scripts.mainFigures.Rmd

---
title: "scripts.mainFigures"
author: "Tuesday Simmons and Alex Styer"
date: "2/20/2020"
output:
  html_document: default
  pdf_document: default
---

The following sripts should run as is assuming you've downloaded all data and markdown files to the same directory. 

```{r, load packages + data + set functions, warning=FALSE}
#load all packages
library("phyloseq")
library("ggplot2")
library("scales")
library("ape")
library("vegan")
library("data.table")
library("RColorBrewer")
library("colorRamps")
library("VennDiagram")
library("beepr")
library("agricolae")
library("seqinr")
library("patchwork")
library("labdsv")
library("reshape2")
library("dplyr")

#import all data
millet_base <- import_biom("feature-table.biom")
beep(2)
mapping <- import_qiime_sample_data("metadata.txt")
millet_base <- merge_phyloseq(millet_base, mapping)
millet_base
#rename columns in taxonomy table
colnames(tax_table(millet_base)) <- c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species")
colnames(tax_table(millet_base))

#a function to get indicator species for levels of a given column in your metadata file
#takes a phyloseq object and returns a dataframe of indicators
get_indicators <- function(physeq, groupfactor){
all.data <- psmelt(physeq)
all.data <- subset(all.data, Abundance > 0)
metadata <- sample_data(physeq)

spec <- subset(all.data, select = c("Sample", "OTU", "Abundance"))
spec <- dcast(spec, Sample ~ OTU, value.var =  "Abundance")
spec <- mutate_at(spec, c(2:ncol(spec)), ~replace(., is.na(.), 0))
rownames(spec) <- spec$Sample
spec <- subset(spec, select = -c(Sample))
spec <- spec[order(rownames(spec)),]

clust <- as.data.frame(metadata[[groupfactor]])
colnames(clust) <- "group"
clust$group.level <- as.integer(clust$group)
rownames(clust) <- metadata$SampleID
clust <- subset(clust, rownames(clust) %in% rownames(spec))
clust <- clust[order(rownames(clust)),]

if (!(identical(rownames(spec), rownames(clust)))) {
  stop("Species and clustering dataframes are misaligned! Check sample orders.")
}

indvals <- indval(spec, clust$group.level, numitr = 1000)
summary(indvals)

#convert all.indvals to a single df for ease of access
pval <- as.data.frame(indvals$pval)
pval$ASV <- rownames(pval)
maxcls <- as.data.frame(indvals$maxcls)
maxcls$ASV <- rownames(maxcls)
indcls <- as.data.frame(indvals$indcls)
indcls$ASV <- rownames(indcls)
all.indvals <- merge(pval, maxcls, by = "ASV")
all.indvals <- merge(all.indvals, indcls, by = "ASV")
colnames(all.indvals) <- c("ASV", "pval", "group.level", "indcls")
all.indvals <- merge(all.indvals, clust, by = "group.level")
all.indvals <- subset(all.indvals, select = c(ASV, group, pval, indcls))
all.indvals <- unique(all.indvals)

#get only significant indicators (p<=0.05) with an indcls >=0.5
sig.indvals <- subset(all.indvals, pval <= 0.05)
sig.indvals <- subset(sig.indvals, indcls >= 0.5)

#return an object with both the subsetted significant indvals and all indvals
indvals <- list("all.indvals" = all.indvals, "sig.indvals" = sig.indvals, "all.data" = indvals)
}

#a function to make a fasta file of just indicator species
#requires original repseqs fasta file and a list of indicators generated by get_indicators
#NOTE--this writes the fasta file out to your working directory
make_fasta <- function(indicator.list, repseqs, outfile){
  #subset repseqs file to indicators only
  indicators <- repseqs[names(repseqs) %in% unique(indicator.list$ASV)]
  #write it out to file to run through fasttree
  indicator.names <- names(indicators)
  indicator.sequences <- getSequence(indicators)
  write.fasta(indicator.sequences, indicator.names, outfile)
}

#make a clean taxonomy table
clean_taxtable <- function(physeq){
  tax.table <- as.data.frame(tax_table(physeq))
  colnames(tax.table) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
  tax.table$Kingdom <- as.character(lapply(tax.table$Kingdom, sub, pattern = "k__", replacement = ""))
  tax.table$Phylum <- as.character(lapply(tax.table$Phylum, sub, pattern = "p__", replacement = ""))
  tax.table$Class <- as.character(lapply(tax.table$Class, sub, pattern = "c__", replacement = ""))
  tax.table$Order <- as.character(lapply(tax.table$Order, sub, pattern = "o__", replacement = ""))
  tax.table$Family <- as.character(lapply(tax.table$Family, sub, pattern = "f__", replacement = ""))
  tax.table$Genus <- as.character(lapply(tax.table$Genus, sub, pattern = "g__", replacement = ""))
  tax.table$Species <- as.character(lapply(tax.table$Species, sub, pattern = "s__", replacement = ""))
  tax.table$Kingdom <- ifelse(tax.table$Kingdom == "", "Unknown", as.character(tax.table$Kingdom))
  tax.table$Phylum <- ifelse(tax.table$Phylum == "", "Unknown", as.character(tax.table$Phylum))
  tax.table$Class <- ifelse(tax.table$Class == "", "Unknown", as.character(tax.table$Class))
  tax.table$Order <- ifelse(tax.table$Order == "", "Unknown", as.character(tax.table$Order))
  tax.table$Family <- ifelse(tax.table$Family == "", "Unknown", as.character(tax.table$Family))
  tax.table$Genus <- ifelse(tax.table$Genus == "", "Unknown", as.character(tax.table$Genus))
  tax.table$Species <- ifelse(tax.table$Species == "", "Unknown", as.character(tax.table$Species))
  tax.table$Kingdom <- ifelse(is.na(tax.table$Kingdom) == TRUE , "Unknown", as.character(tax.table$Kingdom))
  tax.table$Phylum <- ifelse(is.na(tax.table$Phylum) == TRUE , "Unknown", as.character(tax.table$Phylum))
  tax.table$Class <- ifelse(is.na(tax.table$Class) == TRUE , "Unknown", as.character(tax.table$Class))
  tax.table$Order <- ifelse(is.na(tax.table$Order) == TRUE , "Unknown", as.character(tax.table$Order))
  tax.table$Family <- ifelse(is.na(tax.table$Family) == TRUE , "Unknown", as.character(tax.table$Family))
  tax.table$Genus <- ifelse(is.na(tax.table$Genus) == TRUE , "Unknown", as.character(tax.table$Genus))
  tax.table$Species <- ifelse(is.na(tax.table$Species) == TRUE , "Unknown", as.character(tax.table$Species))
  tax.table$ASV <- rownames(tax.table)
  rownames(tax.table) <- NULL
  tax.table <- tax.table
}

theme_set(theme_bw())
#set up color palette
farrowAndBall_palette <- c(
  "#4d5b6a" #Stiffkey Blue
  ,"#6a90b4" #Cook's Blue
  ,"#599ec4" #ST Giles Blue
  ,"#a1c5c8" #Blue Ground
  ,"#7997a1" #Stone Blue
  ,"#427e83" #Vardo
  ,"#84b59c" #Arsenic
  ,"#919f70" #Yeabridge Green
  ,"#686a47" #Bancha
  ,"#c8bd83" #Churlish Green
  ,"#cb9e59" #India Yellow
  ,"#ecc363" #Babouche
  ,"#c57b67" #Red Earth
  ,"#d65f3d" #Charlotte's Locks
  ,"#a04344" #Incarnadine
  ,"#bf7a8f" #Rangwali
  ,"#8d8089" #Brassica
  ,"#50414c" #Pelt
  ,"#e5e0db" #Strong White
  ,"#444546" #Off-Black
)
farrowAndBall_pal <- function() scales::manual_pal(farrowAndBall_palette)
scale_fill_farrowAndBall <- function(...) ggplot2::discrete_scale('fill', 'farrowAndBall', farrowAndBall_pal(), ...)

```

