Contents

1 What curatedMetagenomicData provides

curatedMetagenomicData provides 6 types of data for each dataset:

  1. Species-level taxonomic profiles, expressed as relative abundance from kingdom to strain level
  2. Presence of unique, clade-specific markers
  3. Abundance of unique, clade-specific markers
  4. Abundance of gene families
  5. Metabolic pathway coverage
  6. Metabolic pathway abundance

Types 1-3 are generated by MetaPhlAn2; 4-6 are generated by HUMAnN2.

Currently, curatedMetagenomicData provides:

These datasets are documented in the reference manual.

2 Using curatedMetagenomicData Resources

Because curatedMetagenomicData utilizes Bioconductor’s ExperimentHub platform, data are stored as Amazon S3 buckets and only downloaded individually when requested. Local caching occurs automatically; see ExperimentHub documentation to customize options such as where the cache is stored (by default it is ~/.ExperimentHub). In case you haven’t installed it yet, do:

# BiocInstaller::biocLite("curatedMetagenomicData")
library(curatedMetagenomicData)

Note that newly datasets require installing the development version of Bioconductor.

2.1 Available samples and metadata

The manually curated metadata for all available samples are provided in a single table combined_metadata:

?combined_metadata
View(combined_metadata)
table(combined_metadata$antibiotics_current_use)
## 
##   no  yes 
## 2183  558
table(combined_metadata$disease)
## 
##                                   AD                                AD;AR 
##                                   10                                   16 
##                         AD;AR;asthma                            AD;asthma 
##                                    8                                    4 
##                                   AR                            AR;asthma 
##                                    2                                    8 
##                     CD;schizophrenia                               CDI;NA 
##                                    1                                   14 
##                       CDI;cellulitis                   CDI;osteoarthritis 
##                                    2                                    1 
##                        CDI;pneumonia                    CDI;ureteralstone 
##                                   15                                    1 
##                                  CRC                 CRC;T2D;hypertension 
##                                  258                                    2 
##                      CRC;fatty_liver     CRC;fatty_liver;T2D;hypertension 
##                                    3                                    4 
##         CRC;fatty_liver;hypertension                     CRC;hypertension 
##                                   12                                   12 
##            HBV;HDV;ascites;cirrhosis                    HBV;HDV;cirrhosis 
##                                    2                                    3 
##                     HBV;HE;cirrhosis            HBV;HEV;ascites;cirrhosis 
##                                    2                                    4 
##                    HBV;HEV;cirrhosis                HBV;ascites;cirrhosis 
##                                    3                                   44 
##    HBV;ascites;cirrhosis;schistosoma         HBV;ascites;cirrhosis;wilson 
##                                    1                                    1 
##                        HBV;cirrhosis                        HCV;cirrhosis 
##                                   39                                    1 
##                     HE;HEV;cirrhosis    HEV;ascites;cirrhosis;schistosoma 
##                                    1                                    2 
##                                  IBD                                  IGT 
##                                  148                                   49 
##                                   NK                                 STEC 
##                                    1                                   43 
##                                  T1D                      T1D;coeliac;IBD 
##                                   24                                    3 
##                                  T2D                     T2D;hypertension 
##                                  223                                    1 
##                    T2D;schizophrenia                      abdominalhernia 
##                                    3                                    2 
##                              adenoma                          adenoma;T2D 
##                                   77                                    1 
##             adenoma;T2D;hypertension                  adenoma;fatty_liver 
##                                    2                                   12 
##              adenoma;fatty_liver;T2D adenoma;fatty_liver;T2D;hypertension 
##                                    1                                    1 
##     adenoma;fatty_liver;hypertension                 adenoma;hypertension 
##                                   14                                    8 
##                            arthritis                    ascites;cirrhosis 
##                                    1                                    9 
##        ascites;cirrhosis;pbcirrhosis        ascites;cirrhosis;schistosoma 
##                                    1                                    1 
##             ascites;cirrhosis;wilson                           bronchitis 
##                                    1                                   18 
##                           cellulitis                            cirrhosis 
##                                    6                                    8 
##     coeliac;gestational_diabetes;CMV                                cough 
##                                    2                                    2 
##                             cystitis                          fatty_liver 
##                                    1                                    8 
##                      fatty_liver;T2D         fatty_liver;T2D;hypertension 
##                                    3                                    9 
##             fatty_liver;hypertension                                fever 
##                                   13                                    3 
##                             gangrene                              healthy 
##                                    1                                 3366 
##                            hepatitis                         hypertension 
##                                    3                                  166 
##            infectiousgastroenteritis                       osteoarthritis 
##                                    5                                    1 
##                               otitis                            pneumonia 
##                                  107                                    7 
##                            psoriasis                  psoriasis;arthritis 
##                                   74                                   22 
##                        pyelonefritis                       pyelonephritis 
##                                    2                                    6 
##                       respiratoryinf                        salmonellosis 
##                                   13                                    1 
##                        schizophrenia                               sepsis 
##                                   12                                    1 
##                              skininf                           stomatitis 
##                                    2                                    2 
##                              suspinf                          tonsillitis 
##                                    1                                    3

