R packages

The following code uses several packages particularly tidyr and plyr for organizing data tables and manipulating strings of text.

library(stringr)
library(tidyr)
library(plyr)
library(vegan)

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.5-6

library(knitr)
library(ape)

Data Import, Subsetting, and Reformatting

Prior to assigning functional guild with Funguild, the OTU table must first be reformatted so that it includes the taxonomy assignments from UNITE. To begin I first import the rarefied shared file from mothur.

Data Import

Below ITS data from drill shavings of phloem tissues collected from 47 Juglans nigra trees are imported into R. This includes a rarefied OTU table from mothur and the taxonomic classifications for each OTU. Samples were rarefied to 3400 sequences in mothur prior to import, resulting in the loss of one sample from IN. Taxonomic assignments were made in mothur using the UNITE database.

#Load OTU table (shared file in mothur) and taxonomy files. The first column of the Mothur shared file is titled "label" and is unimportant so
its2.ds.shared.rare<-read.table("Mothur_output/barkits2otutable.rarefied_jnigra17_fa19_gw.shared", header=TRUE, sep="\t", row.names=2)

its2.ds.tax<-read.table("Mothur_output/barkits2taxonomy_jnigra17_fa19_gw.cons.taxonomy", header=TRUE, sep="\t")

Subsetting of OTU Table

To begin, I first subset the Mothur OTU table by removing mothur associated metadata, remove singleton OTUs (OTUs with only 1 sequence following rarefaction), and subset the taxonomy file.

#Here I subset the taxonomy file using the OTU label, transpose the OTU table, and merge the taxonomy info for import into funguild.  

#First I remove Mothur associated metatdata that includes the percent identity threshold and the number of OTUs. 
its2.ds.shared.rare<-(its2.ds.shared.rare[,3:2696])

#I then remove singleton OTUs or OTUs with only one sequence following rarefaction. 
its2.ds.shared.rare.nosingles<-as.data.frame(its2.ds.shared.rare[,(colSums(its2.ds.shared.rare))>1])

its2.ds.otu.lab<-labels(its2.ds.shared.rare.nosingles)

#Below I subset the taxonomy file. This is a universal taxonomy file for multiple compartments. We just want the OTUs specific to the caulosphere.
its2.ds.tax.sub<-subset(its2.ds.tax, its2.ds.tax$OTU %in% its2.ds.otu.lab[[2]])

Reformatting of OTU Table

Prior to running funguild, the OTU table needs to be transposed and the taxonomy strings from UNITE added to meet the input requirements for funguild.

#To add taxonomy information I first transposed the OTU table so that OTU IDs are the row names and sample names are the Column names. 
t.its2.ds.shared.rare.nosingles<-t(its2.ds.shared.rare.nosingles)

#I then create a column called otu that contains the OTU ID. This will be used to merge the OTU table to the taxonomy table.
t.its2.ds.shared.rare.nosingles<-data.frame(otu=labels(t.its2.ds.shared.rare.nosingles)[[1]],t.its2.ds.shared.rare.nosingles)

#I then merge the taxonomy file and the transposed OTU table using the OTU ID column. 
t.its2.ds.otutab.tax<-merge(t.its2.ds.shared.rare.nosingles,its2.ds.tax.sub, by.x="otu",by.y="OTU")

#Then I remove the column titled "size" from the new transposed OTU table with taxonomy information. This column does not really provide any useful information in this context and following funguild assignments the OTU table will be pretty heavily cluttered. 
t.its2.ds.otutab.tax<-t.its2.ds.otutab.tax[,-48]

#I then change the column names of the OTU ID column and the taxonomy column to align with the input requirements for funguild. 
names(t.its2.ds.otutab.tax)[1]<-"OTU_ID"
names(t.its2.ds.otutab.tax)[48]<-"taxonomy"

#write.table(t.its2.ds.otutab.tax, "~/Documents/microbiome2017//funguild.ds.input.file_FA19.txt",sep="\t", row.names=FALSE)

Funguild Assignments Import

