Methods

Reference for publications:

CosmosID Metagenomics Cloud, app.cosmosid.com, CosmosID Inc., www.cosmosid.com

Taxonomic Classification Methods:

The system utilizes a high performance data-mining k-mer algorithm that rapidly disambiguates millions of short sequence reads into the discrete genomes engendering the particular sequences. The pipeline has two separable comparators: the first consists of a pre-computation phase for reference databases and the second is a per-sample computation. The input to the pre-computation phase are databases of reference genomes, virulence markers and antimicrobial resistance markers that are continuously curated by CosmosID scientists. The output of the pre-computational phase is a phylogeny tree of microbes, together with sets of variable length k-mer fingerprints (biomarkers) uniquely associated with distinct branches and leaves of the tree. The second per-sample computational phase searches the hundreds of millions of short sequence reads, or alternatively contigs from draft de novo assemblies, against the fingerprint sets. This query enables the sensitive yet highly precise detection and taxonomic classification of microbial NGS reads. The resulting statistics are analyzed to return the fine-grain taxonomic and relative abundance estimates for the microbial NGS datasets. To exclude false positive identifications the results are filtered using a filtering threshold derived based on internal statistical scores that are determined by analyzing a large number of diverse metagenomes. The same approach is applied to enable the sensitive and accurate detection of genetic markers for virulence and for resistance to antibiotics.

Functional Classification Methods:

Initial QC, adapter trimming and preprocessing of metagenomic sequencing reads are done using BBduk (1). The quality controlled reads are then subjected to a translated search against a comprehensive and non-redundant protein sequence database, UniRef 90. The UniRef90 database, provided by UniProt (2), represents a clustering of all non-redundant protein sequences in UniProt, such that each sequence in a cluster aligns with 90% identity and 80% coverage of the longest sequence in the cluster. The mapping of metagenomic reads to gene sequences are weighted by mapping quality, coverage and gene sequence length to estimate community wide weighted gene family abundances as described by Franzosa et al (3). Gene families are then annotated to MetaCyc (4) reactions (Metabolic Enzymes) to reconstruct and quantify MetaCyc (4) metabolic pathways in the community as described by Franzosa et al (3). Furthermore, the UniRef_90 gene families are also regrouped to GO terms (5) in order to get an overview of GO functions in the community. Lastly, to facilitate comparisons across multiple samples with different sequencing depths, the abundance values are normalized using Total-sum scaling (TSS) normalization to produce "Copies per million" (analogous to TPMs in RNA-Seq) units.

References:
Bushnell, B. (2021). BBDuk Guide - DOE Joint Genome Institute. Retrieved 1 August 2021, from https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbduk-guide/

UniProt: the universal protein knowledgebase. (2016). Nucleic Acids Research, 45(D1), D158-D169. doi: 10.1093/nar/gkw1099

Franzosa, E., McIver, L., Rahnavard, G., Thompson, L., Schirmer, M., & Weingart, G. et al. (2018). Species-level functional profiling of metagenomes and metatranscriptomes. Nature Methods, 15(11), 962-968. doi: 10.1038/s41592-018-0176-y

Caspi, R., Foerster, H., Fulcher, C., Kaipa, P., Krummenacker, M., & Latendresse, M. et al. (2007). The MetaCyc Database of metabolic pathways and enzymes and the BioCyc collection of Pathway/Genome Databases. Nucleic Acids Research, 36(Database), D623-D631. doi: 10.1093/nar/gkm900

Carbon, S., Ireland, A., Mungall, C., Shu, S., Marshall, B., & Lewis, S. (2008). AmiGO: online access to ontology and annotation data. Bioinformatics, 25(2), 288-289. doi: 10.1093/bioinformatics/btn615

16S ASV based Classification Methods

