Evaluating metagenomics and targeted approaches for diagnosis and surveillance of viruses

Sensitivity and limit of detection

We tested simulated post-extraction clinical samples where the input viral composition is known, consisting of a mock community of genomic DNA/RNA from six viruses, two DNA (human mastadenovirus F and human betaherpesvirus 5) and four RNA (mammalian orthoreovirus, human orthopneumovirus, influenza B virus and Zika Virus), at four different concentrations in a constant human DNA and RNA background (Figure 1). The same input was used for all metagenomics approaches: the untargeted Illumina and ONT protocols, and the capture probe enrichment with the Twist Comprehensive Viral Research Panel (VRP) followed by Illumina sequencing (Figure 1). At least two replicates were tested for each technology-concentration pair. We obtained 38.2-81.2 and 38.8-66.9 million reads per sample for untargeted Illumina sequencing and Illumina following the Twist VRP respectively, corresponding to 5.7-12.2 Gb and 5.8-10.0 Gb respectively (Supplementary Table 1). For the untargeted ONT sequencing, we obtained 5.1-12.3 Gb per sample. To improve comparability between methods, we randomly subsampled 5 Gb from each sample across the platforms for analysis.

• Download figure
• Open in new tab

Figure 1: Metagenomic sequencing and experimental outline

A Overview of metagenomic processing pipeline. B Flow chart summarizing experimental design, which involves inputting mock and clinical samples into three metagenomic workflows: Illumina DNA and RNA seq using NEBNext and KAPA kits respectively, ONT DNA and RNA seq using the Rapid PCR barcoding kit and the Rapid Smart-9N method respectively, and finally the targeted DNA- and RNA-based Twist viral research panel, sequenced on the Illumina platform. The resulting data was analyzed using different taxonomic classifiers. Produced with biorender.com.

The Twist VRP was the most sensitive method, as it was the only platform to detect all the expected viruses at 60 genome copies per ml (gc/ml), with coverage over 98.8% for all viruses at 60,000 gc/ml and ranging from 3.7-23.0% at 60 gc/ml (Figure 2A). ONT was less sensitive than Illumina, detecting in at least one of the replicates all six viruses at 60,000 gc/ml, four of six viruses (human betaherpesvirus 5, human mastadenovirus F, orthopneumovirus and Zika virus) at 6000 gc/ml but only two viruses, one double-stranded (ds) DNA (human betaherpesvirus 5) and the other dsRNA virus orthoreovirus at 600 gc/ml and none at 60 gc/ml. The detection of the dsRNA virus orthoreovirus at 600 gc/ml despite not being detected at 6000 gc/ml, represents only four reads in one of the replicates, with no reads detected in the other replicate, likely reflecting stochastic variation. In contrast, untargeted Illumina detected all six viruses at 60,000 and 6,000 gc/ml, five at 600 gc/ml (all apart from human mastadenovirus F) and one at 60 gc/ml (human betaherpesvirus 5) (Error! Reference source not found. 2A). At levels close to the limits of detection, there was sometimes variation between the technical replicates in their ability to detect the viruses (Figure 2). One additional DNA virus (human mastadenovirus F) and one additional RNA virus in one of the repeats (human orthopneumovirus) were detected by Illumina sequencing at 600 gc/ml when additional sequence data was available beyond 5 Gb (9.6 and 10.7 Gb for DNA and 11.1 Gb for RNA) (Supplementary Table 1). Other than this, no additional viruses were detected in the full datasets before subsampling.

• Download figure
• Open in new tab

Figure 2: Sensitivity in mock clinical samples

Coverage and base pairs aligned to the six expected viral species in mock samples, by untargeted Illumina and ONT sequencing and capture probe enrichment with the Twist Biosciences Comprehensive Viral Research Panel followed by Illumina sequencing. A Percentage genome coverage at depth 1x of species in mock community. B log10(bases) aligning to reference genome. Samples where a virus was detected in the full dataset but not the subsampled dataset are indicated with a *. Genome copy numbers refer to an average across the viral species – see Supplementary Table 3. Each point shows the mean of at least two technical replicates – error bars show the range. PCR duplicate reads removed.

