Sequencing Illumina libraries at high accuracy on the ONT MinION using R2C2

Evaluating R2C2 for the sequencing of Illumina RNA-seq libraries

First, we benchmarked the ONT-based R2C2 method for the generation of RNA-seq data from Illumina libraries. We prepared four technical replicate libraries from a single RNA sample in the form of dual indexed paired-end Illumina libraries using the NEBNext Ultra II directional RNA kit with RNA of the human lung carcinoma cell line A549. We pooled and sequenced these libraries with the ONT MinION both directly (1D) and after R2C2 conversion (R2C2), as well as with the Illumina MiSeq.

To establish the effect of R2C2 conversion on the throughput of the ONT MinION when sequencing short Illumina libraries, we processed the raw reads generated by both 1D and R2C2 sequencing runs. Raw read numbers for 1D and R2C2 runs generated from one ONT MinION flowcell were similar at around 11.8 million reads. However, 1D reads were less likely than R2C2 reads to (1) pass filter during base-calling, (2) contain both p5 and p7 Illumina adapter sequences, and (3) be successfully demultiplexed. After preprocessing, only 2.5 million 1D reads (21%) remained compared with around 8 million R2C2 reads (Table 2).

The 1D read numbers we generated for the RNA-seq libraries are similar to published 1D read numbers generated for libraries of similar lengths. A recent large-scale study on GTEx samples (Glinos et al. 2022) sequenced ∼600-nt-long cDNA molecules across dozens of flowcells and generated about 6 million reads per MinION flowcell. The Long-read RNA-seq Genome Annotation Assessment Project (LRGASP) Consortium (Pardo-Palacios et al. 2021) sequenced ∼520-nt-long cDNA molecules and generated about 18 million reads per flowcell. Even the most productive 1D run in these studies, potentially generating up to 20 million raw reads for molecules of this length (Pardo-Palacios et al. 2021), would still generate fewer demultiplexed reads (21% of 20 million or fewer than 5 million) than the R2C2 run we performed here.

To validate the demultiplexing of Illumina library pools from R2C2 data, we compared the ratio of reads assigned to each library in Illumina MiSeq, R2C2, and ONT 1D data based on their combination of i5 and i7 indexes. For all three methods, three technical replicate libraries were pooled at a 4:2:1 ratio. The Illumina MiSeq produced a 4:2.03:1.58 read ratio after demultiplexing. R2C2 produced a 4:1.91:1.34 ratio, and ONT 1D produced a 4:2.5:1.82 ratio. With these results being quite similar, the differences are likely owing to pipetting variability when pooling the libraries for the different sequencing methods. Further, to evaluate our ability to quantitatively pool libraries at different points in the R2C2 workflow, we processed a fourth replicate in parallel and added it at a specific ratio after rolling circle amplification. The fourth replicate represented 40.5% of the R2C2 data, which is slightly more than the 30% of R2C2 DNA it represented in the MinION sequencing run. Finally, 9.71% of the R2C2 reads were not assigned to any index combination and 1.7% of the R2C2 reads were assigned to index combinations not present in the pool, implying only 0.0289% (1.7% × 1.7%) of the R2C2 reads were assigned to the wrong index combination owing to index hopping.

Next, we established the effect of R2C2 conversion on read accuracy compared with ONT 1D and Illumina MiSeq data sets. We aligned all complete p5- and p7-containing and demultiplexed R2C2 (8,066,704) and 1D reads (2,530,950) as well as Illumina MiSeq reads (20,830,560 2 × 300-bp paired-end reads) generated from these RNA-seq libraries using minimap2. We then calculated the median read accuracy, accuracy per base, and read position–dependent accuracy per base (Table 3).