The CosmosID-HUB Microbiome’s 16S workflow implements the DADA2 algorithm(3) as its core engine and utilizes the Nextflow ampliseq pipeline(1) definitions to run it on our cloud infrastructure. Briefly, primer removal is done with Cutadapt (4), and quality trimming parameters are passed to DADA2 to ensure that the median quality score over the length of the read exceeds a certain Phred score threshold. Within DADA2, forward and reverse reads are each trimmed to a uniform length based on the quality of reads in the sample—higher quality data will generally result in longer reads. DADA2 uses machine learning with a parametric error model to learn the error rates for the forward and reverse reads, based on the premise that correct sequences should be more common than any particular error-variant. DADA2 then applies its core sample inference algorithm to the filtered and trimmed data, applying these learned error models. Paired-end reads are then merged if they have at least 12 bases of overlap and are identical across the entire overlap.

The resulting table of sequences and observed frequencies is filtered to remove chimeric sequences (those that exactly match a combination of more-prevalent “parent” sequences). Taxonomy and species-level identification (where possible) are conducted with DADA2’s naive Bayesian classifier, using the Silva version 138 database.

Lastly, the predicted functional potential of the community was profiled using PICRUST2 (5)(6)(7)(8)(9). Briefly, PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) is a tool that predicts functional capabilities and abundances of a microbial community based on the observed amplicon (marker gene) content. Functional capabilities are given by EC classifiers, or MetaCyc ontologies, and these can be aggregated to predict pathways that are likely present in a given sample.

References:

1. Straub, D. et al. Interpretations of Environmental Microbial Community Studies Are Biased by the Selected 16S rRNA (Gene) Amplicon Sequencing Pipeline. Front. Microbiol. 11, 1–18 (2020).
2. Callahan, B. J., McMurdie, P. J. & Holmes, S. P. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 11, 2639–2643 (2017).
3. Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).
4. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17, 10 (2011).
5. Douglas, G. M. et al. PICRUSt2 for prediction of metagenome functions. Nat. Biotechnol. 38, 685–688 (2020).
6. Barbera, P. et al. EPA-ng: Massively Parallel Evolutionary Placement of Genetic Sequences. Syst. Biol. 68, 365–369 (2019).
7. Czech, L., Barbera, P. & Stamatakis, A. Genesis and Gappa: processing, analyzing and visualizing phylogenetic (placement) data. Bioinformatics 36, 3263–3265 (2020).
8. MIRARAB, S., NGUYEN, N. & WARNOW, T. SEPP: SATé-Enabled Phylogenetic Placement. in Biocomputing 2012 247–258 (WORLD SCIENTIFIC, 2011). doi:10.1142/9789814366496_0024.
9. Louca, S. & Doebeli, M. Efficient comparative phylogenetics on large trees. Bioinformatics 34, 1053–1055 (2018).
10. Ye, Y. & Doak, T. G. A Parsimony Approach to Biological Pathway Reconstruction/Inference for Genomes and Metagenomes. PLoS Comput. Biol. 5, e1000465 (2009).
11. Chiarello, M., McCauley, M., Villéger, S. & Jackson, C. R. Ranking the biases: The choice of OTUs vs. ASVs in 16S rRNA amplicon data analysis has stronger effects on diversity measures than rarefaction and OTU identity threshold. PLoS One 17, 1–19 (2022).

I6S OTU based Taxonomic Classification Methods:

For taxonomic profiling based of amplicon data, the CosmosID 16S data analysis pipeline starts with preprocessing of the raw reads from either paired-end or single-end Fastq files through read-trimming to remove adapters as well as reads and bases of low quality. If the reads are in a paired-end format, the forward and reverse overlapping pairs are joined together; the unjoined R1 and R2 reads are then added to the end of the file. The file is then converted to Fasta format and used as input for OTU picking. OTUs are identified against the CosmosID curated 16S database using a closed-reference OTU picker and 97% sequence similarity through the QIIME framework. The final results are then presented in tabular format with the taxonomic names, OTU IDs, frequency, and relative abundance. Results can be downloaded or compared to other 16S samples for visualizations through the CosmosID Comparative Analysis tool.