2.1.1 Read depth of all samples across all studies

combined_metadata also provides technical information for each sample like sequencing platform, read length, and read depth. The following uses combined_metadata to create a boxplot of read depth for each sample in each study, with boxes colored by body site. First, create a ranking of datasets by median read depth:

dsranking <- combined_metadata %>%
  as.data.frame() %>% 
  group_by(dataset_name) %>%
  summarize(mediandepth = median(number_reads) / 1e6) %>%
  mutate(dsorder = rank(mediandepth)) %>%
  arrange(dsorder)
dsranking
## # A tibble: 31 x 3
##    dataset_name    mediandepth dsorder
##    <chr>                 <dbl>   <dbl>
##  1 TettAJ_2016            1.19       1
##  2 SmitsSA_2017           3.81       2
##  3 HanniganGD_2017        5.86       3
##  4 LomanNJ_2013           6.66       4
##  5 OhJ_2014              10.5        5
##  6 ChngKR_2016           14.1        6
##  7 VincentC_2016         15.0        7
##  8 RampelliS_2015        15.2        8
##  9 VatanenT_2016         19.1        9
## 10 LouisS_2016           20.1       10
## # ... with 21 more rows

Create a factor ds that is ordered according to the ranking by median read depth, to show datasets in order from lowest to highest median read depth, then create the box plot.

suppressPackageStartupMessages(library(ggplot2))
combined_metadata %>%
  as.data.frame() %>%
  mutate(ds = factor(combined_metadata$dataset_name, levels=dsranking$dataset_name)) %>%
  ggplot(aes(ds, number_reads / 1e6, fill=body_site)) + 
  geom_boxplot() +
  theme(axis.text.x = element_text(angle=45, hjust=1)) +
  labs(x="Dataset", y="Read Depth (millions)")

For the datasets that provide both stool and oral cavity profiles, notice how much lower the read depth of oral cavity profiles is due to the higher proportions of removed human DNA reads.

3 Accessing datasets

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:

suppressPackageStartupMessages(library(curatedMetagenomicData))
loman.eset = LomanNJ_2013.metaphlan_bugs_list.stool()

However, the above approach lacks additional options and versioning provided by the curatedMetagenomicData function, which returns a list of datasets. In this case the list has length 1:

loman <- curatedMetagenomicData("LomanNJ_2013.metaphlan_bugs_list.stool", dryrun = FALSE)
## Working on LomanNJ_2013.metaphlan_bugs_list.stool
## snapshotDate(): 2018-06-29
## see ?curatedMetagenomicData and browseVignettes('curatedMetagenomicData') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1278'
loman
## List of length 1
## names(1): LomanNJ_2013.metaphlan_bugs_list.stool
loman.eset <- loman[[1]]
loman.eset
## 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:

