A Sample-to-Report Solution for Taxonomic Identification of Cultured Bacteria in the Clinical Setting Based on Nanopore Sequencing

Bacterial isolates all originated from the Institute for Infectious Diseases (IFIK, Bern) Biobank. The IFIK provides the entire spectrum of medical microbiological diagnostic services to the largest Swiss hospital group (Inselgruppe) and other regional hospitals. The diagnostic division of IFIK (clinical microbiology) is ISO/IEC 17025 accredited to perform routine bacterial diagnostics from clinical samples. ATCC strains were obtained from LGC Standards (Wesel, Germany) and were grown on solid medium as recommended by the manufacturer.

Bacterial cultures grown overnight were harvested from agar plates and dissolved in 300 μl of Tris-EDTA (pH 8.0). DNA was extracted with a NucliSens easyMag (bioMérieux, Switzerland) robot according to the manufacturer’s protocol. The 16S rRNA gene PCR was performed with the primer sets 16S_f, 5′-AGAGTTTGATCMTGGCTCAG-3′, and 16S_r, 5′-TACCGCGGCWGCTGGCACRDA-3′, (general bacteria) and mbak_f, 5′-GAGTTTGATCCTGGCTCAGGA-3′, and mbak_r, 5′-TGCACACAGGCCACAAGGGA-3′, (mycobacteria) supplemented with the universal tails 5′-TTTCTGTTGGTGCTGATATTGC-3′ (ONT forward primer), 5′-ACTTGCCTGTCGCTCTATCTTC-3′ (ONT reverse primer), 5′-TGTAAAACGACGGCCAG-3′ (M13f, Sanger forward primer), or 5′-CAGGAAACAGCTATGAC-3′ (M13r, Sanger reverse primer). PCR mixtures (25 μl) for general bacteria and mycobacteria were assembled, respectively, with 1 and 2.5 ng DNA template and 10 μl of a 1.25 and 2.5 μM primer working solution, both with 12.5 μl Q5 master mix. Amplification was performed in a GeneAmp 9700 thermocycler (Thermo Fisher Scientific, Inc., MA, USA) with the following program: 98°C for 1 min; 30 cycles of 98°C for 10 s, 63°C for 15 s, 72°C for 30 s; and 72°C for 2 min. PCR products were purified with CleanNGS beads (CleanNA, Waddinxveen, Netherlands) according to the manufacturer’s instructions with the following modifications: after the washing step, an additional 3-s centrifugation step was introduced, and the purified DNA was eluted in 80 μl of Tris-HCl (0.01 M, pH 8.0). Fragment size of the amplicons was analyzed using the TapeStation D1000 assay (Agilent, Santa Clara, CA, USA), concentrations were measured with the Qubit double-stranded DNA (dsDNA) broad-range (BR) assay (Thermo Fisher Scientific), and the purity of the DNA was analyzed with a NanoDrop spectrophotometer (Thermo Fisher Scientific). Samples with DNA concentrations of <1.05 nM were excluded from the analysis.

A typical library consisted of the pooling of PCR amplicons from 2 to 15 clinical samples and 1 positive control (Mycobacteria intracellulare, amplified with general bacterial primers). Library preparation was performed with the kits EXP-PBC096 and SQL-LSK109 (Oxford Nanopore Technologies, Oxford, UK) using the supplementary reagents NEBNext end repair/dA-tailing module (E7546; New England Biolabs, ON, CA), NEB Blunt/TA ligase master mix (M0367; New England Biolabs), Taq 2× master mix (NEB M0270; New England Biolabs), and CleanNGS beads (CleanNA). All modifications made to the manufacturer’s protocol (PCR barcoding [96] genomic DNA, PBAC96_9069_v109_revK_14Aug2019) are described in the following section (see also Fig. 1A; for a detailed protocol, see Text S1 in the supplemental material). AMPure beads were substituted with CleanNGS beads, and the HulaMixer (Thermo Fisher Scientific) parameters “orbital: 40 rpm, 07 s; reciprocal: 89 deg, 2 s; vibro: 5 deg, 2 s; vertical position” were used. Barcoding PCRs (12 cycles) were set up with 25.2 nmol of template per reaction. Raw barcoded PCR products were quantified with the Qubit dsDNA BR assay and pooled at equal molar proportions. Products containing less than 0.57 pmol DNA were excluded from the analysis. If the total amount of DNA in a pooled library was below 9.23 pmol, “place-holder” (filling) barcoded samples were added to the pooled library to avoid flow cell underloading (see example of calculations and adjustments in Text S1). Place-holder barcoded samples were produced in advance from the same template as the positive controls, with 15 instead of 12 barcoding PCR cycles. Resulting PCR products were quantified with Qubit and stored at −20°C. The pooled library was purified (CleanNGS beads; 50-μl elution volume) and quantified with the Qubit dsDNA BR assay. The purified library pools were diluted to 140 nM before proceeding to the “end preparation” step of the protocol.