Although median read accuracy is a useful and often reported metric to compare error-prone long-read sequencing technologies, it becomes less useful in this study. The sequencing reads we aim to compare are very short—either owing to the short length of the molecules sequenced (1D and R2C2) (Fig. 2A) or to technology limitations (Illumina MiSeq)—and often accurate enough to be unlikely to contain errors at that length, causing many individual sequencing reads to be 100% accurate. This is obvious with read 1 of the Illumina MiSeq having a median accuracy of 100%, which contains little information on the real Illumina MiSeq error rate. Accuracy per base (%), that is, (correct bases of all reads/all bases of all reads) × 100, is a more useful metric to compare accurate short reads. Using this metric, we see that 1D reads are the least accurate with an accuracy per base of 96.59%. R2C2 falls between Illumina MiSeq read 1 (99.47%) and read 2 (98.57%) with an accuracy per base of 98.87%. Further, although R2C2 reads contained more deletion and insertion errors, they contained fewer mismatch errors than both Illumina MiSeq read 1 and read 2.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Sequencing Illumina RNA-seq libraries on the ONT MinION after R2C2 conversion. Insert length distribution (A) and read position–dependent identity to the reference genome (B) of R2C2 and Illumina MiSeq reads generated from the same Illumina library. (C) Comparisons of R2C2 and Illumina MiSeq read-based gene expression and splice junction usage quantification by STAR and kallisto are shown as scatter plots with marginal distributions (log₂ normalized) shown as histograms. (D) Genome browser-style visualization of read alignments to the Actb locus. Mismatches are marked by lines colored by the read base (A, orange; T, green; C, blue; G, purple). Insertions are shown as gaps in the alignments, and deletions are shown as black lines.

Read position–dependent accuracy of 1D, R2C2, and Illumina MiSeq read 1 and read 2 adds further detail to this comparison. In contrast to 1D and R2C2 data, Illumina MiSeq base accuracy decreased with increasing read cycles, particularly in read 2, with R2C2 surpassing Illumina MiSeq accuracy for read 2 lengths over ∼175 bp (Fig. 2B,D). To ensure that our Illumina MiSeq run was not an outlier in terms of accuracy typical of Illumina benchtop sequencers, we performed the same position-dependent accuracy analysis on publicly available Illumina MiSeq, iSeq, and MiniSeq data, which showed the same overall trends (Supplemental Fig. S1).

Next, we aimed to establish whether R2C2 RNA-seq and ONT 1D data could be analyzed using computational tools designed and established for Illumina RNA-seq data. To quantify gene expression levels, we aligned and evaluated the entire demultiplexed R2C2 (8,066,704 reads) and ONT 1D (2,530,950 reads) data sets as well as our Illumina MiSeq data set (20,830,560 read pairs) using the STAR aligner (Dobin et al. 2013) (STARlong executable for R2C2 and ONT1D data), which is routinely used for standard Illumina RNA-seq analysis; 7,365,398 R2C2 reads (91.66%), 1,834,065 ONT 1D reads (72.48%), and 18,649,031 Illumina MiSeq reads (90.08%) mapped uniquely to the human genome. The STAR log files indicated that compared with R2C2 (1) the low overall accuracy of ONT 1D reads and (2) the declining quality of MiSeq reads means STAR aligner is less likely to align them. Although more forgiving aligners like minimap2 exist, they are not intended for spliced short-read alignments and therefore not optimized for this use-case.

Based on these read alignments, STAR determined normalized gene counts for the Illumina MiSeq, R2C2, and ONT 1D data sets. Illumina MiSeq gene counts showed Pearson's r-values of 0.995 and 0.987 compared with R2C2 (Fig. 2C) and ONT 1D, respectively. Additionally, STAR also determined normalized splice junction counts for the three data sets, which provide a higher-resolution view of the transcriptome. Illumina MiSeq splice junction counts showed Pearson's r-values of 0.974 and 0.929 compared with R2C2 (Fig. 2C) and ONT 1D. Finally, we also tested whether ultrafast pseudoalignment-based tools will generate reliable gene expression levels based on R2C2 and ONT 1D reads that feature more insertion and deletion rates compared with standard Illumina data. We used one such tool, kallisto (Bray et al. 2016), and found that gene expression values as determined for Illumina MiSeq had Pearson's r-values of 0.985 and 0.973 compared with R2C2 (Fig. 2C) and ONT 1D.

Overall, this comparison showed that using R2C2 we can convert Illumina RNA-seq libraries into DNA ideally suited for the ONT MinION. Not only does R2C2 generate more reads than regular ONT 1D ligation protocols, but R2C2 reads are also much more accurate. Because they are more accurate, R2C2 reads are also more efficiently demultiplexed and aligned than ONT 1D reads. Further, because they are similar in accuracy to Illumina reads, standard Illumina tools, like STAR and kallisto, can be used to analyze them. The gene expression and splice junction values generated by R2C2 are highly similar to those generated by Illumina MiSeq data from the same libraries.

Evaluating R2C2 for the sequencing of Illumina ChIP-seq libraries