Following functional guild assignment on funguild using the funguild webserver (http://www.stbates.org/guilds/app.php). Funguild assigned 632 OTUs out of 1578 OTUs to a functional guild. Below I import the functional guild assignments into R.

# After running funguild in linux, I then import the funguild output table. 
t.funguild.ds<-read.delim("FunGuild/barkfunguildoutput_jnigra17_fa19.25.11_ajo.guilds.txt",header=TRUE,sep="\t",fill=TRUE,strip.white=TRUE)

#When importing the soil funguild outputs, more rows than what was present in the input may be created. This is because some of the text from the funguild citation wraps inappropriately and ends up in its own row. These rows will be removed in the following chunk.

Funguild Table Reformatting

Prior to creating relative abundance plots, the functional guild output needs to be reformatted. First, compound guild strings need to be broken down into individual guilds and OTUs need to be weighted based on the number of guilds to which they were assigned. Then, OTUs are summed by guild and state.

# When importing the funguild table, some columns containing citation information may get incorrectly wrapped so that they have their own rows. These are first removed with na.omit. 
t.funguild.ds<-na.omit(t.funguild.ds)

#Here I am subsetting the data to pull out the the functional guild information from the funguild output. I then create a data frame that includes just the functional guild information and the abundance by site data. This table is currently a transposed OTU table. 
t.ds.otu.id<-t.funguild.ds[,1]
t.guild.ds<-t.funguild.ds[,52]
t.otutab.ds<-(t.funguild.ds[,2:47])
t.otutab.guild.ds<-data.frame(otu=t.ds.otu.id,Guilds=t.guild.ds,t.otutab.ds)
head(t.otutab.guild.ds)

##       otu                                              Guilds
## 1 Otu0003                                                   -
## 2 Otu0001 Animal Pathogen-Plant Pathogen-Undefined Saprotroph
## 3 Otu0002                                                   -
## 4 Otu0005                                                   -
## 5 Otu0004                                     Wood Saprotroph
## 6 Otu0011                                                   -
##   IN_MCB10_272s_S46_L001_R1_001 IN_MCB11_272s_S7_L001_R1_001
## 1                          1152                         1097
## 2                             4                            1
## 3                             0                            0
## 4                           349                          672
## 5                             0                            0
## 6                             0                            0
##   IN_MCB16_132s_S70_L001_R1_001 IN_MCB17_132s_S78_L001_R1_001
## 1                           223                           281
## 2                             3                             6
## 3                             1                             1
## 4                           222                           161
## 5                             0                             0
## 6                             0                             0
##   IN_MCB24_272s_S94_L001_R1_001 IN_MCB26_130s_S62_L001_R1_001
## 1                           156                          1719
## 2                             5                             0
## 3                             0                             0
## 4                            68                           873
## 5                             0                             0
## 6                             0                             0
##   IN_MCB27_132s_S86_L001_R1_001 IN_MCB28_WTs_S15_L001_R1_001
## 1                          1489                         1089
## 2                             3                            1
## 3                             0                            0
## 4                           332                          343
## 5                             0                            0
## 6                             0                            0
##   IN_MCB2_130s_S38_L001_R1_001 IN_MCB6_55s_S14_L001_R1_001
## 1                          705                        1061
## 2                            6                          22
## 3                            0                           0
## 4                          485                         358
## 5                            0                           0
## 6                            0                           0
##   IN_MCB7_55s_S22_L001_R1_001 IN_MCB9_130s_S54_L001_R1_001
## 1                         687                          969
## 2                           0                            2
## 3                           0                            0
## 4                         391                          409
## 5                           0                            0
## 6                           0                            0
##   IN_WT_MCB_29s_S23_L001_R1 IN_WT_MCB_33s_S31_L001_R1 TN_130_As_S71_L001_R1_001
## 1                       737                      1294                       187
## 2                        29                         1                         1
## 3                         0                         0                         0
## 4                       676                       140                       480
## 5                         0                         0                         0
## 6                         0                         0                      2022
##   TN_130_Bs_S79_L001_R1_001 TN_130_Cs_S87_L001_R1_001 TN_132_As_S95_L001_R1_001
## 1                       478                      1086                        73
## 2                        14                       146                        85
## 3                         0                         0                         0
## 4                       291                        78                       121
## 5                         0                         0                         0
## 6                         0                        88                        28
##   TN_132_Bs_S47_L001_R1_001 TN_132_Cs_S8_L001_R1_001 TN_272_As_S16_L001_R1_001
## 1                        55                       46                         8
## 2                         1                       11                        33
## 3                         0                        0                         0
## 4                         3                       99                         0
## 5                         0                        0                         0
## 6                       972                        4                       380
##   TN_272_Bs_S24_L001_R1_001 TN_272_Cs_S32_L001_R1_001 TN_55_As_S39_L001_R1_001
## 1                        14                         4                      184
## 2                        42                        45                        0
## 3                         0                         0                        0
## 4                         1                         0                        4
## 5                         0                         0                        0
## 6                         1                         0                        0
##   TN_55_Bs_S55_L001_R1_001 TN_55_Cs_S63_L001_R1_001 TN_LS_1s_S40_L001_R1_001
## 1                       75                      671                      271
## 2                        4                        2                        0
## 3                        0                        0                        0
## 4                       12                       40                      278
## 5                        0                        0                        0
## 6                        1                        1                        0
##   TN_LS_2s_S56_L001_R1_001 TN_LS_3s_S64_L001_R1_001 TN_WT_MB_19s_S72_L001_R1
## 1                       90                      285                     1184
## 2                        3                        0                        5
## 3                        0                        2                        0
## 4                      242                      160                      293
## 5                        0                        0                        0
## 6                        0                        2                        4
##   TN_WT_MB_20s_S80_L001_R1 TN_WT_MB_21s_S88_L001_R1
## 1                      633                      256
## 2                        0                        0
## 3                        0                        0
## 4                      288                      291
## 5                        0                        0
## 6                        0                        0
##   WA_BNL17_272s_S77_L001_R1_001 WA_BNL18_272s_S85_L001_R1_001
## 1                             0                             0
## 2                           737                           537
## 3                          1093                           638
## 4                             0                             0
## 5                            38                           282
## 6                             0                             0
##   WA_BNL19_55s_S13_L001_R1_001 WA_BNL20_130s_S29_L001_R1_001
## 1                            0                             0
## 2                         1001                          1740
## 3                          712                           309
## 4                            0                             0
## 5                           12                            74
## 6                            0                             0
##   WA_BNL21_WTs_S93_L001_R1_001 WA_BNL22_WTs_S45_L001_R1_001
## 1                            0                            0
## 2                          947                          854
## 3                          806                          872
## 4                            0                            0
## 5                           24                           44
## 6                            0                            0
##   WA_BNL23_WTs_S6_L001_R1_001 WA_RN10_272s_S69_L001_R1_001
## 1                           0                            0
## 2                        1189                          346
## 3                         215                         1397
## 4                           0                            0
## 5                          76                            0
## 6                           0                            0
##   WA_RN1_55s_S44_L001_R1_001 WA_RN2_55s_S5_L001_R1_001
## 1                          0                         0
## 2                        511                       592
## 3                        249                       783
## 4                          0                         0
## 5                       1033                         0
## 6                          0                         0
##   WA_RN4_130s_S21_L001_R1_001 WA_RN7_132s_S37_L001_R1_001
## 1                           0                           0
## 2                         224                         889
## 3                         232                         647
## 4                           0                           0
## 5                         516                          19
## 6                           0                           0
##   WA_RN8_132s_S53_L001_R1_001 WA_RN9_132s_S61_L001_R1_001
## 1                           0                           0
## 2                         321                         699
## 3                         378                         680
## 4                           0                           0
## 5                        1153                         294
## 6                           0                           0

#Next I use the stringr package to quantify the number of "-" in the guild names. To weight each value, I divide the total count by the number of dashes plus 1. I add 1 because a name with only 1 dash will have two functional guilds. 
library(stringr)
t.otutab.guild.div.ds<-data.frame(Guilds=t.guild.ds,(t.otutab.guild.ds[,3:48]/(str_count(t.otutab.guild.ds$Guilds,pattern="-")+1)))

#Then using separate_rows, I split the functional guild names by dash and create duplicate columns. 
library(tidyr)
t.otutab.guild.div.sep.ds<-separate_rows(t.otutab.guild.div.ds, Guilds, sep="-")

#Then using plyr I take the sum of the data by functional guild and then by site. 
library(plyr)
t.otutab.guild.div.sep.sum.ds<-as.data.frame(ddply(t.otutab.guild.div.sep.ds, .(Guilds),colwise(sum)))

#I then make the row name the functional guild and change the names so that the period between two parts of a name becomes a space. 
rownames(t.otutab.guild.div.sep.sum.ds)<-make.names(t.otutab.guild.div.sep.sum.ds$Guilds)


#I then remove the column that contains guild names since this is a duplicate of the row names. 
otutab.guild.div.sep.sum.ds<-t(t.otutab.guild.div.sep.sum.ds[,2:47])

site<-c("IN","IN","IN","IN","IN","IN","IN","IN","IN","IN","IN","IN","IN","IN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","TN","WA","WA","WA","WA","WA","WA","WA","WA","WA","WA","WA","WA","WA","WA")

#Then I create a data frame that contains the state information alongside the functional guild assignments and the sum functional guilds by state.  
otutab.guild.div.sep.sum.state.ds<-data.frame(state=site,otutab.guild.div.sep.sum.ds)
otutab.guild.div.sep.sum.statesum.ds<-(ddply(otutab.guild.div.sep.sum.state.ds, .(state),colwise(sum)))

#Then I rename the unassigned column, currently delimited with "X" to Unassigned and substitute . for spaces in the funguild names. 
colnames(otutab.guild.div.sep.sum.statesum.ds)[2]<-make.names("Unassigned")

colnames(otutab.guild.div.sep.sum.statesum.ds)<-sub("\\."," ",colnames(otutab.guild.div.sep.sum.statesum.ds))

Funguild Relative Abundance

Relative abundance is then calculated for each functional guild. The first group includes the unassigned and the second group excludes the unassigned.

# Then I calculate the relative abundance of each functional guild with unassigned included.
library(vegan)
funguild.relabund.other.ds<-decostand(otutab.guild.div.sep.sum.statesum.ds[,2:22],method="total")
funguild.relabund.other.ds<-data.frame(state=otutab.guild.div.sep.sum.statesum.ds$state,funguild.relabund.other.ds)
colnames(funguild.relabund.other.ds)<-sub("\\."," ",colnames(funguild.relabund.other.ds))

# Then I calculate the relative abundance of each functional guild. 
library(vegan)
funguild.relabund.ds<-decostand(otutab.guild.div.sep.sum.statesum.ds[,3:22],method="total")
funguild.relabund.ds<-data.frame(state=otutab.guild.div.sep.sum.statesum.ds$state,funguild.relabund.ds)
colnames(funguild.relabund.ds)<-sub("\\."," ",colnames(funguild.relabund.ds))

Relative Abundance Charts

To visualize differences in the relative abundance of different functional guilds between states, I reorganize the OTU table and taxonomy strings to create relative abundance stacked bar charts.

#In this chunk you input the relative abundance of narrow functional groups.
#Input guild data

#Convert to stacked bar plot format using the melt fucntion from the data.table package.
library(data.table)
funguild.melt.other.ds<-melt(funguild.relabund.other.ds, id.vars="state", variable.name="Guild")
ds.colors.n.funguild <- length(unique(funguild.melt.other.ds[,'Guild']))

#Plot in ggplot
library(ggplot2)
funguild.other.ds<-ggplot(funguild.melt.other.ds, aes(x=state, y=value, fill=Guild))+
  geom_bar(stat="identity", color="black" ,show.legend=TRUE)+
   scale_fill_manual(values=c('#999999','#003300','#006600','#666600','#99CC00','#33FF00','#CCFF33','#99FF99','#CCFFFF','#00FFFF','#66CCCC','#3399FF','#333FFF','#006699','#000099','#330066','#6600cc','#9900CC','#CC33FF','#CC99FF','#FFCCFF','#FF66FF','#FF3399','#CC0066','#993366','#993333','#CC6600','#FF6633','#006600','#99CC00','#33FF00','#CCFF33','#99FF99','#00FFFF','#3399FF','#333FFF','#000099','#330066','#6600cc','#CC99FF','#003300','#006600','#666600','#99CC00','#33FF00','#CCFF33','#99FF99','#CCFFFF','#666600' ))+
xlab("State") +
  ylab("Relative Abundance") +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text.x = element_text(size = 14), axis.ticks.x = element_blank(), axis.text.y = element_text(size = 14), axis.line.y.left = element_line(), axis.title.y = element_text(size = 14))
