ConFindr: rapid detection of intraspecies and cross-species contamination in bacterial whole-genome sequence data

ConFindr workflow and implementation

ConFindr determines the presence of contaminating sequencing reads based on the analysis of the set of 53 genes encoding the bacterial ribosomal-protein subunits that are used in the ribosomal multilocus sequence typing scheme (rMLST) (Jolley et al., 2012). The rMLST genes are typically present as single copies and are conserved across the entire bacterial domain, with some exceptions where multiple alleles for a gene exist or no gene exists. ConFindr works on the principle that a genome containing more than one allele for any rMLST gene is contaminated, taking into consideration the known exceptions.

In its first step, ConFindr uses a screening functionality provided in Mash (Ondov et al., 2016) to determine which genera are present in a sample. This screen is done against a custom database derived from the NCBI RefSeq genomes (https://www.ncbi.nlm.nih.gov/refseq/) with one genome representing each species (O’Leary et al., 2015). If more than one genus is detected, ConFindr reports cross-species contamination for the sample and does not proceed further. If only one genus is present, ConFindr creates a genus-specific rMLST database by extracting all rMLST sequences associated with the target genus, excluding genes known to have multiple alleles, and proceeds to attempt to find contamination by searching for multiple alleles of one or more of the benchmark rMLST genes.

To search for multiple alleles, ConFindr begins by using BBDuk (Bushnell, 2014) to extract reads that are likely part of the rMLST gene set. These baited reads are stringently trimmed, again using BBDuk, and then aligned to the rMLST genes using BBMap (Bushnell, 2014). The resulting BAM file is then parsed in order to find ‘Contaminating Single Nucleotide Variants’ (cSNVs)—that is, sites in the pileup where more than one base is present. Since all rMLST genes are known to be present as single copies, the occurrence of multi-base sites in the pileup indicates multiple alleles, and therefore contamination. To be called as a cSNV, at least 2 bases of the minor variant with a Phred score of 20 or greater must be present at that site, and at least 5 percent of bases must support the minor variant (though these parameters can be changed by the user). ConFindr determines that a sample is contaminated if multiple genera are found in the Mash screen step or if three or more cSNVs are found.

In silico dataset creation

To create in silico datasets, we selected complete assemblies from RefSeq for E. coli, S. enterica, and L. monocytogenes (accessions NC_002695.1, NC_003198.1, and NC_003210.1, respectively) and generated variants of these genomes with 100, 500, 1,000, and 2,000 SNVs using a custom script available at https://github.com/lowandrew/MutantCreator. Simulated reads were then created from both variant and base genomes using ART v2.5.8 (commands used to carry this out can be found in Table S2) (Huang et al., 2012) and mixed together in proportions of 0, 1, 5, 10, and 20 percent contamination using scripts found at https://github.com/lowandrew/FastQMixer to a total coverage depth of approximately 60X. Five replicates were created for each mixed read set.

Test datasets

We generated a test dataset of 48 samples comprised of intra- and cross-species mixes of E. coli, S. enterica, and L. monocytogenes isolates with varying levels of relatedness (Table 1). Average nucleotide identity (ANI) was calculated using OrthoANI version 1.4.0 (Lee et al., 2016). This dataset was made both in silico by mixing together reads from previous runs of these isolates using the reformat.sh program of the BBMap package (Bushnell, 2014) to a coverage depth of 80X, as well as by sequencing lab-generated mixes of genomic DNA (gDNA).

To generate WGS data, bacterial isolates were cultured in Brain Heart Infusion (BHI) broth (Oxoid Ltd., Basingstoke, Hampshire, England) for 4 to 6 h at 36 °C, and gDNA was extracted using the Maxwell 16 Cell SEV DNA Purification Kit (Promega, Madison, WI). DNA was quantified using the Quant-it High-Sensitivity DNA Assay Kit (Life Technologies Inc., Burlington, ON). Sequencing libraries were constructed from 1 ng of gDNA using the Nextera XT DNA Sample Preparation Kit (Illumina, Inc., San Diego, CA) and the Nextera XT Index Kit (Illumina, Inc.) according to manufacturers’ instructions. Genomic sequencing was performed on the Illumina MiSeq Platform (Illumina, Inc.) using a 600-cycle MiSeq Reagent kit v3 (Illumina, Inc.).

Nucleotide sequence accession numbers

Raw data have been deposited at DDBJ/EMBL/GenBank under BioProject PRJNA507762. The accession numbers and strain descriptions are listed in Table S1.

Coverage dataset

To evaluate the effect of extremely high coverage depth as well as different sequencing platforms on the performance of ConFindr, we created simulated datasets from RefSeq assemblies (see Table S5) for both the HiSeq and MiSeq platforms using ART v2.5.8 (commands in Table S2) at coverage levels of 50, 100, 200, 300, 400 and 500×, and at contamination levels of 0, 5, 10, and 20 percent.

Genome assemblies and quality metrics

All of the read sets created were also put through the process of de novo assembly. Briefly, reads were quality trimmed using bbduk.sh and error corrected with tadpole.sh (both of the BBMap package (Bushnell, 2014)) and then assembled using SKESA v2.3.0 (Souvorov, Agarwala & Lipman, 2018). Exact commands used to carry this out can be found in Table S2. Genome quality statistics were assessed using QUAST v4.6.3 (Gurevich et al., 2013).

Calculation of number of SNVs between rMLST types

To calculate the number of SNVs between rMLST types within E. coli, S. enterica, and L. monocytogenes, we retrieved all rMLST allele sequences and the list of profiles (accessed at https://pubmlst.org/rmlst/, November 1, 2018). We then extracted sequences for each allele within each rMLST type (1,641 types for L. monocytogenes, 3,062 types for E. coli, and 7,255 for S. enterica). The number of SNVs between every sequence type pair within each species was calculated by aligning each gene in the first type against each gene in the second type using the pairwise2 module in biopython (Cock et al., 2009).

Dataset testing

To detect contamination in the datasets generated, ConFindr v0.4.4 was run on default settings on all samples generated. Kraken v1.0 (Wood & Salzberg, 2014) was run on fastq files that had been trimmed to a quality of 15 with bbduk.sh against the standard Kraken database. Exact commands used to carry this out can be found in Table S2. Strains/species were determined to be present if at least 0.5 percent of classified reads could be assigned to them. CheckM v1.0.11 (Parks et al., 2015) was run with the lineage_wf workflow for each assembly, and samples were called as contaminated if their contamination level was 0.6 percent or greater, as across our test dataset our uncontaminated samples had up to 0.6 percent contamination called by CheckM.

Identification of contaminating SNVs within rMLST genes using ConFindr

As ConFindr is based on finding contaminating SNVs (cSNVs) within the rMLST genes in raw reads, we first looked at the reliability of detection of cSNVs in simulated data with different levels of contamination. Synthetic mutants with 100 to 2000 SNVs relative to reference genomes were generated in silico, and the number of SNVs occurring within rMLST genes in these mutants was calculated. The number of cSNVs found by ConFindr in the in silico datasets at 60 times coverage was compared with the predicted number of SNVs within rMLST genes in the two isolates making up the contaminated sample (Fig. 1). As the relative contamination increased, ConFindr’s estimate of the number of cSNVs in the sample approached the expected number of SNVs. Contamination at 5% was reliably detected when the contaminant had at least 16 SNVs within the rMLST genes, at 10% with at least 7 SNVs and at 20% with at least 3 SNVs within target genes.

Diversity among rMLST sequences types

To illustrate the genetic diversity within the rMLST scheme, we calculated the numbers of SNVs between ribosomal sequence types (rSTs) for L. monocytogenes, S. enterica, and E. coli (Fig. 2). Over 99 percent of all pairs in all three species had three or more SNVs, which is the cutoff chosen in ConFindr as the minimum number of cSNVs that need to be found before a sample will be considered contaminated. Over 80 percent of the sequence types have 16 or more SNVs relative to others. Therefore, ConFindr should almost always be able to detect contamination between two isolates with different rMLST types.

ConFindr detects contamination with more sensitivity that existing tools

We compared ConFindr to existing tools capable of detecting contamination, CheckM and Kraken using both in silico and lab-generated datasets. Mixes were binned based on the Average Nucleotide Identity (ANI) of the two samples being mixed - those with >99 percent ANI, representing very closely related mixes, those with between 98 and 99 percent ANI, representing same-species mixes between strains not as closely related, and those with ANI of less than 98 percent, representing distantly related same-species mixes and cross-species mixes (Table 1).

ConFindr was more sensitive than either CheckM or Kraken for intraspecies contamination detection (Figs. 3A, 3B, 3D, and 3E), and comparable to both for cross-species contamination detection (Figs. 3C and 3F). ConFindr detected contamination successfully in all cases except for one simulated mix of closely related strains at 5 percent contamination (Table S1), while both CheckM and Kraken required either more distance between species or a higher level of contamination for its determination.

To further verify improved sensitivity of ConFindr relative to CheckM or Kraken, we varied the cutoffs used for each approach. Across our test datasets, the contamination percentage reported by CheckM for uncontaminated samples varied from 0 to 0.6 percent. As shown in Fig. 3, at this cutoff CheckM was not as sensitive as ConFindr. Lowering this cutoff further in an attempt to achieve increased sensitivity resulted in a proliferation of false positives, while still not achieving sensitivity equivalent to ConFindr. For example, lowering the cutoff to 0.3 resulted in 13 of 72 contaminated samples being falsely labelled as clean, with 5 of 16 clean samples being falsely labelled as contaminated. Further lowering the cutoff to classify anything with a reported contamination greater than 0 as contaminated resulted in 7 of 72 contaminated samples being labelled as clean, and 10 of 16 clean samples labelled as contaminated. Similar results were found for Kraken—lowering the cutoff used to a value of 0.05 resulted in false negative calls for 12 of 72 contaminated samples, while false positive results are reported for 10 of 16 clean samples. For the same dataset, ConFindr has only 1 false negative, and 0 false positives.

Assembly metrics are insufficient for contamination detection

We assembled the contaminated datasets used in this study to assess metrics such as number of contigs, N50 and total length for contaminated datasets relative to uncontaminated datasets. At 5% contamination, there was an increase in the number of contigs and the total length of the assembly, and a decrease in N50 (Table 2). The relative increase or decrease in these metrics varied depending on the strain used as the contaminant. For example, contamination of L. monocytogenes strain OLF10129 with OLF15140 appeared to have a smaller impact than contamination with the more distantly related isolate OLF09168. Statistics on all assemblies at all contamination levels are available in Table S1.

Effect of sequencing coverage and platform

To verify that ConFindr does not produce false positives at high sequencing depth and works across both the HiSeq and MiSeq platforms, we created a simulated data set from RefSeq genomes at coverage levels of up to 500X. ConFindr found no false positives on any uncontaminated samples, and only failed to detect contamination in one closely related L. monocytogenes mix (50X HiSeq coverage of NC_002973 and NC_012488 at 5 percent contamination). ConFindr also did not call mixes between two closely related E. coli samples (NC_002695 and NC_011353) as contaminated, but this was due to the fact that the rMLST genes between these two samples contained only 2 SNVs, below the threshold of 3 SNVs that ConFindr uses to call contamination—the 2 SNVs present were accurately detected at many of the coverage and contamination levels tested. The full results for this dataset are available in Table S5.

Contamination in SRA data

To evaluate the prevalence of contamination in public databases, we randomly selected 1,500 isolates sequenced on Illumina instruments from the Sequence Read Archive (SRA) −500 for each of E. coli, S. enterica, and L. monocytogenes. A full list of accessions can be found in Table S4. Of the 1,500 samples examined, 78 (5.27%) were determined to be contaminated by ConFindr (Table 3, Table S4). For all species, intraspecies contamination appeared to be more prevalent than cross-species contamination. Assembly metrics (number of contigs, N50, genome assembly size) for contaminated samples were compared to samples determined to be clean within each species (Fig. S1). Metrics for contaminated samples were generally within ranges observed for samples determined to be clean. As with the simulated data, coverage levels did not impact intraspecies contamination detection. Of the 100 samples with the highest level of coverage (average coverage depth 176 to 6,148, Table S4), only three samples were deemed to have intraspecies contamination. In a comparator set of 100 samples with a coverage depth ranging from 40 to 50.5, five samples were determined to be contaminated. In the set of isolates with a high coverage level, cross-species contamination was detected in three samples (ERR2135583, ERR1100941, SRR1206152) compared to two samples (ERR2520947, SRR7298866) in the set with lower coverage. It is notable that cross-species contamination would have been predicted for the latter two samples based on poor assembly metrics. For samples with higher coverage, assembly metrics were within typical ranges and detection of contamination may be a result of an increased sensitivity for cross-species contamination at this coverage level.

Confirmation of contamination in E. coli WGS data

In a recent WGS analysis performed in our laboratories, hybrid assemblies of E. coli samples led to an incorrect serotype determination (Table 4, sample 1, isolate 3). Three strains of E. coli from two samples were sequenced in duplicate or triplicate. Presumptive contamination was identified due to higher number of contigs or larger genome size relative to duplicates from the same sample (Table 4, bold). In one sample, the serotype of an isolate was incorrectly determined (O159:H2). Contamination was confirmed by analysis of samples with ConFindr (Table 4, bold).

Runtime considerations and installation

ConFindr can be installed with a single command via bioconda (Grüning et al., 2018), and completes analysis on a sample in under one minute when using 4 threads and less than 4 GB of RAM. These features make it practical for ConFindr to be installed and run as a standard quality control step in bioinformatics pipelines.