Support UniFrac and pre-calculated dissimilarity in getPERMANOVA (#844)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: mia
 Type: Package
-Version: 1.21.9
+Version: 1.21.10
 Authors@R:
     c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"),
              email = "tuomas.v.borman@utu.fi",

R/addDissimilarity.R

@@ -57,7 +57,7 @@
 #'   links, you can provide a vector with whose length equals to the number of
 #'   rows/columns in \code{x}. Alternatively, you can provide a named vector
 #'   where \code{names} represent names in abundance table and values their
-#'   corresponding node in tree. 
+#'   corresponding node in tree.
 #'
 #'   \item \code{chunkSize}: (JSD) \code{Integer scalar}. Defines the size of
 #'   data  send to the individual worker. Only has an effect, if \code{BPPARAM}
@@ -130,7 +130,7 @@
 #'
 #' Lozupone C, Knight R. ``Unifrac: a new phylogenetic method for comparing
 #' microbial communities.'' Appl Environ Microbiol. 2005 71 (12):8228-35.
-#' 
+#'
 #' McDonald D, Vázquez-Baeza Y, Koslicki D, McClelland J, Reeve N, Xu Z,
 #' Gonzalez A, Knight R. ``Striped UniFrac: enabling microbiome analysis at
 #' unprecedented scale.'' Nat Methods. 2018 15 (11):847-848.
@@ -160,12 +160,12 @@
 #' ### Overlap dissimilarity
 #'
 #' tse <- addDissimilarity(tse, method = "overlap", detection = 0.25)
-#' metadata(tse)[["overlap"]][1:6, 1:6]
+#' metadata(tse)[["overlap"]]
 #'
 #' ### JSD dissimilarity
 #'
 #' tse <- addDissimilarity(tse, method = "jsd")
-#' metadata(tse)[["jsd"]][1:6, 1:6]
+#' metadata(tse)[["jsd"]]
 #'
 #' # Multi Dimensional Scaling applied to JSD dissimilarity matrix
 #' tse <- addMDS(tse, method = "overlap", assay.type = "counts")
@@ -177,15 +177,15 @@
 #' dim(as.matrix(res))
 #'
 #' tse <- addDissimilarity(tse, method = "unifrac", weighted = TRUE)
-#' metadata(tse)[["unifrac"]][1:6, 1:6]
+#' metadata(tse)[["unifrac"]]
 #'
 #' ### Bray dissimilarity
 #'
 #' # Bray is usually applied to relative abundances so we have to apply
 #' # transformation first
 #' tse <- transformAssay(tse, method = "relabundance")
 #' res <- getDissimilarity(tse, method = "bray", assay.type = "relabundance")
-#' as.matrix(res)[1:6, 1:6]
+#' res
 #'
 #' # If applying rarefaction, the input must be count matrix and transformation
 #' # method specified in function call (Note: increase niter)
@@ -194,7 +194,7 @@
 #' }
 #' res <- getDissimilarity(
 #'     tse, method = "euclidean", transf = rclr, niter = 2L)
-#' as.matrix(res)[1:6, 1:6]
+#' res
 #'
 NULL
 
@@ -206,7 +206,7 @@ setMethod(
     #
     res <- getDissimilarity(x, method = method, ...)
     # Add matrix to original SE
-    x <- .add_values_to_metadata(x, names = name, values = as.matrix(res))
+    x <- .add_values_to_metadata(x, names = name, values = res)
     return(x)
     }
 )

R/getPERMANOVA.R

@@ -45,9 +45,9 @@
 #' \code{test.homogeneity = TRUE}. \code{addPERMANOVA} adds these results to
 #' metadata of \code{x}.
 #'
-#' @inheritParams addAlpha
 #' @inheritParams runCCA
-#'
+#' @inheritParams addAlpha
+
 #' @param name \code{Character scalar}. A name for the results that will be
 #' stored to metadata. (Default: \code{"permanova"})
 #'
@@ -105,6 +105,17 @@
 #' # Significance results are similar to PERMANOVA
 #' attr(rda_res, "significance")
 #'
+#' # Run PERMANOVA using a pre-calculated dissimilarity matrix from metadata.
+#' # By doing like this, you can use all dissimilarity methods from
+#' # getDissimilarity() in PERMANOVA.
+#' tse <- addDissimilarity(tse, method = "unifrac", name = "unifrac_dist")
+#' res_precalc <- getPERMANOVA(
+#'     tse,
+#'     dis.name = "unifrac_dist",
+#'     formula = x ~ SampleType,
+#'     permutations = 99
+#' )
+#'
 #' @seealso
 #' For more details on the actual implementation see
 #' \code{\link[vegan:adonis2]{vegan::adonis2}},