funguild.other.ds

library(data.table)
funguild.melt.ds<-melt(funguild.relabund.ds, id.vars="state", variable.name="Guild")
ds.colors.n.funguild <- length(unique(funguild.melt.ds[,'Guild']))

#Plot in ggplot
library(ggplot2)
funguild.ds<-ggplot(funguild.melt.ds, aes(x=state, y=value, fill=Guild))+
  geom_bar(stat="identity", color="black" ,show.legend=TRUE)+
   scale_fill_manual(values=c('#003300','#006600','#666600','#99CC00','#33FF00','#CCFF33','#99FF99','#CCFFFF','#00FFFF','#66CCCC','#3399FF','#333FFF','#006699','#000099','#330066','#6600cc','#9900CC','#CC33FF','#CC99FF','#FFCCFF','#FF66FF','#FF3399','#CC0066','#993366','#993333','#CC6600','#FF6633','#006600','#99CC00','#33FF00','#CCFF33','#99FF99','#00FFFF','#3399FF','#333FFF','#000099','#330066','#6600cc','#CC99FF','#003300','#006600','#666600','#99CC00','#33FF00','#CCFF33','#99FF99','#CCFFFF','#666600' ))+
  ylab("Relative Abundance") +
  theme(panel.border = element_blank(),panel.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text.x = element_text(size = 14), axis.ticks.x = element_blank(), axis.text.y = element_text(size = 14), axis.line.y.left = element_line(), axis.title.y = element_text(size = 14))+
  geom_text(data=NULL,aes(x=0.51, y=1.05,label="Caulosphere: Functional Guilds"),hjust=0,colour="black")+
guides(fill=guide_legend(ncol=2))
funguild.ds

#ggsave(filename = "ITS2_funguild_relativeabundance_ds_FA19_AJO.png", plot(funguild.ds),dpi=300)

Funguild Community Analyses

Mycoparasites

Subset for Mycoparasites and PCOA

I then reformat the OTU table for use in PCoA and funguild richness analyses for mycoparasites.

# First subset the funguild table to include only the fungal parasites
t.mycoparasite.ds<-subset(t.otutab.guild.div.sep.ds,t.otutab.guild.div.sep.ds$Guilds=="Fungal Parasite")

# Then change the row names and then transpose the table so that sites are the row names.
row.names(t.mycoparasite.ds)<-make.names(t.mycoparasite.ds$Guilds,unique = TRUE)
t.mycoparasite.ds<-t.mycoparasite.ds[,-1]
mycoparasite.ds<-as.data.frame((t(t.mycoparasite.ds)))

# Calculate relative abundance using the decostand function from vegan
mycoparasite.ds.relabund<-decostand(mycoparasite.ds,method="total")

# Create distance matrix using the bray method with the vegdist function from vegan
mycoparasite.ds.dist<-vegdist(mycoparasite.ds,method="bray")

#Perform PCOA with pcoa function from ape package
mycoparasite.ds.pcoa.results<-pcoa(mycoparasite.ds.dist)

Plotting of Mycoparasite Principal Coordinate Analysis

Then to plot the PCoA, I use the ggplot function. To use ggplot I need to first pull out the site scores for each sample and generate ellipses.

#Next I import the metadata for the particular habitat of interest. In this case I am working with drill shavings. I will use this metadata to subset the larger OTU file above. 
its2.ds.met<-read.table("Mothur_output/barkits2metadata_jnigra17_su19_ajo.txt",header=TRUE,sep="\t",fill=TRUE,strip.white=TRUE)

#Then we extract PCoA site scores (these will be the x and y coordinates for each sample)
mycoparasite.ds.pcoavec<-as.data.frame(mycoparasite.ds.pcoa.results$vectors)
mycoparasite.ds.pcoasitescores<-data.frame(PC1=mycoparasite.ds.pcoavec$Axis.1, PC2=mycoparasite.ds.pcoavec$Axis.2)

#Create a new dataframe that includes the site scores from above with metadata from your study. 
mycoparasite.ds.pcoagraph<-data.frame(mycoparasite.ds.pcoasitescores,PC1=mycoparasite.ds.pcoasitescores$PC1, PC2=mycoparasite.ds.pcoasitescores$PC2, State=its2.ds.met$State, Clone=its2.ds.met$Clone,group=its2.ds.met$Group)

#This is where you make confidence ellipses. I don't know what all the code means. Just know where to plug in my objects. 
mycoparasite.ds.pcoaellipse<-ordiellipse(mycoparasite.ds.pcoasitescores,mycoparasite.ds.pcoagraph$State, display="sites", kind="sd", draw="none")

df_ell.bark.mycop <- data.frame()
for(g in levels(mycoparasite.ds.pcoagraph$State)){
df_ell.bark.mycop <- rbind(df_ell.bark.mycop, cbind(as.data.frame(with(mycoparasite.ds.pcoagraph[mycoparasite.ds.pcoagraph$State==g,],                                                vegan:::veganCovEllipse(mycoparasite.ds.pcoaellipse[[g]]$cov,mycoparasite.ds.pcoaellipse[[g]]$center,mycoparasite.ds.pcoaellipse[[g]]$scale))) ,State=g))}

#Plot PCoA in ggplot
library(ggplot2)
mycoparasite.ds.pcoa<-ggplot(mycoparasite.ds.pcoagraph, aes(PC1,PC2, colour=State))+
  #geom_text(aes(label=group))+
   geom_path(data=df_ell.bark.mycop, aes(x=PC1, y=PC2, colour=State),size=0.5, linetype=1)+
  geom_point(aes(shape=Clone), size=3.5)+
  xlab("PC1 (27.7%)")+
  ylab("PC2 (14.5%)")+
 theme(axis.title.x=element_text(size=14, face="bold"))+
  theme(axis.title.y=element_text(size=14, face="bold"))+
  theme(axis.text.x=element_text(size=12, face="bold"))+
  theme(axis.text.y=element_text(size=12, face="bold"))+ 
  scale_color_manual(values=c("#006600","#3399FF","#FF9900"))+
  scale_shape_manual(values=c(16,10,8,7,0,2))+
  scale_x_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-0.60,0.60))+
  scale_y_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-.60,.60))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
    panel.border=element_rect(colour="black", size=1, fill=NA))+
  theme(panel.grid.major=element_blank(),
       panel.grid.minor=element_blank(),
       panel.background = element_blank(),
  panel.border=element_rect(color="black", size=1, fill=NA))+
  geom_text(data=NULL,aes(x=-0.60,y=0.60,label="Caulosphere: Mycoparasite"),hjust=0,colour="black")
mycoparasite.ds.pcoa

Mycoparasites PERMANOVA

To formally test how well state and clone explain differences in drill shaving fungal community composition, I perform a PERMANOVA using the adonis function from vegan.

#Below I perform permanovas by state and clone using the adonis function from vegan. 
library(vegan)
mycoparasite.ds.perm.state.inter<-adonis(mycoparasite.ds.relabund~its2.ds.met$State*its2.ds.met$Clone, method="bray",permutations=10000)
mycoparasite.ds.perm.state.inter

## 
## Call:
## adonis(formula = mycoparasite.ds.relabund ~ its2.ds.met$State *      its2.ds.met$Clone, permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                                     Df SumsOfSqs MeanSqs F.Model      R2
## its2.ds.met$State                    2    6.3570  3.1785 17.8715 0.43469
## its2.ds.met$Clone                    4    0.9413  0.2353  1.3232 0.06437
## its2.ds.met$State:its2.ds.met$Clone  8    1.8124  0.2266  1.2738 0.12393
## Residuals                           31    5.5135  0.1779         0.37701
## Total                               45   14.6243                 1.00000
##                                        Pr(>F)    
## its2.ds.met$State                   9.999e-05 ***
## its2.ds.met$Clone                      0.1554    
## its2.ds.met$State:its2.ds.met$Clone    0.1486    
## Residuals                                        
## Total                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mycoparasite.ds.perm.state.add<-adonis(mycoparasite.ds.relabund~its2.ds.met$State+its2.ds.met$Clone, method="bray",permutations=10000)
mycoparasite.ds.perm.state.add