At 60,000 gc/ml, assigned bases ranged from 17,527-217,630 of 5 Gb for Illumina and 2110-134,026 of 5 Gb for ONT (Figure 2B). Both ONT and Illumina untargeted sequencing provided incomplete coverage at all concentrations tested, ranging from 8.4 to 83.9% at 60,000 gc/ml, 0-30.7% at 6000 gc/ml, and 0-8.9% at 600 gc/ml (Figure 2A). Viruses with longer genomes were detected in with greater read numbers, however normalizing for genome length gave similar abundance estimates for each virus, where viral loads were high enough for consistent detection (Supplementary Figure 1). As expected, all technologies displayed levels of PCR duplication ranging from 0-16%, with the Twist VRP showing the highest rates (Supplementary Figure 2A). Including PCR duplicates makes no difference to the conclusions regarding sensitivity (Supplementary Figure 2B).

We also tested the sensitivity of a range of taxonomic classifiers. The classifiers tested and reasons for inclusion are outlined in Table 1. Where no thresholds were applied, all the classifiers had similar sensitivity, although there was some variation in ability to detect viruses at 60-6000 gc/ml for untargeted Illumina sequencing and at 60,000 gc/ml for untargeted ONT sequencing, with Kraken2, Dragen, metaMix-fast and CZ ID being the most sensitive at these viral loads (Figure 3A). MetaMix and MEGAN-LR failed to identify influenza B virus and mammalian orthoreovirus respectively with ONT sequencing at 60,000 gc/ml; both RNA viruses for which fewer than 10 reads were detected by the aligner minimap2. Of the other classifiers, One Codex had substantially lower sensitivity for the Twist VRP data compared to other classifiers, all of which identified almost all the viruses at all concentrations tested (Figure 3A). This may be because the program only reports organisms that reach a set of predetermined abundance thresholds (56), which may not be reached at low viral loads, whilst the other classifiers do not by default use such thresholds. Where viruses were detected, the classifiers provided broadly similar estimates of reads per million, ranging, for example, from 30.3-73.6, 31.0-39.2 and 3695-8164 RPM, for human betaherpesvirus 5 for Illumina, ONT and the Twist VRP respectively (nucleotide-based classifiers only) (Supplementary Figure 3).

• Download figure
• Open in new tab

Figure 3: Sensitivity and specificity of taxonomic classifiers

A Sensitivity to the species in the mock community before and after the application of thresholds, for seven different taxonomic classifiers, by untargeted Illumina and ONT sequencing and capture probe enrichment with the Twist Biosciences Comprehensive Viral Research Panel followed by Illumina sequencing. MEGAN-LR and the One Codex Twist report are only designed for ONT and Twist sequencing respectively so were only run for these platforms. B,C Number of false positive species, defined as a species that is classified as positive but not present in the mock community B False positive species from the raw output of the taxonomic classifiers with no thresholds applied. C Comparison of the numbers of viral positive species identified before and after the application of thresholds. Genome copy numbers refer to an average across the viral species – see Supplementary Table 3. Each bar shows the mean of at least two technical replicates.

Specificity and false positive rates

High precision and low false positive rates are as important as sensitivity in a clinical diagnostic setting, and rational approaches to identifying and reporting contaminants, particularly by non-specialist bioinformaticians, are needed. To evaluate their performance, we compared the number of false positive species identified by a range of commonly used taxonomic classifiers for the mock samples (Table 1). A false positive is defined as any species not present in the mock community. All the classifiers assigned similar numbers of reads to the species in the mock community, except for One Codex, which had lower sensitivity for the Twist VRP data than the other classifiers. (Supplementary Figure 3). However, when no additional thresholds were applied, there was a large variability between the classifiers in terms of the number of species identified by Illumina sequencing (Figure 1B). Most of the false positive species were fungi or bacteria. Kraken2 Illumina’s Dragen Metagenomics Pipeline (which is based on Kraken2) and One Codex, all use kmer methodologies and identified over 1500 false positive species for the untargeted Illumina sequencing. The discrepancy between the number of false positives identified for the Twist VRP data by One Codex at different concentrations may be caused by greater availability of data for the classifier to distinguish between true and false positives at higher read depths (56). By contrast, metaMix and Bracken, which both use Bayesian methods, identified only one false positive viral species at 60,000 gc/ml, (Figure 3A&C). However, both these classifiers were less sensitive at lower genome copy numbers than classifiers such as Kraken2 and CZ ID. In contrast to Illumina, few false positive species, especially viruses, were identified with ONT sequencing. Thus, for ONT the application of thresholds beyond a basic comparison to the negative control may not be required.

