Nanopore Metagenomic Sequencing for Detection and Characterization of SARS-CoV-2 in Clinical Samples

Study population and specimen collection

Clinical specimens were collected from two different populations. First, oropharyngeal swabs were collected in 2 mL of a guanidinium-based inactivation agent (Prestige Diagnostics) as part of a study conducted to estimate SARS-CoV-2 infection prevalence in a UK community from April to June 2020. Swab samples from 2714 individuals from around the greater Oxford area were collected to compare PCR, serology, and Nanopore sequencing for SARS-CoV-2 infected versus uninfected subjects. A set of eight SARS-CoV-2 PCR positives or indeterminate samples from this population were included in the current study. Second, nasopharyngeal swab specimens collected in 3 mL viral transport medium (Yocon Bio-technology Co. Ltd) were obtained from routine SARS-CoV-2 community testing at Vancouver General Hospital (VGH) or the BC Centre for Disease Control (BCCDC) (Vancouver, British Columbia, Canada) (n = 35). RT-PCR testing for COVID-19 was performed for all samples at either the BCCDC Public Health Laboratory or the medical microbiology laboratory at VGH using either the Roche MagNA Pure extraction system (Roche Diagnostics, Laval, Canada) in combination with detection of E-gene and RdRp gene targets, or the Panther Fusion SARS-CoV-2 assay (Hologic Inc., San Diego, CA) detecting two targets in ORF1ab. Primers for the SARS-CoV-2 RNA-dependent RNA polymerase (RdRp) were developed in-house by the BCCDC Public health laboratory and primers for the E gene were based on the World Health Organization RT-qPCR protocol (3). The human RNaseP gene was used as an internal control as suggested by the World Health Organization (https://www.who.int/csr/resources/publications/swineflu/CDCRealtimeRTPCR_SwineH1Assay-2009_20090430.pdf?ua=1). A table containing primers and probes used for these assays can be found in Supplementary Table 1. Additionally, PCR screening for potential variants of concern (VOCs) (Ex. Alpha, Beta, Gamma, Delta variants) was performed on 11 of the positive swabs obtained from VGH that were collected during May 2021. Primers and probes were designed to target the N501Y and E484K mutations (Supplementary Table 2). Swabs were stored at either -80°C for the oropharyngeal swabs or - 20°C for the nasopharyngeal swabs. Specimens were chosen to obtain test performance metrics for Nanopore mNGS across a range of C_t values (Supplementary Figure 1).

RNA Extractions

Prior to extraction, samples were vortexed and 200 uL of each sample was centrifuged at 16,000g for 3 minutes to pellet host cells. 140 uL of supernatant was aspirated and viral RNA was extracted from the supernatant using the QIAmp Viral RNA kit (Qiagen) as previously described (11), and eluted in 30 uL nuclease-free water. Samples were treated with TURBO DNase (Thermo Fisher Scientific) and incubated at 37°C for 30 minutes, followed by concentration and clean-up with the RNA Clean & Concentrator-5 kit (Zymo Research); finally, eluting in 8 uL nuclease-free water.

SISPA Amplification

SISPA amplification was performed as described previously (9-13). Briefly, concentrated RNA was incubated with primer A (100pmol/uL; 5’ - GTTTCCCACTGGAGGATA(N₉) - 3’) and then reverse transcribed using SuperScript IV Reverse Transcriptase (Thermo Fisher Scientific). Second strand synthesis was performed using Sequenase Version 2.0 (Thermo Fisher Scientific), following which, RNase H was performed to digest any remaining RNA. Random amplification was performed on each using AccuTaq LA DNA polymerase (Thermo Fisher Scientific) and SISPA primer B (5’ - GTTTCCCACTGGAGGATA - 3’). This reaction underwent PCR using the following conditions: initial denaturation for 30 seconds at 98°C, followed by 30 cycles of 94°C for 15 seconds, 50°C for 20 seconds, and 68°C for 2 minutes. A final elongation step of 68°C for 10 minutes was added, prior to a final hold at 4°C. Amplified cDNA was purified using a 1:1 ratio of PCR Clean DX beads (Aline Biosciences) and eluted in 50 uL nuclease-free water. Amplified cDNA was quantified using a Qubit 4 Fluorometer (Thermo Fisher Scientific) and fragment lengths were assessed using the TapeStation 2200 automated electrophoresis platform (Agilent).

Library Preparation and MinION Sequencing

Library preparation was performed using ONT’s ligation sequencing kit (SQK-LSK109 or SQK-LSK110). Multiplexing was performed using the native barcoding expansion 96 kit (EXP-NBD196). Library preparation was performed according to the manufacturer’s instructions, with several key modifications. DNA repair and end-prep was performed with 1000 fmol of input cDNA and the incubation times were increased to 30 minutes at 20°C, followed by 30 minutes at 65°C. For the barcoding reaction 200 fmol of input cDNA was incubated with the native barcodes and Blunt/TA Ligase Master Mix (New England Biolabs) for 20 minutes at room temperature (15-25°C), followed by 10 minutes at 65°C to improve barcode ligation efficiency with smaller fragments. Up to four clinical samples (90 fmol/sample) were multiplexed on each minION flowcell, with the addition of a blank viral transport medium negative control sample to each pooled library. Samples were sequenced on FLO-MIN106 flowcells on MinION MK1b sequencing devices for 72 hours using MinKNOW (Version 4.2.8, Oxford Nanopore Technologies) with live basecalling disabled.

Sequence Data Analysis

Raw fast5 files were basecalled using Guppy (Version 5.0.7, Oxford Nanopore Technologies) using the -- device cuda:0 flag to enable GPU basecalling. Output fastq files were uploaded to BugSeq (version 1.1, database version: RefSeq on Jan 28, 2021) for metagenomic classification (22), and results classification results were visualized in Recentrifuge (23). A representative html file containing an example visualization output can be found in the supplemental material (Supplementary Data). In brief, reads were demultiplexed with qcat using default run parameters (enforcing barcodes on both ends, which we have defined as stringent demultiplexing), followed by quality control with prinseq-lite. Reads shorter than 100bp or those deemed low quality (DUST score less than 7) were discarded. Reads were then classified against all of the microbial genomes in RefSeq, as well as the human genome and a library of common contaminants (see 22 for details). Reads classified as SARS-CoV-2 were extracted and used to build a consensus sequence with Medaka. Bases with less than 20X coverage were masked in accordance with public SARS-CoV-2 sequencing guidelines (https://www.aphl.org/programs/preparedness/Crisis-Management/Documents/APHL-SARS-CoV-2-Sequencing.pdf). SARS-CoV-2 lineages were assessed using Pangolin (Version 3.1.5, github.com/cov-lineages/pangolin), and phylogenetic analysis was performed with UShER (24) (Database: GISAID, GenBank, COG-UK and CNCB [2021-07-11]). Phylogenetic trees were constructed using augur (25), rooted at the SARS-CoV-2 reference sequence, and visualized in iTOL (26). Antimicrobial resistance genes were detected by aligning reads against the Resfinder database (27) with minimap2, disabling secondary alignments. Analysis from BugSeq outputs and visualizations were performed in RStudio (R version 4.1.0) and Python, with all code available at https://gitlab.com/bugseq/sars-cov-2-nanopore-mngs-performance (28).

Ethics Approval

This study obtained research ethics board approval from the University of British Columbia (H20-02152). Approval for collection of participant data was obtained by the Central University Research Ethics Committee at the University of Oxford (R69035). Specimens collected as part of routine testing at VGH and the BCCDC were de-identified and only contained a sample ID number, collection date, C_t, and VOC screening result.

Data Availability

Raw FASTQ data has been uploaded to NCBI Bioproject Accession PRJNA752146. Raw reads were mapped against the human reference genome and transcriptome (Ensembl hg38) using minimap2 (2.20) (29) and any human reads were removed.

Sequence data & Sample Descriptions

Amplified cDNA from a total of 43 patient swabs were sequenced on MinION sequencing devices. Of these samples, 38 were either positive or had indeterminate results based on SARS-CoV-2 RT-PCR and 5 samples had negative RT-PCR results. The 38 positive and indeterminate samples had a mean C_t value of 27.6 and ranged from 14.7-38.7 (Supplementary Figure 1). Sample collection dates, sample type, total read counts, as well as dual barcode reads, percent human reads, and SARS-CoV-2 reads per million reads sequenced (RPM) are present in Table 1. On average, negative controls exhibited a 29.7-fold decrease in dual barcode reads compared to the average number of dual barcode specimen reads (Mean dual barcode reads = 20,013, Q1: 158.5, Q3: 17556). SARS-CoV-2 was detected in similar abundance across our six positive control samples (Mean RPM Dual Barcode: 103,521 ± 21,070)

Sensitivity, Specificity, & Limit of Detection

We evaluated the test performance of our mNGS assay for detecting SARS-CoV-2. A sample was considered positive if one or more reads were assigned to SARS-CoV-2. Across all clinical samples, we detect SARS-CoV-2 with 78.4% (95%CI 62.8%-88.6%) sensitivity and 100% specificity (95%CI 56.6%-100%) (Table 2). Previous literature has demonstrated decreased sensitivity of mNGS assays above C_t 30 for other viruses (11,30). To assess the dependence of the mNGS assay on C_t value, we performed a subgroup analysis on samples above and below SARS-CoV-2 C_t 30. For samples with SARS-CoV-2 C_t < 30, sensitivity was 100% (95%CI 84.5%-100%), while for samples with SARS-CoV-2 C_t greater than 30, sensitivity was 50% (95%CI 27.8%-72.0%).

We note that two of 11 negative control samples had a single read assigned to SARS-CoV-2. We investigated these reads (further denoted as read one and two) to identify reasons for false positivity. Both reads had the expected barcode on both ends of the read as identified by BLAST. The first read exhibited 100% identity over the 24 nucleotide barcode on both ends, and the second read had 100% and 83% identity over the 24 nucleotide barcode on both ends. We next search these reads against the NCBI nucleotide database using megaBLAST to assess whether a BugSeq classification error occurred. However, both reads had top hits that exclusively matched SARS-CoV-2 with greater than 95% identity over more than 90% of their total length (923 and 1752 bases, respectively). SARS-CoV-2 was detected, despite strict dual barcode demultiplexing and removal reads with improper barcode insertions. Previous studies have identified barcode crosstalk, ranging from 0.2% to 0.3% of total classified reads, on Nanopore MinION flowcell results (31,32). When we examined the total SARS-CoV-2 read counts for a given flowcell on flowcells with false positive negative controls, we saw that one of those flowcells has the highest total SARS-CoV-2 read count of all flowcells in this study, therefore, we would expect higher levels of barcode crosstalk for that flowcell (Supplementary Figure 2).

We adjusted for barcode crosstalk by controlling for the total number of dual-barcoded SARS-CoV-2 reads on each flowcell. If we assume 0.2% of reads have incorrect barcodes ligated on both ends, and that these misclassified reads are evenly distributed across all barcodes on the flowcell, we can subtract the estimated number of misclassified reads from each sample. This correction yielded an acceptable threshold for classifying specimens as positive or negative. After adjusting for barcode crosstalk in this manner, we find that seven samples and two negative controls with SARS-CoV-2 reads detected would be re-classified as negative, and all negative controls are therefore classified correctly. The overall sensitivity and specificity on clinical samples after adjusting for barcode crosstalk are estimated to be 59.5% (95%CI 43.5%-73.7%) and 100% (56.6%-100%), respectively. Grouping by C_t value, the sensitivity estimates are 95.2% (95%CI 77.3%-99.2%) and 12.5% (95%CI 3.5%-36.0%%) for samples below and above C_t 30, respectively (Table 3).

RT-qPCR/SISPA Correlation, Genome Coverage, & SARS-CoV-2 Phylogeny

We assessed the relationship between SARS-CoV-2 RT-PCR C_t value and SARS-CoV-2 RPM for dual barcode reads, using stringent demultiplexing analysis parameters. SARS-CoV-2 log-RPM showed a strong linear association with RT-qPCR C_t value (R² = 0.71), with lower C_t values having a higher RPM on average (Figure 1). This relationship did not differ by RT-qPCR gene target (E-gene, ORF1ab, or RdRp) (Supplementary Figure 3). SARS-CoV-2 genome coverage depth showed a similar relationship, with decreasing coverage depth across the entire genome being associated with increasing C_t value (Figure 2, Table 4). We produced logistic regression models to assess the probability of attaining greater than 95% genome coverage at 1X, 20X, or 50X depth of coverage. We found that for every one unit increase in C_t value, the odds of recovering a 95% complete genome were 0.765 (95% CI: 0.519, 0.961), 0.263 (95% CI: 0.023, 0.666), or 0.263 (95% CI 0.023, 0.666) on average for coverage depths of 1X, 20X, or 50X, respectively (Figure 3). Interestingly, we did not see any difference in the likelihood of obtaining 95% coverage for 20X or 50X, despite slight differences in coverage depth for both of these thresholds (Table 4; Figure 3).

• Download figure
• Open in new tab

Figure 1:

Log SARS-CoV-2 reads per million reads sequenced across C_t value (E gene or ORF1ab) for all RT-qPCR positive samples. 95% confidence intervals for the linear regression line are shaded in grey. Coefficient of determination = 0.71.

• Download figure
• Open in new tab

Figure 2:

Coverage depth for samples classified as positive by our classifier with log depth of coverage on the y-axis and SARS-CoV-2 reference genome position on the x-axis

• Download figure
• Open in new tab

Figure 3:

Probability of obtaining greater than 95% genome coverage (1 = Yes, 0 = No) for RT-qPCR positive study samples across C_t value for a. 1x, b. 20x, and c. 50x genome coverage. Logistic regression models are represented in blue.

SARS-CoV-2 metagenomic reads were used to reconstruct viral genomes. We produced ten complete (greater than 95% unambiguous bases) and one near-complete consensus genome sequence (greater than 80% unambiguous bases) from our 20 SISPA-positive clinical specimens, masking any bases with less than 20X coverage. Two partial viral genomes were constructed with 20-25% unambiguous bases. Pangolin lineage assignment was successful to all complete or near complete genomes; of these five underwent SARS-CoV-2 VOC PCR testing. All five whole or partial viral genomes were classified as SARS-CoV-2 lineages concordant with PCR results (Table 5). We also detected an additional VOC in a sample that did not undergo VOC PCR testing. We also assessed our complete or near-complete genomes in the context of global SARS-CoV-2 transmission by placing them in a phylogenetic tree containing over two million publicly available SARS-CoV-2 genomes. The ten complete genomes could be placed in the global phylogeny with high confidence (only one maximally parsimonious placement), and the near-complete genome could be placed with lower confidence (nine maximally parsimonious placements). For ten of 11 genomes derived from metagenomic data, the nearest neighbor in this tree was a genome derived from the same province of sample collection, British Columbia. Additionally, for 9/11 study genomes, 80% or more of the nearest 50 genomes were derived from British Columbia; for the remaining two study genomes, 90% or more of the nearest 50 genomes were derived from Canada (Figure 4; Supplementary Figure 4). The UK samples did not yield well covered genomes. Subtrees with nearest neighbors for all study samples are available in the supplementary material (Supplementary Figure 4).

• Download figure
• Open in new tab

Figure 4:

Study samples (marked as P1, P2, etc.) and their nearest three neighbors from all publicly available SARS-CoV-2 sequences

Universal Microbial Detection & Antimicrobial Resistance

We searched the BugSeq metagenomic output of our clinical specimens for alternative respiratory viruses or viral or bacterial co-infections.We did not identify any other pathogenic viruses or atypical bacteria such as Chlamydia pneumoniae or Mycoplasma pneumoniae. We did identify several members of the normal nasopharyngeal microbiota, which when found in the lower respiratory tract, may cause disease; these included two samples with Moraxella catarrhalis, seven samples with Haemophilus influenzae or Haemophilus parainfluenzae, three samples with Neisseria meningitidis, three samples with Staphylococcus aureus, two samples with Streptococcus pneumoniae and two samples with Klebsiella pneumoniae (Supplementary Table 3). These results are consistent with other metagenomic sequencing results from the nasopharynx (30). We searched our data for genes conferring antimicrobial resistance, and identified 10 genes across 6 samples. We found two beta-lactamases in our dataset: blaTEM-234, a class A beta-lactamase which has undetermined spectrum and derived from Escherichia coli in sample P22, as well as blaOXA-85, which confers resistance to amoxicillin and amoxicillin-clavulanate, that derived from Fusobacterium psuedoperiodonticum (P9).