1 Prerequisites

The following knitr options will be used in this vignette to provide the most useful and concise output.

knitr::opts_chunk$set(message = FALSE)

The following packages will be used in this vignette to provide demonstrative examples of what a user might do with HMP16SData.

library(HMP16SData)
library(phyloseq)
library(magrittr)
library(ggplot2)
library(tibble)
library(dplyr)
library(dendextend)
library(circlize)
library(curatedMetagenomicData)
library(gridExtra)
library(cowplot)
library(readr)
library(haven)

Pipe operators from the magrittr package are used in this vignette to provide the most elegant and concise syntax. See the magrittr vignette if the syntax is unclear.

2 Introduction

HMP16SData is a Bioconductor ExperimentData package of the Human Microbiome Project (HMP) 16S rRNA sequencing data for variable regions 1–3 and 3–5. Raw data files are provided in the package as downloaded from the HMP Data Analysis and Coordination Center. Processed data is provided as SummarizedExperiment class objects via ExperimentHub.

HMP16SData can be installed using BiocInstaller as follows.

BiocInstaller::biocLite("HMP16SData")

Once installed, HMP16SData provides two functions to access data – one for variable region 1–3 and another for variable region 3–5. When called, as follows, the functions will download data from an ExperimentHub Amazon S3 (Simple Storage Service) bucket over https or load data from a local cache.

V13()

## class: SummarizedExperiment 
## dim: 43140 2898 
## metadata(2): experimentData phylogeneticTree
## assays(1): 16SrRNA
## rownames(43140): OTU_97.1 OTU_97.10 ... OTU_97.9997 OTU_97.9999
## rowData names(7): CONSENSUS_LINEAGE SUPERKINGDOM ... FAMILY GENUS
## colnames(2898): 700013549 700014386 ... 700114963 700114965
## colData names(7): RSID VISITNO ... HMP_BODY_SUBSITE SRS_SAMPLE_ID

V35()

## class: SummarizedExperiment 
## dim: 45383 4743 
## metadata(2): experimentData phylogeneticTree
## assays(1): 16SrRNA
## rownames(45383): OTU_97.1 OTU_97.10 ... OTU_97.9998 OTU_97.9999
## rowData names(7): CONSENSUS_LINEAGE SUPERKINGDOM ... FAMILY GENUS
## colnames(4743): 700013549 700014386 ... 700114717 700114750
## colData names(7): RSID VISITNO ... HMP_BODY_SUBSITE SRS_SAMPLE_ID

The two data sets are represented as SummarizedExperiment objects, a standard Bioconductor class that is amenable to subsetting and analysis. To maintain brevity, details of the SummarizedExperiment class are not outlined here but the SummarizedExperiment package provides an excellent vignette.

3 Features

V13() %>%
    table_one() %>%
    head()

##   Visit Number    Sex                                  Run Center
## 1          One Female                  Baylor College of Medicine
## 2          One   Male Baylor College of Medicine, Broad Institute
## 3          One   Male Baylor College of Medicine, Broad Institute
## 4          One   Male Baylor College of Medicine, Broad Institute
## 5          One   Male Baylor College of Medicine, Broad Institute
## 6          One   Male Baylor College of Medicine, Broad Institute
##            HMP Body Site HMP Body Subsite
## 1 Gastrointestinal Tract            Stool
## 2 Gastrointestinal Tract            Stool
## 3                   Oral           Saliva
## 4                   Oral    Tongue Dorsum
## 5                   Oral      Hard Palate
## 6                   Oral    Buccal Mucosa

list(V13 = V13(), V35 = V35()) %>%
    table_one() %>%
    kable_one()

V35_stool <-
    V35() %>%
    subset(select = HMP_BODY_SUBSITE == "Stool")

V35_stool

## class: SummarizedExperiment 
## dim: 45383 319 
## metadata(2): experimentData phylogeneticTree
## assays(1): 16SrRNA
## rownames(45383): OTU_97.1 OTU_97.10 ... OTU_97.9998 OTU_97.9999
## rowData names(7): CONSENSUS_LINEAGE SUPERKINGDOM ... FAMILY GENUS
## colnames(319): 700013549 700014386 ... 700114717 700114750
## colData names(7): RSID VISITNO ... HMP_BODY_SUBSITE SRS_SAMPLE_ID

V35_stool_protected <-
    V35_stool %>%
    attach_dbGaP("~/prj_12146.ngc")