To reduce the number of false positive species identified for Illumina sequencing, we imposed more stringent thresholds. Completely disregarding all species with any reads in the negative control may result in a reduction in sensitivity, particularly when there is low-level cross-contamination from high viral load samples into the control. We found that using a combination of reads per million ratio between sample and the corresponding negative control and proportion of microbial reads resulted in optimum sensitivity (91.7%) and specificity (77.4%), which may be useful for classifiers such as Kraken2 and One Codex which require additional thresholds (Figure 1A&C). In contrast, ONT sequencing and classifiers such as metaMix have few false positive reads and can be used with only a comparison to the negative control. More details of the derivation of our thresholds can be found in the supplementary information. Use of protein-based classifiers, including Kaiju (57) and the protein modes of MEGAN-LR, metaMix and CZ ID, did not improve the sensitivity or specificity of classification (Supplementary Figure 4). The false positive viral species that remained after the application of thresholds were mainly viruses that do not infect mammals or birds, making them unlikely to be clinically relevant (Supplementary Figure 5). The remaining false positive viruses were mainly Anelloviridae (often Torque Teno viruses), and viruses that were related to those in the mock community, such as other herpes or adenoviruses. The Anelloviridae, which are very commonly found in human samples, were mostly not found in the controls and are possibly a result of low-level contamination.

Host transcriptomic analysis

Several studies highlight the power of host transcriptomics methods for distinguishing bacterial, viral, and non-infectious causes of illness (27,50–53), although none are being used diagnostically at present. When metagenomics does not identify any pathogens, such analysis could help distinguish between a non-infectious cause of disease and a lack of sensitivity of the metagenomics protocol. Since Illumina RNA sequencing has been extensively used and validated for transcriptomic studies, we compared the estimates of human gene expression provided by the ONT and Twist VRP platforms to those from Illumina. Although the Twist VRP only enriches for viruses, it retains the background, meaning that this analysis remains possible. The number of reads assigned to each human protein-coding gene were positively correlated between Illumina and the other two technologies (correlation coefficients, Spearman’s rho, 0.694 and 0.709 for ONT and the Twist VRP respectively) (Figure 4A-C). Due to the combined DNA and RNA protocol used with the Twist VRP, there were a large number of human genes that were identified as highly expressed by the panel but not untargeted Illumina (Figure 4B). We therefore repeated the analysis, focusing only on reads that mapped across exon-exon junctions, termed henceforth “spliced reads”, which are likely to represent mRNA, resulting in a better agreement between the Illumina and the Twist VRP results (Figure 4D).

• Download figure
• Open in new tab

Figure 4: Host transcriptomic analysis

A-D Read counts per million assigned to each gene in the human genome by untargeted Illumina, untargeted ONT and targeted Illumina sequencing using the Twist Viral Research Panel. Each point represents a gene. A-C raw reads; D only reads that map across splice junctions. E number of genes identified by each pair of technologies. F counts per million of reads by platform. Each panel shows the log2(CPM) as estimated by a different technology. Outliers not shown. All comparisons are statistically significant (p < 0.01) with a pairwise Wilcox test other than those indicated.

While most protein-coding genes were identified by all the technologies (Figure 4E), there was still a substantial minority that were not identified by ONT (Figure 4E). Use of spliced reads for untargeted Illumina and Twist VRP, only resulted in a small drop in the number of protein-coding genes identified, and a larger drop in the non-coding transcripts (Figure 4E). Genes that were identified by all technologies were significantly more highly expressed (Figure 4F) suggesting that low-expressed genes may be less reliably identified by all technologies, particularly ONT.

Turnaround time and cost

Costs and turnaround times from sample to results affect the adoption of metagenomics for routine diagnostics. ONT provides the quickest library preparation method, at just over 5 hours for both DNA and RNA protocols (Figure 5A). Targeted sequencing with the Twist VRP requires overnight hybridization and is the slowest protocol (Figure 5A). The Twist VRP protocol was the cheapest based on 23 samples (+ negative control) and a sequencing depth of 5 Gb, while the untargeted Illumina sequencing was the most expensive (Figure 5B, Supplementary Table 2).

• Download figure
• Open in new tab

Figure 5: Turnaround times and output data volumes

A Time taken for library preparation for the different protocols tested. The Twist panel uses a combined DNA and RNA-Seq protocol. The DNA+RNA bars for the untargeted sequencing indicate the time taken if both protocols are performed by a single operator. B Total cost (including library preparation) to sequence number of samples indicated plus single negative control, to a depth of 5GB. ONT costs are shown with 48- and 72-hour maximum run times per flow cell. C Volume of data output by time for a range of Illumina sequencing kits and ONT sequencing with PromethION flow cells. The Illumina kits produce a set amount of data after the sequencing run is complete – this is shown by pink dots. In ONT sequencing, data is output continuously and the run can be stopped at any time, until the flow cell becomes degraded. PromethION data (green/blue dotted lines) shows the average of our RNA and DNA-Seq runs, passed reads only. Data outputs for Illumina were obtained from the product specification data as of April 2024.