The following creates a list of two ExpressionSet objects providing the BritoIL_2016 and Castro-NallarE_2015 oral cavity taxonomic profiles:

oral <- c("BritoIL_2016.metaphlan_bugs_list.oralcavity",
          "Castro-NallarE_2015.metaphlan_bugs_list.oralcavity")
esl <- curatedMetagenomicData(oral, dryrun = FALSE)
esl
esl[[1]]
esl[[2]]

And the following would provide all stool metaphlan datasets if dryrun = FALSE were set:

curatedMetagenomicData("*metaphlan_bugs_list.stool*", dryrun = TRUE)
## Dry run: see return values for datasets that would be downloaded. Run with `dryrun=FALSE` to actually download these datasets.
##  [1] "AsnicarF_2017.metaphlan_bugs_list.stool"       
##  [2] "BritoIL_2016.metaphlan_bugs_list.stool"        
##  [3] "FengQ_2015.metaphlan_bugs_list.stool"          
##  [4] "HMP_2012.metaphlan_bugs_list.stool"            
##  [5] "HanniganGD_2017.metaphlan_bugs_list.stool"     
##  [6] "Heitz-BuschartA_2016.metaphlan_bugs_list.stool"
##  [7] "Heitz_BuschartA_2016.metaphlan_bugs_list.stool"
##  [8] "KarlssonFH_2013.metaphlan_bugs_list.stool"     
##  [9] "LeChatelierE_2013.metaphlan_bugs_list.stool"   
## [10] "LiJ_2014.metaphlan_bugs_list.stool"            
## [11] "LiJ_2017.metaphlan_bugs_list.stool"            
## [12] "LiuW_2016.metaphlan_bugs_list.stool"           
## [13] "LomanNJ_2013.metaphlan_bugs_list.stool"        
## [14] "LouisS_2016.metaphlan_bugs_list.stool"         
## [15] "NielsenHB_2014.metaphlan_bugs_list.stool"      
## [16] "Obregon-TitoAJ_2015.metaphlan_bugs_list.stool" 
## [17] "Obregon_TitoAJ_2015.metaphlan_bugs_list.stool" 
## [18] "QinJ_2012.metaphlan_bugs_list.stool"           
## [19] "QinN_2014.metaphlan_bugs_list.stool"           
## [20] "RampelliS_2015.metaphlan_bugs_list.stool"      
## [21] "RaymondF_2016.metaphlan_bugs_list.stool"       
## [22] "SchirmerM_2016.metaphlan_bugs_list.stool"      
## [23] "SmitsSA_2017.metaphlan_bugs_list.stool"        
## [24] "VatanenT_2016.metaphlan_bugs_list.stool"       
## [25] "VincentC_2016.metaphlan_bugs_list.stool"       
## [26] "VogtmannE_2016.metaphlan_bugs_list.stool"      
## [27] "XieH_2016.metaphlan_bugs_list.stool"           
## [28] "YuJ_2015.metaphlan_bugs_list.stool"            
## [29] "ZellerG_2014.metaphlan_bugs_list.stool"

3.1 Merging multiple datasets

The following merges the two oral cavity datasets downloaded above into a single ExpressionSet.

eset <- mergeData(esl)
eset
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 1914 features, 172 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: BritoIL_2016.metaphlan_bugs_list.oralcavity:M1.1.SA
##     BritoIL_2016.metaphlan_bugs_list.oralcavity:M1.10.SA ...
##     Castro-NallarE_2015.metaphlan_bugs_list.oralcavity:ES_080 (172
##     total)
##   varLabels: subjectID body_site ... studyID (21 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:

This works for any number of datasets. The function will not stop you from merging different data types (e.g. metaphlan bugs lists with gene families), but you probably don’t want to do that.

4 Using ExpressionSet Objects

All datasets are represented as ExpressionSet objects because of the integrative nature of the class and its ability to bind data and metadata. There are three main functions, from the Biobase package, that provide access to experiment-level metadata, subject-level metadata, and the data itself.

