curatedMetagenomicData

Lucas Schiffer, MPH1* and Levi Waldron, PhD2**

1Section of Computational Biomedicine, Boston University School of Medicine, Boston, MA, U.S.A.
2Graduate School of Public Health and Health Policy, City University of New York, New York, NY, U.S.A.

*schifferl@bu.edu
**levi.waldron@sph.cuny.edu

Abstract

The curatedMetagenomicData package provides standardized, curated human microbiome data for novel analyses. It includes gene families, marker abundance, marker presence, pathway abundance, pathway coverage, and relative abundance for samples collected from different body sites. The bacterial, fungal, and archaeal taxonomic abundances for each sample were calculated with MetaPhlAn3, and metabolic functional potential was calculated with HUMAnN3. The manually curated sample metadata and standardized metagenomic data are available as (Tree)SummarizedExperiment objects.

Package

curatedMetagenomicData 3.13.0

Contents

1 What curatedMetagenomicData provides

curatedMetagenomicData provides processed data from whole-metagenome shotgun metagenomics, with manually-curated metadata, as integrated and documented Bioconductor TreeSummarizedExperiment objects or TSV flat text exports. It provides 6 types of data for each dataset:

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

Types 1-3 are generated by MetaPhlAn3; 4-6 are generated by HUMAnN3 using the UniRef90 database.

Currently there are:

• 22588 samples from 93 datasets; see the list of available studies
• 141 fields of specimen metadata from original papers, supplementary files, and websites, with manual curation and automated syntax-checking to standardize annotations

2 Additional documentation and resources

See:

1. Available Studies
2. Our Pipeline
3. Changes in cMD 3
4. Reference for cMD Functions
5. The command-line tool
6. Example analyses in R and Python, Docker image, free Cloud instance for teaching/learning

3 Installation

To install curatedMetagenomicData from Bioconductor, use BiocManager as follows.

BiocManager::install("curatedMetagenomicData")

To install curatedMetagenomicData from GitHub, use BiocManager as follows.

BiocManager::install("waldronlab/curatedMetagenomicData", dependencies = TRUE, build_vignettes = TRUE)

Most users should simply install curatedMetagenomicData from Bioconductor.

4 R Packages

To demonstrate the functionality of curatedMetagenomicData, the dplyr and DT packages are needed.

library(dplyr)
library(DT)

5 Sample Metadata

The curatedMetagenomicData package contains a data.frame, sampleMetadata, of manually curated sample metadata to help users understand the nature of studies and samples available prior to returning resources. Beyond this, it serves two purposes: 1) to define study_name, which is used with the curatedMetagenomicData() function to query and return resources, and 2) to define sample_id, which is used with the returnSamples() function to return samples across studies.

To demonstrate, the first ten rows and columns (without any NA values) of sampleMetadata for the AsnicarF_2017 study are shown in the table below.

sampleMetadata |>
    filter(study_name == "AsnicarF_2017") |>
    select(where(~ !any(is.na(.x)))) |>
    slice(1:10) |>
    select(1:10) |>
    datatable(options = list(dom = "t"), extensions = "Responsive")

6 Data Access

There are three main ways to access data resources in curatedMetagenomicData.

1. The curatedMetagenomicData() function to search for and return resources.
2. The returnSamples() function to return samples across studies.
3. Through the curatedMetagenomicDataTerminal command-line interface.

6.1 curatedMetagenomicData()

To access curated metagenomic data, users will use the curatedMetagenomicData() function both to query and return resources. The first argument pattern is a regular expression pattern to look for in the titles of resources available in curatedMetagenomicData; "" will return all resources. The title of each resource is a three part string with “.” as a delimiter – the fields are runDate, studyName, and dataType. The runDate is the date we created the resource and can mostly be ignored by users because if there is more than one date corresponding to a resource, the most recent one is selected automatically – it would be used if a specific runDate was needed.

Multiple resources can be queried or returned with a single call to curatedMetagenomicData(), but only the titles of resources are returned by default.

