From 989206b9b082cdd2ad775c92809b0d7faa77853a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Th=C3=A9otime=20Pralas?= <151254073+thpralas@users.noreply.github.com> Date: Wed, 8 May 2024 16:17:19 +0300 Subject: [PATCH] Replace getExperiment* and testExperiment* functions with getCrossAssociation (#511) --- DESCRIPTION | 2 +- NAMESPACE | 2 + NEWS | 1 + R/deprecate.R | 80 +++ ...ossAssociation.R => getCrossAssociation.R} | 118 +--- man/deprecate.Rd | 27 + ...sAssociation.Rd => getCrossAssociation.Rd} | 87 ++- pkgdown/_pkgdown.yml | 2 +- tests/testthat/test-5getCrossAssociation.R | 566 ++++++++++++++++++ .../test-5getExperimentCrossAssociation.R | 540 ----------------- 10 files changed, 745 insertions(+), 680 deletions(-) rename R/{getExperimentCrossAssociation.R => getCrossAssociation.R} (93%) rename man/{getExperimentCrossAssociation.Rd => getCrossAssociation.Rd} (77%) create mode 100644 tests/testthat/test-5getCrossAssociation.R delete mode 100644 tests/testthat/test-5getExperimentCrossAssociation.R diff --git a/DESCRIPTION b/DESCRIPTION index da6cce0f5..d4fdd9f96 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mia Type: Package -Version: 1.13.1 +Version: 1.13.2 Authors@R: c(person(given = "Felix G.M.", family = "Ernst", role = c("aut"), email = "felix.gm.ernst@outlook.com", diff --git a/NAMESPACE b/NAMESPACE index fd90803e8..0c47fba0c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -32,6 +32,7 @@ export(estimateFaith) export(estimateRichness) export(full_join) export(getBestDMNFit) +export(getCrossAssociation) export(getDMN) export(getExperimentCrossAssociation) export(getExperimentCrossCorrelation) @@ -145,6 +146,7 @@ exportMethods(estimateFaith) exportMethods(estimateRichness) exportMethods(full_join) exportMethods(getBestDMNFit) +exportMethods(getCrossAssociation) exportMethods(getDMN) exportMethods(getExperimentCrossAssociation) exportMethods(getExperimentCrossCorrelation) diff --git a/NEWS b/NEWS index 7eb4e182a..0e675243e 100644 --- a/NEWS +++ b/NEWS @@ -115,3 +115,4 @@ Changes in version 1.11.x Changes in version 1.13.x + Added new functions getMediation and addMediation ++ replace getExperiment* and testExperiment* functions with getCrossAssociation diff --git a/R/deprecate.R b/R/deprecate.R index 83fdf012d..cf42cd9e0 100644 --- a/R/deprecate.R +++ b/R/deprecate.R @@ -62,6 +62,86 @@ setMethod("taxonomyTree", signature = c(x = "SummarizedExperiment"), #' @rdname deprecate #' @export +setGeneric("getExperimentCrossAssociation", signature = c("x"), + function(x, ...) + standardGeneric("getExperimentCrossAssociation")) + +#' @rdname deprecate +#' @export +setMethod("getExperimentCrossAssociation", + signature = c(x = "MultiAssayExperiment"), + function(x, ...){ + .Deprecated(msg = paste0("'getExperimentCrossAssociation' is ", + "deprecated. Use ", + "'getCrossAssociation' instead.")) + getCrossAssociation(x, ...) + } +) + +#' @rdname deprecate +#' @export +setMethod("getExperimentCrossAssociation", signature = "SummarizedExperiment", + function(x, ...){ + .Deprecated(msg = paste0("'getExperimentCrossAssociation' is ", + "deprecated. Use ", + "'getCrossAssociation' instead.")) + getCrossAssociation(x, ...) + } +) + +#' @rdname deprecate +#' @export +setGeneric("testExperimentCrossAssociation", signature = c("x"), + function(x, ...) + standardGeneric("testExperimentCrossAssociation")) + +#' @rdname deprecate +#' @export +setMethod("testExperimentCrossAssociation", signature = c(x = "ANY"), + function(x, ...){ + .Deprecated(msg = paste0("'testExperimentCrossAssociation' is ", + "deprecated. Use ", + "'getCrossAssociation' instead.")) + getCrossAssociation(x, test_significance = TRUE, ...) + } +) + +#' @rdname deprecate +#' @export +setGeneric("testExperimentCrossCorrelation", signature = c("x"), + function(x, ...) + standardGeneric("testExperimentCrossCorrelation")) + +#' @rdname deprecate +#' @export +setMethod("testExperimentCrossCorrelation", signature = c(x = "ANY"), + function(x, ...){ + .Deprecated(msg = paste0("'testExperimentCrossCorrelation' is ", + "deprecated. Use ", + "'getCrossAssociation' instead.")) + getCrossAssociation(x, test_significance = TRUE, ...) + } +) + +#' @rdname deprecate +#' @export +setGeneric("getExperimentCrossCorrelation", signature = c("x"), + function(x, ...) + standardGeneric("getExperimentCrossCorrelation")) + +#' @rdname deprecate +#' @export +setMethod("getExperimentCrossCorrelation", signature = c(x = "ANY"), + function(x, ...){ + .Deprecated(msg = paste0("'getExperimentCrossCorrelation' is ", + "deprecated. Use ", + "'getCrossAssociation' instead.")) + getCrossAssociation(x, ...) + } +) + +#' @rdname deprecate +#' @export setGeneric("mergeFeaturesByPrevalence", signature = "x", function(x, ...) standardGeneric("mergeFeaturesByPrevalence")) diff --git a/R/getExperimentCrossAssociation.R b/R/getCrossAssociation.R similarity index 93% rename from R/getExperimentCrossAssociation.R rename to R/getCrossAssociation.R index 2d816003f..1dca1b6e7 100644 --- a/R/getExperimentCrossAssociation.R +++ b/R/getCrossAssociation.R @@ -89,6 +89,7 @@ #' #' @param test_significance A single boolean value for selecting whether to test #' statistical significance of associations. +#' (By default: \code{test_significance = FALSE}) #' #' @param show_warnings A single boolean value for selecting whether to show warnings #' that might occur when correlations and p-values are calculated. @@ -111,10 +112,10 @@ #' } #' #' @details -#' These functions calculates associations between features of two experiments. -#' \code{getExperimentCrossAssociation} calculates only associations by default. -#' \code{testExperimentCrossAssociation} calculates also significance of -#' associations. +#' The function \code{getCrossAssociation} calculates associations between +#' features of two experiments. By default, it not only computes associations +#' but also tests their significance. If desired, setting +#' \code{test_significance} to FALSE disables significance calculation. #' #' We recommend the non-parametric Kendall's tau as the default method for association #' analysis. Kendall's tau has desirable statistical properties and robustness at lower @@ -122,13 +123,13 @@ #' running times are critical. #' #' @return -#' These functions return associations in table or matrix format. In table format, +#' This function returns associations in table or matrix format. In table format, #' returned value is a data frame that includes features and associations #' (and p-values) in columns. In matrix format, returned value is a one matrix #' when only associations are calculated. If also significances are tested, then #' returned value is a list of matrices. #' -#' @name getExperimentCrossAssociation +#' @name getCrossAssociation #' @export #' #' @author Leo Lahti and Tuomas Borman. Contact: \url{microbiome.github.io} @@ -147,7 +148,7 @@ #' mae[[1]] <- transformAssay(mae[[1]], method = "rclr") #' #' # Calculate cross-correlations -#' result <- getExperimentCrossAssociation(mae, method = "pearson", assay.type2 = "nmr") +#' result <- getCrossAssociation(mae, method = "pearson", assay.type2 = "nmr") #' # Show first 5 entries #' head(result, 5) #' @@ -156,47 +157,44 @@ #' # Transform data #' altExp(mae[[1]], "Phylum") <- transformAssay(altExp(mae[[1]], "Phylum"), method = "relabundance") #' # When mode = "matrix", the return value is a matrix -#' result <- getExperimentCrossAssociation(mae, experiment2 = 2, +#' result <- getCrossAssociation(mae, experiment2 = 2, #' assay.type1 = "relabundance", assay.type2 = "nmr", #' altexp1 = "Phylum", #' method = "pearson", mode = "matrix") #' # Show first 5 entries #' head(result, 5) #' -#' # testExperimentCorrelation additionally returns significances +#' # If test_significance = TRUE, then getCrossAssociation additionally returns +#' # significances #' # filter_self_correlations = TRUE filters self correlations #' # p_adj_threshold can be used to filter those features that do not #' # have any correlations whose p-value is lower than the threshold -#' result <- testExperimentCrossAssociation(mae[[1]], experiment2 = mae[[1]], method = "pearson", +#' result <- getCrossAssociation(mae[[1]], experiment2 = mae[[1]], method = "pearson", #' filter_self_correlations = TRUE, -#' p_adj_threshold = 0.05) +#' p_adj_threshold = 0.05, +#' test_significance = TRUE) #' # Show first 5 entries #' head(result, 5) -#' -#' # getExperimentCrossAssociation also returns significances when -#' # test_significance = TRUE -#' # Warnings can be suppressed by using show_warnings = FALSE -#' result <- getExperimentCrossAssociation(mae[[1]], experiment2 = mae[[2]], method = "pearson", -#' assay.type2 = "nmr", -#' mode = "matrix", test_significance = TRUE, -#' show_warnings = FALSE) #' #' # Returned value is a list of matrices #' names(result) #' #' # Calculate Bray-Curtis dissimilarity between samples. If dataset includes #' # paired samples, you can use paired = TRUE. -#' result <- getExperimentCrossAssociation(mae[[1]], mae[[1]], MARGIN = 2, paired = FALSE, -#' association_FUN = vegan::vegdist, method = "bray") +#' result <- getCrossAssociation(mae[[1]], mae[[1]], MARGIN = 2, paired = FALSE, +#' association_FUN = vegan::vegdist, +#' method = "bray") #' #' #' # If experiments are equal and measure is symmetric (e.g., taxa1 vs taxa2 == taxa2 vs taxa1), #' # it is possible to speed-up calculations by calculating association only for unique #' # variable-pairs. Use "symmetric" to choose whether to measure association for only #' # other half of of variable-pairs. -#' result <- getExperimentCrossAssociation(mae, experiment1 = "microbiota", experiment2 = "microbiota", -#' assay.type1 = "counts", assay.type2 = "counts", -#' symmetric = TRUE) +#' result <- getCrossAssociation(mae, experiment1 = "microbiota", +#' experiment2 = "microbiota", +#' assay.type1 = "counts", +#' assay.type2 = "counts", +#' symmetric = TRUE) #' #' # For big data sets, the calculations might take a long time. #' # To speed them up, you can take a random sample from the data. @@ -205,7 +203,7 @@ #' sample_size <- 0.3 #' tse <- mae[[1]] #' tse_sub <- tse[ sample( seq_len( nrow(tse) ), sample_size * nrow(tse) ), ] -#' result <- testExperimentCrossAssociation(tse_sub) +#' result <- getCrossAssociation(tse_sub) #' #' # It is also possible to choose variables from colData and calculate association #' # between assay and sample metadata or between variables of sample metadata @@ -213,22 +211,21 @@ #' # colData_variable works similarly to assay.type. Instead of fetching an assay #' # named assay.type from assay slot, it fetches a column named colData_variable #' # from colData. -#' result <- getExperimentCrossAssociation(mae[[1]], assay.type1 = "counts", -#' colData_variable2 = c("shannon", "coverage")) +#' result <- getCrossAssociation(mae[[1]], assay.type1 = "counts", +#' colData_variable2 = c("shannon", "coverage"), +#' test_significance = TRUE) #' NULL -#' @rdname getExperimentCrossAssociation -#' @aliases getExperimentCrossCorrelation +#' @rdname getCrossAssociation #' @export -setGeneric("getExperimentCrossAssociation", signature = c("x"), +setGeneric("getCrossAssociation", signature = c("x"), function(x, ...) - standardGeneric("getExperimentCrossAssociation")) + standardGeneric("getCrossAssociation")) -#' @rdname getExperimentCrossAssociation -#' @aliases getExperimentCrossCorrelation +#' @rdname getCrossAssociation #' @export -setMethod("getExperimentCrossAssociation", signature = c(x = "MultiAssayExperiment"), +setMethod("getCrossAssociation", signature = c(x = "MultiAssayExperiment"), function(x, experiment1 = 1, experiment2 = 2, @@ -277,12 +274,11 @@ setMethod("getExperimentCrossAssociation", signature = c(x = "MultiAssayExperime } ) -#' @rdname getExperimentCrossAssociation -#' @aliases getExperimentCrossCorrelation +#' @rdname getCrossAssociation #' @importFrom MultiAssayExperiment MultiAssayExperiment ExperimentList #' @importFrom SingleCellExperiment altExps #' @export -setMethod("getExperimentCrossAssociation", signature = "SummarizedExperiment", +setMethod("getCrossAssociation", signature = "SummarizedExperiment", function(x, experiment2 = x, ...){ ############################## INPUT CHECK ############################# # If y is SE or TreeSE object @@ -341,54 +337,6 @@ setMethod("getExperimentCrossAssociation", signature = "SummarizedExperiment", } ) -#' @rdname getExperimentCrossAssociation -#' @aliases testExperimentCrossCorrelation -#' @export -setGeneric("testExperimentCrossAssociation", signature = c("x"), - function(x, ...) - standardGeneric("testExperimentCrossAssociation")) - -#' @rdname getExperimentCrossAssociation -#' @aliases testExperimentCrossCorrelation -#' @export -setMethod("testExperimentCrossAssociation", signature = c(x = "ANY"), - function(x, ...){ - getExperimentCrossAssociation(x, test_significance = TRUE, ...) - } -) - -############################# Methods for aliases ############################## -#' @rdname getExperimentCrossAssociation -#' @aliases testExperimentCrossAssociation -#' @export -setGeneric("testExperimentCrossCorrelation", signature = c("x"), - function(x, ...) - standardGeneric("testExperimentCrossCorrelation")) - -#' @rdname getExperimentCrossAssociation -#' @aliases testExperimentCrossAssociation -#' @export -setMethod("testExperimentCrossCorrelation", signature = c(x = "ANY"), - function(x, ...){ - getExperimentCrossAssociation(x, test_significance = TRUE, ...) - } -) - -#' @rdname getExperimentCrossAssociation -#' @aliases getExperimentCrossAssociation -#' @export -setGeneric("getExperimentCrossCorrelation", signature = c("x"), - function(x, ...) - standardGeneric("getExperimentCrossCorrelation")) - -#' @rdname getExperimentCrossAssociation -#' @aliases getExperimentCrossAssociation -#' @export -setMethod("getExperimentCrossCorrelation", signature = c(x = "ANY"), - function(x, ...){ - getExperimentCrossAssociation(x, ...) - } -) ############################## MAIN FUNCTIONALITY ############################## # This function includes all the main functionality. #' @importFrom S4Vectors unfactor diff --git a/man/deprecate.Rd b/man/deprecate.Rd index cb3d4f16e..01744dda0 100644 --- a/man/deprecate.Rd +++ b/man/deprecate.Rd @@ -8,6 +8,15 @@ \alias{addTaxonomyTree,SummarizedExperiment-method} \alias{taxonomyTree} \alias{taxonomyTree,SummarizedExperiment-method} +\alias{getExperimentCrossAssociation} +\alias{getExperimentCrossAssociation,MultiAssayExperiment-method} +\alias{getExperimentCrossAssociation,SummarizedExperiment-method} +\alias{testExperimentCrossAssociation} +\alias{testExperimentCrossAssociation,ANY-method} +\alias{testExperimentCrossCorrelation} +\alias{testExperimentCrossCorrelation,ANY-method} +\alias{getExperimentCrossCorrelation} +\alias{getExperimentCrossCorrelation,ANY-method} \alias{mergeFeaturesByPrevalence} \alias{mergeFeaturesByPrevalence,SummarizedExperiment-method} \alias{loadFromBiom} @@ -41,6 +50,24 @@ taxonomyTree(x, ...) \S4method{taxonomyTree}{SummarizedExperiment}(x, ...) +getExperimentCrossAssociation(x, ...) + +\S4method{getExperimentCrossAssociation}{MultiAssayExperiment}(x, ...) + +\S4method{getExperimentCrossAssociation}{SummarizedExperiment}(x, ...) + +testExperimentCrossAssociation(x, ...) + +\S4method{testExperimentCrossAssociation}{ANY}(x, ...) + +testExperimentCrossCorrelation(x, ...) + +\S4method{testExperimentCrossCorrelation}{ANY}(x, ...) + +getExperimentCrossCorrelation(x, ...) + +\S4method{getExperimentCrossCorrelation}{ANY}(x, ...) + mergeFeaturesByPrevalence(x, ...) \S4method{mergeFeaturesByPrevalence}{SummarizedExperiment}(x, ...) diff --git a/man/getExperimentCrossAssociation.Rd b/man/getCrossAssociation.Rd similarity index 77% rename from man/getExperimentCrossAssociation.Rd rename to man/getCrossAssociation.Rd index 19c8410cd..963d73264 100644 --- a/man/getExperimentCrossAssociation.Rd +++ b/man/getCrossAssociation.Rd @@ -1,20 +1,14 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/getExperimentCrossAssociation.R -\name{getExperimentCrossAssociation} -\alias{getExperimentCrossAssociation} -\alias{getExperimentCrossCorrelation} -\alias{getExperimentCrossAssociation,MultiAssayExperiment-method} -\alias{getExperimentCrossAssociation,SummarizedExperiment-method} -\alias{testExperimentCrossAssociation} -\alias{testExperimentCrossCorrelation} -\alias{testExperimentCrossAssociation,ANY-method} -\alias{testExperimentCrossCorrelation,ANY-method} -\alias{getExperimentCrossCorrelation,ANY-method} +% Please edit documentation in R/getCrossAssociation.R +\name{getCrossAssociation} +\alias{getCrossAssociation} +\alias{getCrossAssociation,MultiAssayExperiment-method} +\alias{getCrossAssociation,SummarizedExperiment-method} \title{Calculate correlations between features of two experiments.} \usage{ -getExperimentCrossAssociation(x, ...) +getCrossAssociation(x, ...) -\S4method{getExperimentCrossAssociation}{MultiAssayExperiment}( +\S4method{getCrossAssociation}{MultiAssayExperiment}( x, experiment1 = 1, experiment2 = 2, @@ -41,19 +35,7 @@ getExperimentCrossAssociation(x, ...) ... ) -\S4method{getExperimentCrossAssociation}{SummarizedExperiment}(x, experiment2 = x, ...) - -testExperimentCrossAssociation(x, ...) - -\S4method{testExperimentCrossAssociation}{ANY}(x, ...) - -testExperimentCrossCorrelation(x, ...) - -\S4method{testExperimentCrossCorrelation}{ANY}(x, ...) - -getExperimentCrossCorrelation(x, ...) - -\S4method{getExperimentCrossCorrelation}{ANY}(x, ...) +\S4method{getCrossAssociation}{SummarizedExperiment}(x, experiment2 = x, ...) } \arguments{ \item{x}{A @@ -157,7 +139,8 @@ between experiment itself is tested, i.e., when assays are identical. about progress of calculation.} \item{test_significance}{A single boolean value for selecting whether to test -statistical significance of associations.} +statistical significance of associations. +(By default: \code{test_significance = FALSE})} \item{show_warnings}{A single boolean value for selecting whether to show warnings that might occur when correlations and p-values are calculated.} @@ -167,7 +150,7 @@ that might occur when correlations and p-values are calculated.} when \code{MARGIN = 1}. (By default: \code{paired = FALSE})} } \value{ -These functions return associations in table or matrix format. In table format, +This function returns associations in table or matrix format. In table format, returned value is a data frame that includes features and associations (and p-values) in columns. In matrix format, returned value is a one matrix when only associations are calculated. If also significances are tested, then @@ -177,10 +160,10 @@ returned value is a list of matrices. Calculate correlations between features of two experiments. } \details{ -These functions calculates associations between features of two experiments. -\code{getExperimentCrossAssociation} calculates only associations by default. -\code{testExperimentCrossAssociation} calculates also significance of -associations. +The function \code{getCrossAssociation} calculates associations between +features of two experiments. By default, it not only computes associations +but also tests their significance. If desired, setting +\code{test_significance} to FALSE disables significance calculation. We recommend the non-parametric Kendall's tau as the default method for association analysis. Kendall's tau has desirable statistical properties and robustness at lower @@ -201,7 +184,7 @@ mae[[1]] <- mae[[1]][rowSds(assay(mae[[1]])) > 0, ] mae[[1]] <- transformAssay(mae[[1]], method = "rclr") # Calculate cross-correlations -result <- getExperimentCrossAssociation(mae, method = "pearson", assay.type2 = "nmr") +result <- getCrossAssociation(mae, method = "pearson", assay.type2 = "nmr") # Show first 5 entries head(result, 5) @@ -210,47 +193,44 @@ altExp(mae[[1]], "Phylum") <- agglomerateByRank(mae[[1]], rank = "Phylum") # Transform data altExp(mae[[1]], "Phylum") <- transformAssay(altExp(mae[[1]], "Phylum"), method = "relabundance") # When mode = "matrix", the return value is a matrix -result <- getExperimentCrossAssociation(mae, experiment2 = 2, +result <- getCrossAssociation(mae, experiment2 = 2, assay.type1 = "relabundance", assay.type2 = "nmr", altexp1 = "Phylum", method = "pearson", mode = "matrix") # Show first 5 entries head(result, 5) -# testExperimentCorrelation additionally returns significances +# If test_significance = TRUE, then getCrossAssociation additionally returns +# significances # filter_self_correlations = TRUE filters self correlations # p_adj_threshold can be used to filter those features that do not # have any correlations whose p-value is lower than the threshold -result <- testExperimentCrossAssociation(mae[[1]], experiment2 = mae[[1]], method = "pearson", +result <- getCrossAssociation(mae[[1]], experiment2 = mae[[1]], method = "pearson", filter_self_correlations = TRUE, - p_adj_threshold = 0.05) + p_adj_threshold = 0.05, + test_significance = TRUE) # Show first 5 entries head(result, 5) - -# getExperimentCrossAssociation also returns significances when -# test_significance = TRUE -# Warnings can be suppressed by using show_warnings = FALSE -result <- getExperimentCrossAssociation(mae[[1]], experiment2 = mae[[2]], method = "pearson", - assay.type2 = "nmr", - mode = "matrix", test_significance = TRUE, - show_warnings = FALSE) # Returned value is a list of matrices names(result) # Calculate Bray-Curtis dissimilarity between samples. If dataset includes # paired samples, you can use paired = TRUE. -result <- getExperimentCrossAssociation(mae[[1]], mae[[1]], MARGIN = 2, paired = FALSE, - association_FUN = vegan::vegdist, method = "bray") +result <- getCrossAssociation(mae[[1]], mae[[1]], MARGIN = 2, paired = FALSE, + association_FUN = vegan::vegdist, + method = "bray") # If experiments are equal and measure is symmetric (e.g., taxa1 vs taxa2 == taxa2 vs taxa1), # it is possible to speed-up calculations by calculating association only for unique # variable-pairs. Use "symmetric" to choose whether to measure association for only # other half of of variable-pairs. -result <- getExperimentCrossAssociation(mae, experiment1 = "microbiota", experiment2 = "microbiota", - assay.type1 = "counts", assay.type2 = "counts", - symmetric = TRUE) +result <- getCrossAssociation(mae, experiment1 = "microbiota", + experiment2 = "microbiota", + assay.type1 = "counts", + assay.type2 = "counts", + symmetric = TRUE) # For big data sets, the calculations might take a long time. # To speed them up, you can take a random sample from the data. @@ -259,7 +239,7 @@ result <- getExperimentCrossAssociation(mae, experiment1 = "microbiota", experim sample_size <- 0.3 tse <- mae[[1]] tse_sub <- tse[ sample( seq_len( nrow(tse) ), sample_size * nrow(tse) ), ] -result <- testExperimentCrossAssociation(tse_sub) +result <- getCrossAssociation(tse_sub) # It is also possible to choose variables from colData and calculate association # between assay and sample metadata or between variables of sample metadata @@ -267,8 +247,9 @@ mae[[1]] <- estimateDiversity(mae[[1]]) # colData_variable works similarly to assay.type. Instead of fetching an assay # named assay.type from assay slot, it fetches a column named colData_variable # from colData. -result <- getExperimentCrossAssociation(mae[[1]], assay.type1 = "counts", - colData_variable2 = c("shannon", "coverage")) +result <- getCrossAssociation(mae[[1]], assay.type1 = "counts", + colData_variable2 = c("shannon", "coverage"), + test_significance = TRUE) } \author{ diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index bec02517e..88d0c7669 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -46,8 +46,8 @@ reference: - getPrevalence - getPrevalentAbundance - getTopTaxa + - getCrossAssociation - getMediation - - getExperimentCrossCorrelation - hierarchy-tree - subtitle: Dominance - contents: diff --git a/tests/testthat/test-5getCrossAssociation.R b/tests/testthat/test-5getCrossAssociation.R new file mode 100644 index 000000000..7a0d9ac7e --- /dev/null +++ b/tests/testthat/test-5getCrossAssociation.R @@ -0,0 +1,566 @@ + +context("getCrossAssociation") + +test_that("getCrossAssociation", { + + # Try 5 times to fetch the data + for(i in seq_len(5) ){ + mae <- tryCatch( + { + # Try to fetch the data + microbiomeDataSets::peerj32() + }, + error = function(cond) { + # If it was not possible to fetch the data, give FALSE + return(NULL) + } + ) + # Break if mae has the data + if( !is.null(mae) ){ + break + } + } + # Run tests if the data fetch was successful + if( !is.null(mae) ){ + ############################### Test input ############################### + expect_error(getCrossAssociation(mae, + experiment1 = 3, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 3, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + altexp1 = 1, + altexp2 = NULL, + method = "spearman", mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 3, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + altexp1 = FALSE, + altexp2 = NULL, + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 3, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + altexp2 = "test", + altexp1 = NULL, + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = TRUE, + assay.type1 = "counts", + assay.type2 = "counts", + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "test", + assay.type2 = "counts", + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = 1, + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = FALSE, + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = TRUE, + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = 1, + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 2, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = TRUE, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 0.1, + cor_threshold = 2, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 0.1, + cor_threshold = TRUE, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 0.1, + cor_threshold = NULL, + sort = 1, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 0.1, + cor_threshold = NULL, + sort = TRUE, + filter_self_correlations = 1, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 0.1, + cor_threshold = NULL, + sort = TRUE, + filter_self_correlations = TRUE, + verbose = 1)) + expect_error(getCrossAssociation(mae[[1]], + assay(mae1[[2]]), + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 0.1, + cor_threshold = NULL, + sort = TRUE, + filter_self_correlations = TRUE, + verbose = 1)) + expect_error(getCrossAssociation(mae[[1]], + NULL, + experiment1 = 1, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + method = "pearson", + mode = "matrix", + p_adj_method = "fdr", + p_adj_threshold = 0.1, + cor_threshold = NULL, + sort = TRUE, + filter_self_correlations = TRUE, + verbose = 1)) + expect_error(getCrossAssociation(mae, + experiment1 = 3, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + colData_variable1 = FALSE, + colData_variable2 = NULL, + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 3, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + colData_variable1 = NULL, + colData_variable2 = 1, + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + expect_error(getCrossAssociation(mae, + experiment1 = 3, + experiment2 = 2, + assay.type1 = "counts", + assay.type2 = "counts", + colData_variable1 = "test", + colData_variable2 = NULL, + method = "spearman", + mode = "table", + p_adj_method = "fdr", + p_adj_threshold = 0.05, + cor_threshold = NULL, + sort = FALSE, + filter_self_correlations = FALSE, + verbose = TRUE)) + ############################# Test input end ############################# + # Test that association is calculated correctly with numeric data + # Result from + # d1 <- t(assay(mae[[1]])) + # d2 <- t(assay(mae[[2]])) + # cc <- microbiome::associate(d1, d2, method='pearson', mode= "table") + cor_compare <- c(-0.67473156, 0.19980484, -0.19942102, -0.17346641, -0.16954081, + -0.15477367, -0.15279005, 0.08792788, 0.08443186) + p_adj_compare <- c(0.001247967, 0.784472862, 0.785332288, 0.830362548, + 0.836148425, 0.856762552, 0.859203260, 0.938444366, 0.942610008) + # Calculate correlation + cor <- getCrossAssociation(mae, + method = "pearson", + p_adj_threshold = NULL, + show_warnings = FALSE, + test_significance = TRUE) + # Take only specific taxa and lipids + df <- cor[cor$Var1 %in% c("Fusobacteria", "Campylobacter", "Actinomycetaceae") & + cor$Var2 %in% c("PE(48:7)", "TG(50:0)", "SM(d18:1/18:0)"), ] + # Sort the data, so that lowest p-values are first + df <- df[order(df$p_adj), ] + # Correlation values and p-values should be the same + expect_equal(round(df$cor, 7), round(cor_compare, 7)) + expect_equal(round(df$p_adj, 7), round(p_adj_compare, 7)) + + # Test that association is calculated correctly with factor data + # Create a dummy data + assay1 <- matrix(rep(c("A", "B", "B"), 20*30/3), + nrow = 20, ncol = 30) + assay2 <- matrix(rep(c("A", "B", "A", "A", "B", "B"), 20*30/6), + nrow = 20, ncol = 30) + # Reference + # ref <- c() + # for(i in 1:20){ + # ref <- c(ref, GoodmanKruskal::GKtau(assay1[i, ], assay2[i, ])$tauxy) + # } + ref <- c(0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, + 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00) + # Calculate values for 20 feature-pairs + result <- c() + for(i in 1:20){ + result <- c(result, .calculate_gktau(assay1[i, ], assay2[i, ])$estimate) + } + # Values should be the same + expect_equal(round(result, 4), round(ref, 4)) + + mae_sub <- mae[1:10, 1:10] + # Test that output is in correct type + expect_true( is.data.frame( + getCrossAssociation(mae_sub, p_adj_threshold = NULL, + show_warnings = FALSE, test_significance = TRUE)) ) + expect_true( is.data.frame(getCrossAssociation(mae_sub, + show_warnings = FALSE)) ) + # There should not be any p-values that are under 0 + expect_true( is.null( + getCrossAssociation(mae_sub, p_adj_threshold = 0, + show_warnings = FALSE, + test_significance = TRUE)) ) + # Test that output is in correct type + expect_true( is.list( + getCrossAssociation(mae_sub, mode = "matrix", + p_adj_threshold = NULL, + show_warnings = FALSE, + test_significance = TRUE)) ) + + expect_true( is.matrix(getCrossAssociation(mae_sub, + mode = "matrix", + show_warnings = FALSE)) ) + + # There should not be any p-values that are under 0 + expect_true( is.null( + getCrossAssociation(mae_sub, + p_adj_threshold = 0, + mode = "matrix", + show_warnings = FALSE, + test_significance = TRUE)) ) + + # When correlation between same assay is calculated, calculation is made faster + # by not calculating duplicates + expect_error(getCrossAssociation(mae, experiment1 = 1, experiment2 = 1, + show_warnings = FALSE, + symmetric = "TRUE", + test_significance = TRUE)) + expect_error(getCrossAssociation(mae, experiment1 = 1, experiment2 = 1, + show_warnings = FALSE, + symmetric = 1, + test_significance = TRUE)) + expect_error(getCrossAssociation(mae, experiment1 = 1, experiment2 = 1, + show_warnings = FALSE, + symmetric = NULL, + test_significance = TRUE)) + expect_error(getCrossAssociation(mae, experiment1 = 1, experiment2 = 1, + show_warnings = FALSE, + symmetric = c(TRUE, TRUE), + test_significance = TRUE)) + + time <- system.time( + cor <- getCrossAssociation(mae, experiment1 = 1, experiment2 = 1, + show_warnings = FALSE, + symmetric = TRUE, + test_significance = TRUE) + ) + time2 <- system.time( + cor2 <- getCrossAssociation(mae, experiment1 = 1, experiment2 = 1, + show_warnings = FALSE, + test_significance = TRUE) + ) + # Get random variables and test that their duplicates are equal + for(i in 1:10 ){ + random_var1 <- sample(cor$Var1, 1) + random_var2 <- sample(cor$Var1, 1) + expect_equal(as.numeric(cor[cor$Var1 == random_var1 & cor$Var2 == random_var2, c("cor", "pval", "p_adj")]), + as.numeric(cor[cor$Var1 == random_var2 & cor$Var2 == random_var1, c("cor", "pval", "p_adj")])) + } + expect_equal(cor, cor2) + # Test that symmetric = TRUE was faster + expect_true(time[3] < time2[3]) + # Test that paired samples work correctly + tse1 <- mae[[1]] + tse2 <- mae[[1]] + # Convert assay to have random values + mat <- matrix(sample(0:100, nrow(tse2)*ncol(tse2), replace = TRUE), + nrow = nrow(tse2), ncol = ncol(tse2)) + colnames(mat) <- colnames(tse2) + rownames(mat) <- rownames(tse2) + assay(tse2) <- mat + # Calculate with paired samples + cor_paired <- getCrossAssociation(tse1, + experiment2 = tse2, + paired = TRUE, + MARGIN = 2, + show_warnings = FALSE, + test_significance = TRUE) + # Calculate all pairs + cor <- getCrossAssociation(tse1, + experiment2 = tse2, + MARGIN = 2, + show_warnings = FALSE, + test_significance = TRUE) + # Take only pairs that are paired + cor <- cor[cor$Var1 == cor$Var2, ] + rownames(cor) <- NULL + + # Should be equal + expect_equal(cor[, c("cor", "pval")], + cor_paired[, c("cor", "pval")]) + + # Test that result does not depend on names (if there are equal names) + tse <- mae[[1]] + rownames(tse)[1:10] <- rep("Unknown", 10) + cor_table <- getCrossAssociation(tse, show_warnings = FALSE, + test_significance = TRUE) + cor_table_ref <- getCrossAssociation(mae[[1]], show_warnings = FALSE, + test_significance = TRUE) + expect_equal(cor_table[ , 3:5], cor_table_ref[ , 3:5]) + mat <- getCrossAssociation(tse, mode = "matrix", show_warnings = FALSE) + expect_true( is.matrix(mat) ) + expect_true(nrow(mat) == nrow(tse) && ncol(mat) == nrow(tse)) + mat <- getCrossAssociation(tse, mode = "matrix", show_warnings = FALSE, + cor_threshold = 0.8, + filter_self_correlation = TRUE) + expect_true(nrow(mat) < nrow(tse) && ncol(mat) < nrow(tse)) + + + # Test user's own function + expect_true( is.data.frame(getCrossAssociation(tse, method = "canberra", + mode = "table", + show_warnings = T, + association_FUN = stats::dist) ) ) + + expect_true( is.matrix( getCrossAssociation(tse, method = "bray", + show_warnings = FALSE, + mode = "matrix", + association_FUN = vegan::vegdist, + test_significance = TRUE) ) ) + expect_error( getCrossAssociation(tse, method = "bray", + show_warnings = FALSE, + mode = "matrix", + association_FUN = DelayedMatrixStats::rowSums2, + test_significance = TRUE) ) + + # Test that output has right columns + expect_equal(colnames(getCrossAssociation(tse, + show_warnings = FALSE)), + c("Var1", "Var2", "cor")) + expect_equal(colnames(getCrossAssociation(tse, show_warnings = FALSE, + test_significance = TRUE)), + c("Var1", "Var2", "cor", "pval", "p_adj")) + + # Test that the table have same information with different levels + tab1 <- getCrossAssociation(tse, show_warnings = FALSE) + tab1_levels1 <- levels(tab1$Var1) + tab1_levels2 <- levels(tab1$Var2) + tab1$Var1 <- as.character(tab1$Var1) + tab1$Var2 <- as.character(tab1$Var2) + tab2 <- getCrossAssociation(tse, show_warnings = FALSE, sort = TRUE) + tab2_levels1 <- levels(tab2$Var1) + tab2_levels2 <- levels(tab2$Var2) + tab2$Var1 <- as.character(tab2$Var1) + tab2$Var2 <- as.character(tab2$Var2) + expect_equal(tab1, tab2) + expect_true( !all(tab1_levels1 == tab2_levels1) ) + expect_true( !all(tab1_levels2 == tab2_levels2) ) + + # Test altexps + altExps(tse) <- splitByRanks(tse) + # Test that output has right columns + expect_equal(getCrossAssociation(tse, tse, show_warnings = FALSE, + altexp1 = 1, altexp2 = "Phylum"), + getCrossAssociation(altExps(tse)[[1]], altExp(tse, "Phylum"), + show_warnings = FALSE)) + expect_equal(getCrossAssociation(tse, tse, show_warnings = FALSE, + altexp1 = "Family", + altexp2 = NULL), + getCrossAssociation(altExp(tse, "Family"), tse, + show_warnings = FALSE)) + + # Test colData_variable + # Check that all the correct names are included + indices <- c("shannon", "gini_simpson") + tse <- estimateDiversity(tse, index = indices) + res <- getCrossAssociation(tse, tse, + assay.type1 = "counts", + colData_variable2 = indices) + unique_var1 <- unfactor(unique(res$Var1)) + unique_var2 <- unfactor(unique(res$Var2)) + rownames <- rownames(tse) + + expect_true( all(rownames %in% unique_var1) && all(unique_var1 %in% rownames) && + all(indices %in% unique_var2) && all(unique_var2 %in% indices) ) + # Check that assay.type is disabled + res2 <- getCrossAssociation(tse, assay.type1 = "counts", + assay.type2 = "counts", + colData_variable2 = indices) + expect_equal(res, res2) + + colData(tse)[, "test"] <- rep("a") + expect_error( + getCrossAssociation(tse, colData_variable2 = c("shannon", "test"))) + +} +}) diff --git a/tests/testthat/test-5getExperimentCrossAssociation.R b/tests/testthat/test-5getExperimentCrossAssociation.R deleted file mode 100644 index d8c46ea4e..000000000 --- a/tests/testthat/test-5getExperimentCrossAssociation.R +++ /dev/null @@ -1,540 +0,0 @@ - -context("getExperimentCrossAssociation") - -test_that("getExperimentCrossAssociation", { - - # Try 5 times to fetch the data - for(i in seq_len(5) ){ - mae <- tryCatch( - { - # Try to fetch the data - microbiomeDataSets::peerj32() - }, - error = function(cond) { - # If it was not possible to fetch the data, give FALSE - return(NULL) - } - ) - # Break if mae has the data - if( !is.null(mae) ){ - break - } - } - # Run tests if the data fetch was successful - if( !is.null(mae) ){ - ############################### Test input ############################### - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 3, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 3, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - altexp1 = 1, - altexp2 = NULL, - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 3, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - altexp1 = FALSE, - altexp2 = NULL, - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 3, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - altexp2 = "test", - altexp1 = NULL, - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = TRUE, - assay.type1 = "counts", - assay.type2 = "counts", - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "test", - assay.type2 = "counts", - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = 1, - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = FALSE, - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = TRUE, - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = 1, - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 2, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = TRUE, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 0.1, - cor_threshold = 2, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 0.1, - cor_threshold = TRUE, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 0.1, - cor_threshold = NULL, - sort = 1, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 0.1, - cor_threshold = NULL, - sort = TRUE, - filter_self_correlations = 1, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 0.1, - cor_threshold = NULL, - sort = TRUE, - filter_self_correlations = TRUE, - verbose = 1)) - expect_error(getExperimentCrossAssociation(mae[[1]], - assay(mae1[[2]]), - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 0.1, - cor_threshold = NULL, - sort = TRUE, - filter_self_correlations = TRUE, - verbose = 1)) - expect_error(getExperimentCrossAssociation(mae[[1]], - NULL, - experiment1 = 1, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - method = "pearson", - mode = "matrix", - p_adj_method = "fdr", - p_adj_threshold = 0.1, - cor_threshold = NULL, - sort = TRUE, - filter_self_correlations = TRUE, - verbose = 1)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 3, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - colData_variable1 = FALSE, - colData_variable2 = NULL, - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 3, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - colData_variable1 = NULL, - colData_variable2 = 1, - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - expect_error(getExperimentCrossAssociation(mae, - experiment1 = 3, - experiment2 = 2, - assay.type1 = "counts", - assay.type2 = "counts", - colData_variable1 = "test", - colData_variable2 = NULL, - method = "spearman", - mode = "table", - p_adj_method = "fdr", - p_adj_threshold = 0.05, - cor_threshold = NULL, - sort = FALSE, - filter_self_correlations = FALSE, - verbose = TRUE)) - ############################# Test input end ############################# - # Test that association is calculated correctly with numeric data - # Result from - # d1 <- t(assay(mae[[1]])) - # d2 <- t(assay(mae[[2]])) - # cc <- microbiome::associate(d1, d2, method='pearson', mode= "table") - cor_compare <- c(-0.67473156, 0.19980484, -0.19942102, -0.17346641, -0.16954081, - -0.15477367, -0.15279005, 0.08792788, 0.08443186) - p_adj_compare <- c(0.001247967, 0.784472862, 0.785332288, 0.830362548, - 0.836148425, 0.856762552, 0.859203260, 0.938444366, 0.942610008) - # Calculate correlation - cor <- testExperimentCrossAssociation(mae, method = "pearson", - p_adj_threshold = NULL, show_warnings = FALSE) - # Take only specific taxa and lipids - df <- cor[cor$Var1 %in% c("Fusobacteria", "Campylobacter", "Actinomycetaceae") & - cor$Var2 %in% c("PE(48:7)", "TG(50:0)", "SM(d18:1/18:0)"), ] - # Sort the data, so that lowest p-values are first - df <- df[order(df$p_adj), ] - # Correlation values and p-values should be the same - expect_equal(round(df$cor, 7), round(cor_compare, 7)) - expect_equal(round(df$p_adj, 7), round(p_adj_compare, 7)) - - # Test that association is calculated correctly with factor data - # Create a dummy data - assay1 <- matrix(rep(c("A", "B", "B"), 20*30/3), - nrow = 20, ncol = 30) - assay2 <- matrix(rep(c("A", "B", "A", "A", "B", "B"), 20*30/6), - nrow = 20, ncol = 30) - # Reference - # ref <- c() - # for(i in 1:20){ - # ref <- c(ref, GoodmanKruskal::GKtau(assay1[i, ], assay2[i, ])$tauxy) - # } - ref <- c(0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, - 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00, 0.25, 1.00) - # Calculate values for 20 feature-pairs - result <- c() - for(i in 1:20){ - result <- c(result, .calculate_gktau(assay1[i, ], assay2[i, ])$estimate) - } - # Values should be the same - expect_equal(round(result, 4), round(ref, 4)) - - mae_sub <- mae[1:10, 1:10] - # Test that output is in correct type - expect_true( is.data.frame( - testExperimentCrossAssociation(mae_sub, p_adj_threshold = NULL, show_warnings = FALSE)) ) - expect_true( is.data.frame( - getExperimentCrossAssociation(mae_sub, test_significance = TRUE, - p_adj_threshold = NULL, show_warnings = FALSE)) ) - expect_true( is.data.frame(getExperimentCrossAssociation(mae_sub, show_warnings = FALSE)) ) - # There should not be any p-values that are under 0 - expect_true( is.null( - testExperimentCrossAssociation(mae_sub, p_adj_threshold = 0, show_warnings = FALSE)) ) - # Test that output is in correct type - expect_true( is.list( - testExperimentCrossAssociation(mae_sub, mode = "matrix", - p_adj_threshold = NULL, show_warnings = FALSE)) ) - expect_true( is.list( - getExperimentCrossAssociation(mae_sub, test_significance = TRUE, - mode = "matrix", - p_adj_threshold = NULL, show_warnings = FALSE)) ) - expect_true( is.matrix(getExperimentCrossAssociation(mae_sub, mode = "matrix", - show_warnings = FALSE)) ) - - # There should not be any p-values that are under 0 - expect_true( is.null( - testExperimentCrossAssociation(mae_sub, p_adj_threshold = 0, mode = "matrix", show_warnings = FALSE)) ) - - # When correlation between same assay is calculated, calculation is made faster - # by not calculating duplicates - expect_error(testExperimentCrossAssociation(mae, experiment1 = 1, experiment2 = 1, - show_warnings = FALSE, symmetric = "TRUE")) - expect_error(testExperimentCrossAssociation(mae, experiment1 = 1, experiment2 = 1, - show_warnings = FALSE, symmetric = 1)) - expect_error(testExperimentCrossAssociation(mae, experiment1 = 1, experiment2 = 1, - show_warnings = FALSE, symmetric = NULL)) - expect_error(testExperimentCrossAssociation(mae, experiment1 = 1, experiment2 = 1, - show_warnings = FALSE, symmetric = c(TRUE, TRUE))) - - time <- system.time( - cor <- testExperimentCrossAssociation(mae, experiment1 = 1, experiment2 = 1, - show_warnings = FALSE, symmetric = TRUE) - ) - time2 <- system.time( - cor2 <- testExperimentCrossAssociation(mae, experiment1 = 1, experiment2 = 1, - show_warnings = FALSE) - ) - # Get random variables and test that their duplicates are equal - for(i in 1:10 ){ - random_var1 <- sample(cor$Var1, 1) - random_var2 <- sample(cor$Var1, 1) - expect_equal(as.numeric(cor[cor$Var1 == random_var1 & cor$Var2 == random_var2, c("cor", "pval", "p_adj")]), - as.numeric(cor[cor$Var1 == random_var2 & cor$Var2 == random_var1, c("cor", "pval", "p_adj")])) - } - expect_equal(cor, cor2) - # Test that symmetric = TRUE was faster - expect_true(time[3] < time2[3]) - # Test that paired samples work correctly - tse1 <- mae[[1]] - tse2 <- mae[[1]] - # Convert assay to have random values - mat <- matrix(sample(0:100, nrow(tse2)*ncol(tse2), replace = TRUE), - nrow = nrow(tse2), ncol = ncol(tse2)) - colnames(mat) <- colnames(tse2) - rownames(mat) <- rownames(tse2) - assay(tse2) <- mat - # Calculate with paired samples - cor_paired <- testExperimentCrossAssociation(tse1, - experiment2 = tse2, - paired = TRUE, - MARGIN = 2, - show_warnings = FALSE) - # Calculate all pairs - cor <- testExperimentCrossAssociation(tse1, - experiment2 = tse2, - MARGIN = 2, - show_warnings = FALSE) - # Take only pairs that are paired - cor <- cor[cor$Var1 == cor$Var2, ] - rownames(cor) <- NULL - - # Should be equal - expect_equal(cor[, c("cor", "pval")], - cor_paired[, c("cor", "pval")]) - - # Test that result does not depend on names (if there are equal names) - tse <- mae[[1]] - rownames(tse)[1:10] <- rep("Unknown", 10) - cor_table <- testExperimentCrossAssociation(tse, show_warnings = FALSE) - cor_table_ref <- testExperimentCrossAssociation(mae[[1]], show_warnings = FALSE) - expect_equal(cor_table[ , 3:5], cor_table_ref[ , 3:5]) - mat <- getExperimentCrossAssociation(tse, mode = "matrix", show_warnings = FALSE) - expect_true( is.matrix(mat) ) - expect_true(nrow(mat) == nrow(tse) && ncol(mat) == nrow(tse)) - mat <- getExperimentCrossAssociation(tse, mode = "matrix", show_warnings = FALSE, - cor_threshold = 0.8, filter_self_correlation = TRUE) - expect_true(nrow(mat) < nrow(tse) && ncol(mat) < nrow(tse)) - - - # Test user's own function - expect_true( is.data.frame(getExperimentCrossAssociation(tse, method = "canberra", - mode = "table", - show_warnings = T, - association_FUN = stats::dist) ) ) - - expect_true( is.matrix( testExperimentCrossAssociation(tse, method = "bray", - show_warnings = FALSE, - mode = "matrix", - association_FUN = vegan::vegdist) ) ) - expect_error( testExperimentCrossAssociation(tse, method = "bray", - show_warnings = FALSE, - mode = "matrix", - association_FUN = DelayedMatrixStats::rowSums2) ) - - # Test that output has right columns - expect_equal(colnames(getExperimentCrossAssociation(tse, show_warnings = FALSE)), c("Var1", "Var2", "cor")) - expect_equal(colnames(getExperimentCrossAssociation(tse, show_warnings = FALSE, - test_significance = TRUE)), - c("Var1", "Var2", "cor", "pval", "p_adj")) - - # Test that the table have same information with different levels - tab1 <- getExperimentCrossAssociation(tse, show_warnings = FALSE) - tab1_levels1 <- levels(tab1$Var1) - tab1_levels2 <- levels(tab1$Var2) - tab1$Var1 <- as.character(tab1$Var1) - tab1$Var2 <- as.character(tab1$Var2) - tab2 <- getExperimentCrossAssociation(tse, show_warnings = FALSE, sort = TRUE) - tab2_levels1 <- levels(tab2$Var1) - tab2_levels2 <- levels(tab2$Var2) - tab2$Var1 <- as.character(tab2$Var1) - tab2$Var2 <- as.character(tab2$Var2) - expect_equal(tab1, tab2) - expect_true( !all(tab1_levels1 == tab2_levels1) ) - expect_true( !all(tab1_levels2 == tab2_levels2) ) - - # Test altexps - altExps(tse) <- splitByRanks(tse) - # Test that output has right columns - expect_equal(getExperimentCrossAssociation(tse, tse, show_warnings = FALSE, - altexp1 = 1, altexp2 = "Phylum"), - getExperimentCrossAssociation(altExps(tse)[[1]], altExp(tse, "Phylum"), - show_warnings = FALSE)) - expect_equal(getExperimentCrossAssociation(tse, tse, show_warnings = FALSE, - altexp1 = "Family", altexp2 = NULL), - getExperimentCrossAssociation(altExp(tse, "Family"), tse, - show_warnings = FALSE)) - - # Test colData_variable - # Check that all the correct names are included - indices <- c("shannon", "gini_simpson") - tse <- estimateDiversity(tse, index = indices) - res <- getExperimentCrossAssociation(tse, tse, - assay.type1 = "counts", - colData_variable2 = indices) - unique_var1 <- unfactor(unique(res$Var1)) - unique_var2 <- unfactor(unique(res$Var2)) - rownames <- rownames(tse) - - expect_true( all(rownames %in% unique_var1) && all(unique_var1 %in% rownames) && - all(indices %in% unique_var2) && all(unique_var2 %in% indices) ) - # Check that assay.type is disabled - res2 <- getExperimentCrossAssociation(tse, - assay.type1 = "counts", - assay.type2 = "counts", - colData_variable2 = indices) - expect_equal(res, res2) - - colData(tse)[, "test"] <- rep("a") - expect_error( - getExperimentCrossAssociation(tse, - colData_variable2 = c("shannon", "test"))) - -} -})