From 0be90c2516f0975f793ac19512263c580639d4b8 Mon Sep 17 00:00:00 2001
From: Jakob Russel <cmf483@science.domain>
Date: Thu, 2 Jan 2020 14:49:27 +0100
Subject: [PATCH] add functions

---
 NAMESPACE             |  4 ++++
 R/adonis_OmegaSq.R    | 29 +++++++++++++++++++++++
 R/rcurve.R            | 54 +++++++++++++++++++++++++++++++++++++++++++
 README.md             |  8 +++++++
 man/adonis_OmegaSq.Rd | 19 +++++++++++++++
 man/rcurve.Rd         | 23 ++++++++++++++++++
 6 files changed, 137 insertions(+)
 create mode 100644 R/adonis_OmegaSq.R
 create mode 100644 R/rcurve.R
 create mode 100644 man/adonis_OmegaSq.Rd
 create mode 100644 man/rcurve.Rd

diff --git a/NAMESPACE b/NAMESPACE
index b9ab347..4af906a 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/adonis_OmegaSq.R b/R/adonis_OmegaSq.R
new file mode 100644
index 0000000..97ba99f
--- /dev/null
+++ b/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)
+}
diff --git a/R/rcurve.R b/R/rcurve.R
new file mode 100644
index 0000000..a962775
--- /dev/null
+++ b/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)
+    }
+    
+}
diff --git a/README.md b/README.md
index 9ad5d10..adbb802 100644
--- a/README.md
+++ b/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. 
diff --git a/man/adonis_OmegaSq.Rd b/man/adonis_OmegaSq.Rd
new file mode 100644
index 0000000..208b12e
--- /dev/null
+++ b/man/adonis_OmegaSq.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/adonis_OmegaSq.R
+\name{adonis_OmegaSq}
+\alias{adonis_OmegaSq}
+\title{Calculate (partial) Omega-squared (effect-size calculation) for PERMANOVA and add it to the input object}
+\usage{
+adonis_OmegaSq(adonisOutput, partial = TRUE)
+}
+\arguments{
+\item{adonisOutput}{An adonis object}
+
+\item{partial}{Should partial omega-squared be calculated (sample size adjusted). Default TRUE}
+}
+\value{
+Original adonis object with the (partial) Omega-squared values added
+}
+\description{
+Calculate (partial) Omega-squared (effect-size calculation) for PERMANOVA and add it to the input object
+}
diff --git a/man/rcurve.Rd b/man/rcurve.Rd
new file mode 100644
index 0000000..81ce9ee
--- /dev/null
+++ b/man/rcurve.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rcurve.R
+\name{rcurve}
+\alias{rcurve}
+\title{Rarefaction curve on phyloseq object
+Theoretical richness with vegan's rarefy function}
+\usage{
+rcurve(physeq, subsamp = 10^c(1:5), trim = TRUE,
+  add_sample_data = TRUE)
+}
+\arguments{
+\item{physeq}{phyloseq object}
+
+\item{subsamp}{Vector of number of reads to subsample}
+
+\item{trim}{Remove richness estimations from subsamples larger than the library size}
+
+\item{add_sample_data}{Add sample data to data.frame}
+}
+\description{
+Rarefaction curve on phyloseq object
+Theoretical richness with vegan's rarefy function
+}