Next, we tested the ability of R2C2 for the quality control of Illumina ChIP-seq libraries. To do this, we converted a previously generated ChIP-seq library targeting the H3K4me3 histone modification in a Glycine max (soybean) sample. The H3K4me3 library and its corresponding control Input library have previously been sequenced on an Illumina NovaSeq 6000 to a depth of 8,413,865 and 32,377,813 2 × 150-bp paired-end reads, respectively (Table 4). Based on their alignment, the sequenced molecule libraries had an insert length of 390 bp (H3K4me3) and 312 bp (Input) (Table 4).

Because the H3K4me3 and Input libraries were prepared with only a single index distinguishing them, we converted the libraries separately with R2C2 using distinct DNA splints that contained unique index sequences. This added an extra level of indexing to minimize concerns of potential index cross talk. We splint-indexed and pooled the H3K4me3 and Input ChIP-seq Illumina libraries and sequenced the pool on a single ONT MinION flowcell. We then demultiplexed the resulting R2C2 reads, assigning 2,493,021 and 1,530,914 reads (1.6:1) to the H3K4me3 and Input libraries (Table 4), respectively, a ratio that corresponded well with the 1.35:1 ratio at which they were pooled before sequencing. Importantly, the demultiplexing script scored only 163,489 (3.9%) reads as “undetermined” and assigned only 4014 (0.1%) reads to a combination of indexes not present in the library. This indicated that the extra level of indexing was highly successful in minimizing index hopping.

The demultiplexed R2C2 reads showed median read accuracy of 99.23% (H3K4me3) and 98.8% (Input) as well as median read length of 556 bp (H3K4me3) and 459 bp (Input) (Table 4). Molecules sequenced by R2C2 were therefore longer than molecules sequenced by the Illumina NovaSeq 6000 (Fig. 3A). The difference between the technologies is likely owing to the bias of the Illumina NovaSeq toward shorter molecules.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Sequencing ChIP-seq libraries on the ONT MinION after R2C2 conversion. (A) Insert length distribution of R2C2 and Illumina NovaSeq 6000 reads generated from the same Illumina library. (B) Percentage of reads in the R2C2, subsampled Illumina, and full Illumina data sets overlapping with H3K4me3 peaks generated from the full Illumina H3K4me3 data set using MACS2. (C) The comparison of the number of R2C2 and subsampled Illumina reads overlapping with H3K4me3 peaks is shown as scatter plots with marginal distributions shown as histograms. Pearson's r is shown at the bottom right. (D) Genome annotation, H3K4me3 peak areas, and read coverage histograms are shown for a section of the Gmax genome.

To test whether R2C2 reads could replace the same number of Illumina reads, we subsampled the Illumina sequencing data to the depth of the R2C2 data for both samples. We then aligned both Illumina NovaSeq 6000, subsampled Illumina NovaSeq 6000, and R2C2 reads to the G. max genome (Gmax_508_v4.0) (Valliyodan et al. 2019). For alignment, we chose the short-read preset of the minimap2 (Li 2018) aligner for both the Illumina and R2C2 data. We then called peaks on the full H3K4me3 Illumina NovaSeq 6000 data set using MACS2 and tested whether both the subsampled Illumina NovaSeq 6000 and R2C2 data could be used to evaluate the success of a ChIP experiment. Visual inspection of the data using the Phytozome JBrowse genome browser (Goodstein et al. 2012) and our own tools (Fig. 3D) showed that the subsampled Illumina NovaSeq 6000 and R2C2 data both showed the same enrichment patterns as the full Illumina NovaSeq 6000 data. A systematic analysis showed that 84% of the R2C2 reads and 69% of the subsampled Illumina reads overlap with an H3K4me3 peak identified on the full Illumina data, whereas only 18% and 11% of the respective Input reads do so (Fig. 3B).