However, directly comparing the costs of each protocol to obtain at least 5 Gb of sequence data does not account for differences in their sensitivity. Since the Twist VRP approach is at least 10-100x more sensitive than untargeted Illumina and 100-1000x more than untargeted ONT (Figure 2), increases of orders of magnitude in sequencing depth would be required to bring the sensitivity of the untargeted protocols in line with that of the Twist VRP. Even up to twice as much sequence data (10.6 Gb and 11.3 Gb at 60 gc/ml) did not increase the sensitivity of untargeted Illumina and ONT respectively to anything near to the Twist VRP. For untargeted Illumina, greater sequencing depth also amplifies the detection of contaminants, making interpretation more difficult. Achieving increased sensitivity using ONT sequencing would require long sequencing runs and a reduction in the number of samples sequenced per flow cell, significantly increasing costs and turnaround times. This means that targeted metagenomics methods such as Twist VRP are by far the quickest and most cost-effective of the protocols for detection of low viral loads (60-600 gc/ml). Similarly, since untargeted Illumina is more sensitive than this untargeted ONT protocol, it will be quicker and cheaper to reach the sequencing depths required to detect intermediate viral loads (600-6000 gc/ml) using Illumina.

Sequencing costs and turnaround times will also be influenced by the number of samples. For fewer than six samples, including controls, ONT is the cheapest and fastest alternative where microbial load is likely to be high and genomic sequences are achievable with lower sequencing depth, for example 5 Gb of data, per sample (Figure 5A&B). ONT also provides access to the sequencing data in real time, allowing preliminary analysis of the results before the run is completed, which can be advantageous for samples with high viral loads. However, if more samples are processed in parallel or a higher sequencing depth is required to improve sensitivity to a level comparable with untargeted Illumina, longer sequencing runs will be needed (Figure 5B&C). When the total volume of data required per run is higher than around 30 Gb, it may be faster to use Illumina sequencing (Figure 5C). However, it remains cheaper to use ONT with 23 sample runs (24 including negative control, 120 Gb) (Figure 5B). Because of the Twist VRP’s improved sensitivity, lower sequencing depths are required per sample, allowing the use of smaller Illumina sequencers and cheaper kits with shorter sequencing times (Figure 5C). However, for fewer samples the Twist VRP method may be much more expensive, since the kit optimal cost per sample is based on the pooling of 7 (8 including negative control) samples per hybridization (Figure 5B, Supplementary Table 2).

Mock clinical samples

Mock samples were prepared to represent high-biomass samples (e.g. blood and tissue) with a clinically relevant spectrum of viral loads ranging from 60 to 60,000 gc/ml. This was achieved by performing serial dilutions of a commercial genetic material mix – the ATCC Virome Nucleic Acid Mix (ATCC, MSA-1008) (Table 2) in a background of either human DNA, RNA or a DNA+RNA mix at a concentration of 40 ng/µl. Mock samples were produced in a single large batch to reduce variation and stored as 10 µl single-use aliquots. In parallel, a batch of negative controls containing either DNA, RNA or DNA+RNA without any viral nucleic acids were also generated. DNA and samples were spiked with Lambda phage DNA and RNA samples were spiked with MS2 Bacteriophage RNA. DNA + RNA mixed samples were spiked with both standards. Details of mock sample composition can be found in Supplementary Table 3.

Untargeted Illumina sequencing

Untargeted Illumina DNA and RNA metagenomic sequencing of the mock clinical samples was performed as previously described (90). Two technical replicates per mock sample were performed. DNA samples underwent human CpG-methylated DNA depletion using the NEBNext^® Microbiome DNA Enrichment Kit (New England Biolabs; E2612L) followed by library preparation using the NEBNext® Ultra™ II FS DNA Library Prep Kit for Illumina (New England Biolabs; E7805L). RNA samples underwent ribosomal RNA (rRNA) depletion followed by library preparation using KAPA RNA HyperPrep kit with RiboErase HMR (Roche; KK8561).