colData(V35_stool_protected)

V13_tonsils <-
    V13() %>%
    subset(select = HMP_BODY_SUBSITE == "Palatine Tonsils")

V13_stool <-
    V13() %>%
    subset(select = HMP_BODY_SUBSITE == "Stool")

V13_tonsils_phyloseq <-
    as_phyloseq(V13_tonsils)

V13_stool_phyloseq <-
    as_phyloseq(V13_stool)

sample_samples <- function(x, size) {
    sampled_names <-
        sample_names(x) %>%
        sample(size)

    prune_samples(sampled_names, x)
}

V13_tonsils_phyloseq %<>%
    sample_samples(25)

V13_stool_phyloseq %<>%
    sample_samples(25)

sample_data(V13_tonsils_phyloseq)$Study <- "Tonsils"

sample_data(V13_stool_phyloseq)$Study <- "Stool"

V13_phyloseq <-
    merge_phyloseq(V13_tonsils_phyloseq, V13_stool_phyloseq)

V13_phyloseq %<>%
    taxa_sums() %>%
    is_greater_than(0) %>%
    prune_taxa(V13_phyloseq)

richness_measures <-
    c("Observed", "Shannon", "Simpson")

V13_phyloseq %>%
    plot_richness(x = "Study", color = "Study", measures = richness_measures) +
    stat_boxplot(geom ="errorbar") +
    geom_boxplot() +
    theme_bw() +
    theme(axis.title.x = element_blank(), legend.position = "none")

V13_dendrogram <-
    distance(V13_phyloseq, method = "bray") %>%
    hclust() %>%
    as.dendrogram()

V13_sample_data <-
    sample_data(V13_phyloseq) %>%
    data.frame()

V13_sample_data %<>%
    rownames_to_column(var = "PSN") %>%
    mutate(labels_col = if_else(Study == "Stool", "#F8766D", "#00BFC4")) %>%
    mutate(leaves_col = if_else(Study == "Stool", "#F8766D", "#00BFC4")) %>%
    mutate(leaves_pch = if_else(Study == "Stool", 16, 17))

V13_sample_order <-
    labels(V13_dendrogram) %>%
    match(V13_sample_data$PSN)

labels_col <- V13_sample_data$labels_col[V13_sample_order]
leaves_col <- V13_sample_data$leaves_col[V13_sample_order]
leaves_pch <- V13_sample_data$leaves_pch[V13_sample_order]

V13_dendrogram %<>%
    set("labels_col", labels_col) %>%
    set("leaves_col", leaves_col) %>%
    set("leaves_pch", leaves_pch)

V13_dendrogram %>%
    circlize_dendrogram(labels_track_height = 0.2)

V13_ordination <-
    ordinate(V13_phyloseq, method = "PCoA", distance = "bray")

V13_phyloseq %>%
    plot_ordination(V13_ordination, color = "Study", shape = "Study") +
    theme_bw() +
    theme(legend.position = "bottom")

V13_ordination %>%
    plot_scree() +
    theme_bw()

4 Phylum-level Comparison to Metagenomic Shotgun Sequencing

In addition to 16S rRNA sequencing, the HMP Study conducted whole metagenome shotgun (MGX) sequencing. These profiles, along with thousands of profiles from other studies, are available in the curatedMetagenomicData package. Here, a phylum-level relative abundance comparison of the 16S and MGX samples is made to illustrate comparing these data sets.

First a V35_stool_phyloseq object is constructed and contains the 16S variable region 3–5 data for stool samples. Then the MGX stool samples are obtained and coerced to a phyloseq object as follows.

V35_stool_phyloseq <-
    V35_stool %>%
    as_phyloseq()

MGX_stool_phyloseq <-
    HMP_2012.metaphlan_bugs_list.stool(cmdversion = "20170526") %>%
    ExpressionSet2phyloseq()

The curatedMetagenomicData package provides taxonomic relative abundance for MGX data (with an option to estimate counts by multiplying by read depth) from MetaPhlAn2. MetaPhlAn2 directly estimates relative abundance at every taxonomic level based on clade-specific markers, and these estimates are better than summing lower-level taxonomic abundances.

Note, this comparison becomes complicated because phyloseq does not currently support row and column matching and reordering. Instead, we use the phyloseq::psmelt() method to generate data.frame objects for further manipulation.

