Microbial Community Diversity Analysis Tutorial with Phyloseq
This is a demo of how to import amplicon microbiome data into R using Phyloseq and run some basic analyses to understand microbial community diversity and composition accross your samples. More demos of this package are available from the authors here. This script was created with Rmarkdown.
Author: Michelle Berry
Updated: April 14, 2016
===================================================================
In this tutorial, we are working with illumina 16s data that has already been processed into an OTU and taxonomy table from the mothur pipeline. Phyloseq has a variety of import options if you processed your raw sequence data with a different pipeline.
The samples were collected from the Western basin of Lake Erie between May and November 2014 at three different locations. The goal of this dataset was to understand how the bacterial community in Lake Erie shifts during toxic algal blooms caused predominantly by a genus of cyanobacteria called Microcystis.
In this tutorial, we will learn how to import an OTU table and sample metadata into R with the Phyloseq package. We will perform some basic exploratory analyses, examining the taxonomic composition of our samples, and visualizing the dissimilarity between our samples in a low-dimensional space using ordinations. Lastly, we will estimate the alpha diversity (richness and evenness) of our samples.
Libraries
#Load libraries
library(ggplot2)
library(vegan)
library(dplyr)
library(scales)
library(grid)
library(reshape2)
library(phyloseq)# Set working directory
setwd("~/chabs/miseq_may2015/analysis/")
# Source code files
# miseqR.R can be found in this repository
source("~/git_repos/MicrobeMiseq/R/miseqR.R")
# Set plotting theme
theme_set(theme_bw())Data import
First, we will import the mothur shared file, consensus taxonomy file, and our sample metadata and store them in one phyloseq object. By storing all of our data structures together in one object we can easily interface between each of the structures. For example, as we will see later, we can use criteria in the sample metadata to select certain samples from the OTU table.
# Assign variables for imported data
sharedfile = "mothur/chabs.shared"
taxfile = "mothur/chabs-silva.taxonomy"
mapfile = "other/habs_metadata_cleaned.csv"
# Import mothur data
mothur_data <- import_mothur(mothur_shared_file = sharedfile,
mothur_constaxonomy_file = taxfile)
# Import sample metadata
map <- read.csv(mapfile)The sample metadata is just a basic csv with columns for sample attributes. Here is a preview of what the sample metadata looks like. As you can see, there is one column called SampleID with the names of each of the samples. The remaining columns contain information on the environmental or sampling conditions related to each sample.
head(map)## SampleID Date Station Shore Fraction Type Samplenum Date_year
## 1 E0010_CNA 6/16 nearshore2 nearshore CNA sample E0010 6/16/14
## 2 E0012_CNA 6/16 nearshore1 nearshore CNA sample E0012 6/16/14
## 3 E0014_CNA 6/16 offshore offshore CNA sample E0014 6/16/14
## 4 E0016_CNA 6/30 nearshore2 nearshore CNA sample E0016 6/30/14
## 5 E0018_CNA 6/30 nearshore1 nearshore CNA sample E0018 6/30/14
## 6 E0020_CNA 6/30 offshore offshore CNA sample E0020 6/30/14
## Month Days pH SpCond StationDepth Secchi Temp Turbidity ParMC DissMC
## 1 June 20 8.080 290.6 5.1 0.9 22.4 10.10 0 NA
## 2 June 20 8.390 300.4 6.1 0.9 23.3 11.90 0 NA
## 3 June 20 8.308 281.7 8.5 4.0 22.5 1.34 0 NA
## 4 June 34 8.647 291.5 4.9 1.3 24.7 4.81 0 NA
## 5 June 34 8.435 290.0 6.0 1.4 24.5 4.01 0 NA
## 6 June 34 8.609 311.8 8.0 2.0 24.8 1.79 0 NA
## Phycocyanin Chla DOC SRP TP TDP POC PON PP N.P Nitrate
## 1 0.31 10.72 3.26025 4.16 36.1 12.3 0.68 0.15 23.84 6.15 482.5
## 2 1.06 34.37 3.40300 1.11 56.8 10.5 1.72 0.32 46.22 6.89 545.0
## 3 0.44 12.26 3.24075 0.24 19.7 4.4 0.84 0.13 15.31 8.31 475.0
## 4 1.69 19.94 3.84425 0.65 26.3 6.7 0.96 0.15 19.69 7.63 472.0
## 5 0.63 7.88 3.02160 0.17 19.4 4.6 0.59 0.12 14.83 7.92 444.5
## 6 1.59 12.93 3.38725 0.27 23.7 5.7 0.78 0.13 18.05 7.24 611.5
## Ammonia H2O2 CDOM
## 1 3.51 217.1213 4.212388
## 2 2.14 180.3269 4.252438
## 3 10.41 224.5438 3.523810
## 4 8.64 190.3620 3.454660
## 5 16.18 167.6386 4.595779
## 6 20.65 199.2388 6.792312We convert this dataframe into phyloseq format with a simple constructor. The only formatting required to merge the sample data into a phyloseq object is that the rownames must match the sample names in your shared and taxonomy files.
map <- sample_data(map)
# Assign rownames to be Sample ID's
rownames(map) <- map$SampleIDWe need to merge our metadata into our phyloseq object.
# Merge mothurdata object with sample metadata
moth_merge <- merge_phyloseq(mothur_data, map)
moth_merge## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 19656 taxa and 56 samples ]
## sample_data() Sample Data: [ 56 samples by 32 sample variables ]
## tax_table() Taxonomy Table: [ 19656 taxa by 6 taxonomic ranks ]Now we have a phyloseq object called moth.merge. If we wanted to, we could also add a phylogenetic tree or a fasta with OTU representative sequences into this object. At anytime, we can print out the data structures stored in a phyloseq object to quickly view its contents.
Before we move on with analysis, we need to do some basic reformatting and filtering.
What are the column names of our taxonomy file?
colnames(tax_table(moth_merge))## [1] "Rank1" "Rank2" "Rank3" "Rank4" "Rank5" "Rank6"These taxonomy names are not helpful, so let’s rename them
colnames(tax_table(moth_merge)) <- c("Kingdom", "Phylum", "Class",
"Order", "Family", "Genus")Now, let’s filter out samples we don’t want to include in our analysis such as the extraction and pcr blanks (We can look at these later to see what’s there)
Note: there is a column in my metadata named “Type”
The “.” that appears in the prune_taxa command is used to refer to the data object that we are piping in (the phyloseq object with non-samples removed)
moth_sub <- moth_merge %>%
subset_samples(Type == "sample") %>%
prune_taxa(taxa_sums(.) > 0, .)Now we will filter out Eukaryotes, Archaea, chloroplasts and mitochondria, because we only intended to amplify bacterial sequences. You may have done this filtering already in mothur, but it’s good to check you don’t have anything lurking in the taxonomy table. I like to keep these organisms in my dataset when running mothur because they are easy enough to remove with Phyloseq and sometimes I’m interested in exploring them.
erie <- moth_sub %>%
subset_taxa(
Kingdom == "Bacteria" &
Family != "mitochondria" &
Class != "Chloroplast"
)
erie## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 7522 taxa and 56 samples ]
## sample_data() Sample Data: [ 56 samples by 32 sample variables ]
## tax_table() Taxonomy Table: [ 7522 taxa by 6 taxonomic ranks ]Sample summary
As a first analysis, we will look at the distribution of read counts from our samples
# Make a data frame with a column for the read counts of each sample
sample_sum_df <- data.frame(sum = sample_sums(erie))
# Histogram of sample read counts
ggplot(sample_sum_df, aes(x = sum)) +
geom_histogram(color = "black", fill = "indianred", binwidth = 2500) +
ggtitle("Distribution of sample sequencing depth") +
xlab("Read counts") +
theme(axis.title.y = element_blank())
# mean, max and min of sample read counts
smin <- min(sample_sums(erie))
smean <- mean(sample_sums(erie))
smax <- max(sample_sums(erie))The minimum sample read count is 1.563110^{4}
The mean sample read count is 3.0975310^{4}
The max sample read count is 4.777110^{4}
Stacked barplots
Let’s make a stacked barplot of Phyla to get a sense of the community composition in these samples.
Since this is not a quantitative analysis, and since we have more Phyla in this dataset than we can reasonably distinguish colors (43!), we will prune out low abundance taxa and only include Phyla that contribute more than 2% of the relative abundance of each sample. Depending on your dataset and the taxonomic level you are depicting, you can adjust this prune parameter. In later analyses, we will of course included these taxa, but for now they will just clutter our plot.
# melt to long format (for ggploting)
# prune out phyla below 2% in each sample
erie_phylum <- erie %>%
tax_glom(taxrank = "Phylum") %>% # agglomerate at phylum level
transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
psmelt() %>% # Melt to long format
filter(Abundance > 0.02) %>% # Filter out low abundance taxa
arrange(Phylum) # Sort data frame alphabetically by phylum# Set colors for plotting
phylum_colors <- c(
"#CBD588", "#5F7FC7", "orange","#DA5724", "#508578", "#CD9BCD",
"#AD6F3B", "#673770","#D14285", "#652926", "#C84248",
"#8569D5", "#5E738F","#D1A33D", "#8A7C64", "#599861"
)
# Plot
ggplot(erie_phylum, aes(x = Date, y = Abundance, fill = Phylum)) +
facet_grid(Station~.) +
geom_bar(stat = "identity") +
scale_fill_manual(values = phylum_colors) +
scale_x_discrete(
breaks = c("7/8", "8/4", "9/2", "10/6"),
labels = c("Jul", "Aug", "Sep", "Oct"),
drop = FALSE
) +
# Remove x axis title
theme(axis.title.x = element_blank()) +
#
guides(fill = guide_legend(reverse = TRUE, keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance (Phyla > 2%) \n") +
ggtitle("Phylum Composition of Lake Erie \n Bacterial Communities by Sampling Site") 
This plot was created using facets to seperate samples along the y axis by sampling station. This is a great feature of ggplot.
Notice that each sample doesn’t fully add up to 1. This reflects the rare phyla that were removed. If you want your plot to look like everything adds up, you can add position = “fill” to the geom_bar() command.
Unconstrained Ordinations
One of the best exploratory analyses for amplicon data is unconstrained ordinations. Here we will look at ordinations of our full community samples. We will use the scale_reads() function in miseqR.R to scale to the smallest library size, which is the default. If you want to scale to another depth, you can do so by setting the “n” argument.
# Scale reads to even depth
erie_scale <- erie %>%
scale_reads(round = "round")
# Fix month levels in sample_data
sample_data(erie_scale)$Month <- factor(
sample_data(erie_scale)$Month,
levels = c("June", "July", "August", "September", "October")
)
# Ordinate
erie_pcoa <- ordinate(
physeq = erie_scale,
method = "PCoA",
distance = "bray"
)
# Plot
plot_ordination(
physeq = erie_scale,
ordination = erie_pcoa,
color = "Month",
shape = "Station",
title = "PCoA of Lake Erie bacterial Communities"
) +
scale_color_manual(values = c("#a65628", "red", "#ffae19",
"#4daf4a", "#1919ff", "darkorchid3", "magenta")
) +
geom_point(aes(color = Month), alpha = 0.7, size = 4) +
geom_point(colour = "grey90", size = 1.5) 
Let’s try an NMDS instead. For NMDS plots it’s important to set a seed since the starting positions of samples in the alogrithm is random.
Important: if you calculate your bray-curtis distance metric “in-line” it will perform a square root transformation and Wisconsin double standardization. If you don’t want this, you can calculate your bray-curtis distance separately
set.seed(1)
# Ordinate
erie_nmds <- ordinate(
physeq = erie_scale,
method = "NMDS",
distance = "bray"
)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1579025
## Run 1 stress 0.1819437
## Run 2 stress 0.240828
## Run 3 stress 0.1866267
## Run 4 stress 0.185961
## Run 5 stress 0.2000213
## Run 6 stress 0.1710585
## Run 7 stress 0.1711768
## Run 8 stress 0.2248862
## Run 9 stress 0.1580173
## ... procrustes: rmse 0.005865558 max resid 0.02795903
## Run 10 stress 0.1579805
## ... procrustes: rmse 0.00415466 max resid 0.02854737
## Run 11 stress 0.210445
## Run 12 stress 0.1579419
## ... procrustes: rmse 0.002487412 max resid 0.01392694
## Run 13 stress 0.1579215
## ... procrustes: rmse 0.004302403 max resid 0.02680986
## Run 14 stress 0.2153886
## Run 15 stress 0.2408162
## Run 16 stress 0.2032362
## Run 17 stress 0.1710288
## Run 18 stress 0.214908
## Run 19 stress 0.1819041
## Run 20 stress 0.1711611# Plot
plot_ordination(
physeq = erie_scale,
ordination = erie_nmds,
color = "Month",
shape = "Station",
title = "NMDS of Lake Erie bacterial Communities"
) +
scale_color_manual(values = c("#a65628", "red", "#ffae19",
"#4daf4a", "#1919ff", "darkorchid3", "magenta")
) +
geom_point(aes(color = Month), alpha = 0.7, size = 4) +
geom_point(colour = "grey90", size = 1.5) 
NMDS plots attempt to show ordinal distances between samples as accurately as possible in two dimensions. It is important to report the stress of these plots, because a high stress value means that the algorithm had a hard time representing the distances between samples in 2 dimensions. The stress of this plot was OK - it was .148 (generally anything below .2 is considered acceptable). However, the PCoA for this data was able to show a lot of variation in just two dimensions, and it shows the temporal trends in this dataset better, so we will stick with that plot.
Permanova
Here is an example of how to run a permanova test using the adonis function in vegan. In this example we are testing the hypothesis that the three stations we collected samples from have different centroids
set.seed(1)
# Calculate bray curtis distance matrix
erie_bray <- phyloseq::distance(erie_scale, method = "bray")
# make a data frame from the sample_data
sampledf <- data.frame(sample_data(erie))
# Adonis test
adonis(erie_bray ~ Station, data = sampledf)##
## Call:
## adonis(formula = erie_bray ~ Station, data = sampledf)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Station 2 0.6754 0.33772 2.7916 0.09531 0.003 **
## Residuals 53 6.4118 0.12098 0.90469
## Total 55 7.0872 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Homogeneity of dispersion test
beta <- betadisper(erie_bray, sampledf$Station)
permutest(beta)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.025388 0.012694 2.3781 999 0.114
## Residuals 53 0.282915 0.005338This output tells us that our adonis test is significant so we can reject the null hypothesis that our three sites have the same centroid.
Additionally, our betadisper results are not significant, meaning we cannot reject the null hypothesis that our groups have the same dispersions. This means we can be more confident that our adonis result is a real result, and not due to differences in group dispersions
There is a lot more analysis that can be done here. We could use a distance metric other than Bray-curtis, we could test different grouping variables, or we could create a more complex permanova by testing a model that combines multiple variables. Unfortunately, there are currently no post-hoc tests developed for adonis.
Constrained Ordinations
Above we used unconstrained ordinations (PCoA, NMDS) to show relationships between samples in low dimensions. We can use a constrained ordination to see how environmental variables are associated with these changes in community composition. We constrain the ordination axes to linear combinations of environmental variables. We then plot the environmental scores onto the ordination
# Remove data points with missing metadata
erie_not_na <- erie_scale %>%
subset_samples(
!is.na(Phycocyanin) &
!is.na(SRP) &
!is.na(pH) &
!is.na(ParMC) &
!is.na(H2O2)
)
bray_not_na <- phyloseq::distance(physeq = erie_not_na, method = "bray")
# CAP ordinate
cap_ord <- ordinate(
physeq = erie_not_na,
method = "CAP",
distance = bray_not_na,
formula = ~ ParMC + Nitrate + SRP + Phycocyanin + Ammonia + pH + H2O2
)
# CAP plot
cap_plot <- plot_ordination(
physeq = erie_not_na,
ordination = cap_ord,
color = "Month",
axes = c(1,2)
) +
aes(shape = Station) +
geom_point(aes(colour = Month), alpha = 0.4, size = 4) +
geom_point(colour = "grey90", size = 1.5) +
scale_color_manual(values = c("#a65628", "red", "#ffae19", "#4daf4a",
"#1919ff", "darkorchid3", "magenta")
)
# Now add the environmental variables as arrows
arrowmat <- vegan::scores(cap_ord, display = "bp")
# Add labels, make a data.frame
arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)
# Define the arrow aesthetic mapping
arrow_map <- aes(xend = CAP1,
yend = CAP2,
x = 0,
y = 0,
shape = NULL,
color = NULL,
label = labels)
label_map <- aes(x = 1.3 * CAP1,
y = 1.3 * CAP2,
shape = NULL,
color = NULL,
label = labels)
arrowhead = arrow(length = unit(0.02, "npc"))
# Make a new graphic
cap_plot +
geom_segment(
mapping = arrow_map,
size = .5,
data = arrowdf,
color = "gray",
arrow = arrowhead
) +
geom_text(
mapping = label_map,
size = 4,
data = arrowdf,
show.legend = FALSE
)
Do a permutational ANOVA on constrained axes used in ordination
anova(cap_ord)## Permutation test for capscale under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = distance ~ ParMC + Nitrate + SRP + Phycocyanin + Ammonia + pH + H2O2, data = data)
## Df SumOfSqs F Pr(>F)
## Model 7 3.1215 5.6409 0.001 ***
## Residual 45 3.5574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Alpha Diversity
Estimating alpha diversity of microbial communities is problematic no matter what you do. My best stab at it is to subsample the libraries with replacement to estimate the species abundance of the real population while standardizing sampling effort.
min_lib <- min(sample_sums(erie))We will subsample to 1.563110^{4}, the minimum number of reads. We will repeat this 100 times and average the diversity estimates from each trial.
# Initialize matrices to store richness and evenness estimates
nsamp = nsamples(erie)
trials = 100
richness <- matrix(nrow = nsamp, ncol = trials)
row.names(richness) <- sample_names(erie)
evenness <- matrix(nrow = nsamp, ncol = trials)
row.names(evenness) <- sample_names(erie)
# It is always important to set a seed when you subsample so your result is replicable
set.seed(3)
for (i in 1:100) {
# Subsample
r <- rarefy_even_depth(erie, sample.size = min_lib, verbose = FALSE, replace = TRUE)
# Calculate richness
rich <- as.numeric(as.matrix(estimate_richness(r, measures = "Observed")))
richness[ ,i] <- rich
# Calculate evenness
even <- as.numeric(as.matrix(estimate_richness(r, measures = "InvSimpson")))
evenness[ ,i] <- even
}Let’s calculate the mean and standard deviation per sample for observed richness and inverse simpson’s index and store those values in a dataframe.
# Create a new dataframe to hold the means and standard deviations of richness estimates
SampleID <- row.names(richness)
mean <- apply(richness, 1, mean)
sd <- apply(richness, 1, sd)
measure <- rep("Richness", nsamp)
rich_stats <- data.frame(SampleID, mean, sd, measure)
# Create a new dataframe to hold the means and standard deviations of evenness estimates
SampleID <- row.names(evenness)
mean <- apply(evenness, 1, mean)
sd <- apply(evenness, 1, sd)
measure <- rep("Inverse Simpson", nsamp)
even_stats <- data.frame(SampleID, mean, sd, measure)Now we will combine our estimates for richness and evenness into one dataframe
alpha <- rbind(rich_stats, even_stats)Let’s add the sample metadata into this dataframe using the merge() command
s <- data.frame(sample_data(erie))
alphadiv <- merge(alpha, s, by = "SampleID") Lastly, we will reorder some factors in this dataset before plotting them
alphadiv <- order_dates(alphadiv)Finally, we will plot the two alpha diversity measures in a timeseries using a facet
ggplot(alphadiv, aes(x = Date, y = mean, color = Station, group = Station, shape = Station)) +
geom_point(size = 2) +
geom_line(size = 0.8) +
facet_wrap(~measure, ncol = 1, scales = "free") +
scale_color_manual(values = c("#E96446", "#302F3D", "#87CEFA")) +
scale_x_discrete(
breaks = c("7/8", "8/4", "9/2", "10/6"),
labels = c("Jul", "Aug", "Sep", "Oct"),
drop = FALSE
) +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
Session info
sessionInfo()## R version 3.2.2 (2015-08-14)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.9.5 (Mavericks)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] phyloseq_1.12.2 reshape2_1.4.1 scales_0.4.0 dplyr_0.4.3
## [5] vegan_2.3-1 lattice_0.20-33 permute_0.8-4 ggplot2_2.1.0
##
## loaded via a namespace (and not attached):
## [1] Biobase_2.28.0 splines_3.2.2
## [3] foreach_1.4.3 Formula_1.2-1
## [5] assertthat_0.1 stats4_3.2.2
## [7] latticeExtra_0.6-26 yaml_2.1.13
## [9] RSQLite_1.0.0 chron_2.3-47
## [11] digest_0.6.9 GenomicRanges_1.20.8
## [13] RColorBrewer_1.1-2 XVector_0.8.0
## [15] colorspace_1.2-6 htmltools_0.2.6
## [17] Matrix_1.2-2 plyr_1.8.3
## [19] DESeq2_1.8.2 XML_3.98-1.3
## [21] genefilter_1.50.0 zlibbioc_1.14.0
## [23] xtable_1.7-4 BiocParallel_1.2.22
## [25] annotate_1.46.1 mgcv_1.8-7
## [27] IRanges_2.2.9 lazyeval_0.1.10
## [29] nnet_7.3-11 BiocGenerics_0.14.0
## [31] proto_0.3-10 survival_2.38-3
## [33] RJSONIO_1.3-0 magrittr_1.5
## [35] evaluate_0.8 nlme_3.1-122
## [37] MASS_7.3-44 RcppArmadillo_0.6.100.0.0
## [39] foreign_0.8-66 tools_3.2.2
## [41] data.table_1.9.6 formatR_1.2.1
## [43] stringr_1.0.0 S4Vectors_0.6.6
## [45] locfit_1.5-9.1 munsell_0.4.3
## [47] cluster_2.0.3 AnnotationDbi_1.30.1
## [49] lambda.r_1.1.7 Biostrings_2.36.4
## [51] ade4_1.7-2 GenomeInfoDb_1.4.3
## [53] futile.logger_1.4.1 iterators_1.0.8
## [55] biom_0.3.12 igraph_1.0.1
## [57] labeling_0.3 rmarkdown_0.9.5
## [59] multtest_2.24.0 gtable_0.2.0
## [61] codetools_0.2-14 DBI_0.3.1
## [63] R6_2.1.1 gridExtra_2.0.0
## [65] knitr_1.11 Hmisc_3.17-0
## [67] futile.options_1.0.0 ape_3.3
## [69] stringi_1.0-1 parallel_3.2.2
## [71] Rcpp_0.12.3 geneplotter_1.46.0
## [73] rpart_4.1-10 acepack_1.3-3.3