```{r, Figure 1, include=TRUE}
#subset main data file to only include samples from the multi-species field experiment
milletfield <- subset_samples(millet_base, Experiment=="Drought Levels")
#remove any metadata columns that do not pertain to any samples in the subsetted data
milletfield@sam_data <- subset(milletfield@sam_data, select=-c(SideTreatment, RootType, Divison, SubDivision, Timepoint, Plant, Root, Notes))
###filtering
#Removed ASVs without at least 1 read in at least 1 sample
mf_filter <- filter_taxa(milletfield, function(x) sum(x >= 1) >= 1, TRUE)
mf_filter
#Removed samples with no reads
mf_filter <- prune_samples(sample_sums(mf_filter) > 0, mf_filter)
mf_filter #removed 3 samples
#set order for SampleType
mf_filter@sam_data$SampleType <- factor(mf_filter@sam_data$SampleType,levels=c("Soil","Rhizo","Root"))

#################################
########### FIGURE 1A ###########
#################################
richness = estimate_richness(mf_filter,measure="Shannon")
s <- data.frame(sample_data(mf_filter))
alphadiv <- cbind(richness, s)
alphadiv <- data.frame(alphadiv)
alphadiv$Species <- factor(alphadiv$Species,levels=c("Foxtail Millet","Pearl Millet","Proso Millet","Japanese Millet"))
alphadiv$Treatment <- factor(alphadiv$Treatment,levels=c("Severe Drought","Moderate Drought","Control"))
alphadiv$SampleType <- factor(alphadiv$SampleType,levels=c("Soil","Rhizo","Root"))

#statistical difference between categories
interac <- with(alphadiv, interaction(SampleType, Treatment))
anova_interac <- aov(data=alphadiv, formula = Shannon~interac)
summary(anova_interac)
out <- HSD.test(anova_interac, "interac", group=TRUE)
plot(out)

#set the order of the three treatments along the x axis
alphadiv$Treatment <- factor(alphadiv$Treatment, levels=c("Control","Moderate Drought","Severe Drought"))

#taking the stats data from the out table and putting it into a dataframe
vars <- data.frame(expand.grid(levels(mf_filter@sam_data$Treatment),levels(mf_filter@sam_data$SampleType)))
colnames(vars) <- c("Treatment","SampleType")
dat <- data.frame(x= vars$Treatment, y = c(6.7, 6.9, 6.6, 6.7, 6.3,6.2, 5.8, 5.5, 5.2), vars, labs=c("a","a","a","ab","ab","bc","c","c","c"))

alpha1 <- ggplot(data=alphadiv,aes(x=Treatment,y=Shannon)) + 
  geom_boxplot(aes(fill=SampleType)) + 
  facet_grid(~SampleType,scales="free",space="free",
             labeller=labeller(SampleType=c("Soil"="Bulk Soil","Rhizo"="Rhizosphere","Root"="Root Endosphere"))) +
  scale_fill_manual(breaks = c("Soil","Rhizo","Root"),
                    values=c("#c57b67","#c8bd83","#84b59c"))+
  guides(fill=FALSE) +
  ylim(3,7) +
  geom_text(data = dat, size = 5, aes(x = x, y = y, label = labs)) +
  scale_x_discrete("Treatment", breaks = c("Control","Moderate Drought","Severe Drought"),
                   labels=c("Control","Moderate\nDrought","Severe\nDrought")) +
  theme(axis.text.x=element_text(size=14,color="black",angle=45,hjust=1.0,vjust=1.0), 
        axis.text.y=element_text(size=16,color="black"), 
        axis.title=element_text(size=16,face="bold"),
        text=element_text(size=16)) 
alpha1

#################################
########### FIGURE 1B ###########
#################################
#Transform
mf_trans  = transform_sample_counts(mf_filter, function(x) 100 * x / sum(x) )
#Agglomerate taxa
mf_glomClass <- tax_glom(mf_trans, taxrank = "Class")
mf_glomClass@tax_table[is.na(mf_glomClass@tax_table)] <- "Unknown"
#Need to check number of classes for next step
mf_glomClass #114 Classes

#Class level sort based on taxa abundance
## By Class
#label taxa that are NOT in the top 9 as "other"
mf_othersClass <- names(sort(taxa_sums(mf_glomClass), decreasing = TRUE)[10:114])
mf_merge <- merge_taxa(mf_glomClass, mf_othersClass,1)
mf_meltClass <- psmelt(mf_merge)
#convert Class to a character vector from a factor 
str(mf_meltClass$Class)
mf_meltClass$Class <- as.character(mf_meltClass$Class)
#replace all NAs with Other
mf_meltClass$Class[is.na(mf_meltClass$Class)] <- "Other"
# convert Class back to a factor
mf_meltClass$Class <- as.factor(mf_meltClass$Class)
levels(factor(mf_meltClass$Class))
#sort taxa
mf_meltClass$Class <- factor(mf_meltClass$Class, levels = c("Other", 
                                                            "c__Deltaproteobacteria" ,
                                                            "c__[Saprospirae]"  ,      
                                                            "c__Acidobacteriia",
                                                            "c__Alphaproteobacteria" , 
                                                            "c__Bacilli",
                                                            "c__Betaproteobacteria" ,  
                                                            "c__Cytophagia"   ,        
                                                            "c__Gammaproteobacteria",
                                                            "c__Actinobacteria" ))
#sort sample type
mf_meltClass$SampleType <- factor(mf_meltClass$SampleType, levels = c("Soil", "Rhizo","Root"))
#renaming sample types in the facet names
sampletype_names <- c(
  'Soil'='Bulk\nSoil',
  'Rhizo'="Rhizosphere",
  'Root'='Root\nEndophyte'
)

rel1 <- ggplot(mf_meltClass, aes(x = Treatment, y = Abundance,  fill=Class, order=Class))+
  scale_fill_manual(values=farrowAndBall_palette)+
  geom_bar(stat = "identity", position = "fill") +
  facet_grid(. ~ SampleType,labeller=as_labeller(sampletype_names)) +
  theme(axis.text.x=element_text(size=14,color="black",angle=45,hjust=1.0), axis.text.y=element_text(size=16,color="black"), 
        axis.title=element_text(size=16,face="bold"),text=element_text(size=16),legend.title=element_blank())  
rel1

#################################
########### FIGURE 1C ###########
#################################
ord1c <- plot_ordination(mf_filter, ordinate(mf_filter, "MDS"), axes=1:2, color = "SampleType") + 
  geom_point(size = 4)+ 
  scale_color_manual(values=c("#c57b67","#c8bd83","#84b59c"),
                     name = "Sample Type",
                     breaks=c("Soil","Rhizo","Root"),
                     labels=c("Bulk Soil","Rhizosphere","Root Endosphere"))+
    theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
        axis.title=element_text(size=14),text=element_text(size=14))
ord1c

#################################
########### FIGURE 1D ###########
#################################
mf_filter_root <- subset_samples(mf_filter, SampleType=="Root")

ord1d <- plot_ordination(mf_filter_root, ordinate(mf_filter_root, "CAP", "bray", ~Treatment+Species), axes=1:2, color = "Treatment", shape="Species") + 
  geom_point(size = 4)+ 
  scale_color_manual(values=c("#84b59c","#c8bd83","#c57b67"))+
  theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
        axis.title=element_text(size=14),text=element_text(size=14))
ord1d

#################################
########### PERMANOVA ###########
#################################

#PERMANOVA for all samples
set.seed(1)
# Calculate bray curtis distance matrix
all_dist <- phyloseq::distance(mf_filter, method = "bray")
# make a data frame
all_data <- data.frame(sample_data(mf_filter))
# Then use an adonis test to look at your experimental variables.
adonis_all <- adonis(all_dist ~ SampleType + Treatment + Species, data = all_data)
adonis_all

## PERMANOVA for root samples
set.seed(1)
# Subset phyloseq object, calculate bray curtis distance matrix, make data frame
root_dist <- phyloseq::distance(mf_filter_root, method = "bray")
root_data <- data.frame(sample_data(mf_filter_root))
# Then use an adonis test to look at your experimental variables.
adonis_root <- adonis(root_dist ~ Treatment + Species, data = root_data)
adonis_root
# get test results as a table (necessary for next step)
adonis_root_tab <- adonis_root$aov.tab
adonis_root_tab_rows <- rownames(adonis_root_tab)
# Get a 'for-loop' to get your percent-variances attributable to each factor.
for(i in 1:4){ #1:4 is how many rows you have in table
  a1 <- c("The percent of variance attributable to the factor")
  b2 <- toString(adonis_root_tab_rows[i])
  c3 <- c("is")
  d4 <- toString(round(adonis_root_tab$SumsOfSqs[i]/adonis_root_tab$SumsOfSqs[4], 3)) #1st number should be "total" cell, 2nd number is how many digits it rounds to
  print(paste(a1,b2,c3,d4))
  
#################################
########### PLOT ALL ############
#################################
#requires library("patchwork")
Figure1 <- (alpha1 | ord1c) / (rel1 | ord1d)
#Plot this directly in console to be able to view the complete figure in the plots viewer
Figure1
}
```