MGX_stool_melted <-
    MGX_stool_phyloseq %>%
    subset_taxa(is.na(Class) & !is.na(Phylum)) %>%
    psmelt()

The 16S data sets do not contain summarized counts for higher taxonomic levels, so we use the phyloseq::tax_glom() method to agglomerate at phylum level.

V35_stool_melted <-
    V35_stool_phyloseq %>%
    tax_glom(taxrank = "PHYLUM") %>%
    psmelt()

There is a column SRS_SAMPLE_ID present in the 16S variable region 3–5 data that is explicitly for mapping to the MGX samples and the matching identifiers are in the NCBI_accession column of the MGX sample data. The intersection of the two vectors represents the matching samples that are present in both data sets. This intersection, shown below as SRS_SAMPLE_IDS, can then be used to filter both the 16S and MGX samples to include only the samples in common.

Along with the filter step shown below, standardization of sample identifiers to the SRS numbers and conversion of taxonomic counts to relative abundance is done. The conversion to relative abundance is necessary for comparability across the 16S and MGX samples.

In either case, data is first grouped by samples and the percent composition of each phylum relative to the others in the sample is calculated by taking the count abundance of the phylum divided by the sum of all abundance counts in the sample. Samples are then sorted by phylum and descending relative abundance before being grouped by phylum – the grouping by phylum allows for the assignment a phylum rank that is then used to sort the phyla by descending relative abundance. Finally, only the sample, phylum, and relative abundance columns are kept once the order has been established and all groupings can be removed.

SRS_SAMPLE_IDS <-
    intersect(V35_stool_melted$SRS_SAMPLE_ID, MGX_stool_melted$NCBI_accession)

V35_stool_melted %<>%
    filter(SRS_SAMPLE_ID %in% SRS_SAMPLE_IDS) %>%
    rename(Phylum = PHYLUM) %>%
    mutate(Sample = SRS_SAMPLE_ID) %>%
    group_by(Sample) %>%
    mutate(`Relative Abundance` = Abundance / sum(Abundance)) %>%
    arrange(Phylum, -`Relative Abundance`) %>%
    group_by(Phylum) %>%
    mutate(phylum_rank = sum(`Relative Abundance`)) %>%
    arrange(-phylum_rank) %>%
    select(Sample, Phylum, `Relative Abundance`) %>%
    ungroup()

MGX_stool_melted %<>%
    filter(NCBI_accession %in% SRS_SAMPLE_IDS) %>%
    group_by(Sample) %>%
    mutate(`Relative Abundance` = Abundance / sum(Abundance)) %>%
    arrange(Phylum, -`Relative Abundance`) %>%
    group_by(Phylum) %>%
    mutate(phylum_rank = sum(`Relative Abundance`)) %>%
    arrange(-phylum_rank) %>%
    select(Sample, Phylum, `Relative Abundance`) %>%
    ungroup()

Provided that the phyla are ordered by relative abundance, the top phyla from each data set can be obtained and intersected to give the top eight phyla present in both data sets. The top eight phyla can then be used to filter 16S and MGX samples to include only the desired phyla.

V35_top_phyla <-
    V35_stool_melted %$%
    unique(Phylum) %>%
    as.character()

MGX_top_phyla <-
    MGX_stool_melted %$%
    unique(Phylum) %>%
    as.character()

top_eight_phyla <-
    intersect(V35_top_phyla, MGX_top_phyla) %>%
    extract(1:8)

V35_stool_melted %<>%
    filter(Phylum %in% top_eight_phyla)

MGX_stool_melted %<>%
    filter(Phylum %in% top_eight_phyla)

To achieve ordering of samples by the relative abundance of the top phylum when plotting, a vector of the unique sample identifiers is constructed and will be used as the levels of a factor.

sample_order <-
    V35_stool_melted %$%
    unique(Sample)

A vector of unique phyla is also constructed and will be used as the levels of a factor when plotting. If this were not done, phyla would simply be sorted alphabetically.

phylum_order <-
    V35_stool_melted %$%
    unique(Phylum)

Color blindness affects a significant portion of the population, but, using an intelligent color pallet, figures can be designed to be friendly to those with deuteranopia, protanopia, and tritanopia. The following colors are derived from Wong, B. Points of view: Color blindness. Nat. Methods 8, 441 (2011).

bang_wong_colors <-
    c("#CC79A7", "#D55E00", "#0072B2", "#F0E442", "#009E73", "#56B4E9",
      "#E69F00", "#000000")

