RpNGS: an automated platform for pathogen identification and monitoring in clinical metagenomics data

RpNGS software overview

RpNGS was written in R programming language (https://www.R-project.org) as a modular Shiny app (Chang et al., 2024), an R package for constructing interactive web applications. The web interface of this application was developed using shiny, shinyFiles, shinycssloaders, shinydashboard, shinyWidgets, and shinyBS libraries. This interface runs in a web browser and provides dynamic changes lab data and clinical metadata through interactive JavaScript tables powered by the DT package. In addition, RpNGS also delivers a set of sophisticated and well-designed interactive visualizations based on the plotly and leaflet package.

Analysis workflow

A typical clinical metagenomic next-generation sequencing bioinformatics pipeline usually comprises quality control, host reads removal, taxonomic categorization and validation (Miller et al., 2019). The RpNGS analysis methodology automates metagenomic processing in a series of replicable stages by click “Confirm and Analyze” button after loading raw FASTQ files of one sequencing run into working directory and fill the related lab data.

Given these processes are conducted by R scripts that connect with external tools within a Conda environment (Fig. S1). The preparation of analysis workflow involves the Conda environment set-up and reference databases. We firstly constructed a Fastp conda environment and installed Fastp, Bowtie2, Kraken2, and Bracken programs before executing analysis pipeline (Lu et al., 2017; Chen et al., 2018; Wood, Lu & Langmead, 2019; Bush et al., 2020). The Bowtie2 index of each microbes was built by the bowtie2-build command with default parameters, while the reference database (PlusPFP.tar.gz) of kraken2 and bracken and human bowtie2 index were download from their websites (https://benlangmead.github.io/aws-indexes/k2 and https://bowtie-bio.sourceforge.net/bowtie2/index.shtml, respectively).

The backend computational pipeline can be conducted within an R environment or delivered on a server for automated processing. The analytic pipeline utilizes Fastp for getting high-quality clean FASTQ file with read length >36 bp. Considering the genomic size disparity between humans and microorganisms, the majority of metagenomic read sets pertained to human DNA, even in samples with a similar number of cells. To better focus on the microbe’s identification, the host reads should be removed before run classification analysis (Bush et al., 2020). Therefore, the cleaned FASTQ file are aligned to the hg37dec_v0.1 reference genome using Bowtie2 with -U, –very-sensitive parameters to obtain the filtered data. Reads categorization and abundance estimate is independently created by Kranken2uniq with –minimum-hit-groups 3 and Bracken with -r 50 -l S against a reference database. After that, we mapping taxonomy list to our pathogen list to acquire the potential pathogens.

In order to exclude the background influences such as microbial nucleic acid residing in reagents, sampling and lab settings, RpNGS construct a z-score for each species. The taxon with z-score >1 was classified as actually detected, which suggests the abundance of species in sample is higher than the controls. The z-score approach was first described in Wilson et al. (2018) and is applied in the IDseq and Pavian metagenomic platform (Breitwieser & Salzberg, 2020) to reflect the significance of relative abundance estimations in a sample as compared to the water controls. The z-score value of each taxon in one sample is calculated as follow: ${\displaystyle z=\frac{x-\mu }{\sigma }}$

where x stands for rpm of taxon in samples, μ is the average of rpm of taxon in control samples, and σ means the of rpm of taxon in control samples.

Reads for candidate pathogen are independently extracted from the filtered FASTQ file of each sample. To validate the pathogen detection results, we run Bowtie2 to map the pathogen specific reads to its matching reference genome for coverage detection, which illustrate by Gviz, Rsamtools and GenomicAlignments R packages in the third tab of RpNGS. The pathogen with greater mapped reads but less coverage will be removed from candidate pathogen list.

Report workflow

R packages flextable and officer were used to generate a pathogen report each sample by exporting a reactive table from shiny into a pre-existing word template.

External benchmarks-datasets and metrics

To evaluate the efficacy of in-house pathogen identification pipeline in RpNGS, we benchmarked it to previously published metagenomic pathogen identification results (Diao et al., 2023). A total of 23 respiratory pathogens at varied concentrations were individually mixed to create 14 microbial mixtures (samples S1–S13 and one negative control), which were then disseminated to 122 laboratories as part of a multicenter mNGS quality assessment study (see Table S1). The raw FASTQ files from these samples were deposited in the Genome Warehouse of the National Genomics Data Center under project PRJCA015554. For benchmarking, we evaluated a total of 28 raw FASTQ files generated by lab005 using IDseq, with the host filter set to “human” and the background set to “none” (see Table S2). The IDseq was selected because it is well-known cloud-based, open-source bioinformatics platform developed for clinical metagenomic study. The sequence analysis procedure of IDseq platform involves three primary steps: host filtering and QC, assembly-based alignment, and taxonomy reporting and visualization. The computed metrics including the precision (P), the recall (R), and the F-score (F1) were calculated based on below formula: P=TP/(TP+FP), R= TP/(TP+FN), and F1=2/((1/R) + (1/P)), where TP are true positives, FP stands for false positives, FN represents false negatives.

Application for pathogen identification

For the comparison with IDseq, we additionally downloaded real clinical metagenomics sequencing files of CHRF0000, CHRF0094 and CHRF0002 from the NCBI under BioProject PRJNA516582 and re-analysis on our in-house system. Samples CHRF0000, CHRF0094 and CHRF0002 in the original study, which contained 91 CSF samples and six water control (Saha et al., 2019), are a water control, Streptococcus pneumoniae, chikungunya virus infection, respectively.

Key functionalities

The summary page (Fig. 1) permits users to visually assess the achieved mNGS test. The interactive visualization choice is a map, pie chart, and bar plot to illustrate the sales volume among locations, proportion of each sample type and samples size distribution between months in one year, respectively. RpNGS employ primary data table to display detail information of processed mNGS test for feasible searching and double checking the clinical reports.

In the RpNGS second tab (Fig. S2), the trained experimenter should update the information of each batch including flow cell ID, sample ID, nucleic acids concentration after extraction and library preparation processes, adaptor ID, and file name of sequencing data. Then start the process step by click the process button. There are six steps inside the process pipeline include detecting and copying the sequencing data from sequencer to server, quality control, host reads removal, and classification, abundance estimate, mapped reads extraction. To assist with distinguishing reads from microbes existing in the reagent and lab’s environment, RpNGS produced z-scores of taxons in each sample by comparison relative abundance estimations of species to water-only or other control sample collections. The progress log will show the status of five critical operations during data analyze.

Multiple detected pathogens are presented in the candidate pathogen list obtained by analysis process (Liu et al., 2021). As a clinical level pathogen detection test, it is important for subset pathogens from the candidate pathogen list depending on their clinical significance. However, identifying whether the detected microbe is the causing pathogen requires a full examination integrating the patient’s symptoms, epidemiological context, and other relevant aspects, that should be done by certified clinician, consular or medical doctor. Results are displaced in the third tab (as in Fig. S3) based on user selection of a number of criteria, including z-scores, mapped reads, average coverage for a particular pathogen, gender, age, sample types, clinically relevant pathogen types, and anti-infection treatment, and interface allow user to select the pathogen for a particular patient from the list of microbes. The user should input the flow cell ID for getting test samples in a status table with three columns, flow cell ID, sample ID and status columns. The clinical information and raw microorganisms list for the sample are activated based on the selection in the report status table. To assist with pathogen detection, the background color of microbes indicated in all putative pathogen candidates extracted from 132 clinical metagenomics reports, generated by Shaanxi Lifegen Co., Ltd, are filled with green. The alignment view and related data of specific species will emerge by choosing one of these colored pathogens in raw microorganism list table. Finally, the pathogen table will be formed by repeated selections in the raw microbe’s table. An important feature is the ability to export experiment data, patient basic information, confirmed pathogens and related aspects of corresponding pathogen in a Word document as clinical pathogen test report for each sample.

Pathogens identification pipeline validation

The presence of diverse amounts and types of bacteria, fungi, DNA viruses, RNA viruses, and human cells in the S1–S12 samples, majority of laboratories performed DNA and RNA extraction and sequencing individually for all samples, including NC and S13. The DNA sequencing data were largely utilized used to calculate precision, recall, and F1-score based on the detection of bacteria, fungi, and DNA viruses in each sample, while the RNA sequencing data were used to evaluate the corresponding metrics for RNA viruses. The detection of low-abundance taxa by IDseq may contribute to a drop in precision compared to the precision acquired by RpNGS (Fig. 2).

We also note that low precision (<0.1) among all instruments that caused by significant number of false positive microbes based on precision formula. As shown in Fig. S4, the recall of DNA sequencing data for samples S11, S12, S3, S4, S6, and S9 analyzed by RpNGS are lower than those acquired by IDseq_NR and IDseq_NT. This mismatch is likely due to Bracken’s default filtering threshold (Lu et al., 2017), which requires the read count to exceed 10. Compared to IDseq, RpNGS did not detect certain pathogens, as shown in Table S3. For example, in sample S11, Kraken2uniq identified only six reads for Haemophilus influenzae, which is comparable with the number of reads allocated to Haemophilus influenzae by IDseq_NT. These results indicate that the threshold of taxonomic reads also critical for DNA-based pathogen detection, special for bacteria. For the RNA sequencing data of samples, RpNGS exhibited an equivalent recall rate for RNA virus detection compared to IDseq.

To further remove false positive caused by factors such as reagents and environmental contamination, we employed Z-score normalization to correct the microbial abundance in the samples. With the removal of contaminant-derived false-positive microorganisms, the precision of RpNGS and IDseq in detecting bacteria, fungus, DNA viruses, and RNA viruses improved while maintaining recall rates (Figs. S5, S6). Additionally, this method greatly boosted the F-score of IDseq comparing it with non-correction (Figs. S7, S8).

Use case: RpNGS for pathogen detection in cases of pediatric meningitis

The RpNGS pipeline was developed for microbial species identification within a human genetic background. However, finding the true causative pathogens of a disease remains problematic due to several influencing factors, including the choice of metagenomic shotgun sequencing analytic methodologies, categorization databases, microbial genetic similarity, and environmental contamination.

To boost pathogen detection, we integrate many components into a single panel, including read categorization tables, clinical metadata, pathogen lists, and associated genome coverage data. To illustrate RpNGS’s capability in clinical pathogen identification, we re-evaluated three cerebrospinal fluid (CSF) samples from a study investigating pediatric meningitis in Bangladesh, which were also analyzed using IDseq.

Figure 3 illustrates the production of a clinical report following data analysis section. Clinical information for sample CHRF0002 is given in an editable table, allowing users to alter patient and sample details (Fig. 3A). Users can pick pathogens from the raw microbial list depending on clinical significance (Fig. 3B). Microbes with the greatest mapping reads in the raw table, indicated in green according to the specified pathogen list, undergo genome coverage visualization to assess read alignment with reference genomes. For instance, the genome coverage of Streptococcus pneumoniae reached 81.52%, increasing the evidence for its pathogenic role (Figs. 3C, 3D).

Additionally, we observed a considerably larger number of reads mapping to Streptococcus pneumoniae in CHRF0002 and Chikungunya virus in CHRF0094 compared to other detected microorganisms (Fig. 4). Based on these data, we identified Streptococcus pneumoniae and Chikungunya virus as the primary pathogens for these two samples, which fits with prior studies (Fig. S9, Table S4) (Saha et al., 2019; Kalantar et al., 2020). These results highlight the amount of categorized reads as critical evidence for infection detection.