```{r, Figure 2, include=TRUE}
#subset main data file to only include samples from the multi-species field experiment
sorghumfield <- subset_samples(millet_base, Experiment=="Spatial Organization")
#remove any metadata columns that do not pertain to any samples in the subsetted data
sorghumfield@sam_data <- subset(sorghumfield@sam_data, select=-c(Species, Latin, SideTreatment, RootType, Root, Notes))
###filtering
#Removed ASVs without at least 1 read in at least 1 sample
so_filter <- filter_taxa(sorghumfield, function(x) sum(x >= 1) >= 1, TRUE)
so_filter
#Removed samples with no reads
so_filter <- prune_samples(sample_sums(so_filter) > 0, so_filter)
so_filter #removed 1 sample

#set order for SampleType
so_filter@sam_data$SampleType <- factor(so_filter@sam_data$SampleType,levels=c("Soil","Rhizo","Root"))

#################################
########### FIGURE 2A ###########
#################################
richness2 = estimate_richness(so_filter,measure="Shannon")
s2 <- data.frame(sample_data(so_filter))
alphadiv2 <- cbind(richness2, s2)
alphadiv2 <- data.frame(alphadiv2)
alphadiv2$Treatment <- factor(alphadiv2$Treatment,levels=c("Control","Drought"))
alphadiv2$SampleType <- factor(alphadiv2$SampleType,levels=c("Soil","Rhizo","Root"))
alphadiv2$Divison <- factor(alphadiv2$Divison,levels=c("Soil","AllRhizo","SingleRhizo","AllRoot","SingleRoot"))
alphadiv2$SubDivision <- factor(alphadiv2$SubDivision,levels=c("Soil","AllRhizo","SingleRhizo","AllRoot","SingleRoot","SubRootTop","SubRootMiddle","SubRootBottom"))
alphadiv2$Timepoint <- as.character(alphadiv2$Timepoint)
alphadiv2$Timepoint <- factor(alphadiv2$Timepoint,levels=c("1","2","3","9","11"))

#statistical differnece between categories
interac2 <- with(alphadiv2, interaction(SubDivision))
anova_interac2 <- aov(data=alphadiv2, formula = Shannon~interac2)
summary(anova_interac2)
out2 <- HSD.test(anova_interac2, "interac2", group=TRUE)
plot(out2)

#taking the stats data from the out table and putting it into a dataframe
vars2 <- data.frame(expand.grid(levels(so_filter@sam_data$SubDivision)))
colnames(vars2) <- c("SubDivision")
vars2$SampleType <- c("Rhizo","Root","Rhizo","Root","Soil","Root","Root","Root")
out2[["groups"]]$SubDivision <- paste(row.names(out2[["groups"]]))
vars2 <- merge(vars2, out2[["groups"]], by="SubDivision")
dat2 <- data.frame(x= vars2$SubDivision, 
                   y = c(7,6.3,7.1,6.6,6.8,6.8,6.8,6.4), 
                   vars2)
dat2$SampleType <- factor(dat2$SampleType,levels=c("Soil","Rhizo","Root"))


alpha2 <- ggplot(data=alphadiv2,aes(x=SubDivision,y=Shannon)) + 
  geom_boxplot(aes(fill=SampleType)) + 
  geom_text(data = dat2, size = 5, aes(x = x, y = y, label = groups)) +
  facet_grid(~SampleType,
             scales="free",
             space="free",
             labeller=labeller(SampleType=c("Soil"="Bulk Soil","Rhizo"="Rhizosphere","Root"="Root Endosphere"))) +
  scale_fill_manual(breaks = c("Soil","Rhizo","Root"),
                    values=c("#c57b67","#c8bd83","#84b59c"))+
  guides(fill=FALSE) +
  ylim(3,7.5) +
  theme(axis.text.x=element_text(size=14,color="black",angle=45,hjust=1.0,vjust=1.0), 
        axis.text.y=element_text(size=16,color="black"), 
        axis.title=element_text(size=16,face="bold"),
        text=element_text(size=16)) 
alpha2

#################################
########### FIGURE 2B ###########
#################################
#subset to root samples only
so_filter_roots <- subset_samples(so_filter, SampleType=="Root")
#subset to the later time points only
so_f_r_late <- subset_samples(so_filter_roots, Timepoint%in%c("9","11"))
#subset to remove AllRoot samples
so_f_r_l_sub <- subset_samples(so_f_r_late, Divison!="AllRoot")

#Transform
so_trans_rlsub  = transform_sample_counts(so_f_r_l_sub, function(x) 100 * x / sum(x) )
#Agglomerate taxa
so_glomClass_rlsub <- tax_glom(so_trans_rlsub, taxrank = "Class")
so_glomClass_rlsub@tax_table[is.na(so_glomClass_rlsub@tax_table)] <- "Unknown"
#Need to check number of classes for next step
so_glomClass_rlsub #121 Classes

#Class level sort based on taxa abundance
#label taxa that are NOT in the top 9 as "other"
so_othersClass_rlsub <- names(sort(taxa_sums(so_glomClass_rlsub), decreasing = TRUE)[10:121])
so_merge <- merge_taxa(so_glomClass_rlsub, so_othersClass_rlsub,1)
so_meltClass_rlsub <- psmelt(so_merge)
#convert Class to a character vector from a factor 
str(so_meltClass_rlsub$Class)
so_meltClass_rlsub$Class <- as.character(so_meltClass_rlsub$Class)
#replace all NAs with Other
so_meltClass_rlsub$Class[is.na(so_meltClass_rlsub$Class)] <- "Other"
# convert Class back to a factor
so_meltClass_rlsub$Class <- as.factor(so_meltClass_rlsub$Class)
levels(factor(so_meltClass_rlsub$Class))
#sort taxa
so_meltClass_rlsub$Class <- factor(so_meltClass_rlsub$Class, levels = c("Other", 
                                                            "c__Gemmatimonadetes" ,
                                                            "c__[Saprospirae]"  ,      
                                                            "c__Acidobacteriia",
                                                            "c__Alphaproteobacteria" , 
                                                            "c__Bacilli",
                                                            "c__Betaproteobacteria" ,
                                                            "c__Deltaproteobacteria",
                                                            "c__Thermoleophilia"   ,        
                                                            "c__Gammaproteobacteria",
                                                            "c__Actinobacteria" ))
#sort sample type
so_meltClass_rlsub$SubDivision <- factor(so_meltClass_rlsub$SubDivision, levels = c("AllRoot","SingleRoot","SubRootTop","SubRootMiddle","SubRootBottom"))

rel2 <- ggplot(so_meltClass_rlsub, aes(x = Treatment, y = Abundance,  fill=Class, order=Class))+
  scale_fill_manual(values=farrowAndBall_palette)+
  geom_bar(stat = "identity", position = "fill") +
  facet_grid(. ~ SubDivision, scales="free", space="free") +
  theme(axis.text.x=element_text(size=14,color="black",angle=45,hjust=1.0), axis.text.y=element_text(size=16,color="black"), 
        axis.title=element_text(size=16,face="bold"),text=element_text(size=16),legend.title=element_blank())  
rel2

#################################
########### FIGURE 2C ###########
#################################
ord2 <- plot_ordination(so_f_r_late, ordinate(so_f_r_late, "MDS",distance="bray"),axes=1:2, color="SubDivision", shape="Treatment") +
  geom_point(size = 4)+
  scale_color_manual(values=c("#a04344","#4d5b6a","#84b59c","#686a47","#ecc363"))+
  theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
        axis.title=element_text(size=14),text=element_text(size=14))
ord2

#################################
########### PLOT ALL ############
#################################
#requires library("patchwork")
layout2 <- "
AABB
CCBB"
Figure2 <- alpha2 + rel2 + ord2 + plot_layout(design = layout2)
#Plot this directly in console to be able to view the complete figure in the plots viewer
Figure2
```