ONT sequencing was performed on a GridION X5 instrument (Oxford Nanopore Technologies) with real-time basecalling enabled (ont-guppy-for-gridion v.1.4.3-1 and v.3.0.3-1; fast basecalling mode). Sequencing runs were terminated after production of 1 million reads or when sequencing rates dropped below 20 reads per second. Purified PCR products were submitted to Sanger sequencing at Microsynth (Balgach, Switzerland).

(i) LORCAN pipeline description. LORCAN was developed to facilitate reproducible ONT sequencing-based marker gene analysis in diagnostics facilities. The pipeline, written in Perl 5, R, and BASH, runs on Linux servers or workstations. The code is publicly available (19) and is based on publicly available, third-party dependencies (see Table S1 in the supplemental material). Major steps of the workflow are described in the following section (numbers correspond to the steps in Fig. 1B). In step 1, basecalled reads are demultiplexed and adapters trimmed (Porechop [20], parameters: –format fasta, –discard_unassigned, –require_two_barcodes). In step 2, reads are filtered by length, keeping only those with lengths of −20 to +100 bases (lower boundary adjustable) around the modal sequence length (custom Perl and R scripts) (Fig. 1B). In step 3, reads are mapped to a nonredundant reference database (minimap2 [21]; see database preparation below). In step 4, reads are extracted, binned by taxonomic level (here species), and remapped to the reference sequence that obtained the highest number of mapped reads among all sequences of the corresponding species (minimap2, SAMtools [22], SeqKit [23]). In step 5, consensus sequences are derived using a 50% majority rule consensus. In step 6, the 10 closest reference sequences are selected by sequence similarity to the consensus sequence (blastn, BLAST+ [24]). In step 7, phylogenetic trees for each consensus sequence with its 10 closest references are created (MAFFT [25] with parameters -maxiterate 1000 –localpair; gBlocks [26] with parameters −t = d; and IQ-Tree [27] with parameters -m GTR+I+G -bb 1000 -czb). Parameters of all software are also provided in the LORCAN GitHub repository.

