Do you cov me? Effect of coverage reduction on metagenome shotgun sequencing studies

Samples description and DNA extraction

The following samples were used in the present work: the mock community DNA sample “20 Strain Staggered Mix Genomic Material” ATCC^® MSA-1003^TM (short name: A1), two biological medicines (B1 and B2), two horse fecal samples (F1 and F2), three food samples (M1, M2, and M3), and two human fecal samples (V1 and V2).

Biological medicines were two different lots of live attenuated MPRV vaccine, widely used for immunization against measles, mumps, rubella and chickenpox in infants. Lyophilised vaccines were resuspended in 500 μl sterile water for injection and DNA extracted from 250 μl using Maxwell^® 16 Instrument and the Maxwell^® 16 Tissue DNA Purification Kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. The vaccine composition declared by the producer is the following: live attenuated viruses: 1) Measles (ssRNA) Swartz strain, cultured in embryo chicken cell cultures; Mumps (ssRNA) strain RIT 4385, derived from the Jeryl Linn strain, cultured in embryo chicken cell cultures; Rubella (ssRNA) Wistar RA 27/3 strain, grown in human diploid cells (MRC-5); Varicella (dsDNA) OKA strain grown in human diploid cells (MRC-5).

Horse feces from two individuals were processed as follows: 100 mg of starting material stored in 70% ethanol were used for DNA extraction using the QIAamp PowerFecal DNA Kit (QIAGEN GmbH, Hilden, Germany), according to the manufacturer's instructions.

Food samples were raw materials of animal and plant origin, used to industrially prepare bouillon cubes. DNA extractions from those three samples were performed starting from 2 grams of material each, using the DNeasy mericon Food Kit (QIAGEN GmbH, Hilden, Germany), according to the manufacturer's instructions. The declared sample composition was Agaricus bisporus for M1, spice (Piper nigrum) for M2 and mix of animal extracts for M3.

The mock community declared components are: 0.18% Acinetobacter baumannii (ATCC 17978), 0.02% Actinomyces odontolyticus (ATCC 17982), 1.80% Bacillus cereus (ATCC 10987), 0.02% Bacteroides vulgatus (ATCC 8482), 0.02% Bifidobacterium adolescentis (ATCC 15703), 1.80% Clostridium beijerinckii (ATCC 35702), 0.18% Cutibacterium acnes (ATCC 11828), 0.02% Deinococcus radiodurans (ATCC BAA-816), 0.02% Enterococcus faecalis (ATCC 47077), 18.0% Escherichia coli (ATCC 700926), 0.18% Helicobacter pylori (ATCC 700392), 0.18% Lactobacillus gasseri (ATCC 33323), 0.18% Neisseria meningitidis (ATCC BAA-335), 18.0% Porphyromonas gingivalis (ATCC 33277), 1.80% Pseudomonas aeruginosa (ATCC 9027), 18.0% Rhodobacter sphaeroides (ATCC 17029), 1.80% Staphylococcus aureus (ATCC BAA-1556), 18.0% Staphylococcus epidermidis (ATCC 12228), 1.80% Streptococcus agalactiae (ATCC BAA-611), 18.0% Streptococcus mutans (ATCC 700610).

DNA purity and concentration were estimated using a NanoDrop Spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA) and Qubit 2.0 fluorimeter (Invitrogen, Carlsbad, CA, USA).

Human fecal samples V1 and V2 derive from a study investigating the virome composition of feces of uncontacted Amerindians²⁶. Data are publicly available on Sequence Read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra/). The two samples with the highest sequencing depth were chosen; accession numbers are SRR6287060 and SRR6287079, respectively.

Whole metagenome DNA library construction and sequencing

DNA library preparations were performed according to manufacturer’s protocol, using the kit Ovation^® Ultralow System V4 1–96 (Nugen, San Carlos, CA). Library prep monitoring and validation were performed both by Qubit 2.0 fluorimeter (Invitrogen, Carlsbad, CA, USA) and Agilent 2100 Bioanalyzer DNA High Sensitivity Analysis kit (Agilent Technologies, Santa Clara, CA). Obtained DNA concentrations were as follows: A1 8 ng/µl (total amount = 640 ng), B1 10.7 ng/µl (total amount = 535 ng), B2 9.41 ng/µl (total amount = 470.5 ng), F1 42.3 ng/µl (total amount = 4,230 ng), F2 22.6 ng/µl (total amount = 2,260 ng), M1 16.6 ng/µl (total amount = 1,494 ng), M2 1.87 ng/µl (total amount = 168.3 ng), M3 16 ng/µl (total amount = 640 ng).

