add functions

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(UniFrac.multi)
+export(adonis_OmegaSq)
 export(comdist.par)
 export(comdistnt.par)
 export(community_rrna)
@@ -10,6 +11,7 @@ export(proportionality)
 export(rarefy_rrna)
 export(rarefy_rrna.matrix)
 export(rarefy_rrna.phyloseq)
+export(rcurve)
 export(ses.UniFrac)
 export(ses.comdist)
 export(ses.comdist2)
@@ -25,11 +27,13 @@ import(foreach)
 import(parallel)
 import(phyloseq)
 import(picante)
+import(utils)
 import(vegan)
 importFrom(abind,abind)
 importFrom(stats,as.dist)
 importFrom(stats,cov)
 importFrom(stats,dnorm)
+importFrom(stats,na.omit)
 importFrom(stats,p.adjust)
 importFrom(stats,pbeta)
 importFrom(stats,quantile)

R/adonis_OmegaSq.R

@@ -0,0 +1,29 @@
+#' Calculate (partial) Omega-squared (effect-size calculation) for PERMANOVA and add it to the input object
+#' 
+#' @param adonisOutput An adonis object
+#' @param partial Should partial omega-squared be calculated (sample size adjusted). Default TRUE 
+#' @return Original adonis object with the (partial) Omega-squared values added
+#' @import vegan
+#' @export
+adonis_OmegaSq <- function(adonisOutput, partial = TRUE){
+    if(!is(adonisOutput, "adonis")) stop("Input should be an adonis object")
+    aov_tab <- adonisOutput$aov.tab
+    heading <- attr(aov_tab, "heading")
+    MS_res <- aov_tab[rownames(aov_tab) == "Residuals", "MeanSqs"]
+    SS_tot <- aov_tab[rownames(aov_tab) == "Total", "SumsOfSqs"]
+    N <- aov_tab[rownames(aov_tab) == "Total", "Df"] + 1
+    if(partial){
+        omega <- apply(aov_tab, 1, function(x) (x["Df"]*(x["MeanSqs"]-MS_res))/(x["Df"]*x["MeanSqs"]+(N-x["Df"])*MS_res))
+        aov_tab$parOmegaSq <- c(omega[1:(length(omega)-2)], NA, NA)
+        cn_order <- c("Df", "SumsOfSqs", "MeanSqs", "F.Model", "R2", "parOmegaSq", "Pr(>F)")
+    } else {
+        omega <- apply(aov_tab, 1, function(x) (x["SumsOfSqs"]-x["Df"]*MS_res)/(SS_tot+MS_res))
+        aov_tab$OmegaSq <- c(omega[1:(length(omega)-2)], NA, NA)
+        cn_order <- c("Df", "SumsOfSqs", "MeanSqs", "F.Model", "R2", "OmegaSq", "Pr(>F)")
+    }
+    aov_tab <- aov_tab[, cn_order]
+    attr(aov_tab, "names") <- cn_order
+    attr(aov_tab, "heading") <- heading
+    adonisOutput$aov.tab <- aov_tab
+    return(adonisOutput)
+}

R/rcurve.R

@@ -0,0 +1,54 @@
+#' Rarefaction curve on phyloseq object
+#' Theoretical richness with vegan's rarefy function
+#' 
+#' @param physeq phyloseq object
+#' @param subsamp Vector of number of reads to subsample
+#' @param trim Remove richness estimations from subsamples larger than the library size
+#' @param add_sample_data Add sample data to data.frame
+#' @import vegan
+#' @import utils
+#' @importFrom stats na.omit
+#' @export
+rcurve <- function(physeq, subsamp = 10^c(1:5), trim = TRUE, add_sample_data = TRUE){
+
+    otu <- otu_table(physeq)
+    if(!taxa_are_rows(physeq)){
+        otu <- t(otu)
+    }
+    colS <- colSums(otu)
+
+    pb <- txtProgressBar(min = 0, max = length(subsamp), style = 3)
+    rars <- list()
+    for(i in seq_along(subsamp)){
+        setTxtProgressBar(pb, i)
+        rars[[i]] <- vegan::rarefy(otu, sample = subsamp[i], MARGIN = 2)
+    }
+    mat <- do.call(cbind, rars)
+
+    if(trim){
+        mat_bool <- sapply(subsamp, function(i) sapply(colS, function(j) i <= j))
+        mat_new <- mat * mat_bool
+        mat_new[mat_new == 0] <- NA
+    } else {
+        mat_new <- mat
+    }
+    colnames(mat_new) <- subsamp
+
+    # To a data.frame
+    df <- as.data.frame.table(mat_new)
+    colnames(df) <- c("Sample", "Reads", "Richness")
+
+    if(trim){
+        df <- na.omit(df)
+    }
+
+    # Add sample data
+    if(add_sample_data){
+        samp <- sample_data(physeq)
+        df2 <- merge(df, samp, by = "Sample", by.y = "row.names")
+        return(df2)
+    } else {
+        return(df)
+    }
+
+}

README.md

@@ -12,6 +12,14 @@ MicEco: Various functions for analysis for microbial community data
 ### Citation
 [![DOI](https://zenodo.org/badge/83547545.svg)](https://zenodo.org/badge/latestdoi/83547545)
 
+#### adonis_OmegaSq
+
+Calculate the unbiased effect size estimation (partial) omega-squared for adonis (PERMANOVA) models
+
+#### rcurve
+
+Rarefaction curve (theoretical and fast) from a phyloseq object. Output ready for plotting in ggplot2
+
 #### UniFrac.multi
 
 With unrooted phylogenies UniFrac sets the root randomly on the tree.