From 4d3f3b33ba75e39f78e1c912948466ef61d9780f Mon Sep 17 00:00:00 2001
From: "Thomas A. Christensen II" <25492070+MillironX@users.noreply.github.com>
Date: Thu, 16 May 2019 23:06:06 -0600
Subject: [PATCH] Added automated rarefaction selection

---
 main.sh                  | 22 ++++++++++++++++------
 manipulatefeaturetable.R | 22 +++++++++++++++++-----
 2 files changed, 33 insertions(+), 11 deletions(-)

diff --git a/main.sh b/main.sh
index b8d5f8e..1c12951 100644
--- a/main.sh
+++ b/main.sh
@@ -44,12 +44,13 @@ echo "--^-- X: Finding metadata columns..."
 Rscript ./fetchmetadata.R
 echo "--^-- X: Finding metadata columns...Done!"
 
+# Our minimum taxa count is 11123 - this will be needed for rarefaction
+MINRAREFACTION=$(<rarefaction.min.txt)
+MAXRAREFACTION=$(<rarefaction.max.txt)
+
 # Start up QIIME
 source activate qiime2
 
-# Our minimum taxa count is 11123 - this will be needed for rarefaction
-RAREFACTION=11123
-
 # Convert the feature table into BIOM format
 echo "--^-- X: Importing data..."
 biom convert \
@@ -79,7 +80,7 @@ rm -r "core-metrics-results"
 # This is one of the few QIIME commands that can use multithreading
 qiime diversity core-metrics \
  --i-table feature-table.qza \
- --p-sampling-depth "$RAREFACTION" \
+ --p-sampling-depth "$MINRAREFACTION" \
  --m-metadata-file metadata.tsv \
  --p-n-jobs 4 \
  --output-dir core-metrics-results \
@@ -99,8 +100,17 @@ qiime taxa barplot \
  --o-visualization visualizations/barplot.qzv
 echo "--^-- X: Generating barplot...Done!"
 
- # Run alpha-diversity group significance: this will automatically include all the columns
- # Evenness first
+echo "--^-- X: Plotting rarefaction curve..."
+# Create a rarefaction curve to make sure the magic of rarefaction is valid
+qiime diversity alpha-rarefaction \
+ --i-table feature-table.qza \
+ --p-maxdepth "$MAXRAREFACTION" \
+ --m-metadata-file metadata.tsv \
+ --o-visualization visualizations/rarefaction-curve.qzv
+echo "--^-- X: Plotting rarefaction curve...Done!"
+
+# Run alpha-diversity group significance: this will automatically include all the columns
+# Evenness first
 echo "--^-- X: Finding alpha-group-significance..."
 qiime diversity alpha-group-significance \
  --i-alpha-diversity core-metrics-results/evenness_vector.qza \
diff --git a/manipulatefeaturetable.R b/manipulatefeaturetable.R
index 64c09c0..6dc6000 100644
--- a/manipulatefeaturetable.R
+++ b/manipulatefeaturetable.R
@@ -1,10 +1,10 @@
 # Read the inital feature table in
 feature_table <- read.table("metaxa-feature-table.tsv",
-							header = TRUE,
-                            sep = "\t",
-                            quote = "",
-                            strip.white = TRUE,
-							check.names = FALSE)
+							              header=TRUE,
+                            sep="\t",
+                            quote="",
+                            strip.white=TRUE,
+							              check.names=FALSE)
 
 # Get the dimensions of the table
 numSamples <- ncol(feature_table) - 1
@@ -40,6 +40,18 @@ for (i in 1:numFeatures){
 feature_table$'#SampleId' <- ids
 feature_table <- feature_table[c(numSamples+2, 1:(numSamples+1))]
 
+# Find minimum and maximum rarefaction values
+numCounts <- vector(length=numSamples)
+for (i in 2:(numSamples+1)) {
+  numCounts[i-1] <- sum(feature_table[,i])
+}
+minRarefaction <- min(numCounts)
+maxRarefaction <- max(numCounts)
+write.table(minRarefaction, "rarefaction.min.txt",
+            row.names=FALSE, col.names=FALSE, quote=FALSE)
+write.table(maxRarefaction, "rarefaction.max.txt",
+            row.names=FALSE, col.names=FALSE, quote=FALSE)
+
 # Write the file out
 write.table(feature_table,
             file="feature-table.tsv",