From e7653ed3a51bac4aa8d0f87b7e7fa4a6816db0e2 Mon Sep 17 00:00:00 2001 From: Daena Rys Date: Sat, 3 Jan 2026 08:20:10 +0200 Subject: [PATCH] up --- .Rbuildignore | 2 + DESCRIPTION | 16 +- NAMESPACE | 24 +++ R/AllGenerics.R | 25 +++ R/getTtest.R | 149 +++++++++++++ R/getWilcoxon.R | 144 +++++++++++++ R/utils.R | 322 +++++++++++++++++++++++++++++ README.Rmd | 97 +++------ README.md | 148 +++---------- man/figures/README-pressure-1.png | Bin 17097 -> 0 bytes man/getTtest.Rd | 99 +++++++++ man/getWilcoxon.Rd | 97 +++++++++ tests/testthat.R | 7 +- tests/testthat/test-example_test.R | 3 - tests/testthat/test-getTtest.R | 144 +++++++++++++ tests/testthat/test-getWilcoxon.R | 130 ++++++++++++ tests/testthat/test-utils.R | 149 +++++++++++++ vignettes/daa.Rmd | 134 +++++------- 18 files changed, 1406 insertions(+), 284 deletions(-) create mode 100644 R/AllGenerics.R create mode 100644 R/getTtest.R create mode 100644 R/getWilcoxon.R create mode 100644 R/utils.R delete mode 100644 man/figures/README-pressure-1.png create mode 100644 man/getTtest.Rd create mode 100644 man/getWilcoxon.Rd delete mode 100644 tests/testthat/test-example_test.R create mode 100644 tests/testthat/test-getTtest.R create mode 100644 tests/testthat/test-getWilcoxon.R create mode 100644 tests/testthat/test-utils.R diff --git a/.Rbuildignore b/.Rbuildignore index 40e9f50..4fc5fa0 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -4,3 +4,5 @@ ^README\.Rmd$ ^\.github$ ^codecov\.yml$ +^pkgdown$ +^README\.html$ diff --git a/DESCRIPTION b/DESCRIPTION index e961c29..f2c3927 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,17 +15,23 @@ BugReports: https://github.com/microbiome/daa/issues/ biocViews: Software, Microbiome Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 Depends: - R (>= 4.5.0) + R (>= 4.5.0), + mia Imports: - dplyr + dplyr, + methods, + purrr, + rstatix, + S4Vectors, + SingleCellExperiment, + SummarizedExperiment Suggests: BiocStyle, knitr, - RefManageR, + MultiAssayExperiment, rmarkdown, - sessioninfo, testthat Config/testthat/edition: 3 VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 6ae9268..e5bf65f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,2 +1,26 @@ # Generated by roxygen2: do not edit by hand +export(addTtest) +export(addWilcoxon) +export(getTtest) +export(getWilcoxon) +exportMethods(addTtest) +exportMethods(addWilcoxon) +exportMethods(getTtest) +exportMethods(getWilcoxon) +importFrom(S4Vectors,DataFrame) +importFrom(SingleCellExperiment,altExp) +importFrom(SummarizedExperiment,assay) +importFrom(SummarizedExperiment,colData) +importFrom(SummarizedExperiment,rowData) +importFrom(dplyr,across) +importFrom(dplyr,all_of) +importFrom(dplyr,arrange) +importFrom(dplyr,group_split) +importFrom(methods,callNextMethod) +importFrom(mia,meltSE) +importFrom(purrr,map_df) +importFrom(rstatix,adjust_pvalue) +importFrom(rstatix,t_test) +importFrom(rstatix,wilcox_test) +importFrom(stats,as.formula) diff --git a/R/AllGenerics.R b/R/AllGenerics.R new file mode 100644 index 0000000..1a48ae6 --- /dev/null +++ b/R/AllGenerics.R @@ -0,0 +1,25 @@ +# This file contains all the generics + +#' @rdname getWilcoxon +#' @export +setGeneric("getWilcoxon", signature = "x", function(x, ...) { + standardGeneric("getWilcoxon") +}) + +#' @rdname getWilcoxon +#' @export +setGeneric("addWilcoxon", signature = "x", function(x, ...) { + standardGeneric("addWilcoxon") +}) + +#' @rdname getTtest +#' @export +setGeneric("getTtest", signature = "x", function(x, ...) { + standardGeneric("getTtest") +}) + +#' @rdname getTtest +#' @export +setGeneric("addTtest", signature = "x", function(x, ...) { + standardGeneric("addTtest") +}) diff --git a/R/getTtest.R b/R/getTtest.R new file mode 100644 index 0000000..0455782 --- /dev/null +++ b/R/getTtest.R @@ -0,0 +1,149 @@ +#' @name +#' getTtest +#' +#' @title +#' Perform t-test +#' +#' @description +#' These functions perform t-test to compare values between two groups. +#' +#' @details +#' By default, Welch's t-test is used which does not assume equal variances. +#' Set \code{var.equal = TRUE} for Student's t-test. +#' +#' Specify exactly one of \code{assay.type}, \code{row.var}, or \code{col.var} +#' to define the values to test. +#' +#' @return +#' \code{getTtest} returns a \code{DataFrame} with test results. +#' \code{addTtest} returns \code{x} with results added to \code{rowData(x)}. +#' +#' @param x A \code{SummarizedExperiment} object. +#' +#' @param assay.type \code{Character scalar} or \code{NULL}. Specifies assay to +#' test. Tests are run per feature. (Default: \code{NULL}) +#' +#' @param row.var \code{Character scalar} or \code{NULL}. Specifies variable +#' from \code{rowData(x)} to test. (Default: \code{NULL}) +#' +#' @param col.var \code{Character scalar} or \code{NULL}. Specifies variable +#' from \code{colData(x)} to test. (Default: \code{NULL}) +#' +#' @param formula \code{formula}. A formula specifying the grouping variable, +#' e.g., \code{~ SampleType}. The RHS specifies the comparison groups. +#' For >2 levels, pairwise comparisons are performed. +#' +#' @param split.by \code{Character vector} or \code{NULL}. Columns to split by. +#' Tests are run separately for each combination. (Default: \code{NULL}) +#' +#' @param pair.by \code{Character scalar} or \code{NULL}. Column for pairing +#' samples in paired tests. (Default: \code{NULL}) +#' +#' @param features \code{Character vector} or \code{NULL}. Specific features +#' to test when using \code{assay.type}. (Default: \code{NULL}) +#' +#' @param var.equal \code{Logical scalar}. Assume equal variances? +#' (Default: \code{FALSE}) +#' +#' @param p.adjust.method \code{Character scalar}. Method for p-value +#' adjustment. (Default: \code{"fdr"}) +#' +#' @param name \code{Character scalar}. Column name prefix for results. +#' (Default: \code{"ttest"}) +#' +#' @param ... Additional arguments passed to \code{rstatix::t_test}. +#' +#' @examples +#' data(GlobalPatterns, package = "mia") +#' tse <- GlobalPatterns +#' tse <- tse[1:50, tse$SampleType %in% c("Feces", "Skin")] +#' +#' # Test assay values (per feature) +#' res <- getTtest(tse, assay.type = "counts", formula = ~SampleType) +#' +#' # Test colData variable (e.g., alpha diversity) +#' tse <- mia::addAlpha(tse, index = "shannon_diversity") +#' res <- getTtest(tse, col.var = "shannon_diversity", formula = ~SampleType) +#' +#' @seealso +#' \code{\link[rstatix:t_test]{rstatix::t_test}}, +#' \code{\link[daa:getWilcoxon]{getWilcoxon}} +#' +#' @export +NULL + +#' @rdname getTtest +#' @export +#' @importFrom SingleCellExperiment altExp +#' @importFrom methods callNextMethod +setMethod("getTtest", "SingleCellExperiment", function (x, ... ){ + x <- .check_and_get_altExp(x, ...) + res <- callNextMethod(x, ...) + return(res) +}) + +#' @rdname getTtest +#' @export +#' @importFrom SummarizedExperiment assay colData rowData +#' @importFrom rstatix t_test adjust_pvalue +#' @importFrom S4Vectors DataFrame +setMethod( + "getTtest", "SummarizedExperiment", + function (x, assay.type = NULL, row.var = NULL, col.var = NULL, + formula, split.by = NULL, pair.by = NULL, features = NULL, + var.equal = FALSE, p.adjust.method = "fdr", ... ){ + ############################# Input check ############################## + group <- .check_input( + x, assay.type, row.var, col.var, formula, + split.by, pair.by, features + ) + if( !.is_a_bool(var.equal) ){ + stop("'var.equal' must be TRUE or FALSE.", call. = FALSE) + } + ########################### Input check end ############################ + # Get y variable name (the column to test) + y <- c(assay.type, row.var, col.var) + # Get data based on source + df <- .get_data( + x, assay.type, row.var, col.var, group, + split.by, pair.by, features + ) + # Run tests + paired <- !is.null(pair.by) + res <- .run_ttest( + df, y, group, split.by, paired, var.equal, + p.adjust.method, features, ... + ) + return(res) + } +) + +#' @rdname getTtest +#' @export +setMethod( + "addTtest", "SummarizedExperiment", + function (x, name = "ttest", ... ){ + if( !.is_non_empty_string(name) ){ + stop("'name' must be a single character value.", call. = FALSE) + } + res <- getTtest(x, ...) + x <- .add_values_to_colData(x, list(res$p.adj), paste0(name, "_padj"), + MARGIN = 1 + ) + return(x) + } +) + +################################################################################ +# Internal function using .calc_daa engine +################################################################################ + +#' @importFrom rstatix t_test +.run_ttest <- function (df, y, group, split.by, paired, var.equal, + p.adjust.method, features = NULL, ... ){ + .calc_daa( + df = df, y = y, group = group, split.by = split.by, + paired = paired, FUN = rstatix::t_test, + p.adjust.method = p.adjust.method, features = features, var.equal = var.equal, ... + ) +} diff --git a/R/getWilcoxon.R b/R/getWilcoxon.R new file mode 100644 index 0000000..5dfae5e --- /dev/null +++ b/R/getWilcoxon.R @@ -0,0 +1,144 @@ +#' @name +#' getWilcoxon +#' +#' @title +#' Perform Wilcoxon rank-sum test +#' +#' @description +#' These functions perform Wilcoxon rank-sum test (Mann-Whitney U) to compare +#' values between two groups. +#' +#' @details +#' The Wilcoxon rank-sum test is a non-parametric test used to compare two +#' independent groups. It is suitable when data does not meet normality +#' assumptions. +#' +#' Specify exactly one of \code{assay.type}, \code{row.var}, or \code{col.var} +#' to define the values to test. +#' +#' @return +#' \code{getWilcoxon} returns a \code{DataFrame} with test results. +#' \code{addWilcoxon} returns \code{x} with results added to \code{rowData(x)}. +#' +#' @param x A \code{SummarizedExperiment} object. +#' +#' @param assay.type \code{Character scalar} or \code{NULL}. Specifies assay to +#' test. Tests are run per feature. (Default: \code{NULL}) +#' +#' @param row.var \code{Character scalar} or \code{NULL}. Specifies variable +#' from \code{rowData(x)} to test. (Default: \code{NULL}) +#' +#' @param col.var \code{Character scalar} or \code{NULL}. Specifies variable +#' from \code{colData(x)} to test. (Default: \code{NULL}) +#' +#' @param formula \code{formula}. A formula specifying the grouping variable, +#' e.g., \code{~ SampleType}. The RHS specifies the comparison groups. +#' For >2 levels, pairwise comparisons are performed. +#' +#' @param split.by \code{Character vector} or \code{NULL}. Columns to split by. +#' Tests are run separately for each combination. (Default: \code{NULL}) +#' +#' @param pair.by \code{Character scalar} or \code{NULL}. Column for pairing +#' samples in paired tests. (Default: \code{NULL}) +#' +#' @param features \code{Character vector} or \code{NULL}. Specific features +#' to test when using \code{assay.type}. (Default: \code{NULL}) +#' +#' @param p.adjust.method \code{Character scalar}. Method for p-value +#' adjustment. (Default: \code{"fdr"}) +#' +#' @param name \code{Character scalar}. Column name prefix for results. +#' (Default: \code{"wilcoxon"}) +#' +#' @param ... Additional arguments passed to \code{rstatix::wilcox_test}. +#' +#' @examples +#' data(GlobalPatterns, package = "mia") +#' tse <- GlobalPatterns +#' tse <- tse[1:50, tse$SampleType %in% c("Feces", "Skin")] +#' +#' # Test assay values (per feature) +#' res <- getWilcoxon(tse, assay.type = "counts", formula = ~SampleType) +#' +#' # Test colData variable (e.g., alpha diversity) +#' tse <- mia::addAlpha(tse, index = "shannon_diversity") +#' res <- getWilcoxon(tse, col.var = "shannon_diversity", formula = ~SampleType) +#' +#' @seealso +#' \code{\link[rstatix:wilcox_test]{rstatix::wilcox_test}}, +#' \code{\link[daa:getTtest]{getTtest}} +#' +#' @export +NULL + +#' @rdname getWilcoxon +#' @export +#' @importFrom SingleCellExperiment altExp +#' @importFrom methods callNextMethod +setMethod("getWilcoxon", "SingleCellExperiment", function(x, ...) { + x <- .check_and_get_altExp(x, ...) + res <- callNextMethod(x, ...) + return(res) +}) + +#' @rdname getWilcoxon +#' @export +#' @importFrom SummarizedExperiment assay colData rowData +#' @importFrom rstatix wilcox_test adjust_pvalue +#' @importFrom S4Vectors DataFrame +setMethod( + "getWilcoxon", "SummarizedExperiment", + function(x, assay.type = NULL, row.var = NULL, col.var = NULL, + formula, split.by = NULL, pair.by = NULL, features = NULL, + p.adjust.method = "fdr", ...) { + ############################# Input check ############################## + group <- .check_input( + x, assay.type, row.var, col.var, formula, + split.by, pair.by, features + ) + ########################### Input check end ############################ + # Get y variable name (the column to test) + y <- c(assay.type, row.var, col.var) + # Get data based on source + df <- .get_data( + x, assay.type, row.var, col.var, group, + split.by, pair.by, features + ) + # Run tests + paired <- !is.null(pair.by) + res <- .run_wilcoxon( + df, y, group, split.by, paired, + p.adjust.method, features, ... + ) + return(res) + } +) + +#' @rdname getWilcoxon +#' @export +setMethod( + "addWilcoxon", "SummarizedExperiment", + function(x, name = "wilcoxon", ...) { + if( !.is_non_empty_string(name)) { + stop("'name' must be a single character value.", call. = FALSE) + } + res <- getWilcoxon(x, ...) + x <- .add_values_to_colData(x, list(res$p.adj), paste0(name, "_padj"), + MARGIN = 1 + ) + return(x) + } +) + +################################################################################ +# Internal function using .calc_daa engine +################################################################################ + +#' @importFrom rstatix wilcox_test +.run_wilcoxon <- function(df, y, group, split.by, paired, p.adjust.method, features = NULL, ...) { + .calc_daa( + df = df, y = y, group = group, split.by = split.by, + paired = paired, FUN = rstatix::wilcox_test, + p.adjust.method = p.adjust.method, features = features, ... + ) +} diff --git a/R/utils.R b/R/utils.R new file mode 100644 index 0000000..661c970 --- /dev/null +++ b/R/utils.R @@ -0,0 +1,322 @@ +################################################################################ +# Internal methods imported from mia +################################################################################ + +.is_a_bool <- mia:::.is_a_bool +.is_non_empty_string <- mia:::.is_non_empty_string +.is_non_empty_character <- mia:::.is_non_empty_character +.check_assay_present <- mia:::.check_assay_present +.check_and_get_altExp <- mia:::.check_and_get_altExp +.add_values_to_colData <- mia:::.add_values_to_colData + +################################################################################ +# Helper function for metadata variable validation +################################################################################ + +# This function checks whether variable can be found from colData or rowData. +.check_metadata_variable <- function(tse, var, row = FALSE, col = FALSE, multiple = FALSE, + var.name = .get_name_in_parent(var)) { + if( !.is_a_bool(multiple)) { + stop("'multiple' must be TRUE or FALSE.", call. = FALSE) + } + # If the variable is not NULL + if( !is.null(var)) { + # It must be a string and found from colData/rowData + is_string <- ifelse(multiple, is.character(var), .is_non_empty_string(var)) + check_values <- c() + check_values <- c(check_values, if( col) colnames(colData(tse))) + check_values <- c(check_values, if( row) colnames(rowData(tse))) + var_found <- all(var %in% check_values) + if( !(is_string && var_found)) { + stop("'", var.name, "' must be ", ifelse(multiple, "", "a single "), + "character value from the following options: '", + paste0(check_values, collapse = "', '"), "'", + call. = FALSE + ) + } + } + return(NULL) +} + +.get_name_in_parent <- function(var) { + deparse(substitute(var)) +} + +#' Check inputs for statistical tests +#' @keywords internal +#' @noRd +#' @importFrom SummarizedExperiment colData rowData +.check_input <- function(x, assay.type, row.var, col.var, formula, + split.by, pair.by, features) { + # Either assay.type, row.var or col.var must be specified + if( sum(c(is.null(assay.type), is.null(row.var), is.null(col.var))) != 2L) { + stop("Please specify either 'assay.type', 'row.var', or 'col.var'.", + call. = FALSE + ) + } + # features cannot be specified if row.var or col.var is specified + if( is.null(assay.type) && !is.null(features)) { + stop("'features' can only be specified when 'assay.type' is specified.", + call. = FALSE + ) + } + # As features points to rownames, the TSE must have rownames + if( !is.null(features) && is.null(rownames(x))) { + stop("'x' must have rownames.", call. = FALSE) + } + if( !(is.null(features) || + (is.character(features) && all(features %in% rownames(x))))) { + stop("'features' must be NULL or character vector specifying rownames.", + call. = FALSE + ) + } + # If assay was specified, check that it is correct + if( !is.null(assay.type)) { + .check_assay_present(assay.type, x) + } + # Check row.var exists in rowData + if( !is.null(row.var)) { + if( !.is_non_empty_string(row.var)) { + stop("'row.var' must be a single character value.", call. = FALSE) + } + if( !row.var %in% colnames(rowData(x))) { + stop("'", row.var, "' not found in rowData(x).", call. = FALSE) + } + } + # Check col.var exists in colData + if( !is.null(col.var)) { + if( !.is_non_empty_string(col.var)) { + stop("'col.var' must be a single character value.", call. = FALSE) + } + if( !col.var %in% colnames(colData(x))) { + stop("'", col.var, "' not found in colData(x).", call. = FALSE) + } + } + # Check formula + group <- .get_rhs_var(formula) + # Check group variable exists in appropriate metadata + .check_metadata_var(x, group, assay.type, row.var, col.var) + # Check split.by variables + if( !is.null(split.by)) { + if( !is.character(split.by)) { + stop("'split.by' must be a character vector.", call. = FALSE) + } + for( var in split.by) { + .check_metadata_var(x, var, assay.type, row.var, col.var) + } + } + # Check pair.by variable + if( !is.null(pair.by)) { + if( !.is_non_empty_string(pair.by)) { + stop("'pair.by' must be a single character value.", call. = FALSE) + } + .check_metadata_var(x, pair.by, assay.type, row.var, col.var) + } + # Return group for use in caller + group +} + +#' Check metadata variable exists in appropriate location +#' @keywords internal +#' @noRd +#' @importFrom SummarizedExperiment colData rowData +.check_metadata_var <- function(x, var, assay.type, row.var, col.var) { + # When row.var is used, check rowData + # When assay.type or col.var is used, check colData + .check_metadata_variable( + tse = x, + var = var, + row = !is.null(row.var), + col = !is.null(assay.type) || !is.null(col.var) + ) + invisible(NULL) +} + +#' Extract RHS variable from formula +#' @keywords internal +#' @noRd +.get_rhs_var <- function(formula) { + if( !inherits(formula, "formula")) { + stop("'formula' must be a formula.", call. = FALSE) + } + # Get RHS of formula (handles both ~ group and value ~ group) + rhs <- as.character(formula)[length(as.character(formula))] + if( length(rhs) != 1 || rhs == "") { + stop("Formula must specify a grouping variable.", call. = FALSE) + } + rhs +} + +#' Check group has at least 2 levels +#' @keywords internal +#' @noRd +.check_group <- function(df, group) { + if( !group %in% names(df)) { + stop("'", group, "' not found in data.", call. = FALSE) + } + vals <- unique(df[[group]]) + vals <- vals[!is.na(vals)] + if( length(vals) < 2) { + stop("'group' must have at least 2 levels. Found ", length(vals), ".", + call. = FALSE + ) + } + invisible(NULL) +} + +#' Get data for testing based on source +#' @keywords internal +#' @noRd +#' @importFrom SummarizedExperiment colData rowData +#' @importFrom mia meltSE +.get_data <- function(x, assay.type, row.var, col.var, + group, split.by, pair.by, features) { + # Collect all variable names needed + all_vars <- c(group, split.by, pair.by) + all_vars <- unique(all_vars[!sapply(all_vars, is.null)]) + + # Get data based on source + if( !is.null(assay.type)) { + # Automatically detect which variables are in rowData vs colData + row_vars <- vapply(all_vars, function(v) { + v %in% colnames(rowData(x)) + }, logical(1L)) + col_vars <- all_vars[!row_vars] + row_vars <- all_vars[row_vars] + # Use meltSE with smart routing + df <- meltSE( + x, + assay.type = assay.type, + add.col = c(col.var, pair.by, col_vars), + add.row = c(row.var, row_vars) + ) + } + # If row.var was specified, get the data from rowData + if( !is.null(row.var)) { + df <- rowData(x)[, c(row.var, all_vars), drop = FALSE] + } + # If col.var was specified, get the data from colData + if( !is.null(col.var)) { + df <- colData(x)[, c(col.var, pair.by, all_vars), drop = FALSE] + } + + df <- as.data.frame(df) + attr(df, "pair.by") <- pair.by + .check_group(df, group) + + # Validate dependent variable is numeric + y_var <- if( !is.null(assay.type)) assay.type else if( !is.null(row.var)) row.var else col.var + if( !is.numeric(df[[y_var]])) { + stop("The dependent variable must be numeric.", call. = FALSE) + } + + df +} + + +#' daa dispatcher +#' @keywords internal +#' @noRd +#' @importFrom rstatix adjust_pvalue +#' @importFrom S4Vectors DataFrame +#' @importFrom dplyr across all_of group_split arrange +#' @importFrom purrr map_df +#' @importFrom stats as.formula +.calc_daa <- function(df, y, group, split.by, paired, FUN, + p.adjust.method = "fdr", features = NULL, ...) { + # Identify pairing variable + pair.by <- attr(df, "pair.by") + + # Combine grouping vars: FeatureID (from meltSE) + split.by + grouping_vars <- c(split.by) + if( "FeatureID" %in% names(df)) { + grouping_vars <- c("FeatureID", grouping_vars) + } + + # Build formula: y ~ group + formula <- as.formula(paste0(y, " ~ ", group)) + + # Run tests per group (FeatureID + split.by) + if( length(grouping_vars) > 0) { + res <- df |> + group_split(across(all_of(grouping_vars))) |> + map_df(function(sub_df) { + # Ensure correct alignment for paired tests + if( paired && !is.null(pair.by)) { + sub_df <- sub_df |> arrange(across(all_of(pair.by))) + } + + # Run test with error handling + test_res <- tryCatch( + { + FUN(sub_df, formula, paired = paired, ...) + }, + error = function(e) { + # Print warning to terminal + msg <- conditionMessage(e) + warning("Statistical test failed for a feature: ", msg, + call. = FALSE + ) + # Create empty/NA result with error message + out <- data.frame( + .y. = y, + group1 = NA_character_, + group2 = NA_character_, + n1 = NA_integer_, + n2 = NA_integer_, + statistic = NA_real_, + p = NA_real_, + warning = msg + ) + return(out) + } + ) + + # Attach grouping info + group_info <- sub_df[1, grouping_vars, drop = FALSE] + res_with_groups <- cbind(test_res, group_info) + return(res_with_groups) + }) + } else { + # Single test (no FeatureID/split.by) + if( paired && !is.null(pair.by)) { + df <- df |> arrange(across(all_of(pair.by))) + } + + res <- tryCatch( + { + FUN(df, formula, paired = paired, ...) + }, + error = function(e) { + msg <- conditionMessage(e) + warning("Statistical test failed: ", msg, call. = FALSE) + out <- data.frame( + .y. = y, + group1 = NA_character_, + group2 = NA_character_, + n1 = NA_integer_, + n2 = NA_integer_, + statistic = NA_real_, + p = NA_real_, + warning = msg + ) + return(out) + } + ) + } + + # P-value adjustment + if( "p" %in% colnames(res) && any(!is.na(res$p))) { + res <- res |> adjust_pvalue(method = p.adjust.method) + } else { + res$p.adj <- NA_real_ + } + + # Subset to relevant features AFTER statistical calculations + if( !is.null(features) && "FeatureID" %in% colnames(res)) { + res <- res[res$FeatureID %in% features, , drop = FALSE] + } + + res <- DataFrame(res) + return(res) +} diff --git a/README.Rmd b/README.Rmd index 5f5b812..6e5340d 100644 --- a/README.Rmd +++ b/README.Rmd @@ -13,31 +13,31 @@ knitr::opts_chunk$set( ) ``` -# daa +# daa - Differential Abundance Analysis -[![GitHub issues](https://img.shields.io/github/issues/microbiome/daa)](https://github.com/microbiome/daa/issues) -[![GitHub pulls](https://img.shields.io/github/issues-pr/microbiome/daa)](https://github.com/microbiome/daa/pulls) [![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) -[![Bioc release status](http://www.bioconductor.org/shields/build/release/bioc/daa.svg)](https://bioconductor.org/checkResults/release/bioc-LATEST/daa) -[![Bioc devel status](http://www.bioconductor.org/shields/build/devel/bioc/daa.svg)](https://bioconductor.org/checkResults/devel/bioc-LATEST/daa) -[![Bioc downloads rank](https://bioconductor.org/shields/downloads/release/daa.svg)](http://bioconductor.org/packages/stats/bioc/daa/) -[![Bioc support](https://bioconductor.org/shields/posts/daa.svg)](https://support.bioconductor.org/tag/daa) -[![Bioc history](https://bioconductor.org/shields/years-in-bioc/daa.svg)](https://bioconductor.org/packages/release/bioc/html/daa.html#since) -[![Bioc last commit](https://bioconductor.org/shields/lastcommit/devel/bioc/daa.svg)](http://bioconductor.org/checkResults/devel/bioc-LATEST/daa/) -[![Bioc dependencies](https://bioconductor.org/shields/dependencies/release/daa.svg)](https://bioconductor.org/packages/release/bioc/html/daa.html#since) -[![check-bioc](https://github.com/microbiome/daa/actions/workflows/check-bioc.yml/badge.svg)](https://github.com/microbiome/daa/actions/workflows/check-bioc.yml) +[![Bioc-release](http://bioconductor.org/shields/build/release/bioc/daa.svg)](http://bioconductor.org/packages/release/bioc/html/daa.html) [![Codecov test coverage](https://codecov.io/gh/microbiome/daa/graph/badge.svg)](https://app.codecov.io/gh/microbiome/daa) -The goal of `daa` is to provide method for differential abundance analysis -(DAA). +## Using the package -## Installation instructions +This package provides functions for differential abundance analysis (DAA) +of microbiome data. It works with `TreeSummarizedExperiment` and +`SummarizedExperiment` objects from the miaverse ecosystem. -Get the latest stable `R` release from [CRAN](http://cran.r-project.org/). Then -install `daa` from [Bioconductor](http://bioconductor.org/) using the following -code: +**Available methods:** + +- `getWilcoxon()` / `addWilcoxon()` - Wilcoxon rank-sum test +- `getTtest()` / `addTtest()` - t-test (Welch's or Student's) + +For more information, see the [function reference](https://microbiome.github.io/daa/reference/index.html) +and the [Orchestrating Microbiome Analysis (OMA)](https://microbiome.github.io/OMA) online book. + +## Installation + +### Bioc-release ```{r 'install', eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) { @@ -47,65 +47,24 @@ if (!requireNamespace("BiocManager", quietly = TRUE)) { BiocManager::install("daa") ``` -And the development version from [GitHub](https://github.com/microbiome/daa) -with: +### Bioc-devel ```{r 'install_dev', eval = FALSE} -BiocManager::install("microbiome/daa") -``` - -## Example - -This is a basic example which shows you how to solve a common problem: - -```{r example, eval = requireNamespace('daa')} -library("daa") -## basic example code -``` - -What is special about using `README.Rmd` instead of just `README.md`? You can include R chunks like so: - -```{r cars} -summary(cars) -``` - -You'll still need to render `README.Rmd` regularly, to keep `README.md` up-to-date. - -You can also embed plots, for example: +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager") +} -```{r pressure, echo = FALSE} -plot(pressure) +BiocManager::install(version = "devel") +BiocManager::install("daa") ``` -In that case, don't forget to commit and push the resulting figure files, so they display on GitHub! - -## Citation +## Contributing -Below is the citation output from using `citation('daa')` in R. Please -run this yourself to check for any updates on how to cite __daa__. - -```{r 'citation', eval = requireNamespace('daa')} -print(citation("daa"), bibtex = TRUE) -``` - -Please note that the `daa` was only made possible thanks to many other R and -bioinformatics software authors, which are cited either in the vignettes and/or -the paper(s) describing this package. +Contributions are welcome in the form of feedback, issues, and pull requests. +See [contributor guidelines](CONTRIBUTING.md). ## Code of Conduct -Please note that the `daa` project is released with a -[Contributor Code of Conduct](http://bioconductor.org/about/code-of-conduct/) +Please note that the daa project is released with a +[Contributor Code of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms. - -## Development tools - -* Continuous code testing is possible thanks to [GitHub actions](https://www.tidyverse.org/blog/2020/04/usethis-1-6-0/) through `r BiocStyle::CRANpkg('usethis')`, `r BiocStyle::CRANpkg('remotes')`, and `r BiocStyle::CRANpkg('rcmdcheck')` customized to use [Bioconductor's docker containers](https://www.bioconductor.org/help/docker/) and `r BiocStyle::Biocpkg('BiocCheck')`. -* Code coverage assessment is possible thanks to [codecov](https://codecov.io/gh) and `r BiocStyle::CRANpkg('covr')`. -* The [documentation website](http://microbiome.github.io/daa) is automatically updated thanks to `r BiocStyle::CRANpkg('pkgdown')`. -* The code is styled automatically thanks to `r BiocStyle::CRANpkg('styler')`. -* The documentation is formatted thanks to `r BiocStyle::CRANpkg('devtools')` and `r BiocStyle::CRANpkg('roxygen2')`. - -For more details, check the `dev` directory. - -This package was developed using `r BiocStyle::Biocpkg('biocthis')`. diff --git a/README.md b/README.md index 6bfef92..b4bc61f 100644 --- a/README.md +++ b/README.md @@ -1,43 +1,36 @@ -# daa +# daa - Differential Abundance Analysis -[![GitHub -issues](https://img.shields.io/github/issues/microbiome/daa)](https://github.com/microbiome/daa/issues) -[![GitHub -pulls](https://img.shields.io/github/issues-pr/microbiome/daa)](https://github.com/microbiome/daa/pulls) [![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental) -[![Bioc release -status](http://www.bioconductor.org/shields/build/release/bioc/daa.svg)](https://bioconductor.org/checkResults/release/bioc-LATEST/daa) -[![Bioc devel -status](http://www.bioconductor.org/shields/build/devel/bioc/daa.svg)](https://bioconductor.org/checkResults/devel/bioc-LATEST/daa) -[![Bioc downloads -rank](https://bioconductor.org/shields/downloads/release/daa.svg)](http://bioconductor.org/packages/stats/bioc/daa/) -[![Bioc -support](https://bioconductor.org/shields/posts/daa.svg)](https://support.bioconductor.org/tag/daa) -[![Bioc -history](https://bioconductor.org/shields/years-in-bioc/daa.svg)](https://bioconductor.org/packages/release/bioc/html/daa.html#since) -[![Bioc last -commit](https://bioconductor.org/shields/lastcommit/devel/bioc/daa.svg)](http://bioconductor.org/checkResults/devel/bioc-LATEST/daa/) -[![Bioc -dependencies](https://bioconductor.org/shields/dependencies/release/daa.svg)](https://bioconductor.org/packages/release/bioc/html/daa.html#since) -[![check-bioc](https://github.com/microbiome/daa/actions/workflows/check-bioc.yml/badge.svg)](https://github.com/microbiome/daa/actions/workflows/check-bioc.yml) +[![Bioc-release](http://bioconductor.org/shields/build/release/bioc/daa.svg)](http://bioconductor.org/packages/release/bioc/html/daa.html) [![Codecov test coverage](https://codecov.io/gh/microbiome/daa/graph/badge.svg)](https://app.codecov.io/gh/microbiome/daa) -The goal of `daa` is to provide method for differential abundance -analysis (DAA). +## Using the package -## Installation instructions +This package provides functions for differential abundance analysis +(DAA) of microbiome data. It works with `TreeSummarizedExperiment` and +`SummarizedExperiment` objects from the miaverse ecosystem. -Get the latest stable `R` release from -[CRAN](http://cran.r-project.org/). Then install `daa` from -[Bioconductor](http://bioconductor.org/) using the following code: +**Available methods:** + +- `getWilcoxon()` / `addWilcoxon()` - Wilcoxon rank-sum test +- `getTtest()` / `addTtest()` - t-test (Welch’s or Student’s) + +For more information, see the [function +reference](https://microbiome.github.io/daa/reference/index.html) and +the [Orchestrating Microbiome Analysis +(OMA)](https://microbiome.github.io/OMA) online book. + +## Installation + +### Bioc-release ``` r if (!requireNamespace("BiocManager", quietly = TRUE)) { @@ -47,103 +40,24 @@ if (!requireNamespace("BiocManager", quietly = TRUE)) { BiocManager::install("daa") ``` -And the development version from -[GitHub](https://github.com/microbiome/daa) with: - -``` r -BiocManager::install("microbiome/daa") -``` - -## Example - -This is a basic example which shows you how to solve a common problem: +### Bioc-devel ``` r -library("daa") -## basic example code -``` - -What is special about using `README.Rmd` instead of just `README.md`? -You can include R chunks like so: +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager") +} -``` r -summary(cars) -#> speed dist -#> Min. : 4.0 Min. : 2.00 -#> 1st Qu.:12.0 1st Qu.: 26.00 -#> Median :15.0 Median : 36.00 -#> Mean :15.4 Mean : 42.98 -#> 3rd Qu.:19.0 3rd Qu.: 56.00 -#> Max. :25.0 Max. :120.00 +BiocManager::install(version = "devel") +BiocManager::install("daa") ``` -You’ll still need to render `README.Rmd` regularly, to keep `README.md` -up-to-date. - -You can also embed plots, for example: - - - -In that case, don’t forget to commit and push the resulting figure -files, so they display on GitHub! - -## Citation - -Below is the citation output from using `citation('daa')` in R. Please -run this yourself to check for any updates on how to cite **daa**. - -``` r -print(citation("daa"), bibtex = TRUE) -#> To cite package 'daa' in publications use: -#> -#> Borman T, Lahti L, Muluh M (2025). _daa: Differential abundance -#> analysis_. R package version 0.99.0, -#> . -#> -#> A BibTeX entry for LaTeX users is -#> -#> @Manual{, -#> title = {daa: Differential abundance analysis}, -#> author = {Tuomas Borman and Leo Lahti and Muluh Muluh}, -#> year = {2025}, -#> note = {R package version 0.99.0}, -#> url = {https://github.com/microbiome/daa}, -#> } -``` +## Contributing -Please note that the `daa` was only made possible thanks to many other R -and bioinformatics software authors, which are cited either in the -vignettes and/or the paper(s) describing this package. +Contributions are welcome in the form of feedback, issues, and pull +requests. See [contributor guidelines](CONTRIBUTING.md). ## Code of Conduct -Please note that the `daa` project is released with a [Contributor Code -of Conduct](http://bioconductor.org/about/code-of-conduct/) By -contributing to this project, you agree to abide by its terms. - -## Development tools - -- Continuous code testing is possible thanks to [GitHub - actions](https://www.tidyverse.org/blog/2020/04/usethis-1-6-0/) - through *[usethis](https://CRAN.R-project.org/package=usethis)*, - *[remotes](https://CRAN.R-project.org/package=remotes)*, and - *[rcmdcheck](https://CRAN.R-project.org/package=rcmdcheck)* customized - to use [Bioconductor’s docker - containers](https://www.bioconductor.org/help/docker/) and - *[BiocCheck](https://bioconductor.org/packages/3.22/BiocCheck)*. -- Code coverage assessment is possible thanks to - [codecov](https://codecov.io/gh) and - *[covr](https://CRAN.R-project.org/package=covr)*. -- The [documentation website](http://microbiome.github.io/daa) is - automatically updated thanks to - *[pkgdown](https://CRAN.R-project.org/package=pkgdown)*. -- The code is styled automatically thanks to - *[styler](https://CRAN.R-project.org/package=styler)*. -- The documentation is formatted thanks to - *[devtools](https://CRAN.R-project.org/package=devtools)* and - *[roxygen2](https://CRAN.R-project.org/package=roxygen2)*. - -For more details, check the `dev` directory. - -This package was developed using -*[biocthis](https://bioconductor.org/packages/3.22/biocthis)*. +Please note that the daa project is released with a [Contributor Code of +Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). +By contributing to this project, you agree to abide by its terms. diff --git a/man/figures/README-pressure-1.png b/man/figures/README-pressure-1.png deleted file mode 100644 index e208311082fc7bbdf9503283d73351b89a7f71fe..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 17097 zcmeIa2T+tz7bV!3P(a!Mq97ompW%)oBO=Z*H@r2P z?Wl<#RX;s8$}7uN70~ozR;87fQkk4zX(sxWr-7=SDIuqL2WJQ(9WK%yoS8w92ws6? z-ce(UQ-hBLI&CNS^bcdwo^YYX3G>pI-i-0$D>b-{I(sr!p%NpJ zVJ+&kSBk;+OqSKBh0Kc4J$sg8+8N)|)8n!};^O3_TkDvskRMZ5@U$;b#9iXfwO?&L z_48jIQA)qj$A>@y><`w4BPfM@&rg=6BOY^daS=eF)V$WCf%wp+?i6|#i=bR~Ji{CJ zryvleH+$eCUYqfsAFm<|J-xgzleL3QO(wew?T3el10x!5JnwFs^7h|`K!O~$CTq13 zW{cJij*hzr2e)tCvY4u?d@V^DPL&#Ho<#xiPgiD&(<$cF%$CYlX1cwUS-$VJoy*_f z4ZOndbmz}$ywvO0&42zRG`xfOZ}8ZR@1qU<(nNi*Veth81wVfLz#gp2X{krOq23Nk zURm<|PFE+RsHphnO`Dy?%a^~4b*n=)FYU3I9Vu!#kP{LT(*ETxZBbECKw#jXa0)&j ztTQunV3N8z)}OCrswFc`_eveE zG1cIKcyg0GymeFyI#>Q>|NGoVi0zGynK~Ox$TurD%+yykD|xl5p15S9wHMjuC^UHu zXYXx@IMe6EtpML^4No{OZJH}`d({xGsQU%L?7Z|@ICVV16Ap>>S!P6*kols26 z@6R{7V*uEe1NaR~Yv+NHA&*oCt<+51*Q7EZ{ukh zsER|%I&X{;BjFWbz%p|ph(FQYZ3v|84j%OX<8L3y>HmB&UW4+T+1;c5l~6L}llm>S zgi~XTUwyZLVzgx0zmXMQki!Oiu_+7cwQbzh$S99I^+KY~n_mUNHP2d$b*hQ>i)|Ae z_97)|belE?>)w_U}W@Gz>~x}-n8QrAp@FVMLHQJatP$JEY@FBOGh_! zr_5p8JKkI*L{3FVD>2qW`9BbDusUV#KT~<)bEE5G(fWP#(&81)7w2Z%M}uCCl}yM# zDH8Dm>Ng?&zxqL1+>Zo9girmi^|3d!_Hf$qt;GOYjGaC zKV@y}e*Au*p6c*oXm|`^p0AjZJ~ObtzYpHT|Ih(*3E0QGaLl4W3x_fOe$7{Dg+AYI zf=ViPIjx4=r?Yk1B!@)=J z0k2cFf)HC47M6g3fL;A46_}F!A4-6<5?KFBV28-2QC$=MmoLWmpUz&Fc=zwaDA7rL z_M0!q)6SUlk=0x5;ZJS`@M3$P(+ctJ7i#q*B!X%7&cDpE0SXu@1qDDSD4HE-nD2%X zpp1Eb^CfmQQph!Rx3h`B<_6J*3iIA*16qgZWc8!M8WV`+X9Lt(=&1{bk%a0iOBvIm zExkH4&uk3qcH=Kj5yA9cDK`V2^gQP0{V}mYO6KK9{=}yk?EcNwquI(M+v&T>M31+u zitk*LprCZn+#D>fUq;}jNe_#sx8{A>7aWbJ`Kqr%{DlFgcQ8>6ll(PjdlnM#}*WzIz_tLSnc_;^pA z3K}2Mr<#g&UPm})+c)yP*$A&r57XL;0k?7*xH#msyw+?f*e zuytV4h{1eC%_6NG$vVP-Of@EY#q@@EzXC4Yf3E+c@dR7T%*NK!-OZ|La1ku#N#Rf1_0DaY@&f%No5`zB+)e*{5{ZhNs4TG5)~E!*Z`|Pr6{WRyi5$l} zGxPRZn1o%$e>u;sJV9SsJh-2x4|Orsq6Vq&^XU}>H_Su`hPw3p*T)oPOW{Yq%;Y;CSy$gYmsGn>0 z>keg1g?qNFm!DyyISuhr?oqZNDXKN7hwqw1Fi$Ikuay+icsn-D*ZeQBBy1m7glJXeU zl;pv}Re0RIgW+V07vByKdHKu;qjWz9T_#n@)aVEeSP^r>I zw=V`3iz;JlcC6}?$ThwrKD#A8qiO#wtMQ~W+&}9wjaMLAhu4Oy+ST9}xyIKd?jQ!D zV^jxtTA8dUC`e=Gu0WJnPI@IXf{XA`b^@sjTk6Q6K)(6u?~9m2i#7G0@!uUI{gtWk z?6llmo98pccD-Yc0yL#oXJBdPyJ=y*gR{KU8Q2L)*SfFtP%6<%>*`&X(g7q?^lZSK zaza7Kyg#$kiq><;(G5|}#WV2!S&K*wrKck%^w&o*s&Q}YZmY%6ecfvi|BH%KwsOS& zmaub8?+ksNa^k2O!O9Rqa#$r?hZCEhv{y28v9W(B|o}Fzx@*lY#GO z0h~H0cIKBWGm4K67bkN)S@*k9avoy3hj!|^gwiNJF$Pn=AwH_g*%#u40hu#3z_nG? zj&gc`Z&geAPQ}qh;#jFr>`S#;N5{=ODe~i= zG8kG+fiGZCLf-d&@9>M;G;_-*O@36>+DgRQpm}Nf*>`H`u2Q@p3!C21Ihy@kIvhIg@3h?sicY}S0R~gyY->d^{pb%fie4(~enTC~ zG&hsbY=(oQL)e6#VCJsJg(HVzUAOoiCqyLV3(^%`&%c?tZ9%bj)T&(fa}H=jBa zGgH+0f%)!&Y*ZDaM0dfs-SXrdvQo08BM;YxnNhhg7;I-}=XoJ1$S6^EYUT!~+-<0- zKK6*_wdOx5t+(~;Rkn+wT-9 z{#&=8NDVYS#1^274c3ABlGOb4`>X*+dhD6KefzdMnMYqwFNaOKYFFsB-CgbtR@|XU zGJf#`c1m=GL7WfuX$~}AuatY&um);y=M$*;)~;z3$DGN2t@oZ?k+jLbLHflV@-1w3 zcD6rT+Nx4VQxgpe+nzKbcJ1=!>1l(ty7vq-w449z$y*D{j?=`{iWO5{+EWN8dYmD9 zCVJdp@e)OKhz*$iTCPO~$t3y|9DKAnaTLbhSmU}mUS{&=DJ^YraWU>-o!&)%_t}}x zc0i*2hrc;?rIwC7>sh@m=6R(LMeVE0-fg^C=!6Drco?%l9+m?MK|Yb=J75xH?5r#- z@of5Qo0~p=e_!|W^Ruvu$+!ZoTEEA3QDrpv3-#1Ee8Rh-+~jw9TgRr3TKAss!9AtK zL{mMS$zIYm_1wD3PIx+Bn2aKivDto_E;;S5<}FxNESdnqV5M1V;4>f2=UBxM*O4cm z^7-@Ud)v8e7yi9Bb*DNn>s0XQ_i(UioIb96jsC5BmH4` zh~NmFy_}^9W-{mzLiJ|AH<$O;l_il&kIxsw2hJt&Am2U(1xc!?sFbt2JYL_6cA;qP68xsKW}2wJie3Fkh}!q#LRI=j z2K}L07BzrEIJ`F}s%L%}z3RDCBEaUwv+ByKM}J);+qTFjPOtKF>iNpMC~uH@!pnv} z&kiZR7Z~~ZjQIr4f?jLdBqMGxKbW%ewGkV25TeDvHX)E3misIHN0=GRFkI8#e!Dw` z&+qte@B>yI@BM!1NtPqKHq0F)-3wK7Ys+L@IdAU3n?B3R3!Q!PF}~xONORpY4%5Dg zkL#q5o}yME;-YFy>R=R2LmwcTR3s)ow6ntiy@S_oT5qjxff40jZ84eWyZ`&T)y(uh zD!F%}d10?pgBtOHzevLI=^c44kuAAe2GXZM=`VdeoIuV67YM}pCc2t1Nlq?Av7FjnBr(NK_qfSgLOf-XCWtLO%NT~?Qp_<3l?Rb#P zZ`#H@caz9sW zf{gD-uUQUScv9rMPX_FT?<#JOg1FWk2G)CYC;` z;p1hP_jx|F3DA?sRV>SQnD7HJ2^4izq0M|)pjhu0$qi-nUf zc(*DyL2(lx6X?d*V{M}TdH4#%|7APosgQoWj&iL!iqDX{rdDHr{{nYmVr<%x>k5X= zJ~b;xK(Y3(l*Z!0e#sIhbmIMQq|g6o`Fj=@JP73M0I>Sb1a@6*?criwbVrBmnxoP+ zaKRI(Wu<3FETe9f_0Ha&oQzDZ({fLP$KJFr&IY#r(g;X3{=#n)IsUfC(DT_&{y-q^ zU!xHgKCMI>AgZe|V~G5r-0w>K_)th2D<`KzXFThS_ny2fC~RgoH<9n#qb0*0$ml*K zxa?glxUP?|LZMI!ij>rF5N*VSgfWHBNJxPY)Yje6!NkeQscGWv-B4Ls`RpL`ArSb* zZ4C?z1k<**|27Fwpw6-D0QTk7?}dhjZ2(<^ktDPfbzbLpaeklf?MmPGS3pH=4W|eU z3{*_zr*W}hVDTYvx9Mo2FhXNQ72YVIAHuCunLLh%3b#JcxoAVS4M|8yoU!|>L4;%x z&|CbhM0Hzc%|~ym(7eNt`(k zN6RADQsnd*0It2l`lF}DD&E(he4(V*tUz?vrflgnmP=2yQsOrZAG^$iFh{&S&%1vM z3K6G~@gqg@VkIhm3;p%qB%f-pcyc+P$IwBTZ|Z@)IrK-_1WFe`5${ z`1uvJ<6y<%@f3~kj8CO!=d@vhO^9q;Ur4}9BV+XvolPm5d^-xhIxl&v!x`|1D?XPZ z&D%~s>axs}FV_Q}tdw&N4C&l=B79+@0*KmNe~drbwU7T80}uA@e}ozsu-$(MPs%=? zVd(uSVUOKsXhfvrzz1&f@H{ok8xZkcpm$!%*gvyLlsPW)oA#m9qLWoNDGW;@^bu0o zK-p$Dss!crHPCVxMfg$<9;|DoVRnrCyP5BvTW6(Ift&Lb(gelM_Eml5K)werISm%D z!vChqh%V6mQ)N2#X=blqlC@0dRhgfJ#vRD8{2o1@oPs-jdFrgekIA~$snVuAo)#-4 zXbh_l#Nk8yjiSyi90Q)vrJq)$Nr&Iq#NDQvtk$PI!qQexs-UzAYile=b)Amrw`;T@ z5CVB1-}vmt&>(s=sSvBb&uZ`!j=YC8%S=B+TKw7d_5{;nAOf3Z_4{zE-q5Y%P||LQ z>%k`i%ITBw38s=?1sC1xPQFGa5gk_{N{*d&Fkc;o|~9ukJ{+==9Q;i5%=XL zm!GJz=}i?t*iO}LZ*N~Gd{6|miG{W(74-EG8NaAB0I}tIogr{t1@GL-VK0gR`eIN@ zZK=wcn<4T18?2`r#+sWK5)6Hn8At-20BvFo;DD&8C>fjH_Ro!}cO=;7?e$KR7> zvfY*nbSh{(2zIqCiWYQGdKHeWnc@n~d(%NvC+Y{&8-51VkKGl zkXUg5Bu5)#G2g!3AthBudb-HT$-$Td_5l=ryXE!o1%vxH&fXFEzb?2m>+M67@`fxA zKyOGfL$p|pDTUjD%YFOZ{rmS}VF5&C0LxA3nV1Sej|M121Ox<@(+ww}`*nJDmYtoQ zTS9oL`iPG`5^~-6M)Ru5W`g%aU*`8$!CT{%a4U2P26>b|1PwGFUIQ0!a`ump>*T5|T2~+tsthEcI#E%-Ui)+wFY_zE zfIx0;ZcA(HlJ(cGU)5l0BO@dK^~=P>L?UwT52wdFj~+c5V+Hrk1MufG{!|ep5f2Q; zAOyUVn~Q7r@DNCa`(VqFNF*plM9+^+VCw-y4pEelSjsrrgZW6Y&!0Y}3OKd4wSD;g z=@#%F6|(Cdmoq3j!Ct;(^{>B=c(wn3ZdzPC;pHV1Ke_(kMG~` z{P}a_(X>}pRTUTz{Ck0$xtN=u&t)TEU;z%bn>lPL&VQ^ z1nviuyGjrQ1`$35Hp^lt(XH+o8`HjX4Ie)k)Esflnr!H7V+c|b>;qL#y2wPWlX=j~ zC`1`BKaet`3ZG@Z0TyTreEGd^)w!13MNq$_4S)XuNVFaw26Ea#7IiW`J=~lCeM%;E z8SyBfE`AFJ$T&7Oc7T}?d;0X{olm$U4D=D#Ggj7OyO~Cngc~FTm$xf1kbqh5A0FC{ z6eWR@L0sG)4qrb%!-a=?e8MNXs>=jP`ttRwun)FErqz3whMu+Lzsa|&`QJ8 zvraab$s<@?+Zama|LeY_!`cw+Rg;6?|z(^t@BA#O&b0}RC_U8iX^8`>($1ANs;2>WG6XW6G0riyGC_}{S zP_N>Ha!>rhIvRv&eN9d55A&b)4k?fQ%d-bp>)bHYfK(L}FgQ6mN!$T3GXkdf!n{Wf zrjk24<%?a-mynPE)AQQv=98C~zc`+c^uZl&&Y)2H1o%XtHv}Pwj*4t5A$n|Zp2D@GjhY*l_RtEEX0WfJ6>!?QqeT}j=O$ZJ87vQ2(Fcy#H&`m@=+meDttqemEXnxz_?E5hwP845Lfr=q@t~Jos#TWt zDK|MW@hdqPr$G7n2L89(hqHJ|U&6!FUYIMOe9xpJsgS_;qTSp94FgOO^t|^bYn||k zAJ+i0TUxGx=#v}Bm(7(`Q0N&>jE#M7WTYHzi`Nk9e_;;z9dJbe*bGx%miO=70}f3A z&@a@-k~>9-T;1Q_4~i@A<833x!=&V7PtZIE5gt3|Y(gHsG}jUa{IbFG zz|!J?@#gQpO~GOpr+c7Q3=Rs4lMM#BD#P#G1GMIm4Go_AD?7^>e(?>xO~Ft!OqE_# zv=R6?c!3?+0QYEUXvkpo8r*kIPxj7EPsZp)z+ zu*ByZsw6F7XVijS^0e91s7L2@L>f71^Hxd2NxG`>|_&^$#SFl z3i-IXF!t26Gz#tyvlt9!Zmt{z$h||x&^fSKKHlD7ztJE$f&Hzp9xGqSm6DdOwjLYk z?X?4?0T7N#l)nE41%%vjX>xr0NTx&}h#q`=d_%DHV9KD!Z`k1OAcIwd50f%$2sgBEJB(lAk#x1^D9N$ns+BJE-qGrqZbp^cGNyvJ?KcY-pzfy5Ji}xOe3qO`8 zT^h!f5Hefjgc_G!IR_xI*qPwBgM36uIR(~pmsTn^hbT<}4Pp;8HS+WF%uuI)gQ@lO z^|Q0GWCu1uAOerxnG1gT`t@soP)l7&s?nl@gM;8~$#QQ-#@Qmg{{xV9baG|mE}Mio z6YP!0PdPX^o<4nwTLdKx=)0hcw81@KJUu-<*R9EHbA;Y+peVN)y`pN9Fqo)_*=Z{^eKVIq22K(R>^=K^bZ4wta zsbPG?7@CCw%;XH`RP|a?a%^G(KD=gpbkYPBbKaU1&@=&Gf##L+hG?ZuF~EL+cuu?q zHTK_iDlJ>KfPaF{DyY9eDF*@8`%qQo6guZ#mqTs`e+b-eu-oFD>JRr|-gUm#+g&Qc(KihykY9vD~o z8W}!FP~!0*=HMBdjD?YMqFrn#DKo%TK2$EC1Tc^uqbyZnYNfuLc{w=^H9^%6ySlp# z&(fk&Q}uo?D<`LVu`HX#;8uHD>GgT`aL@EU9SbKMcuwqm60^y_bBYVuk=;4 zEET4O;CN-_gPYqCkgx%VeJv&Bi`%RwWm3E3+HzX$6&8b_$L+c?x-}rpW(%DIi<_F9H0w&d`{02qpfP}1*4ETK7Iy#fDVO*WR}N?(gZzZ_yEq5g6fpS9 z4*9)94}bw=WC@C*C(pRKxt~3I-}?1YNl6K4NHfZ0fK(pNHx0^xUjSLz^lCuz`V1Zk zg1Y`wE{KQv`ua;12M}J}I%f-UadCimg9EajCVxP7)`n}aDk}B@Nf+?`+f>5*AT56a z5Ku6<3aStA{Bl;^Dq-*A9}ARM7zA!zHJ}O28p6L@ZgKB_d}_q=VA9)x+ji1@XRZaX zOi3)fxw#q80KlxoaCa8>R|hG0t+jTIK^LGBI3g%PL9*T7oNxnW#d87iv`yf~2=m^a zbb<=76hIGHPWsi9N%(+h3~mA0wzj?=N+o>B0C>5WR`+iKFbCT%0O+p@%zC_HRhrso zVQER@x9D|=JykHfbWxwRzkmNOE*{vGt^h7^amWPD zIuRTLW2FH;HhfE-@UkicxwI5eU6|(c@1T&Aj-&!+lhc1=0pOw9Ziek}ALqFa5)%~_ z6<6nPKng&a$cR4e_{_$|HPGG;8&y>AQv_m91$b3f#?^frAfce}0_i17G&AR+=h2pa zrRA`n7pn&rI)!g|uO5If_yJW;E?W*a;I>4Rf)&pcBl=%~)dK_`i>!#cTmgWCHfyqu z!9k~^tts$qu4KQoCQ+@B%|)qu+@ zC{!U3F?HWm4uP3L9Rhbi zOKaNjWS^}^%ePO`2Zf{=)E$&T_e4UQeDBoVMxes7|iKlxeI*he}Y?Rs^liOuLz2V9VvNl}RM^ z2Tn>oy|e+6W-5&jk^kW3mY=gEl(^s6U>C`L(a|v$?_e{pWawoQ^n~FnOQ8$yR4uFI zPT2k%szBNC$<8^uUp5O|-lc+GN$F4~W=C9e-dM(C38VM(gBS1@w|WU!-o*^=6?0xl zJlc|en_%x~r0GT|$vf!SJv~91>(SWtZiaC)8%bT|{KAKqNj)yYDY8Xm;-6}UpI#7Od|tHP2u4! z|{fY(I5+slN9!@`4C;IC=18)ETUo^I-Hf^#1zKBRBm>kH?3d*Tm!o z>!{Ouztv3pWrmIM6eU}Px92;0A9Jz$-Fs(iV6~9fC(T)X=rh?dU9-00qj+HvMb&Pl zlnmcH?Y4@TayVGTxvndRDv%ubDvb+uPG}{mZ6fx!yN_dQ`fTn~j6c+TSpLF*{y|Oq z>d}Z+UnbEZZM>_qx6A#6DkrGaGq@^-+Hk0)d-OXxvR_b{d)7_%==}tQET#h)pXQmE zr%;E$WShCR#z}3--IY`(8T6%iXVve}o@>QEQ|Psm1hlof8pCN8^NmTP^QbSQ_+)%U zg+`Xcj;r$WAK?Rb*SsaTjWTW|j_)7;9a9vW%C4kj?b7viLC?C*-RE_!wtrkF6xnMy zU}qE z`s7Fc&T_)bsiEYjB$j4IB4ka(P6$@#n;;GI6lK`ll$3vKMS3WZs6{A3Uk8hEV`k+qGGnOF1Cf$A@hAxmi_)SjGQ*h;`oR>!)cFm`Lr^( z8hFSuE+6q)d1%^TNbO()b)j%%Z-Y;WFbm7L(EiWd(;aq(F9zQ|Twz=+rx69^XlRyE zFfZ|#v% z<)|@{nk|`uPTP2!P}IxewYUMQ&xQ@oXAWzf(J$E3`Ac_fXR3A)VhGEz@37VGn0ynW zqo}oO-gp(tca?QIL*pK&`(=htGHG%s&xE-b(5E;US!_+!CLDYiq!?y)KRLmvRk@YC zu-{x%1IPa5)tY_2d&d#5i%!!mGS}WZ+L|1U$!zPMJSvU^3`xk#nQ$3}@2&Pac0~mi zaI$y9^rVB7q|6PI`{gfOE39S+pi%6t_4sH(CjUu4{9=PJ+% z*VSPQ^M)uh58UoY6TikZ#+w0bhqmv_vy2xBdc}eW;J`S|UU#uvg~m^237_QU<}g}s zSR0ZQ>U~N=RNBdwp5vZC$y8GdcwK_(l({jZd<6}2r?Kcj6GVbFULW0Ud%D%yR9AcsidX&Os zb5*D@Z#puRPpbXGLu)I&bd*S~;`E1?ve7KmI)AtFZU+yN(9yLLa4z4Gbdma;sl;^K zy2*X9_}=g%OXcstVXw19=Z*YEc-!c|c3$Z3HQ1|A)z-Y}ZnesH-%4!kyt_xMU+G4i zs}*iP!iAh@!dAe#^PF@_kK`p@O6qs(@KR#1X3m&a!r$qlA;K4f48uo0<8MWyp9Z_h za;%r1_YWpc(2pFvJ#KJXKB3VLc(<{4-j5$yV(mfS*FEsNA8{6Yka19Px;}yFHDmlL z|Mhyf&e^VX=;h+1u!rB72t>TZ_QRsG5rN-gzs-M^IT&nq4hoD3{D}P98c)nd9M6s+ zDzlu-YLwYUw(+Oc#>B*ZOZ;}c8;Irp3wwuIc=wLtb5wRr_BJ<~|Mh)*dRMaTXaj$Ft^h-9An0Fne!pyQx>)-dDPGz)3h(!Ky}kGjBwJ z?U0?GTbZBNX3J!M+g2?v3d}r;jgw||vEbhQbn@9FA7bpG?+R%)foSE%3wA>yFJA&W z=KYriRjHAYVG?JEkR5YBy>jE;%>?@B=GX<2v=C%qc$7o8AWx7k8N1j_UJJcn(BzBV z64e(+dSspo+Sw7+7S3F2f%n_!fWm+WEF3Ea}7^i+6N3v7iw};W1>?_ez=wM z4d5g~cD(4FsdmE_6ZL8SpuZA0`#$mkz=v(*SrWPgK8u%TCpvYC{P71_k``9xv^n{? zVK+_{OsPEcqfLL(YZ`2YrHqbi6RcWTWNrQGcadupINV+|NJ&gheLBi7n0Qm8a;+vz zliff=Og9U8-|++{SuD7P*HyGmuVVhJZm^^;@;Rhdr>tQ4+2m6T9UV34C}@JMYWXtna~hEPfE-iFj< zR|@m^MAG%H`mBDOghdJKhxc1U>%r0GEfEP4EJ_CL^_WOK3{yhGzthy$OQ8 zeyhE$UZ`OxNupQl3za=`_aU>|!%`c1`=aLtY3D`cH?Y>eS%WQSqn#r;OU>De8B@4a zd;6k%*t#z2`O0f$AIs2uixvq~3wTd5$kw&}!fB)axq{-B2rKcG#&^jnR)Rt)@phud zt@ML8c$_)GxticNB4b~|#GTetwrpD;43871aJlTZu~k}!>kF~eM!(z;?7ui%JJQto#47PMU96Y!ZP zCwv_H*SC4dNG|Yk>w`*ozqYMQJ1+lv#seZ`UNn3>X3 z*WO$lxioZUNhJ8TwfauWP-k7u3yMaq+PzXOQ;!P+lNqnUUWJ)n0gKX73f`>`v(XqS zbO)z)MIs4yFu~A9gd%QHEx)s%y(y)RiLmUA-GTsv{%LN8X7PlB;M49_q{4j3tIw{5 zMkn4$=ftvJxC!O##`-i|wB_J~7j8YaHgqI_RBnzOEB^By>?YE%bljU|!y#F_cDxuH zPCnf--y&u91rrPU#!u)B&ukWx4Hc)pc2^XvO!&?w(iAg^_=$0CQ5~Bt=DpH1`PwnO zZ+C|M!AD-#z)z6ja#7bugaV1dvUx|0arSc)@$v+0YJIEu@ z9E5xWeWm}$-g4YAemNS2%-_nB)0h7HP@~fS@!f6zkA0ttDW=y0DLDl}MNiUH_~0Au z6{>HHX>@Ro)xtr`mO1BA898;bWlTrOt64G7{32?y;_ZGgBLhv>8gBS?YG3?>F4)$) zT9uN|D$Yd!!%G{V_N7z7X|wYB?P8stX>PcUB=d{=PWX|BsuwlXUUP_P7d=+i4JICU zr^OZ$i97R+tMApl=Fk~OB^GvO#+KS8uFFSNdR`8~A~QXuzA7_?n;j_Pp7nV{9(#o^ zQg0hi&+7Gkn8#pmyRDgr)=#(cww)#JMEkaBZ!XTZ^ROhd*Wyl>KQDFR(jUl^O^qz} z8#TgjD7}9QdG*cRZ@GYJ>7!Gxh-Y+R|D-m$Wo)RpJ7*ZRHGSS)!KU6-bj2)9j!aVd zz9{M=EF{S)CWL}Goxj|n2o*OlUYlx;7SN5PoH;WsWolhwO_2KFQnoIboRXBHJK||F z-D>D}VrCY(w$!TrdF#t?NS2oPvgQ&t34QZs(TokYZu1vk1|8*N)#1_aJEQ8TO|csD zz2&Qx+KGHRcS`QMtO~&8#dZu_)DbOr$IWTm1w~q34KIZ0j%x~ab+`>fLPLS^UoY*r z!C_#HaMG6jCX5>AL83kP$ynCcJ?dee%JrkY!3}>tPCKM9mykk?-$j$7L^0aYs|S8(h{Ms(>VTR;8h6ve!R?je=9$rIub*j%v2yPUYGXS?<8~&) z*IYJ;xbmuTwh1?M1i4hjW)6?iH6B3h4FgIRWfD>rBGZ%5%kWlDkJ7M$d7~o%(KqR} ztF{MCqo|E_@@+K93=_qkG&6SkqDWZz)pxdGRwl|{^*<+2p3m1nNbtP2{eLNI>1Zbn zr5tFE%n<*vpf0hs()nc2T$H$SOjbW(?D*EYE-E~&ayLUz^PvcC{EkWWY3b1)GNkCV zAK%XEye=8kN2tNo$n}x4jAD!v-C-E9oVgyP97NtP&K-+^pSNjBcfel{FuDA@0K2W{ zKGo^ewIfrwEmo8iW>b+>g)b?;`pgxhJz$ABt=UN|QJ9CU|8^^>_1GyLvP9Ao+vxQV zyt(tKONftlPdL6riNQH!=T@&qF0p{KFkA2ZaV|F-B;biX*r~@-vb3A?KM;OX9Plxh zUH1$Pgu!~k>3QY@b!Sb#R3WE&4p%=_If=dp-^DiDBkOLht0n0s z>{)Uygf&6SOHI+AdhOnA9sEi%53U`?Nyoi_>KrH~OtoYWZ+7WRZ>X*tqg}Y&9-nTW zpCB#kb#(YO?D}lxjx<%$OPK(AbRE1X{IwWWN4954E^iE-O?sN~YPNw7X{=O>U`YZ@`J;OLxqQHVzp_mX9orxZez_V7r4 zS(duy3~TR`&n1jxHKS|EP39ki@763p#VzLVi1>45aa1}iISOPJ*@M#}5Uyq#&so%) z1ac=aiB-LR<>xNXA?MHk9kT@A{NMR!8P1^`m3L;?@!i3%lR+Rdk_xX&UcL+XFN@__ A2 levels, pairwise comparisons are performed.} + +\item{split.by}{\code{Character vector} or \code{NULL}. Columns to split by. +Tests are run separately for each combination. (Default: \code{NULL})} + +\item{pair.by}{\code{Character scalar} or \code{NULL}. Column for pairing +samples in paired tests. (Default: \code{NULL})} + +\item{features}{\code{Character vector} or \code{NULL}. Specific features +to test when using \code{assay.type}. (Default: \code{NULL})} + +\item{var.equal}{\code{Logical scalar}. Assume equal variances? +(Default: \code{FALSE})} + +\item{p.adjust.method}{\code{Character scalar}. Method for p-value +adjustment. (Default: \code{"fdr"})} + +\item{name}{\code{Character scalar}. Column name prefix for results. +(Default: \code{"ttest"})} +} +\value{ +\code{getTtest} returns a \code{DataFrame} with test results. +\code{addTtest} returns \code{x} with results added to \code{rowData(x)}. +} +\description{ +These functions perform t-test to compare values between two groups. +} +\details{ +By default, Welch's t-test is used which does not assume equal variances. +Set \code{var.equal = TRUE} for Student's t-test. + +Specify exactly one of \code{assay.type}, \code{row.var}, or \code{col.var} +to define the values to test. +} +\examples{ +data(GlobalPatterns, package = "mia") +tse <- GlobalPatterns +tse <- tse[1:50, tse$SampleType \%in\% c("Feces", "Skin")] + +# Test assay values (per feature) +res <- getTtest(tse, assay.type = "counts", formula = ~SampleType) + +# Test colData variable (e.g., alpha diversity) +tse <- mia::addAlpha(tse, index = "shannon_diversity") +res <- getTtest(tse, col.var = "shannon_diversity", formula = ~SampleType) + +} +\seealso{ +\code{\link[rstatix:t_test]{rstatix::t_test}}, +\code{\link[daa:getWilcoxon]{getWilcoxon}} +} diff --git a/man/getWilcoxon.Rd b/man/getWilcoxon.Rd new file mode 100644 index 0000000..8e053e2 --- /dev/null +++ b/man/getWilcoxon.Rd @@ -0,0 +1,97 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/getWilcoxon.R +\name{getWilcoxon} +\alias{getWilcoxon} +\alias{addWilcoxon} +\alias{getWilcoxon,SingleCellExperiment-method} +\alias{getWilcoxon,SummarizedExperiment-method} +\alias{addWilcoxon,SummarizedExperiment-method} +\title{Perform Wilcoxon rank-sum test} +\usage{ +getWilcoxon(x, ...) + +addWilcoxon(x, ...) + +\S4method{getWilcoxon}{SingleCellExperiment}(x, ...) + +\S4method{getWilcoxon}{SummarizedExperiment}( + x, + assay.type = NULL, + row.var = NULL, + col.var = NULL, + formula, + split.by = NULL, + pair.by = NULL, + features = NULL, + p.adjust.method = "fdr", + ... +) + +\S4method{addWilcoxon}{SummarizedExperiment}(x, name = "wilcoxon", ...) +} +\arguments{ +\item{x}{A \code{SummarizedExperiment} object.} + +\item{...}{Additional arguments passed to \code{rstatix::wilcox_test}.} + +\item{assay.type}{\code{Character scalar} or \code{NULL}. Specifies assay to +test. Tests are run per feature. (Default: \code{NULL})} + +\item{row.var}{\code{Character scalar} or \code{NULL}. Specifies variable +from \code{rowData(x)} to test. (Default: \code{NULL})} + +\item{col.var}{\code{Character scalar} or \code{NULL}. Specifies variable +from \code{colData(x)} to test. (Default: \code{NULL})} + +\item{formula}{\code{formula}. A formula specifying the grouping variable, +e.g., \code{~ SampleType}. The RHS specifies the comparison groups. +For >2 levels, pairwise comparisons are performed.} + +\item{split.by}{\code{Character vector} or \code{NULL}. Columns to split by. +Tests are run separately for each combination. (Default: \code{NULL})} + +\item{pair.by}{\code{Character scalar} or \code{NULL}. Column for pairing +samples in paired tests. (Default: \code{NULL})} + +\item{features}{\code{Character vector} or \code{NULL}. Specific features +to test when using \code{assay.type}. (Default: \code{NULL})} + +\item{p.adjust.method}{\code{Character scalar}. Method for p-value +adjustment. (Default: \code{"fdr"})} + +\item{name}{\code{Character scalar}. Column name prefix for results. +(Default: \code{"wilcoxon"})} +} +\value{ +\code{getWilcoxon} returns a \code{DataFrame} with test results. +\code{addWilcoxon} returns \code{x} with results added to \code{rowData(x)}. +} +\description{ +These functions perform Wilcoxon rank-sum test (Mann-Whitney U) to compare +values between two groups. +} +\details{ +The Wilcoxon rank-sum test is a non-parametric test used to compare two +independent groups. It is suitable when data does not meet normality +assumptions. + +Specify exactly one of \code{assay.type}, \code{row.var}, or \code{col.var} +to define the values to test. +} +\examples{ +data(GlobalPatterns, package = "mia") +tse <- GlobalPatterns +tse <- tse[1:50, tse$SampleType \%in\% c("Feces", "Skin")] + +# Test assay values (per feature) +res <- getWilcoxon(tse, assay.type = "counts", formula = ~SampleType) + +# Test colData variable (e.g., alpha diversity) +tse <- mia::addAlpha(tse, index = "shannon_diversity") +res <- getWilcoxon(tse, col.var = "shannon_diversity", formula = ~SampleType) + +} +\seealso{ +\code{\link[rstatix:wilcox_test]{rstatix::wilcox_test}}, +\code{\link[daa:getTtest]{getTtest}} +} diff --git a/tests/testthat.R b/tests/testthat.R index 54620d4..f2759bb 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -1,10 +1,9 @@ # This file is part of the standard setup for testthat. # It is recommended that you do not modify it. # -# Where should you do additional test configuration? -# Learn more about the roles of various files in: -# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview -# * https://testthat.r-lib.org/articles/special-files.html +# Where should you do your testing? +# Learn more about the recommended testthat workflow here: +# https://r-pkgs.org/testing-design.html#sec-tests-files-overview library(testthat) library(daa) diff --git a/tests/testthat/test-example_test.R b/tests/testthat/test-example_test.R deleted file mode 100644 index 94b301e..0000000 --- a/tests/testthat/test-example_test.R +++ /dev/null @@ -1,3 +0,0 @@ -test_that("multiplication works", { - expect_equal(2 * 2, 4) -}) diff --git a/tests/testthat/test-getTtest.R b/tests/testthat/test-getTtest.R new file mode 100644 index 0000000..c2e7b7f --- /dev/null +++ b/tests/testthat/test-getTtest.R @@ -0,0 +1,144 @@ +# Tests for getTtest + +# Use HintikkaXOData - has good variance for statistical tests +# Need getWithColData to get colData from MAE +data(HintikkaXOData, package = "mia") +mae <- HintikkaXOData +tse <- MultiAssayExperiment::getWithColData(mae, "microbiota") +tse <- mia::transformAssay(tse, method = "relabundance") +tse <- tse[1:10, ] + +################################################################################ +# assay.type tests (per-feature) +################################################################################ + +test_that("getTtest with assay.type returns DataFrame", { + result <- getTtest(tse, + assay.type = "relabundance", + formula = ~Fat + ) + + expect_s4_class(result, "DataFrame") + expect_true("p" %in% names(result)) + expect_true("p.adj" %in% names(result)) +}) + +test_that("getTtest with assay.type returns one row per feature", { + result <- getTtest(tse, + assay.type = "relabundance", + formula = ~Fat + ) + expect_equal(nrow(result), nrow(tse)) +}) + +test_that("getTtest respects features argument", { + feat_subset <- rownames(tse)[1:5] + result <- getTtest(tse, + assay.type = "relabundance", + formula = ~Fat, features = feat_subset + ) + expect_equal(nrow(result), 5) +}) + +################################################################################ +# col.var tests (colData variable) +################################################################################ + +test_that("getTtest with col.var works", { + tse_alpha <- mia::addAlpha(tse, index = "shannon_diversity") + result <- getTtest(tse_alpha, + col.var = "shannon_diversity", + formula = ~Fat + ) + + expect_s4_class(result, "DataFrame") + expect_equal(nrow(result), 1) # Single test, not per-feature +}) + +################################################################################ +# Validation tests +################################################################################ + +test_that("getTtest fails without data source", { + expect_error(getTtest(tse, formula = ~Fat), "either") +}) + +test_that("getTtest fails with invalid formula", { + expect_error( + getTtest(tse, assay.type = "relabundance", formula = "not a formula"), + "must be a formula" + ) +}) + +test_that("getTtest var.equal parameter works", { + result_welch <- getTtest(tse, + assay.type = "relabundance", + formula = ~Fat, var.equal = FALSE + ) + result_student <- getTtest(tse, + assay.type = "relabundance", + formula = ~Fat, var.equal = TRUE + ) + expect_true(is.numeric(result_welch$p)) + expect_true(is.numeric(result_student$p)) +}) + +################################################################################ +# Robustness tests +################################################################################ + +test_that("getTtest handles zero variance with warnings", { + # Create a feature with zero variance to trigger error in t-test + tse_sparse <- tse + assay(tse_sparse, "relabundance")[1, ] <- 1 # All same values + + expect_warning( + result <- getTtest(tse_sparse, + assay.type = "relabundance", + formula = ~Fat + ), + "failed" + ) + expect_true(any(is.na(result$p))) + expect_true("warning" %in% colnames(result)) +}) + +test_that("getTtest works with un-sorted paired data", { + # Create paired data and shuffle it + tse_paired <- tse[, 1:10] + tse_paired$Subject <- rep(1:5, each = 2) + tse_paired$Time <- rep(c("T1", "T2"), 5) + + # Shuffle + idx <- sample(ncol(tse_paired)) + tse_shuffled <- tse_paired[, idx] + + # This should pass because we sort internally by Subject + expect_silent( + result <- getTtest(tse_shuffled, + assay.type = "relabundance", + formula = ~Time, pair.by = "Subject" + ) + ) + expect_s4_class(result, "DataFrame") +}) + +test_that("getTtest p.adjust.method parameter works", { + result <- getTtest(tse, + assay.type = "relabundance", + formula = ~Fat, p.adjust.method = "bonferroni" + ) + expect_true("p.adj" %in% names(result)) +}) + +################################################################################ +# addTtest tests +################################################################################ + +test_that("addTtest adds p.adj to rowData", { + result <- addTtest(tse, + assay.type = "relabundance", + formula = ~Fat + ) + expect_true("ttest_padj" %in% colnames(rowData(result))) +}) diff --git a/tests/testthat/test-getWilcoxon.R b/tests/testthat/test-getWilcoxon.R new file mode 100644 index 0000000..0c81d7e --- /dev/null +++ b/tests/testthat/test-getWilcoxon.R @@ -0,0 +1,130 @@ +# Tests for getWilcoxon + +# Use HintikkaXOData - has good variance for statistical tests +data(HintikkaXOData, package = "mia") +mae <- HintikkaXOData +tse <- MultiAssayExperiment::getWithColData(mae, "microbiota") +tse <- mia::transformAssay(tse, method = "relabundance") +tse <- tse[1:10, ] + +################################################################################ +# assay.type tests (per-feature) +################################################################################ + +test_that("getWilcoxon with assay.type returns DataFrame", { + result <- getWilcoxon(tse, + assay.type = "relabundance", + formula = ~Fat + ) + + expect_s4_class(result, "DataFrame") + expect_true("p" %in% names(result)) + expect_true("p.adj" %in% names(result)) +}) + +test_that("getWilcoxon with assay.type returns one row per feature", { + result <- getWilcoxon(tse, + assay.type = "relabundance", + formula = ~Fat + ) + expect_equal(nrow(result), nrow(tse)) +}) + +test_that("getWilcoxon respects features argument", { + feat_subset <- rownames(tse)[1:5] + result <- getWilcoxon(tse, + assay.type = "relabundance", + formula = ~Fat, features = feat_subset + ) + expect_equal(nrow(result), 5) +}) + +################################################################################ +# col.var tests (colData variable) +################################################################################ + +test_that("getWilcoxon with col.var works", { + tse_alpha <- mia::addAlpha(tse, index = "shannon_diversity") + result <- getWilcoxon(tse_alpha, + col.var = "shannon_diversity", + formula = ~Fat + ) + + expect_s4_class(result, "DataFrame") + expect_equal(nrow(result), 1) # Single test, not per-feature +}) + +################################################################################ +# Robustness tests +################################################################################ + +test_that("getWilcoxon handles sparse data with warnings", { + # Create a feature with all NAs to trigger error + tse_sparse <- tse + assay(tse_sparse, "relabundance")[1, ] <- NA_real_ + + expect_warning( + result <- getWilcoxon(tse_sparse, + assay.type = "relabundance", + formula = ~Fat + ), + "failed" + ) + expect_true(any(is.na(result$p))) + expect_true("warning" %in% colnames(result)) +}) + +test_that("getWilcoxon works with un-sorted paired data", { + # Create paired data and shuffle it + tse_paired <- tse[, 1:10] + tse_paired$Subject <- rep(1:5, each = 2) + tse_paired$Time <- rep(c("T1", "T2"), 5) + + # Shuffle + idx <- sample(ncol(tse_paired)) + tse_shuffled <- tse_paired[, idx] + + # This should pass because we sort internally by Subject + expect_silent( + result <- getWilcoxon(tse_shuffled, + assay.type = "relabundance", + formula = ~Time, pair.by = "Subject" + ) + ) + expect_s4_class(result, "DataFrame") +}) + +################################################################################ +# Validation tests +################################################################################ + +test_that("getWilcoxon fails without data source", { + expect_error(getWilcoxon(tse, formula = ~Fat), "either") +}) + +test_that("getWilcoxon fails with invalid formula", { + expect_error( + getWilcoxon(tse, assay.type = "relabundance", formula = "not a formula"), + "must be a formula" + ) +}) + +test_that("getWilcoxon p.adjust.method parameter works", { + result <- getWilcoxon(tse, + assay.type = "relabundance", + formula = ~Fat, p.adjust.method = "bonferroni" + ) + expect_true("p.adj" %in% names(result)) +}) + +################################################################################ +# addWilcoxon tests +################################################################################ + +test_that("addWilcoxon adds p.adj to rowData", { + result <- addWilcoxon(tse, + assay.type = "relabundance", + formula = ~Fat + ) + expect_true("wilcoxon_padj" %in% colnames(rowData(result))) +}) diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R new file mode 100644 index 0000000..1c4e2ec --- /dev/null +++ b/tests/testthat/test-utils.R @@ -0,0 +1,149 @@ +# Tests for utils.R functions + +# Use HintikkaXOData for reliable tests +# Need getWithColData to get colData from MAE +data(HintikkaXOData, package = "mia") +mae <- HintikkaXOData +tse <- MultiAssayExperiment::getWithColData(mae, "microbiota") +tse <- mia::transformAssay(tse, method = "relabundance") +tse <- tse[1:10, ] + +################################################################################ +# .check_input tests +################################################################################ + +test_that(".check_input requires exactly one data source", { + expect_error( + .check_input(tse, NULL, NULL, NULL, ~Fat, NULL, NULL, NULL), + "either" + ) + expect_error( + .check_input( + tse, "relabundance", "var", NULL, ~Fat, + NULL, NULL, NULL + ), + "either" + ) +}) + +test_that(".check_input validates features only with assay.type", { + expect_error( + .check_input(tse, NULL, NULL, "Fat", ~Fat, NULL, NULL, + features = "Taxa1" + ), + "features" + ) +}) + +test_that(".check_input validates features exist", { + expect_error( + .check_input(tse, "relabundance", NULL, NULL, ~Fat, NULL, NULL, + features = c("nonexistent") + ), + "features" + ) +}) + +test_that(".check_input validates formula", { + expect_error( + .check_input( + tse, "relabundance", NULL, NULL, "not a formula", + NULL, NULL, NULL + ), + "must be a formula" + ) +}) + +test_that(".check_input validates formula variable exists", { + expect_error( + .check_input( + tse, "relabundance", NULL, NULL, ~nonexistent, + NULL, NULL, NULL + ), + "must be a single character value from the following options" + ) +}) + +test_that(".check_input validates split.by exists", { + expect_error( + .check_input( + tse, "relabundance", NULL, NULL, ~Fat, + "nonexistent", NULL, NULL + ), + "must be a single character value from the following options" + ) +}) + +test_that(".check_input validates pair.by exists", { + expect_error( + .check_input( + tse, "relabundance", NULL, NULL, ~Fat, NULL, + "nonexistent", NULL + ), + "must be a single character value from the following options" + ) +}) + +test_that(".check_input passes with valid inputs", { + expect_silent( + .check_input( + tse, "relabundance", NULL, NULL, ~Fat, + NULL, NULL, NULL + ) + ) +}) + +################################################################################ +# .get_rhs_var tests +################################################################################ + +test_that(".get_rhs_var extracts group from formula", { + result <- .get_rhs_var(~Fat) + expect_equal(result, "Fat") +}) + +test_that(".get_rhs_var fails with non-formula", { + expect_error(.get_rhs_var("not a formula"), "must be a formula") +}) + +################################################################################ +# .check_group tests +################################################################################ + +test_that(".check_group validates at least 2 levels", { + df <- data.frame(group = c("A"), value = 1) + expect_error(.check_group(df, "group"), "at least 2 levels") +}) + +test_that(".check_group passes with valid 2-level group", { + df <- data.frame(group = c("A", "B"), value = 1:2) + expect_silent(.check_group(df, "group")) +}) + +################################################################################ +# .get_data tests +################################################################################ + +test_that(".get_data with assay.type returns correct format", { + df <- .get_data( + tse, "relabundance", NULL, NULL, "Fat", + NULL, NULL, NULL + ) + + expect_s3_class(df, "data.frame") + expect_true("relabundance" %in% names(df)) + expect_true("Fat" %in% names(df)) + expect_true("FeatureID" %in% names(df)) +}) + +test_that(".get_data with col.var returns correct format", { + tse_alpha <- mia::addAlpha(tse, index = "shannon_diversity") + df <- .get_data( + tse_alpha, NULL, NULL, "shannon_diversity", + "Fat", NULL, NULL, NULL + ) + + expect_s3_class(df, "data.frame") + expect_true("shannon_diversity" %in% names(df)) + expect_true("Fat" %in% names(df)) +}) diff --git a/vignettes/daa.Rmd b/vignettes/daa.Rmd index 86bb835..1348a11 100644 --- a/vignettes/daa.Rmd +++ b/vignettes/daa.Rmd @@ -1,19 +1,9 @@ --- title: "Introduction to daa" -author: - - name: Your name - affiliation: - - Your institution - email: your.email@somewhere.com output: BiocStyle::html_document: - self_contained: yes toc: true toc_float: true - toc_depth: 2 - code_folding: show -date: "`r doc_date()`" -package: "`r pkg_ver('daa')`" vignette: > %\VignetteIndexEntry{Introduction to daa} %\VignetteEngine{knitr::rmarkdown} @@ -22,113 +12,85 @@ vignette: > ```{r setup, include = FALSE} knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>", - crop = NULL ## Related to https://stat.ethz.ch/pipermail/bioc-devel/2020-April/016656.html + collapse = TRUE, + comment = "#>" ) ``` +# Introduction -```{r vignetteSetup, echo=FALSE, message=FALSE, warning = FALSE} -## Bib setup -library("RefManageR") - -## Write bibliography information -bib <- c( - R = citation(), - BiocStyle = citation("BiocStyle")[1], - knitr = citation("knitr")[1], - RefManageR = citation("RefManageR")[1], - rmarkdown = citation("rmarkdown")[1], - sessioninfo = citation("sessioninfo")[1], - testthat = citation("testthat")[1], - daa = citation("daa")[1] -) -``` - -# Basics - -## Install `daa` +The **daa** package provides methods for differential abundance analysis (DAA) +of microbiome data. It works with `TreeSummarizedExperiment` and +`SummarizedExperiment` objects from the miaverse ecosystem. -`R` is an open-source statistical environment which can be easily modified to enhance its functionality via packages. `r Biocpkg("daa")` is a `R` package available via the [Bioconductor](http://bioconductor.org) repository for packages. `R` can be installed on any operating system from [CRAN](https://cran.r-project.org/) after which you can install `r Biocpkg("daa")` by using the following commands in your `R` session: +## Installation -```{r "install", eval = FALSE} +```{r install, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) { - install.packages("BiocManager") + install.packages("BiocManager") } BiocManager::install("daa") - -## Check that you have a valid Bioconductor installation -BiocManager::valid() ``` -## Required knowledge +# Available Methods -`r Biocpkg("daa")` is based on many other packages and in particular in those that have implemented the infrastructure needed for dealing with RNA-seq data (EDIT!). That is, packages like `r Biocpkg("SummarizedExperiment")` (EDIT!). +## Wilcoxon rank-sum test -If you are asking yourself the question "Where do I start using Bioconductor?" you might be interested in [this blog post](http://lcolladotor.github.io/2014/10/16/startbioc/#.VkOKbq6rRuU). +The `getWilcoxon()` function performs nonparametric Wilcoxon rank-sum tests +for each feature, comparing two groups. -## Asking for help +```{r wilcoxon, message = FALSE} +library(daa) +library(mia) -As package developers, we try to explain clearly how to use our packages and in which order to use the functions. But `R` and `Bioconductor` have a steep learning curve so it is critical to learn where to ask for help. The blog post quoted above mentions some but we would like to highlight the [Bioconductor support site](https://support.bioconductor.org/) as the main resource for getting help: remember to use the `daa` tag and check [the older posts](https://support.bioconductor.org/tag/daa/). Other alternatives are available such as creating GitHub issues and tweeting. However, please note that if you want to receive help you should adhere to the [posting guidelines](http://www.bioconductor.org/help/support/posting-guide/). It is particularly critical that you provide a small reproducible example and your session information so package developers can track down the source of the error. +# Load example data (HintikkaXOData is used in unit tests) +data(HintikkaXOData, package = "mia") +tse <- MultiAssayExperiment::getWithColData(HintikkaXOData, "microbiota") -## Citing `daa` +# Standardize assay data +tse <- transformAssay(tse, method = "relabundance") -We hope that `r Biocpkg("daa")` will be useful for your research. Please use the following information to cite the package and the overall approach. Thank you! +# Subset for demonstration +tse <- tse[1:50, ] -```{r "citation"} -## Citation info -citation("daa") +# Run Wilcoxon test comparing groups in 'Fat' metadata column +res <- getWilcoxon(tse, assay.type = "relabundance", formula = ~Fat) +head(res) ``` -# Quick start to using `daa` +Use `addWilcoxon()` to store results in `rowData()`: -```{r "start", message=FALSE} -library("daa") +```{r add_wilcoxon} +tse <- addWilcoxon(tse, assay.type = "relabundance", formula = ~Fat) +head(rowData(tse)[, "wilcoxon_padj", drop = FALSE]) ``` -Edit this as you see fit =) - -Here is an example of you can cite your package inside the vignette: - -* `r Biocpkg("daa")` `r Citep(bib[["daa"]])` +## t-test +The `getTtest()` function performs t-tests (Welch's by default). - -# Reproducibility - -The `r Biocpkg("daa")` package `r Citep(bib[["daa"]])` was made possible thanks to: - -* R `r Citep(bib[["R"]])` -* `r Biocpkg("BiocStyle")` `r Citep(bib[["BiocStyle"]])` -* `r CRANpkg("knitr")` `r Citep(bib[["knitr"]])` -* `r CRANpkg("RefManageR")` `r Citep(bib[["RefManageR"]])` -* `r CRANpkg("rmarkdown")` `r Citep(bib[["rmarkdown"]])` -* `r CRANpkg("sessioninfo")` `r Citep(bib[["sessioninfo"]])` -* `r CRANpkg("testthat")` `r Citep(bib[["testthat"]])` - -This package was developed using `r BiocStyle::Biocpkg("biocthis")`. - -`R` session information. - -```{r reproduce3, echo=FALSE} -## Session info -library("sessioninfo") -options(width = 120) -session_info() +```{r ttest} +res <- getTtest(tse, assay.type = "relabundance", formula = ~Fat) +head(res) ``` +## Testing metadata variables +You can also test variables from `colData` or `rowData`. For example, testing +alpha diversity indices: -# Bibliography +```{r metadata} +# Add alpha diversity to colData +tse <- mia::addAlpha(tse, index = "shannon_diversity") -This vignette was generated using `r Biocpkg("BiocStyle")` `r Citep(bib[["BiocStyle"]])` -with `r CRANpkg("knitr")` `r Citep(bib[["knitr"]])` and `r CRANpkg("rmarkdown")` `r Citep(bib[["rmarkdown"]])` running behind the scenes. +# Test diversity differences between groups +res <- getWilcoxon(tse, col.var = "shannon_diversity", formula = ~Fat) +head(res) +``` -Citations made with `r CRANpkg("RefManageR")` `r Citep(bib[["RefManageR"]])`. +# Session Info -```{r vignetteBiblio, results = "asis", echo = FALSE, warning = FALSE, message = FALSE} -## Print bibliography -PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html")) +```{r session} +sessionInfo() ```