To investigate this discrepancy in percentage of reads overlapping with H3K4me3 peaks, especially for the H3K4me3 library, we focused on differences between the R2C2 and Illumina sequencing reads. The most obvious difference is the read length with the Illumina reads originating from much shorter molecules (or library inserts). Indeed, when we recalculated this read percentage for Illumina reads originating from inserts >450 nt, it increased to 76%. Next, we analyzed the GC content of Illumina and R2C2 reads and found that—in contrast to all other experiments in this paper (Supplemental Fig. S2)—Illumina reads had a lower GC content than R2C2 reads (39% vs. 42%). To see whether the difference in insert length and GC content together would explain the discrepancy in percentage of reads overlapping with H3K4me3, we again recalculated this read percentage only for Illumina reads originating from inserts >450 nt and with a GC content >39%, that is, reads derived from long and GC-rich molecules. Here, we found that this read percentage increased to 83.2%, virtually matching the R2C2 percentage. Ultimately, this suggested that R2C2 sampled longer and slightly more GC-rich molecules from the ChIP-seq libraries. Although it is not clear why the longer molecules are more likely to overlap with H3K4me3 peaks, these peaks happen to be more GC-rich than the rest of the genome (40% vs. 30%), explaining why more GC-rich molecules are more likely to overlap with H3K4me3 peaks.

