New version

DESCRIPTION

@@ -1,10 +1,11 @@
 Package: DAtest
 Title: Comparing Differential Abundance methods
-Version: 2.6.1
+Version: 2.6.2
 Authors@R: person("Jakob", "Russel", email = "[email protected]", role = c("aut", "cre"))
 Description: What the title says.
 Depends: R (>= 3.2.5)
 Imports: foreach, doParallel, doSNOW, snow, MASS, pROC, ggplot2, cowplot, lme4, nlme, stats, methods, utils, statmod, pscl, samr, minqa, nloptr 
+Suggests: lsmeans, venneuler
 License: GPL (>= 3)
 Encoding: UTF-8
 LazyData: true

NAMESPACE

@@ -3,10 +3,13 @@
 S3method(plot,DA)
 S3method(print,DA)
 S3method(summary,DA)
+export(DA.TukeyHSD)
 export(DA.adx)
 export(DA.anc)
+export(DA.anova)
 export(DA.aov)
 export(DA.bay)
+export(DA.drop1)
 export(DA.ds2)
 export(DA.ere)
 export(DA.ere2)
@@ -22,6 +25,7 @@ export(DA.lli2)
 export(DA.llm)
 export(DA.llm2)
 export(DA.lrm)
+export(DA.lsmeans)
 export(DA.ltt)
 export(DA.ltt2)
 export(DA.msf)
@@ -65,4 +69,3 @@ importFrom(lme4,glmer)
 importFrom(lme4,glmer.nb)
 importFrom(pROC,roc)
 importFrom(parallel,detectCores)
-importFrom(venneuler,venneuler)

R/DA.poi.R