## 
## Call:
## adonis(formula = mycoparasite.ds.relabund ~ its2.ds.met$State +      its2.ds.met$Clone, permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                   Df SumsOfSqs MeanSqs F.Model      R2    Pr(>F)    
## its2.ds.met$State  2    6.3570  3.1785 16.9211 0.43469 9.999e-05 ***
## its2.ds.met$Clone  4    0.9413  0.2353  1.2528 0.06437    0.2039    
## Residuals         39    7.3259  0.1878         0.50094              
## Total             45   14.6243                 1.00000              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mycoparasite.ds.perm.state.sing<-adonis(mycoparasite.ds.relabund~its2.ds.met$State, method="bray",permutations=10000)
mycoparasite.ds.perm.state.sing

## 
## Call:
## adonis(formula = mycoparasite.ds.relabund ~ its2.ds.met$State,      permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                   Df SumsOfSqs MeanSqs F.Model      R2    Pr(>F)    
## its2.ds.met$State  2    6.3570  3.1785  16.532 0.43469 9.999e-05 ***
## Residuals         43    8.2673  0.1923         0.56531              
## Total             45   14.6243                 1.00000              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mycoparasites Richness

Observed OTU richness is first calculated for each sample. Then due to the unbalanced nature of our study design, we use a type 3 ANOVA.

#Below I calculate observed species richness using the specnumber function from vegan. I then create a new data frame to include sample metadata. 
library(vegan)
mycoparasite.ds.rich<-specnumber(mycoparasite.ds)
mycoparasite.ds.rich.tabmet<-data.frame(State=its2.ds.met$State, clone=its2.ds.met$Clone, sample=its2.ds.met$Group,Richness=mycoparasite.ds.rich)

#I first test for a significant interaction between state and clone. Due to the unbalanced design we use a type 3 ANOVA. 
library(car)

## Loading required package: carData

mycoparasite.ds.anova.typ3.rich<-aov((Richness)~State*clone, data=mycoparasite.ds.rich.tabmet)
mycoparasite.ds.typ3aov.rich<-Anova(mycoparasite.ds.anova.typ3.rich,type="III")
mycoparasite.ds.typ3aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 320.33  1 43.7139 2.181e-07 ***
## State        61.21  2  4.1764   0.02477 *  
## clone        22.43  4  0.7652   0.55602    
## State:clone  46.70  8  0.7966   0.60988    
## Residuals   227.17 31                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(mycoparasite.ds.typ3aov.rich)

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	320.33333	1	43.7138665	0.0000002
State	61.20833	2	4.1763573	0.0247741
clone	22.42857	4	0.7651714	0.5560195
State:clone	46.69911	8	0.7965915	0.6098813
Residuals	227.16667	31	NA	NA

#The next line generates plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(mycoparasite.ds.anova.typ3.rich)

#There was no significant interaction between state and clone and thus we drop the interaction term.  
mycoparasite.ds.anova.add.rich<-aov((Richness)~State+clone, data=mycoparasite.ds.rich.tabmet)
mycoparasite.ds.add.aov.rich<-Anova(mycoparasite.ds.anova.add.rich,type="III")
mycoparasite.ds.add.aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 798.34  1 113.6879 4.033e-13 ***
## State       403.00  2  28.6944 2.174e-08 ***
## clone        16.44  4   0.5851    0.6753    
## Residuals   273.87 39                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Create a kable
mycoparasite.ds.add.aov.rich.tukey<-TukeyHSD(mycoparasite.ds.anova.add.rich)
mycoparasite.ds.add.aov.rich.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = (Richness) ~ State + clone, data = mycoparasite.ds.rich.tabmet)
## 
## $State
##            diff       lwr       upr     p adj
## TN-IN -0.484127 -2.784739  1.816485 0.8656333
## WA-IN -6.857143 -9.297311 -4.416975 0.0000001
## WA-TN -6.373016 -8.673628 -4.072404 0.0000001
## 
## $clone
##                diff       lwr      upr     p adj
## 132-130 -0.96263228 -4.644621 2.719356 0.9437163
## 272-130 -0.51818783 -4.200177 3.163801 0.9942482
## 55-130  -1.64285714 -5.431592 2.145878 0.7283919
## WT-130  -0.06448413 -3.523110 3.394142 0.9999981
## 272-132  0.44444444 -3.127609 4.016498 0.9964227
## 55-132  -0.68022487 -4.362214 3.001764 0.9838969
## WT-132   0.89814815 -2.443202 4.239498 0.9380731
## 55-272  -1.12466931 -4.806658 2.557319 0.9048954
## WT-272   0.45370370 -2.887647 3.795054 0.9949876
## WT-55    1.57837302 -1.880253 5.036999 0.6897683

mycoparasite.ds.add.aov.rich.tukey.table<-rbind(mycoparasite.ds.add.aov.rich.tukey$State,mycoparasite.ds.add.aov.rich.tukey$clone)
kable(mycoparasite.ds.add.aov.rich.tukey.table)

	diff	lwr	upr	p adj
TN-IN	-0.4841270	-2.784739	1.816485	0.8656333
WA-IN	-6.8571429	-9.297311	-4.416975	0.0000001
WA-TN	-6.3730159	-8.673628	-4.072404	0.0000001
132-130	-0.9626323	-4.644621	2.719356	0.9437163
272-130	-0.5181878	-4.200177	3.163801	0.9942482
55-130	-1.6428571	-5.431592	2.145878	0.7283919
WT-130	-0.0644841	-3.523110	3.394142	0.9999981
272-132	0.4444444	-3.127609	4.016498	0.9964227
55-132	-0.6802249	-4.362214	3.001764	0.9838969
WT-132	0.8981481	-2.443202	4.239498	0.9380731
55-272	-1.1246693	-4.806658	2.557319	0.9048954
WT-272	0.4537037	-2.887647	3.795054	0.9949876
WT-55	1.5783730	-1.880253	5.036999	0.6897683

#The below two lines generate plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(mycoparasite.ds.anova.add.rich)

#There was no significant interaction between state and clone and thus we drop the interaction term.  
mycoparasite.ds.anova.sing.rich<-aov((Richness)~State, data=mycoparasite.ds.rich.tabmet)
mycoparasite.ds.sing.aov.rich<-Anova(mycoparasite.ds.anova.sing.rich,type="III")
mycoparasite.ds.sing.aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##              Sum Sq Df F value  Pr(>F)    
## (Intercept) 1672.07  1 247.670 < 2e-16 ***
## State        424.13  2  31.412 3.9e-09 ***
## Residuals    290.30 43                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(mycoparasite.ds.sing.aov.rich)

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	1672.0714	1	247.67027	0
State	424.1332	2	31.41169	0
Residuals	290.3016	43	NA	NA

#Perform Tukey test
mycoparasite.ds.anova.sing.rich.tukey<-TukeyHSD(mycoparasite.ds.anova.sing.rich)
mycoparasite.ds.anova.sing.rich.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = (Richness) ~ State, data = mycoparasite.ds.rich.tabmet)
## 
## $State
##            diff       lwr       upr     p adj
## TN-IN -0.484127 -2.731699  1.763445 0.8606096
## WA-IN -6.857143 -9.241053 -4.473232 0.0000000
## WA-TN -6.373016 -8.620588 -4.125444 0.0000001

mycoparasite.ds.anova.sing.rich.tukey.table<-rbind(mycoparasite.ds.anova.sing.rich.tukey$State)
kable(mycoparasite.ds.anova.sing.rich.tukey.table)

	diff	lwr	upr	p adj
TN-IN	-0.484127	-2.731699	1.763445	0.8606096
WA-IN	-6.857143	-9.241053	-4.473232	0.0000000
WA-TN	-6.373016	-8.620588	-4.125444	0.0000001