To compare whether the R2C2 and subsampled Illumina NovaSeq 6000 data sets are also similar quantitatively, we counted how many reads for each of the data sets fell into each H3K4me3 peak we identified using the full Illumina NovaSeq 6000 data set and MACS2. We found that the peak depths are correlated (Pearson's r = 0.776) (Fig. 3C). This correlation is increased to r = 0.866 when this analysis was performed with the longer/more GC-rich subsample of Illumina reads but remained lower than what we observed with the RNA-seq data. This means that although R2C2 can be used to evaluate whether a ChIP-seq experiment successfully enriched targeted chromatin, in this particular experiment, R2C2 sampled a different population of molecules than the Illumina NovaSeq 6000, thereby complicating quantitative comparisons.

Evaluating R2C2 for the sequencing of size-selected Illumina Tn5 libraries

In contrast to the other parts of the paper, which represent head-to-head comparisons between R2C2- and Illumina-based sequencing of the same short-read libraries, here, we tested whether the ability of R2C2 to sequence “medium-length” molecules >600 nt could aid in small genome assembly tasks. Illumina library preparation methods like Tn5-based tagmentation can generate library molecules >600 nt, which are too long to be sequenced efficiently by Illumina sequencers but can be efficiently processed and sequenced using R2C2. To generate these medium-length molecules for the purpose of genome assembly, we chose to size-select a Tn5-based Illumina library for molecules between 800 and 1200 bp in length, corresponding to genomic DNA inserts of ∼600–1000 bp. We then R2C2-converted and sequenced this size-selected library on the ONT MinION.

For this test, we chose to sequence the 1.2-Mb genome of the Wolbachia bacterial endosymbiont of D. melanogaster and prepared Tn5 libraries from DNA extracted from Wolbachia-containing D. melanogaster S2 cells. We generated a total of 3,338,280 R2C2 consensus reads with a median length of 680 bp. Out of these reads, we assembled 879,303 reads that did not align to the D. melanogaster genome. We used miniasm (Li 2016) for this assembly task and polished the resulting assembly using Medaka (v.1.4.4; https://github.com/nanoporetech/medaka). The resulting assembly contained 95 contigs, which covered 97.2% of the Wolbachia genome (Fig. 4A,B), and had a NGA50 of 29,963 bp and 8.5/5.6 mismatches/indels per 100 kb of sequence.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Comparing R2C2 and Illumina based assemblies of a small genome. Illumina 2 × 150 reads were assembled in 134 contigs using Meraculous. R2C2 reads were assembled using miniasm into 95 contigs. The alignments of the contigs of both assemblies—(A) Illumina and (B) R2C2—are shown as dot plots generated by MUMmer (Kurtz et al. 2004). Both approaches failed to assemble a section of the Wolbachia genome that contains pseudogenes and a transposable element near coordinate 500,000.

We also generated an assembly from Illumina NextSeq 2 × 150 bp reads from a non-size-selected Tn5 library of the same cell line. From 2,552,018 2 × 150-bp Illumina reads, we extracted 779,206 reads that did not align to the D. melanogaster genome and assembled those reads using Meraculous (Chapman et al. 2011). The resulting assembly contained 136 contigs that covered 91.6% of the Wolbachia genome (Fig. 4) and had a NGA50 of 23,217 bp and 0.5/0.6 mismatches/indels per 100 kb of sequence. Neither assembly had misassemblies as determined by QUAST (Gurevich et al. 2013).

Comparing Illumina and R2C2 assemblies of the Wolbachia genome (NC_002978.6) showed R2C2 can generate more contiguous and complete assemblies from the same library type (Fig. 4A,B). However, systematic errors produced by the ONT MinION cannot be fully removed by the R2C2 consensus process or Medaka polishing. The assembly we generated, therefore, does have more mismatches and indel errors than its Illumina counterpart. This ultimately suggests that when limited to a single Tn5 library owing to sample constraints, R2C2 can be a valuable addition to an assembly effort, but depending on use-case, further polishing with Illumina data might be required to achieve the desired base accuracy.

Evaluating R2C2 for the sequencing of target-enriched Illumina Tn5 libraries

We tested the ability of R2C2 to evaluate target-enriched Tn5 libraries and benchmark our ability to detect germline variants in the resulting data. To this end, we generated dual-indexed Tn5 libraries from genomic DNA of two cancer cell lines (NCI-H1650 and NCI-H1975) with known mutations in the EGFR gene. We pooled these libraries and enriched the pool for a panel of cancer genes based on the Stanford solid tumor STAMP panel (Newman et al. 2014) using a Twist Bioscience oligos panel and reagents (Supplemental Table S2). We performed this enrichment experiment once, without optimization, and using custom blocking oligos, therefore expecting enrichment to be far from optimal. To compare R2C2 and Illumina MiSeq, we sequenced these enriched Tn5 libraries on (1) a multiplexed Illumina MiSeq 2 × 300-bp paired-end run and (2) an ONT MinION after R2C2 conversion.

The multiplexed MiSeq run generated 7,430,624 read pairs for the NCI-H1650 library and 1,142,187 read pairs for the NCI-H1975 library. The ONT MinION run generated 3,825,657 R2C2 reads after C3POa processing. Demultiplexing then assigned 2,057,155 (53.7%) R2C2 reads to the NCI-H1650 library and 1,021,758 (26.7%) R2C2 reads to NCI-H1975. Although 537,997 (14.1%) R2C2 reads were not assigned to any sample, only 5.4% of reads were assigned to one of the two combinations of Illumina indexes not included in the pool, implying that only 0.29% (5.4% × 5.4%) of reads were assigned to the wrong sample in our dual indexed library.

After demultiplexing, we compared the insert length and target enrichment across samples and methods. We did so by merging the Illumina MiSeq read pairs using bbmerge (Bushnell et al. 2017). As with the ChIP-seq experiment, R2C2 data showed longer insert lengths than the Illumina MiSeq data, with the R2C2 insert length more closely resembling the actual length of the input library (Fig. 5A,D; Supplemental Fig. S3). We aligned the reads of different samples and methods to the human genome using the short-read preset of minimap2 and determined the percentage of reads overlapped with a target region and the coverage for each region. For NCI-H1650, 15.8% of R2C2 reads and 14.4% of Illumina MiSeq reads overlapped with a target region, producing a median coverage of 128 (fifth percentile: 28; 95th percentile: 310) for R2C2 and 558 (fifth percentile: 134; 95th percentile: 1220) for Illumina MiSeq. For NCI-H1975, 18.5% of R2C2 reads and 16.8% of Illumina MiSeq reads overlapped with a target region with a median coverage of 69 (fifth percentile: 13; 95th percentile: 166) for R2C2 and 110 (fifth percentile: 23; 95th percentile: 225) for Illumina MiSeq. The per-base coverage of the R2C2 and Illumina MiSeq data sets was very well correlated within samples, with NCI-H1650 showing a Pearson's r = 0.91 and NCI-H1975 showing a Pearson's r = 0.89 (Fig. 5B,E).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Evaluating target-enriched Tn5 libraries with R2C2. (A,D) Insert length of library molecules sequenced by Illumina or R2C2 approaches. (B,E) Comparison of per-base coverage in the Illumina and R2C2 data sets. Marginal distributions are log₂ normalized. (C,F) Alignment-based read position–dependent accuracy shown for the indicated sequencing reads and methods. (G,H) Sequencing coverage plot of the target-enriched Tn5 libraries for R2C2 and Illumina results at Chromosome 7: 55,134,584–55,211,629, which covers a part of the EGFR gene. Top panel shows the annotation of one EGFR isoform. The x-axis of the coverage plot is the base pair position, and the y-axis is the total number of reads at each position. The dotted lines indicate zoomed-in views of exons that contain the 15-bp deletion in NCI-H1650 (left) and the C-to-T and T-to-G point mutations in NCI-H1975 (right). Both samples’ Illumina reads and the R2C2 read alignments of the selected regions are shown. The mismatches are colored based on the read base (A, orange; T, green; C, blue; G, purple).

Next, we used the read alignments to determine per-base accuracy levels for all samples and method combinations. The NCI-H1975 sample, which also produced fewer reads than expected on the Illumina MiSeq, produced reads at a lower-than-expected accuracy. Read alignments suggested that the average per-base accuracies for read 1 and read 2 in NCI-H1975 were 96.81% and 98.26% compared with 98.37% and 97.88% for NCI-H1650. As expected, the per-base accuracy was highly position dependent and declined with increasing sequencing cycle number (Fig. 5C,F). Furthermore, the actual accuracy of the MiSeq reads is likely even lower owing to alignments not being extended once the read and genome are too dissimilar. The accuracy of the R2C2 reads in both NCI-H1975 and NCI-H1650 was similar and stable throughout the reads at 98.40% and 98.28%, meaning that, in this case, the R2C2 reads had a higher per-base accuracy than the combined MiSeq reads.

Visualizing Illumina MiSeq and the R2C2 read alignments showed that both methods successfully enriched for (Fig. 5G) and detected the 15-bp heterozygous deletion in the EGFR gene in the NCI-H1650 cell line and the C-to-T heterozygous variants in the EGFR gene in the NCI-H1975 cell line (Fig. 5H). To systematically evaluate the germline variant detection ability of Illumina MiSeq and R2C2 reads, we used DeepVariant (Poplin et al. 2018) for calling germline variants based on the Illumina MiSeq data and used PEPPER-Margin-DeepVariant (Shafin et al. 2021), a variant caller designed for ONT data sets, for calling germline variants in the R2C2 sequencing results. Because of the poor sequencing performance of the Illumina MiSeq for the NCI-H1975 library, we only performed this analysis on NCI-H1650. For NCI-H1650, Illumina/DeepVariant detected 119 variants in the enriched genomic regions when using a QUAL cut-off of 33.3 or more. R2C2/PEPPER-Margin-DeepVariant detected 122 variants in the enriched genomic regions when using a QUAL score of 3.8 or more, including 117 of the 119 Illumina/DeepVariant calls. When we used Illumina/DeepVariant variants as a ground truth, the R2C2/PEPPER-Margin-DeepVariant method achieved 95.9% precision and 98.3% recall.

When we visualized the reads on which the false-positive and false-negative R2C2/PEPPER-Margin-DeepVariant variant calls were made (Supplemental Fig. S4), we found that the false-positive variants were supported by fewer than half of the R2C2 reads. Moreover, when we colored the reads based on the direction of their raw reads, we found that false-positive variants were supported only by reads originating from one raw read direction. We hypothesized that if we oriented reads using the direction of their raw reads—instead of using the p5 and p7 adapters on their ends—before variant calling, it would more closely resemble regular ONT reads and provide more useful information to PEPPER-Margin-DeepVariant. Indeed, when reanalyzing the reoriented reads and using a QUAL score of nine or more, PEPPER-Margin-DeepVariant detected 116 variants that were all present in the Illumina/DeepVariant calls. This means that reorienting the reads before variant calling eliminated all false positives in the R2C2/PEPPER-Margin-DeepVariant variant calls. Reflecting known systematic errors of ONT sequencers, two of the three false negatives missing from the R2C2/PEPPER-Margin-DeepVariant variant calls were a deletion (TA → T) next to a 13-nt A homopolymer at Chr 17: 7,667,260 and a variant (G → C) next to an 8-nt C homopolymer at Chr 12: 120,994,314. The third missing variant, a G → A call at Chr 5: 112,839,666, had a 46% frequency in both Illumina and R2C2 reads and was initially identified as a candidate by PEPPER-Margin-DeepVariant but was ultimately scored as a “RefCall,” not a variant. Overall, reorienting the reads by raw read direction before running PEPPER-Margin-DeepVariant increased precision to 100% while achieving a recall of 97.4%.

This showed that R2C2 can accurately quantify what percentage of molecules in an enriched Tn5 Illumina library overlap with a target region. Despite showing longer insert lengths than the Illumina MiSeq data set, the R2C2 data set showed per-base coverage that was highly correlated with the Illumina MiSeq data. In this experiment, R2C2 actually showed a higher average per-base accuracy than the Illumina MiSeq. After reorienting R2C2 reads, variants called based on the R2C2 and Illumina MiSeq data were very similar. This shows the promise of variant calling based on ONT data but also highlights that extra care has to be taken when preparing data for use in neural network–based variant callers like DeepVariant.

Real-time analysis of Illumina library metrics using PLNK

To enable the real-time monitoring of sequencing runs and the rapid evaluation of metrics of libraries sequenced in those runs, we created the computational pipeline Processing Live Nanopore Experiments (PLNK). PLNK controls real-time base-calling, raw read processing into R2C2 consensus reads, demultiplexing of R2C2 reads, and the alignment of demultiplexed R2C2 reads to a genome. Based on the resulting alignments and the user-defined regions of interest, PLNK then determines the on-target percentage and resulting target coverage for each demultiplexed sample. PLNK runs alongside a MinION sequencing run, tracking the creation of new FAST5 files and processing them individually in the order they are generated. To do this, PLNK controls several external tools: guppy5 for base-calling, C3POa for R2C2 consensus generation, a separate Python script for demultiplexing (based on splint sequences and Illumina indexes), and mappy (minimap2 Python library) for aligning reads to a provided genome (Fig. 6A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

Real-time characterization of Illumina sequencing libraries. (A) Diagram of PLNK functionality; FAST5 files processed in the order they are produced. PLNK controls guppy5 for base-calling, C3POa for consensus calling, and mappy for alignment, as well as calculates metrics based on those alignments. (B–D) Simulation of real-time analysis for enriched Tn5 (B), ChIP-seq (C), and RNA-seq (D) libraries. For each time point, panels from top to bottom show (1) the number of FAST5 files that are produced and processed, (2) the number of demultiplexed reads produced by guppy5/C3POa/demultiplexing, (3) the percentage of reads associated with each library in the sequenced pool, (4) the percentage of reads overlapping with target regions, and (5) the median read coverage of bases in the target regions.

To test whether our pipeline could keep up with ONT MinION data generation and provide real-time analysis, we simulated ONT MinION runs using FAST5 files from previously completed sequencing experiments, our Tn5, ChIP-seq, and RNA-seq data. We used the FAST5 files’ metadata to determine the time intervals at which files were generated by the MinKnow software and copied the FAST5 files to a new output directory at those intervals. We then started PLNK to monitor the generation and control the processing of FAST5 files in this new output directory. First, we simulated the real-time analysis of the target-enriched Tn5 data. Using a desktop computer and limiting PLNK to the use of eight CPU threads and two Nvidia RTX2070 GPUs, the pipeline processed sequencing data at the same rate that a single MinION produced FAST5 files. Importantly, both the library composition (percentage of demultiplexed reads assigned to either sample; NCI-H1650 and NCI-1975), as well as the percentage of reads on-target, stabilized after less than an hour and agreed very well with the numbers generated from the whole data set (Fig. 6B). Additionally, throughout the run, PLNK reported the overall coverage of target regions in real-time.

When we simulated the analysis of ChIP-seq and RNA-seq experiments, PLNK kept up with ChIP-seq but not with the RNA-seq experiment (Fig. 6C,D). Because the RNA-seq experiment produced the largest amount of data in the study, this was not unexpected. In both cases, however, library composition and on-target percent both stabilized within the first hour of sequencing and reflected the number derived from the complete data set. This means that the library composition and quality of target-enriched Tn5 libraries (as measured by reads overlapping target areas), ChIP-seq libraries (as measured by reads overlapping with peak areas, promoters, or gene bodies—depending on targeted histone mark), and RNA-seq libraries (as measured by reads overlapping with exons) can be determined with minimal sequencing time.

The bottleneck for analysis in our desktop computer setup seemed to be the guppy5-based base-calling using the slower yet most accurate “sup” base-calling configuration. Although we could use a faster, less accurate setting to keep up with even the fastest data-producing experiments, using the most accurate model means the data can be used for in-depth analysis once the run has completed and PLNK has processed all the files, without the need to re-base-call the raw data.

Overall, this suggests that PLNK can be used to monitor ONT sequencing runs in real-time. This makes it possible to stop ONT sequencing runs when the goal of an experiment is achieved. For the rapid evaluation of library pools, this could be 1 h into a run once library composition and quality metrics have stabilized. For run monitoring, this could be several hours into a run once a specific coverage of defined target regions is reached. In both cases, a run can be stopped, allowing the ONT MinION flowcell to be flushed, stored, and ultimately reused.