@@ -131,10 +142,9 @@ setMethod("getPERMANOVA", "SingleCellExperiment", function(x,  ...){
 #' @rdname getPERMANOVA
 setMethod("getPERMANOVA", "SummarizedExperiment",
     function(
-        x, assay.type = "counts", formula = NULL, col.var = NULL, ...){
+        x, assay.type = NULL, dis.name = NULL, formula = NULL, col.var = NULL,
+         ...){
         ############################# Input check ##############################
-        # Assay must be present
-        .check_assay_present(assay.type, x)
         # Formula must be either correctly specified or not specified
         if( !(is.null(formula) || is(formula, "formula")) ){
             stop("'formula' must be formula or NULL.", call. = FALSE)
@@ -145,14 +155,32 @@ setMethod("getPERMANOVA", "SummarizedExperiment",
             stop("'col.var' must specify column from colData(x) or be NULL.",
                 call. = FALSE)
         }
+        if( !is.null(formula) && !is.null(col.var) ){
+            stop("Specify either 'formula' or 'col.var'.", call. = FALSE)
+        }
+        # User can specify either abundance matrix or dissimilarity matrix from
+        # metadata.
+        if( sum(c(is.null(assay.type), is.null(dis.name))) != 1L ){
+            stop("Either 'assay.type' or 'dis.name' must be specified.",
+                 call. = FALSE)
+        }
+        if( !is.null(assay.type) ){
+            # Get assay
+            .check_assay_present(assay.type, x)
+            mat <- assay(x, assay.type)
+        } else{
+            # Get dissimilarity matrix
+            .check_metadata_present(dis.name, x)
+            mat <- metadata(x)[[dis.name]]
+        }
         ########################### Input check end ############################
-        # Get abundance table, formula and related sample metadata as DF
-        mat <- assay(x, assay.type)
+        # Get formula and related sample metadata as DF
         temp <- .get_formula_and_covariates(x, formula, col.var)
         formula <- temp[["formula"]]
         covariates <- temp[["variables"]]
-        # Calculate PERMANOVA with matrix method
-        res <- getPERMANOVA(mat, formula = formula, data = covariates, ...)
+        # Calculate PERMANOVA with matrix/dist method
+        res <- getPERMANOVA(
+            mat, formula = formula, data = covariates, ...)
         return(res)
     }
 )
@@ -161,17 +189,19 @@ setMethod("getPERMANOVA", "SummarizedExperiment",
 #' @rdname getPERMANOVA
 setMethod("getPERMANOVA", "ANY", function(
         x, formula, data, method = "bray", test.homogeneity = TRUE, ...){
-    if( !is.matrix(x) ){
-        stop("'x' must be matrix.", call. = FALSE)
+    is_dist <- inherits(x, "dist")
+    if( !(is.matrix(x) || is_dist) ){
+        stop("'x' must be a matrix or a 'dist' object.", call. = FALSE)
     }
     if( !is(formula, "formula") ){
         stop("'formula' must be formula or NULL.", call. = FALSE)
     }
     if( !(is.data.frame(data) || is.matrix(data) || is(data, "DFrame")) ){
-        stop("'data' must be data.frame or coarcible to one.", call. = FALSE)
+        stop("'data' must be data.frame or coercible to one.", call. = FALSE)
     }
-    if( ncol(x) != nrow(data) ){
-        stop("Number of columns in 'x' should match with number of rows in ",
+    n_samples <- if( is_dist ) attr(x, "Size") else ncol(x)
+    if( n_samples != nrow(data) ){
+        stop("Number of samples in 'x' should match with number of rows in ",
             "'data'.", call. = FALSE)
     }
     if( !.is_a_string(method) ){
@@ -200,7 +230,7 @@ setMethod("addPERMANOVA", "SummarizedExperiment",
     }
     # Calculate permanova
     res <- getPERMANOVA(x, ...)
-    # Addresults to metadata
+    # Add results to metadata
     x <- .add_values_to_metadata(x, name, res)
     return(x)
     }
@@ -209,7 +239,7 @@ setMethod("addPERMANOVA", "SummarizedExperiment",
 ################################ HELP FUNCTIONS ################################
 
 # This function is internal function to perform PERMANOVA from abundance
-# matrix, formula and sample metadata table.
+# matrix or distance matrix, formula and sample metadata table.
 #' @importFrom vegan adonis2
 .calculate_permanova <- function(
         x, formula, data, by = "margin", na.action = na.fail, ...){
@@ -219,8 +249,10 @@ setMethod("addPERMANOVA", "SummarizedExperiment",
     }
     #
     # Get abundance data into correct orientation. Samples must be in rows and
-    # features in columns. Also ensure that the abundance table is matrix.
-    x <- as.matrix(t(x))
+    # features in columns.
+    if( is.matrix(x) ){
+        x <- x |> t()
+    }
     # Covariate table must be data.frame
     data <- data.frame(data, check.names = FALSE)
     # Instead of letting na.action pass through, give informative error
@@ -230,7 +262,7 @@ setMethod("addPERMANOVA", "SummarizedExperiment",
             "to remove samples with missing values.", call. = FALSE)
     }
     # This next step ensures that the formula points correctly to abundance
-    # table
+    # table or distance matrix
     formula <- as.character(formula)
     formula[[2]] <- "x"
     formula <- as.formula(

R/runCCA.R

@@ -674,10 +674,10 @@ setMethod("addRDA", "SingleCellExperiment",
             paste0(names(cols)[!cols], collapse = "','"), "'", call. = FALSE)
         variables <- variables[, cols, drop = FALSE]
     }
-    #
     # Calculate dissimilarity matrix
-    mat <- t(mat)
-    diss_mat <- vegdist(mat, method = method, ...)
+    if( is.matrix(mat) ){
+        mat <- vegdist(mat |> t(), method = method, ...)
+    }
     # For all variables run the analysis
     homogeneity <- lapply(colnames(variables), function(x){
         # Get variable values
@@ -688,7 +688,7 @@ setMethod("addRDA", "SingleCellExperiment",
         # "missing observations due to 'group' removed"
         suppressWarnings(
             suppressMessages(
-                betadisper_res <- betadisper(diss_mat, group = var)
+                betadisper_res <- betadisper(mat, group = var)
             )
         )
         # Run significance test

tests/testthat/test-5calculateOverlap.R

@@ -2,7 +2,7 @@
 context("Overlap")
 
 test_that("Overlap", {
-    
+
     # Get data
     data(esophagus, package="mia")
     tse <- esophagus
@@ -11,7 +11,7 @@ test_that("Overlap", {
         getDissimilarity(tse, method = "overlap", assay.type = "relabundance", detection = 0.15))
     expect_error(getDissimilarity(tse, method = "overlap", detection = "TEST"))
     expect_error(getDissimilarity(tse, method = "overlap", detection = TRUE))
-    
+
     # Calculate overlap
     tse <- transformAssay(tse, method = "relabundance")
     result <- getDissimilarity(tse, method = "overlap", assay.type = "relabundance", detection = 0.15)
@@ -20,24 +20,18 @@ test_that("Overlap", {
     expect_true(class(result) == "dist")
     # Test values
     result <- as.matrix(result)
-    
+
     # Reference
     reference <- matrix(c(0.0000000, 0.3799961, 0.4111838, 0.3799961, 0.0000000,
                           0.3634972, 0.4111838, 0.3634972, 0.0000000), nrow=3)
     expect_equal(round(unname(result), 7), round(reference), 7)
     # Test with different detection threshold
     result <- getDissimilarity(tse, method = "overlap", assay.type = "relabundance", detection = 0)
     result <- as.matrix(result)
-    
+
     # Reference
-    reference <- matrix(c(0.0000000, 0.8811552, 0.9038734, 0.8811552, 0.0000000, 
+    reference <- matrix(c(0.0000000, 0.8811552, 0.9038734, 0.8811552, 0.0000000,
                           0.8390008, 0.9038734, 0.8390008, 0.00000000), nrow=3)
     expect_equal(round(unname(result), 7), round(reference, 7))
-
-    # Test .calculate_overlap
-    expect_equal( metadata(addDissimilarity(tse, "overlap"))[["overlap"]], 
-                  as.matrix(getDissimilarity(tse, method = "overlap")) )
-    expect_equal( metadata(addDissimilarity(tse, method = "overlap", assay.type = "relabundance", detection = 0.08))[["overlap"]], 
-                  as.matrix(getDissimilarity(tse, method = "overlap", assay.type = "relabundance", detection = 0.08)) )
 })