Cluster generation was then performed on Illumina cBot and flowcell HiSeq SBS V4 (250 cycle), and sequenced on HiSeq2500 Illumina sequencer producing 125bp paired-end reads.

Samples F1 and F2 were loaded on flowcell HiSeq Rapid SBS Kit v2 (500 cycles) producing 250bp paired-end reads. The estimated library insert sizes were: 539 bp (A1), 531 bp (B1), 536 bp (B2), 620 bp (F1), 620 bp (F2), 342 bp (M1), 178 bp (M2), 496 bp (M3). Samples were sequenced in different runs and pooled with other libraries of similar insert sizes.

The CASAVA Illumina Pipeline version 1.8.2 was used for base-calling and de-multiplexing. Adapters were masked using cutadapt²⁷. Masked and low quality bases were filtered using erne-filter version 1.4.6.²⁸.

Bioinformatics analysis

The bioinformatics analysis performed in the present work are summarized in Figure 1; a standard pipeline for reproducing the main steps of analysis is available on GitHub (http://www.doi.org/10.5281/zenodo.2593798).

Figure 1. Workflow of the main bioinformatics analysis performed in the present work.

Since different read lengths among samples may constitute an additional confounder in analysis, 250 bp long reads belonging to F1, F2, V1 and V2 were trimmed to a length of 125bp using fastx-toolkit version 0.0.13 (http://hannonlab.cshl.edu/fastx_toolkit/) before analysis.

Reduction in coverage was simulated by randomly sampling a fixed number of reads from the full set of reads. Subsamples of 10,000, 25,000, 50,000, 100,000, 250,000, 500,000 and 1,000,000 reads were extracted from the raw reads using seqtk version 1.3. To estimate the variability due to random effects, subsampling was replicated five times for each simulated depth and 99% confidence limits were estimated and plotted.

To classify the largest possible number of prokaryotes, eukaryotes and viruses, reads were classified against the complete NCBI nt database using kraken2, version 2.0.6²⁹. The nt database was converted to kraken2 format using the built-in kraken2-build script with default parameters. Among the most significant parameters, kmer size for the database is by default set to 35 and the minimizer length to 31. A simplified representation of species composition was obtained using Krona³⁰. To obtain accurate species abundances Bracken, version 2.2³¹ was used on species supported by at least 10 reads; since the reads used in all the experiments were 125bp, the bracken database was built using 125bp kmers.

The threshold for declaring a species as present was set according to results of a performance analysis on the mock community (A1) for which species presence and abundance was known. Performance was assessed using F1-score, calculated as 2*TP/(2*TP+FP+FN), as previously reported³². F1-score is a measure used in performance analysis when the number of true negatives is extremely high or unknown.

The effect of the selected database on reads classification was assessed on the mock community full sequencing experiment, by observing the variation in present and absent species when using different databases. In addition to the nt database built explicitly for this study, the standard kraken2 database, the minikraken2 v1 database and the minikraken2_v2 were used. The standard kraken2 database contains complete genomes in RefSeq for the bacterial, archaeal, and viral domains, along with the human genome and a collection of known vectors (UniVec_Core), and the minikraken2 v1 database contains RefSeq bacteria, archaea, and viral libraries, and the minikraken2_v2 database contains RefSeq bacteria, archaea, and viral libraries and the GRCh38 human genome.

Bracken database was built for 125bp kmers for the standard database. Minikraken2 instead is distributed as a prebuilt database, from which it is not possible to build the bracken database, but for which bracken databases with kmers distribution of 100bp, 150bp and 200 kmers are available. kmers 100bp and 150bp were tested, since they are the closest to the read length used in this study.

Observed number of taxa, Shannon’s diversity index³³ and Pielou’s index³⁴ were estimated using the R package vegan version 2.4.2³⁵ or base R, version 3.3.3³⁶ functions. The number of observed taxa was computed as the number of species passing the detection threshold.

Shannon diversity index is estimated as

Where N is the total number of species and p_i is the frequency of the species i.

Pielou’s evenness index is estimated as

Where H is Shannon’s diversity index and S is the total number of observed species. The value ln S corresponds to the maximum possible value of H, observed when all species have the same frequency.

The effect of sequencing depth on beta-diversity was assessed using the procrustes and protest functions, implemented in vegan.

Assembly of the metagenome was performed using megahit version 1.1.2³⁷ with default parameters, with kmer sizes varying as follows: 21, 29, 39, 59, 79, 99, 119, 141. Reconstructed contigs were binned at the species level using kraken2, and only contigs assigned to species above the detection threshold were used for further analysis. Completeness of the assemblies of each species was assessed using BUSCO³⁸. For each species, the proportion of the reconstructed genes was measured as the proportion of genes that were fully reconstructed, plus the proportion of genes that were partially reconstructed. For each sample, results were then averaged over detected species to provide the average proportion of reconstructed genes. BUSCO analysis was performed on prokaryotic database for all the samples with the exception of M1 (predominantly composed by fungi) for which the fungal database was used.

Unless otherwise specified, all the analysis were performed using R 3.3.3³⁶.

Determination of detection threshold

The mock community sample “20 Strain Staggered Mix Genomic Material” (ATCC^® MSA-1003^TM) was used as a reference to control performance of sequencing and classification procedures at various depth. The community includes a total of 20 bacterial species, of which 5 have a frequency of 0.02%, 5 a frequency of 0.18%, 5 a frequency of 1.8% and 5 a frequency of 18%.

Results of the performance analysis on the mock dataset are shown in Table 1. The highest F1 score (0.8) was obtained when applying a 0.1% threshold. Using this threshold, 14 species were correctly identified while 6 of them were not detected. Five out of the 6 undetected species had a nominal frequency of 0.02%; the sixth undetected species was Helicobacter pylori, with a nominal frequency of 0.18%, for which we recorded a frequency of 0.096%, below the 0.1% threshold. The only false positive was Shigella flexneri, a species highly related to Escherichia coli³⁹, that was observed at a frequency of 0.128%. Based on these results we used a threshold of 0.1% for declaring a species as present in a sample.

Sample composition

Summary statistics for the samples included in the study are shown in Table 2.

The number of reads obtained in the samples selected for the present study ranged from slightly more than 1 million (sample V1) to more than 12 million (sample F1). The number of species identified in each sample ranged from 4 in sample B1 to 138 in sample M2. Figure 2 summarizes the composition of each sample at the Phylum level. Viruses are aggregated at the division level. Only phyla more abundant than 1% were plotted. Reads that were either unclassified or assigned to rare phyla were aggregated under the name “Unknown/Other”. Samples B1, B2 and M3 where mainly composed of Chordata, sample M1 was mostly composed of Basidiomycota, and sample V2 was mainly composed of Viruses. Samples F1, F2 and, to a lesser extent, M2 were characterized by a large proportion of reads unclassified or assigned to rare phyla. For a more detailed view of raw taxonomy composition, interactive html Chrona are available for download on Open Science Framework (https://osf.io/y7c39/), under the project “Do you cov me”, DOI: 10.17605/OSF.IO/Y7C39.

Figure 2. Phylum composition of the samples.

Only phyla represented by at least 1% of the reads are shown. Viruses are presented at division level. Unclassified reads and reads assigned to rare phyla are aggregated under the name “Unknown/Other”.

Effect of the choice of database

The effect of the selected database on reads classification was assessed only on the mock community. Results are shown in Table 3. All the Spearman correlation coefficients were >0.9 (not shown). Estimated abundances according to the minikraken2 v1 database were very similar to those obtained according to the minikraken2 v2 database (irrespective of the kmer length which slightly altered the results). The main differences between the two were observed for species not declared in the mock community (i.e. false positives, indicated by ND in Abundance column), such as Homo sapiens, Bacillus anthracis and Salmonella enterica.

Homo sapiens contamination was detected with the Minikraken2 v2 (the only Minikraken2 containing human sequences) and using the standard database. No contamination was detected using the nt database, for which Homo sapiens was recorded with a frequency of 0.086% and was therefore below the detection threshold of 0.1%. Bacillus anthracis was only detected when using the Minikraken2 v2 database. Salmonella enterica was detected in all the experiments involving the Minikraken2 databases, but it was the only species for which substantial variation in the estimated frequencies was observed, ranging from 0.697% to 1.573% according to Minikraken2 v1 with 150 kmers and Minikraken2 V2 with 100 kmers, respectively.

False positives were generally low: the total contribution of false positives ranged from 0.168% for nt to 2.628% for Minikraken2 v2 100 kmers. However, none of the databases was immune to false positives. The nt database showed only one false positive, Shigella flexneri, while minikraken and standard databases showed more than one false positive each.

The estimated abundances of the declared species were in excellent agreement across databases, with some exceptions. For Escherichia coli, with declared abundance of 18%, estimated abundances ranged from 20.043% to 25.954% when using nt and the standard database, respectively. Staphylococcus epidermidis, with declared abundance of 18%, was estimated at 13.112% by the standard database and 16.701% by nt. Finally, Clostridium beijerinckii, declared at 1.8% was estimated at 2.524% by the standard Kraken databse and at 4.69% by nt database.

Species abundance

The effect of reducing sequencing depth on the accuracy of taxonomical classification was assessed by using the mock community, given its known composition. Expected and observed abundancies of the 20 mock species maintained a high correlation (r=0.94) even when decreasing sequencing to 10,000 reads (Figure 3). However, decreasing sequencing depth caused an increase in uncertainty, as shown by the broader confidence intervals for lower depths, particularly for rare species.

Figure 3. Observed and expected abundance of bacterial species present in the mock community “20 Strain Staggered Mix Genomic Material” (ATCC^® MSA-1003^TM) at varying sequencing depths.

In red, species identified at frequency lower than the selected threshold of 0.1% and arbitrarily plotted at 0.002%. Error bars represent 95% confidence intervals obtained from five resampling experiments. Both axes are plotted in log scale to facilitate visualization of rare species.

We also measured the correlation in species abundance between the full and reduced datasets in all the samples (Figure 4). The correlation between the two quantities was in general high and improved at increasing sequencing depths. The average Spearman’s correlation coefficient between full and reduced samples ranged from 0.71 in the 25,000 reads subsample to 0.94 in the 1,000,000 reads subsample.

Figure 4. Correlation of species abundance estimated using full and reduced datasets.

Data for all the five subsampled replicates are plotted. Each point (colored by sample of origin) represents a given species. Both axes are plotted in log scale to facilitate visualization of rare species.

Diversity analysis

We further evaluated the impact of reducing sequencing depth on several diversity measures, such as the observed number of taxa, Shannon’s diversity index and Pielou’s diversity index (Figure 5).

Figure 5.

Effect of reduction of sequencing depth on: A) Observed number of taxa, B) Shannon’s diversity index, C) Pielou’s diversity index, and D) Total length of de novo assembly. In all panels X axis is in log scale and Y axis is in linear scale with the exception of panel F, in which both axes are in log scale. Shaded areas represent the confidence limits of resampling experiments. The number of reads used for the analysis is shown from the smallest (10,000) on the left, to the full dataset on the right.