Using the sample_order and phylum_order vectors constructed above, stacked phylum-level relative abundance bar plots sorted by decreasing relative abundance of the top phylum and stratified by the top eight phyla can be made for 16S and MGX samples. The two plots are made separately and serialized so that they can be arranged in a single figure for comparison.

V35_plot <-
    V35_stool_melted %>%
    mutate(Sample = factor(Sample, sample_order)) %>%
    mutate(Phylum = factor(Phylum, phylum_order)) %>%
    ggplot(aes(Sample, `Relative Abundance`, fill = Phylum)) +
    geom_bar(stat = "identity", position = "fill", width = 1) +
    scale_fill_manual(values = bang_wong_colors) +
    theme_minimal() +
    theme(axis.text.x = element_blank(), axis.title.x = element_blank(),
          panel.grid = element_blank(), legend.position = "none",
          legend.title = element_blank()) +
    ggtitle("Phylum-Level Relative Abundance", "16S Stool Samples")

MGX_plot <-
    MGX_stool_melted %>%
    mutate(Sample = factor(Sample, sample_order)) %>%
    mutate(Phylum = factor(Phylum, phylum_order)) %>%
    ggplot(aes(Sample, `Relative Abundance`, fill = Phylum)) +
    geom_bar(stat = "identity", position = "fill", width = 1) +
    scale_fill_manual(values = bang_wong_colors) +
    theme_minimal() +
    theme(axis.text.x = element_blank(), axis.title.x = element_blank(),
          panel.grid = element_blank(), legend.position = "none",
          legend.title = element_blank()) +
    ggtitle("", "MGX Stool Samples")

In the plots above, legends are removed to reduce redundancy, but a legend is still necessary and can be serialized as its own plot using the get_legend method from the cowplot package.

plot_legend <- {
        MGX_plot +
            theme(legend.position = "bottom")
    } %>%
    get_legend()

Finally, the grid.arrange method from the gridExtra package is used to arrange, scale, and plot the three plots in a single figure.

grid.arrange(V35_plot, MGX_plot, plot_legend, ncol = 1, heights = c(3, 3, 1))

Notably, the figure illustrates the Bacteroidetes/Firmicutes gradient with reasonable agreement between the 16S and MGX samples.

5 Exporting Data to CSV, SAS, SPSS, and STATA Formats

To our knowledge, R and Bioconductor provide the most and best methods for the analysis of microbiome data. However, we realize they are not the only analysis environments and wish to provide methods to export the data from HMP16SData to CSV, SAS, SPSS, and STATA formats. As we do not use these other languages regularly, we are unaware of how to represent phylogenitic trees in them and will not attempt to export the trees here.

V13_participant_data <-
  V13() %>%
  colData() %>%
  as.data.frame() %>%
  rownames_to_column(var = "PSN")

V13_OTU_counts <-
  V13() %>%
  assay() %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column(var = "PSN")

V13_data <-
  merge.data.frame(V13_participant_data, V13_OTU_counts, by = "PSN")

V13_dictionary <-
  V13() %>%
  rowData() %>%
  t.data.frame() %>%
  as.data.frame()

colnames(V13_data) <-
    colnames(V13_data) %>%
    gsub(pattern = "\\.", replacement = "_", x = .)

colnames(V13_dictionary) <-
    colnames(V13_dictionary) %>%
    gsub(pattern = "\\.", replacement = "_", x = .)

write_csv(V13_data, "~/V13_data.csv")
write_csv(V13_dictionary, "~/V13_dictionary.csv")

write_sas(V13_data, "~/V13_data.sas7bdat")
write_sas(V13_dictionary, "~/V13_dictionary.sas7bdat")

write_sav(V13_data, "~/V13_data.sav")
write_sav(V13_dictionary, "~/V13_dictionary.sav")

write_dta(V13_data, "~/V13_data.dta", version = 13)
write_dta(V13_dictionary, "~/V13_dictionary.dta", version = 13)

6 Session Information

sessionInfo()

## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.7-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] haven_1.1.2                   readr_1.1.1                  
##  [3] cowplot_0.9.3                 gridExtra_2.3                
##  [5] curatedMetagenomicData_1.10.2 bindrcpp_0.2.2               
##  [7] ExperimentHub_1.6.0           AnnotationHub_2.12.0         
##  [9] circlize_0.4.4                dendextend_1.8.0             
## [11] dplyr_0.7.6                   tibble_1.4.2                 
## [13] ggplot2_3.0.0                 magrittr_1.5                 
## [15] phyloseq_1.24.2               HMP16SData_1.0.1             
## [17] SummarizedExperiment_1.10.1   DelayedArray_0.6.5           
## [19] BiocParallel_1.14.2           matrixStats_0.54.0           
## [21] Biobase_2.40.0                GenomicRanges_1.32.6         
## [23] GenomeInfoDb_1.16.0           IRanges_2.14.11              
## [25] S4Vectors_0.18.3              BiocGenerics_0.26.0          
## [27] BiocStyle_2.8.2              
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.2               plyr_1.8.4                   
##   [3] igraph_1.2.2                  lazyeval_0.2.1               
##   [5] splines_3.5.1                 digest_0.6.16                
##   [7] foreach_1.4.4                 BiocInstaller_1.30.0         
##   [9] htmltools_0.3.6               viridis_0.5.1                
##  [11] memoise_1.1.0                 cluster_2.0.7-1              
##  [13] Biostrings_2.48.0             colorspace_1.3-2             
##  [15] blob_1.1.1                    rvest_0.3.2                  
##  [17] xfun_0.3                      crayon_1.3.4                 
##  [19] RCurl_1.95-4.11               jsonlite_1.5                 
##  [21] bindr_0.1.1                   survival_2.42-6              
##  [23] iterators_1.0.10              ape_5.1                      
##  [25] glue_1.3.0                    kableExtra_0.9.0             
##  [27] gtable_0.2.0                  zlibbioc_1.26.0              
##  [29] XVector_0.20.0                kernlab_0.9-27               
##  [31] Rhdf5lib_1.2.1                shape_1.4.4                  
##  [33] prabclus_2.2-6                DEoptimR_1.0-8               
##  [35] scales_1.0.0                  mvtnorm_1.0-8                
##  [37] DBI_1.0.0                     Rcpp_0.12.18                 
##  [39] viridisLite_0.3.0             xtable_1.8-2                 
##  [41] bit_1.1-14                    mclust_5.4.1                 
##  [43] httr_1.3.1                    fpc_2.1-11.1                 
##  [45] modeltools_0.2-22             pkgconfig_2.0.2              
##  [47] flexmix_2.3-14                nnet_7.3-12                  
##  [49] tidyselect_0.2.4              labeling_0.3                 
##  [51] rlang_0.2.2                   reshape2_1.4.3               
##  [53] later_0.7.3                   AnnotationDbi_1.42.1         
##  [55] munsell_0.5.0                 tools_3.5.1                  
##  [57] RSQLite_2.1.1                 ade4_1.7-11                  
##  [59] evaluate_0.11                 biomformat_1.8.0             
##  [61] stringr_1.3.1                 yaml_2.2.0                   
##  [63] knitr_1.20                    bit64_0.9-7                  
##  [65] robustbase_0.93-2             purrr_0.2.5                  
##  [67] nlme_3.1-137                  whisker_0.3-2                
##  [69] mime_0.5                      xml2_1.2.0                   
##  [71] compiler_3.5.1                rstudioapi_0.7               
##  [73] curl_3.2                      interactiveDisplayBase_1.18.0
##  [75] stringi_1.2.4                 highr_0.7                    
##  [77] forcats_0.3.0                 lattice_0.20-35              
##  [79] trimcluster_0.1-2.1           Matrix_1.2-14                
##  [81] vegan_2.5-2                   permute_0.9-4                
##  [83] multtest_2.36.0               pillar_1.3.0                 
##  [85] GlobalOptions_0.1.0           data.table_1.11.4            
##  [87] bitops_1.0-6                  httpuv_1.4.5                 
##  [89] R6_2.2.2                      bookdown_0.7                 
##  [91] promises_1.0.1                codetools_0.2-15             
##  [93] MASS_7.3-50                   assertthat_0.2.0             
##  [95] rhdf5_2.24.0                  rprojroot_1.3-2              
##  [97] withr_2.1.2                   GenomeInfoDbData_1.1.0       
##  [99] diptest_0.75-7                mgcv_1.8-24                  
## [101] hms_0.4.2                     grid_3.5.1                   
## [103] tidyr_0.8.1                   class_7.3-14                 
## [105] rmarkdown_1.10                shiny_1.1.0