RAPiD: a rapid and accurate plant pathogen identification pipeline for on-site nanopore sequencing

Pipeline description

The RAPiD pipeline (Fig. 1) (https://github.com/SteveKnobloch/RAPiD_pipeline) was developed in Nextflow (Di Tommaso et al., 2017) and containerised with Docker (Merkel, 2014). In the first step, basecalled reads are quality filtered using Porechop (Wick, 2017) to remove adapters, followed by processing in NanoFilt and NanoLyse (De Coster et al., 2018) to remove low quality (q < 8) and small (<1,000 bp) reads, as well as reads matching the Lambda phage control genome (GenBank accession J02459.1). Quality filtered reads are then aligned to an indexed reference database (see section on reference database construction) using minimap2 version 2.21 (Li, 2018) with the -ax map-ont flag and without secondary alignments. Alignments are then filtered using SAMtools (Li et al., 2009) to remove supplemental alignments and alignments with a DP alignment score (AS) below 2,000, less than 200 minimizers on the chain (cm) or a gap-compressed per-base sequence divergence (de) above 0.075. These values were determined suitable to filter the majority of non-target sequences at the species level (see results section on cut-off value estimation). Filtered alignments are then summarised by species or subspecies rank and a mean normalised alignment score ${\displaystyle (\frac{AS}{sequencelength}\times 50)}$ and mean per base identity score ( $\left(1-de\right)\times 100$) calculated for each taxonomic group. The results are output into an alignment and report summary. For simplifying interpretation of the report, a confidence value is provided for each taxonomic group based on the number of hits and the mean normalised alignment score, with at least two hits and a score above 75 considered “high”, two hits and a score between 70 and 75 considered “medium” and one hit or a score below or equal 70 considered “low”. Additional features of the web-based application are live base calling using guppy (available from https://community.nanoporetech.com) and automatic expansion of the reference database with user-supplied reference genomes (Fig. 1).

Cut-off value estimation

A reduced reference index size is necessary to improve computational speed and memory requirements of the pipeline, but enables false positive matches due to a lack of off-target reference genomes and presence of highly conserved genes across taxa. To reduce the rate of false positive matches, suitable alignment score cut-off values were determined based on the alignment of publicly available whole genome nanopore sequences, consisting of 14 pairs of closely related species, against a reference index containing only the genome of one species of each pair (see Table S1). Due to differences in sequence quality between datasets, a wide range of species were chosen for which publicly available Nanopore data was available. These included species of the genera Pseudomonas, Burkholderia, Spiroplasma, Lactobacillus, Xanthomonas, Clostridium, Vibrio, Staphylococcus, Puccinia, Fusarium, Phytophthora, Bipolaris, Colletotrichum and Penicillium. Raw genomic reads for each BioSample were downloaded from the NCBI SRA using the sratoolkit (Leinonen, Sugawara & Shumway, 2011). Reads smaller than 1,000 bp were removed with NanoFilt and each resulting dataset normalised to 100,000 reads using the reformat.sh script of the BBMap tool (Bushnell, 2014). A reference index was built in minimap2 using the corresponding representative genomes on NCBI GenBank for one species of each pair. All raw reads were then mapped against the index as detailed above.

Reference database construction and contaminant removal

To construct a non-redundant database for minimap2 alignment within the RAPiD pipeline, representative genomes of 190 plant pathogenic bacteria (13 genera), fungi (76 genera) and oomycetes (eight genera) (Table S2) were downloaded with genome_updater version 0.2.5 (available at https://github.com/pirovc/genome_updater) from NCBI GenBank (accessed on 05.05.2022). Sequences smaller than 2,000 bp were removed. Mitochondrial genomes were removed by performing a BLASTN (Altschul et al., 1990; Camacho et al., 2009) query against the NCBI mitochondrion dataset (https://ftp.ncbi.nlm.nih.gov/refseq/release/mitochondrion/, accessed 15.05.2022). In total, 661 contigs with a percent identity above 80% to a mitochondrial genome, a coverage of at least 500 bp and a subject size below 200,000 bp were removed. Cross-kingdom contamination was checked using Conterminator (Steinegger & Salzberg, 2020). In brief, all-against-all alignments in Conterminator were performed with all representative plant genomes and the representative human genome from NCBI RefSeq (accessed on 01.06.2022). In total, 2,476 sequences that matched a different kingdom than the original accession were removed. Additional, manual inspection of long sequences, that did not get flagged in Conterminator due to their size, was performed using BLASTN against the NCBI nt database. Extensive removal of contaminants was important as representative genomes such as Pseudoperonospora cubensis (GenBank accession AHJF00000000.1) contained sequences larger than 49,000 bp of Cucumis sativus and Acinetobacter calcoaceticus which would otherwise flag as positives in an uninfected sample. A list of all removed sequences is compiled in Table S3. The filtered database was then indexed using minimap2 with the parameter −I 2G. For sub-species identification, the same procedure for index creation was performed, apart from that all genomes corresponding to the 190 plant pathogenic bacteria, fungi and oomycetes were downloaded from NCBI GenBank.

Mock community construction

Two in silico mock communities and one biological mock community were constructed to compare the classification results of the RAPiD pipeline against other available tools. The first in silico community consisted of publicly available nanopore sequences from 15 individual bacteria, fungi and oomycetes (see Table S4) downloaded with sratools from the NCBI SRA database. Each dataset was quality filtered before analysis, using Porechop to remove adapters, NanoFilt to remove reads smaller than 1,000 bp and below a quality score of eight, NanoLyse to remove reads aligning to the Lambda reference genome and BBMap to subset each dataset to 5,000 reads. The second in silico mock community consisted of nanopore sequences from 11 pathovars or strains in the Pseudomonas syringae group (see Table S4) to test for sub-species identification. Quality filtering and subsetting was performed as mentioned above, apart for dataset SRR9943113, which only contained 1,219 reads after quality filtering and dataset SRR10431754, which consisted of four different strains and was thus subset to 20,000 reads. The biological mock community consisted of nine Fusarium species (see Table S4). DNA was extracted from each strain using a modified CTAB-extraction, quantified using a Qubit dsDNA HS assay (Thermo-Fisher Scientific, Waltham, MA, USA) and pooled in equal concentrations by adding 200 ng gDNA of each organism to a tube. A sequence library was constructed using the LSK110 kit (Oxford Nanopore Technologies, Oxford, UK) according to the manufacturer’s instructions and sequenced on a R9.4.1 Flongle flow cell connected to a MinION Mk1C device (Oxford Nanopore Technologies, Oxford, UK). The sequencing run generated 131.17 k reads with an N50 of 14.61 kbp. Raw sequences were basecalled with guppy version 6.1.7 and quality filtering was performed as for the other mock communities. No subsetting of read numbers was conducted, leaving a total of 85,532 reads for analysis.

Comparison of RAPiD to other metagenomic classification tools

The three mock communities were analysed with: (1) the RAPiD pipeline; (2) BugSeq (Fan, Huang & Chorlton, 2021), a web-based and commercial tool; (3) WIMP (Juul et al., 2015), the classification pipeline of Oxford Nanopore Technologies which uses centrifuge (Kim et al., 2016) as classification engine; (4) Kraken2 version 2.1.2 (Wood, Lu & Langmead, 2019), a short-read classifier; and (5) an approach using minimap2 alignment against the entire NCBI nt database followed by MEGAN (MEGAN6 Community Edition) lowest common ancestor calculation (Huson et al., 2016). For analysis with BugSeq (accessed February 2022), each mock community was uploaded and analysed against the NCBI nt database. Classification with WIMP was performed by uploading each mock community through the EPI2ME agent (Oxford Nanopore Technologies, Oxford, UK) and choosing Fastq WIMP analysis against the database v2021.11.26. Classification with Kraken2 was performed with default parameters against the precompiled database PlusPF and the memory reduced database PlusPF-8 retrieved from https://benlangmead.github.io/aws-indexes/k2 (accessed July 2023). For the minimap2-MEGAN pipeline, first the NCBI nt database was downloaded (http://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz, accessed February 2022) and indexed using minimap2 with default parameters. Alignment of each mock community was then performed against the index with minimap2 and the -ax map-ont flag. The alignment file was subsequently subject to lowest common ancestor (LCA) calculation using the MEGAN sam2rma script with the -lg -alg longReads flags against the MEGAN mapping file megan-nucl-Jan201.db (https://software-ab.cs.uni-tuebingen.de/download/megan6/welcome.html, accessed January 2022).

Read tables were then summarised by counting reads at species or sub-species level. Reads that matched a target species or sub-species present in the mock community were classified as true positives. In contrast, reads that matched a species or sub-species not present in the mock community were classified as false positives. For each mock community, the precision, recall and F1 score was calculated along with the number of true and false taxa detected with each tool. Precision was calculated as the ratio between true positive reads and all reads classified at the species or sub-species level. Recall was calculated as the ratio between true positive reads and all reads analysed. The F1-score was calculated as $2\ast \left(precision\ast recall\right)/\left(precision+recall\right)$. For the Pseudomonas syringae mock community, reads not assigned to the correct pathovar or strain were considered false positives.

The Kraken2 pipeline with the PlusPF database and the minimap2-MEGAN approach were run on one node of an HPC (five nodes, each with 64 cores (AMD EPYC 7301) and 1 TB RAM). The RAPiD pipeline and Kraken2 with the PlusPF-8 database were run on a laptop (Dell Latitude 5420, Intel Core i5-1135G7 and 16 GB RAM). With the exception of the web-based tools, runtime, CPU-time and maximum memory requirements were recorded for each analysis.

In-field sample collection, DNA extraction and sequencing

The following protocol was developed to rapidly process samples for sequencing outside of the laboratory without the use of dangerous chemicals in the DNA extraction procedure or specialised equipment, such as centrifuges, thermocyclers and heat blocks. The set-up for in-field sample processing and sequencing is depicted in Fig. 2. Samples of wheat, tomato and apple leaves showing signs of infection were collected and processed individually as follows. Approximately 100 mg of sample was cut into small pieces with a sterile scalpel, transferred to a 1.5 ml Eppendorf tube and immediately macerated in 800 µl of CTAB-Buffer (2% CTAB, 2% PVP, 100 mM Tris, 1.4M NaCl, 20 mM Na₂EDTA, pH 8) with 14 µl of 20 g l⁻¹ Proteinase K using a hand-held pestle grinder (DWK Life Sciences Kontes). The sample was then kept at approximately 65 °C for 15 min using a thermal flask and glass thermometer. The lysate was then subject to DNA purification with a S-TECH device (DNAiTECH) according to the manufacturer’s instructions for soil samples and resuspended in 140 µl of elution buffer (10 mM Tris-HCl pH 8.5). The resuspended DNA was then normalised to approximately 500 ng in 23.5 µl elution buffer with a pre-made, modified concentration of MAG-BIND TotalPure NGS magnetic beads (Omega Bio-tek). In short, 75 µl of MAG-BIND TotalPure NGS magnetic beads was thoroughly resuspended and 3.5 µl transferred to a new tube. The remaining reagent was spun-down to pellet the beads and 66.5 µl of the supernatant was added to the previously aliquoted 3.5 µl of beads. This resulted in 70 ul of binding reagent with a reduced number of beads for a 0.5x clean-up and concentration of the DNA. In previous experiments, this concentration of beads bound and subsequently released approximately 400–600 ng of DNA in 23.5 µl of elution buffer from various plant samples (data not shown). All subsequent clean-up steps were performed according to the manufacturer’s instructions. Sequence libraries were constructed using the LSK114 kit according to the manufacturer’s instructions with omission of the first magnetic bead clean-up step. The temperature-controlled steps were achieved by adding hot or cold water to a thermal flask with a glass thermometer. Libraries were sequenced on R10.4.1 Flongle flow cells on a battery-operated MinION Mk1C device (Oxford Nanopore Technologies, Oxford, UK). Basecalling was performed in simplex high accuracy (HAC) mode. RAPiD was run on a laptop with 16 Gb RAM. The workflow from sample to DNA library took less than 2 h per sample. The total runtime for each sequencing run was dependent on the battery pack and was approximately 5.5 h. The DNA concentration for each sample after DNA extraction was analysed in the laboratory after the field trial using a Qubit dsDNA HS assay and is listed in Table S5.

Data availability

Raw sequence data of the Fusarium spp. mock community and the three infected plant datasets sequenced in the field are stored in the NCBI SRA under BioProject PRJNA1063731.

Cut-off value estimation

The values for alignment metrics of target and non-target species against the reference database are shown in Fig. 3. The metrics most discriminative of a true positive match were the DP alignment score (AS), number of minimizers on the chain (cm) and the gap-compressed per-base sequence divergence (de). The chaining score (s1) largely overlapped with cm. Despite not having a reference genome in the database, each non-target species matched a false reference genome. Different cut-off values were compared to exclude most false positives while retaining most true positive reads. Different cut-off values were investigated to find a setting that retained the majority of target reads while excluding over 99% of off-target reads. Choosing values of 2,000, 200 and 0.075 for AS, cm and de, respectively, included 57.98 ± 23.95% of all mapped target reads while excluding 99.18 ± 1.93% of all off-target reads. Upon closer examination, many remaining false positive reads matched conserved genes of the mitochondrial genome, the ribosomal RNA genes or internal transcribed spacer regions. Other alignments metrics such as mapping quality (MAPQ) and total number of mismatches and gaps (NM) were less suitable for species-level discrimination (Fig. 3).

Mock community analysis

The RAPiD pipeline had the highest precision for each of the mock communities, with 0.99, 0.94 and 0.98 for the mixed mock, P. syringae mock and Fusarium mock, respectively (Table 1). In addition, it detected all taxa in the mock correctly. High precision, however, came at the cost of lower recall of 0.50, 0.50 and 0.37, respectively. RAPiD had the highest F1 scores of 0.65 for the P. syringae mock and 0.53 for the Fusarium mock. For the mixed mock, RAPiD had a lower F1 score than all pipelines apart from WIMP. BugSeq correctly assigned all 15 species of the mixed mock with a precision of 0.94 and recall of 0.80. However, it did not detect five pathovars of the P. syringae mock and seven species of the Fusarium mock, leading to a lower precision of 0.33 and 0.44, respectively. WIMP detected seven of the species in the mixed mock, six pathovars in the P. syringae mock and none of the correct species in the Fusarium mock. The precision was 0.68, 0.51 and 0 for each of the mock communities, respectively. The Kraken2 pipeline performed similarly between both the precompiled, large and memory-reduced databases, detecting 10 species in the mixed mock and seven pathovars in the P. syringae mock, leading to a precision of between 0.26 and 0.88 and recall of 0.22 and 0.61, respectively. Similar to WIMP, none of the Fusarium species were detected with Kraken2. The minimap2-MEGAN pipeline performed best on the mixed mock community, detecting all species with a precision of 0.95. On the P. syringae and Fusarium mocks it detected 6 and 5 taxa with a precision of 0.24 and 0.31, respectively.

The number of taxa wrongly assigned was highest for WIMP and Kraken2, which assigned reads to over 1,000 species not in the mock community. The pipeline with the lowest number of wrongly assigned taxa was the minimap2-MEGAN approach with 7 to 25 taxa depending on the mock community, followed by the RAPiD pipeline with 9 to 82 wrongly assigned taxa. Runtime was shortest for Kraken2 with the 8GB database, followed by Kraken2 with the full reference database and then WIMP, RAPiD and BugSeq. The memory requirements were lowest for Kraken2 with the 8GB database, followed by RAPiD, the minimap2-MEGAN approach and Kraken2 with the full reference database.

Pathogen detection in the field

Three infected plant samples were processed and their DNA sequenced in the field. These included a leaf of a common wheat plant with wheat leaf rust (Puccinia triticina), a leaf of a tomato plant with leaf blight (Alternaria spp.) and a leaf of an apple tree with apple scab (Venturia inaequalis) (Fig. 2). The instrument set-up for DNA extraction and sequencing is shown in Fig. 2. DNA concentrations of the three plant samples ranged from 4.5 to 12.2 ng µl⁻¹ (Table S5). The time from sample collection to starting the sequencing run took less than 2 h. Sequencing generated 68.7 to 192.6 Mbp of raw sequences or 37,565 to 117,610 reads per sample (Table S5). After sequencing commenced, the first sequence match for P. triticina on the wheat sample was detected within 30 min. In total, 303 sequences were assigned to P. triticina out of 16,027 reads that passed quality filtering and the average normalized alignment score for the matches was 77.48 with a mean per base identity of 95.37% (Fig. 2). For the tomato sample, two hits matched Alternaria alternata out of 27,853 quality filtered reads. The normalized alignment score for these matches was 89.44 with a mean per base identity of 97.60%. The first hit was detected within 3 h of sequencing. For the apple sample, three sequences were assigned to V. inaequalis out of 46,039 reads after quality filtering. The average normalized alignment score was 87.19 with a mean per base identity of 96.63%. The first hit was detected within 1.5 h after sequencing commenced. No false positive hits were recorded for any of the three samples. A BLASTN search of all matching sequences against the nt database (accessed 01.07.24) confirmed the presence of Puccinia triticina and Alternaria alternata. However, the apple sample sequences showed closer sequence similarity to Venturia effusa, albeit with low percent sequence identity, suggesting a different pathogen currently not present in the RAPiD database. All sequence hits are listed in Table S6.