To access the experiment-level metadata the experimentData() function is used to return a MIAME (Minimum Information About a Microarray Experiment) object.

experimentData( loman.eset )
## Experiment data
##   Experimenter name: Loman NJ, Constantinidou C, Christner M, Rohde H, Chan JZ, Quick J, Weir JC, Quince C, Smith GP, Betley JR, Aepfelbacher M, Pallen MJ 
##   Laboratory: [1] Institute of Microbiology and Infection, University of Birmingham, Birmingham, England. 
##   Contact information:  
##   Title: A culture-independent sequence-based metagenomics approach to the investigation of an outbreak of Shiga-toxigenic Escherichia coli O104:H4. 
##   URL:  
##   PMIDs: 23571589 
## 
##   Abstract: A 308 word abstract is available. Use 'abstract' method.
##   notes:
##    Sequencing platform:     
##       IlluminaHiSeq

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 disease
## OBK1122    OBK1122     stool                    <NA>            STEC    STEC
## OBK1196    OBK1196     stool                    <NA>            STEC    STEC
## OBK1253    OBK1253     stool                    <NA>            STEC    STEC
## OBK2535    OBK2535     stool                    <NA>            STEC    STEC
## OBK2535b  OBK2535b     stool                    <NA>            STEC    STEC
## OBK2638    OBK2638     stool                    <NA>            STEC    STEC
##          age_category country non_westernized DNA_extraction_kit
## OBK1122         adult     DEU              no             Qiagen
## OBK1196         adult     DEU              no             Qiagen
## OBK1253         adult     DEU              no             Qiagen
## OBK2535         adult     DEU              no             Qiagen
## OBK2535b        adult     DEU              no             Qiagen
## OBK2638         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

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

Bioconductor provides further documentation of the ExpressionSet class and has published an excellent introduction.

4.1 Estimating Absolute Raw Count Data

Absolute raw count data can be estimated from the relative count data by multiplying the columns of the ExpressionSet data by the number of reads for each sample, as found in the pData column “number_reads”. For demo purposes you could (but don’t have to!) do this manually by dividing by 100 and multiplying by the number of reads:

loman.counts = sweep(exprs( loman.eset ), 2, loman.eset$number_reads / 100, "*")
loman.counts = round(loman.counts)
loman.counts[1:6, 1:5]
##                                             OBK1122  OBK1196  OBK1253
## k__Bacteria                                  464060 30181380 65704818
## k__Bacteria|p__Bacteroidetes                 387753  9180464 52177530
## k__Bacteria|p__Firmicutes                     70297 20614459 12673040
## k__Bacteria|p__Proteobacteria                  6011    19023   638769
## k__Bacteria|p__Bacteroidetes|c__Bacteroidia  387753  9180464 52177530
## k__Bacteria|p__Firmicutes|c__Clostridia       69291 10028427  7399861
##                                              OBK2535 OBK2535b
## k__Bacteria                                 43335945 22171164
## k__Bacteria|p__Bacteroidetes                 6602986  1564633
## k__Bacteria|p__Firmicutes                   11210077 10052545
## k__Bacteria|p__Proteobacteria               25519054 10523706
## k__Bacteria|p__Bacteroidetes|c__Bacteroidia  6602986  1564633
## k__Bacteria|p__Firmicutes|c__Clostridia      7194991  3111004

or just set the counts argument in curatedMetagenomicData() to TRUE:

loman.eset2 = curatedMetagenomicData("LomanNJ_2013.metaphlan_bugs_list.stool",
                                     counts = TRUE, dryrun = FALSE)[[1]]
## Working on LomanNJ_2013.metaphlan_bugs_list.stool
## snapshotDate(): 2018-06-29
## see ?curatedMetagenomicData and browseVignettes('curatedMetagenomicData') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1278'
all.equal(exprs(loman.eset2), loman.counts)
## [1] TRUE