Samples F1, F2, M2 and V1 had more than 50 taxa, and all the remaining samples had less than 25 (Table 2 and Figure 5A). Downsampling only produced significant differences in the four samples with high number of observed taxa (panel A). In samples F1 and F2, intermediate levels of downsampling (e.g. 100,000 reads) caused an increase in the number of species exceeding the 0.1% abundance threshold.

Shannon’s diversity index (panel B) is a widely used method to assess biological diversity of ecological and microbiological communities. The effect of sequencing depth on Shannon’s diversity index is negligible for most samples, with the exception of the samples with the richest species composition (F1, F2, and M2) for which downsampling led to a significant variation in the estimate.

Pielou’s index (panel C) is a measure of the species’ distribution evenness. Values close to 1 denote equifrequent species, and lower values denote uneven distribution of species relative abundance. The effect of the number of reads on Pielou’s index is negligible for all samples.

Effect of sequencing depth on beta-diversity using procrustes analysis is shown in Figure 6. Procrustes analysis is shown between the full dataset (starting point of the arrows) and each replicate of each of the reduced dataset (arrival points of the arrowheads). The smallest dataset was excluded for the analysis because it did not have sufficient information for procrustes analysis. The relative positions between the full dataset and the reduced datasets tended to remain similar across different subsampling, with some notable exception, such as V2 in the 250,000 reads sample and A1 in the 500,000 sample. F1 and F2, always clustered together. Correlation was good between all the matrices, but it was excellent (r>0.9) only for the comparison of the full dataset with the 1,000,000 reads subsample.