```{r, Figure 3, include=TRUE}
#Indicator analyses are below along with scripts to make supporting files for tree-building
#Also makes dataframes that were used to make annotation files for iTOL
#iTOL was used to plot trees and change aesthetics
all <- import_biom("feature-table.biom")
metadata <- import_qiime_sample_data("metadata.txt")
tree <- read_tree("tree.nwk")
all <- merge_phyloseq(all, metadata, tree)
all.root <- subset_samples(all, SampleType == "Root")
metadata <- read.delim("metadata.txt")
taxonomy <- read.delim("taxonomy.tsv")
repseqs <- read.fasta(file = "rep-seqs.fasta")

#this will merge all ASVs that have matching taxonomy assignments at the genus level
#it will retain, but not merge, ASVs whose assigned genera are either NA, "", " ", "\t", or "g__"
#takes a while...
all.genglom <- tax_glom(all.root, taxrank=rank_names(all)[6], 
                NArm = FALSE, 
                bad_empty = c(NA, "", " ", "\t", "g__"))

#remove ultralow abundance ASVs to prevent false postives in indicator analysis
#we're interested in ASVs that are reliable indicators
#an ASV with 10 (out of 10K!) reads in one or two samples in a category could be a significant indicator simply by virtue of not being found at all in other categories
#removes ASVs with fewer than 25 total reads and must be present in at least 1% of all samples (i.e. 4 samples in this case)
filter <- phyloseq::genefilter_sample(all.genglom, filterfun_sample(function(x) x >= 25), 
                                       A = 0.01*nsamples(all.genglom)-1)
all.genglom.filtered <- prune_taxa(filter, all.genglom)

#remove any samples with fewer than 2000 reads, then transform to relative abundance
all.genglom.filtered <- prune_samples(sample_sums(all.genglom.filtered)>=2000, all.genglom.filtered)
all.genglom.filtered <- transform_sample_counts(all.genglom.filtered, function(x) x/sum(x))
tax.table <- clean_taxtable(all.genglom.filtered)

spatialOrg <- subset_samples(all.genglom.filtered, Experiment == "Spatial Organization")
spatialOrg_roots <- subset_samples(spatialOrg, SampleType == "Root")

#run indicator analysis on each subset of data
SO.tips <- subset_samples(spatialOrg_roots, SubDivision == "SubRootBottom")
SO.mids <- subset_samples(spatialOrg_roots, SubDivision == "SubRootMiddle")
SO.bases <- subset_samples(spatialOrg_roots, SubDivision == "SubRootTop")
SO.tipIndicators <- get_indicators(SO.tips, groupfactor = "Treatment_all")
SO.midIndicators <- get_indicators(SO.mids, groupfactor = "Treatment_all")
SO.baseIndicators <- get_indicators(SO.bases, groupfactor = "Treatment_all")

#subset only  significant indicators
SO.tipIndicators$sig.indvals$subsection <- "tip"
SO.midIndicators$sig.indvals$subsection <- "mid"
SO.baseIndicators$sig.indvals$subsection <- "base"

#combine all subsections into the same df
SO.subIndicators <- rbind(SO.tipIndicators$sig.indvals, SO.midIndicators$sig.indvals, SO.baseIndicators$sig.indvals)
SO.subASVs <- unique(subset(SO.subIndicators, select = ASV))
#make_fasta(SO.subASVs, repseqs, "SO.subIndicators.fasta")

#make a csv to help build aestethic mappings for iTOL
SO.subIndicators.iTOL <- SO.subIndicators
SO.subIndicators.iTOL$group <- as.numeric(SO.subIndicators.iTOL$group)
SO.subIndicators.iTOL$group[SO.subIndicators.iTOL$group == 2] <- -1
SO.subIndicators.iTOL <- subset(SO.subIndicators.iTOL, select = c(ASV, group, subsection))
SO.subIndicators.iTOL <- recast(SO.subIndicators.iTOL, ASV~subsection, measure.var = "group")
SO.subIndicators.iTOL[is.na(SO.subIndicators.iTOL)] <- 0
SO.subIndicators.iTOL <- merge(SO.subIndicators.iTOL, tax.table, by = "ASV")
#write.csv(SO.subIndicators.iTOL, row.names = FALSE, quote = FALSE, file = "iTOL-files/SO-subIndicators-iTOL.csv")
```

