Instructors names and contact information

1. Levi Waldron, Graduate School of Public Health and Health Policy, City University of New York, New York, NY. email: Levi.Waldron@sph.cuny.edu

2. Ni Zhao, Department of Biostatistics, Johns Hopkins University, Baltimore, MD. email: nzhao10@jhu.edu

3. Ekaterina Smirnova, Department of Biostatistics, Virginia Commonwealth University, Richmond, VA. email: ekaterina.smirnova@vcuhealth.org

Workshop Description

The composition of microbial species in a human body is essential for maintaining human health, and it is associated with a number of diseases including obesity, bowel inflammatory disease, and bacterial vaginosis. Over the last decade, microbiome data analysis almost entirely shifted towards using samples taken directly from various sites of human body and to explore a large number of microbes using 16S or whole metagenome sequencing. This technological shift has led to a radical change in the data collected and led to the creation of the Human Microbiome Project (HMP) in 2008. With the growth and success of the microbiomics field, the size and complexity, and availability of the microbiome data in any given experiment have increased exponentially. Publicly-available data sets from amplicon sequencing of two 16S ribosomal RNA variable regions, with extensive controlled-access participant data, provide a reference for ongoing microbiome studies. However, utilization of these data sets can be hindered by the complex bioinformatic steps required to access, import, decrypt, and merge the various components in formats suitable for ecological and statistical analysis.

The main goal of this workshop is to provide a much-needed tutorial on downloading, understanding, and analyzing the publicly-available HMP data. We will describe the two packages, HMP16SData and HMP2Data, that provide the comprehensive tutorial on the analysis of Phase 1 and 2 of the HMP project, respectively. These packages provide count data for both 16S ribosomal RNA variable regions, integrated with phylogeny, taxonomy, public participant data, and controlled participant data for authorized researchers, using standard integrative Bioconductor data objects. As such, these packages provide the first comprehensive publicly-available tool for working with microbiome data collected at different body sites (e.g., fecal, nasal, vaginal) coupled with other omics approaches (e.g., cytokines, whole metagenome), and is conducted longitudinally on multiple subjects at multiple visits. By removing bioinformatic hurdles of data access and management, HMP16SData and HMP2Data enable researchers with only basic R skills to quickly analyze HMP data. This allows studying the dynamics of microbial composition over the disease progression, discovering disease-specific microbial biomarkers, and understanding the role of microbiome in connection to other omics measurements.

This workshop will consist of the series on the instructors-led live demos presented together with the brief lecture materials and statistical methods review when necessary. Completely worked out examples are available through the HMP16SData and HMP2Data package vignettes.

Pre-requisites

• Basic knowledge of R syntax
• Familiarity with the microbiome data
• Familiarity with alpha and beta diversity in the context of microbiome data analysis
• Familiarity with principal components and principal coordinates analysis

Background reading:

• Schiffer L, Azhar R, Shepherd L, Ramos M, Geistlinger L, Huttenhower C, Dowd JB, Segata N, Waldron L (2019). “HMP16SData: Efficient Access to the Human Microbiome Project through Bioconductor.” American Journal of Epidemiology. doi: 10.1093/aje/kwz006.

Time outline

Workshop goals and objectives

The main purpose of this workshop is to help researchers interested in public microbiome data to start using the HMP16SData and HMP2Data packages.

Learning goals

• understand the goals of HMP projects and the studies available through the HMP DAC portal
• download and use HMP data packages
• understand the difference between phase 1 and phase 2 projects
• understand principles of -omics data integration and visualization

Learning objectives

• install HMP data packages and access vignettes
• download HMP data to produce analyses illustrated in package vignettes
• create sample summary statistics
• create alpha and beta diversity plots
• understand longitudinal patterns and examine microbiome diversity changes
  
• integrate 16S and cytokines data using co-inertia analyses techniques
• learn capturing patterns across -omics data sets of different structure

Integrative Human Microbiome Project (iHMP) https://hmpdacc.org/ihmp/

• One of the largest open data resources for studying longitudinal microbiome changes and connection to disease

