Assessment of metrics in next-generation sequencing experiments for use in core-genome multilocus sequence type

Introduction

With the reduction in the cost of next-generation sequencing (NGS), whole-genome sequencing (WGS)–based methods are being widely used in genomic epidemiology to characterize bacterial pathogens and perform strain typing (Deng, Bakker & Hendriksen, 2016; Fratamico et al., 2016; Lindsey et al., 2016). Multilocus sequence type (MLST) genotyping (Maiden et al., 1998) has been used for many years for cross-laboratory comparison and outbreak investigation among closely-related strains. Core-genome MLST (cgMLST), an advanced version of MLST genotyping, is a genome-wide gene-by-gene comparison approach (Maiden et al., 2013) that has been successfully used for detecting disease clusters and investigating outbreaks (Barkley, Gosciminski & Miller, 2016; De Been et al., 2015; Jackson et al., 2016). Several websites and databases, such as PubMLST.org (Jolley, Bray & Maiden, 2018) and Pathogenwatch (https://pathogen.watch/), that are funded by large companies and governments have been using cgMLST. Because of the increasing significance of cgMLST in the field of epidemiology, evaluating its robustness is crucial. Segerman (2020) has reviewed sequencing technologies and assembly methods for the bacterial surveillance and the RefSeq Genome Database. He found that Illumina sequencers were the mostly used sequencing platforms and SPAdes, SKESA and CLC were the most popular assemblers. Based on Segerman’s (2020) findings, we designed a metric “number of core genes unrecalled” to find out the minimum sequencing depth/coverage for SPAdes, SKESA and CLC at read lengths with 150 bp and 250 bp, which were common in Illumina platforms, to recover the most completely “core gene alleles” (i.e., not only gene locus but also nucleotide sequence of such gene locus needed to be the same). The idea of metric “core gene unrecalled” was from the benchmarking metrics of genome assemblies (i.e., Contiguity, Correctness and Completeness) suggested by Molina-Mora et al. (2020) with the scale from genome level down to gene level. Also, because the genes order within a genome does not influence the generated cgMLST profile, we only consider the correctness and completeness of core genes. Therefore, our designed metric “number of core genes unrecalled” could fully reflect the quality of cgMLST profiles. Since the sequencing read length, read depth, and assembler might substantially affect cgMLST results, we investigated the effect of these factors on cgMLST results. In this study, we simulated different read depths of different lengths from four common food-borne pathogens, namely, Escherichia coli, Listeria monocytogenes, and Salmonella enterica, and performed assembling by using different assemblers to determine the minimum read depths required under different situations (i.e., different combinations of read lengths, read depths, and assemblers). The minimum read depths determined in this study might help researchers in estimating the depths before conducting cgMLST studies.

Methods and Materials

To evaluate the minimum read depth required for recalling genes, we simulated read sets with different read depths from complete reference genomes downloaded from NCBI. Three food-borne pathogens were tested: E. coli, L. monocytogenes, and S. enterica. Different assemblers and read lengths were included in the evaluation. The experiments were repeated three times to ensure the robustness of the results.

Results

We evaluated 14 sequencing depths or coverages (10 ×, 20 ×, 30 ×, 40 ×, 50 ×, 60 ×, 70 ×, 80 ×, 90 ×, 100 ×, 200 ×, 300 ×, 400 ×, and 500 ×) for determining the assembly quality. The number of unrecalled genes from the reference genomes of S. enterica LT2, E. coli IAI39, and L. monocytogenes EGD-e represented the assembly quality. Three commonly used assemblers, namely SPAdes, CLC, and SKESA, were applied to run the tests. Two read lengths, 150 bp and 250 bp, representing the widely used sequencing read lengths in Illumina HiSeq and Illumina MiSeq platforms, respectively, were evaluated.

Minimum sequencing depth achieving stable number of core genes unrecalled for different assemblers by using different read lengths

As shown in Fig. 1, a sequencing coverage of 60 × might be a safe choice irrespective of the assembler and read length. We observed that the SPAdes assembler required 30 × read depths (irrespective of whether the read length of 150 or 250 bp was used) to achieve minimum depth of the stable core genes unrecalled compared with CLC and SKESA that required read depths of 40 × ∼60 ×. Because the read lengths of 150 and 250 bp are mainly used in Illumina platforms, we evaluated the minimum sequencing coverage required for these two read lengths. As shown in Fig. 1, sequencing coverages of at least 60 × and 50 × were required for the read lengths of 150 and 250 bp, respectively, for assembly to achieve the stable core genes unrecalled irrespective of the assemblers used (i.e., SPAdes, CLC, and SKESA). Regarding assemblers, we observed that SPAdes was not considerably affected by the read length and required a read depth of only 30 × to recover reference genes. However, CLC and SKESA required a read depth of at least 40 ×–50 × and 50 ×–60 ×, respectively, to achieve assembly quality similar to that obtained using SPAdes.

Discussion

In our evaluation, we applied an index “genes called” to represent the assembly quality. Because read sets were simulated from a complete genome, the number of “genes called” can directly represent the quality of assemblies. Because the number of genes called can indicate the completeness of a “pan genome”, which covers “core genes”, we used genes called as our evaluation index. The three most important factors in NGS were evaluated: read depth, read length, and assembler. We found that an assembler was the most crucial factor that affected the quality of assemblies, especially at a low read depth. For low read depths (20 ×–30 ×), SPAdes outperformed CLC and SKESA with an error rate of <2.0%, although the performance of CLC was close to that of SPAdes. Compared with SPAdes and CLC, SKESA usually required 40 ×–50 × to reach an error rate of <2.0%. Although SPAdes demonstrated the highest performance, its running time considerably increased with the read depth, especially when the depth was >100 ×. No difference in results was observed between long reads (250 bp) and short reads (150 bp) for SPAdes and CLC; however, SKESA required a larger depth to assemble short reads to reach an error rate of <2.0%. In addition, to investigate if some similarity sharing among unrecalled genes at even high sequencing depth, we analyzed the unrecalled core genes of S. enterica, E. coli and L. monocytogenes at depth 100 ×, and no commonality among these genes was found. The unrecalled core genes at depth 100 × from Table 1 are listed in Table S2 (the union of the triple repeat are listed). In summary, we recommend sequencing at a read depth of 30 ×–50 × and a read length of 250 bp by using SPAdes as the assembler to maintain a balance of cost/pay ratio in cgMLST.

Supplemental Information