Figure 6. Procrustes analysis.

Arrows start from the NMDS matrix of the full sample and the arrowhead ends at the NMDS matrices of the five subsampled replicates for each sample. The correlation between the matrices are shown at the bottom right of each plot. All correlations were significant (p<=0.001).

Completeness of de novo assembly

We investigated the effect of coverage reduction on the completeness of de novo assembly. We reconstructed the metagenome of the full and reduced datasets and compared the completeness of the reconstructed genomes. Results are summarized in Figure 5 (panel D). As expected, the size of the assembly was strongly influenced by the sequencing depth. Assembly size for the full dataset ranged from less than 1 Mb (V2) to nearly 100 Mb (F1 and F2). A decrease in the sequencing depth led to a steady decrease in assembly size in all samples. At 1,000,000 reads the size ranged from slightly more than 100 kb (V2) to slightly more than 10Mb (A1 and M1).

BUSCO analysis³⁸ was used as an additional measure to assess the completeness of the reconstructed metagenome.

First, we assessed the performance on the A1 mock community for the full set of reads (4,969,245 reads) and for the largest subset, (1,000,000 reads). Genomes of species present at 0.02% and 0.18% were not reconstructed, while genomes of species present at 1.8% or 18% were reconstructed. The proportion of BUSCO genes reconstructed using the whole set of reads ranged 59%–99% for species present at 1.8% and 93%–99% in the most abundant species (Figure 7). At 1,000,000 reads the proportion of reconstructed BUSCO genes dropped to 0.7%–3% for species present at 1.8% while it ranged 93%–99% for the species present at 18%.