• Novel data – first scientific results just published in the special collection of Nature and its sister journals (https://www.nature.com/articles/s41586-019-1238-8)

• Human Microbiome Project Data Portal (https://portal.hmpdacc.org/)
• Open data – but not straightforward to download and work with

• Need Aspera client to even download the data – not every lab/researcher has these expertise available
• Still need multiple processing steps after data is downloaded to:

1. map taxonomy ids from the database

2. merge with meta data – not available on the DAC portal

3. construct phylogenetic tree

4. merge omics data.

HMP2Data package

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

suppressPackageStartupMessages({
  library(phyloseq)
  library(SummarizedExperiment)
  library(MultiAssayExperiment)
  library(dplyr)
  library(ggplot2)
  library(UpSetR)
  library(ade4)
  library(vegan)
  library(MiRKAT)
})

HMP2Data package is currently under review as a Bioconductor package. Current (development) version can be installed can be installed from John Stansfield’s GitHub repository https://github.com/jstansfield0/HMP2Data as follows. This is the version used for the workshop.

BiocManager::install("jstansfield0/HMP2Data")

The stable version is available at the dozmorovlab GitHub repository https://github.com/dozmorovlab/HMP2Data

BiocManager::install("dozmorovlab/HMP2Data")

Load HMP2Data package

library(HMP2Data)

Once installed, HMP2Data provides functions to access the data from the three major HMP2 studies.

data("momspi16S_mtx")

data("momspi16S_tax")
# Check if Greengene IDs match between the 16S and taxonomy data
# all.equal(rownames(momspi16S_mtx), rownames(momspi16S_tax)) # Should be TRUE

data("momspi16S_samp")

# Check if sample names match between the 16S and sample data
# all.equal(colnames(momspi16S_mtx), rownames(momspi16S_samp)) # Should be TRUE

momspi16S_phyloseq <- momspi16S()
momspi16S_phyloseq
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 7665 taxa and 9107 samples ]
#> sample_data() Sample Data:       [ 9107 samples by 13 sample variables ]
#> tax_table()   Taxonomy Table:    [ 7665 taxa by 7 taxonomic ranks ]

data("momspiCyto_mtx")
momspiCyto_mtx[1:5, 1:3]
#>         EP004835_K10_MVAX EP004835_K20_MVAX EP004835_K30_MVAX
#> Eotaxin             68.05             60.08             37.84
#> FGF                265.45             -2.00             -2.00
#> G-CSF              624.79           2470.70             -2.00
#> GM-CSF              51.39             -2.00             42.01
#> IFN-g              164.58            193.22            184.53

data("momspiCyto_samp")
colnames(momspiCyto_samp)
#>  [1] "file_id"          "md5"              "size"            
#>  [4] "urls"             "sample_id"        "file_name"       
#>  [7] "subject_id"       "sample_body_site" "visit_number"    
#> [10] "subject_gender"   "subject_race"     "study_full_name" 
#> [13] "project_name"

momspiCyto <- momspiCytokines()
momspiCyto
#> class: SummarizedExperiment 
#> dim: 29 1396 
#> metadata(0):
#> assays(1): ''
#> rownames(29): Eotaxin FGF ... FGF basic IL-17
#> rowData names(1): cytokine
#> colnames(1396): EP004835_K10_MVAX EP004835_K20_MVAX ...
#>   EP996091_K40_MVAX EP996091_K60_MVAX
#> colData names(13): file_id md5 ... study_full_name project_name

data("IBD16S_mtx")

data("IBD16S_tax")
colnames(IBD16S_tax)
#> [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"

data("IBD16S_samp")
colnames(IBD16S_samp) %>% head()
#> [1] "Project"       "sample_id"     "subject_id"    "site_sub_coll"
#> [5] "data_type"     "week_num"

IBD <- IBD16S()
IBD
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 982 taxa and 178 samples ]
#> sample_data() Sample Data:       [ 178 samples by 490 sample variables ]
#> tax_table()   Taxonomy Table:    [ 982 taxa by 6 taxonomic ranks ]

data("T2D16S_mtx")

data("T2D16S_tax")
colnames(T2D16S_tax)
#> [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"

data("T2D16S_samp")
colnames(T2D16S_samp)
#>  [1] "file_id"          "md5"              "size"            
#>  [4] "urls"             "sample_id"        "file_name"       
#>  [7] "subject_id"       "sample_body_site" "visit_number"    
#> [10] "subject_gender"   "subject_race"     "study_full_name" 
#> [13] "project_name"

T2D <- T2D16S()
T2D
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 12062 taxa and 2208 samples ]
#> sample_data() Sample Data:       [ 2208 samples by 13 sample variables ]
#> tax_table()   Taxonomy Table:    [ 12062 taxa by 7 taxonomic ranks ]

Features

list("MOMS-PI 16S" = momspi16S_phyloseq, "MOMS-PI Cytokines" = momspiCyto, "IBD 16S" = IBD, "T2D 16S" = T2D) %>% table_two()

# make data.frame for plotting
plot_visits <- data.frame(study = c(rep("MOMS-PI Cytokines", nrow(momspiCyto_samp)),
                     rep("IBD 16S", nrow(IBD16S_samp)),
                     rep("T2D 16S", nrow(T2D16S_samp))),
          visits = c(momspiCyto_samp$visit_number,
                     IBD16S_samp$visit_number,
                     T2D16S_samp$visit_number))
