R Pipeline

biom convert -i $PATH/otus/otu_table_mc2_w_tax_no_ChloroMitoUnID.biom -o $PATH/otus/otu_table_mc2_w_tax_no_ChloroMitoUnID-JSON.biom --to-json --head-key taxonomy

library('phyloseq')

library('ggplot2')

library('vegan')

#$PATH represents the complete file path for the input files below

BCsoil<-import_biom("$PATH/otus/otu_table_mc2_w_tax_no_ChloroMitoUnID-JSON.biom","$PATH/otus/rep_set.tre", parseFunction=parse_taxonomy_default)

BCsoil_Map<-read.table("$PATH/CampusSoil-mapping-forR.txt", header=TRUE)

BCsoil_SampleData<-sample_data(as.data.frame(BCsoil_Map, check.rows=TRUE, row.names=sort(sample_names(BCsoil)), stringAsFactors=TRUE))

#Creating the phyloseq object

BCsoil<-merge_phyloseq(BCsoil, BCsoil_SampleData); BCsoil

#Assigning names to taxonomy levels

rank_names(BCsoil)

colnames(tax_table(BCsoil)) = c("Kingdom","Phylum","Class","Order","Family","Genus","Species")

#Alpha diversity metrics can also included "ACE", "InvSimpson"

BCsoil.Alpha<-estimate_richness(BCsoil, measures=c("Observed", "Chao1", "Shannon", "Simpson") ); BCsoil.Alpha

plot1<-plot_richness(BCsoil,  measures=alpha_meas) 

plot2<-plot_richness(BCsoil, "Season", color="Season", measures=alpha_meas) 

plot2 + geom_boxplot(data=plot1$data, aes(x= Season, y=value, color=NULL), alpha=0.1) + geom_point(size=10)

#Create distance matrix

BCsoil_WeightedUniFrac<-UniFrac(BCsoil, weighted=TRUE)

#Plot PCoA

plot_ordination(BCsoil, ordinate(BCsoil, "PCoA"), color="Season") + geom_point(size=10)

#Calculate statistic if the season explains the difference between samples

metadata<-as(sample_data(BCsoil.Rare), "data.frame")

adonis(distance(BCsoil.Rare, method="bray")~Season, data=metadata)

#RAREFACTION - not a universal practice - but an attempt to normalize data when samples can be unevenly sequenced - could strengthen comparative analyses

BCsoil.Rare<-rarefy_even_depth(BCsoil, sample.size = min(20000))

#Community proportion bars

BCsoil.Rare_1<-prune_taxa(taxa_sums(BCsoil.Rare)>1,BCsoil.Rare) #Select for OTUs with over 1 seqs

#plot_bar(BCsoil.Rare_1) #takes a while to run

BCsoil.Rare_1000<-prune_taxa(taxa_sums(BCsoil.Rare)>1000,BCsoil.Rare) #Select for OTUs with over 1000 seqs

plot_bar(BCsoil.Rare_1000, fill='Phylum') 

plot_tree(BCsoil.Rare_1000) 

BCsoil.Rare.Proteobacteria <-subset_taxa(BCsoil.Rare, Phylum == "p__Proteobacteria" ) 

PhylumPlot5<-plot_bar(BCsoil.Rare.Proteobacteria, fill='Order')

PhylumPlot5 + facet_wrap(~Order) + theme(legend.position = "none") 

set.seed(42) #Provides continuity for running the script and getting the same output repeatedly with the same data input

calculate_rarefaction_curves <- function(psdata, measures, depths) {

  require('plyr') # ldply

  require('reshape2') # melt

  estimate_rarified_richness <- function(psdata, measures, depth) {

    if(max(sample_sums(psdata)) < depth) return()

    psdata <- prune_samples(sample_sums(psdata) >= depth, psdata)

    rarified_psdata <- rarefy_even_depth(psdata, depth, verbose = FALSE)

    alpha_diversity <- estimate_richness(rarified_psdata, measures = measures)

    # as.matrix forces the use of melt.array, which includes the Sample names (rownames)

    molten_alpha_diversity <- melt(as.matrix(alpha_diversity), varnames = c('Sample', 'Measure'), value.name = 'Alpha_diversity')

    molten_alpha_diversity

  }

  names(depths) <- depths # this enables automatic addition of the Depth to the output by ldply

  rarefaction_curve_data <- ldply(depths, estimate_rarified_richness, psdata = psdata, measures = measures, .id = 'Depth', .progress = ifelse(interactive(), 'text', 'none'))

  # convert Depth from factor to numeric

  rarefaction_curve_data$Depth <- as.numeric(levels(rarefaction_curve_data$Depth))[rarefaction_curve_data$Depth]

  rarefaction_curve_data

}

#From here on psdata = BCsoil

rarefaction_curve_data <- calculate_rarefaction_curves(psdata, c('Observed', 'Shannon'), rep(c(1, 10, 100, 1000, 1:100 * 10000), each = 10))

summary(rarefaction_curve_data)

rarefaction_curve_data_summary <- ddply(rarefaction_curve_data, c('Depth', 'Sample', 'Measure'), summarise, Alpha_diversity_mean = mean(Alpha_diversity), Alpha_diversity_sd = sd(Alpha_diversity))

rarefaction_curve_data_summary_verbose <- merge(rarefaction_curve_data_summary, data.frame(sample_data(psdata)), by.x = 'Sample', by.y = 'row.names')

#Create Plot

ggplot(

  data = rarefaction_curve_data_summary_verbose,

  mapping = aes(

    x = Depth,

    y = Alpha_diversity_mean,

    ymin = Alpha_diversity_mean - Alpha_diversity_sd,

    ymax = Alpha_diversity_mean + Alpha_diversity_sd,

    colour = Season,

    group = Sample

  )

) + geom_line(

) + geom_pointrange(

) + facet_wrap(

  facets = ~ Measure,

  scales = 'free_y'

)