#The below two lines generate plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(mycoparasite.ds.anova.sing.rich)

#I then create a box plot for the richness data by state using ggplot. 
library(ggplot2)
library(ggpubr)

## Loading required package: magrittr

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:tidyr':
## 
##     extract

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:ape':
## 
##     rotate

## The following object is masked from 'package:plyr':
## 
##     mutate

#Below I create a boxplot for richness values.
mycoparasite.ds.rich.plot<-ggboxplot(mycoparasite.ds.rich.tabmet, x="State", y="Richness", outlier.shape=NA)+
  geom_jitter(data=mycoparasite.ds.rich.tabmet,position=position_jitter(0.2), aes(shape=clone),size=3)+
  scale_shape_manual(values=c(16,10,8,7,0,2))+
  #geom_text(data=richness.summarized, aes(x=State, y = 25 + Richness, label=Group, fontface="bold"))+
  geom_text(data=NULL,aes(x=1, y=17.5, label="A"))+
  geom_text(data=NULL,aes(x=2, y=17.5, label="A"))+
  geom_text(data=NULL,aes(x=3, y=7.5, label="B"))+
  geom_text(data=NULL,aes(x=0.5,y=20,label="Caulosphere: Mycoparasite"),hjust=0)+
  theme(axis.line.x=element_blank(), 
        axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.x=element_blank())
mycoparasite.ds.rich.plot

Plant Pathogens

Subset for Plant Pathogens and PCOA

I then reformat the OTU table for use in PCoA and funguild richness analyses.

# First subset the funguild table to include only the fungal parasites
t.plant.pathogen.ds<-subset(t.otutab.guild.div.sep.ds,t.otutab.guild.div.sep.ds$Guilds=="Plant Pathogen")

# Then change the row names and then transpose the table so that sites are the row names.
row.names(t.plant.pathogen.ds)<-make.names(t.plant.pathogen.ds$Guilds,unique = TRUE)
t.plant.pathogen.ds<-t.plant.pathogen.ds[,-1]
plant.pathogen.ds<-as.data.frame((t(t.plant.pathogen.ds)))

# Calculate relative abundance using the decostand function from vegan
plant.pathogen.ds.relabund<-decostand(plant.pathogen.ds,method="total")

# Create distance matrix using the bray method with the vegdist function from vegan
plant.pathogen.ds.dist<-vegdist(plant.pathogen.ds,method="bray")

#Perform PCOA with pcoa function from ape package
plant.pathogen.ds.pcoa.results<-pcoa(plant.pathogen.ds.dist)

Plotting of plant.pathogen Principal Coordinate Analysis

Then to plot the PCoA, I use the ggplot function. To use ggplot I need to first pull out the site scores for each sample and generate ellipses.

#Next I import the metadata for the particular habitat of interest. In this case I am working with drill shavings. I will use this metadata to subset the larger OTU file above. 
its2.ds.met<-read.table("Mothur_output/barkits2metadata_jnigra17_su19_ajo.txt",header=TRUE,sep="\t",fill=TRUE,strip.white=TRUE)

#Then we extract PCoA site scores (these will be the x and y coordinates for each sample)
plant.pathogen.ds.pcoavec<-as.data.frame(plant.pathogen.ds.pcoa.results$vectors)
plant.pathogen.ds.pcoasitescores<-data.frame(PC1=plant.pathogen.ds.pcoavec$Axis.1, PC2=plant.pathogen.ds.pcoavec$Axis.2)

#Create a new dataframe that includes the site scores from above with metadata from your study. 
plant.pathogen.ds.pcoagraph<-data.frame(plant.pathogen.ds.pcoasitescores,PC1=plant.pathogen.ds.pcoasitescores$PC1, PC2=plant.pathogen.ds.pcoasitescores$PC2, State=its2.ds.met$State, Clone=its2.ds.met$Clone,group=its2.ds.met$Group)

#This is where you make confidence ellipses. I don't know what all the code means. Just know where to plug in my objects. 
plant.pathogen.ds.pcoaellipse<-ordiellipse(plant.pathogen.ds.pcoasitescores,plant.pathogen.ds.pcoagraph$State, display="sites", kind="sd", draw="none")

df_ell.bark.pp <- data.frame()
for(g in levels(plant.pathogen.ds.pcoagraph$State)){
df_ell.bark.pp <- rbind(df_ell.bark.pp, cbind(as.data.frame(with(plant.pathogen.ds.pcoagraph[plant.pathogen.ds.pcoagraph$State==g,],                                                vegan:::veganCovEllipse(plant.pathogen.ds.pcoaellipse[[g]]$cov,plant.pathogen.ds.pcoaellipse[[g]]$center,plant.pathogen.ds.pcoaellipse[[g]]$scale))) ,State=g))}

#Plot PCoA in ggplot
library(ggplot2)
plant.pathogen.ds.pcoa<-ggplot(plant.pathogen.ds.pcoagraph, aes(PC1,PC2, colour=State))+
  #geom_text(aes(label=group))+
  geom_path(data=df_ell.bark.pp, aes(x=PC1, y=PC2, colour=State),size=0.5, linetype=1)+
  geom_point(aes(shape=Clone), size=3.5)+
  xlab("PC1 (30.4%)")+
  ylab("PC2 (12.0%)")+
 theme(axis.title.x=element_text(size=14, face="bold"))+
  theme(axis.title.y=element_text(size=14, face="bold"))+
  theme(axis.text.x=element_text(size=12, face="bold"))+
  theme(axis.text.y=element_text(size=12, face="bold"))+ 
  scale_color_manual(values=c("#006600","#3399FF","#FF9900"))+
  scale_shape_manual(values=c(16,10,8,7,0,2))+
  scale_x_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-0.60,0.60))+
  scale_y_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-.60,.60))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
    panel.border=element_rect(colour="black", size=1, fill=NA))+
  theme(panel.grid.major=element_blank(),
       panel.grid.minor=element_blank(),
       panel.background = element_blank(),
  panel.border=element_rect(color="black", size=1, fill=NA))+
  geom_text(data=NULL,aes(x=-0.60,y=0.60,label="Caulosphere: Plant Pathogen"),hjust=0,colour="black")
plant.pathogen.ds.pcoa

#ggsave(filename = "ITS2_drillshavings_pcoa.png", plot(its2.ds.pcoa),dpi=300)

Plant pathogens PERMANOVA

To formally test how well state and clone explain differences in drill shaving fungal community composition, I perform a PERMANOVA using the adonis function from vegan.

#Below I perform permanovas by state and clone using the adonis function from vegan. 
library(vegan)
plant.pathogen.ds.perm.state.inter<-adonis(plant.pathogen.ds.relabund~its2.ds.met$State*its2.ds.met$Clone, method="bray",permutations=10000)
plant.pathogen.ds.perm.state.inter

## 
## Call:
## adonis(formula = plant.pathogen.ds.relabund ~ its2.ds.met$State *      its2.ds.met$Clone, permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                                     Df SumsOfSqs MeanSqs F.Model      R2
## its2.ds.met$State                    2    6.5120  3.2560 15.9590 0.39827
## its2.ds.met$Clone                    4    1.4093  0.3523  1.7269 0.08619
## its2.ds.met$State:its2.ds.met$Clone  8    2.1047  0.2631  1.2895 0.12872
## Residuals                           31    6.3247  0.2040         0.38682
## Total                               45   16.3506                 1.00000
##                                        Pr(>F)    
## its2.ds.met$State                   9.999e-05 ***
## its2.ds.met$Clone                      0.0201 *  
## its2.ds.met$State:its2.ds.met$Clone    0.1012    
## Residuals                                        
## Total                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plant.pathogen.ds.perm.state.add<-adonis(plant.pathogen.ds.relabund~its2.ds.met$State+its2.ds.met$Clone, method="bray",permutations=10000)
plant.pathogen.ds.perm.state.add

