Do you cov me? Effect of coverage reduction on species identification and genome reconstruction in complex biological matrices by metagenome shotgun high-throughput sequencing

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.

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.

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.

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 by 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”.

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 very good and improved at increasing sequencing depths. Pearson’s correlation coefficient between the two datasets was >0.95 in all subsampling replicates at all sequencing depths, with the exception of samples F1 and F2, for which the correlation reached 0.95 only for subsamples of 500,000 or more reads.

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 observed species, due to the increased 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 moderate.

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 F). 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 6). 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 6. 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 7 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 7. 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.