```{r, Figure 4, include=TRUE}
#subset main data file to only include samples from the split pot experiments
sp <- subset_samples(millet_base, Experiment=="Split Pot")
#remove any metadata columns that do not pertain to any samples in the subsetted data
colnames(sp@sam_data)
sp@sam_data <- subset(sp@sam_data, select=-c(Latin, RootType, Divison, SubDivision, Plant, Root, Experiment, Notes))
#Removed ASVs without at least 1 read in at least 1 sample
sp_filter <- filter_taxa(sp, function(x) sum(x >= 1) >= 1, TRUE)
sp_filter
#Removed samples with no reads
sp_filter <- prune_samples(sample_sums(sp_filter) > 0, sp_filter)
sp_filter #removed 1 sample

#################################
########### FIGURE 4A ###########
#################################
sp_filterM <- subset_samples(sp_filter, Species=="Japanese Millet")
sp_filterMRt <- subset_samples(sp_filterM, SampleType=="Root")

# ord4a <- plot_ordination(sp_filterMRt, ordinate(sp_filterMRt, "MDS",distance="bray"),axes=1:2, color="SideTreatment") +
#   geom_point(size = 4)+
#   scale_color_manual(values=c("#a04344","#4d5b6a","#c57b67","#6a90b4"))+
#   theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"),
#         axis.title=element_text(size=14),text=element_text(size=14))
# ord4a

ord4a <- plot_ordination(sp_filterMRt, ordinate(sp_filterMRt, "CAP",distance="bray", ~SideTreatment),axes=1:2, color="SideTreatment") +
  geom_point(size = 4)+
  scale_color_manual(values=c("#a04344","#4d5b6a","#c57b67","#6a90b4"))+
  theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
        axis.title=element_text(size=14),text=element_text(size=14))
ord4a

#################################
########### FIGURE 4B ###########
#################################

sp_filterMRh <- subset_samples(sp_filterM, SampleType=="Rhizo")

# ord4b <- plot_ordination(sp_filterMRh, ordinate(sp_filterMRh, "MDS",distance="bray"),axes=c(1,3), color="SideTreatment") +
#   geom_point(size = 4)+
#   scale_color_manual(values=c("#a04344","#4d5b6a","#c57b67","#6a90b4"))+
#   theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
#         axis.title=element_text(size=14),text=element_text(size=14))
# ord4b

ord4b <- plot_ordination(sp_filterMRh, ordinate(sp_filterMRh, "CAP",distance="bray", ~SideTreatment),axes=c(1,2), color="SideTreatment") +
  geom_point(size = 4)+
  scale_color_manual(values=c("#a04344","#4d5b6a","#c57b67","#6a90b4"))+
  theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
        axis.title=element_text(size=14),text=element_text(size=14))
ord4b

#################################
########### FIGURE 4C ###########
#################################

#subset to remove soil samples
sp_filterMRR <- subset_samples(sp_filterM, SampleType!="Soil")

#Transform
sp_filterMRR_trans  = transform_sample_counts(sp_filterMRR, function(x) 100 * x / sum(x) )
#Agglomerate taxa
spMRR_glomClass <- tax_glom(sp_filterMRR_trans, taxrank = "Class")
spMRR_glomClass@tax_table[is.na(spMRR_glomClass@tax_table)] <- "Unknown"
#Need to check number of classes for next step
spMRR_glomClass #114 Classes

#Class level sort based on taxa abundance
## By Class
#label taxa that are NOT in the top 9 as "other"
spMRR_othersClass <- names(sort(taxa_sums(spMRR_glomClass), decreasing = TRUE)[10:104])
spMRR_merge <- merge_taxa(spMRR_glomClass, spMRR_othersClass,1)
spMRR_meltClass <- psmelt(spMRR_merge)
#convert Class to a character vector from a factor 
str(spMRR_meltClass$Class)
spMRR_meltClass$Class <- as.character(spMRR_meltClass$Class)
#replace all NAs with Other
spMRR_meltClass$Class[is.na(spMRR_meltClass$Class)] <- "Other"
# convert Class back to a factor
spMRR_meltClass$Class <- as.factor(spMRR_meltClass$Class)
levels(factor(spMRR_meltClass$Class))
#sort taxa
spMRR_meltClass$Class <- factor(spMRR_meltClass$Class, levels = c("Other", 
                                                                  "c__[Saprospirae]"  ,      
                                                                  "c__Acidobacteriia",
                                                                  "c__Alphaproteobacteria" , 
                                                                  "c__Bacilli",
                                                                  "c__Betaproteobacteria" ,
                                                                  "c__Solibacteres"    ,    
                                                                  "c__Sphingobacteriia",
                                                                  "c__Gammaproteobacteria",
                                                                  "c__Actinobacteria" ))
#sort SideTreatment
spMRR_meltClass$SideTreatment <- factor(spMRR_meltClass$SideTreatment, levels = c("FullWater","SplitWater","SplitDrought","FullDrought"))

rel4c <- ggplot(spMRR_meltClass, aes(x = SideTreatment, y = Abundance,  fill=Class, order=Class))+
  scale_fill_manual(values=farrowAndBall_palette)+
  geom_bar(stat = "identity", position = "fill") +
  facet_grid(~SampleType, scales="free", space="free", labeller=as_labeller(c("Rhizo"="Rhizosphere","Root"="Root"))) +
  theme(axis.text.x=element_text(size=14,color="black",angle=45,hjust=1.0), axis.text.y=element_text(size=16,color="black"), 
        axis.title=element_text(size=16,face="bold"),text=element_text(size=16),legend.title=element_blank())  
rel4c

#################################
########### PLOT ALL ############
#################################
#requires library("patchwork")
Figure4 <- ord4a + ord4b + rel4c
#Plot this directly in console to be able to view the complete figure in the plots viewer
Figure4
```

