diff --git a/DESCRIPTION b/DESCRIPTION index 9b0a8b7..f7112b8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,7 +7,7 @@ Authors@R: person("Jakob", "Russel", Description: Runs multiple differential abundance/expressions methods, with the same standardized input and output, to compare their results and test which are better at finding spiked features. -Depends: R (>= 3.5) +Depends: R (>= 3.4) biocViews: Sequencing, RNASeq, Microbiome, @@ -50,3 +50,4 @@ License: GPL (>= 3) | file LICENSE Encoding: UTF-8 LazyData: true RoxygenNote: 6.1.1 +BugReports: https://github.com/Russel88/DAtest/issues diff --git a/NAMESPACE b/NAMESPACE index 9729b18..87645cf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -64,7 +64,7 @@ export(norm_alr) export(norm_clr) export(powerDA) export(preDA) -export(prune.tests.DA) +export(pruneTests) export(runtimeDA) export(spikein) export(testDA) diff --git a/NEWS b/NEWS deleted file mode 100644 index 12210a6..0000000 --- a/NEWS +++ /dev/null @@ -1,4 +0,0 @@ -CHANGES IN VERSION 2.8.0 -------------------------- - - o Compliant with Bioconductor \ No newline at end of file diff --git a/R/DA.adx.R b/R/DA.adx.R index 5701a05..5e21b5f 100644 --- a/R/DA.adx.R +++ b/R/DA.adx.R @@ -5,6 +5,15 @@ #' @param predictor The predictor of interest. Factor, OR if \code{data} is a \code{phyloseq} object the name of the variable in \code{sample_data(data)} in quotation #' @param ... Additional arguments for the \code{aldex} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running ALDEx2 +#' res <- DA.adx(data = mat, predictor = pred) #' @export DA.adx <- function(data, predictor, ...){ diff --git a/R/DA.aoa.R b/R/DA.aoa.R index 20a17ca..2fc881b 100644 --- a/R/DA.aoa.R +++ b/R/DA.aoa.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{aov} function #' @param ... Additional arguments for the \code{aov} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running ANOVA on each feature +#' res <- DA.aoa(data = mat, predictor = pred) #' @export DA.aoa <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 1, allResults = FALSE, ...){ @@ -43,7 +52,7 @@ DA.aoa <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 1, all } # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") # ALR transformation diff --git a/R/DA.aoc.R b/R/DA.aoc.R index d1a4ba7..10b5f38 100644 --- a/R/DA.aoc.R +++ b/R/DA.aoc.R @@ -9,6 +9,15 @@ #' @param allResults If TRUE will return raw results from the \code{aov} function #' @param ... Additional arguments for the \code{aov} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running ANOVA on each feature +#' res <- DA.aoc(data = mat, predictor = pred) #' @export DA.aoc <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 1, allResults = FALSE, ...){ @@ -41,7 +50,7 @@ DA.aoc <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 1, all } # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") # ALR transformation diff --git a/R/DA.aov.R b/R/DA.aov.R index a544a1f..a2ca696 100644 --- a/R/DA.aov.R +++ b/R/DA.aov.R @@ -9,6 +9,15 @@ #' @param allResults If TRUE will return raw results from the \code{aov} function #' @param ... Additional arguments for the \code{aov} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running ANOVA on each feature +#' res <- DA.aov(data = mat, predictor = pred) #' @export DA.aov <- function(data, predictor, covars = NULL, relative = TRUE, p.adj = "fdr", allResults = FALSE, ...){ diff --git a/R/DA.bay.R b/R/DA.bay.R index 3dacdff..a53561a 100644 --- a/R/DA.bay.R +++ b/R/DA.bay.R @@ -6,6 +6,15 @@ #' @param allResults If TRUE will return raw results from the \code{getLikelihoods} function #' @param ... Additional arguments to the \code{getPriors.NB} and \code{getLikelihoods} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(100, size = 0.1, mu = 500), nrow = 50, ncol = 10) +#' rownames(mat) <- 1:50 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running baySeq +#' res <- DA.bay(data = mat, predictor = pred) #' @export DA.bay <- function(data, predictor, allResults = FALSE, ...){ diff --git a/R/DA.ds2.R b/R/DA.ds2.R index 043b6b8..3c93143 100644 --- a/R/DA.ds2.R +++ b/R/DA.ds2.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{DESeq} function #' @param ... Additional arguments for the \code{DESeq} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(200, size = 0.1, mu = 500), nrow = 20, ncol = 10) +#' rownames(mat) <- 1:20 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running DESeq2 +#' res <- DA.ds2(data = mat, predictor = pred) #' @export DA.ds2 <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", coeff = 2, coeff.ref = 1, allResults = FALSE, ...){ diff --git a/R/DA.ds2x.R b/R/DA.ds2x.R index b923090..8ce2316 100644 --- a/R/DA.ds2x.R +++ b/R/DA.ds2x.R @@ -12,6 +12,15 @@ #' @param allResults If TRUE will return raw results from the \code{DESeq} function #' @param ... Additional arguments for the \code{DESeq} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(6) +#' mat <- matrix(rnbinom(200, size = 0.2, mu = 500), nrow = 20, ncol = 10) +#' rownames(mat) <- 1:20 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running DESeq2 +#' res <- DA.ds2x(data = mat, predictor = pred) #' @export DA.ds2x <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", coeff = 2, coeff.ref = 1, allResults = FALSE, ...){ diff --git a/R/DA.ere.R b/R/DA.ere.R index fc8d574..0e0fe4f 100644 --- a/R/DA.ere.R +++ b/R/DA.ere.R @@ -6,6 +6,15 @@ #' @param p.adj Character. P-value adjustment. Default "fdr". See \code{p.adjust} for details #' @param ... Additional arguments for the \code{calcNormFactors}, \code{estimateCommonDisp}, \code{estimateTagwiseDisp} and \code{exactTest} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running edgeR +#' res <- DA.ere(data = mat, predictor = pred) #' @export DA.ere <- function(data, predictor, p.adj = "fdr", ...){ diff --git a/R/DA.ere2.R b/R/DA.ere2.R index 3597c11..33027c7 100644 --- a/R/DA.ere2.R +++ b/R/DA.ere2.R @@ -6,6 +6,15 @@ #' @param p.adj Character. P-value adjustment. Default "fdr". See \code{p.adjust} for details #' @param ... Additional arguments for the \code{calcNormFactors}, \code{estimateCommonDisp}, \code{estimateTagwiseDisp} and \code{exactTest} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running edgeR +#' res <- DA.ere2(data = mat, predictor = pred) #' @export DA.ere2 <- function(data, predictor, p.adj = "fdr", ...){ diff --git a/R/DA.erq.R b/R/DA.erq.R index 081ed23..608127a 100644 --- a/R/DA.erq.R +++ b/R/DA.erq.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{glmQLFTest} function #' @param ... Additional arguments for the \code{calcNormFactors}, \code{estimateDisp}, \code{glmQLFit} and \code{glmQLFTest} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running edgeR +#' res <- DA.erq(data = mat, predictor = pred) #' @export DA.erq <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", coeff = 2, allResults = FALSE, ...){ diff --git a/R/DA.erq2.R b/R/DA.erq2.R index 9e20d04..0fd4dc7 100644 --- a/R/DA.erq2.R +++ b/R/DA.erq2.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{glmQLFTest} function #' @param ... Additional arguments for the \code{calcNormFactors}, \code{estimateDisp}, \code{glmQLFit} and \code{glmQLFTest} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running edgeR +#' res <- DA.erq2(data = mat, predictor = pred) #' @export DA.erq2 <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", coeff = 2, allResults = FALSE, ...){ diff --git a/R/DA.fri.R b/R/DA.fri.R index 8488117..9c9dd7f 100644 --- a/R/DA.fri.R +++ b/R/DA.fri.R @@ -9,6 +9,16 @@ #' @param allResults If TRUE will return raw results from the \code{friedman.test} function #' @param ... Additional arguments for the \code{friedman.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table, predictor, and paired variable +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' subject <- rep(1:5, 3) +#' +#' # Running Friedman test on each feature +#' res <- DA.fri(data = mat, predictor = pred, paired = subject) #' @export DA.fri <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr", allResults = FALSE, ...){ @@ -44,9 +54,9 @@ DA.fri <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr if(allResults){ return(reslist) } else { - res <- data.frame(statistic = vapply(reslist, function(x) x$statistic), - parameter = vapply(reslist, function(x) x$parameter), - pval = vapply(reslist, function(x) x$p.value)) + res <- data.frame(statistic = sapply(reslist, function(x) x$statistic), + parameter = sapply(reslist, function(x) x$parameter), + pval = sapply(reslist, function(x) x$p.value)) res$pval.adj <- p.adjust(res$pval, method = p.adj) res$Feature <- gsub(".Friedman.*","",rownames(res)) diff --git a/R/DA.kru.R b/R/DA.kru.R index a702f05..42201ab 100644 --- a/R/DA.kru.R +++ b/R/DA.kru.R @@ -8,6 +8,15 @@ #' @param allResults If TRUE will return raw results from the \code{kruskal.test} function #' @param ... Additional arguments for the \code{kruskal.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running Kruskal-Wallis on each feature +#' res <- DA.kru(data = mat, predictor = pred) #' @export DA.kru <- function(data, predictor, relative = TRUE, p.adj = "fdr", allResults = FALSE, ...){ @@ -41,7 +50,7 @@ DA.kru <- function(data, predictor, relative = TRUE, p.adj = "fdr", allResults = if(allResults){ return(tests) } else { - res <- data.frame(pval = vapply(tests, function(x) x$p.value)) + res <- data.frame(pval = sapply(tests, function(x) x$p.value)) res$pval.adj <- p.adjust(res$pval, method = p.adj) res$Feature <- rownames(res) res$Method <- "Kruskal-Wallis (kru)" diff --git a/R/DA.lao.R b/R/DA.lao.R index 4e39e26..8e60251 100644 --- a/R/DA.lao.R +++ b/R/DA.lao.R @@ -10,6 +10,15 @@ #' @param allResults If TRUE will return raw results from the \code{aov} function #' @param ... Additional arguments for the \code{aov} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running ANOVA on each feature +#' res <- DA.lao(data = mat, predictor = pred) #' @export DA.lao <- function(data, predictor, covars = NULL, relative = TRUE, p.adj = "fdr", delta = 1, allResults = FALSE, ...){ diff --git a/R/DA.lao2.R b/R/DA.lao2.R index 5b87164..698e282 100644 --- a/R/DA.lao2.R +++ b/R/DA.lao2.R @@ -9,6 +9,15 @@ #' @param allResults If TRUE will return raw results from the \code{aov} function #' @param ... Additional arguments for the \code{aov} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running ANOVA on each feature +#' res <- DA.lao2(data = mat, predictor = pred) #' @export DA.lao2 <- function(data, predictor, covars = NULL, p.adj = "fdr", delta = 0.001, allResults = FALSE, ...){ diff --git a/R/DA.lia.R b/R/DA.lia.R index e96ecc0..cd08260 100644 --- a/R/DA.lia.R +++ b/R/DA.lia.R @@ -12,6 +12,15 @@ #' @param allResults If TRUE will return raw results from the \code{eBayes} function #' @param ... Additional arguments for the \code{eBayes} and \code{lmFit} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running limma +#' res <- DA.lia(data = mat, predictor = pred) #' @export DA.lia <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", delta = 1, coeff = 2, allResults = FALSE, ...){ @@ -46,7 +55,7 @@ DA.lia <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL } # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") # ALR transformation diff --git a/R/DA.lic.R b/R/DA.lic.R index 2ec34ee..0cb1447 100644 --- a/R/DA.lic.R +++ b/R/DA.lic.R @@ -12,6 +12,15 @@ #' @param allResults If TRUE will return raw results from the \code{eBayes} function #' @param ... Additional arguments for the \code{eBayes} and \code{lmFit} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running limma +#' res <- DA.lic(data = mat, predictor = pred) #' @export DA.lic <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", delta = 1, coeff = 2, allResults = FALSE, ...){ @@ -46,7 +55,7 @@ DA.lic <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL } # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") # ALR transformation diff --git a/R/DA.lim.R b/R/DA.lim.R index 32f3f6a..213d2b6 100644 --- a/R/DA.lim.R +++ b/R/DA.lim.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{eBayes} function #' @param ... Additional arguments for the \code{eBayes} and \code{lmFit} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running limma +#' res <- DA.lim(data = mat, predictor = pred) #' @export DA.lim <- function(data, predictor, paired = NULL, covars = NULL, relative = TRUE, out.all = NULL, p.adj = "fdr", coeff = 2, allResults = FALSE, ...){ diff --git a/R/DA.lli.R b/R/DA.lli.R index 5999e14..67f363f 100644 --- a/R/DA.lli.R +++ b/R/DA.lli.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{eBayes} function #' @param ... Additional arguments for the \code{eBayes} and \code{lmFit} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running limma +#' res <- DA.lli(data = mat, predictor = pred) #' @export DA.lli <- function(data, predictor, paired = NULL, covars = NULL, relative = TRUE, out.all = NULL, p.adj = "fdr", delta = 1, coeff = 2, allResults = FALSE, ...){ diff --git a/R/DA.lli2.R b/R/DA.lli2.R index e55bda8..b9d5d64 100644 --- a/R/DA.lli2.R +++ b/R/DA.lli2.R @@ -12,6 +12,15 @@ #' @param allResults If TRUE will return raw results from the \code{eBayes} function #' @param ... Additional arguments for the \code{eBayes} and \code{lmFit} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running limma +#' res <- DA.lli2(data = mat, predictor = pred) #' @export DA.lli2 <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", delta = 0.001, coeff = 2, allResults = FALSE, ...){ diff --git a/R/DA.llm.R b/R/DA.llm.R index 57aa079..ddc0948 100644 --- a/R/DA.llm.R +++ b/R/DA.llm.R @@ -14,6 +14,15 @@ #' @param allResults If TRUE will return raw results from the \code{lm}/\code{lme} function #' @param ... Additional arguments for the \code{lm}/\code{lme} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running linear model on each feature +#' res <- DA.llm(data = mat, predictor = pred) #' @import nlme #' @export diff --git a/R/DA.llm2.R b/R/DA.llm2.R index f5e5412..04ac672 100644 --- a/R/DA.llm2.R +++ b/R/DA.llm2.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{lm}/\code{lme} function #' @param ... Additional arguments for the \code{lm}/\code{lme} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running linear model on each feature +#' res <- DA.llm2(data = mat, predictor = pred) #' @import nlme #' @export diff --git a/R/DA.lma.R b/R/DA.lma.R index 66f8566..7d3ccc6 100644 --- a/R/DA.lma.R +++ b/R/DA.lma.R @@ -14,6 +14,15 @@ #' @param allResults If TRUE will return raw results from the \code{lm}/\code{lme} function #' @param ... Additional arguments for the \code{lm}/\code{lme} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running linear model on each feature +#' res <- DA.lma(data = mat, predictor = pred) #' @import nlme #' @export @@ -44,7 +53,7 @@ DA.lma <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL count_table <- as.data.frame.matrix(count_table) # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") # ALR transformation diff --git a/R/DA.lmc.R b/R/DA.lmc.R index 365386f..48d19db 100644 --- a/R/DA.lmc.R +++ b/R/DA.lmc.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{lm}/\code{lme} function #' @param ... Additional arguments for the \code{lm}/\code{lme} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running linear model on each feature +#' res <- DA.lmc(data = mat, predictor = pred) #' @import nlme #' @export @@ -43,10 +52,10 @@ DA.lmc <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL count_table <- as.data.frame.matrix(count_table) # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") - # ALR transformation + # CLR transformation count_table <- norm_clr(count_table) # Define design diff --git a/R/DA.lrm.R b/R/DA.lrm.R index f5b48c0..d751a49 100644 --- a/R/DA.lrm.R +++ b/R/DA.lrm.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{lm}/\code{lme} function #' @param ... Additional arguments for the \code{lm}/\code{lme} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running linear model on each feature +#' res <- DA.lrm(data = mat, predictor = pred) #' @import nlme #' @export diff --git a/R/DA.ltt.R b/R/DA.ltt.R index c505ff2..6f60b26 100644 --- a/R/DA.ltt.R +++ b/R/DA.ltt.R @@ -12,6 +12,15 @@ #' @param allResults If TRUE will return raw results from the \code{t.test} function #' @param ... Additional arguments for the \code{t.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running t-test on each feature +#' res <- DA.ltt(data = mat, predictor = pred) #' @export DA.ltt <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr", delta = 1, testStat = function(case,control){log2((mean(case)+0.001)/(mean(control)+0.001))}, testStat.pair = function(case,control){log2(mean((case+0.001)/(control+0.001)))},allResults = FALSE, ...){ diff --git a/R/DA.ltt2.R b/R/DA.ltt2.R index b33c046..87d72ff 100644 --- a/R/DA.ltt2.R +++ b/R/DA.ltt2.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{t.test} function #' @param ... Additional arguments for the \code{t.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running t-test on each feature +#' res <- DA.ltt2(data = mat, predictor = pred) #' @export DA.ltt2 <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 0.001, testStat = function(case,control){log2((mean(exp(case)))/(mean(exp(control))))}, testStat.pair = function(case,control){log2(mean((exp(case))/(exp(control))))},allResults = FALSE, ...){ diff --git a/R/DA.msf.R b/R/DA.msf.R index 3e34aaf..0436536 100644 --- a/R/DA.msf.R +++ b/R/DA.msf.R @@ -8,6 +8,15 @@ #' @param allResults If TRUE will return raw results from the \code{fitFeatureModel} function #' @param ... Additional arguments for the \code{fitFeatureModel} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running MetagenomeSeq feature model +#' res <- DA.msf(data = mat, predictor = pred) #' @export DA.msf <- function(data, predictor, p.adj = "fdr", allResults = FALSE, ...){ diff --git a/R/DA.mva.R b/R/DA.mva.R index e9a5fb7..99fb6bd 100644 --- a/R/DA.mva.R +++ b/R/DA.mva.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{mvabund} function #' @param ... Additional arguments for the \code{manyglm} and \code{summary.manyglm} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(100, size = 0.1, mu = 500), nrow = 10, ncol = 10) +#' rownames(mat) <- 1:10 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running mvabund manyglm +#' res <- DA.mva(data = mat, predictor = pred) #' @export DA.mva <- function(data, predictor, paired = NULL, covars = NULL, relative = TRUE, p.adj = "fdr", coeff = 2, coeff.ref = 1, resamp = "montecarlo", allResults = FALSE, ...){ diff --git a/R/DA.neb.R b/R/DA.neb.R index 6d10236..21394df 100644 --- a/R/DA.neb.R +++ b/R/DA.neb.R @@ -15,6 +15,15 @@ #' @param allResults If TRUE will return raw results from the \code{glm.nb}/\code{glmer.nb} function #' @param ... Additional arguments for the \code{glm.nb}/\code{glmer.nb} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running Negative Binomial regression on each feature +#' res <- DA.neb(data = mat, predictor = pred) #' @import MASS #' @importFrom lme4 glmer.nb glmer #' @export @@ -198,7 +207,7 @@ DA.neb <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU rownames(res) <- rownames(count_table) } res$pval.adj <- p.adjust(res$pval, method = p.adj) - res$Feature <- rownames(res) + res$Feature <- rownames(count_table) res$Method <- "Negbinom GLM (neb)" if(nrow(res) > 1){ diff --git a/R/DA.pea.R b/R/DA.pea.R index 330480c..fad496e 100644 --- a/R/DA.pea.R +++ b/R/DA.pea.R @@ -7,6 +7,15 @@ #' @param p.adj Character. P-value adjustment. Default "fdr". See \code{p.adjust} for details #' @param ... Additional arguments for the \code{cor.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 50, ncol = 20) +#' rownames(mat) <- 1:50 +#' pred <- rnorm(20) +#' +#' # Running Pearson correlation on each feature +#' res <- DA.pea(data = mat, predictor = pred) #' @export DA.pea <- function(data, predictor, relative = TRUE, p.adj = "fdr", ...){ @@ -36,11 +45,11 @@ DA.pea <- function(data, predictor, relative = TRUE, p.adj = "fdr", ...){ peas <- apply(count.rel, 1, pea) # Collect results - res <- data.frame(pval = vapply(peas, function(x) x$p.value)) + res <- data.frame(pval = sapply(peas, function(x) x$p.value)) res$pval.adj <- p.adjust(res$pval, method = p.adj) - res$cor <- vapply(peas, function(x) x$estimate) - res$`cor_2.5%` <- vapply(peas, function(x) x$conf.int[1]) - res$`cor_97.5%` <- vapply(peas, function(x) x$conf.int[2]) + res$cor <- sapply(peas, function(x) x$estimate) + res$`cor_2.5%` <- sapply(peas, function(x) x$conf.int[1]) + res$`cor_97.5%` <- sapply(peas, function(x) x$conf.int[2]) res$Feature <- rownames(res) res$Method <- "Pearson (pea)" diff --git a/R/DA.per.R b/R/DA.per.R index 7025f04..1cdb8ac 100644 --- a/R/DA.per.R +++ b/R/DA.per.R @@ -18,6 +18,15 @@ #' @param noOfIterations Integer. Iterations for permutations. Default 10000 #' @param margin Numeric. Margin for when to stop iterations if p-value is high and unlikely to become low #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running permutation test on each feature +#' res <- DA.per(data = mat, predictor = pred) #' @export DA.per <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr", testStat = function(case,control){log2((mean(case)+1)/(mean(control)+1))}, testStat.pair = function(case,control){log2(mean((case+1)/(control+1)))}, noOfIterations = 10000, margin = 50){ diff --git a/R/DA.phyloseq.R b/R/DA.phyloseq.R index ab70749..c9d4062 100644 --- a/R/DA.phyloseq.R +++ b/R/DA.phyloseq.R @@ -1,5 +1,7 @@ #' Extract data from a \code{phyloseq} object to be used in \code{DAtest} #' +#' Internal function +#' #' @param data A \code{phyloseq} object #' @param predictor The \code{predictor} of interest #' @param paired Factor for paired/blocked experimental designs. diff --git a/R/DA.poi.R b/R/DA.poi.R index 7136fae..72f4767 100644 --- a/R/DA.poi.R +++ b/R/DA.poi.R @@ -15,6 +15,15 @@ #' @param allResults If TRUE will return raw results from the \code{glm}/\code{glmer} function #' @param ... Additional arguments for the \code{glm}/\code{glmer} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running Poisson regression on each feature +#' res <- DA.poi(data = mat, predictor = pred) #' @importFrom lme4 glmer #' @export diff --git a/R/DA.qpo.R b/R/DA.qpo.R index bd9bf24..7054e6e 100644 --- a/R/DA.qpo.R +++ b/R/DA.qpo.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{glm} function #' @param ... Additional arguments for the \code{glm} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running Quasi-Poisson regression on each feature +#' res <- DA.qpo(data = mat, predictor = pred) #' @export DA.qpo <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NULL, p.adj = "fdr", coeff = 2, coeff.ref = 1, allResults = FALSE, ...){ diff --git a/R/DA.qua.R b/R/DA.qua.R index 9ab5858..944a306 100644 --- a/R/DA.qua.R +++ b/R/DA.qua.R @@ -9,6 +9,16 @@ #' @param allResults If TRUE will return raw results from the \code{quade.test} function #' @param ... Additional arguments for the \code{quade.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table, predictor, and paired variable +#' set.seed(4) +#' mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' subject <- rep(1:5, 3) +#' +#' # Running Quade test on each feature +#' res <- DA.qua(data = mat, predictor = pred, paired = subject) #' @export DA.qua <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr", allResults = FALSE, ...){ @@ -44,10 +54,10 @@ DA.qua <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr if(allResults){ return(reslist) } else { - res <- data.frame(statistic = vapply(reslist, function(x) x$statistic), - num.df = vapply(reslist, function(x) x$parameter[1]), - denom.df = vapply(reslist, function(x) x$parameter[2]), - pval = vapply(reslist, function(x) x$p.value)) + res <- data.frame(statistic = sapply(reslist, function(x) x$statistic), + num.df = sapply(reslist, function(x) x$parameter[1]), + denom.df = sapply(reslist, function(x) x$parameter[2]), + pval = sapply(reslist, function(x) x$p.value)) res[is.nan(res$statistic),"pval"] <- 1 res$pval.adj <- p.adjust(res$pval, method = p.adj) diff --git a/R/DA.sam.R b/R/DA.sam.R index b8aed16..ba5f7c5 100644 --- a/R/DA.sam.R +++ b/R/DA.sam.R @@ -8,6 +8,15 @@ #' @param allResults If TRUE will return raw results from the \code{SAMseq} function #' @param ... Additional arguments for the \code{SAMseq} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running SamSeq +#' res <- DA.sam(data = mat, predictor = pred) #' @export DA.sam <- function(data, predictor, paired = NULL, fdr.output = 0.05, allResults = FALSE, ...){ diff --git a/R/DA.spe.R b/R/DA.spe.R index 1f84861..e1e0137 100644 --- a/R/DA.spe.R +++ b/R/DA.spe.R @@ -7,6 +7,15 @@ #' @param p.adj Character. P-value adjustment. Default "fdr". See \code{p.adjust} for details #' @param ... Additional arguments for the \code{cor.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 50, ncol = 20) +#' rownames(mat) <- 1:50 +#' pred <- rnorm(20) +#' +#' # Running Spearman correlation on each feature +#' res <- DA.spe(data = mat, predictor = pred) #' @export DA.spe <- function(data, predictor, relative = TRUE, p.adj = "fdr", ...){ @@ -36,9 +45,9 @@ DA.spe <- function(data, predictor, relative = TRUE, p.adj = "fdr", ...){ spes <- apply(count.rel, 1, spe) # Collect results - res <- data.frame(pval = vapply(spes, function(x) x$p.value)) + res <- data.frame(pval = sapply(spes, function(x) x$p.value)) res$pval.adj <- p.adjust(res$pval, method = p.adj) - res$rho <- vapply(spes, function(x) x$estimate) + res$rho <- sapply(spes, function(x) x$estimate) res$Feature <- rownames(res) res$Method <- "Spearman (spe)" diff --git a/R/DA.tta.R b/R/DA.tta.R index 4646317..4a25efc 100644 --- a/R/DA.tta.R +++ b/R/DA.tta.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{t.test} function #' @param ... Additional arguments for the \code{t.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running t-test on each feature +#' res <- DA.tta(data = mat, predictor = pred) #' @export DA.tta <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 1, testStat = function(case,control){mean(case)-mean(control)}, testStat.pair = function(case,control){mean(case-control)}, allResults = FALSE, ...){ @@ -43,7 +52,7 @@ DA.tta <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 1, tes } # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") # ALR transformation diff --git a/R/DA.ttc.R b/R/DA.ttc.R index 3792d1e..39d6ca5 100644 --- a/R/DA.ttc.R +++ b/R/DA.ttc.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{t.test} function #' @param ... Additional arguments for the \code{t.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running t-test on each feature +#' res <- DA.ttc(data = mat, predictor = pred) #' @export DA.ttc <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 1, testStat = function(case,control){mean(case)-mean(control)}, testStat.pair = function(case,control){mean(case-control)}, allResults = FALSE, ...){ @@ -41,7 +50,7 @@ DA.ttc <- function(data, predictor, paired = NULL, p.adj = "fdr", delta = 1, tes } # Zero-correction - count_table <- apply(count_table, 2, function(y) vapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) + count_table <- apply(count_table, 2, function(y) sapply(y,function(x) ifelse(x==0,delta,(1-(sum(y==0)*delta)/sum(y))*x))) if(any(count_table <= 0)) stop("Zero-correction failed. Dataset likely contains too many zeroes") # CLR transformation diff --git a/R/DA.ttt.R b/R/DA.ttt.R index 9653626..e925091 100644 --- a/R/DA.ttt.R +++ b/R/DA.ttt.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{t.test} function #' @param ... Additional arguments for the \code{t.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running t-test +#' res <- DA.ttt(data = mat, predictor = pred) #' @export DA.ttt <- function(data, predictor,paired = NULL, relative = TRUE, p.adj = "fdr", testStat = function(case,control){log2((mean(case)+0.001)/(mean(control)+0.001))}, testStat.pair = function(case,control){log2(mean((case+0.001)/(control+0.001)))},allResults = FALSE, ...){ diff --git a/R/DA.vli.R b/R/DA.vli.R index 0c5279c..72f1e30 100644 --- a/R/DA.vli.R +++ b/R/DA.vli.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{eBayes} function #' @param ... Additional arguments for the \code{voom}, \code{eBayes} and \code{lmFit} functions #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running limma +#' res <- DA.vli(data = mat, predictor = pred) #' @export DA.vli <- function(data, predictor, paired = NULL, covars = NULL, out.all = NULL, p.adj = "fdr", coeff = 2, allResults = FALSE, ...){ diff --git a/R/DA.wil.R b/R/DA.wil.R index 6877f1f..9832f58 100644 --- a/R/DA.wil.R +++ b/R/DA.wil.R @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{wilcox.test} function #' @param ... Additional arguments for the \code{wilcox.test} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running Wilcoxon test on each feature +#' res <- DA.wil(data = mat, predictor = pred) #' @export DA.wil <- function(data, predictor, paired = NULL, relative = TRUE, p.adj = "fdr", testStat = function(case,control){log2((mean(case)+0.001)/(mean(control)+0.001))}, testStat.pair = function(case,control){log2(mean((case+0.001)/(control+0.001)))}, allResults = FALSE, ...){ diff --git a/R/DA.zig.R b/R/DA.zig.R index e4039b3..4da9a50 100644 --- a/R/DA.zig.R +++ b/R/DA.zig.R @@ -1,4 +1,4 @@ -#' MetageonomeSeq ZIG +#' MetagenomeSeq ZIG #' #' Implementation of Metagenome zero-inflated gaussian model for \code{DAtest} #' @param data Either a matrix with counts/abundances, OR a \code{phyloseq} object. If a matrix/data.frame is provided rows should be taxa/genes/proteins and columns samples @@ -11,6 +11,15 @@ #' @param allResults If TRUE will return raw results from the \code{fitZig} function #' @param ... Additional arguments for the \code{fitZig} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running MetagenomeSeq Zero-inflated Gaussian +#' res <- DA.zig(data = mat, predictor = pred) #' @export DA.zig <- function(data, predictor, paired = NULL, covars = NULL, p.adj = "fdr", by = 2, eff = 0.5, allResults = FALSE, ...){ diff --git a/R/DA.znb.R b/R/DA.znb.R index 6de3f16..a1d5758 100644 --- a/R/DA.znb.R +++ b/R/DA.znb.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{zeroinfl} function #' @param ... Additional arguments for the \code{zeroinfl} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running Zero-inflated Negative Binomial regression on each feature +#' res <- DA.znb(data = mat, predictor = pred) #' @export DA.znb <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NULL, p.adj = "fdr", coeff = 2, coeff.ref = 1, allResults = FALSE, ...){ @@ -154,7 +163,7 @@ DA.znb <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NU } res$pval.adj <- p.adjust(res$pval, method = p.adj) - res$Feature <- rownames(res) + res$Feature <- rownames(count_table) res$Method <- "ZI-NegBin GLM (znb)" if(nrow(res) > 1){ diff --git a/R/DA.zpo.R b/R/DA.zpo.R index e30be8f..146dd24 100644 --- a/R/DA.zpo.R +++ b/R/DA.zpo.R @@ -13,6 +13,15 @@ #' @param allResults If TRUE will return raw results from the \code{zeroinfl} function #' @param ... Additional arguments for the \code{zeroinfl} function #' @return A data.frame with with results. +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running Zero-inflated Poisson regression on each feature +#' res <- DA.zpo(data = mat, predictor = pred) #' @export DA.zpo <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NULL, p.adj = "fdr", coeff = 2, coeff.ref = 1, allResults = FALSE, ...){ @@ -153,7 +162,7 @@ DA.zpo <- function(data, predictor, covars = NULL, relative = TRUE, out.all = NU rownames(res) <- rownames(count_table) } res$pval.adj <- p.adjust(res$pval, method = p.adj) - res$Feature <- rownames(res) + res$Feature <- rownames(count_table) res$Method <- "ZI-Poisson GLM (zpo)" if(nrow(res) > 1){ diff --git a/R/DA.zzz.R b/R/DA.zzz.R index f01a6b4..6927e1b 100644 --- a/R/DA.zzz.R +++ b/R/DA.zzz.R @@ -8,6 +8,39 @@ #' @param p.adj Character. P-value adjustment. Default "fdr". See \code{p.adjust} for details #' @param FUN Function to apply to data. Should take input in the following order: \code{count_table} (data.frame, samples are columns), \code{predictor} (vector), \code{paired} (factor), \code{covars} (named list with vector). Output should be a dataframe with at least the following columns: \code{Feature}, \code{pval} and \code{Method}. #' @return A data.frame with results from the user-defined method +#' @examples +#' # Creating random count_table and predictor +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Define function for t-test +#' myfun <- function(count_table, predictor, paired, covars){ +#' +#' # Relative abundance +#' rel <- apply(count_table, 2, function(x) x/sum(x)) +#' +#' # t-test function +#' # Wrapping this function in tryCatch(..., error = function(e){NA}) +#' # ensures that our main function won't fail if t.test fails on some features +#' tfun <- function(x){ +#' tryCatch(t.test(x ~ predictor)$p.value, error = function(e){NA}) +#' } +#' +#' # P-values for each feature +#' pvals <- apply(rel, 1, tfun) +#' +#' # Collect and return data +#' df <- data.frame(Feature = rownames(count_table), +#' pval = pvals) +#' df$pval.adj <- p.adjust(df$pval, method = "fdr") +#' df$Method <- "My own t-test" +#' return(df) +#' } +#' +#' # Running the test +#' res <- DA.zzz(data = mat, predictor = pred, FUN = myfun) #' @export DA.zzz <- function(data, predictor, paired = NULL, covars = NULL, p.adj = "fdr", FUN = NULL){ diff --git a/R/add.tax.DA.R b/R/add.tax.DA.R index f6767b0..320a5d5 100644 --- a/R/add.tax.DA.R +++ b/R/add.tax.DA.R @@ -3,9 +3,8 @@ #' Internal function #' @param data phyloseq object #' @param res data.frame with results -#' @return A \code{res} merged with the tax.table from the phyloseq object +#' @return \code{res} merged with the tax.table from the phyloseq object #' @export - add.tax.DA <- function(data, res){ loadNamespace("phyloseq") diff --git a/R/allDA.R b/R/allDA.R index 8d6fe2f..fa825fe 100644 --- a/R/allDA.R +++ b/R/allDA.R @@ -24,7 +24,37 @@ #' \item details - A dataframe with details from the run #' \item results - A complete list of output from all the methods. Example: Get wilcoxon results from 2. run as such: \code{$results[[2]]["wil"]} #' } +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(500, size = 0.1, mu = 500), nrow = 50, ncol = 10) +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) #' +#' # Running allDA to compare methods +#' # This example uses 1 core (cores = 1). +#' # Remove the cores argument to get it as high (and thereby fast) as possible. +#' res <- allDA(data = mat, predictor = pred, cores = 1) +#' +#' # View adjusted p-values from all methods +#' print(res$adj) +#' +#' # View estimates from all methods +#' print(res$est) +#' +#' \donttest{ +#' # Include a paired variable for dependent/blocked samples +#' subject <- rep(1:5, 2) +#' res <- allDA(data = mat, predictor = pred, paired = subject) +#' +#' # Include covariates +#' covar1 <- rnorm(10) +#' covar2 <- rep(c("A","B"), 5) +#' res <- allDA(data = mat, predictor = pred, +#' covars = list(FirstCovar = covar1, CallItWhatYouWant = covar2)) +#' +#' # Data is absolute abundance +#' res <- allDA(data = mat, predictor = pred, relative = FALSE) +#' } #' @export allDA <- function(data, predictor, paired = NULL, covars = NULL, @@ -79,7 +109,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL, if(!isTRUE(all(unlist(count_table) == floor(unlist(count_table))))) decimal <- TRUE if(any(count_table == 0)) zeroes <- TRUE tests <- unique(tests) - if(!"zzz" %in% tests) tests <- prune.tests.DA(tests, predictor, paired, covars, relative, decimal, zeroes) + if(!"zzz" %in% tests) tests <- pruneTests(tests, predictor, paired, covars, relative, decimal, zeroes) if(length(tests) == 0) stop("No tests to run!") if(verbose) message(paste("Running on",cores,"cores")) @@ -211,7 +241,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL, names(results) <- tests # Handle failed tests - results <- results[!vapply(results,is.null)] + results <- results[!sapply(results,is.null)] if(length(names(results)) != length(tests)){ if(length(tests) - length(names(results)) == 1){ if(verbose) message(paste(paste(tests[!tests %in% names(results)],collapse = ", "),"was excluded due to failure")) @@ -250,7 +280,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL, # Raw p-values Pval.raw <- lapply(results,function(x) tryCatch(as.data.frame(x[,c("Feature","pval")]), error = function(e) NULL)) - Pval.raw <- Pval.raw[!vapply(Pval.raw,is.null)] + Pval.raw <- Pval.raw[!sapply(Pval.raw,is.null)] if(length(Pval.raw) > 0){ df.raw <- suppressWarnings(Reduce(function(x,y) merge(x, y, by= "Feature", all.x = TRUE, all.y = TRUE), Pval.raw)) colnames(df.raw)[2:ncol(df.raw)] <- names(Pval.raw) @@ -261,7 +291,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL, # Adjusted p-values Pval.adj <- lapply(results,function(x) tryCatch(as.data.frame(x[,c("Feature","pval.adj")]), error = function(e) NULL)) - Pval.adj <- Pval.adj[!vapply(Pval.adj,is.null)] + Pval.adj <- Pval.adj[!sapply(Pval.adj,is.null)] if(length(Pval.adj) > 0){ df.adj <- suppressWarnings(Reduce(function(x,y) merge(x, y, by= "Feature", all.x = TRUE, all.y = TRUE), Pval.adj)) colnames(df.adj)[2:ncol(df.adj)] <- names(Pval.adj) @@ -331,7 +361,7 @@ allDA <- function(data, predictor, paired = NULL, covars = NULL, } else return(NULL) } else return(NULL) } - list.est <- list.est[!vapply(list.est, is.null)] + list.est <- list.est[!sapply(list.est, is.null)] if(length(list.est) > 0){ df.est <- Reduce(function(x,y) merge(x, y, by= "Feature", all.x = TRUE, all.y = TRUE), list.est) if(class(data) == "phyloseq") df.est <- add.tax.DA(data, df.est) diff --git a/R/featurePlot.R b/R/featurePlot.R index 5a03dad..b9f4a4d 100644 --- a/R/featurePlot.R +++ b/R/featurePlot.R @@ -15,8 +15,15 @@ #' @param delta Pseudocount for log10 normalization #' @param covar.quant Quantiles for cutting quantitative \code{covars} #' @return A ggplot +#' @examples +#' # Create random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(500, size = 0.1, mu = 500), nrow = 50, ncol = 10) +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' rownames(mat) <- 1:50 +#' +#' featurePlot(mat, pred, feature = "5") #' @export - featurePlot <- function(data, predictor, paired = NULL, covars = NULL, feature = NULL, relative = TRUE, logScale = FALSE, delta = 0.001, covar.quant = c(0,1/3,2/3,1)){ stopifnot(exists("data"),exists("predictor")) diff --git a/R/norm.R b/R/norm.R index 3118d29..bdbefa5 100644 --- a/R/norm.R +++ b/R/norm.R @@ -2,6 +2,7 @@ #' #' @param x numeric vector #' @return The geometric mean +#' @examples gm_mean(1:10) #' @export gm_mean = function(x){ if(any(x < 0, na.rm = TRUE)){ @@ -14,6 +15,7 @@ gm_mean = function(x){ #' #' @param x numeric matrix. Samples are columns #' @return A CLR normalized count_table +#' @examples norm_clr(matrix(1:100, nrow = 10, ncol = 10)) #' @export norm_clr <- function(x){ gm <- apply(x, 2, function(y) gm_mean(y)) @@ -25,6 +27,7 @@ norm_clr <- function(x){ #' @param x numeric matrix. Samples are columns #' @param ref reference feature #' @return An ALR normalized count_table +#' @examples norm_alr(matrix(1:100, nrow = 10, ncol = 10)) #' @export norm_alr <- function(x, ref = nrow(x)){ return(apply(x, 2, function(y) log(y/y[ref]))[-ref,]) diff --git a/R/posthocs.R b/R/posthocs.R index a1be8c7..899fcf8 100644 --- a/R/posthocs.R +++ b/R/posthocs.R @@ -6,6 +6,17 @@ #' @param p.adj P-value adjustment method. See \code{p.adjust} for details. Default "fdr" #' @param ... Additional arguments for \code{drop1} function #' @return A data.frame with output from drop1 and adjusted p.values for each predictor and feature +#' @examples +#' # Creating random count_table, predictor, and covariate +#' set.seed(5) +#' mat <- matrix(rnbinom(1500, size = 0.5, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' covar <- rnorm(15) +#' +#' # Running linear model and then drop1 on each feature +#' res <- DA.lmc(mat, pred, covars = list(Something = covar), allResults = TRUE) +#' res.drop <- DA.drop1(res) #' @export DA.drop1 <- function(results, test = "Chisq", p.adj = "fdr", ...){ @@ -132,6 +143,17 @@ DA.drop1 <- function(results, test = "Chisq", p.adj = "fdr", ...){ #' @param p.adj P-value adjustment method. See \code{p.adjust for details}. Default "fdr" #' @param ... Additional arguments for \code{anova} function #' @return A data.frame with output from anova and adjusted p.values for each predictor and feature +#' @examples +#' # Creating random count_table, predictor, and covariate +#' set.seed(5) +#' mat <- matrix(rnbinom(1500, size = 0.5, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' covar <- rnorm(15) +#' +#' # Running linear model and then anova on each feature +#' res <- DA.lmc(mat, pred, covars = list(Something = covar), allResults = TRUE) +#' res.ano <- DA.anova(res) #' @export DA.anova <- function(results, p.adj = "fdr", ...){ @@ -191,6 +213,16 @@ DA.anova <- function(results, p.adj = "fdr", ...){ #' @param p.adj P-value adjustment method. See \code{p.adjust} for details #' @param ... Additional arguments for \code{TukeyHSD} function #' @return A data.frame with output from TukeyHSD and adjusted p.values for each predictor and feature +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(3000, size = 0.1, mu = 500), nrow = 100, ncol = 30) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 10), rep("B", 10), rep("C", 10)) +#' +#' # Run anova and then TukeyHSD on each predictor +#' res <- DA.aov(data = mat, predictor = pred, allResults = TRUE) +#' res.tuk <- DA.TukeyHSD(res) #' @export DA.TukeyHSD <- function(results, variable = "predictor", p.adj = "fdr", ...){ @@ -237,6 +269,16 @@ DA.TukeyHSD <- function(results, variable = "predictor", p.adj = "fdr", ...){ #' @param p.adj P-value adjustment method. See \code{p.adjust} for details #' @param ... Additional arguments for \code{lsmeans} function #' @return A data.frame with output from lsmeans::pairs and adjusted p.values for each predictor and feature +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(1500, size = 0.5, mu = 500), nrow = 100, ncol = 15) +#' rownames(mat) <- 1:100 +#' pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +#' +#' # Running linear model and then pairwise comparisons on each feature +#' res <- DA.lmc(mat, pred, allResults = TRUE) +#' res.lsm <- DA.lsmeans(res) #' @export DA.lsmeans <- function(results, variable = "predictor", predictor = NULL, covars = NULL, p.adj = "fdr", ...){ diff --git a/R/powerDA.R b/R/powerDA.R index a2a26e6..b4cb40e 100644 --- a/R/powerDA.R +++ b/R/powerDA.R @@ -23,6 +23,33 @@ #' @import snow doSNOW foreach utils #' @importFrom parallel detectCores #' @importFrom pROC roc +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(1000, size = 0.5, mu = 500), nrow = 50, ncol = 20) +#' rownames(mat) <- 1:50 +#' pred <- c(rep("Control", 10), rep("Treatment", 10)) +#' +#' # Running powerDA on Wilcoxon test to test it with different effect sizes +#' # This example uses 1 core (cores = 1). +#' # Remove the cores argument to get it as high (and thereby fast) as possible. +#' res <- powerDA(data = mat, predictor = pred, test = "wil", cores = 1) +#' summary(res) +#' +#' \donttest{ +#' # Include a paired variable for dependent/blocked samples +#' subject <- rep(1:10, 2) +#' res <- powerDA(data = mat, predictor = pred, paired = subject, test = "ttt") +#' +#' # Include covariates +#' covar1 <- rnorm(20) +#' covar2 <- rep(c("A","B"), 10) +#' res <- powerDA(data = mat, predictor = pred, +#' covars = list(FirstCovar = covar1, CallItWhatYouWant = covar2), test = "lrm") +#' +#' # Data is absolute abundance +#' res <- powerDA(data = mat, predictor = pred, relative = FALSE, test = "ttt") +#' } #' @export powerDA <- function(data, predictor, paired = NULL, covars = NULL, test = NULL, effectSizes = c(2,4,8,16,32), alpha.p = 0.05, alpha.q = 0.1, p.adj = "fdr", R = 5, relative = TRUE, k = NULL, cores = (detectCores()-1), args = list(), out.all = NULL, core.check = TRUE, verbose = TRUE){ diff --git a/R/preDA.R b/R/preDA.R index 063f62f..77edc49 100644 --- a/R/preDA.R +++ b/R/preDA.R @@ -6,6 +6,16 @@ #' @param min.reads Minimum number of total reads the features should have. Default 0 #' @param min.abundance Minimum mean relative abundance features should have. Default 0 #' @return Similar to input, but with features not reaching the criteria given grouped as "Others" +#' @examples +#' # Creating random count_table with many low abundant +#' set.seed(4) +#' mat <- matrix(rnbinom(1000, size = 0.05, mu = 500), nrow = 100, ncol = 10) + +#' # Keep only those present in at least 3 samples +#' res <- preDA(mat, min.samples = 3) +#' +#' # Last feature is now called 'Others' and is a sum of all features present in less than 3 samples +#' rownames(res)[nrow(res)] #' @export preDA <- function(data, min.samples = 0, min.reads = 0, min.abundance = 0){ diff --git a/R/prune.tests.DA.R b/R/pruneTests.R similarity index 98% rename from R/prune.tests.DA.R rename to R/pruneTests.R index ffa0c3f..16ba65d 100644 --- a/R/prune.tests.DA.R +++ b/R/pruneTests.R @@ -10,7 +10,7 @@ #' @return A character with the tests argument pruned according to the parameters #' @export -prune.tests.DA <- function(tests, predictor, paired, covars, relative, decimal, zeroes){ +pruneTests <- function(tests, predictor, paired, covars, relative, decimal, zeroes){ # Prune test argument if packages are not installed if(!"baySeq" %in% rownames(installed.packages())) tests <- tests[tests != "bay"] diff --git a/R/runtimeDA.R b/R/runtimeDA.R index d8f7cc2..79421dd 100644 --- a/R/runtimeDA.R +++ b/R/runtimeDA.R @@ -17,6 +17,18 @@ #' @param ... Additional arguments for the \code{testDA} function #' @return A data.frame with estimated runtimes for 1 run #' @importFrom parallel detectCores +#' @examples +#' # Creating large random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(150000, size = 0.5, mu = 500), nrow = 10000, ncol = 10) +#' rownames(mat) <- 1:10000 +#' pred <- c(rep("A", 5), rep("B", 5)) +#' +#' # Use runtimeDA to predict total runtime for all features +#' # This example uses 1 core (cores = 1). +#' # Remove the cores argument to get it as high (and thereby fast) as possible. +#' # Also, in this example only a subset of tests are run. +#' runtimeDA(mat, pred, cores = 1, tests = c("ttt","wil"), tests.slow = c("neb")) #' @export runtimeDA <- function(data, predictor, paired = NULL, covars = NULL, subsamples = c(500,1000,1500,2000), subsamples.slow = c(100,150,200,250), tests = c("sam", "qua", "fri", "vli", "qpo", "pea", "wil", "ttt", "ltt", "ltt2","ere", "ere2", "msf", "zig", "lim", "lli", "lli2", "aov", "lao", "lao2", "kru", "lrm", "llm", "llm2", "spe", "aoa", "aoc", "tta", "ttc", "lma", "lmc", "lia", "lic"), diff --git a/R/summary.DAPower.R b/R/summary.DAPower.R index e05149c..429b53b 100644 --- a/R/summary.DAPower.R +++ b/R/summary.DAPower.R @@ -21,10 +21,10 @@ summary.DAPower <- function(object, decimals = 2, ...){ } else { if(x$Method[1] %in% c("ALDEx2 t-test (adx)","ALDEx2 wilcox (adx)")){ # Find medians - output.summary.sdr <- aggregate(Power ~ Method + EffectSize, data = x, FUN = function(y) median) - output.summary.auc <- aggregate(AUC ~ Method + EffectSize, data = x, FUN = function(y) median) - output.summary.fpr <- aggregate(FPR ~ Method + EffectSize, data = x, FUN = function(y) median) - output.summary.fdr <- aggregate(FDR ~ Method + EffectSize, data = x, FUN = function(y) median) + output.summary.sdr <- aggregate(Power ~ Method + EffectSize, data = x, FUN = median) + output.summary.auc <- aggregate(AUC ~ Method + EffectSize, data = x, FUN = median) + output.summary.fpr <- aggregate(FPR ~ Method + EffectSize, data = x, FUN = median) + output.summary.fdr <- aggregate(FDR ~ Method + EffectSize, data = x, FUN = median) # Merge df <- cbind(output.summary.sdr,output.summary.auc[,3],output.summary.fpr[,3],output.summary.fdr[,3]) @@ -32,10 +32,10 @@ summary.DAPower <- function(object, decimals = 2, ...){ colnames(df) <- c("Method","EffectSize","Power","AUC","FPR","FDR") } else { # Find medians - output.summary.sdr <- aggregate(Power ~ EffectSize, data = x, FUN = function(y) median) - output.summary.auc <- aggregate(AUC ~ EffectSize, data = x, FUN = function(y) median) - output.summary.fpr <- aggregate(FPR ~ EffectSize, data = x, FUN = function(y) median) - output.summary.fdr <- aggregate(FDR ~ EffectSize, data = x, FUN = function(y) median) + output.summary.sdr <- aggregate(Power ~ EffectSize, data = x, FUN = median) + output.summary.auc <- aggregate(AUC ~ EffectSize, data = x, FUN = median) + output.summary.fpr <- aggregate(FPR ~ EffectSize, data = x, FUN = median) + output.summary.fdr <- aggregate(FDR ~ EffectSize, data = x, FUN = median) # Merge df <- merge(merge(merge(output.summary.sdr,output.summary.auc, by = "EffectSize"), diff --git a/R/testDA.R b/R/testDA.R index 015b436..0618492 100644 --- a/R/testDA.R +++ b/R/testDA.R @@ -26,6 +26,35 @@ #' \item details - A dataframe with details from the run #' \item run.times - A dataframe with average run times of the different methods #' } +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(1000, size = 0.5, mu = 500), nrow = 50, ncol = 20) +#' rownames(mat) <- 1:50 +#' pred <- c(rep("Control", 10), rep("Treatment", 10)) +#' +#' # Running testDA to find the best method +#' # This example only repeats the test 1 time (R = 1). +#' # R should be increased to at least 20 for real test. +#' # It also uses 1 core (cores = 1). +#' # Remove this argument to get it as high (and thereby fast) as possible. +#' res <- testDA(data = mat, predictor = pred, cores = 1, R = 1) +#' summary(res) +#' +#' \donttest{ +#' # Include a paired variable for dependent/blocked samples +#' subject <- rep(1:10, 2) +#' res <- testDA(data = mat, predictor = pred, paired = subject) +#' +#' # Include covariates +#' covar1 <- rnorm(20) +#' covar2 <- rep(c("A","B"), 10) +#' res <- testDA(data = mat, predictor = pred, +#' covars = list(FirstCovar = covar1, CallItWhatYouWant = covar2)) +#' +#' # Data is absolute abundance +#' res <- testDA(data = mat, predictor = pred, relative = FALSE) +#' } #' #' @import stats snow doSNOW foreach utils doParallel #' @importFrom parallel detectCores @@ -92,7 +121,7 @@ testDA <- function(data, predictor, paired = NULL, covars = NULL, R = 20, if(!isTRUE(all(unlist(count_table) == floor(unlist(count_table))))) decimal <- TRUE if(any(count_table == 0)) zeroes <- TRUE tests <- unique(tests) - if(!"zzz" %in% tests) tests <- prune.tests.DA(tests, predictor, paired, covars, relative, decimal, zeroes) + if(!"zzz" %in% tests) tests <- pruneTests(tests, predictor, paired, covars, relative, decimal, zeroes) tests.par <- paste0(unlist(lapply(seq_len(R), function(x) rep(x,length(tests)))),"_",rep(tests,R)) if(length(tests) == 0) stop("No tests to run!") @@ -273,7 +302,7 @@ testDA <- function(data, predictor, paired = NULL, covars = NULL, R = 20, names(run.times) <- tests.par # Handle failed tests - results <- results[!vapply(results,is.null)] + results <- results[!sapply(results,is.null)] cat("\n") if(length(unique(gsub(".*_","",names(results)))) != length(tests)){ @@ -344,20 +373,20 @@ testDA <- function(data, predictor, paired = NULL, covars = NULL, R = 20, names(res.sub) <- newnames # Confusion matrix - totalPos <- vapply(res.sub,function(x) nrow(x[x$pval <= 0.05,])) - totalNeg <- vapply(res.sub,function(x) nrow(x[x$pval > 0.05,])) + totalPos <- sapply(res.sub,function(x) nrow(x[x$pval <= 0.05,])) + totalNeg <- sapply(res.sub,function(x) nrow(x[x$pval > 0.05,])) trueNeg <- totalNeg #if effectSize == 1 truePos <- 0 #if effectSize == 1 falseNeg <- 0 #if effectSize == 1 if(effectSize != 1){ - truePos <- vapply(res.sub, function(x) sum(x[x$pval <= 0.05,"Feature"] %in% spikeds[[r]][[2]])) - falseNeg <- vapply(res.sub, function(x) sum(x[x$pval > 0.05,"Feature"] %in% spikeds[[r]][[2]])) + truePos <- sapply(res.sub, function(x) sum(x[x$pval <= 0.05,"Feature"] %in% spikeds[[r]][[2]])) + falseNeg <- sapply(res.sub, function(x) sum(x[x$pval > 0.05,"Feature"] %in% spikeds[[r]][[2]])) } falsePos <- totalPos - truePos trueNeg <- totalNeg - falseNeg # FPR - fprs <- vapply(seq_along(res.sub), function(x) { + fprs <- sapply(seq_along(res.sub), function(x) { if((falsePos[x] + trueNeg[x]) != 0){ falsePos[x] / (falsePos[x] + trueNeg[x]) } else {0}}) @@ -367,12 +396,12 @@ testDA <- function(data, predictor, paired = NULL, covars = NULL, R = 20, # True positive for adjusted p-values truePos.adj <- 0 #if effectSize == 1 if(effectSize != 1){ - truePos.adj <- vapply(res.sub, function(x) sum(x[x$pval.adj <= alpha,"Feature"] %in% spikeds[[r]][[2]])) + truePos.adj <- sapply(res.sub, function(x) sum(x[x$pval.adj <= alpha,"Feature"] %in% spikeds[[r]][[2]])) } - sdrs <- vapply(seq_along(res.sub), function(x) truePos.adj[x] / sum(k)) + sdrs <- sapply(seq_along(res.sub), function(x) truePos.adj[x] / sum(k)) # AUC - aucs <- vapply(seq_along(res.sub), function(x) { + aucs <- sapply(seq_along(res.sub), function(x) { if(effectSize != 1){ test_roc <- NULL tryCatch( @@ -389,17 +418,17 @@ testDA <- function(data, predictor, paired = NULL, covars = NULL, R = 20, }) # Confusion matrix adjusted - totalPos.adj <- vapply(res.sub, function(x) nrow(x[x$pval.adj <= alpha,])) - truePos.adj <- vapply(res.sub, function(x) sum(x[x$pval.adj <= alpha,"Feature"] %in% spikeds[[r]][[2]])) + totalPos.adj <- sapply(res.sub, function(x) nrow(x[x$pval.adj <= alpha,])) + truePos.adj <- sapply(res.sub, function(x) sum(x[x$pval.adj <= alpha,"Feature"] %in% spikeds[[r]][[2]])) falsePos.adj <- totalPos.adj - truePos.adj - fdrs <- vapply(seq_along(res.sub), function(x) { + fdrs <- sapply(seq_along(res.sub), function(x) { if(totalPos.adj[x] != 0){ falsePos.adj[x] / totalPos.adj[x] } else {0}}) # Combine and return - df.combined <- data.frame(Method = vapply(res.sub, function(x) x$Method[1]), + df.combined <- data.frame(Method = sapply(res.sub, function(x) x$Method[1]), AUC = aucs, FPR = fprs, FDR = fdrs, diff --git a/R/vennDA.R b/R/vennDA.R index 9fbaa1f..34e492e 100644 --- a/R/vennDA.R +++ b/R/vennDA.R @@ -11,6 +11,21 @@ #' @param pkg Use either "eulerr" package (default) or "venneuler" for drawing diagrams. #' @param ... Additional arguments for plotting #' @return If output TRUE then a data.frame with Features detected by the different methods +#' @examples +#' # Creating random count_table and predictor +#' set.seed(5) +#' mat <- matrix(rnbinom(500, size = 0.1, mu = 500), nrow = 50, ncol = 10) +#' pred <- c(rep("Control", 5), rep("Treatment", 5)) +#' +#' # Running allDA to compare methods +#' # This example uses 1 core (cores = 1). +#' # Remove the cores argument to get it as high (and thereby fast) as possible. +#' res <- allDA(data = mat, predictor = pred, cores = 1) +#' +#' # Plot venn diagram comparing significant features from znb and zpo +#' # znb and zpo only have significant features due to high false positive rates in this example +#' # split = TRUE splits the significant features in positive and negative estimates +#' vennDA(res, tests = c("znb","zpo"), split = TRUE) #' @export vennDA <- function(x, tests = NULL, alpha = 0.1, split = FALSE, output = FALSE, pkg = "eulerr", ...){ diff --git a/man/DA.TukeyHSD.Rd b/man/DA.TukeyHSD.Rd index 81c9beb..e59aa66 100644 --- a/man/DA.TukeyHSD.Rd +++ b/man/DA.TukeyHSD.Rd @@ -21,3 +21,14 @@ A data.frame with output from TukeyHSD and adjusted p.values for each predictor \description{ Works on "aov", "lao", "lao2", "aoc", "aoa" } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(3000, size = 0.1, mu = 500), nrow = 100, ncol = 30) +rownames(mat) <- 1:100 +pred <- c(rep("A", 10), rep("B", 10), rep("C", 10)) + +# Run anova and then TukeyHSD on each predictor +res <- DA.aov(data = mat, predictor = pred, allResults = TRUE) +res.tuk <- DA.TukeyHSD(res) +} diff --git a/man/DA.adx.Rd b/man/DA.adx.Rd index 19e9431..e64ffe9 100644 --- a/man/DA.adx.Rd +++ b/man/DA.adx.Rd @@ -19,3 +19,13 @@ A data.frame with with results. \description{ Implementation of \code{aldex} for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running ALDEx2 +res <- DA.adx(data = mat, predictor = pred) +} diff --git a/man/DA.anova.Rd b/man/DA.anova.Rd index ee50d03..9182271 100644 --- a/man/DA.anova.Rd +++ b/man/DA.anova.Rd @@ -19,3 +19,15 @@ A data.frame with output from anova and adjusted p.values for each predictor and \description{ Works on "lrm", "llm", "llm2", "lma", "lmc". Non-paired "neb" } +\examples{ +# Creating random count_table, predictor, and covariate +set.seed(5) +mat <- matrix(rnbinom(1500, size = 0.5, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +covar <- rnorm(15) + +# Running linear model and then anova on each feature +res <- DA.lmc(mat, pred, covars = list(Something = covar), allResults = TRUE) +res.ano <- DA.anova(res) +} diff --git a/man/DA.aoa.Rd b/man/DA.aoa.Rd index 1e7942c..792684d 100644 --- a/man/DA.aoa.Rd +++ b/man/DA.aoa.Rd @@ -31,3 +31,13 @@ Apply ANOVA on multiple features with one \code{predictor}. \details{ Note: Last feature in the data is used as reference for the log-ratio transformation. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running ANOVA on each feature +res <- DA.aoa(data = mat, predictor = pred) +} diff --git a/man/DA.aoc.Rd b/man/DA.aoc.Rd index f9a16fa..459a6fa 100644 --- a/man/DA.aoc.Rd +++ b/man/DA.aoc.Rd @@ -28,3 +28,13 @@ A data.frame with with results. \description{ Apply ANOVA on multiple features with one \code{predictor}. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running ANOVA on each feature +res <- DA.aoc(data = mat, predictor = pred) +} diff --git a/man/DA.aov.Rd b/man/DA.aov.Rd index 28fa935..dc39203 100644 --- a/man/DA.aov.Rd +++ b/man/DA.aov.Rd @@ -28,3 +28,13 @@ A data.frame with with results. \description{ Run ANOVA on multiple features with \code{predictor} as independent variable } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running ANOVA on each feature +res <- DA.aov(data = mat, predictor = pred) +} diff --git a/man/DA.bay.Rd b/man/DA.bay.Rd index b43bea1..ef118cc 100644 --- a/man/DA.bay.Rd +++ b/man/DA.bay.Rd @@ -21,3 +21,13 @@ A data.frame with with results. \description{ Implementation of baySeq for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(100, size = 0.1, mu = 500), nrow = 50, ncol = 10) +rownames(mat) <- 1:50 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running baySeq +res <- DA.bay(data = mat, predictor = pred) +} diff --git a/man/DA.drop1.Rd b/man/DA.drop1.Rd index cbffcef..fd3b29b 100644 --- a/man/DA.drop1.Rd +++ b/man/DA.drop1.Rd @@ -21,3 +21,15 @@ A data.frame with output from drop1 and adjusted p.values for each predictor and \description{ Works on "zpo", "znb", "qpo", "neb", "poi". Non-paired "lrm", "llm", "llm2", "lma", "lmc" } +\examples{ +# Creating random count_table, predictor, and covariate +set.seed(5) +mat <- matrix(rnbinom(1500, size = 0.5, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +covar <- rnorm(15) + +# Running linear model and then drop1 on each feature +res <- DA.lmc(mat, pred, covars = list(Something = covar), allResults = TRUE) +res.drop <- DA.drop1(res) +} diff --git a/man/DA.ds2.Rd b/man/DA.ds2.Rd index 6476a86..b207da0 100644 --- a/man/DA.ds2.Rd +++ b/man/DA.ds2.Rd @@ -35,3 +35,13 @@ A data.frame with with results. Implementation of DESeq2 for \code{DAtest} Manual geometric means calculated to avoid errors, see https://github.com/joey711/phyloseq/issues/387 } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(200, size = 0.1, mu = 500), nrow = 20, ncol = 10) +rownames(mat) <- 1:20 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running DESeq2 +res <- DA.ds2(data = mat, predictor = pred) +} diff --git a/man/DA.ds2x.Rd b/man/DA.ds2x.Rd index b853846..38e0d9a 100644 --- a/man/DA.ds2x.Rd +++ b/man/DA.ds2x.Rd @@ -35,3 +35,13 @@ A data.frame with with results. \description{ Implementation of DESeq2 for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(6) +mat <- matrix(rnbinom(200, size = 0.2, mu = 500), nrow = 20, ncol = 10) +rownames(mat) <- 1:20 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running DESeq2 +res <- DA.ds2x(data = mat, predictor = pred) +} diff --git a/man/DA.ere.Rd b/man/DA.ere.Rd index f88312e..356d2ef 100644 --- a/man/DA.ere.Rd +++ b/man/DA.ere.Rd @@ -21,3 +21,13 @@ A data.frame with with results. \description{ Implementation of edgeR exact test for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running edgeR +res <- DA.ere(data = mat, predictor = pred) +} diff --git a/man/DA.ere2.Rd b/man/DA.ere2.Rd index b82a644..e417230 100644 --- a/man/DA.ere2.Rd +++ b/man/DA.ere2.Rd @@ -21,3 +21,13 @@ A data.frame with with results. \description{ Implementation of edgeR exact test for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running edgeR +res <- DA.ere2(data = mat, predictor = pred) +} diff --git a/man/DA.erq.Rd b/man/DA.erq.Rd index afb8d70..b4cc46e 100644 --- a/man/DA.erq.Rd +++ b/man/DA.erq.Rd @@ -32,3 +32,13 @@ A data.frame with with results. \description{ Implementation of edgeR quasi-likelihood test for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running edgeR +res <- DA.erq(data = mat, predictor = pred) +} diff --git a/man/DA.erq2.Rd b/man/DA.erq2.Rd index deb0f95..d04a2b3 100644 --- a/man/DA.erq2.Rd +++ b/man/DA.erq2.Rd @@ -33,3 +33,13 @@ A data.frame with with results. \description{ Implementation of edgeR quasi-likelihood test for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running edgeR +res <- DA.erq2(data = mat, predictor = pred) +} diff --git a/man/DA.fri.Rd b/man/DA.fri.Rd index 3a08301..bc2569f 100644 --- a/man/DA.fri.Rd +++ b/man/DA.fri.Rd @@ -28,3 +28,14 @@ A data.frame with with results. \description{ Apply friedman test to multiple features with one \code{predictor} } +\examples{ +# Creating random count_table, predictor, and paired variable +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +subject <- rep(1:5, 3) + +# Running Friedman test on each feature +res <- DA.fri(data = mat, predictor = pred, paired = subject) +} diff --git a/man/DA.kru.Rd b/man/DA.kru.Rd index 61055ad..6a11b77 100644 --- a/man/DA.kru.Rd +++ b/man/DA.kru.Rd @@ -26,3 +26,13 @@ A data.frame with with results. \description{ Apply kruskal-wallis test on multiple features with one \code{predictor} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running Kruskal-Wallis on each feature +res <- DA.kru(data = mat, predictor = pred) +} diff --git a/man/DA.lao.Rd b/man/DA.lao.Rd index 40b20d1..fd40033 100644 --- a/man/DA.lao.Rd +++ b/man/DA.lao.Rd @@ -30,3 +30,13 @@ A data.frame with with results. \description{ Apply ANOVA on multiple features with one \code{predictor}, with log transformation of counts before normalization. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running ANOVA on each feature +res <- DA.lao(data = mat, predictor = pred) +} diff --git a/man/DA.lao2.Rd b/man/DA.lao2.Rd index c5cb75f..a4018c5 100644 --- a/man/DA.lao2.Rd +++ b/man/DA.lao2.Rd @@ -28,3 +28,13 @@ A data.frame with with results. \description{ Apply ANOVA to multiple features with one \code{predictor}, with log transformation of relative abundances. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running ANOVA on each feature +res <- DA.lao2(data = mat, predictor = pred) +} diff --git a/man/DA.lia.Rd b/man/DA.lia.Rd index 7db5444..a207ad1 100644 --- a/man/DA.lia.Rd +++ b/man/DA.lia.Rd @@ -34,3 +34,13 @@ A data.frame with with results. \description{ Note: Last feature in the data is used as reference for the log-ratio transformation. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running limma +res <- DA.lia(data = mat, predictor = pred) +} diff --git a/man/DA.lic.Rd b/man/DA.lic.Rd index 369d85a..8c09dc3 100644 --- a/man/DA.lic.Rd +++ b/man/DA.lic.Rd @@ -34,3 +34,13 @@ A data.frame with with results. \description{ Note: Last feature in the data is used as reference for the log-ratio transformation. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running limma +res <- DA.lic(data = mat, predictor = pred) +} diff --git a/man/DA.lim.Rd b/man/DA.lim.Rd index ee05d7f..96c29e9 100644 --- a/man/DA.lim.Rd +++ b/man/DA.lim.Rd @@ -35,3 +35,13 @@ A data.frame with with results. \description{ LIMMA } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running limma +res <- DA.lim(data = mat, predictor = pred) +} diff --git a/man/DA.lli.Rd b/man/DA.lli.Rd index fa91f37..cf74828 100644 --- a/man/DA.lli.Rd +++ b/man/DA.lli.Rd @@ -37,3 +37,13 @@ A data.frame with with results. \description{ Limma with log transformation for counts before normalization. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running limma +res <- DA.lli(data = mat, predictor = pred) +} diff --git a/man/DA.lli2.Rd b/man/DA.lli2.Rd index 98423d8..09ea5c0 100644 --- a/man/DA.lli2.Rd +++ b/man/DA.lli2.Rd @@ -35,3 +35,13 @@ A data.frame with with results. \description{ Limma with log transformation of relative abundances. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running limma +res <- DA.lli2(data = mat, predictor = pred) +} diff --git a/man/DA.llm.Rd b/man/DA.llm.Rd index cd0c612..598b8d1 100644 --- a/man/DA.llm.Rd +++ b/man/DA.llm.Rd @@ -38,3 +38,13 @@ A data.frame with with results. Apply linear regression for multiple features with one \code{predictor}, with log transformation of counts before normalization. Mixed-effect model is used when a \code{paired} argument is included, with the \code{paired} variable as a random intercept. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running linear model on each feature +res <- DA.llm(data = mat, predictor = pred) +} diff --git a/man/DA.llm2.Rd b/man/DA.llm2.Rd index cfe77dc..71eab2f 100644 --- a/man/DA.llm2.Rd +++ b/man/DA.llm2.Rd @@ -36,3 +36,13 @@ A data.frame with with results. Apply linear regression on multiple features with one \code{predictor}, with log transformation of relative abundances. Mixed-effect model is used when a \code{paired} argument is included, with the \code{paired} variable as a random intercept. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running linear model on each feature +res <- DA.llm2(data = mat, predictor = pred) +} diff --git a/man/DA.lma.Rd b/man/DA.lma.Rd index e413fc5..ddf14c3 100644 --- a/man/DA.lma.Rd +++ b/man/DA.lma.Rd @@ -36,3 +36,13 @@ Apply linear regression for multiple features with one \code{predictor}. Mixed-effect model is used when a \code{paired} argument is included, with the \code{paired} variable as a random intercept. Note: Last feature in the data is used as reference for the log-ratio transformation. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running linear model on each feature +res <- DA.lma(data = mat, predictor = pred) +} diff --git a/man/DA.lmc.Rd b/man/DA.lmc.Rd index d511724..ac5c1b8 100644 --- a/man/DA.lmc.Rd +++ b/man/DA.lmc.Rd @@ -35,3 +35,13 @@ A data.frame with with results. Apply linear regression for multiple features with one \code{predictor}. Mixed-effect model is used when a \code{paired} argument is included, with the \code{paired} variable as a random intercept. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running linear model on each feature +res <- DA.lmc(data = mat, predictor = pred) +} diff --git a/man/DA.lrm.Rd b/man/DA.lrm.Rd index 333da69..c8550ee 100644 --- a/man/DA.lrm.Rd +++ b/man/DA.lrm.Rd @@ -36,3 +36,13 @@ A data.frame with with results. Apply linear regression to mulitple features with one \code{predictor} Mixed-effect model is used when a \code{paired} argument is present, with the \code{paired} variable as random intercept } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running linear model on each feature +res <- DA.lrm(data = mat, predictor = pred) +} diff --git a/man/DA.lsmeans.Rd b/man/DA.lsmeans.Rd index 8d17098..2d870aa 100644 --- a/man/DA.lsmeans.Rd +++ b/man/DA.lsmeans.Rd @@ -29,3 +29,14 @@ Pairwise testing on predictor and covars. Works on "poi", "neb", "lrm", "lma", " \details{ Require the \code{lsmeans} package } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(1500, size = 0.5, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) + +# Running linear model and then pairwise comparisons on each feature +res <- DA.lmc(mat, pred, allResults = TRUE) +res.lsm <- DA.lsmeans(res) +} diff --git a/man/DA.ltt.Rd b/man/DA.ltt.Rd index bac39c5..8735b40 100644 --- a/man/DA.ltt.Rd +++ b/man/DA.ltt.Rd @@ -37,3 +37,13 @@ A data.frame with with results. \description{ Apply welch t-test to multiple features with one \code{predictor}, with log transformaiton of counts before normalization } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running t-test on each feature +res <- DA.ltt(data = mat, predictor = pred) +} diff --git a/man/DA.ltt2.Rd b/man/DA.ltt2.Rd index 4603edf..65ddae1 100644 --- a/man/DA.ltt2.Rd +++ b/man/DA.ltt2.Rd @@ -35,3 +35,13 @@ A data.frame with with results. \description{ Apply welch t-test to multiple features and one \code{predictor}, and with log transformed relative abundances } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running t-test on each feature +res <- DA.ltt2(data = mat, predictor = pred) +} diff --git a/man/DA.msf.Rd b/man/DA.msf.Rd index 32a979f..b08e7ee 100644 --- a/man/DA.msf.Rd +++ b/man/DA.msf.Rd @@ -24,3 +24,13 @@ A data.frame with with results. Implemented as in: https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-016-0208-8 } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running MetagenomeSeq feature model +res <- DA.msf(data = mat, predictor = pred) +} diff --git a/man/DA.mva.Rd b/man/DA.mva.Rd index 8a301ad..7204a7d 100644 --- a/man/DA.mva.Rd +++ b/man/DA.mva.Rd @@ -37,3 +37,13 @@ A data.frame with with results. \description{ Implementation of mvabund manyglm for \code{DAtest}. With negative binomial family and an offset of log(LibrarySize) when relative = TRUE } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(100, size = 0.1, mu = 500), nrow = 10, ncol = 10) +rownames(mat) <- 1:10 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running mvabund manyglm +res <- DA.mva(data = mat, predictor = pred) +} diff --git a/man/DA.neb.Rd b/man/DA.neb.Rd index 7eb3e5a..f742c1d 100644 --- a/man/DA.neb.Rd +++ b/man/DA.neb.Rd @@ -39,3 +39,13 @@ Apply negative binomial generalized linear model for multiple features with one With \code{log(librarySize)} as offset if \code{relative=TRUE}. Mixed-effect model is used when a \code{paired} argument is included, with the \code{paired} variable as a random intercept. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running Negative Binomial regression on each feature +res <- DA.neb(data = mat, predictor = pred) +} diff --git a/man/DA.pea.Rd b/man/DA.pea.Rd index 63214e7..0962b3d 100644 --- a/man/DA.pea.Rd +++ b/man/DA.pea.Rd @@ -23,3 +23,13 @@ A data.frame with with results. \description{ Apply pearson correlation between multiple features and one \code{predictor} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 50, ncol = 20) +rownames(mat) <- 1:50 +pred <- rnorm(20) + +# Running Pearson correlation on each feature +res <- DA.pea(data = mat, predictor = pred) +} diff --git a/man/DA.per.Rd b/man/DA.per.Rd index 6089838..b73b1ae 100644 --- a/man/DA.per.Rd +++ b/man/DA.per.Rd @@ -42,3 +42,13 @@ P-values are now two-sided, and test statistic is a simple log fold change A paired permutation test is implemented specifically for this package. The test is similar to the original, but with a different test statistic and permutation scheme. The permutations are constrained in the paired version such that the predictor is only permuted within each level of the paired argument (e.g. subjects). The test statistic first finds the log-ratio between the two predictor levels (e.g. case and control) for each level of the paired argument and the final statistic is the mean of these log-ratios. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running permutation test on each feature +res <- DA.per(data = mat, predictor = pred) +} diff --git a/man/DA.phyloseq.Rd b/man/DA.phyloseq.Rd index c1f2cae..579398a 100644 --- a/man/DA.phyloseq.Rd +++ b/man/DA.phyloseq.Rd @@ -19,5 +19,5 @@ DA.phyloseq(data, predictor, paired = NULL, covars = NULL) A list with data extracted from the \code{phyloseq} object } \description{ -Extract data from a \code{phyloseq} object to be used in \code{DAtest} +Internal function } diff --git a/man/DA.poi.Rd b/man/DA.poi.Rd index a964833..7580b9b 100644 --- a/man/DA.poi.Rd +++ b/man/DA.poi.Rd @@ -39,3 +39,13 @@ Apply poisson generalized linear models on multiple features with one \code{pred With \code{log(librarySize)} as offset. Mixed-effect model is used when a \code{paired} argument is included, with the \code{paired} variable as a random intercept. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running Poisson regression on each feature +res <- DA.poi(data = mat, predictor = pred) +} diff --git a/man/DA.qpo.Rd b/man/DA.qpo.Rd index b897105..32cbab3 100644 --- a/man/DA.qpo.Rd +++ b/man/DA.qpo.Rd @@ -36,3 +36,13 @@ A data.frame with with results. Apply quasi-poisson generalized linear model for multiple features with one \code{predictor} With \code{log(librarySize)} as offset if \code{relative=TRUE}. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running Quasi-Poisson regression on each feature +res <- DA.qpo(data = mat, predictor = pred) +} diff --git a/man/DA.qua.Rd b/man/DA.qua.Rd index 3d35c14..a78c131 100644 --- a/man/DA.qua.Rd +++ b/man/DA.qua.Rd @@ -28,3 +28,14 @@ A data.frame with with results. \description{ Apply \code{quade.test} to multiple features with one \code{predictor} } +\examples{ +# Creating random count_table, predictor, and paired variable +set.seed(4) +mat <- matrix(rnbinom(1500, size = 0.1, mu = 500), nrow = 100, ncol = 15) +rownames(mat) <- 1:100 +pred <- c(rep("A", 5), rep("B", 5), rep("C", 5)) +subject <- rep(1:5, 3) + +# Running Quade test on each feature +res <- DA.qua(data = mat, predictor = pred, paired = subject) +} diff --git a/man/DA.sam.Rd b/man/DA.sam.Rd index 5dcf44f..b28e0bf 100644 --- a/man/DA.sam.Rd +++ b/man/DA.sam.Rd @@ -26,3 +26,13 @@ A data.frame with with results. \description{ SAMSeq implementation for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running SamSeq +res <- DA.sam(data = mat, predictor = pred) +} diff --git a/man/DA.spe.Rd b/man/DA.spe.Rd index c157929..a6986e9 100644 --- a/man/DA.spe.Rd +++ b/man/DA.spe.Rd @@ -23,3 +23,13 @@ A data.frame with with results. \description{ Apply Spearman correlation to multiple features with one \code{predictor} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 50, ncol = 20) +rownames(mat) <- 1:50 +pred <- rnorm(20) + +# Running Spearman correlation on each feature +res <- DA.spe(data = mat, predictor = pred) +} diff --git a/man/DA.tta.Rd b/man/DA.tta.Rd index 391a547..be208a1 100644 --- a/man/DA.tta.Rd +++ b/man/DA.tta.Rd @@ -37,3 +37,13 @@ Apply welch t-test to multiple features with one \code{predictor} \details{ Note: Last feature in the data is used as reference for the log-ratio transformation. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running t-test on each feature +res <- DA.tta(data = mat, predictor = pred) +} diff --git a/man/DA.ttc.Rd b/man/DA.ttc.Rd index 7fa9a8f..b7fe048 100644 --- a/man/DA.ttc.Rd +++ b/man/DA.ttc.Rd @@ -34,3 +34,13 @@ A data.frame with with results. \description{ Apply welch t-test to multiple features with one \code{predictor} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running t-test on each feature +res <- DA.ttc(data = mat, predictor = pred) +} diff --git a/man/DA.ttt.Rd b/man/DA.ttt.Rd index 9cff4e7..c26fd44 100644 --- a/man/DA.ttt.Rd +++ b/man/DA.ttt.Rd @@ -35,3 +35,13 @@ A data.frame with with results. \description{ Apply welch t-test to multiple features with one \code{predictor} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running t-test +res <- DA.ttt(data = mat, predictor = pred) +} diff --git a/man/DA.vli.Rd b/man/DA.vli.Rd index 0fd1848..74234f3 100644 --- a/man/DA.vli.Rd +++ b/man/DA.vli.Rd @@ -32,3 +32,13 @@ A data.frame with with results. \description{ Limma is run after variance stabilization with \code{voom} using \code{edgeR} normalized library sizes } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running limma +res <- DA.vli(data = mat, predictor = pred) +} diff --git a/man/DA.wil.Rd b/man/DA.wil.Rd index 036e660..9ea2e55 100644 --- a/man/DA.wil.Rd +++ b/man/DA.wil.Rd @@ -35,3 +35,13 @@ A data.frame with with results. \description{ Apply wilcoxon test for multiple features with one \code{predictor} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running Wilcoxon test on each feature +res <- DA.wil(data = mat, predictor = pred) +} diff --git a/man/DA.zig.Rd b/man/DA.zig.Rd index a16d6ab..5ca8623 100644 --- a/man/DA.zig.Rd +++ b/man/DA.zig.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/DA.zig.R \name{DA.zig} \alias{DA.zig} -\title{MetageonomeSeq ZIG} +\title{MetagenomeSeq ZIG} \usage{ DA.zig(data, predictor, paired = NULL, covars = NULL, p.adj = "fdr", by = 2, eff = 0.5, allResults = FALSE, ...) @@ -32,3 +32,13 @@ A data.frame with with results. \description{ Implementation of Metagenome zero-inflated gaussian model for \code{DAtest} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running MetagenomeSeq Zero-inflated Gaussian +res <- DA.zig(data = mat, predictor = pred) +} diff --git a/man/DA.znb.Rd b/man/DA.znb.Rd index ee35a89..90c126f 100644 --- a/man/DA.znb.Rd +++ b/man/DA.znb.Rd @@ -36,3 +36,13 @@ A data.frame with with results. Apply zero-inflated negative binomial generalized linear model to multiple features, with one independent variable With \code{log(librarySize)} as offset. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running Zero-inflated Negative Binomial regression on each feature +res <- DA.znb(data = mat, predictor = pred) +} diff --git a/man/DA.zpo.Rd b/man/DA.zpo.Rd index f50c5f6..8c3e3e6 100644 --- a/man/DA.zpo.Rd +++ b/man/DA.zpo.Rd @@ -36,3 +36,13 @@ A data.frame with with results. Apply zero-inflated poisson generalized linear model for multiple features, with one \code{predictor} as independent variable. With \code{log(librarySize)} as offset when \code{relative=TRUE}. } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running Zero-inflated Poisson regression on each feature +res <- DA.zpo(data = mat, predictor = pred) +} diff --git a/man/DA.zzz.Rd b/man/DA.zzz.Rd index 070f1e3..045e53a 100644 --- a/man/DA.zzz.Rd +++ b/man/DA.zzz.Rd @@ -26,3 +26,37 @@ A data.frame with results from the user-defined method \description{ Apply a user-defined function to all features of a count table. For implemetation in \code{testDA} and \code{allDA} } +\examples{ +# Creating random count_table and predictor +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.1, mu = 500), nrow = 100, ncol = 10) +rownames(mat) <- 1:100 +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Define function for t-test +myfun <- function(count_table, predictor, paired, covars){ + +# Relative abundance +rel <- apply(count_table, 2, function(x) x/sum(x)) + +# t-test function +# Wrapping this function in tryCatch(..., error = function(e){NA}) +# ensures that our main function won't fail if t.test fails on some features +tfun <- function(x){ + tryCatch(t.test(x ~ predictor)$p.value, error = function(e){NA}) +} + +# P-values for each feature +pvals <- apply(rel, 1, tfun) + +# Collect and return data +df <- data.frame(Feature = rownames(count_table), + pval = pvals) +df$pval.adj <- p.adjust(df$pval, method = "fdr") +df$Method <- "My own t-test" +return(df) +} + +# Running the test +res <- DA.zzz(data = mat, predictor = pred, FUN = myfun) +} diff --git a/man/add.tax.DA.Rd b/man/add.tax.DA.Rd index 7641966..a2ff32b 100644 --- a/man/add.tax.DA.Rd +++ b/man/add.tax.DA.Rd @@ -12,7 +12,7 @@ add.tax.DA(data, res) \item{res}{data.frame with results} } \value{ -A \code{res} merged with the tax.table from the phyloseq object +\code{res} merged with the tax.table from the phyloseq object } \description{ Internal function diff --git a/man/allDA.Rd b/man/allDA.Rd index 1a69a88..3823b9a 100644 --- a/man/allDA.Rd +++ b/man/allDA.Rd @@ -56,3 +56,35 @@ Run many differential abundance and expression tests at a time, to easily compar \details{ mva is excluded by default, as it is slow. } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(500, size = 0.1, mu = 500), nrow = 50, ncol = 10) +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running allDA to compare methods +# This example uses 1 core (cores = 1). +# Remove the cores argument to get it as high (and thereby fast) as possible. +res <- allDA(data = mat, predictor = pred, cores = 1) + +# View adjusted p-values from all methods +print(res$adj) + +# View estimates from all methods +print(res$est) + +\donttest{ +# Include a paired variable for dependent/blocked samples +subject <- rep(1:5, 2) +res <- allDA(data = mat, predictor = pred, paired = subject) + +# Include covariates +covar1 <- rnorm(10) +covar2 <- rep(c("A","B"), 5) +res <- allDA(data = mat, predictor = pred, + covars = list(FirstCovar = covar1, CallItWhatYouWant = covar2)) + +# Data is absolute abundance +res <- allDA(data = mat, predictor = pred, relative = FALSE) +} +} diff --git a/man/featurePlot.Rd b/man/featurePlot.Rd index e5de922..d4bd41f 100644 --- a/man/featurePlot.Rd +++ b/man/featurePlot.Rd @@ -38,3 +38,12 @@ Boxplot for categorical variables, points and smooth line for quantitative varia If a \code{paired} variable is supplied, it is always plotted as points with lines grouped by the \code{paired} variable If \code{covars} are supplied data is split in facets. Quantitative covars are cut in intervals according to the quantiles given in \code{covar.quant} } +\examples{ +# Create random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(500, size = 0.1, mu = 500), nrow = 50, ncol = 10) +pred <- c(rep("Control", 5), rep("Treatment", 5)) +rownames(mat) <- 1:50 + +featurePlot(mat, pred, feature = "5") +} diff --git a/man/gm_mean.Rd b/man/gm_mean.Rd index 7887ab1..c9c1c06 100644 --- a/man/gm_mean.Rd +++ b/man/gm_mean.Rd @@ -15,3 +15,6 @@ The geometric mean \description{ Geometric means } +\examples{ +gm_mean(1:10) +} diff --git a/man/norm_alr.Rd b/man/norm_alr.Rd index f91fdde..da7da2b 100644 --- a/man/norm_alr.Rd +++ b/man/norm_alr.Rd @@ -17,3 +17,6 @@ An ALR normalized count_table \description{ Additive log-ratio normalization } +\examples{ +norm_alr(matrix(1:100, nrow = 10, ncol = 10)) +} diff --git a/man/norm_clr.Rd b/man/norm_clr.Rd index 71b3aa1..214593e 100644 --- a/man/norm_clr.Rd +++ b/man/norm_clr.Rd @@ -15,3 +15,6 @@ A CLR normalized count_table \description{ Centered log-ratio normalization } +\examples{ +norm_clr(matrix(1:100, nrow = 10, ncol = 10)) +} diff --git a/man/powerDA.Rd b/man/powerDA.Rd index f762eb4..d15d103 100644 --- a/man/powerDA.Rd +++ b/man/powerDA.Rd @@ -54,3 +54,31 @@ Estimating (empirical) statistical power for a specific differential abundance a \details{ Currently implemented methods: see \code{testDA} } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(1000, size = 0.5, mu = 500), nrow = 50, ncol = 20) +rownames(mat) <- 1:50 +pred <- c(rep("Control", 10), rep("Treatment", 10)) + +# Running powerDA on Wilcoxon test to test it with different effect sizes +# This example uses 1 core (cores = 1). +# Remove the cores argument to get it as high (and thereby fast) as possible. +res <- powerDA(data = mat, predictor = pred, test = "wil", cores = 1) +summary(res) + +\donttest{ +# Include a paired variable for dependent/blocked samples +subject <- rep(1:10, 2) +res <- powerDA(data = mat, predictor = pred, paired = subject, test = "ttt") + +# Include covariates +covar1 <- rnorm(20) +covar2 <- rep(c("A","B"), 10) +res <- powerDA(data = mat, predictor = pred, + covars = list(FirstCovar = covar1, CallItWhatYouWant = covar2), test = "lrm") + +# Data is absolute abundance +res <- powerDA(data = mat, predictor = pred, relative = FALSE, test = "ttt") +} +} diff --git a/man/preDA.Rd b/man/preDA.Rd index c6f190d..c1c75bd 100644 --- a/man/preDA.Rd +++ b/man/preDA.Rd @@ -21,3 +21,13 @@ Similar to input, but with features not reaching the criteria given grouped as " \description{ Pre-process the count table before running \code{DAtest} functions } +\examples{ +# Creating random count_table with many low abundant +set.seed(4) +mat <- matrix(rnbinom(1000, size = 0.05, mu = 500), nrow = 100, ncol = 10) +# Keep only those present in at least 3 samples +res <- preDA(mat, min.samples = 3) + +# Last feature is now called 'Others' and is a sum of all features present in less than 3 samples +rownames(res)[nrow(res)] +} diff --git a/man/prune.tests.DA.Rd b/man/pruneTests.Rd similarity index 80% rename from man/prune.tests.DA.Rd rename to man/pruneTests.Rd index 735dba6..5cca0bd 100644 --- a/man/prune.tests.DA.Rd +++ b/man/pruneTests.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prune.tests.DA.R -\name{prune.tests.DA} -\alias{prune.tests.DA} +% Please edit documentation in R/pruneTests.R +\name{pruneTests} +\alias{pruneTests} \title{Prune tests argument for \code{testDA}} \usage{ -prune.tests.DA(tests, predictor, paired, covars, relative, decimal, zeroes) +pruneTests(tests, predictor, paired, covars, relative, decimal, zeroes) } \arguments{ \item{tests}{A character vector with names of tests} diff --git a/man/runtimeDA.Rd b/man/runtimeDA.Rd index 84ac6c2..686a772 100644 --- a/man/runtimeDA.Rd +++ b/man/runtimeDA.Rd @@ -45,3 +45,16 @@ Outputs the estimated times for running each method 1 time. With cores=1 the run Runtime of all methods are expected to scale linearly with the number of features, except "anc" and "bay" which are modelled with a 2. order polynomial. } +\examples{ +# Creating large random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(150000, size = 0.5, mu = 500), nrow = 10000, ncol = 10) +rownames(mat) <- 1:10000 +pred <- c(rep("A", 5), rep("B", 5)) + +# Use runtimeDA to predict total runtime for all features +# This example uses 1 core (cores = 1). +# Remove the cores argument to get it as high (and thereby fast) as possible. +# Also, in this example only a subset of tests are run. +runtimeDA(mat, pred, cores = 1, tests = c("ttt","wil"), tests.slow = c("neb")) +} diff --git a/man/testDA.Rd b/man/testDA.Rd index 2db0bbd..16714ef 100644 --- a/man/testDA.Rd +++ b/man/testDA.Rd @@ -62,3 +62,34 @@ Calculating Power, False Discovery Rates, False Positive Rates and AUC (Area Und \details{ mva is excluded by default, as it is slow. } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(1000, size = 0.5, mu = 500), nrow = 50, ncol = 20) +rownames(mat) <- 1:50 +pred <- c(rep("Control", 10), rep("Treatment", 10)) + +# Running testDA to find the best method +# This example only repeats the test 1 time (R = 1). +# R should be increased to at least 20 for real test. +# It also uses 1 core (cores = 1). +# Remove this argument to get it as high (and thereby fast) as possible. +res <- testDA(data = mat, predictor = pred, cores = 1, R = 1) +summary(res) + +\donttest{ +# Include a paired variable for dependent/blocked samples +subject <- rep(1:10, 2) +res <- testDA(data = mat, predictor = pred, paired = subject) + +# Include covariates +covar1 <- rnorm(20) +covar2 <- rep(c("A","B"), 10) +res <- testDA(data = mat, predictor = pred, + covars = list(FirstCovar = covar1, CallItWhatYouWant = covar2)) + +# Data is absolute abundance +res <- testDA(data = mat, predictor = pred, relative = FALSE) +} + +} diff --git a/man/vennDA.Rd b/man/vennDA.Rd index 37817fa..a0ae003 100644 --- a/man/vennDA.Rd +++ b/man/vennDA.Rd @@ -31,3 +31,19 @@ Plot a Venn (Euler) diagram of features found by different methods. \details{ Require the eulerr package unless output is TRUE. } +\examples{ +# Creating random count_table and predictor +set.seed(5) +mat <- matrix(rnbinom(500, size = 0.1, mu = 500), nrow = 50, ncol = 10) +pred <- c(rep("Control", 5), rep("Treatment", 5)) + +# Running allDA to compare methods +# This example uses 1 core (cores = 1). +# Remove the cores argument to get it as high (and thereby fast) as possible. +res <- allDA(data = mat, predictor = pred, cores = 1) + +# Plot venn diagram comparing significant features from znb and zpo +# znb and zpo only have significant features due to high false positive rates in this example +# split = TRUE splits the significant features in positive and negative estimates +vennDA(res, tests = c("znb","zpo"), split = TRUE) +}