From 95bc38991eff46a9119fff90cb73eb58dabc6ed9 Mon Sep 17 00:00:00 2001 From: Daena Rys Date: Tue, 8 Sep 2026 12:51:38 +0300 Subject: [PATCH 1/5] Support UniFrac and pre-calculated dissimilarity in getPERMANOVA (#844) - Add dis.name to getPERMANOVA,SummarizedExperiment to retrieve pre-calculated dissimilarity matrices from metadata(x). - Automatically delegate to mia::getDissimilarity() when method %in% c('unifrac', 'overlap', 'jsd'). - Allow getPERMANOVA,ANY to accept dist objects directly and align sample metadata rows by sample names. - Update .calculate_homogeneity() in runCCA.R to accept dist and symmetric matrices without re-running vegdist(). - Add unit tests covering on-the-fly UniFrac, direct dist inputs, dis.name usage, and input validations. - Synchronize documentation and examples in man/getPERMANOVA.Rd. --- R/getPERMANOVA.R | 103 ++++++++++++++++++++++++----- R/runCCA.R | 11 ++- man/getPERMANOVA.Rd | 32 ++++++++- tests/testthat/test-getPERMANOVA.R | 97 ++++++++++++++++++++++++++- 4 files changed, 220 insertions(+), 23 deletions(-) diff --git a/R/getPERMANOVA.R b/R/getPERMANOVA.R index 7c5a5c786..57d43a996 100644 --- a/R/getPERMANOVA.R +++ b/R/getPERMANOVA.R @@ -54,6 +54,9 @@ #' @param method \code{Character scalar}. A dissimilarity metric used in #' PERMANOVA and group dispersion calculation. (Default: \code{"bray"}) #' +#' @param dis.name \code{Character scalar}. Specifies the name of dissimilarity +#' matrix from \code{metadata(x)} used in calculation. (Default: \code{NULL}) +#' #' @param test.homogeneity \code{Logical scalar}. Should the homogeneity of #' group dispersions be evaluated? (Default: \code{TRUE}) #' @@ -105,6 +108,25 @@ #' # Significance results are similar to PERMANOVA #' attr(rda_res, "significance") #' +#' # Perform PERMANOVA with UniFrac directly on TreeSE +#' tse <- addPERMANOVA( +#' tse, +#' assay.type = "counts", +#' method = "unifrac", +#' formula = x ~ SampleType, +#' name = "unifrac_permanova", +#' permutations = 99 +#' ) +#' +#' # Run PERMANOVA using a pre-calculated dissimilarity matrix from metadata +#' 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}}, @@ -128,10 +150,9 @@ setMethod("getPERMANOVA", "SingleCellExperiment", function(x, ...){ #' @rdname getPERMANOVA setMethod("getPERMANOVA", "SummarizedExperiment", function( - x, assay.type = "counts", formula = NULL, col.var = NULL, ...){ + x, assay.type = "counts", formula = NULL, col.var = NULL, + dis.name = NULL, method = "bray", ...){ ############################# 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) @@ -142,14 +163,45 @@ 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) + } + if( !.is_a_string(method) ){ + stop("'method' must be a single character value.", call. = FALSE) + } + # User can specify either abundance matrix or dissimilarity matrix from + # metadata. + if( !is.null(dis.name) ){ + if( !.is_a_string(dis.name) ){ + stop("'dis.name' must be a single character value.", + call. = FALSE) + } + .check_metadata_present(dis.name, x) + mat <- metadata(x)[[dis.name]] + if( is.matrix(mat) ){ + mat <- as.dist(mat) + } + if( !inherits(mat, "dist") ){ + stop("'", dis.name, "' in metadata(x) must be a 'dist' object or ", + "a symmetric distance matrix.", call. = FALSE) + } + } else{ + .check_assay_present(assay.type, x) + if( method %in% c("unifrac", "overlap", "jsd") ){ + mat <- getDissimilarity( + x, method = method, assay.type = assay.type, ...) + } else{ + mat <- assay(x, assay.type) + } + } ########################### 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, method = method, ...) return(res) } ) @@ -158,18 +210,31 @@ 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) + } + n_samples <- if( is_dist ) attr(x, "Size") else ncol(x) + if( n_samples != nrow(data) ){ + if( is_dist ){ + stop("Number of samples in 'x' should match with number of rows in ", + "'data'.", call. = FALSE) + } else { + stop("Number of columns in 'x' should match with number of rows in ", + "'data'.", call. = FALSE) + } } - if( ncol(x) != nrow(data) ){ - stop("Number of columns in 'x' should match with number of rows in ", - "'data'.", call. = FALSE) + # Align rows of data if names match + sample_names <- if( is_dist ) labels(x) else colnames(x) + if( !is.null(sample_names) && !is.null(rownames(data)) && + setequal(sample_names, rownames(data)) ){ + data <- data[sample_names, , drop = FALSE] } if( !.is_a_string(method) ){ stop("'method' must be a single character value.", call. = FALSE) @@ -197,7 +262,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) } @@ -206,7 +271,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, ...){ @@ -216,8 +281,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( !inherits(x, "dist") ){ + x <- as.matrix(t(x)) + } # Covariate table must be data.frame data <- data.frame(data, check.names = FALSE) # Instead of letting na.action pass through, give informative error @@ -227,7 +294,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( diff --git a/R/runCCA.R b/R/runCCA.R index af2ba2446..26b068d7e 100644 --- a/R/runCCA.R +++ b/R/runCCA.R @@ -674,10 +674,15 @@ 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( inherits(mat, "dist") ){ + diss_mat <- mat + } else if( is.matrix(mat) && isSymmetric(unname(mat)) ){ + diss_mat <- as.dist(mat) + } else { + mat <- t(mat) + diss_mat <- vegdist(mat, method = method, ...) + } # For all variables run the analysis homogeneity <- lapply(colnames(variables), function(x){ # Get variable values diff --git a/man/getPERMANOVA.Rd b/man/getPERMANOVA.Rd index ccb3c4574..8358c7b5b 100644 --- a/man/getPERMANOVA.Rd +++ b/man/getPERMANOVA.Rd @@ -15,7 +15,15 @@ addPERMANOVA(x, ...) \S4method{getPERMANOVA}{SingleCellExperiment}(x, ...) -\S4method{getPERMANOVA}{SummarizedExperiment}(x, assay.type = "counts", formula = NULL, col.var = NULL, ...) +\S4method{getPERMANOVA}{SummarizedExperiment}( + x, + assay.type = "counts", + formula = NULL, + col.var = NULL, + dis.name = NULL, + method = "bray", + ... +) \S4method{getPERMANOVA}{ANY}(x, formula, data, method = "bray", test.homogeneity = TRUE, ...) @@ -68,6 +76,9 @@ including covariates defined by \code{formula}.} \item{method}{\code{Character scalar}. A dissimilarity metric used in PERMANOVA and group dispersion calculation. (Default: \code{"bray"})} +\item{dis.name}{\code{Character scalar}. Specifies the name of dissimilarity +matrix from \code{metadata(x)} used in calculation. (Default: \code{NULL})} + \item{test.homogeneity}{\code{Logical scalar}. Should the homogeneity of group dispersions be evaluated? (Default: \code{TRUE})} @@ -140,6 +151,25 @@ rda_res <- getRDA( # Significance results are similar to PERMANOVA attr(rda_res, "significance") +# Perform PERMANOVA with UniFrac directly on TreeSE +tse <- addPERMANOVA( + tse, + assay.type = "counts", + method = "unifrac", + formula = x ~ SampleType, + name = "unifrac_permanova", + permutations = 99 +) + +# Run PERMANOVA using a pre-calculated dissimilarity matrix from metadata +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 diff --git a/tests/testthat/test-getPERMANOVA.R b/tests/testthat/test-getPERMANOVA.R index 0e586264d..1986f57b7 100644 --- a/tests/testthat/test-getPERMANOVA.R +++ b/tests/testthat/test-getPERMANOVA.R @@ -189,7 +189,7 @@ test_that("getPERMANOVA matches direct calculations", { permutations = 99 ) homogeneity_direct <- vegan::betadisper( - vegdist(t(assay(tse, "relabundance"))), + vegan::vegdist(t(assay(tse, "relabundance"))), group = tse$SampleType ) @@ -209,3 +209,98 @@ test_that("getPERMANOVA matches direct calculations", { expect_equal( res[[2]][[2]][[1]][[1]]$distances, homogeneity_direct$distances) }) + +test_that("getPERMANOVA works with UniFrac and dist object directly", { + # Test on-the-fly UniFrac calculation on TreeSummarizedExperiment + set.seed(42) + res_unifrac <- getPERMANOVA( + tse, assay.type = "counts", + method = "unifrac", + formula = x ~ SampleType, + test.homogeneity = TRUE, + permutations = 99 + ) + expect_type(res_unifrac, "list") + expect_s3_class(res_unifrac$permanova, "anova.cca") + expect_s3_class(res_unifrac$homogeneity, "data.frame") + + # Test passing a dist object directly + d <- getDissimilarity(tse, method = "unifrac") + expect_s3_class(d, "dist") + setseed_wrap <- set.seed(42) + res_dist <- suppressWarnings(getPERMANOVA( + d, formula = x ~ SampleType, data = colData(tse), + test.homogeneity = TRUE, permutations = 99 + )) + expect_type(res_dist, "list") + expect_s3_class(res_dist$permanova, "anova.cca") + expect_equal(res_unifrac$permanova$SumOfSqs, res_dist$permanova$SumOfSqs) + expect_equal(res_unifrac$permanova$R2, res_dist$permanova$R2) + expect_equal(res_unifrac$permanova$F, res_dist$permanova$F) +}) + +test_that("getPERMANOVA and addPERMANOVA work with pre-calculated dissimilarity (dis.name)", { + # Pre-calculate dissimilarity and store in metadata + tse_diss <- addDissimilarity(tse, method = "unifrac", name = "unifrac_dist") + expect_true("unifrac_dist" %in% names(metadata(tse_diss))) + + # Run getPERMANOVA using dis.name + set.seed(42) + res_precalc <- getPERMANOVA( + tse_diss, dis.name = "unifrac_dist", + formula = x ~ SampleType, permutations = 99 + ) + expect_s3_class(res_precalc$permanova, "anova.cca") + + d <- getDissimilarity(tse, method = "unifrac") + set.seed(42) + res_dist <- suppressWarnings(getPERMANOVA( + d, formula = x ~ SampleType, data = colData(tse), permutations = 99 + )) + expect_equal(res_precalc$permanova$SumOfSqs, res_dist$permanova$SumOfSqs) + expect_equal(res_precalc$permanova$F, res_dist$permanova$F) + + # Test addPERMANOVA with dis.name + tse_meta <- addPERMANOVA( + tse_diss, dis.name = "unifrac_dist", + formula = x ~ SampleType, name = "perm_precalc", permutations = 99 + ) + expect_true("perm_precalc" %in% names(metadata(tse_meta))) + expect_s3_class(metadata(tse_meta)[["perm_precalc"]]$permanova, "anova.cca") + + # Test addPERMANOVA with method = unifrac on-the-fly + tse_meta2 <- addPERMANOVA( + tse, method = "unifrac", + formula = x ~ SampleType, name = "perm_unifrac", permutations = 99 + ) + expect_true("perm_unifrac" %in% names(metadata(tse_meta2))) + expect_s3_class(metadata(tse_meta2)[["perm_unifrac"]]$permanova, "anova.cca") +}) + +test_that("getPERMANOVA input validations for dis.name and dist objects", { + # Nonexistent dis.name in metadata + expect_error( + getPERMANOVA(tse, dis.name = "nonexistent", formula = x ~ SampleType), + "'dis.name' must be a valid name of metadata" + ) + + # Invalid dis.name type + expect_error( + getPERMANOVA(tse, dis.name = 123, formula = x ~ SampleType), + "'dis.name' must be a single character value" + ) + + # Dist object with sample size mismatch + d <- getDissimilarity(tse, method = "bray") + expect_error( + getPERMANOVA(d, formula = x ~ SampleType, data = colData(tse)[1:10, ]), + "Number of samples in 'x' should match" + ) + + # Invalid x (neither matrix nor dist) + expect_error( + getPERMANOVA(1:10, formula = x ~ SampleType, data = colData(tse)), + "'x' must be a matrix or a 'dist' object" + ) +}) + From 3d8e0b17c7a8ab2aed67072160471cde83687dce Mon Sep 17 00:00:00 2001 From: Tuomas Borman Date: Fri, 25 Sep 2026 09:50:32 +0300 Subject: [PATCH 2/5] up --- DESCRIPTION | 2 +- R/addDissimilarity.R | 6 +- R/getPERMANOVA.R | 79 +++++------------- R/runCCA.R | 11 +-- tests/testthat/test-5getDissimilarity.R | 6 +- tests/testthat/test-getPERMANOVA.R | 102 +++++++----------------- 6 files changed, 61 insertions(+), 145 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a4c7e7952..a6547382d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mia Type: Package -Version: 1.21.7 +Version: 1.21.9 Authors@R: c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"), email = "tuomas.v.borman@utu.fi", diff --git a/R/addDissimilarity.R b/R/addDissimilarity.R index 844d80478..323b6ec42 100644 --- a/R/addDissimilarity.R +++ b/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. @@ -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) } ) diff --git a/R/getPERMANOVA.R b/R/getPERMANOVA.R index 57d43a996..78c8e4528 100644 --- a/R/getPERMANOVA.R +++ b/R/getPERMANOVA.R @@ -45,18 +45,15 @@ #' \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"}) #' #' @param method \code{Character scalar}. A dissimilarity metric used in #' PERMANOVA and group dispersion calculation. (Default: \code{"bray"}) #' -#' @param dis.name \code{Character scalar}. Specifies the name of dissimilarity -#' matrix from \code{metadata(x)} used in calculation. (Default: \code{NULL}) -#' #' @param test.homogeneity \code{Logical scalar}. Should the homogeneity of #' group dispersions be evaluated? (Default: \code{TRUE}) #' @@ -108,17 +105,9 @@ #' # Significance results are similar to PERMANOVA #' attr(rda_res, "significance") #' -#' # Perform PERMANOVA with UniFrac directly on TreeSE -#' tse <- addPERMANOVA( -#' tse, -#' assay.type = "counts", -#' method = "unifrac", -#' formula = x ~ SampleType, -#' name = "unifrac_permanova", -#' permutations = 99 -#' ) -#' -#' # Run PERMANOVA using a pre-calculated dissimilarity matrix from metadata +#' # 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, @@ -150,8 +139,8 @@ setMethod("getPERMANOVA", "SingleCellExperiment", function(x, ...){ #' @rdname getPERMANOVA setMethod("getPERMANOVA", "SummarizedExperiment", function( - x, assay.type = "counts", formula = NULL, col.var = NULL, - dis.name = NULL, method = "bray", ...){ + x, assay.type = NULL, dis.name = NULL, formula = NULL, col.var = NULL, + ...){ ############################# Input check ############################## # Formula must be either correctly specified or not specified if( !(is.null(formula) || is(formula, "formula")) ){ @@ -166,33 +155,20 @@ setMethod("getPERMANOVA", "SummarizedExperiment", if( !is.null(formula) && !is.null(col.var) ){ stop("Specify either 'formula' or 'col.var'.", call. = FALSE) } - if( !.is_a_string(method) ){ - stop("'method' must be a single character value.", call. = FALSE) - } # User can specify either abundance matrix or dissimilarity matrix from # metadata. - if( !is.null(dis.name) ){ - if( !.is_a_string(dis.name) ){ - stop("'dis.name' must be a single character value.", - call. = FALSE) - } + 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]] - if( is.matrix(mat) ){ - mat <- as.dist(mat) - } - if( !inherits(mat, "dist") ){ - stop("'", dis.name, "' in metadata(x) must be a 'dist' object or ", - "a symmetric distance matrix.", call. = FALSE) - } - } else{ - .check_assay_present(assay.type, x) - if( method %in% c("unifrac", "overlap", "jsd") ){ - mat <- getDissimilarity( - x, method = method, assay.type = assay.type, ...) - } else{ - mat <- assay(x, assay.type) - } } ########################### Input check end ############################ # Get formula and related sample metadata as DF @@ -201,7 +177,7 @@ setMethod("getPERMANOVA", "SummarizedExperiment", covariates <- temp[["variables"]] # Calculate PERMANOVA with matrix/dist method res <- getPERMANOVA( - mat, formula = formula, data = covariates, method = method, ...) + mat, formula = formula, data = covariates, ...) return(res) } ) @@ -222,19 +198,8 @@ setMethod("getPERMANOVA", "ANY", function( } n_samples <- if( is_dist ) attr(x, "Size") else ncol(x) if( n_samples != nrow(data) ){ - if( is_dist ){ - stop("Number of samples in 'x' should match with number of rows in ", - "'data'.", call. = FALSE) - } else { - stop("Number of columns in 'x' should match with number of rows in ", - "'data'.", call. = FALSE) - } - } - # Align rows of data if names match - sample_names <- if( is_dist ) labels(x) else colnames(x) - if( !is.null(sample_names) && !is.null(rownames(data)) && - setequal(sample_names, rownames(data)) ){ - data <- data[sample_names, , drop = FALSE] + stop("Number of samples in 'x' should match with number of rows in ", + "'data'.", call. = FALSE) } if( !.is_a_string(method) ){ stop("'method' must be a single character value.", call. = FALSE) @@ -282,8 +247,8 @@ setMethod("addPERMANOVA", "SummarizedExperiment", # # Get abundance data into correct orientation. Samples must be in rows and # features in columns. - if( !inherits(x, "dist") ){ - x <- as.matrix(t(x)) + if( is.matrix(x) ){ + x <- x |> t() } # Covariate table must be data.frame data <- data.frame(data, check.names = FALSE) diff --git a/R/runCCA.R b/R/runCCA.R index 26b068d7e..c06486f4d 100644 --- a/R/runCCA.R +++ b/R/runCCA.R @@ -675,13 +675,8 @@ setMethod("addRDA", "SingleCellExperiment", variables <- variables[, cols, drop = FALSE] } # Calculate dissimilarity matrix - if( inherits(mat, "dist") ){ - diss_mat <- mat - } else if( is.matrix(mat) && isSymmetric(unname(mat)) ){ - diss_mat <- as.dist(mat) - } else { - 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){ @@ -693,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 diff --git a/tests/testthat/test-5getDissimilarity.R b/tests/testthat/test-5getDissimilarity.R index bd4752b58..6dafeab94 100644 --- a/tests/testthat/test-5getDissimilarity.R +++ b/tests/testthat/test-5getDissimilarity.R @@ -18,13 +18,13 @@ test_that("Dissimilarity calculation", { } set.seed(123) res1 <- vegan::avgdist(t(mat), distfun = vegdist, dmethod = "euclidean", - sample = min(colSums2(mat)), + sample = min(colSums2(mat)), iterations = 10, transf = clr) |> expect_warning() set.seed(123) - res2 <- getDissimilarity(tse_sub, method = "euclidean", niter = 10, + res2 <- getDissimilarity(tse_sub, method = "euclidean", niter = 10, transf = clr) |> expect_warning() expect_equal(as.matrix(res1), as.matrix(res2)) - + # Test unifrac with named vector as node labels rownames(tse) <- paste0("taxa", seq_len(nrow(tse))) nodes <- rowLinks(tse)[["nodeLab"]] diff --git a/tests/testthat/test-getPERMANOVA.R b/tests/testthat/test-getPERMANOVA.R index 1986f57b7..9cdc35a58 100644 --- a/tests/testthat/test-getPERMANOVA.R +++ b/tests/testthat/test-getPERMANOVA.R @@ -13,7 +13,7 @@ test_that("getPERMANOVA works on SummarizedExperiment", { permutations = 99 ) expect_s3_class(res, "anova.cca") - + # Test PERMANOVA with homogeneity check enabled res <- getPERMANOVA( tse, assay.type = "relabundance", @@ -26,7 +26,7 @@ test_that("getPERMANOVA works on SummarizedExperiment", { expect_true("homogeneity" %in% names(res)) expect_s3_class(res$permanova, "anova.cca") expect_s3_class(res$homogeneity, "data.frame") - + # Test full results with nested structure validation for detailed outputs res <- getPERMANOVA( tse, assay.type = "relabundance", @@ -70,7 +70,7 @@ test_that("getPERMANOVA input validations", { getPERMANOVA( assay(tse[, 1:10]), formula = x ~ SampleType, data = colData(tse) ), - "Number of columns in 'x' should match with number of rows in 'data'" + "Number of samples in 'x' should match with number of rows in 'data'" ) # Check that invalid homogeneity test settings raise an error expect_error( @@ -81,37 +81,37 @@ test_that("getPERMANOVA input validations", { ), "'test.homogeneity' must be TRUE or FALSE" ) - + # Invalid assay type (non-existent type) expect_error( getPERMANOVA(tse, assay.type = "nonexistent", formula = x ~ SampleType) ) - + # Invalid variable expect_error( getPERMANOVA(tse, assay.type = "relabundance", col.var = "invalid") ) - + # Incorrect permutations expect_error( getPERMANOVA( - tse, assay.type = "relabundance", formula = x ~ SampleType, + tse, assay.type = "relabundance", formula = x ~ SampleType, permutations = -10 ) ) - + # Incorrect 'by' parameter expect_error( getPERMANOVA( - tse, assay.type = "relabundance", formula = x ~ SampleType, + tse, assay.type = "relabundance", formula = x ~ SampleType, by = "invalid_option" ) ) - + # Missing or incorrect 'homogeneity.test' option expect_error( getPERMANOVA( - tse, assay.type = "relabundance", formula = x ~ SampleType, + tse, assay.type = "relabundance", formula = x ~ SampleType, homogeneity.test = "unsupported_test" ) ) @@ -136,7 +136,7 @@ test_that("getPERMANOVA 'by' and 'homogeneity.test' options", { formula = x ~ SampleType, by = "terms" ) expect_s3_class(res_by_terms$permanova, "anova.cca") - + # Test homogeneity test with ANOVA option res_anova <- getPERMANOVA( tse, assay.type = "relabundance", @@ -144,7 +144,7 @@ test_that("getPERMANOVA 'by' and 'homogeneity.test' options", { homogeneity.test = "anova" ) expect_s3_class(res_anova$homogeneity, "data.frame") - + # Test homogeneity test with Tukey HSD and full results res_tukey <- getPERMANOVA( tse, assay.type = "relabundance", @@ -158,8 +158,8 @@ test_that("getPERMANOVA 'by' and 'homogeneity.test' options", { test_that("getPERMANOVA handles edge cases", { # Test handling of a missing formula (default behavior) - expect_no_error(getPERMANOVA(tse, assay.type = "relabundance")) - + expect_error(getPERMANOVA(tse, assay.type = "relabundance")) + # Test for error when permutations count is zero expect_error( getPERMANOVA( @@ -167,10 +167,10 @@ test_that("getPERMANOVA handles edge cases", { formula = x ~ SampleType, permutations = 0 ) ) - + # Test for error when input matrix is NULL expect_error(getPERMANOVA(NULL, formula = x ~ SampleType)) - + # Test for warning when only one level exists in the factor variable tse_subset <- tse[, tse$SampleType == "Soil"] expect_warning( @@ -192,7 +192,7 @@ test_that("getPERMANOVA matches direct calculations", { vegan::vegdist(t(assay(tse, "relabundance"))), group = tse$SampleType ) - + # Run the getPERMANOVA function and compare results res <- getPERMANOVA( tse, assay.type = "relabundance", @@ -201,80 +201,36 @@ test_that("getPERMANOVA matches direct calculations", { full = TRUE, permutations = 99 ) - + # Verify permanova results match expect_equal(res$permanova$aov.tab, permanova_direct$aov.tab) - + # Verify homogeneity results match expect_equal( res[[2]][[2]][[1]][[1]]$distances, homogeneity_direct$distances) }) -test_that("getPERMANOVA works with UniFrac and dist object directly", { - # Test on-the-fly UniFrac calculation on TreeSummarizedExperiment - set.seed(42) - res_unifrac <- getPERMANOVA( - tse, assay.type = "counts", - method = "unifrac", - formula = x ~ SampleType, - test.homogeneity = TRUE, - permutations = 99 - ) - expect_type(res_unifrac, "list") - expect_s3_class(res_unifrac$permanova, "anova.cca") - expect_s3_class(res_unifrac$homogeneity, "data.frame") - - # Test passing a dist object directly - d <- getDissimilarity(tse, method = "unifrac") - expect_s3_class(d, "dist") - setseed_wrap <- set.seed(42) - res_dist <- suppressWarnings(getPERMANOVA( - d, formula = x ~ SampleType, data = colData(tse), - test.homogeneity = TRUE, permutations = 99 - )) - expect_type(res_dist, "list") - expect_s3_class(res_dist$permanova, "anova.cca") - expect_equal(res_unifrac$permanova$SumOfSqs, res_dist$permanova$SumOfSqs) - expect_equal(res_unifrac$permanova$R2, res_dist$permanova$R2) - expect_equal(res_unifrac$permanova$F, res_dist$permanova$F) -}) - test_that("getPERMANOVA and addPERMANOVA work with pre-calculated dissimilarity (dis.name)", { # Pre-calculate dissimilarity and store in metadata - tse_diss <- addDissimilarity(tse, method = "unifrac", name = "unifrac_dist") - expect_true("unifrac_dist" %in% names(metadata(tse_diss))) + tse <- addDissimilarity(tse, method = "unifrac", name = "unifrac_dist") # Run getPERMANOVA using dis.name set.seed(42) res_precalc <- getPERMANOVA( - tse_diss, dis.name = "unifrac_dist", + tse, dis.name = "unifrac_dist", formula = x ~ SampleType, permutations = 99 ) expect_s3_class(res_precalc$permanova, "anova.cca") - d <- getDissimilarity(tse, method = "unifrac") set.seed(42) - res_dist <- suppressWarnings(getPERMANOVA( - d, formula = x ~ SampleType, data = colData(tse), permutations = 99 - )) + res_dist <- getPERMANOVA( + x = metadata(tse)[["unifrac_dist"]], + formula = x ~ SampleType, + data = colData(tse)[, "SampleType", drop = FALSE], + permutations = 99 + ) expect_equal(res_precalc$permanova$SumOfSqs, res_dist$permanova$SumOfSqs) expect_equal(res_precalc$permanova$F, res_dist$permanova$F) - - # Test addPERMANOVA with dis.name - tse_meta <- addPERMANOVA( - tse_diss, dis.name = "unifrac_dist", - formula = x ~ SampleType, name = "perm_precalc", permutations = 99 - ) - expect_true("perm_precalc" %in% names(metadata(tse_meta))) - expect_s3_class(metadata(tse_meta)[["perm_precalc"]]$permanova, "anova.cca") - - # Test addPERMANOVA with method = unifrac on-the-fly - tse_meta2 <- addPERMANOVA( - tse, method = "unifrac", - formula = x ~ SampleType, name = "perm_unifrac", permutations = 99 - ) - expect_true("perm_unifrac" %in% names(metadata(tse_meta2))) - expect_s3_class(metadata(tse_meta2)[["perm_unifrac"]]$permanova, "anova.cca") }) test_that("getPERMANOVA input validations for dis.name and dist objects", { @@ -287,7 +243,7 @@ test_that("getPERMANOVA input validations for dis.name and dist objects", { # Invalid dis.name type expect_error( getPERMANOVA(tse, dis.name = 123, formula = x ~ SampleType), - "'dis.name' must be a single character value" + "'dis.name' must be a single non-empty character value." ) # Dist object with sample size mismatch From 9307a6fe2063dcf279911024bdf247ebb7a6601c Mon Sep 17 00:00:00 2001 From: Tuomas Borman Date: Fri, 25 Sep 2026 09:53:17 +0300 Subject: [PATCH 3/5] up --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index a6547382d..7d6959702 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mia Type: Package -Version: 1.21.9 +Version: 1.21.8 Authors@R: c(person(given = "Tuomas", family = "Borman", role = c("aut", "cre"), email = "tuomas.v.borman@utu.fi", From a94bbe617aa57ddef5c2345a624797b1a2555b09 Mon Sep 17 00:00:00 2001 From: Tuomas Borman Date: Fri, 25 Sep 2026 12:38:20 +0300 Subject: [PATCH 4/5] up --- R/addDissimilarity.R | 10 +++++----- man/addAlpha.Rd | 5 +++++ man/getDissimilarity.Rd | 20 ++++++++++++-------- man/getPERMANOVA.Rd | 29 ++++++++++------------------- man/getRPCA.Rd | 4 ++++ 5 files changed, 36 insertions(+), 32 deletions(-) diff --git a/R/addDissimilarity.R b/R/addDissimilarity.R index 323b6ec42..e2c3bfbb4 100644 --- a/R/addDissimilarity.R +++ b/R/addDissimilarity.R @@ -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,7 +177,7 @@ #' dim(as.matrix(res)) #' #' tse <- addDissimilarity(tse, method = "unifrac", weighted = TRUE) -#' metadata(tse)[["unifrac"]][1:6, 1:6] +#' metadata(tse)[["unifrac"]] #' #' ### Bray dissimilarity #' @@ -185,7 +185,7 @@ #' # 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 diff --git a/man/addAlpha.Rd b/man/addAlpha.Rd index acc06a0f2..723620d39 100644 --- a/man/addAlpha.Rd +++ b/man/addAlpha.Rd @@ -416,6 +416,11 @@ Efficient computation of Faith's phylogenetic diversity with applications in characterizing microbiomes. \emph{Genome Res.} 31(11):2131-2137. doi: 10.1101/gr.275777.121 +Armstrong G. et al. (2021) +Efficient computation of Faith's phylogenetic diversity with applications +in characterizing microbiomes. +\emph{Genome Res.} 31(11):2131-2137. doi: 10.1101/gr.275777.121 + Beisel J-N. et al. (2003) A Comparative Analysis of Diversity Index Sensitivity. \emph{Internal Rev. Hydrobiol.} 88(1):3-15. diff --git a/man/getDissimilarity.Rd b/man/getDissimilarity.Rd index 0769b0e70..7b6fa8822 100644 --- a/man/getDissimilarity.Rd +++ b/man/getDissimilarity.Rd @@ -132,9 +132,8 @@ similarity is. When using relative abundances, overlap value 1 means that all the abundances of features are equal between two samples, and 0 means that samples have completely different relative abundances. -Unifrac is calculated with \code{ -\link[ecodive:unweighted_unifrac]{ecodive:unweighted_unifrac()}} -or \code{\link[ecodive:weighted_unifrac]{ecodive:weighted_unifrac()}}. +Unifrac is calculated with a C++ implementation of the Striped Unifrac +algorithm (McDonald et al. 2021). If rarefaction is enabled, \code{\link[vegan:avgdist]{vegan:avgdist()}} is utilized. @@ -163,12 +162,12 @@ tse <- GlobalPatterns ### 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") @@ -180,7 +179,7 @@ res <- getDissimilarity(tse, method = "unifrac", weighted = FALSE) dim(as.matrix(res)) tse <- addDissimilarity(tse, method = "unifrac", weighted = TRUE) -metadata(tse)[["unifrac"]][1:6, 1:6] +metadata(tse)[["unifrac"]] ### Bray dissimilarity @@ -188,7 +187,7 @@ metadata(tse)[["unifrac"]][1:6, 1:6] # 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) @@ -197,7 +196,7 @@ rclr <- function(x){ } res <- getDissimilarity( tse, method = "euclidean", transf = rclr, niter = 2L) -as.matrix(res)[1:6, 1:6] +res } \references{ @@ -216,6 +215,11 @@ microbial communities.'' Appl Environ Microbiol. 2007 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. +doi: 10.1038/s41592-018-0187-8. + For JSD dissimilarity: Jensen-Shannon Divergence and Hilbert space embedding. Bent Fuglede and Flemming Topsoe University of Copenhagen, diff --git a/man/getPERMANOVA.Rd b/man/getPERMANOVA.Rd index 8358c7b5b..da7fb4f51 100644 --- a/man/getPERMANOVA.Rd +++ b/man/getPERMANOVA.Rd @@ -17,11 +17,10 @@ addPERMANOVA(x, ...) \S4method{getPERMANOVA}{SummarizedExperiment}( x, - assay.type = "counts", + assay.type = NULL, + dis.name = NULL, formula = NULL, col.var = NULL, - dis.name = NULL, - method = "bray", ... ) @@ -30,7 +29,7 @@ addPERMANOVA(x, ...) \S4method{addPERMANOVA}{SummarizedExperiment}(x, name = "permanova", ...) } \arguments{ -\item{x}{a \code{\link[SummarizedExperiment]{SummarizedExperiment}} object.} +\item{x}{\code{\link[TreeSummarizedExperiment:TreeSummarizedExperiment-class]{TreeSummarizedExperiment}}.} \item{...}{additional arguments passed to \code{vegan::adonis2}. \itemize{ @@ -56,7 +55,10 @@ Options include 'permanova' }} \item{assay.type}{\code{Character scalar}. Specifies the name of assay -used in calculation. (Default: \code{"counts"})} +used in calculation. (Default: \code{NULL})} + +\item{dis.name}{\code{Character scalar}. Specifies the name of dissimilarity +matrix from \code{metadata} slot used in calculation. (Default: \code{NULL})} \item{formula}{\code{formula}. If \code{x} is a \code{\link[SummarizedExperiment:SummarizedExperiment-class]{SummarizedExperiment}} @@ -76,9 +78,6 @@ including covariates defined by \code{formula}.} \item{method}{\code{Character scalar}. A dissimilarity metric used in PERMANOVA and group dispersion calculation. (Default: \code{"bray"})} -\item{dis.name}{\code{Character scalar}. Specifies the name of dissimilarity -matrix from \code{metadata(x)} used in calculation. (Default: \code{NULL})} - \item{test.homogeneity}{\code{Logical scalar}. Should the homogeneity of group dispersions be evaluated? (Default: \code{TRUE})} @@ -151,17 +150,9 @@ rda_res <- getRDA( # Significance results are similar to PERMANOVA attr(rda_res, "significance") -# Perform PERMANOVA with UniFrac directly on TreeSE -tse <- addPERMANOVA( - tse, - assay.type = "counts", - method = "unifrac", - formula = x ~ SampleType, - name = "unifrac_permanova", - permutations = 99 -) - -# Run PERMANOVA using a pre-calculated dissimilarity matrix from metadata +# 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, diff --git a/man/getRPCA.Rd b/man/getRPCA.Rd index 0a7a73a33..78fb5ea77 100644 --- a/man/getRPCA.Rd +++ b/man/getRPCA.Rd @@ -125,6 +125,10 @@ reconstructed matrix, and reconstruction error) stored in attributes. data("ibdmdb") mae <- ibdmdb +# Apply filtering +mae[[1]] <- filterRPCAInput(mae[[1]], assay.type = "mgx") +mae[[2]] <- filterRPCAInput(mae[[2]], assay.type = "mtx") + # Apply data transformations. With impute=FALSE, missing values are preserved # and not imputed. mae[[1]] <- transformAssay( From 19d2c8e428f7dfb93dd9367255150937b4e65e60 Mon Sep 17 00:00:00 2001 From: Tuomas Borman Date: Fri, 25 Sep 2026 13:33:36 +0300 Subject: [PATCH 5/5] up --- tests/testthat/test-5calculateOverlap.R | 16 +++++----------- tests/testthat/test-getPERMANOVA.R | 2 +- 2 files changed, 6 insertions(+), 12 deletions(-) diff --git a/tests/testthat/test-5calculateOverlap.R b/tests/testthat/test-5calculateOverlap.R index 9e0c360c9..c200567ac 100644 --- a/tests/testthat/test-5calculateOverlap.R +++ b/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,7 +20,7 @@ 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) @@ -28,16 +28,10 @@ test_that("Overlap", { # 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)) ) }) diff --git a/tests/testthat/test-getPERMANOVA.R b/tests/testthat/test-getPERMANOVA.R index 9cdc35a58..a50197ee2 100644 --- a/tests/testthat/test-getPERMANOVA.R +++ b/tests/testthat/test-getPERMANOVA.R @@ -158,7 +158,7 @@ test_that("getPERMANOVA 'by' and 'homogeneity.test' options", { test_that("getPERMANOVA handles edge cases", { # Test handling of a missing formula (default behavior) - expect_error(getPERMANOVA(tse, assay.type = "relabundance")) + expect_no_error(getPERMANOVA(tse, assay.type = "relabundance")) # Test for error when permutations count is zero expect_error(