```{r, Figure 5, include=TRUE}
#subset main data file to only include samples from the split pot experiments
sr <- subset_samples(millet_base, Experiment=="Live/Dead Root")
#remove any metadata columns that do not pertain to any samples in the subsetted data
colnames(sr@sam_data)
sr@sam_data <- subset(sr@sam_data, select=-c(Species, Latin, SideTreatment, Divison, SubDivision, Timepoint, Plant, Root, Experiment, Notes))
#Removed ASVs without at least 1 read in at least 1 sample
sr_filter <- filter_taxa(sr, function(x) sum(x >= 1) >= 1, TRUE)
sr_filter
#Removed samples with no reads
sr_filter <- prune_samples(sample_sums(sr_filter) > 0, sr_filter)
sr_filter #removed 0 sample

#################################
########### FIGURE 5A ###########
#################################

sr_filterR <- subset_samples(sr_filter, SampleType=="Root")

# ord5 <- plot_ordination(sr_filterR, ordinate(sr_filterR, "MDS",distance="bray"),axes=c(1,2), color="RootType", shape="Treatment") +
#   geom_point(size = 4)+
#   scale_color_manual(values=c("#a04344","#4d5b6a"))+
#   #"#c57b67","#6a90b4"))+
#   theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
#         axis.title=element_text(size=14),text=element_text(size=14))
# ord5

ord5 <- plot_ordination(sr_filterR, ordinate(sr_filterR, "CAP",distance="bray", ~Treatment+RootType),axes=c(1,2), color="RootType", shape="Treatment") +
  geom_point(size = 4)+
  scale_color_manual(values=c("#a04344","#4d5b6a"))+
  #"#c57b67","#6a90b4"))+
  theme(axis.text.x=element_text(size=14,color="black",angle=90), axis.text.y=element_text(size=14,color="black"), 
        axis.title=element_text(size=14),text=element_text(size=14))
ord5

#################################
########### FIGURE 5B ###########
#################################
#Transform
sr_filterR_trans  = transform_sample_counts(sr_filterR, function(x) 100 * x / sum(x) )
#Agglomerate taxa
srR_glomClass <- tax_glom(sr_filterR_trans, taxrank = "Class")
srR_glomClass@tax_table[is.na(srR_glomClass@tax_table)] <- "Unknown"
#Need to check number of classes for next step
srR_glomClass #99 Classes

#Class level sort based on taxa abundance
#label taxa that are NOT in the top 9 as "other"
srR_othersClass <- names(sort(taxa_sums(srR_glomClass), decreasing = TRUE)[10:99])
srR_merge <- merge_taxa(srR_glomClass, srR_othersClass,1)
srR_meltClass <- psmelt(srR_merge)
#convert Class to a character vector from a factor 
str(srR_meltClass$Class)
srR_meltClass$Class <- as.character(srR_meltClass$Class)
#replace all NAs with Other
srR_meltClass$Class[is.na(srR_meltClass$Class)] <- "Other"
# convert Class back to a factor
srR_meltClass$Class <- as.factor(srR_meltClass$Class)
levels(factor(srR_meltClass$Class))
#sort taxa
srR_meltClass$Class <- factor(srR_meltClass$Class, levels = c("Other", 
                                                              "c__[Saprospirae]"  ,      
                                                              "c__Flavobacteriia",
                                                              "c__Alphaproteobacteria" , 
                                                              "c__Sphingobacteriia" ,
                                                              "c__Betaproteobacteria" ,
                                                              "c__Deltaproteobacteria",
                                                              "c__Cytophagia"  ,        
                                                              "c__Gammaproteobacteria",
                                                              "c__Actinobacteria" ))
#sort SideTreatment
srR_meltClass$Treatment <- factor(srR_meltClass$Treatment, levels = c("Watered","Drought"))
srR_meltClass$RootType <- factor(srR_meltClass$RootType, levels = c("Living","Dead"))

rel5 <- ggplot(srR_meltClass, aes(x = Treatment, y = Abundance,  fill=Class, order=Class))+
  scale_fill_manual(values=farrowAndBall_palette)+
  geom_bar(stat = "identity", position = "fill") +
  facet_grid(~RootType, scales="free", space="free") +
  theme(axis.text.x=element_text(size=14,color="black",angle=45,hjust=1.0), axis.text.y=element_text(size=16,color="black"), 
        axis.title=element_text(size=16,face="bold"),text=element_text(size=16),legend.title=element_blank())  
rel5

## PERMANOVA for root samples
set.seed(1)

# Subset phyloseq object, calculate bray curtis distance matrix, make data frame
root_dist <- phyloseq::distance(sr_filterR, method = "bray")
root_data <- data.frame(sample_data(sr_filterR))

# Then use an adonis test to look at your experimental variables.
adonis_root <- adonis(root_dist ~ Treatment + RootType, data = root_data)
adonis_root

# get test results as a table (necessary for next step)
adonis_root_tab <- adonis_root$aov.tab
adonis_root_tab_rows <- rownames(adonis_root_tab)

# Get a 'for-loop' to get your percent-variances attributable to each factor.
for(i in 1:4){ #1:4 is how many rows you have in table
  a1 <- c("The percent of variance attributable to the factor")
  b2 <- toString(adonis_root_tab_rows[i])
  c3 <- c("is")
  d4 <- toString(round(adonis_root_tab$SumsOfSqs[i]/adonis_root_tab$SumsOfSqs[4], 3)) #1st number should be "total" cell, 2nd number is how many digits it rounds to
  print(paste(a1,b2,c3,d4))
}

#################################
########### PLOT ALL ############
#################################
#requires library("patchwork")
Figure5 <- ord5 + rel5
#Plot this directly in console to be able to view the complete figure in the plots viewer
Figure5

```