p2 <- ggplot(plot_visits, aes(x = visits, fill = study)) + 
  geom_histogram(position = "dodge", alpha = 0.7, bins = 30, color = "#00BFC4") + xlim(c(0, 40)) +
  theme(legend.position = c(0.8, 0.8))  + 
  xlab("Visit number") + ylab("Count")
p2
#> Warning: Removed 677 rows containing non-finite values (stat_bin).
#> Warning: Removed 4 rows containing missing values (geom_bar).

UpsetR plots

# set up data.frame for UpSetR
momspi_upset <- aggregate(momspi16S_samp$sample_body_site, by = list(momspi16S_samp$subject_id), table)
tmp <- as.matrix(momspi_upset[, -1])
tmp <- (tmp > 0) *1
momspi_upset <- data.frame(patient = momspi_upset$Group.1, tmp)

# plot
upset(momspi_upset, order.by = 'freq', matrix.color = "blue", text.scale = 2)

# set up data.frame for UpSetR
momspiCyto_upset <- aggregate(momspiCyto_samp$sample_body_site, by = list(momspiCyto_samp$subject_id), table)
tmp <- as.matrix(momspiCyto_upset[, -1])
tmp <- (tmp > 0) *1
momspiCyto_upset <- data.frame(patient = momspiCyto_upset$Group.1, tmp)

# plot
upset(momspiCyto_upset, order.by = 'freq', matrix.color = "blue", text.scale = 2)

# set up data.frame for UpSetR
T2D_upset <- aggregate(T2D16S_samp$sample_body_site, by = list(T2D16S_samp$subject_id), table)
tmp <- as.matrix(T2D_upset[, -1])
tmp <- (tmp > 0) *1
T2D_upset <- data.frame(patient = T2D_upset$Group.1, tmp)

# plot
upset(T2D_upset, order.by = 'freq', matrix.color = "blue", text.scale = 2)

Alpha Diversity Definitions

Alpha diversity is a within-sample diversity measurement that models the richness and/or eveness of the microbial community. Different alpha diversities emphasize different aspect of the community. Commonly used alpha diversities include:

Of course, the alpha diversity measurement is impacted by the sequencing depth (total number of reads per sample). Rarefying (through vegan:rarefy) is necessary before calculating alpha diversity.

Comparing Alpha Diversities between Groups

We can use “plot_richness” from phyloseq package to plot the alpha diversities.

IBD %<>%
  taxa_sums() %>%
  is_greater_than(0) %>%
  prune_taxa(IBD)
sample_data(IBD)$diagnosis = factor(sample_data(IBD)$diagnosis)

p = plot_richness(IBD, x="diagnosis", color="diagnosis", measures=c("Observed", "InvSimpson", "Shannon", "Chao1")) + geom_jitter()
p + geom_boxplot(data = p$data, aes(x = diagnosis, y = value, color = NULL), alpha = 0.1)
#> Warning: Removed 534 rows containing missing values (geom_errorbar).

Once we summarize the whole community through a one-dimensional summary statistics (alpha diversity), all the traditional methods for hypothesis testing are applicable.

alpha = estimate_richness(IBD, measures=c("Observed", "InvSimpson", "Shannon"))
alpha_df = data.frame(sample_data(IBD)[,c("sample_id","subject_id")], alpha)
rownames(alpha_df) = NULL
kable(head(alpha_df),
      digits = 3,caption = "Some Alpha Diversities of the IBD dataset") %>% kable_styling(bootstrap_options = c("striped", "hover"), full_width = F,position = "center") 

Beta diversity analysis

par(mfrow = c(1,2))
jac <- as.matrix(distance(t(otu_table(IBD)), method = "jaccard"))
mod1 <- vegan::betadisper(as.dist(jac), as.factor(IBD16S_samp$diagnosis))
plot(mod1, ellipse = TRUE, hull = FALSE, main="Jaccard Distance", xlab="PC 1", ylab="PC 2", sub=NULL)
box()

bc = as.matrix(ecodist::bcdist(t(otu_table(IBD))))
mod2 <- vegan::betadisper(as.dist(bc), as.factor(IBD16S_samp$diagnosis))
plot(mod2, ellipse = TRUE, hull = FALSE, main="Bray-Curtis Dissimilarity", xlab="PC 1", ylab="PC 2", sub=NULL, xlim = c(-0.2, 0.2))
legend("bottomright", legend=c("CD", "nonIBD", "UC"), pch = 1:3, col=c("Black", "red", "Green"), cex=1.3, box.lty=0)
box()

IBD0 = subset_samples(IBD,  week_num == 0&diagnosis %in% c("nonIBD", "CD"))
IBD0 %<>%
  taxa_sums() %>%
  is_greater_than(0) %>%
  prune_taxa(IBD0)

