Make phyloseq object

This tutorial is useful for analysis of output files from (Mothur), (QIIME or QIIME2) or any tool that gives a biom file as output. There is also a simple way to read comma seperated (*.csv) files.

Simple comma seperated files:

library(microbiome)


otu.file <-
    system.file("extdata/qiita1629_otu_table.csv",
        package='microbiome')

tax.file <- system.file("extdata/qiita1629_taxonomy_table.csv",
        package='microbiome')

meta.file <- system.file("extdata/qiita1629_mapping_subset.csv",
        package='microbiome')

pseq.csv <- read_phyloseq(
          otu.file=otu.file, 
          taxonomy.file=tax.file, 
          metadata.file=meta.file, type = "simple")

Biom file:

# Read the biom file
biom.file <- 
  system.file("extdata/qiita1629.biom", 
              package = "microbiome")

# Read the mapping/metadata file
 meta.file <- 
  system.file("extdata/qiita1629_mapping.csv", 
              package = "microbiome")
# Make phyloseq object
pseq.biom <- read_phyloseq(otu.file = biom.file, 
                         metadata.file = meta.file, 
                         taxonomy.file = NULL, type = "biom")

Mothur shared OTUs and Consensus Taxonomy:

otu.file <- system.file(
 "extdata/Baxter_FITs_Microbiome_2016_fit.final.tx.1.subsample.shared",
    package='microbiome')

tax.file <- system.file(
 "extdata/Baxter_FITs_Microbiome_2016_fit.final.tx.1.cons.taxonomy",
    package='microbiome')

meta.file <- system.file(
 "extdata/Baxter_FITs_Microbiome_2016_mapping.csv",
    package='microbiome')
 
pseq.mothur <- read_phyloseq(otu.file=otu.file,
        taxonomy.file =tax.file,
        metadata.file=meta.file, type = "mothur")
print(pseq.mothur)

Now, we proceed to core microbiota analysis.

Core microbiota analysis

Here the data from Caporaso, J. Gregory, et al. “Moving pictures of the human microbiome.” Genome biology 12.5 (2011): R50. will be used which is stored as example in jeevanuDB

# install
# install.packages("devtools")
# devtools::install_github("microsud/jeevanuDB")

# check the data 
library(jeevanuDB)
ps <- moving_pictures
table(meta(ps)$sample_type, meta(ps)$host_subject_id)
##         
##           F4  M3
##   saliva 135 373
##   skin   268 723
##   stool  131 336
# Filter the data to include only gut samples from M3 subject
ps.m3 <- subset_samples(ps, sample_type == "stool" & host_subject_id == "M3") 
print(ps.m3)
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3644 taxa and 336 samples ]
## sample_data() Sample Data:       [ 336 samples by 89 sample variables ]
## tax_table()   Taxonomy Table:    [ 3644 taxa by 7 taxonomic ranks ]
# keep only taxa with positive sums
ps.m3 <- prune_taxa(taxa_sums(ps.m3) > 0, ps.m3)
print(ps.m3)
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 611 taxa and 336 samples ]
## sample_data() Sample Data:       [ 336 samples by 89 sample variables ]
## tax_table()   Taxonomy Table:    [ 611 taxa by 7 taxonomic ranks ]
# Calculate compositional version of the data
# (relative abundances)
ps.m3.rel <- microbiome::transform(ps.m3, "compositional")

Output of deblur/dada2 will most likely have seqs as rownames instead of OTU ids or taxa names

taxa_names(ps.m3.rel)[1:3]
## [1] "TACGGAGGATTCAAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGTTTGATAAGTTAGAGGTGAAATACCGGTGCTTAACACCGGAACTGCCTCTAATACTGTTGAACTAGAGAGTAGTTGCGGTAGGCGGAATGTATGG"
## [2] "TACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGTCAGATGTGAAATCCCCGAGCTCAACTTGGGAACTGCGTTTGAAACTACAAGACTAGAATATGTCAGAGGGGGGTAGAATTCCACG"
## [3] "TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGACGGCATAACAAGTCTGATGTGAAAGGCTGGGGCTTAACCCCGGGACTGCATTGGAAACTGTTAAGCTTGAGTGCCGGAGGGGTAAGCGGAATTCCTAG"

We can change it to ASVIDs

library(Biostrings)
dna <- Biostrings::DNAStringSet(taxa_names(ps.m3.rel))
names(dna) <- taxa_names(ps.m3.rel)
ps.m3.rel <- merge_phyloseq(ps.m3.rel, dna)
taxa_names(ps.m3.rel) <- paste0("ASV", seq(ntaxa(ps.m3.rel)))
# now check again
taxa_names(ps.m3.rel)[1:3]
## [1] "ASV1" "ASV2" "ASV3"

Core microbiota analysis

If you only need the names of the core taxa, do as follows. This returns the taxa that exceed the given prevalence and detection thresholds.

core.taxa.standard <- core_members(ps.m3.rel, detection = 0.0001, prevalence = 50/100)

core.taxa.standard

A full phyloseq object of the core microbiota is obtained as follows:

pseq.core <- core(ps.m3.rel, detection = 0.0001, prevalence = .5)

Retrieving the associated taxa names from the phyloseq object:

core.taxa <- taxa(pseq.core)
class(core.taxa)
## [1] "character"
# get the taxonomy data
tax.mat <- tax_table(pseq.core)
tax.df <- as.data.frame(tax.mat)

# add the OTus to last column
tax.df$OTU <- rownames(tax.df)

# select taxonomy of only 
# those OTUs that are core memebers based on the thresholds that were used.
core.taxa.class <- dplyr::filter(tax.df, rownames(tax.df) %in% core.taxa)
knitr::kable(head(core.taxa.class))
Domain Phylum Class Order Family Genus Species OTU
ASV20 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Ruminococcaceae g__Faecalibacterium s__prausnitzii ASV20
ASV34 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Ruminococcaceae g__Oscillospira s__ ASV34
ASV35 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__[Odoribacteraceae] g__Odoribacter s__ ASV35
ASV55 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Lachnospiraceae g__ s__ ASV55
ASV66 k__Bacteria p__Firmicutes c__Clostridia o__Clostridiales f__Lachnospiraceae NA NA ASV66
ASV88 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Rikenellaceae g__ s__ ASV88

Core visualization

Core line plots

Determine core microbiota across various abundance/prevalence thresholds with the blanket analysis (Salonen et al. CMI, 2012) based on various signal and prevalences.

# With compositional (relative) abundances
det <- c(0, 0.1, 0.5, 2, 5, 20)/100
prevalences <- seq(.05, 1, .05)

plot_core(ps.m3.rel, prevalences = prevalences, 
          detections = det, plot.type = "lineplot") + 
  xlab("Relative Abundance (%)") + 
  theme_bw()

Core heatmaps

This visualization method has been used for instance in Intestinal microbiome landscaping: Insight in community assemblage and implications for microbial modulation strategies. Shetty et al. FEMS Microbiology Reviews fuw045, 2017.

Note that you can order the taxa on the heatmap with the order.taxa argument.

# Core with compositionals:
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-4), log10(.2), length = 10)

# Also define gray color palette
gray <- gray(seq(0,1,length=5))
p1 <- plot_core(ps.m3.rel, 
                plot.type = "heatmap", 
                colours = gray,
                prevalences = prevalences, 
                detections = detections, min.prevalence = .5) +
    xlab("Detection Threshold (Relative Abundance (%))")
p1 <- p1 + theme_bw() + ylab("ASVs")
p1

Using viridis color palette

library(viridis)
print(p1 + scale_fill_viridis())
## Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

As it can be seen, we see only OTu IDs and this may not be useful to interpret the data. We need to repreoccess this figure to include taxonomic information. We can do this as follows:

library(RColorBrewer)
library(knitr)
# get the data used for plotting 
df <- p1$data 

# get the list of OTUs
list <- df$Taxa 

# check the OTU ids
# print(list) 

# get the taxonomy data
tax <- as.data.frame(tax_table(ps.m3.rel))

# add the ASVs to last column
tax$ASV <- rownames(tax)

# select taxonomy of only 
# those OTUs that are used in the plot
tax2 <- dplyr::filter(tax, rownames(tax) %in% list) 

# head(tax2)

# We will merege all the column into one except the Doamin as all is bacteria in this case
tax.unit <- tidyr::unite(tax2, Taxa_level,c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species", "ASV"), sep = "_;", remove = TRUE)

tax.unit$Taxa_level <- gsub(pattern="[a-z]__",replacement="", tax.unit$Taxa_level)

# add this new information into the plot data df

df$Taxa <- tax.unit$Taxa_level

# you can see now we have the taxonomic information
knitr::kable(head(df))
Taxa DetectionThreshold Prevalence
Bacteria_;Firmicutes_;Clostridia_;Clostridiales_;Ruminococcaceae_;Faecalibacterium_;prausnitzii_;ASV20 1e-04 0.7797619
Bacteria_;Firmicutes_;Clostridia_;Clostridiales_;Ruminococcaceae_;Oscillospira_;_;ASV34 1e-04 0.7886905
Bacteria_;Bacteroidetes_;Bacteroidia_;Bacteroidales_;[Odoribacteraceae];Odoribacter;_;ASV35 1e-04 0.8869048
Bacteria_;Firmicutes_;Clostridia_;Clostridiales_;Lachnospiraceae_;;;ASV55 1e-04 0.7797619
Bacteria_;Firmicutes_;Clostridia_;Clostridiales_;Lachnospiraceae_;NA_;NA_;ASV66 1e-04 0.5744048
Bacteria_;Bacteroidetes_;Bacteroidia_;Bacteroidales_;Rikenellaceae_;;;ASV88 1e-04 0.5476190 
# replace the data in the plot object
p1$data <- df

plot(p1 + theme(axis.text.y = element_text(face="italic")))

Genus level

ps.m3.rel.gen <- aggregate_taxa(ps.m3.rel, "Genus")
library(RColorBrewer)
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-4), log10(.2), length = 10)

p1 <- plot_core(ps.m3.rel.gen, 
                plot.type = "heatmap", 
                colours = rev(brewer.pal(5, "RdBu")),
                prevalences = prevalences, 
                detections = detections, min.prevalence = .5) +
    xlab("Detection Threshold (Relative Abundance (%))")
p1 <- p1 + theme_bw() + ylab("ASVs")
p1