curatedMetagenomicData("AsnicarF_20.+")
## 2021-03-31.AsnicarF_2017.gene_families
## 2021-03-31.AsnicarF_2017.marker_abundance
## 2021-03-31.AsnicarF_2017.marker_presence
## 2021-03-31.AsnicarF_2017.pathway_abundance
## 2021-03-31.AsnicarF_2017.pathway_coverage
## 2021-03-31.AsnicarF_2017.relative_abundance
## 2021-10-14.AsnicarF_2017.gene_families
## 2021-10-14.AsnicarF_2017.marker_abundance
## 2021-10-14.AsnicarF_2017.marker_presence
## 2021-10-14.AsnicarF_2017.pathway_abundance
## 2021-10-14.AsnicarF_2017.pathway_coverage
## 2021-10-14.AsnicarF_2017.relative_abundance
## 2021-03-31.AsnicarF_2021.gene_families
## 2021-03-31.AsnicarF_2021.marker_abundance
## 2021-03-31.AsnicarF_2021.marker_presence
## 2021-03-31.AsnicarF_2021.pathway_abundance
## 2021-03-31.AsnicarF_2021.pathway_coverage
## 2021-03-31.AsnicarF_2021.relative_abundance

When the dryrun argument is set to FALSE, a list of SummarizedExperiment and/or TreeSummarizedExperiment objects is returned. The rownames argument determines the type of rownames to use for relative_abundance resources: either "long" (the default), "short" (species name), or "NCBI" (NCBI Taxonomy ID). When a single resource is requested, a single element list is returned.

curatedMetagenomicData("AsnicarF_2017.relative_abundance", dryrun = FALSE, rownames = "short")
## $`2021-10-14.AsnicarF_2017.relative_abundance`
## class: TreeSummarizedExperiment 
## dim: 296 24 
## metadata(1): agglomerated_by_rank
## assays(1): relative_abundance
## rownames(296): Escherichia coli Bifidobacterium bifidum ...
##   Streptococcus gordonii Abiotrophia sp. HMSC24B09
## rowData names(7): superkingdom phylum ... genus species
## colnames(24): MV_FEI1_t1Q14 MV_FEI2_t1Q14 ... MV_MIM5_t2M14
##   MV_MIM5_t3F15
## colData names(22): study_name subject_id ... lactating curator
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (296 rows)
## rowTree: 1 phylo tree(s) (10430 leaves)
## colLinks: NULL
## colTree: NULL

When the counts argument is set to TRUE, relative abundance proportions are multiplied by read depth and rounded to the nearest integer prior to being returned. Also, when multiple resources are requested, the list will contain named elements corresponding to each SummarizedExperiment and/or TreeSummarizedExperiment object.

curatedMetagenomicData("AsnicarF_20.+.relative_abundance", dryrun = FALSE, counts = TRUE, rownames = "short")
## $`2021-10-14.AsnicarF_2017.relative_abundance`
## class: TreeSummarizedExperiment 
## dim: 296 24 
## metadata(1): agglomerated_by_rank
## assays(1): relative_abundance
## rownames(296): Escherichia coli Bifidobacterium bifidum ...
##   Streptococcus gordonii Abiotrophia sp. HMSC24B09
## rowData names(7): superkingdom phylum ... genus species
## colnames(24): MV_FEI1_t1Q14 MV_FEI2_t1Q14 ... MV_MIM5_t2M14
##   MV_MIM5_t3F15
## colData names(22): study_name subject_id ... lactating curator
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (296 rows)
## rowTree: 1 phylo tree(s) (10430 leaves)
## colLinks: NULL
## colTree: NULL
## 
## $`2021-03-31.AsnicarF_2021.relative_abundance`
## class: TreeSummarizedExperiment 
## dim: 633 1098 
## metadata(1): agglomerated_by_rank
## assays(1): relative_abundance
## rownames(633): Phocaeicola vulgatus Bacteroides stercoris ...
##   Pyramidobacter sp. C12-8 Brevibacterium aurantiacum
## rowData names(7): superkingdom phylum ... genus species
## colnames(1098): SAMEA7041133 SAMEA7041134 ... SAMEA7045952 SAMEA7045953
## colData names(24): study_name subject_id ... family treatment
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (633 rows)
## rowTree: 1 phylo tree(s) (10430 leaves)
## colLinks: NULL
## colTree: NULL

6.1.1 mergeData()

To merge the list elements returned from the curatedMetagenomicData() function into a single SummarizedExperiment or TreeSummarizedExperiment object, users will use the mergeData() function, provided elements are of the same dataType.