## 
## Call:
## adonis(formula = plant.pathogen.ds.relabund ~ its2.ds.met$State +      its2.ds.met$Clone, permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                   Df SumsOfSqs MeanSqs F.Model      R2    Pr(>F)    
## its2.ds.met$State  2    6.5120  3.2560 15.0644 0.39827 9.999e-05 ***
## its2.ds.met$Clone  4    1.4093  0.3523  1.6301 0.08619    0.0319 *  
## Residuals         39    8.4294  0.2161         0.51554              
## Total             45   16.3506                 1.00000              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plant pathogens Richness

Observed OTU richness is first calculated for each sample. Then due to the unbalanced nature of our study design, we use a type 3 ANOVA.

#Below I calculate observed species richness using the specnumber function from vegan. I then create a new data frame to include sample metadata. 
library(vegan)
plant.pathogen.ds.rich<-specnumber(plant.pathogen.ds)
plant.pathogen.ds.rich.tabmet<-data.frame(State=its2.ds.met$State, clone=its2.ds.met$Clone, sample=its2.ds.met$Group,Richness=plant.pathogen.ds.rich)

#I first test for a significant interaction between state and clone. Due to the unbalanced design we use a type 3 ANOVA. 
library(car)
plant.pathogen.ds.anova.typ3.rich<-aov((Richness)~State*clone, data=plant.pathogen.ds.rich.tabmet)
plant.pathogen.ds.typ3aov.rich<-Anova(plant.pathogen.ds.anova.typ3.rich,type="III")
plant.pathogen.ds.typ3aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 1452.00  1 64.7965 4.333e-09 ***
## State         51.00  2  1.1380   0.33350    
## clone        219.21  4  2.4456   0.06722 .  
## State:clone  306.17  8  1.7079   0.13604    
## Residuals    694.67 31                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(plant.pathogen.ds.typ3aov.rich)

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	1452.0000	1	64.796545	0.0000000
State	51.0000	2	1.137956	0.3334969
clone	219.2143	4	2.445649	0.0672210
State:clone	306.1719	8	1.707893	0.1360412
Residuals	694.6667	31	NA	NA

#The next line generates plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(plant.pathogen.ds.anova.typ3.rich)

#There was no significant interaction between state and clone and thus we drop the interaction term.  
plant.pathogen.ds.anova.add.rich<-aov((Richness)~State+clone, data=plant.pathogen.ds.rich.tabmet)
plant.pathogen.ds.add.aov.rich<-Anova(plant.pathogen.ds.anova.add.rich,type="III")
plant.pathogen.ds.add.aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 3368.2  1 131.2500 4.749e-14 ***
## State        274.9  2   5.3557  0.008808 ** 
## clone         52.1  4   0.5075  0.730491    
## Residuals   1000.8 39                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Create a kable
plant.pathogen.ds.add.aov.rich.tukey<-TukeyHSD(plant.pathogen.ds.anova.add.rich)
plant.pathogen.ds.add.aov.rich.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = (Richness) ~ State + clone, data = plant.pathogen.ds.rich.tabmet)
## 
## $State
##            diff        lwr       upr     p adj
## TN-IN  2.690476  -1.707537  7.088490 0.3064943
## WA-IN -3.500000  -8.164798  1.164798 0.1738793
## WA-TN -6.190476 -10.588490 -1.792463 0.0040325
## 
## $clone
##               diff       lwr      upr     p adj
## 132-130  0.6537698 -6.384981 7.692521 0.9988569
## 272-130 -0.5684524 -7.607203 6.470298 0.9993409
## 55-130  -1.8125000 -9.055315 5.430315 0.9516684
## WT-130   1.2470238 -5.364731 7.858779 0.9826098
## 272-132 -1.2222222 -8.050813 5.606369 0.9856906
## 55-132  -2.4662698 -9.505021 4.572481 0.8528807
## WT-132   0.5932540 -5.794308 6.980816 0.9988572
## 55-272  -1.2440476 -8.282798 5.794703 0.9863501
## WT-272   1.8154762 -4.572086 8.203038 0.9251428
## WT-55    3.0595238 -3.552231 9.671279 0.6786989

plant.pathogen.ds.add.aov.rich.tukey.table<-rbind(plant.pathogen.ds.add.aov.rich.tukey$State,plant.pathogen.ds.add.aov.rich.tukey$clone)
kable(plant.pathogen.ds.add.aov.rich.tukey.table)

	diff	lwr	upr	p adj
TN-IN	2.6904762	-1.707537	7.088490	0.3064943
WA-IN	-3.5000000	-8.164798	1.164798	0.1738793
WA-TN	-6.1904762	-10.588490	-1.792463	0.0040325
132-130	0.6537698	-6.384981	7.692521	0.9988569
272-130	-0.5684524	-7.607203	6.470298	0.9993409
55-130	-1.8125000	-9.055315	5.430315	0.9516684
WT-130	1.2470238	-5.364731	7.858779	0.9826098
272-132	-1.2222222	-8.050814	5.606369	0.9856906
55-132	-2.4662698	-9.505021	4.572481	0.8528807
WT-132	0.5932540	-5.794308	6.980816	0.9988572
55-272	-1.2440476	-8.282798	5.794703	0.9863501
WT-272	1.8154762	-4.572086	8.203038	0.9251428
WT-55	3.0595238	-3.552231	9.671279	0.6786989

#The below two lines generate plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(plant.pathogen.ds.anova.add.rich)

#There was no significant interaction between state and clone and thus we drop the interaction term.  
plant.pathogen.ds.anova.sing.rich<-aov((Richness)~State, data=plant.pathogen.ds.rich.tabmet)
plant.pathogen.ds.sing.aov.rich<-Anova(plant.pathogen.ds.anova.sing.rich,type="III")
plant.pathogen.ds.sing.aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 7825.8  1 319.5932 < 2.2e-16 ***
## State        301.8  2   6.1623  0.004434 ** 
## Residuals   1052.9 43                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(plant.pathogen.ds.sing.aov.rich)

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	7825.7857	1	319.593175	0.0000000
State	301.7888	2	6.162298	0.0044343
Residuals	1052.9286	43	NA	NA

#Perform Tukey test
plant.pathogen.ds.anova.sing.rich.tukey<-TukeyHSD(plant.pathogen.ds.anova.sing.rich)
plant.pathogen.ds.anova.sing.rich.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = (Richness) ~ State, data = plant.pathogen.ds.rich.tabmet)
## 
## $State
##            diff        lwr       upr     p adj
## TN-IN  2.690476  -1.589964  6.970916 0.2891179
## WA-IN -3.500000  -8.040092  1.040092 0.1593213
## WA-TN -6.190476 -10.470916 -1.910036 0.0029946

plant.pathogen.ds.anova.sing.rich.tukey.table<-rbind(plant.pathogen.ds.anova.sing.rich.tukey$State)
kable(plant.pathogen.ds.anova.sing.rich.tukey.table)

	diff	lwr	upr	p adj
TN-IN	2.690476	-1.589964	6.970916	0.2891179
WA-IN	-3.500000	-8.040092	1.040092	0.1593213
WA-TN	-6.190476	-10.470916	-1.910036	0.0029946

#The below two lines generate plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(plant.pathogen.ds.anova.sing.rich)

#I then create a box plot for the richness data by state using ggplot. 

library(ggplot2)
library(ggpubr)
#Below I create a boxplot for richness values.
plant.pathogen.ds.rich.plot<-ggboxplot(plant.pathogen.ds.rich.tabmet, x="State", y="Richness", outlier.shape=NA)+
  geom_jitter(data=plant.pathogen.ds.rich.tabmet,position=position_jitter(0.2), aes(shape=clone),size=3)+
  scale_shape_manual(values=c(16,10,8,7,0,2))+
  #geom_text(data=richness.summarized, aes(x=State, y = 25 + Richness, label=Group, fontface="bold"))+
  geom_text(data=NULL,aes(x=1, y=37.5, label="AB"))+
  geom_text(data=NULL,aes(x=2, y=37.5, label="A"))+
  geom_text(data=NULL,aes(x=3, y=29.5, label="B"))+
  geom_text(data=NULL,aes(x=0.5,y=40,label="Caulosphere: Plant Pathogen"),hjust=0)+
  theme(axis.line.x=element_blank(), 
        axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.x=element_blank())