5 E. coli prevalence

Here’s a direct, exploratory analysis of E. coli prevalence in the Loman dataset using the ExpressionSet object. More elegant solutions will be provided later using subsetting methods provided by the phyloseq package, but for users familiar with grep() and the ExpressionSet object, such manual methods may suffice.

First, which E. coli-related taxa are available?

grep("coli", rownames(loman.eset), value=TRUE)
## [1] "k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacteriales|f__Enterobacteriaceae|g__Escherichia|s__Escherichia_coli"                                 
## [2] "k__Bacteria|p__Proteobacteria|c__Gammaproteobacteria|o__Enterobacteriales|f__Enterobacteriaceae|g__Escherichia|s__Escherichia_coli|t__Escherichia_coli_unclassified"
## [3] "k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Anaerotruncus|s__Anaerotruncus_colihominis"                                          
## [4] "k__Bacteria|p__Firmicutes|c__Clostridia|o__Clostridiales|f__Ruminococcaceae|g__Anaerotruncus|s__Anaerotruncus_colihominis|t__GCF_000154565"

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

This could be plotted as a histogram:

hist( x, xlab = "Relative Abundance", main="Prevalence of E. Coli",
      breaks="FD")

6 Taxonomy-Aware Analysis using phyloseq

For 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:

suppressPackageStartupMessages(library(phyloseq))
loman.pseq = ExpressionSet2phyloseq( loman.eset )

Note the simplified row names of the OTU table, showing only the most detailed level of the taxonomy. This results from the default argument simplify=TRUE, which is convenient and lossless because taxonomic information is now attainable by tax_table(loman.pseq2)

6.1 phylogenetic trees and UniFrac distances

Set phylogenetictree = TRUE to include a phylogenetic tree in the phyloseq object:

loman.tree <- ExpressionSet2phyloseq( loman.eset, phylogenetictree = TRUE)
wt = UniFrac(loman.tree, weighted=TRUE, normalized=FALSE, 
             parallel=FALSE, fast=TRUE)
plot(hclust(wt), main="Weighted UniFrac distances")

6.2 Components of a phyloseq object

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 disease
## OBK1122    OBK1122     stool                    <NA>            STEC    STEC
## OBK1196    OBK1196     stool                    <NA>            STEC    STEC
## OBK1253    OBK1253     stool                    <NA>            STEC    STEC
## OBK2535    OBK2535     stool                    <NA>            STEC    STEC
## OBK2535b  OBK2535b     stool                    <NA>            STEC    STEC
## OBK2638    OBK2638     stool                    <NA>            STEC    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

6.3 Subsetting / Pruning

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))
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 535 taxa and 43 samples ]
## sample_data() Sample Data:       [ 43 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 535 taxa by 8 taxonomic ranks ]

To keep only phylum-level data (not class or finer, and not kingdom-level):

subset_taxa( loman.pseq, is.na(Class) & !is.na(Phylum))
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 8 taxa and 43 samples ]
## sample_data() Sample Data:       [ 43 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 8 taxa by 8 taxonomic ranks ]

Or to keep only Bacteroidetes phylum. Note that taxa names have been shortened from the rownames of the ExpressionSet object, for nicer plotting.

loman.bd = subset_taxa( loman.pseq, Phylum == "Bacteroidetes")
head( taxa_names( loman.bd ) )
## [1] "p__Bacteroidetes"      "c__Bacteroidia"        "o__Bacteroidales"     
## [4] "f__Bacteroidaceae"     "f__Porphyromonadaceae" "f__Rikenellaceae"

6.4 Advanced Pruning

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 ]

Note that phyloseq also provides topk() for selecting the most abundant k taxa, and other functions for advanced pruning of taxa.

6.5 Taxonomy Heatmap

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 strain-level OTUs. 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.

loman.filt = subset_taxa(loman.pseq, keepotu & !is.na(Strain))
plot_heatmap(loman.filt, method="PCoA", distance="bray")