All pre-PCR steps were carried out under an MSC class II cabinet and moved to a post-PCR area following amplification. Libraries were quantified with high sensitivity dSDNA kit (Invitrogen; Q33231) on an Invitrogen Qubit 4 Fluorometer and the average peak sizes for libraries were checked using high sensitivity D1000 screentapes (Agilent 5067-5584) on a Tapestation 4200. Samples were sequenced in equimolar pools using a NextSeq 2000 or a NovaSeq 6000 300 cycle kit (2 x 150 bp) depending on the number of samples processed. A minimum output of 5 Gb per sample was obtained (Supplementary Table 1).

ONT sequencing

ONT sequencing was performed using PCR-based protocols and Q20+ chemistry (Version 14 kits). Two technical replicates per mock sample were performed. DNA samples underwent human CpG-methylated DNA depletion using the NEBNext Microbiome DNA enrichment kit (New England Biolabs) prior library preparation using the Rapid PCR Barcoding kit 24 V 14 (SQK-RPB114.24) according to manufacturer’s instructions. RNA sequencing was performed using Rapid-Smart 9N(91). First, for annealing of the tagged random oligonucleotide, 10 µl of RNA was mixed with 1 μl of 2 μM RLB RT 9N oligo (TTTTTCGTGCGCCGCTTCAACNNNNNNNNN) and 1 μl 10 mM dNTPs. Mix was incubated for 5 min at 65°C, then cooled on ice. For cDNA synthesis and generation of double-tagged cDNA, 4 μlSuperScript IV First-strand Buffer, 1 μL 0.1 M DTT, 1 μl RNase OUT (Thermo Fisher Scientific), 1 μl 2 μM RLB TSO (GCTAATCATTGCTTTTTCGTGCGCCGCTTCAACATrGrGrG), and 1 μL SuperScript IV (Thermo Fisher Scientific) was mixed with the 12 μl annealed RNA. Reaction was incubated for 90 min at 42°C followed by 10 min at 70°C. 5 µl of double-tagged cDNA were used as input for the PCR step in the Rapid PCR Barcoding kit 24 V 14 (SQK-RPB114.24). From this step onwards, manufacturer’s instructions were followed.

All pre-PCR steps were carried out under an MSC class II cabinet and moved to a post-PCR area following amplification. Sequencing was performed using PromethION Flow cells (R.10.4.1) on a P2 solo device connected to a GridION. Real-time basecalling was performed in MinKnow Version 23.07.5 using the high-accuracy model. Samples were sequenced until a minimum output of 5 Gb per sample was obtained.

Targeted Illumina sequencing with Twist Comprehensive Viral Research Panel

Targeted Illumina sequencing was performed on samples with a combined DNA + RNA background using the Twist Comprehensive Viral Research Panel (Twist Biosciences, 103550) following the Twist Biosciences Total Nucleic Acids Library Preparation EF Kit 2.0 for Viral Pathogen Detection and Characterization protocol. Two technical replicates per mock sample were processed other than for 60 and 600 gc/ml and the negative control, where four replicates were performed. Additional replicates were included to thoroughly test for potential cross-contamination and to assess potential sensitivity loss in low copy number samples when combined with high copy number samples in hybridisation-capture reactions.

First, cDNA synthesis was performed using ProtoScript II First strand synthesis kit (New England Biolabs, E6560) followed by the NEBNext Ultra Non-Directional Second Strand Synthesis module (New England Biolabs, E6111). 25ng of the double-stranded cDNA and dsDNA mix were used as input for adapter ligation, indexing and pre-capture amplification using the Twist Library preparation EF Kit 2.0 (Twist Biosciences, 104207 + 100573). All pre-PCR steps were carried out under an MSC class II cabinet until the indexing step was complete.

Following pre-capture amplification, indexed samples were pooled, for a total of 7 samples plus a negative control per hybridisation reaction. Hybridisation was performed overnight for 16 hours. Hybridisation targets were then captured with Streptavidin Binding Beads. At this step, samples were washed using the Twist Wash Buffers (Twist Biosciences, 104178) instead of the washing buffers V2 as per recommendation of the manufacturer. Post-capture amplification was performed on the enriched libraries (8 cycles). Final enriched libraries were quantified with Qubit high sensitivity kit and average peaks obtained with high sensitivity D1000 tapes. Samples were sequenced in equimolar pools using a NextSeq 2000 or a NovaSeq 6000 300 cycle kit (2 x 150 bp) depending on the number of samples processed. A minimum output of 5 Gb per sample was obtained (Supplementary Table 1).

Databases for taxonomic classification