(ii) Database preparation. Reference databases used by LORCAN are nonredundant and assay specific. Detailed instructions for database creation are provided online (https://github.com/aramette/LORCAN/). In short, the reference database (in this study, leBIBI SSU-rDNA-mk37_stringent, https://umr5558-bibiserv.univ-lyon1.fr/BIBIDOCNEW/db-BIBI.html [28]) was trimmed to the region of interest (amplified region minus primers) and dereplicated (mothur [29]), and sequence names were simplified (custom Perl scripts). The names of identical sequences are saved to a file during the dereplication step. The resulting nonredundant database is then used to generate a custom BLAST database, which is used in LORCAN pipeline.

(iii) Sanger sequence analyses. Forward and reverse sequences were assembled into consensus sequences using SeqMan Pro (DNA Star, Madison, WI, USA), primers were trimmed manually, and ambiguous bases were resolved based on visual inspection of the chromatograms. Consensus sequences were taxonomically classified using the online tool leBIBI QBPP (28, 30).

(iv) SNV discrimination and performance with mixed samples. Amplicons produced from pure samples were quantified (Qubit dsDNA BR assay). Mixtures of pure amplicons were produced at defined ratios before library preparation to produce libraries of heterogeneous (“mixed”) samples. Artificial read mixtures were also produced in silico by mixing reads originating from pure amplicon samples. Those reads were obtained from the LORCAN output directories (output file 1_fasta/BC*.fasta produced by step 2) (Fig. 1B) and sampled using Seqtk subseq (v.1.3-r106) (https://github.com/lh3/seqtk) to produce different proportions of original, pure amplicons. Reads from mixed amplicon samples were fed back into LORCAN, and detected species compositions were extracted from the resulting LORCAN reports. Sequence identities between the paired Mycobacterium species were determined based on pairwise alignment of the amplified region using MultAlin v.5.4.1 (http://multalin.toulouse.inra.fr/multalin/ [31]).

(v) Influence of database completeness on consensus accuracy. Amplicons from a set of seven ATCC reference strains were ONT sequenced and analyzed with LORCAN using the full nonredundant leBIBI 16S rRNA database, restricted to the region amplified by the general bacterial primer set. The resulting top consensus sequences were extracted and combined with the above-mentioned database. The resulting sequence data set was aligned (MAFFT v.7.313, FFT-NS-1, progressive method), and pairwise distances were calculated (mothur v.1.40.5, dist.seqs, calc = eachgap, countends = F, cutoff = 0.20). For each consensus sequence, 10 subsets of sequences with minimal distances below thresholds ranging from 0 to 0.1 were extracted (Seqtk subseq), and minimal distances between each data set and the corresponding consensus sequence were analyzed. The seven ATCC read sets were reanalyzed with LORCAN and the corresponding database subsets to produce consensus sequences. Top consensus sequences from each combination of sample and subsetted database were extracted, combined with the consensus sequences generated with the full database, and aligned (MAFFT v.7.313, L-INS-I, iterative refinement method (<16) with local pairwise alignment information). Pairwise distances were analyzed as described above, and distances between the consensus sequences generated with the full and the subsetted databases were extracted.

All reads and consensus sequences corresponding to the data presented in Table 1 and the LORCAN-derived consensus sequences used as references in Fig. 3 were deposited in the European Nucleotide Archive under project accession number PRJEB34167 or made available as supplementary multi-FASTA files.

To test the ability of LORCAN to resolve mixed samples, artificial mixtures were created by mixing either amplicons (Fig. 2A) or reads produced from pure samples (Fig. 2B and C; see also Fig. S1 and S2 in the supplemental material). The taxonomic identity of all involved strains was successfully recovered by LORCAN. The slightly lower amplicon length of Pseudomonas aeruginosa compared to those of Staphylococcus aureus and Enterococcus faecalis resulted in a slight underrepresentation of the latter in the mixtures (Fig. 2B) due to the narrow size window chosen for read size selection (the lower boundary of the size window around the modal read length is adjustable in the LORCAN command line). The mixture of two Mycobacterium species (97.6% sequence identity in the amplified region) (Fig. 2C) was accurately reproduced.

We analyzed the influence of reference database completeness on the resulting consensus quality and accuracy by creating incomplete reference databases, from which we excluded reference sequences if they were too close to the ideal reference sequence, and then performed LORCAN analysis with each of these truncated databases in turn. The genetic distances of the closest reference sequences in the reference database strongly influenced the accuracy of the resulting consensus sequences. For instance, Enterococcus faecalis showed the lowest consensus accuracy at 95% database identity (Fig. 3). This was caused by gaps in the closest reference sequence available. For databases with closest identities of ≤94%, the reference sequence with the identified gaps was absent and consensus quality increased again (see Fig. S3 and S4 in the supplemental material). Classification at the species level was, however, virtually unaffected in pure amplicons. The Eggerthella lenta data set contained a contamination of Pseudomonas stutzeri reads (0.8% of all reads), which did not influence classification when reference sequences that enabled the mapping of Eggerthella lenta reads were available. In the absence of sufficiently close reference sequences, the sample was misidentified (Fig. 3A). Information provided in the LORCAN report did, however, reveal that the Pseudomonas stutzeri consensus sequence was only based on 20 out of 850 reads, which therefore indicated a likely case of suboptimal taxonomic classification.

The comparison of 78 LORCAN-generated consensus sequences from 14 sequencing runs (including 49 clinical samples and 15 ATCC reference strains) to their corresponding Sanger sequences revealed an average sequence identity of 99.6% ± 0.6 (standard deviation). The positive control (originating from the same pool of amplicons) that was systematically sequenced in these 14 runs showed an average identity of 99.8% ± 0.2 to its corresponding Sanger sequence. All reference strains were correctly identified at the species level by LORCAN. Identification by leBIBI QBPP resulted in assignment of the expected species (lowest patristic distance) or the placement of the expected species in the proximal cluster of the query sequence (in the phylogenetic tree) in all but two cases. In these cases, the analyzed strains were placed in close neighborhood of the expected species in the phylogenetic tree produced by leBIBI QBPP (Table 1; see also Fig. S5 in the supplemental material).

Costs per sample of the Sanger method were the lowest across different sequencing technologies (Fig. 4), provided the analyzed amplicons are pure and short enough to be covered by a single sequence at sufficient quality. Among the analyzed next-generation sequencing (NGS) methods, nanopore sequencing was by far the most cost-effective option particularly at throughputs of 24 to 48 samples. The high costs per sample for Illumina are mainly caused by the nonreusable sequencing cartridges (the full costs apply regardless of the number of processed samples) and the comparably high prices of the library preparation kits.

We studied the influence of the read size fraction (relative to the modal read length) and the number of input reads on LORCAN consensus quality. In short, optimal results were obtained when reads shorter than 20 bases below the modal read length were excluded from the analysis (see Fig. S6 in the supplemental material). Further, we found 100 reads to be sufficient for the generation of high-quality consensus sequences (see Fig. S7, S8, and S9 in the supplemental material). The required number of input reads may vary with the taxonomic complexity of the analyzed samples and the resolution required by the operator. From a theoretical viewpoint (Fig. 1B, step 2), a total of 3,000 size-selected reads may allow for the creation of high-quality consensus sequences and reliable species identification for species contributing ≥3.3% of those 3,000 selected reads (i.e., when setting a minimum reference mapping depth of 100 reads in LORCAN, which corresponds to the minimum number of reads recommended for reliable consensus creation [see Fig. S7]). In most cases, however, even when a sample may consist of amplicons derived from a unique species, not all reads are assigned to the target species (e.g., due to read errors and/or the presence of highly similar sequences associated with other species in the reference database). Furthermore, demultiplexing and size selection could result in significant reduction of available reads. For illustrative purposes, during our last 11 sequencing runs consisting of 89 samples (including place-holder samples; see Library preparation in Materials and Methods), an average of 639,944 ± 267,704 basecalled reads were produced, while multiplexing on average 8 ± 3 barcoded samples per sequencing run. Read demultiplexing produced thereafter an average of 46,571 ± 22,129 reads per library (i.e., in 58% of all reads, both index sequences have been identified and assigned to the same barcode). This comparably high read loss resulted from the stringent demultiplexing parameters used (detection of both 5′ and 3′ barcodes required, exclusion of reads with internal barcodes), which may effectively prevent cross talk between libraries (32). Subsequent size selection (read length of −20 to +100 bp around the modal sequence length) resulted in an average of 43,265 ± 21,305 reads per barcode, which were available for further processing. Samples producing more than 3,000 reads of the expected amplicon size were further down-sampled (adjustable LORCAN parameter), resulting in an average number of used reads of 3,008 ± 6 reads per sample. All samples, controls, and place holders processed in these 11 sequencing runs were successfully taxonomically identified. Although species identification could have been achieved with a lower number of reads per sample, sequence production was fast (i.e., approximately 1 to 2 h for 1 million reads), and even if flow cells may have been reused up to four times, the maximal sequencing capacity of the flow cells was never utilized (see Table S2 in the supplemental material).