Metagenomic Information Recovery from Human Stool Samples Is Influenced by Sequencing Depth and Profiling Method

Abstract

Sequencing of the 16S rRNA gene (16S) has long been a go-to method for microbiome characterization due to its accessibility and lower cost compared to shotgun metagenomic sequencing (SMS). However, 16S sequencing rarely provides species-level resolution and cannot provide direct assessment of other taxa (e.g., viruses and fungi) or functional gene content. Shallow shotgun metagenomic sequencing (SSMS) has emerged as an approach to bridge the gap between 16S sequencing and deep metagenomic sequencing. SSMS is cost-competitive with 16S sequencing, while also providing species-level resolution and functional gene content insights. In the present study, we evaluated the effects of sequencing depth on marker gene-mapping- and alignment-based annotation of bacteria in healthy human stool samples. The number of identified taxa decreased with lower sequencing depths, particularly with the marker gene-mapping-based approach. Other annotations, including viruses and pathways, also showed a depth-dependent effect on feature recovery. These results refine the understanding of the suitability and shortcomings of SSMS, as well as annotation tools for metagenomic analyses in human stool samples. Results may also translate to other sample types and may open the opportunity to explore the effect of sequencing depth and annotation method.

Keywords: alignment; marker gene; microbiome; shallow sequencing; shotgun metagenomic sequencing; virome.

Figure 1

Boxplots of the α-diversity of ten stool samples from healthy subjects at various sequencing depths, including 5 Gb (16.67 M reads), 3 Gb (10.00 M reads), 1 Gb (3.00 M reads), 0.75 Gb 2.50 M reads), 0.5 Gb (1.65 M reads), 0.25 Gb (0.85 M reads), and 0.1 Gb (0.34 M reads). The 16S sequencing analysis at the genus level was performed for comparison. (A) Number of observed genera and Shannon diversity at the genus level for the marker gene-mapping and alignment methods. (B) Number of observed species and Shannon diversity at the species level for the marker gene-mapping and alignment methods. Outliers are shown.

Figure 2

PCoA plots of the β-diversity (weighted Bray–Curtis distances) of ten stool samples from healthy subjects at various sequencing depths including 5 Gb (16.67 M reads), 3 Gb (10.00 M reads), 1 Gb (3.00 M reads), 0.75 Gb (2.50 M reads), 0.5 Gb (1.65 M reads), 0.25 Gb (0.85 M reads), and 0.1 Gb (0.34 M reads). The 16S rRNA (16S) V4 hypervariable region sequencing analysis at the genus level was performed for comparison. (A) SMS data annotated using the marker gene mapping and alignment methods, in comparison with the 16S sequencing information described above at the genus level. (B) SMS data annotated using the marker gene-mapping method at the genus level. (C) SMS data annotated using the alignment method at the genus level. (D) SMS data annotated using the marker gene mapping and alignment methods at the species level. (E) SMS data annotated using the marker gene-mapping method at the species level. (F) SMS data annotated using the alignment method at the species level. Clustering of data by subjects is highlighted using circles or ellipses.

Figure 3

Line charts of the average relative abundances (%) (genus level) of selected stool taxa based upon SMS data annotated using marker gene-mapping- and alignment-based methods. Data were annotated at various sequencing depths including 5 Gb (16.67 M reads), 3 Gb (10.00 M reads), 1 Gb (3.00 M reads), 0.75 Gb (2.50 M reads), 0.5 Gb (1.65 M reads), 0.25 Gb (0.85 M reads), and 0.1 Gb (0.35 M reads). The 16S sequencing data were included for comparison. (A) Bifidobacterium SMS data annotated using marker gene-mapping- and alignment-based methods. (B) Alistipes SMS data annotated using marker gene-mapping- and alignment-based methods. (C) Roseburia SMS data annotated using marker gene-mapping- and alignment-based methods. (D) Lactobacillus SMS data annotated using marker gene-mapping- and alignment-based methods. Standard error is shown by error bars.

Figure 4

Correlation plots of the relative abundances of bacterial genera annotated using marker gene-mapping and alignment methods. (A) Correlations performed with data obtained at 5 Gb (16.67 M reads) vs. 3 Gb (10.00 M reads). (B) Correlations performed with data obtained at 5 Gb vs. 1 Gb (3.00 M reads). (C) Correlations performed with data obtained at 5 Gb vs. 0.75 Gb (2.50 M reads). (D) Correlations performed with data obtained at 5 Gb vs. 0.5 Gb (1.65 M reads). (E) Correlations performed with data obtained at 5 Gb vs. 0.25 (0.85 M reads). (F) Correlations performed with data obtained at 5 Gb vs. 0.1 Gb (0.35 M reads). (G) Correlations performed with data obtained at 5 Gb vs. 16S sequencing data.

Figure 5

Viruses identified in healthy individuals at various sequencing depths including 5 Gb (16.67 M reads), 3 Gb (10.00 M reads), 1 Gb (3.00 M reads), 0.75 Gb (2.50 M reads), 0.5 Gb (1.65 M reads), 0.25 Gb (0.85 M reads), and 0.1 Gb (0.35 M reads). (A) Boxplots of the total number of viruses identified in ten healthy individuals at the various sequencing depths described. Figure also shows barplots of normalized counts (%), consisting of total read counts per virus at the various sequencing depths, divided by the genome size (bp). (B) Arthrobacter phage Mendel (taxid 2484218) total read counts/genome size (bp) (%). (C) CrAssphage (taxid 2212563) total read counts/genome size (bp) (%). (D) Faecalibacterium phage FP oengus (taxid 2070188) total read counts/genome size (bp) (%). (E) Lactococcus phage 16802 (taxid 2029659) total read counts/genome size (bp) (%). (F) Poophage MBI 2015a (taxid 1926504) total read counts/genome size (bp) (%).

Figure 6

Boxplots of the total number of functional pathways identified in ten stool samples across various sequencing depths including 5 Gb (16.67 M reads), 3 Gb (10.00 M reads), 1 Gb (3.00 M reads), 0.75 Gb (2.50 M reads), 0.5 Gb (1.65 M reads), 0.25 Gb (0.85 M reads), and 0.1 Gb (0.34 M reads). (A) Number of UniProt pathways across the various sequencing depths. (B) Number of MetaCyc pathways across the various sequencing depths.