Since database composition has been shown to have a significant impact on the results of metagenomics (92), a common set of sequences was used to build the databases where possible. For the tools where it was possible to create a custom database (Kraken2 (58), Bracken (59), Dragen Metagenomics Pipeline (60), EPI2ME labs wf-metagenomics (61), metaMix (63), MEGAN-LR (62) and Kaiju (57)), a database was created based on the bacterial (complete genomes only), viral, fungal, protozoa and human nucleotide from Refseq (downloaded 6^th June 2023). Databases were built using the default parameters, other than for MEGAN-LR, where the recommended settings for ONT data described in (39) were used. A common set of taxonomy files downloaded from NCBI (31^st July 2023) were also used. Unplaced contigs were removed from the parasites and fungal nucleotide sequences prior to building the databases to reduce human contamination present in some of the reference sequences. It is not currently possible for the user to alter the databases for CZ ID (64) or One Codex (65), so the inbuilt databases were used.

Read preprocessing and taxonomic classification

Reads were randomly subsampled from the raw output fastq files, using seqtk (93) sample for the Illumina data and a custom python script for the ONT data, to obtain 5 Gb for each sample across all the technologies.

Kraken2, Bracken and Kaiju were run through the nf-core Taxprofiler pipeline (94), which aims to provide a reproducible best-practice workflow for metagenomics analysis. As recommended, read preprocessing involving adaptor trimming and complexity filtering with fastp (95) was performed for Illumina but not ONT sequencing (96). Host removal was performed for both platforms by alignment to the human genome (version Ch38).

The reads obtained following preprocessing and host removal from the Taxprofiler pipeline were used as input to MEGAN-LR, run through the PB-metagenomics tools pipeline(97), with the adjustments for ONT sequencing recommended in (39).

For Illumina data processed with metaMix, a separate preprocessing pipeline was used for a more thorough removal of host reads. This involves read trimming using TrimGalore (98), followed by removal of human DNA/RNA and ribosomal RNA using alignment with both Bowtie2 (99) and BLAST (100,101). For the other classifiers, the time saved in classification was shorter than the time taken for the longer host removal pipeline, so a single alignment step is sufficient. For the ONT data, the output of the preprocessing and host removal from Taxprofiler was used as input. Reads were then aligned to the reference database with BLAST (nucleotide mode) and DIAMOND (102) (protein mode) before input to the metaMix R package. metaMix-fast is the first two steps of the metaMix R package, before the time-consuming MCMC step.

Raw reads were uploaded to CZ ID metagenomics workflow through the online interface. For CZ ID, the fields nr_count or nt_count were used for protein and nucleotide analyses respectively. Raw reads were also uploaded to the One Codex platform. The reads field from the reads field used for the analysis. For the Twist panel data, the Twist Comprehensive Viral Research Panel report was also run. Initially, this report failed to identify Reovirus, but this issue has since been rectified by One Codex.

Raw reads were also used as input into either Illumina’s Dragen Metagenomics Pipeline or ONT’s EPI2ME labs wf-metagenomics as appropriate. The same Kraken2 custom database as used for running Kraken2 through nf-core Taxprofiler was used, and both tools were run using the command-line interface.

Alignment and sensitivity analysis

Reads were aligned to reference genomes downloaded from ACC using Bowtie2 (99) for the Illumina data and Minimap2 (103) for the ONT data, using the “very-sensitive” mode and the default parameters respectively. PCR duplicate reads were removed before calculating coverage and depth using samtools (104).

Specificity analysis

To standardize the results between classifiers for comparison, taxonomic ranks were identified, organisms were classified as bacteria, viruses, fungi or other eukaryotes, and all reads assigned to taxonomic levels below species were assigned to the relevant species, using custom R scripts and the taxonomizr package (105). Where an organism was detected by both DNA and RNA sequencing, the result with the higher number of reads was retained. All analysis was performed in terms of reads rather than base pairs since not all classifiers output assignments by read, making it impossible to calculate base pair assignments for the ONT data. Read per million ratios and proportion of microbial reads were calculated and used to identify positive species as described in the supplementary information. False positive species were defined as species that were identified by the classifiers and were not present in the mock community or the positive controls. False positive viral species were classified according to host using the Virus-Host DB (106).

Host transcriptomic analysis

Genes and transcripts were quantified using Kallisto(107) with human genome GRCh38.p14 downloaded from Gencode (108). Analysis was conducted in R using the tximport (109) and rtracklayer (110) packages. Spliced reads were identified by alignment to the human genome using STAR (111), using the presence of the CIGAR string to identify gapped alignments.

Plots

Plots were produced in R using Tidyverse (112) packages or using Biorender.com.