Figure 7. Proportion of BUSCO genes reconstructed in the full set of reads, and in a subset of 1,000,000 reads.

Error bars represent 95% confidence intervals based on five subsampling experiments. 1.8%: Species with nominal abundance of 1.8%; 18%: Species with nominal abundance of 18%; Full: results using full set of reads; Reduced; results using a subsample of 1,000,000 reads.

In addition, we plot in Figure 8 the proportion of reconstructed genes in full (X axis) and reduced (Y axis) datasets obtained by randomly sampling 1,000,000 reads. The proportion of reconstructed BUSCO genes is very low even in the full samples, indicating that in general the sequencing depth is still too low to obtain an accurate reconstruction of the metagenome. Only in samples A1 and M1, the average proportion of BUSCO genes reconstructed in the full sample was greater than 10%. Reducing sequencing depth to 1,000,000 reads significantly lowered the proportion of reconstructed genes in all the samples, as testified by the fact that all the points and their confidence limits lie below the diagonal.

Figure 8. Completeness of the BUSCO genes in the full dataset (X axis) and in the largest of the reduced datasets (consisting of 1,000,000 reads, Y axis); error bars are based on the five replicate experiments performed for each sample.

The plot is in log-log scale.

Underlying data

Raw reads generated in the present study are available at NCBI Sequence Read Archive.

Sample A1 is available under accession number SRP174028: https://identifiers.org/insdc.sra/SRP174028.

Samples F1 and F2 are available under accession number SRP163102: https://identifiers.org/insdc.sra/SRP163102.

Samples B1 and B2 are available under accession number SRP163096: https://identifiers.org/insdc.sra/SRP163096;

and samples M1, M2 and M3 are available under accession number SRP163007: https://identifiers.org/insdc.sra/SRP163007.

Extended data

Open Science Framework: Do you cov me. https://doi.org/10.17605/OSF.IO/Y7C39⁴².

This project contains the raw html graphs, produced using Krona.