@@ -89,7 +89,7 @@ DA.poi <- function(data, predictor, paired = NULL, covars = NULL, relative = TRU
       pois <- function(x){
         fit <- NULL
         tryCatch(
-          fit <- lme4::glmer(x ~ predictor + offset(log(log(libSize))) + (1|paired),family="poisson", ...), 
+          fit <- lme4::glmer(x ~ predictor + offset(log(libSize)) + (1|paired),family="poisson", ...), 
           error = function(x) fit <- NULL)
         if(!is.null(fit)) {
           if(nrow(coef(summary(fit))) > 1) {

R/posthocs.R

@@ -0,0 +1,254 @@
+#' Run drop1 on all features from DAtest results with allResults = TRUE
+#'
+#' Works on "zpo", "znb", "qpo", "neb", "poi". Non-paired "lrm", "llm", "llm2"
+#' @param results Output from a DA."test" function with allResults = TRUE
+#' @param test Which test to use to calculate p-values. See ?drop1 for details
+#' @param p.adj P-value adjustment method. See ?p.adjust for details
+#' @param ... Additional arguments for drop1 function
+#' @export
+DA.drop1 <- function(results, test = "Chisq", p.adj = "fdr", ...){
+
+  # Class
+  k <- 1
+  while(class(results[[k]])[1] == "NULL"){
+    k<- k+1
+  } 
+  xclass <- class(results[[k]])
+
+  if(!any(c("lm","glm","zeroinfl","negbin","glmerMod") %in% xclass)) stop(paste("Class should be one of lm, glm, zeroinfl, negbin or glmerMod and not:",xclass))
+
+  xres <- lapply(results, function(x) tryCatch(drop1(x, test = test, ...),error = function(e) NA))
+
+  if(all(xclass == "lm")){
+    AIC <- lapply(xres, function(x) x[,4])
+    pv <-  lapply(xres, function(x) x[,5])
+
+    AIC <- do.call(rbind,AIC[lapply(AIC, length) > 1])
+    pv <- do.call(rbind,pv[lapply(pv, length) > 1])
+    pva <- apply(pv, 2, function(x) p.adjust(x, method=p.adj))
+
+    colnames(AIC) <- paste0("AIC_",rownames(drop1(results[[1]], ...)))
+    colnames(pv) <- paste0("pval_",rownames(drop1(results[[1]], ...)))
+    colnames(pva) <- paste0("pval.adj_",rownames(drop1(results[[1]], ...)))
+
+    res <- cbind(AIC,pv,pva)
+
+  }
+
+  if(xclass[1] == "glm"){
+
+    if(is.na(results[[k]]$aic)){
+
+      Dev <- lapply(xres, function(x) x[,2])
+      pv <-  lapply(xres, function(x) x[,4])
+
+      Dev <- do.call(rbind,Dev[lapply(Dev, length) > 1])
+      pv <- do.call(rbind,pv[lapply(pv, length) > 1])
+      pva <- apply(pv, 2, function(x) p.adjust(x, method=p.adj))
+
+      colnames(Dev) <- paste0("Deviance_",rownames(drop1(results[[1]], ...)))
+      colnames(pv) <- paste0("pval_",rownames(drop1(results[[1]], ...)))
+      colnames(pva) <- paste0("pval.adj_",rownames(drop1(results[[1]], ...)))
+
+      res <- cbind(Dev,pv,pva)
+
+    } else {
+
+      AIC <- lapply(xres, function(x) x[,3])
+      LRT <- lapply(xres, function(x) x[,4])
+      pv <-  lapply(xres, function(x) x[,5])
+
+      AIC <- do.call(rbind,AIC[lapply(AIC, length) > 1])
+      LRT <- do.call(rbind,LRT[lapply(LRT, length) > 1])
+      pv <- do.call(rbind,pv[lapply(pv, length) > 1])
+      pva <- apply(pv, 2, function(x) p.adjust(x, method=p.adj))
+
+      colnames(AIC) <- paste0("AIC_",rownames(drop1(results[[1]], ...)))
+      colnames(LRT) <- paste0("LRT",rownames(drop1(results[[1]], ...)))
+      colnames(pv) <- paste0("pval_",rownames(drop1(results[[1]], ...)))
+      colnames(pva) <- paste0("pval.adj_",rownames(drop1(results[[1]], ...)))
+
+      res <- cbind(AIC,LRT,pv,pva)
+
+    }
+
+  }
+
+  if(all(xclass == "zeroinfl" | xclass == "glmerMod")){
+    AIC <- lapply(xres, function(x) x[,2])
+    LRT <- lapply(xres, function(x) x[,3])
+    pv <-  lapply(xres, function(x) x[,4])
+
+    AIC <- do.call(rbind,AIC[lapply(AIC, length) > 1])
+    LRT <- do.call(rbind,LRT[lapply(LRT, length) > 1])
+    pv <- do.call(rbind,pv[lapply(pv, length) > 1])
+    pva <- apply(pv, 2, function(x) p.adjust(x, method=p.adj))
+
+    colnames(AIC) <- paste0("AIC_",rownames(drop1(results[[1]], ...)))
+    colnames(LRT) <- paste0("LRT",rownames(drop1(results[[1]], ...)))
+    colnames(pv) <- paste0("pval_",rownames(drop1(results[[1]], ...)))
+    colnames(pva) <- paste0("pval.adj_",rownames(drop1(results[[1]], ...)))
+
+    res <- cbind(AIC,LRT,pv,pva)
+
+  }
+
+  if(xclass[1] == "negbin"){
+    AIC <- lapply(xres, function(x) x[,3])
+    LRT <- lapply(xres, function(x) x[,4])
+    pv <-  lapply(xres, function(x) x[,5])
+
+    AIC <- do.call(rbind,AIC[lapply(AIC, length) > 1])
+    LRT <- do.call(rbind,LRT[lapply(LRT, length) > 1])
+    pv <- do.call(rbind,pv[lapply(pv, length) > 1])
+    pva <- apply(pv, 2, function(x) p.adjust(x, method=p.adj))
+
+    colnames(AIC) <- paste0("AIC_",rownames(drop1(results[[k]], ...)))
+    colnames(LRT) <- paste0("LRT",rownames(drop1(results[[k]], ...)))
+    colnames(pv) <- paste0("pval_",rownames(drop1(results[[k]], ...)))
+    colnames(pva) <- paste0("pval.adj_",rownames(drop1(results[[k]], ...)))
+
+    res <- cbind(AIC,LRT,pv,pva)
+
+  }
+
+  return(res)
+}
+
+
+#' Run anova on all features from DAtest results with allResults = TRUE
+#'
+#' Works on "lrm", "llm", "llm2". Non-paired "neb"
+#' @param results Output from a DA."test" function with allResults = TRUE
+#' @param p.adj P-value adjustment method. See ?p.adjust for details
+#' @param ... Additional arguments for anova function
+#' @export
+DA.anova <- function(results, p.adj = "fdr", ...){
+
+  # Class
+  k <- 1
+  while(class(results[[k]])[1] == "NULL"){
+    k<- k+1
+  } 
+  xclass <- class(results[[k]])
+
+  if(!any(c("lm","nebgin","lme") %in% xclass)) stop(paste("Class should be one of lm, lme or negbin and not:",xclass))
+
+  if(all(xclass == "lme")){
+
+    pv <-  lapply(results, function(x) tryCatch(anova(x, ...)[,4],error = function(e) NA))
+
+    pv <- do.call(rbind,pv[lapply(pv, length) > 1])
+    pva <- apply(pv, 2, function(x) p.adjust(x, method=p.adj))
+
+    colnames(pv) <- paste0("pval_",rownames(anova(results[[1]], ...)))
+    colnames(pva) <- paste0("pval.adj_",rownames(anova(results[[1]], ...)))
+
+    res <- cbind(pv,pva)
+
+  }
+
+  if(xclass[1] == "negbin" | xclass[1] == "lm"){
+
+    pv <-  lapply(results, function(x) tryCatch(anova(x)[,5],error = function(e) NA))
+
+    pv <- do.call(rbind,pv[lapply(pv, length) > 1])
+    pva <- apply(pv, 2, function(x) p.adjust(x, method=p.adj))
+
+    colnames(pv) <- paste0("pval_",rownames(anova(results[[1]])))
+    colnames(pva) <- paste0("pval.adj_",rownames(anova(results[[1]])))
+
+    res <- cbind(pv,pva)
+
+  }
+
+
+  return(res)
+}
+
+
+#' Run TukeyHSD on all features from DAtest results with allResults = TRUE
+#'
+#' Works on "aov", "lao", "lao2"
+#' @param results Output from a DA."test" function with allResults = TRUE
+#' @param variable Which variable to test. Default predictor. Alternatively, the name of a covar
+#' @param p.adj P-value adjustment method. See ?p.adjust for details
+#' @param ... Additional arguments for TukeyHSD function
+#' @export
+DA.TukeyHSD <- function(results, variable = "predictor", p.adj = "fdr", ...){
+
+  # Class
+  k <- 1
+  while(class(results[[k]])[1] == "NULL"){
+    k<- k+1
+  } 
+  xclass <- class(results[[k]])
+
+  if(xclass[1] != "aov") stop(paste("Class should be aov and not:",xclass))
+  if(!variable %in% attr(results[[k]]$terms,"term.labels")) stop(paste(variable,"not found in the models."))
+
+  pv <-  lapply(results, function(x) tryCatch(as.data.frame(TukeyHSD(x, ...)[variable])[,4],error = function(e) NA))
+
+  pvs <- do.call(rbind,pv[lapply(pv, length) > 1])
+  colnames(pvs) <- rownames(as.data.frame(TukeyHSD(results[[k]], ...)[variable]))
+  pva <- apply(pvs, 2, function(x) p.adjust(x, method=p.adj))
+
+  colnames(pvs) <- paste0("pval_",colnames(pvs))
+  colnames(pva) <- paste0("pval.adj_",colnames(pva))
+
+  res <- cbind(pvs,pva)
+
+  return(res)
+}
+
+
+#' Run lsmeans on all features from DAtest results with allResults = TRUE
+#'
+#' Pairwise testing on predictor and covars. Works on "poi", "neb", "lrm", "llm", "llm2", "qpo", "znb", "zpo".
+#' 
+#' Require the lsmeans package
+#' @param results Output from a DA."test" function with allResults = TRUE
+#' @param variable Which variable to test. Default predictor. Alternatively, the name of a covar
+#' @param predictor If results come from a paired "lrm", "llm" or "llm2" supply the original predictor variable in the form of as a vector
+#' @param covars If results come from a paired "lrm", "llm" or "llm2" supply the original covars in the form of a named list
+#' @param p.adj P-value adjustment method. See ?p.adjust for details
+#' @param ... Additional arguments for lsmeans function
+#' @export
+DA.lsmeans <- function(results, variable = "predictor", predictor = NULL, covars = NULL, p.adj = "fdr", ...){
+
+  library(lsmeans)
+
+  # Class
+  k <- 1
+  while(class(results[[k]])[1] == "NULL"){
+    k<- k+1
+  } 
+
+  if(class(results[[k]])[1] == "lme"){
+    form <<- as.formula(paste("x ~",paste(attr(results[[1]]$terms,"term.labels"), collapse = "+")))
+    if(is.null(predictor)) stop("predictor has to be supplied for a paired lrm, llm and llm2")
+    if(!is.null(covars)){
+      for(i in 1:length(covars)){
+        assign(names(covars)[i],covars[[i]])
+      }
+    }
+  }
+
+  mc <- lapply(results, function(x) tryCatch(summary(pairs(lsmeans(x, variable))),error = function(e) NA))
+  pv <- lapply(mc, function(x) as.data.frame(x)$p.value)
+  est <- lapply(mc, function(x) as.data.frame(x)$estimate)
+
+  pvs <- do.call(rbind,pv[lapply(pv, length) > 1])
+  est <- do.call(rbind,est[lapply(est, length) > 1])
+  colnames(pvs) <- summary(pairs(lsmeans(results[[k]], variable)))$contrast
+  pva <- apply(pvs, 2, function(x) p.adjust(x, method=p.adj))
+
+  colnames(est) <- paste0("estimate_",colnames(pvs))
+  colnames(pvs) <- paste0("pval_",colnames(pvs))
+  colnames(pva) <- paste0("pval.adj_",colnames(pva))
+
+  res <- cbind(est,pvs,pva)
+  if(class(results[[k]])[1] == "lme") rm(form, envir = .GlobalEnv)
+  return(res)
+}
+

R/vennDA.R

@@ -1,15 +1,18 @@
 #' Plot Venn diagram from allDA object
 #'
-#' Plot a Venn (Euler) diagram of features found by different methods
+#' Plot a Venn (Euler) diagram of features found by different methods.
+#' 
+#' Require the venneuler package.
 #' @param x Output from the allDA function
 #' @param tests Character vector with tests to plot (E.g. c("ttt","adx.t","wil"), see names(x$results)). Default none
 #' @param alpha Numeric. q-value threshold for significant features. Default 0.05
 #' @param ... Additional arguments for plotting
 #' @return Nothing
-#' @importFrom  venneuler venneuler
 #' @export
 vennDA <- function(x, tests = NULL, alpha = 0.05, ...){
 
+  library(venneuler)
+
   if(!all(names(x) == c("table","results"))) stop("x is not an allDA object")
 
   plottests <- tests[tests %in% names(x[[2]])]  

R/zzz.R

@@ -1,6 +1,6 @@
 
 .onLoad <- function(libname, pkgname){
-  message("DAtest version 2.6.1")
+  message("DAtest version 2.6.2")
 
   # Fix samr problem
   if(.Platform$OS.type == "windows"){