plant.pathogen.ds.rich.plot

Saprotrophs

Subset for Plant Pathogens and PCOA

I then reformat the OTU table for use in PCoA and funguild richness analyses.

# First subset the funguild table to include only the fungal parasites
t.saprotroph.ds<-subset(t.otutab.guild.div.sep.ds,t.otutab.guild.div.sep.ds$Guilds=="Wood Saprotroph")

# Then change the row names and then transpose the table so that sites are the row names.
row.names(t.saprotroph.ds)<-make.names(t.saprotroph.ds$Guilds,unique = TRUE)
t.saprotroph.ds<-t.saprotroph.ds[,-1]
saprotroph.ds<-as.data.frame((t(t.saprotroph.ds)))

# Calculate relative abundance using the decostand function from vegan
saprotroph.ds.relabund<-decostand(saprotroph.ds,method="total")

# Create distance matrix using the bray method with the vegdist function from vegan
saprotroph.ds.dist<-vegdist(saprotroph.ds,method="bray")

#Perform PCOA with pcoa function from ape package
saprotroph.ds.pcoa.results<-pcoa(saprotroph.ds.dist)

Plotting of saprotroph Principal Coordinate Analysis

Then to plot the PCoA, I use the ggplot function. To use ggplot I need to first pull out the site scores for each sample and generate ellipses.

#Next I import the metadata for the particular habitat of interest. In this case I am working with drill shavings. I will use this metadata to subset the larger OTU file above. 
its2.ds.met<-read.table("Mothur_output/barkits2metadata_jnigra17_su19_ajo.txt",header=TRUE,sep="\t",fill=TRUE,strip.white=TRUE)

#Then we extract PCoA site scores (these will be the x and y coordinates for each sample)
saprotroph.ds.pcoavec<-as.data.frame(saprotroph.ds.pcoa.results$vectors)
saprotroph.ds.pcoasitescores<-data.frame(PC1=saprotroph.ds.pcoavec$Axis.1, PC2=saprotroph.ds.pcoavec$Axis.2)

#Create a new dataframe that includes the site scores from above with metadata from your study. 
saprotroph.ds.pcoagraph<-data.frame(saprotroph.ds.pcoasitescores,PC1=saprotroph.ds.pcoasitescores$PC1, PC2=saprotroph.ds.pcoasitescores$PC2, State=its2.ds.met$State, Clone=its2.ds.met$Clone,group=its2.ds.met$Group)

#This is where you make confidence ellipses. I don't know what all the code means. Just know where to plug in my objects. 
saprotroph.ds.pcoaellipse<-ordiellipse(saprotroph.ds.pcoasitescores,saprotroph.ds.pcoagraph$State, display="sites", kind="sd", draw="none")

df_ell.bark.sap <- data.frame()
for(g in levels(saprotroph.ds.pcoagraph$State)){
df_ell.bark.sap <- rbind(df_ell.bark.sap, cbind(as.data.frame(with(saprotroph.ds.pcoagraph[saprotroph.ds.pcoagraph$State==g,],                                                vegan:::veganCovEllipse(saprotroph.ds.pcoaellipse[[g]]$cov,saprotroph.ds.pcoaellipse[[g]]$center,saprotroph.ds.pcoaellipse[[g]]$scale))) ,State=g))}

#Plot PCoA in ggplot
library(ggplot2)
saprotroph.ds.pcoa<-ggplot(saprotroph.ds.pcoagraph, aes(PC1,PC2, colour=State))+
  #geom_text(aes(label=group))+
    geom_path(data=df_ell.bark.sap, aes(x=PC1, y=PC2, colour=State),size=0.5, linetype=1)+
  geom_point(aes(shape=Clone), size=3.5)+
xlab("PC1 (19.1%)")+
  ylab("PC2 (10.1%)")+
 theme(axis.title.x=element_text(size=14, face="bold"))+
  theme(axis.title.y=element_text(size=14, face="bold"))+
  theme(axis.text.x=element_text(size=12, face="bold"))+
  theme(axis.text.y=element_text(size=12, face="bold"))+ 
  scale_color_manual(values=c("#006600","#3399FF","#FF9900"))+
  scale_shape_manual(values=c(16,10,8,7,0,2))+
  scale_x_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-0.60,0.60))+
  scale_y_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-.60,.60))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
    panel.border=element_rect(colour="black", size=1, fill=NA))+
  theme(panel.grid.major=element_blank(),
       panel.grid.minor=element_blank(),
       panel.background = element_blank(),
  panel.border=element_rect(color="black", size=1, fill=NA))+
  geom_text(data=NULL,aes(x=-0.60,y=0.60,label="Caulosphere: Wood Saprotroph"),hjust=0,colour="black")
saprotroph.ds.pcoa

#ggsave(filename = "ITS2_drillshavings_pcoa.png", plot(its2.ds.pcoa),dpi=300)

Saprotrophs PERMANOVA

To formally test how well state and clone explain differences in drill shaving fungal community composition, I perform a PERMANOVA using the adonis function from vegan.

#Below I perform permanovas by state and clone using the adonis function from vegan. 
library(vegan)
saprotroph.ds.perm.state.inter<-adonis(saprotroph.ds.relabund~its2.ds.met$State*its2.ds.met$Clone, method="bray",permutations=10000)
saprotroph.ds.perm.state.inter

## 
## Call:
## adonis(formula = saprotroph.ds.relabund ~ its2.ds.met$State *      its2.ds.met$Clone, permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                                     Df SumsOfSqs MeanSqs F.Model      R2
## its2.ds.met$State                    2    6.3305  3.1653 12.8313 0.35522
## its2.ds.met$Clone                    4    1.4083  0.3521  1.4272 0.07902
## its2.ds.met$State:its2.ds.met$Clone  8    2.4355  0.3044  1.2341 0.13666
## Residuals                           31    7.6472  0.2467         0.42910
## Total                               45   17.8215                 1.00000
##                                        Pr(>F)    
## its2.ds.met$State                   9.999e-05 ***
## its2.ds.met$Clone                     0.05879 .  
## its2.ds.met$State:its2.ds.met$Clone   0.11429    
## Residuals                                        
## Total                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

saprotroph.ds.perm.state.add<-adonis(saprotroph.ds.relabund~its2.ds.met$State+its2.ds.met$Clone, method="bray",permutations=10000)
saprotroph.ds.perm.state.add

## 
## Call:
## adonis(formula = saprotroph.ds.relabund ~ its2.ds.met$State +      its2.ds.met$Clone, permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                   Df SumsOfSqs MeanSqs F.Model      R2    Pr(>F)    
## its2.ds.met$State  2    6.3305  3.1653 12.2433 0.35522 9.999e-05 ***
## its2.ds.met$Clone  4    1.4083  0.3521  1.3618 0.07902   0.08119 .  
## Residuals         39   10.0827  0.2585         0.56576              
## Total             45   17.8215                 1.00000              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

saprotroph.ds.perm.state.sing<-adonis(saprotroph.ds.relabund~its2.ds.met$State, method="bray",permutations=10000)
saprotroph.ds.perm.state.sing

## 
## Call:
## adonis(formula = saprotroph.ds.relabund ~ its2.ds.met$State,      permutations = 10000, method = "bray") 
## 
## Permutation: free
## Number of permutations: 10000
## 
## Terms added sequentially (first to last)
## 
##                   Df SumsOfSqs MeanSqs F.Model      R2    Pr(>F)    
## its2.ds.met$State  2    6.3305  3.1653  11.845 0.35522 9.999e-05 ***
## Residuals         43   11.4909  0.2672         0.64478              
## Total             45   17.8215                 1.00000              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plant pathogens Richness

Observed OTU richness is first calculated for each sample. Then due to the unbalanced nature of our study design, we use a type 3 ANOVA.

#Below I calculate observed species richness using the specnumber function from vegan. I then create a new data frame to include sample metadata. 
library(vegan)
saprotroph.ds.rich<-specnumber(saprotroph.ds)
saprotroph.ds.rich.tabmet<-data.frame(State=its2.ds.met$State, clone=its2.ds.met$Clone, sample=its2.ds.met$Group,Richness=saprotroph.ds.rich)