curatedMetagenomicData("AsnicarF_20.+.marker_abundance", dryrun = FALSE) |>
    mergeData()
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(.fns = ~replace_na(.x, 0))`.
## Caused by warning:
## ! Using `across()` without supplying `.cols` was deprecated in dplyr 1.1.0.
## ℹ Please supply `.cols` instead.
## class: SummarizedExperiment 
## dim: 76639 1122 
## metadata(0):
## assays(1): marker_abundance
## rownames(76639): 39491__A0A395UYA6__EUBREC_0408
##   39491__C4Z9E9__T1815_12341 ... 78448__A0A087CC86__BPULL_0419
##   356829__A0A087E8C8__BITS_5024
## rowData names(0):
## colnames(1122): MV_FEI1_t1Q14 MV_FEI2_t1Q14 ... SAMEA7045952
##   SAMEA7045953
## colData names(26): study_name subject_id ... family treatment

The mergeData() function works for every dataType and will always return the appropriate data structure (a single SummarizedExperiment or TreeSummarizedExperiment object).

curatedMetagenomicData("AsnicarF_20.+.pathway_abundance", dryrun = FALSE) |>
    mergeData()
## class: SummarizedExperiment 
## dim: 16913 1122 
## metadata(0):
## assays(1): pathway_abundance
## rownames(16913): UNMAPPED UNINTEGRATED ... PWY-6277: superpathway of
##   5-aminoimidazole ribonucleotide
##   biosynthesis|g__Massiliomicrobiota.s__Massiliomicrobiota_timonensis
##   PWY-6151: S-adenosyl-L-methionine cycle
##   I|g__Massiliomicrobiota.s__Massiliomicrobiota_timonensis
## rowData names(0):
## colnames(1122): MV_FEI1_t1Q14 MV_FEI2_t1Q14 ... SAMEA7045952
##   SAMEA7045953
## colData names(26): study_name subject_id ... family treatment

This is useful for analysis across entire studies (e.g. meta-analysis); however, when doing analysis across individual samples (e.g. mega-analysis) the returnSamples() function is preferable.

curatedMetagenomicData("AsnicarF_20.+.relative_abundance", dryrun = FALSE, rownames = "short") |>
    mergeData()
## class: TreeSummarizedExperiment 
## dim: 673 1122 
## metadata(0):
## assays(1): relative_abundance
## rownames(673): Escherichia coli Bifidobacterium bifidum ...
##   Pyramidobacter sp. C12-8 Brevibacterium aurantiacum
## rowData names(7): superkingdom phylum ... genus species
## colnames(1122): MV_FEI1_t1Q14 MV_FEI2_t1Q14 ... SAMEA7045952
##   SAMEA7045953
## colData names(26): study_name subject_id ... family treatment
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (673 rows)
## rowTree: 1 phylo tree(s) (10430 leaves)
## colLinks: NULL
## colTree: NULL

6.2 returnSamples()

The returnSamples() function takes the sampleMetadata data.frame subset to include only desired samples and metadata as input, and returns a single SummarizedExperiment or TreeSummarizedExperiment object that includes only desired samples and metadata. To use this function, filter rows and/or select columns of interest from the sampleMetadata data.frame, maintaining at least one row, and the sample_id and study_name columns. Then provide the subset data.frame as the first argument to the returnSamples() function.

The returnSamples() function requires a second argument dataType (either "gene_families", "marker_abundance", "marker_presence", "pathway_abundance", "pathway_coverage", or "relative_abundance") to be specified. It is often most convenient to subset the sampleMetadata data.frame using dplyr syntax.

sampleMetadata |>
    filter(age >= 18) |>
    filter(!is.na(alcohol)) |>
    filter(body_site == "stool") |>
    select(where(~ !all(is.na(.x)))) |>
    returnSamples("relative_abundance", rownames = "short")
## class: TreeSummarizedExperiment 
## dim: 824 702 
## metadata(0):
## assays(1): relative_abundance
## rownames(824): Prevotella copri Prevotella sp. CAG:520 ...
##   Corynebacterium aurimucosum Corynebacterium coyleae
## rowData names(7): superkingdom phylum ... genus species
## colnames(702): JAS_1 JAS_10 ... YSZC12003_37879 YSZC12003_37880
## colData names(45): study_name subject_id ... inr zigosity
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: a LinkDataFrame (824 rows)
## rowTree: 1 phylo tree(s) (10430 leaves)
## colLinks: NULL
## colTree: NULL

The counts and rownames arguments apply to returnSamples() as well, and can be passed the function. Finally, users should know that any arbitrary columns added to sampleMetadata will be present in the colData of the SummarizedExperiment or TreeSummarizedExperiment object that is returned.

7 Example Analysis

To demonstrate the utility of curatedMetagenomicData, an example analysis is presented below. However, readers should know analysis is generally beyond the scope of curatedMetagenomicData and the analysis presented here is for demonstration alone. It is best to consider the output of curatedMetagenomicData as the input of analysis more than anything else.

7.1 R Packages

To demonstrate the utility of curatedMetagenomicData, the stringr, mia, scater, and vegan packages are needed.

library(stringr)
library(mia)
library(scater)
library(vegan)

7.2 Prepare Data

In our hypothetical study, let’s examine the association of alcohol consumption and stool microbial composition across all annotated samples in curatedMetagenomicData. We will examine the alpha diversity (within subject diversity), beta diversity (between subject diversity), and conclude with a few notes on differential abundance analysis.

7.2.1 Return Samples

First, as above, we use the returnSamples() function to return the relevant samples across all studies available in curatedMetagenomicData. We want adults over the age of 18, for whom alcohol consumption status is known, and we want only stool samples. The select(where... line below simply removes metadata columns which are all NA values – they exist in another study but are all NA once subsetting has been done. Lastly, the "relative_abundance" dataType is requested because it contains the relevant information about microbial composition.

alcoholStudy <-
    filter(sampleMetadata, age >= 18) |>
    filter(!is.na(alcohol)) |>
    filter(body_site == "stool") |>
    select(where(~ !all(is.na(.x)))) |>
    returnSamples("relative_abundance", rownames = "short")

7.2.2 Mutate colData

Most of the values in the sampleMetadata data.frame (which becomes colData) are in snake case (e.g. snake_case) and don’t look nice in plots. Here, the values of the alcohol variable are made into title case using stringr so they will look nice in plots.

colData(alcoholStudy) <-
    colData(alcoholStudy) |>
    as.data.frame() |>
    mutate(alcohol = str_replace_all(alcohol, "no", "No")) |>
    mutate(alcohol = str_replace_all(alcohol, "yes", "Yes")) |>
    DataFrame()

7.2.3 Agglomerate Ranks

Next, the splitByRanks function from mia is used to create alternative experiments for each level of the taxonomic tree (e.g. Genus). This allows for diversity and differential abundance analysis at specific taxonomic levels; with this step complete, our data is ready to analyze.

altExps(alcoholStudy) <-
    splitByRanks(alcoholStudy)

7.3 Alpha Diversity

Alpha diversity is a measure of the within sample diversity of features (relative abundance proportions here) and seeks to quantify the evenness (i.e. are the amounts of different microbes the same) and richness (i.e. are they are large variety of microbial taxa present). The Shannon index (H’) is a commonly used measure of alpha diversity, it’s estimated here using the estimateDiversity() function from the mia package.

To quickly plot the results of alpha diversity estimation, the plotColData() function from the scater package is used along with ggplot2 syntax.

alcoholStudy |>
    estimateDiversity(assay.type = "relative_abundance", index = "shannon") |>
    plotColData(x = "alcohol", y = "shannon", colour_by = "alcohol", shape_by = "alcohol") +
    labs(x = "Alcohol", y = "Alpha Diversity (H')") +
    guides(colour = guide_legend(title = "Alcohol"), shape = guide_legend(title = "Alcohol")) +
    theme(legend.position = "none")

Figure 1: Alpha Diversity – Shannon Index (H’)

The figure suggest that those who consume alcohol have higher Shannon alpha diversity than those who do not consume alcohol; however, the difference does not appear to be significant, at least qualitatively.

7.4 Beta Diversity

Beta diversity is a measure of the between sample diversity of features (relative abundance proportions here) and seeks to quantify the magnitude of differences (or similarity) between every given pair of samples. Below it is assessed by Bray–Curtis Principal Coordinates Analysis (PCoA) and Uniform Manifold Approximation and Projection (UMAP).

7.4.1 Bray–Curtis PCoA

To calculate pairwise Bray–Curtis distance for every sample in our study we will use the runMDS() function from the scater package along with the vegdist() function from the vegan package.

To quickly plot the results of beta diversity analysis, the plotReducedDim() function from the scater package is used along with ggplot2 syntax.

alcoholStudy |>
    runMDS(FUN = vegdist, method = "bray", exprs_values = "relative_abundance", altexp = "genus", name = "BrayCurtis") |>
    plotReducedDim("BrayCurtis", colour_by = "alcohol", shape_by = "alcohol") +
    labs(x = "PCo 1", y = "PCo 2") +
    guides(colour = guide_legend(title = "Alcohol"), shape = guide_legend(title = "Alcohol")) +
    theme(legend.position = c(0.90, 0.85))
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.