San Michele all’Adige (TN), June 25, 2019
curatedMetagenomicData provides 6 types of data for each dataset:
Types 1-3 are generated by MetaPhlAn2; 4-6 are generated by HUMAnN2.
Because curatedMetagenomicData utilizes Bioconductor's ExperimentHub platform, data are stored as Amazon S3 buckets and only downloaded individually when requested.
library(curatedMetagenomicData)
Individual data projects can be fetched via per-dataset functions or the curatedMetagenomicData() function. A function call to a dataset name returns a Bioconductor ExpressionSet object:
loman.eset = LomanNJ_2013.metaphlan_bugs_list.stool()
## snapshotDate(): 2018-10-30
## see ?curatedMetagenomicData and browseVignettes('curatedMetagenomicData') for documentation
## downloading 0 resources
## loading from cache ## '/home/francesco//.ExperimentHub/1278'
## ExpressionSet (storageMode: lockedEnvironment) ## assayData: 736 features, 43 samples ## element names: exprs ## protocolData: none ## phenoData ## sampleNames: OBK1122 OBK1196 ... OBK5066 (43 total) ## varLabels: subjectID body_site ... NCBI_accession (19 total) ## varMetadata: labelDescription ## featureData: none ## experimentData: use 'experimentData(object)' ## pubMedIds: 23571589 ## Annotation:
All datasets are represented as ExpressionSet objects because of the integrative nature of the class and its ability to bind data and metadata.
To access the data itself (in this case relative abundance), the exprs() function returns a variables by samples (rows by columns) numeric matrix. Note the presence of "synthetic" clades at all levels of the taxonomy, starting with kingdom, e.g. k__bacteria here:
exprs( loman.eset )[1:6, 1:5] #first 6 rows and 5 columns
## OBK1122 OBK1196 OBK1253 ## k__Bacteria 100.00000 100.00000 100.00000 ## k__Bacteria|p__Bacteroidetes 83.55662 30.41764 79.41203 ## k__Bacteria|p__Firmicutes 15.14819 68.30191 19.28784 ## k__Bacteria|p__Proteobacteria 1.29520 0.06303 0.97218 ## k__Bacteria|p__Bacteroidetes|c__Bacteroidia 83.55662 30.41764 79.41203 ## k__Bacteria|p__Firmicutes|c__Clostridia 14.93141 33.22720 11.26228 ## OBK2535 OBK2535b ## k__Bacteria 99.39953 99.63084 ## k__Bacteria|p__Bacteroidetes 15.14525 7.03101 ## k__Bacteria|p__Firmicutes 25.71252 45.17325 ## k__Bacteria|p__Proteobacteria 58.53298 47.29051 ## k__Bacteria|p__Bacteroidetes|c__Bacteroidia 15.14525 7.03101 ## k__Bacteria|p__Firmicutes|c__Clostridia 16.50313 13.97996
To access the subject-level metadata the pData() function is used to return a data.frame containing subject-level variables of study.
head( pData( loman.eset ) )
## subjectID body_site antibiotics_current_use study_condition ## OBK1122 OBK1122 stool <NA> STEC ## OBK1196 OBK1196 stool <NA> STEC ## OBK1253 OBK1253 stool <NA> STEC ## OBK2535 OBK2535 stool <NA> STEC ## OBK2535b OBK2535b stool <NA> STEC ## OBK2638 OBK2638 stool <NA> STEC ## disease age_category country non_westernized DNA_extraction_kit ## OBK1122 STEC adult DEU no Qiagen ## OBK1196 STEC adult DEU no Qiagen ## OBK1253 STEC adult DEU no Qiagen ## OBK2535 STEC adult DEU no Qiagen ## OBK2535b STEC adult DEU no Qiagen ## OBK2638 STEC adult DEU no Qiagen ## number_reads number_bases minimum_read_length median_read_length ## OBK1122 464060 69609000 150 150 ## OBK1196 30181380 4527207000 150 150 ## OBK1253 65704818 9855722700 150 150 ## OBK2535 43597736 6539660400 150 150 ## OBK2535b 22253314 3337997100 150 150 ## OBK2638 1105492 165823800 150 150 ## days_after_onset stec_count shigatoxin_2_elisa stool_texture ## OBK1122 NA <NA> <NA> <NA> ## OBK1196 NA <NA> <NA> <NA> ## OBK1253 NA <NA> <NA> <NA> ## OBK2535 3 moderate positive smooth ## OBK2535b NA <NA> <NA> <NA> ## OBK2638 5 high positive bloody ## curator NCBI_accession ## OBK1122 Paolo_Manghi ERR262939;ERR262938 ## OBK1196 Paolo_Manghi ERR262941;ERR262940 ## OBK1253 Paolo_Manghi ERR262942;ERR260476 ## OBK2535 Paolo_Manghi ERR260478;ERR260477 ## OBK2535b Paolo_Manghi ERR260479 ## OBK2638 Paolo_Manghi ERR262943;ERR260480
We can now retrieve all the stool-associated MetaPhlAn2 datasets with the following command:
esl <- curatedMetagenomicData("*metaphlan_bugs_list.stool*", dryrun = FALSE)
The previous command creates a list of 29 ExpressionSet objects. To merge all them in a single ExpressionSet object:
eset <- mergeData(esl) eset
## ExpressionSet (storageMode: lockedEnvironment) ## assayData: 3986 features, 5249 samples ## element names: exprs ## protocolData: none ## phenoData ## sampleNames: ## AsnicarF_2017.metaphlan_bugs_list.stool:MV_FEI1_t1Q14 ## AsnicarF_2017.metaphlan_bugs_list.stool:MV_FEI2_t1Q14 ... ## ZellerG_2014.metaphlan_bugs_list.stool:MMPU99077057ST (5249 ## total) ## varLabels: subjectID body_site ... studyID (74 total) ## varMetadata: labelDescription ## featureData: none ## experimentData: use 'experimentData(object)' ## Annotation:
Here's a direct, exploratory analysis of E. coli prevalence in the Loman dataset.
Create a vector of E. coli relative abundances. This grep call with a "$" at the end selects the only row that ends with "s__Escherichia_coli":
x = exprs( loman.eset )[grep("s__Escherichia_coli$", rownames( loman.eset )), ]
summary( x )
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.0000 0.6047 2.3710 13.6897 11.8291 90.4474
hist( x, xlab = "Relative Abundance", main="Prevalence of E. Coli",
breaks="FD")
phyloseqFor the MetaPhlAn2 bugs datasets (but not other data types), you gain a lot of taxonomy-aware, ecological analysis and plotting by conversion to a phyloseq class object. curatedMetagenomicData provides the ExpressionSet2phyloseq() function to make this easy:
library(phyloseq) loman.pseq = ExpressionSet2phyloseq( loman.eset )
This phyloseq objects contain 3 components, with extractor functions hinted at by its show method:
loman.pseq
## phyloseq-class experiment-level object ## otu_table() OTU Table: [ 736 taxa and 43 samples ] ## sample_data() Sample Data: [ 43 samples by 19 sample variables ] ## tax_table() Taxonomy Table: [ 736 taxa by 8 taxonomic ranks ]
otu_table() returns the same thing as exprs( loman.eset ) did, the Operational Taxanomic Unit (OTU) table. Here are the first 6 rows and 5 columns:
otu_table( loman.pseq )[1:6, 1:5]
## OTU Table: [6 taxa and 5 samples] ## taxa are rows ## OBK1122 OBK1196 OBK1253 OBK2535 OBK2535b ## k__Bacteria 100.00000 100.00000 100.00000 99.39953 99.63084 ## p__Bacteroidetes 83.55662 30.41764 79.41203 15.14525 7.03101 ## p__Firmicutes 15.14819 68.30191 19.28784 25.71252 45.17325 ## p__Proteobacteria 1.29520 0.06303 0.97218 58.53298 47.29051 ## c__Bacteroidia 83.55662 30.41764 79.41203 15.14525 7.03101 ## c__Clostridia 14.93141 33.22720 11.26228 16.50313 13.97996
The same patient or participant data that was available from pData( loman.eset ) is now availble using sample_data() on the phyloseq object:
sample_data( loman.pseq )[1:6, 1:5]
## subjectID body_site antibiotics_current_use study_condition ## OBK1122 OBK1122 stool <NA> STEC ## OBK1196 OBK1196 stool <NA> STEC ## OBK1253 OBK1253 stool <NA> STEC ## OBK2535 OBK2535 stool <NA> STEC ## OBK2535b OBK2535b stool <NA> STEC ## OBK2638 OBK2638 stool <NA> STEC ## disease ## OBK1122 STEC ## OBK1196 STEC ## OBK1253 STEC ## OBK2535 STEC ## OBK2535b STEC ## OBK2638 STEC
But this object also is aware of the taxonomic structure, which will enable the powerful subsetting methods of the phyloseq package.
head( tax_table( loman.pseq ) )
## Taxonomy Table: [6 taxa by 8 taxonomic ranks]: ## Kingdom Phylum Class Order Family ## k__Bacteria "Bacteria" NA NA NA NA ## p__Bacteroidetes "Bacteria" "Bacteroidetes" NA NA NA ## p__Firmicutes "Bacteria" "Firmicutes" NA NA NA ## p__Proteobacteria "Bacteria" "Proteobacteria" NA NA NA ## c__Bacteroidia "Bacteria" "Bacteroidetes" "Bacteroidia" NA NA ## c__Clostridia "Bacteria" "Firmicutes" "Clostridia" NA NA ## Genus Species Strain ## k__Bacteria NA NA NA ## p__Bacteroidetes NA NA NA ## p__Firmicutes NA NA NA ## p__Proteobacteria NA NA NA ## c__Bacteroidia NA NA NA ## c__Clostridia NA NA NA
The process of subsetting begins with the names of taxonomic ranks:
rank_names( loman.pseq )
## [1] "Kingdom" "Phylum" "Class" "Order" "Family" "Genus" "Species" ## [8] "Strain"
Taxa can be filtered by these rank names. For example, to return an object with only species and strains:
subset_taxa( loman.pseq, !is.na(Species) & is.na(Strain))
## phyloseq-class experiment-level object ## otu_table() OTU Table: [ 287 taxa and 43 samples ] ## sample_data() Sample Data: [ 43 samples by 19 sample variables ] ## tax_table() Taxonomy Table: [ 287 taxa by 8 taxonomic ranks ]
Or to keep only the Bacteroidetes phylum with the prune_taxa function
loman.bd = prune_taxa('p__Bacteroidetes', loman.pseq)
head( taxa_names( loman.bd ) )
## [1] "p__Bacteroidetes"
The phyloseq package provides advanced pruning of taxa, such as the following which keeps only taxa that are among the most abundant 5% in at least five samples:
keepotu = genefilter_sample(loman.pseq, filterfun_sample(topp(0.05)), A=5) summary(keepotu)
## Mode FALSE TRUE ## logical 624 112
subset_taxa(loman.pseq, keepotu)
## phyloseq-class experiment-level object ## otu_table() OTU Table: [ 112 taxa and 43 samples ] ## sample_data() Sample Data: [ 43 samples by 19 sample variables ] ## tax_table() Taxonomy Table: [ 112 taxa by 8 taxonomic ranks ]
The phyloseq package provides the plot_heatmap() function to create heatmaps using a variety of built-in dissimilarity metrics for clustering.
Here, we apply the same abundance filter as above, keep only species.
loman.filt = subset_taxa(loman.pseq, keepotu & !is.na(Species) & is.na(Strain))
This function supports a large number of distance and ordination methods, here we use Bray-Curtis dissimilarity for distance and PCoA as the ordination method for organizing the heatmap.
plot_heatmap(loman.filt, method="PCoA", distance="bray")
Here we plot the top 20 most abundant species (not strains), defined by the sum of abundance across all samples in the dataset:
barplot(sort(taxa_sums(loman.filt) / nsamples(loman.filt), TRUE)[1:20],
ylab = "Total counts", las = 2)
The phyloseq package calculates numerous alpha diversity measures. Here we compare three diversity in the species-level data, stratifying by stool texture:
alphas = c("Shannon", "Simpson", "InvSimpson")
plot_richness(loman.filt, "stool_texture", measures = alphas)
Here is a Principal Coordinates Analysis plot of species-level taxa from the Loman dataset, using Bray-Curtis distance:
ordinated_taxa = ordinate(loman.filt, method="PCoA", distance="bray")
plot_ordination(loman.filt, ordinated_taxa, color="stool_texture",
title = "Bray-Curtis Principal Coordinates Analysis")
Using the combined_metadata object we can see all the bodysites and the dataset names included in curatedMetagenomicData
table(combined_metadata$dataset_name)
## ## AsnicarF_2017 BritoIL_2016 Castro_NallarE_2015 ## 24 312 32 ## ChngKR_2016 FengQ_2015 HanniganGD_2017 ## 78 154 82 ## Heitz_BuschartA_2016 HMP_2012 KarlssonFH_2013 ## 53 749 145 ## LeChatelierE_2013 LiJ_2014 LiJ_2017 ## 292 260 196 ## LiuW_2016 LomanNJ_2013 LouisS_2016 ## 110 43 92 ## NielsenHB_2014 Obregon_TitoAJ_2015 OhJ_2014 ## 396 58 291 ## QinJ_2012 QinN_2014 RampelliS_2015 ## 363 237 38 ## RaymondF_2016 SchirmerM_2016 SmitsSA_2017 ## 72 471 40 ## TettAJ_2016 VatanenT_2016 VincentC_2016 ## 97 785 229 ## VogtmannE_2016 XieH_2016 YuJ_2015 ## 110 250 128 ## ZellerG_2014 ## 199
table(combined_metadata$body_site)
## ## milk nasalcavity oralcavity skin stool vagina ## 8 91 651 493 5138 5
Choose a dataset and look if it has different bodysite or different conditions
ds_nielsen <- NielsenHB_2014.metaphlan_bugs_list.stool() table(pData(ds_nielsen)$study_condition)
## ## control IBD ## 248 148
table(pData(ds_nielsen)$bodysite)
## < table of extent 0 >
Plot the heatmap of the 10 most abundant species
ds_nielsen_pseq <- ExpressionSet2phyloseq(ds_nielsen) ds_nielsen_pseq_sp <- subset_taxa(ds_nielsen_pseq, !is.na(Species) & is.na(Strain)) top_10 <- names(sort(taxa_sums(ds_nielsen_pseq_sp), TRUE)[1:10]) ds_nielsen_pseq_sp <- prune_taxa(top_10, ds_nielsen_pseq) plot_heatmap( ds_nielsen_pseq_sp )
## Warning: Transformation introduced infinite values in discrete y-axis
Choose a species and look for its abundance between the bodysites or conditions
relab <- exprs(ds_nielsen)[grep("s__Butyrivibrio_crossotus$", rownames(ds_nielsen)), ]
sample_condt <- pData(ds_nielsen)$study_condition
x <- data.frame(relab, sample_condt)
library(ggplot2)
ggplot(x) +
geom_boxplot(aes(sample_condt, relab))
curatedMetagenomicData includes the three output tables that HUMAnN2 generates when metagenomes are profiled:
LomanNJ_2013.genefamilies_relab.stool()
LomanNJ_2013.pathabundance_relab.stool()
LomanNJ_2013.pathcoverage.stool()
x <- LomanNJ_2013.pathabundance_relab.stool() x=ExpressionSet2phyloseq(x, simplify = FALSE)
vignette(curatedMetagenomicData) is your friend!
Manuscript on Nature Methods www.nature.com/articles/nmeth.4468
Package on Bioconductor doi.org/doi:10.18129/B9.bioc.curatedMetagenomicData