#I first test for a significant interaction between state and clone. Due to the unbalanced design we use a type 3 ANOVA. 
library(car)
saprotroph.ds.anova.typ3.rich<-aov((Richness)~State*clone, data=saprotroph.ds.rich.tabmet)
saprotroph.ds.typ3aov.rich<-Anova(saprotroph.ds.anova.typ3.rich,type="III")
saprotroph.ds.typ3aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 1160.33  1 86.2253 1.824e-10 ***
## State        158.17  2  5.8767  0.006856 ** 
## clone         46.93  4  0.8718  0.491894    
## State:clone   62.01  8  0.5760  0.789424    
## Residuals    417.17 31                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(saprotroph.ds.typ3aov.rich)

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	1160.33333	1	86.2253296	0.0000000
State	158.16667	2	5.8767479	0.0068556
clone	46.92857	4	0.8718252	0.4918936
State:clone	62.01128	8	0.5760137	0.7894244
Residuals	417.16667	31	NA	NA

#The next line generates plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(saprotroph.ds.anova.typ3.rich)

#There was no significant interaction between state and clone and thus we drop the interaction term.  
saprotroph.ds.anova.add.rich<-aov((Richness)~State+clone, data=saprotroph.ds.rich.tabmet)
saprotroph.ds.add.aov.rich<-Anova(saprotroph.ds.anova.add.rich,type="III")
saprotroph.ds.add.aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 2116.82  1 172.287 6.860e-16 ***
## State        560.17  2  22.796 2.772e-07 ***
## clone         53.96  4   1.098    0.3711    
## Residuals    479.18 39                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Create a kable
saprotroph.ds.add.aov.rich.tukey<-TukeyHSD(saprotroph.ds.anova.add.rich)
saprotroph.ds.add.aov.rich.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = (Richness) ~ State + clone, data = saprotroph.ds.rich.tabmet)
## 
## $State
##            diff        lwr       upr     p adj
## TN-IN -8.095238 -11.138385 -5.052091 0.0000003
## WA-IN -6.785714 -10.013459 -3.557969 0.0000250
## WA-TN  1.309524  -1.733623  4.352671 0.5513530
## 
## $clone
##               diff       lwr      upr     p adj
## 132-130 -0.3273810 -5.197751 4.542989 0.9996810
## 272-130 -2.1051587 -6.975529 2.765212 0.7306850
## 55-130  -2.2767857 -7.288355 2.734784 0.6933092
## WT-130   0.2619048 -4.313012 4.836821 0.9998313
## 272-132 -1.7777778 -6.502731 2.947175 0.8176747
## 55-132  -1.9494048 -6.819775 2.920966 0.7819912
## WT-132   0.5892857 -3.830503 5.009075 0.9953277
## 55-272  -0.1716270 -5.041997 4.698743 0.9999756
## WT-272   2.3670635 -2.052726 6.786853 0.5489628
## WT-55    2.5386905 -2.036226 7.113607 0.5145061

saprotroph.ds.add.aov.rich.tukey.table<-rbind(saprotroph.ds.add.aov.rich.tukey$State,saprotroph.ds.add.aov.rich.tukey$clone)
kable(saprotroph.ds.add.aov.rich.tukey.table)

	diff	lwr	upr	p adj
TN-IN	-8.0952381	-11.138385	-5.052091	0.0000003
WA-IN	-6.7857143	-10.013459	-3.557969	0.0000250
WA-TN	1.3095238	-1.733623	4.352671	0.5513530
132-130	-0.3273810	-5.197751	4.542989	0.9996810
272-130	-2.1051587	-6.975529	2.765212	0.7306850
55-130	-2.2767857	-7.288355	2.734784	0.6933092
WT-130	0.2619048	-4.313012	4.836821	0.9998313
272-132	-1.7777778	-6.502731	2.947176	0.8176747
55-132	-1.9494048	-6.819775	2.920966	0.7819912
WT-132	0.5892857	-3.830503	5.009075	0.9953277
55-272	-0.1716270	-5.041997	4.698743	0.9999756
WT-272	2.3670635	-2.052726	6.786853	0.5489628
WT-55	2.5386905	-2.036226	7.113607	0.5145061

#The below two lines generate plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(saprotroph.ds.anova.add.rich)

#There was no significant interaction between state and clone and thus we drop the interaction term.  
saprotroph.ds.anova.sing.rich<-aov((Richness)~State, data=saprotroph.ds.rich.tabmet)
saprotroph.ds.sing.aov.rich<-Anova(saprotroph.ds.anova.sing.rich,type="III")
saprotroph.ds.sing.aov.rich

## Anova Table (Type III tests)
## 
## Response: (Richness)
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 4500.1  1 362.948 < 2.2e-16 ***
## State        564.6  2  22.768 1.805e-07 ***
## Residuals    533.1 43                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kable(saprotroph.ds.sing.aov.rich)

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	4500.0714	1	362.94788	0e+00
State	564.5963	2	22.76842	2e-07
Residuals	533.1429	43	NA	NA

#Perform Tukey test
saprotroph.ds.anova.sing.rich.tukey<-TukeyHSD(saprotroph.ds.anova.sing.rich)
saprotroph.ds.anova.sing.rich.tukey

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = (Richness) ~ State, data = saprotroph.ds.rich.tabmet)
## 
## $State
##            diff        lwr       upr     p adj
## TN-IN -8.095238 -11.141104 -5.049372 0.0000002
## WA-IN -6.785714 -10.016343 -3.555086 0.0000215
## WA-TN  1.309524  -1.736342  4.355390 0.5538189

saprotroph.ds.anova.sing.rich.tukey.table<-rbind(saprotroph.ds.anova.sing.rich.tukey$State)
kable(saprotroph.ds.anova.sing.rich.tukey.table)

	diff	lwr	upr	p adj
TN-IN	-8.095238	-11.141104	-5.049372	0.0000002
WA-IN	-6.785714	-10.016343	-3.555086	0.0000215
WA-TN	1.309524	-1.736342	4.355390	0.5538189

#The below two lines generate plots to evaluate regression assumptions. 
par(mfrow=c(2,2))
plot(saprotroph.ds.anova.sing.rich)

#I then create a box plot for the richness data by state using ggplot. 

library(ggplot2)
library(ggpubr)
#Below I create a boxplot for richness values.
saprotroph.ds.rich.plot<-ggboxplot(saprotroph.ds.rich.tabmet, x="State", y="Richness", outlier.shape=NA)+
  geom_jitter(data=saprotroph.ds.rich.tabmet,position=position_jitter(0.2), aes(shape=clone),size=3)+
  scale_shape_manual(values=c(16,10,8,7,0,2))+
  #geom_text(data=richness.summarized, aes(x=State, y = 25 + Richness, label=Group, fontface="bold"))+
  geom_text(data=NULL,aes(x=1, y=20, label="B"))+
  geom_text(data=NULL,aes(x=2, y=17.5, label="A"))+
  geom_text(data=NULL,aes(x=3, y=17.5, label="B"))+
  geom_text(data=NULL,aes(x=0.5,y=30,label="Caulosphere: Wood Saprotrophs"),hjust=0)+
  theme(axis.line.x=element_blank(),
        axis.title.x=element_blank())
saprotroph.ds.rich.plot

Funguild_ds_FA19_AJO

Funguild_ds_FA19_AJO

R packages

Data Import, Subsetting, and Reformatting

Data Import

Subsetting of OTU Table

Reformatting of OTU Table

Funguild Assignments Import

Funguild Table Reformatting

Funguild Relative Abundance

Relative Abundance Charts

Funguild Community Analyses

Mycoparasites

Subset for Mycoparasites and PCOA

Plotting of Mycoparasite Principal Coordinate Analysis

Mycoparasites PERMANOVA

Mycoparasites Richness

Plant Pathogens

Subset for Plant Pathogens and PCOA

Plotting of plant.pathogen Principal Coordinate Analysis

Plant pathogens PERMANOVA

Plant pathogens Richness

Saprotrophs

Subset for Plant Pathogens and PCOA

Plotting of saprotroph Principal Coordinate Analysis

Saprotrophs PERMANOVA

Plant pathogens Richness