```{r, Figure 6, include=TRUE}
#Indicator analyses are below along with scripts to make supporting files for tree-building
#Also makes dataframes that were used to make annotation files for iTOL
#iTOL was used to plot trees and change aesthetics
droughtLevels <- subset_samples(all.genglom.filtered, Experiment == "Drought Levels")
splitPot <- subset_samples(all.genglom.filtered, Experiment == "Split Pot")
splitPot_millet <- subset_samples(splitPot, Species == "Japanese Millet")
liveDead <- subset_samples(all.genglom.filtered, Experiment == "Live/Dead Root")
spatialOrg <- subset_samples(all.genglom.filtered, Experiment == "Spatial Organization")

#get indicators for each experiment using a function defined in block 1 above
DL.root.indicators <- get_indicators(droughtLevels, groupfactor = "Treatment_all")
SO.root.indicators <- get_indicators(spatialOrg, groupfactor = "Treatment_all")
LD.root.indicators <- get_indicators(liveDead, groupfactor = "Treatment_all")
SPmillet.root.indicators <- get_indicators(splitPot_millet, groupfactor = "Treatment_all")

#combine all significant root indicators into a single data frame
LD.root.indicators$sig.indvals$exp <- "LD"
DL.root.indicators$sig.indvals$exp <- "DL"
SPmillet.root.indicators$sig.indvals$exp <- "SPM"
SO.root.indicators$sig.indvals$exp <- "SO"
root.indicators <- rbind(LD.root.indicators$sig.indvals, DL.root.indicators$sig.indvals, 
                      SPmillet.root.indicators$sig.indvals, SO.root.indicators$sig.indvals)
root.indicators$type <- "root"

#combine all root indicator values, not just significant ones in a single data frame
LD.root.indicators$all.indvals$exp <- "LD"
DL.root.indicators$all.indvals$exp <- "DL"
SPmillet.root.indicators$all.indvals$exp <- "SPM"
SO.root.indicators$all.indvals$exp <- "SO"
all.root.indicators <- rbind(LD.root.indicators$all.indvals, DL.root.indicators$all.indvals, 
                      SPmillet.root.indicators$all.indvals, SO.root.indicators$all.indvals)
all.root.indicators$type <- "root"

all.drought <- subset(all.root.indicators, group == "Drought")
all.water <- subset(all.root.indicators, group == "Watered")

#creates dfs that are useful in the making of annotation files for iTOL
heatmap <- all.root.indicators
heatmap$indcls <- if_else(heatmap$group == "Watered", heatmap$indcls*-1, heatmap$indcls)
heatmap <- subset(heatmap, select = c(ASV, indcls, exp))
heatmap <- recast(heatmap, ASV~exp, measure.var = "indcls")
heatmap[is.na(heatmap)] <- "X"
heatmap <- merge(heatmap, tax.table, by = "ASV")
 
#Below are files that can be written out to your working directory to help make annotation files for iTOL
#write.csv(heatmap, row.names = FALSE, quote = FALSE, file = "figure6-heatmap.csv")
#write.csv(root.indicators, row.names = FALSE, quote = FALSE, file = "figure6-binary.csv")
#Below writes a fasta format file that can be used to build the tree seen in figure 6
#simply align via muscle and feed the alignment into FastTree; commands below:
#muscle3.8.31_i86darwin64 -in file.fasta -out alignment.afa
#FastTree -nt alignment.afa > tree.tre
#make_fasta(heatmap, repseqs, "figure6.fasta")
```