jac0 <- as.matrix(distance(t(otu_table(IBD0)), method = "jaccard"))
bc0 =  as.matrix(ecodist::bcdist(t(otu_table(IBD0))))
CD = as.numeric(sample_data(IBD0)$diagnosis == "CD")
gender = as.numeric(factor(sample_data(IBD0)$subject_gender))

fit = vegan::adonis(bc0 ~ gender + CD, permutations=500)
kable(fit$aov.tab,
      digits = 3,caption = "PERMANOVA result using Bray-Curtis distance") %>% kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)  

 K_bc = MiRKAT::D2K(bc0); K_jac = D2K(jac0)
 y = ifelse(sample_data(IBD0)$diagnosis == "CD",1, 0 )
 out = MiRKAT::MiRKAT(y=y, X = gender, Ks = list(K_bc, K_jac), out_type = "D") 

Co-inertia analysis

Basics:

1. Let X and Y be 16S and cytokines tables respectively
2. Rows: same n women at first visit
3. Columns: p₁ taxa and p₂ cytokines
4. Column weights: Q_X and Q_Y
5. Row weights: D
6. PCA analysis of each table: (X, Q_X, D) and (Y, Q_Y, D)
7. Co-inertia axes: Y^TDX = KΛ^1/2A^T from eigendecomposition of (Y^TDX, Q_X, Q_Y)
8. Plot F_X = XA and F_Y = YK

Taxa normalization:

1. Center 16S data to work with PCA on the covariance matrix Σ_X = Cov(X)
2. To normalize the magnitude, divide each value of X by the total variance: $\sqrt{\mbox{tr}(\Sigma_X)}$
3. Note: step 2 is quivalent to dividing the matrix by $\sqrt{\sum_{k=1}^r \lambda_k}$, where λ_k are the eigevalues of Σ_X and r is the rank of X.

taxa_mtx <- scale(combined_16S_mtx, center = TRUE, scale = FALSE)
#use fast trace computation formula: tr(A^B) = sum(A*B), where '*' operator refers to elemetwise product
taxa_tr <- sum(taxa_mtx*taxa_mtx)/(dim(taxa_mtx)[1]-1)
taxa_mtx <- taxa_mtx/sqrt(taxa_tr)
taxa.pca <- ade4::dudi.pca(taxa_mtx, scannf=FALSE, nf =61,
                     center = FALSE, scale = FALSE)

Cytokines normalization:

1. Center and scale cytokines data to work with PCA on the correlation matrix Σ_X = Cor(Y)
2. To normalize the magnitude, divide each value of Y by the total variance: $\sqrt{\mbox{tr}(\Sigma_Y)}$

cyto_mtx <- scale(combined_Cyto_mtx, center = TRUE, scale = TRUE)
cyto_tr <- sum(cyto_mtx*cyto_mtx)/(dim(cyto_mtx)[1]-1)
cyto_mtx <- cyto_mtx/sqrt(cyto_tr)
cyto.pca <- ade4::dudi.pca(cyto_mtx, scannf=FALSE, nf =61,
                     center = FALSE, scale = FALSE)

Combine the tables using co-inertia

Co-inertia is available through R package ade4. It takes ade4 PCA objects and performes joint eigendecomposition.

coin <- ade4::coinertia(taxa.pca, cyto.pca, scannf = FALSE, nf = 2)

RV coefficient – measure of similarity between 16S and cytokines tables

RV <- coin$RV
RV
#> [1] 0.0396224

Visualization of output

Interpretation of the plots is similar to examining PCA results.

Plot below provides variables projection on common co-inertia axes.

1. Variables projected close to each other reflect similarilty within and across two data sets

2. Variables with larger values on each component (x- and y-axes) have more importance

We start by identifying important taxa (larger component values):

1. visually from the plot, or

taxa.inx <- c(6,  19,  24,  28,  31,  68,  69, 396)

2. by specifying x and y coordinates to pull all variables with coordinates larger (or smaller) than these values. Table 1 (taxa) values can be accessed $co object of coinertia output named coin.

kable(head(coin$co),
      digits = 5)%>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

Similarly, table 2 (here cytokines) variables loadings can be extracted from the $li object of coinertia output named coin.

kable(head(coin$li),
      digits = 5)%>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)

Sample scores plots. Length of the arrows indicates the samples that have larger differences across two data sets.

#> Warning: Ignoring unknown aesthetics: na.rm

#> Warning: Ignoring unknown aesthetics: na.rm

Samples with largest difference across two data sets. Samples with arrow lengths in 0.9 quantile are chosen.

#Taxa with major differences across two sets
large.difs <- rownames(Samp.coin$Dissimilarity[Samp.coin$Dissimilarity$Quantile >= 0.9, ])

Samples with largest differences across two data sets.

small.difs <- rownames(Samp.coin$Dissimilarity[Samp.coin$Dissimilarity$Quantile < 0.1, ])

Examine these taxa and cytokines