
==== Front
PNAS Nexus
PNAS Nexus
pnasnexus
PNAS Nexus
2752-6542
Oxford University Press US

10.1093/pnasnexus/pgae336
pgae336
Biological, Health, and Medical Sciences
AcademicSubjects/MED00010
AcademicSubjects/SCI00010
AcademicSubjects/SOC00010
PNAS_Nexus/app-bio
PNAS_Nexus/biophys-bio
R2C2 + UMI: Combining concatemeric and unique molecular identifier–based consensus sequencing enables ultra-accurate sequencing of amplicons on Oxford Nanopore Technologies sequencers
https://orcid.org/0000-0002-9320-9430
Deng Dori Z Q Department of Molecular, Cellular, and Developmental Biology, University of California Santa Cruz, Santa Cruz, CA 95064, USA

Verhage Jack Department of Biomolecular Engineering, University of California Santa Cruz, Santa Cruz, CA 95064, USA

Neudorf Celine Department of Biomolecular Engineering, University of California Santa Cruz, Santa Cruz, CA 95064, USA

https://orcid.org/0000-0001-6535-2478
Corbett-Detig Russell Department of Biomolecular Engineering, University of California Santa Cruz, Santa Cruz, CA 95064, USA

Mekonen Honey Department of Biomolecular Engineering, University of California Santa Cruz, Santa Cruz, CA 95064, USA

Castaldi Peter J Channing Division of Network Medicine, Brigham and Women's Hospital, Boston, MA 02115, USA
Division of General Internal Medicine and Primary Care, Brigham and Women's Hospital, Boston, MA 02115, USA

https://orcid.org/0000-0002-7626-6222
Vollmers Christopher Department of Biomolecular Engineering, University of California Santa Cruz, Santa Cruz, CA 95064, USA

Yooseph Shibu Editor
To whom correspondence should be addressed: Email: vollmers@ucsc.edu
Present address: Chan Zuckerberg Biohub, San Francisco, CA, USA.

Competing Interest: C.V. has submitted a patent application for the R2C2 method. P.J.C. has received grant support for Bayer and Sanofi, and he has received consulting fees from Verona Pharma.

9 2024
21 8 2024
21 8 2024
3 9 pgae33604 1 2024
29 7 2024
05 9 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of National Academy of Sciences.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

The sequencing of PCR amplicons is a core application of high-throughput sequencing technology. Using unique molecular identifiers (UMIs), individual amplified molecules can be sequenced to very high accuracy on an Illumina sequencer. However, Illumina sequencers have limited read length and are therefore restricted to sequencing amplicons shorter than 600 bp unless using inefficient synthetic long-read approaches. Native long-read sequencers from Pacific Biosciences and Oxford Nanopore Technologies can, using consensus read approaches, match or exceed Illumina quality while achieving much longer read lengths. Using a circularization-based concatemeric consensus sequencing approach (R2C2) paired with UMIs (R2C2 + UMI), we show that we can sequence an ∼550-nt antibody heavy chain (Immunoglobulin heavy chain - IGH) and an ∼1,500-nt 16S amplicons at accuracies up to and exceeding Q50 (<1 error in 100,000 sequenced bases), which exceeds accuracies of UMI-supported Illumina-paired sequencing as well as synthetic long-read approaches.

NIH 10.13039/501100012264 NIGMS 10.13039/100000057 R35GM133569 NHLBI 10.13039/100000050 R01HL124233 U01 HL089897 R01 HL147326 U01 HL089856 NCT00608764 COPD Foundation 10.13039/100008184
==== Body
pmcSignificance Statement

In this study, we benchmark the ultra-accurate sequencing of amplicons on Oxford Nanopore Technologies (ONT) sequencers. We show that using the R2C2 sample prep method, unique molecular identifiers, and the newly developed BC1 computational tool, ONT sequencers can surpass the accuracy of gold-standard Illumina sequencers typically used for amplicon sequencing. Combined with the low capital cost of ONT sequencers, this approach opens high-throughput ultra-accurate amplicon sequencing to most molecular biology labs. Further, we show that this approach can be applied to amplicons like 16S which are over 1 kb long, freeing amplicon sequencing from the read length limitations of short-read sequencers. Finally, the BC1 tool we developed for this application represents an easy-to-use and fast tool for merging reads based on their unique molecular identifiers.

Introduction

High-throughput sequencing of PCR amplicons is an increasingly essential application that enables a large range of biological, medical, and diagnostic analyses and investigations. Amplicon sequencing can provide detailed insights into the molecular diversity at targeted genetic elements. For example, IGH amplicon sequencing can be used to determine the activation state of the adaptive immune system (1), track minimal residual disease (2), or perform basic research into the human immune system (3). Similarly, 16S amplicon sequencing can be used to profile the composition of bacteria in a variety of biosamples in clinical diagnostics or basic research (4). For each, high-sequencing accuracy is crucial because each sequenced molecule can represent a unique isolate within the sample.

Amplicon sequencing approaches have evolved substantially in recent years and present a host of opportunities for improvement. Illumina sequencers provide good accuracy, but the molecule length limitations complicate the sequencing of amplicons longer than 600 bp. For example, a full-length 16S amplicon is ∼1,500 nt in length and, therefore, requires synthetic long-read approaches like LoopSeq. LoopSeq can achieve accuracies around Q42 but has a complex, proprietary library prep and does not generate many full-length molecules (5).

Long-read sequencing technologies, Oxford Nanopore Technologies (ONT) or Pacific Biosciences (PacBio), can sequence longer amplicons. Methods like PacBio CCS, and ONT-based Nano-ampliseq and INC-seq have explored concatemer-based consensus approaches for 16S sequencing (6–8) but generally only reach accuracies below Q30—far below LoopSeq accuracy. ONT reads generated using the new R10 pore chemistry have been used in conjunction with unique molecular identifiers (UMIs) to sequence 16S or full-length rRNA operon amplicons and achieve an accuracy matching LoopSeq (∼Q42) (9). For the same full-length rRNA operons, PacBio CCS and UMIs achieved an accuracy of ∼Q52 (9). PacBio CCS and UMIs have also been used by FLAIRR-seq which analyzed full-length IGH amplicons but did not empirically determine its method's accuracy (10). In both cases, consensus generation was either done by a multistep ad hoc chain of tools or a likely suboptimal approach designed for Illumina reads.

In summary, while there are several powerful approaches currently available for high-accuracy amplicon sequencing, there currently exists no single approach that (ⅰ) is mostly agnostic to amplicon length, (ⅱ) has an easy-to-implement and robust library preparation protocol, (ⅲ) has an easy-to-use and efficient computational analysis pipeline, (ⅳ) achieves accuracies up to and exceeding Q50, and (ⅴ) is cost-effective on the instrument and per-read level.

Here, we established the easy-to-use R2C2 + UMI method. For this method, we combined our established ONT-based R2C2 method (11–16) with commonly used UMI-containing Illumina-style amplicons. To generate R2C2 + UMI reads, we developed the BC1 computational tool (https://github.com/christopher-vollmers/BC1).

We applied the R2C2 + UMI method to two amplicon types: a 550-nt Illumina-style (p5/p7 adapters) IGH amplicon and an ∼1,500-nt Illumina-style 16S amplicon. We compared the resulting R2C2 + UMI data to either a 2 × 300 MiSeq flow cell (Illumina + UMI for IGH amplicon) or publicly available LoopSeq synthetic long-read data (16S amplicon). We show that R2C2 + UMI outperforms the accuracy of both, reaching Q52, thereby establishing new accuracy benchmarks for both IGH and 16S amplicon sequencing.

R2C2 + UMI in combination with the BC1 tool enables unprecedented accuracy in interpreting variation captured by amplicon sequencing techniques.

Results

To achieve the highest possible accuracy for sequencing amplicons, we combined both molecular and computational methods.

First, we used an established protocol (17) to generate Illumina-style amplicons with each original molecule labeled by a UMI. To sequence these Illumina-style amplicons on ONT long-read sequencers, we processed the libraries into high-molecular-weight DNA using the R2C2 method. R2C2 circularizes library molecules using Gibson assembly and a DNA splint complementary to p5 and p7 Illumina adapter sequences (11). It then uses rolling circle amplification (RCA) to generate long, linear concatemers, containing multiple tandem repeats of the original library molecule (Fig. 1).

Fig. 1. Combining concatemer and UMI-based consensus approaches. Illumina-style libraries containing UMIs are converted to R2C2 libraries and sequenced on ONT sequencers. C3POa splits the resulting raw reads into subreads and combines these subreads into high-accuracy R2C2 consensus reads. BC1 parses the UMIs contained in these consensus reads, groups reads with identical or highly similar UMIs and uses R2C2 consensus reads and subreads of these groups to generate R2C2 + UMI consensus reads of even higher accuracy.

Second, after sequencing this concatemeric DNA on ONT sequencers, we used the established computational C3POa (v3.2) and newly developed BC1 (v0.95) tools to generate consensus sequences for each original library molecule. C3POa parses concatemeric raw reads into subreads and generates accurate R2C2 consensus reads from these subreads (Fig. 1).

BC1 parses R2C2 consensus reads using a highly flexible syntax for the locating and parsing of UMI sequences, enabling the detection of fixed bases used as spacers or IUPAC wildcard base codes like B, D, H, V, R, or Y which, in addition to Ns, can be used to optimize UMIs for more indel-prone long-reads. BC1 then uses the R2C2 consensus reads and subreads associated with each UMI to generate and polish R2C2 + UMI consensus reads for each original molecule: R2C2 consensus reads are aligned to each other by abpoa (18) to form a preliminary consensus. Using the subreads, this preliminary consensus is then error corrected using racon (19) and polished using medaka. Because each R2C2 consensus read can be based on multiple subreads (∼4–20 subreads per R2C2 consensus read in this study) and each UMI can be covered by multiple R2C2 consensus reads (e.g. ∼20 R2C2 consensus reads for IGH amplicons in this study), R2C2 + UMI consensus reads can be covered by hundreds of subreads and consequently reach very high accuracies (Fig. 1).

IGH sequencing

We applied the R2C2 + UMI approach to ∼550 nt IGH amplicon libraries which were designed to capture most of the information contained within antibody heavy-chain transcripts. We created these libraries with R2C2, as previously described (1, 3, 17, 20–22). In short, following reverse transcription of human whole-blood RNA, dsDNA was generated using a UMI-containing primer pool for Framing Region 1 and another UMI-containing primer pool for the first exon of the isotype-defining constant regions (Table S1). This dsDNA is then amplified using indexed Illumina p5/p7 primers. This IGH Illumina amplicon library was first sequenced on an Illumina MiSeq and then processed with R2C2 using a protocol we previously developed for Illumina libraries (11). To effectively circularize the Illumina-style amplicon using Gibson assembly, we used a splint compatible with Illumina p5 and p7 sequences. Following circularization, we performed RCA with phi29 and debranched the resulting linear DNA using T7 Endonuclease I.

We then sequenced this R2C2 library on a single PromethION flow cell as part of a pool of libraries in which it represented about 30% of pooled DNA. We processed the raw reads using C3POa resulting in 8,397,298 R2C2 consensus reads, 6,263,929 of which showed p5/p7 adapters on their ends. C3POa further demultiplexed these reads, producing 1,984,245 reads for the library. Because the original Illumina library was only sequenced to read depths of 383,867 MiSeq 2 × 300 reads, we subsampled the demultiplexed R2C2 reads to this read depth for further analysis.

We then identified UMIs and combined similar UMIs (≤1 mismatch) into groups and generated UMI-based consensus reads using the presto toolkit (23) for Illumina data and the BC1 tool for R2C2 data. These separate workflows resulted in 23,785 Illumina + UMI consensus reads and 28,511 R2C2 + UMI consensus reads. The read coverage of these UMI-based consensus reads was highly similar between the R2C2 and Illumina technologies (Fig. 2A).

Fig. 2. IGH amplicon sequencing. Several characteristics are shown for an UMI-containing IGH amplicon library sequenced either on a Illumina MiSeq (Illumina + UMI) or using R2C2 on an ONT PromethION (R2C2 + UMI). A) Read coverage distribution of UMIs (log2-scaled y-axis). B, C) Isotype and V-segment usage (log2(read numbers + 1)) with Pearson's r calculated using scipy.stats python module. D) SHM location and levels, E) CDR3 length distribution, F) CDR3 abundance and co-occurrence between the two methods, and G) estimated per base accuracy are shown for R2C2 + UMI and Illumina + UMI approaches. To estimate per base accuracy for V and J segments in the presence of SHM, R2C2 + UMI reads and Illumina + UMI reads were independently processed with IgBlast. Then mismatches and Indels in R2C2 + UMI and Illumina + UMI for the same molecule—as identified by UMI—were compared with each other. A mismatch or Indel in one method but not the other was scored as a sequencing error. For the C segments, sequencing errors were detected separately for each method by aligning each molecule to the C segment reference sequences and simply counting mismatches and Indels. The gap in per base accuracy between V and J segments in the plot represents the approximate size of the CDR3 region which is a product of somatic rearrangement and for which therefore no reference sequence exists.

For subsequent analysis, we used IgBlast (24, 25) to process R2C2 + UMI and Illumina + UMI reads. IgBlast processing yielded 27,300 R2C2 + UMI and 22,620 Illumina + UMI reads that were matched to antibody gene segments. To focus only on highly accurate reads, we restricted further analysis on reads that had four or more Illumina or R2C2 reads covering their UMI—16,261 Illumina + UMI and 15,743 R2C2 + UMI reads.

Next, we used a custom script to add isotype information to the otherwise comprehensive list of information contained within IgBlast output. Then, we compared key features of IGH transcripts as determined by the Illumina + UMI and R2C2 + UMI methods. Those metrics included (ⅰ) isotype information, i.e. does an IGH transcript encode for an IgM, IgG1, IgA2, etc., antibody that have different characteristics like membrane permeability and receptor binding, (ⅱ) V-segment usage, i.e. which of the ∼40 V segment in the human germline genome was used by a B cell during somatic recombination to create a functional IGH gene, (ⅲ) somatic hypermutation (SHM), i.e. what positions in that V segment were mutated to fine-tune the affinity of the antibody the IGH transcript encodes, (ⅳ) complementarity determining region 3 (CDR3) length, i.e. how long is the semi-random CDR3 sequence which largely defines antibody specificity and results from the somatic rearrangement of V, D, and J segments, and (ⅴ) clonal composition, i.e. what are the exact sequences of these CDR3 sequences and what is their frequencies.

First, we found that isotype (Fig. 2B) and V-segment usage (Fig. 2C) were highly similar between methods. Further, we found that the SHM locations and levels within those V segments (Fig. 2D) and the length distribution of CDR3s created by somatic recombination were highly similar between methods (Fig. 2E).

Additionally, the Illumina + UMI and R2C2 + UMI methods were used to sequence aliquots of the same IGH library containing copies of the same original molecules which, based on the read coverage distribution of UMIs (Fig. 2A), were sequenced deeply by both methods.

First, this implied that the Illumina + UMI and R2C2 + UMI methods should identify identical CDR3 sequences. Indeed, the methods detected 9,864 shared CDR3s (Fig. 2F); however, R2C2 + UMI detected 1,036 unique CDR3s and Illumina + UMI detected 1,604 unique UMIs. CDR3s unique to any method suggests sequencing errors in either method changed an existing CDR3 sequence to an artifactual CDR3 sequence.

Second, this implied that the Illumina + UMI and R2C2 + UMI methods should identify identical UMIs. This allowed us to estimate Illumina + UMI and R2C2 + UMI error rates which otherwise would have been impossible, because IGH genes can be modified by SHM and therefore lack an accurate reference. By comparing the Illumina + UMI and R2C2 + UMI sequences for the same molecule—as identified by the same UMI—as well as an unmutated reference, we showed that R2C2 + UMI consensus reads had an average per-base accuracies of Q51.7 (<1 error: 100,000 sequenced bases) compared with Q37.7 for Illumina + UMI. The comparably low Illumina + UMI average per-base accuracy was surprising so we determined per-base accuracy of both methods throughout the amplicon. Illumina + UMI had comparable accuracy to R2C2 + UMI at the amplicon ends; however, R2C2 + UMI maintained this accuracy throughout the entire amplicon, while Illumina accuracy dips to Q30 in the center of the amplicon (Fig. 2G). This is likely a consequence of Illumina MiSeq 2 × 300 paired reads, featuring declining accuracies toward their ends. This declining raw accuracy seems to be severe enough that even combining four or more Illumina reads cannot create a highly accurate consensus read. The analysis of a different IGH library containing more molecules showed that this issue is even larger in amplicons with lower per molecule coverage. While R2C2 + UMI accuracy is reduced in this low coverage library as expected, Illumina + UMI coverage declines even more toward the middle of the amplicon (Fig. S1).

Overall this analysis indicated that the R2C2 preparation did not distort library composition. Further, it suggested that R2C2, with the same number of reads, generated equivalent metrics from IGH libraries as Illumina MiSeq 2 × 300 reads. Importantly, R2C2 + UMI consensus reads were 100-times more accurate than Illumina + UMI consensus reads at the same coverage in the center of ∼550 nt amplicons.

16S sequencing

Next, we evaluated the performance of R2C2 + UMI for 16S amplicons. 16S amplicons are too long to be sequenced with Illumina paired-end read approaches, because of both read and cluster generation length limitations. Because of this, we compared R2C2 + UMI reads covering full-length ∼1,500 nt 16S amplicons generated from the ZymoBIOMICS Microbial Community DNA Standard (Zymo Research D6305) to uncorrected ONT, PacBio, and LoopSeq synthetic long-read data generated from the same standard (5, 8).

To this end, we generated 16S dsDNA using previously published primers (26) modified to contain UMIs (Table S1) and then amplified this dsDNA using indexed p5/p7 Illumina adapters. To test whether BC1 can accommodate a wide range of UMI designs we specifically designed these UMIs to not feature homopolymers longer than three nucleotides by using a BDHV nucleotide pattern. This UMI design is likely incompatible with Illumina sequencing which requires nucleotide diversity at each position but will make it possible to identify and exclude UMI-containing insertion and deletion errors which are more common in ONT and PacBio data.

We sequenced these amplicons for about 18 h on a PromethION flow cell resulting in 2,545,474 total R2C2 reads—<25% of this flow cell ultimate total output of ∼15 million reads. 2,319,680 of these R2C2 reads contained p5/p7 adapters and the expected combination of sample indexes. BC1 parsed the custom designed UMIs and consolidated these R2C2 reads into 444,217 reads with unique UMIs—R2C2 + UMI reads. About 84,154 of these R2C2 + UMI reads were covered by more than 1 R2C2 read. We then used a sliding window approach, as previously described (5) to identify and remove chimeric reads and filtered the remaining R2C2 + UMI reads for read length between 1,400 and 1,600 nt, resulting in 77,800 multiread, nonchimeric, and length-filtered R2C2 + UMI reads.

A total of 36,277 of these filtered R2C2 + UMI reads had a coverage of ≥12, which we subsequently compared with 17,862 LoopSeq reads. We also included ONT 1D reads derived from R2C2 subreads and PacBio reads (8) generated from the same ZYMOBIOMICS standard to provide context to this comparison.

To compare read accuracy, we aligned reads to 16S gene reference sequences for the ZYMOBIOMICS standard using minimap2 and then determined errors for each read.

First, we counted the total number of errors produced by the four technologies. Across 10,000 randomly sampled ∼1,500 nt reads, R2C2 + UMI contained 509 errors, LoopSeq contained 1,570 errors, PacBio CCS reads contained 12,470 errors, and ONT 1D reads contained 90,110 errors.

Second, we calculated the percentage of error-free reads produced by the four technologies. About 94.95% of R2C2 + UMI reads were free of errors, followed by LoopSeq (92.54%), PacBio (40.44%), and ONT 1D reads (4.72%; Fig. 3A). The overall number of errors per 10,000 reads as well as the percentage of error-free reads meant that overall accuracies were Q44.8 for R2C2 + UMI, Q42.6 for LoopSeq, Q32 for PacBio, and Q24 for ONT 1D reads.

Fig. 3. 16S amplicon sequencing. A) Distribution of errors per read is shown for ONT raw reads, PacBio CCS reads, LoopSeq reads, and R2C2 + UMI reads. B) Estimated per base accuracy of R2C2 + UMI and LoopSeq reads. For A and B, sequencing errors were determined by aligning reads to the 16S gene reference sequences present in the ZymoBIOMICS Microbial Community DNA Standard. Alignments were then parsed to catalog numbers and positions of mismatches and Indels. C) Species composition covering 16S genes present in the ZymoBIOMICS Microbial Community DNA Standard is shown.

Third, to see whether the previously measured overall accuracies were constant throughout the amplicon, we measured position-dependent per-base accuracies of the two most accurate methods—R2C2 + UMI and LoopSeq. LoopSeq accuracy remained between ∼Q39 and Q44 throughout most of the amplicon but dipped to as low as Q30 in the first and last 50nt—possibly a result of their (proprietary) library preparation. In contrast, R2C2 + UMI accuracy never dipped below Q40, instead remaining between ∼Q43 and Q48 throughout most of the amplicon (Fig. 3B).

Overall, this suggested that using ∼20% of a PromethION flow cell run-time, R2C2 + UMI generated 2× the reads at 3× the accuracy of a publicly available LoopSeq dataset, all while capturing the overall diversity of the input sample (Fig. 3C).

Discussion

When sequencing amplicons, achieving the highest possible sequencing accuracy removes experimental uncertainty from analysis and enables more meaningful interpretation of the resulting data. This applies to the IGH and 16S amplicons sequenced here but extends to other applications like the targeted sequencing of mixed cancer/normal samples to look for rare mutations. Illumina sequencing using UMIs has been the gold-standard for amplicon sequencing (27). However, Illumina sequencers are in fact ill-suited for amplicon sequencing, struggling with read length and library complexity. Despite this, there are a lot of protocols that produce Illumina-style amplicon libraries which would represent a significant resource if the limitations of Illumina sequencers could be overcome. We have previously shown that, using R2C2, <100 ng of Illumina-style libraries can be converted and sequenced accurately and cost-effectively on ONT sequencers, immediately overcoming restrictive library complexity and length limitations (11).

In this study, we showed that using R2C2 to convert Illumina style, UMI-containing amplicon libraries and sequencing these converted libraries with R10.4 flowcells (R2C2 + UMI) outperforms sequencing those same amplicon libraries on an Illumina MiSeq 2 × 300 (Illumina + UMI). At the same number of reads, R2C2 + UMI achieved similar consensus accuracy at the ends, but 100-fold higher consensus accuracy in the center of a ∼550 nt IGH amplicon. This highlighted that while Illumina sequencers likely perform very well with amplicons <150 nt, they struggle considerably with amplicons that push their nominal read lengths.

While effective amplicon design for Illumina sequencers is therefore limited by their read length dependent accuracy, R2C2 + UMI can sequence longer amplicons with position independent and very high accuracy.

This includes amplicons like the ∼1,500 nt 16S amplicons for which Illumina sequencers require synthetic long-read approaches like LoopSeq. Notably, R2C2 + UMI also achieved 3-fold higher accuracy than synthetic long-reads generated using the LoopSeq method when applied to 16S amplicons. We also compared R2C2 + UMI to public PacBio CCS data. However, the data were produced on PacBio SequelII without UMIs. PacBio data, produced on the more accurate PacBio Revio sequencer using the new Kinnex protocol to increase throughput and containing UMIs to further increase accuracy would likely be a highly competitive approach in terms of cost (disregarding instrument cost) and accuracy.

In summary, by achieving Q51.7 sequencing accuracies for IGH amplicons and ∼Q44.8 sequencing accuracies for 16S amplicons R2C2 + UMI has set new accuracy benchmarks for these applications.

To achieve this level of accuracy, we developed and applied the easy-to-setup-and-use BC1 computational tool with two major considerations: useability and scalability. To increase useability, BC1 was designed to require only a single command to find and parse UMIs, then to create very high-accuracy consensus reads for each UMI. To enable scalability, BC1 was designed to require only a small amount of RAM, and to be multithreaded, with the slowest part of the pipeline—polishing with medaka—being optional.

BC1 optionally uses medaka for the polishing of read sequences which, albeit slow, achieves the highest accuracies. It is however important to point out that although R2C2 + UMI/BC1 allows us to create consensus reads based on high numbers of subreads coming from both strands, R2C2 + UMI/BC1 cannot correct for highly systematic errors produced by ONT sequencers like those caused by long homopolymer tracks (>10 nucleotides). Amplicons containing these kinds of homopolymer tracks may not be an ideal application for R2C2 + UMI or really any current long-read based method.

This leaves cost as the final consideration when deciding to generate and analyze R2C2 + UMI data instead of Illumina + UMI data. R2C2 reads, when generated using a PromethION flow cell, have comparable cost to Illumina MiSeq reads—a $900 PromethION flow cell generates between 8 and 15 million R2C2 reads, an $2,500 2 × 300 MiSeq flow cell generates about 30 million paired-end reads. The more capital intensive Illumina NextSeq 2,000 now also offers ∼$3,000 2 × 300 flow cells which produce up to 100 million paired-end reads. NextSeq 2,000 reads are therefore ∼3 times cheaper than PromethION R2C2 reads; however, taking full advantage of this throughput and cost advantage will require careful planning and pooling of many libraries.

Importantly, PromethION flow cells can now be processed on the P2Solo which carries very low capital cost, making it possible for most molecular biology labs to perform these experiments in their own labs. Additionally, because C3POa and BC1 workflows can be performed on consumer-grade computer workstations, analysis of the resulting data can be performed locally as well.

As with most UMI-based approaches, R2C2 + UMI will require trial and error for every new amplicon and sample type. First, to determine the number of unique molecules in a sample that make it through the amplicon generation and R2C2 workflow and second, what R2C2 reads per unique molecule cutoff to use for BC1 to achieve optimal balance between accuracy and throughput. In this study, we found that in both IGH and 16S libraries, BC1 consensus read accuracy increased steadily until about 10 R2C2 reads per unique molecule (Fig. S2). The cutoffs of ≥4 R2C2 reads for IGH and ≥12 R2C2 for 16S were therefore trade-offs. For the deeply sequenced IGH R2C2 + UMI data, we chose a ≥4 R2C2 read cutoff because it allowed us to retain most of the ∼12,000 unique molecules in the sample without markedly decreasing overall accuracy. For 16S R2C2 + UMI data, we chose a ≥12 R2C2 read cutoff because it allowed us to retain very high overall accuracy; however, it caused us to exclude ∼50% of unique molecules from analysis. Ultimately, we showed that R2C2 read cutoffs can be chosen to achieve the throughput and accuracy requirements of a given amplicon experiment.

Overall, our results show that the R2C2 + UMI approach is an ideal choice for sequencing amplicons, combining the highest possible accuracy with long-read length, and low cost. R2C2 + UMI will therefore enable detailed biological insights by empowering ultra-accurate amplicon sequencing for a diverse array of researchers and research questions.

Materials and methods

IGH

Samples

Two de-identified blood RNA samples from the COPDGene Study were analyzed as part of this project. COPDGene is a multicenter observational study to identify genetic and multiomic biomarkers of chronic obstructive pulmonary disease (COPD). Details of this study have been published previously (28). COPDGene was approved by the institutional review boards at all participating study centers, and informed consent was obtained from all study subjects.

IGH amplicon generation

Two hundred nanograms of RNA were extracted from whole human blood and used as input for reverse transcription using Smartscribe RT (Takara Clontech) with primers targeting the IGH constant regions (IDT; Table S1; 72°C for 3 min, 42°C for 60 min and 70°C for 5 min). A nested set of constant region-specific primers and a set of primers V-segment-specific primers, both of which contain UMIs and partial p5 or p7 sequences were then used for a two-cycle PCR using Phusion polymerase (Thermo Fisher) to uniquely label each original IGH transcript (2 cycles of 98°C for 2 min, 52°C for 2 min, and 72°C for 4 min). The labeled molecules were then amplified using primers attaching complete p5 and p7 Illumina adapters to the molecules using KAPA HiFi HotStart ReadyMix (Roche; 95°C for 3 min followed by 25 cycles of 98°C for 20 s, 67°C for 15 s, and 72°C for 1 min, and last elongation at 72°C for 5 min).

IGH amplicon sequencing

UMI-labeled IGH amplicons were sequenced on multiplexed 2 × 300 flow cells on the Illumina MiSeq or converted to R2C2 libraries using the Illumina but with Nanopore workflow, as described previously (11) and then sequenced on the ONT PromethION using LSK114 ligation kits and R10.4 flow cells.

IGH amplicon analysis

Illumina data

Paired-end 2 × 300 Illumina reads were processed into full-length amplicon sequences using the pResto group of tools with the exception of UMI clustering, which was performed using our own IGH processing pipeline (https://github.com/christopher-vollmers/IGH_pipeline) to keep more consistent with the UMI clustering of the R2C2 data. For detailed commands, see Supplementary material.

R2C2 data

ONT reads were basecalled using the R10.4 SUP model of guppy(v6). Raw reads were then processed and demultiplexed into R2C2 consensus reads and their subreads using C3POa. These R2C2 consensus reads and subreads were then processed into R2C2 + UMI consensus reads using the following BC1 command:

python3 bc1.py -i IGH_C3POa.fasta -u 5.0:12.GACAG.NNNNNNNNNNNNNN.,3.0:12.GACAG.NNNNNNNNNNNNNN. -o Sample1.bd1 -s R2C2_subreads.fastq -f -t 60 -m

Error analysis

R2C2 + UMI and Illumina + UMI reads were processed using IgBlast (v1.14). IgBlast assigned each read to V, D, and J segments (as downloaded from IMGT) and determined mismatches, and indels between the read and assigned segment. Separately, R2C2 + UMI and Illumina + UMI reads were aligned to constant region exons using mappy to determine isotype and mismatches and indels in that region.

For error analysis within the constant regions, counting mismatches/indels were simply counted for each read. For error analysis within V and J segments—to take SHM into account, which is indistinguishable from sequencing error—R2C2 + UMI and Illumina + UMI reads were matched into pairs based on their identical UMIs. For each pair, mismatch and indel positions as determined by IgBlast for R2C2 + UMI and Illumina + UMI reads were compared between the reads. If one of the reads contained a mismatch or indel at a specific position and the other read did not, a sequencing error for the first technology was declared.

To adjust for read depth when visualizing position-dependent error rates, we calculated error rates for each method at each position and then simulated the number of errors likely to occur in exactly 100,000 bases at those error rates.

16S

16S amplicon generation

One hundred nanograms of ZymoBIOMICS Microbial Community DNA Standard (Zymo Research D6305) were used as template for two-cycle PCR (Phusion polymerase; Thermo Fisher) with 16S primers capturing the majority of the 16S gene and containing UMIs and partial p5 or p7 adapters (Table S1). The resulting UMI-labeled 16S dsDNA was then amplified using indexed p5/p7 primers and KAPA HiFi HotStart ReadyMix (Roche; (95°C for 3 min followed by 25 cycles of 98°C for 20 s, 67°C for 15 s, and 72°C for 1 min, and last elongation at 72°C for 5 min).

16S amplicon sequencing

UMI-labeled 16S amplicons were then converted to R2C2 libraries using the Illumina but with Nanopore workflow described previously with the modification that the splint was generated through the simple annealing of two oligos (Table S1) instead of primer extension and the sequencing using LSK114 ligation kits on R10.4 flow cells on the ONT P2Solo.

16S amplicon analysis

ONT reads were basecalled using the R10.4 SUP model of guppy(v6). Raw reads were then processed and into R2C2 consensus reads and their subreads using C3POa. These R2C2 consensus reads and subreads were then processed into R2C2 + UMI consensus reads using the following BC1 command:

python3 bc1.py -i 16S_C3POa.fasta -o 16S.BC1_consensus -s R2C2_Subreads.fastq -f -t 60 -u 5.38:50.ACAG.BDHVBDHVBDHV.AG,3.38:50.ACAG.BDHVBDHVBDHV.CG

PacBio reads from Callahan et al. (8) were downloaded from the SRA at SRR9089357.

LoopSeq reads from Callahan et al. (5) were downloaded from the SRA at SRR12148161.

Both PacBio and LoopSeq reads were filtered for perfect matches to the primers used in their respective studies and the primers were then trimmed.

Error analysis

Reads were aligned to the 16S genes extracted from the complete reference genomes of the bacteria present in the ZymoBIOMICS Microbial Community DNA Standard (Zymo Research D6305) using minimap2. Total errors were extracted from the resulting SAM file. Position-dependent errors were determined by converting the resulting SAM file into pairwise alignments using sam2pairwise (29) which were then parsed to catalog errors and their positions.

To adjust for read depth when visualizing position-dependent error rates, we calculated error rates for each method at each position and then simulated the number of errors likely to occur in exactly 100,000 bases at those error rates.

General

Samtools (30), minimap2 (31), numpy (32), scipy (33), matplotlib (34), and sam2pairwise (29) were used throughout the study.

Supplementary Material

pgae336_Supplementary_Data

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

This work was supported by the National Institutes of Health (NIH)/National Institute of General Medical Sciences (NIGMS) grant R35GM133569 to C.V. This work was supported by National Institutes of Health (NIH)/National Heart, Lung, and Blood Institute (NHLBI) R01HL124233, U01 HL089897, R01 HL147326, and U01 HL089856. The COPDGene study (NCT00608764) is also supported by the COPD Foundation through contributions made to an Industry Advisory Board that has included AstraZeneca, Bayer Pharmaceuticals, Boehringer Ingelheim, Genentech, GlaxoSmithKline, Novartis, Pfizer, and Sunovion.

Preprints

This manuscript was posted on a preprint: https://doi.org/10.1101/2023.08.19.553937.

Data Availability

Fully processed and filtered IGH amplicon data are available from the Sequence Read Archive (SRA) at PRJNA1005293. Fully processed and filtered 16S amplicon data are available from the SRA at PRJNA1005263. 16S sequence references for the ZymoBIOMICS Microbial Community DNA Standard (Zymo Research D6305) are provided as Supplementary Data File S1. The BC1 tool is available at github at https://github.com/christopher-vollmers/BC1. Tools for trimming sequencing primers and filtering chimeras used in the study can also be found in that repository.
==== Refs
References

1 Vollmers  C, et al  2015. Monitoring pharmacologically induced immunosuppression by immune repertoire sequencing to detect acute allograft rejection in heart transplant patients: a proof-of-concept diagnostic accuracy study. PLoS Med.  12 :e1001890.26466143
2 Logan  AC, et al  2011. High-throughput VDJ sequencing for quantification of minimal residual disease in chronic lymphocytic leukemia and immune reconstitution assessment. Proc Natl Acad Sci U S A. 108 :21194–21199.22160699
3 Horns  F, et al  2016. Lineage tracing of human B cells reveals the in vivo landscape of human antibody class switching. Elife. 5 :e16578.27481325
4 Kim  J, Atkinson  C, Miller  MJ, Kim  KH, Jin  Y-S. 2023. Microbiome engineering using probiotic yeast: and the secreted human lysozyme lead to changes in the gut microbiome and metabolome of mice. Microbiol Spectr. 11 (4 ):e0078023.37436157
5 Callahan  BJ, Grinevich  D, Thakur  S, Balamotis  MA, Yehezkel  TB. 2021. Ultra-accurate microbial amplicon sequencing with synthetic long reads. Microbiome. 9 :130.34090540
6 Li  C, et al  2016. INC-Seq: accurate single molecule reads using nanopore sequencing. Gigascience. 5 :34.27485345
7 Calus  ST, Ijaz  UZ, Pinto  AJ. 2018. NanoAmpli-Seq: a workflow for amplicon sequencing for mixed microbial communities on the nanopore sequencing platform. Gigascience. 7 :giy140.30476081
8 Callahan  BJ, et al  2019. High-throughput amplicon sequencing of the full-length 16S rRNA gene with single-nucleotide resolution. Nucleic Acids Res. 47 :e103.31269198
9 Karst  SM, et al  2020. Enabling high-accuracy long-read amplicon sequences using unique molecular identifiers with Nanopore or PacBio sequencing. bioRxiv 645903. 10.1101/645903.
10 Ford  EE, et al  2023. FLAIRR-seq: a method for single-molecule resolution of near full-length antibody H chain repertoires. J Immunol. 210 :1607–1619.37027017
11 Zee  A, et al  2022. Sequencing Illumina libraries at high accuracy on the ONT MinION using R2C2. Genome Res. 32 :2092–2106.36351772
12 Volden  R, Vollmers  C. 2022. Single-cell isoform analysis in human immune cells. Genome Biol. 23 :47.35130954
13 Volden  R, et al  2018. Improving nanopore read accuracy with the R2C2 method enables the sequencing of highly multiplexed full-length single-cell cDNA. Proc Natl Acad Sci U S A.  115 :9726–9731.30201725
14 Cole  C, Byrne  A, Adams  M, Volden  R, Vollmers  C. 2020. Complete characterization of the human immune cell transcriptome using accurate full-length cDNA sequencing. Genome Res. 30 :589–601.32312742
15 Byrne  A, et al  2019. Depletion of hemoglobin transcripts and long read sequencing improves the transcriptome annotation of the polar bear (Ursus maritimus). Front Genet. 10 :643. 10.3389/fgene.2019.00643 31379921
16 Vollmers  AC, Mekonen  HE, Campos  S, Carpenter  S, Vollmers  C. 2021. Generation of an isoform-level transcriptome atlas of macrophage activation. J Biol Chem. 296 :100784. 10.1016/j.jbc.2021.100784 34000296
17 Vollmers  C, Sit  RV, Weinstein  JA, Dekker  CL, Quake  SR. 2013. Genetic measurement of memory B-cell recall using antibody repertoire sequencing. Proc Natl Acad Sci U S A.  110 :13463–13468.23898164
18 Gao  Y, et al  2021. abPOA: an SIMD-based C library for fast partial order alignment using adaptive band. Bioinformatics. 37 (15 ):2209–2211. 10.1093/bioinformatics/btaa963 33165528
19 Vaser  R, Sović  I, Nagarajan  N, Šikić  M. 2017. Fast and accurate de novo genome assembly from long uncorrected reads. Genome Res. 27 :737–746.28100585
20 Cole  C, Volden  R, Dharmadhikari  S, Scelfo-Dalbey  C, Vollmers  C. 2016. Highly accurate sequencing of full-length immune repertoire amplicons using Tn5-enabled and molecular identifier–guided amplicon assembly. J Immunol. 196 :2902–2907.26856699
21 Horns  F, Vollmers  C, Dekker  CL, Quake  SR. 2019. Signatures of selection in the human antibody repertoire: selective sweeps, competing subclones, and neutral drift. Proc Natl Acad Sci U S A.  116 :1261–1266.30622180
22 de Bourcy  CFA, et al  2017. Phylogenetic analysis of the human antibody repertoire reveals quantitative signatures of immune senescence and aging. Proc Natl Acad Sci U S A. 114 :1105–1110.28096374
23 Vander Heiden  JA, et al  2014. pRESTO: a toolkit for processing high-throughput sequencing raw reads of lymphocyte receptor repertoires. Bioinformatics. 30 :1930–1932.24618469
24 Ye  J, Ma  N, Madden  TL, Ostell  JM. 2013. IgBLAST: an immunoglobulin variable domain sequence analysis tool. Nucleic Acids Res. 41 :W34–W40.23671333
25 Lefranc  M-P, et al  2004. IMGT-ONTOLOGY for immunogenetics and immunoinformatics. In Silico Biol. 4 :17–29.15089751
26 Matsuo  Y, et al  2021. Full-length 16S rRNA gene amplicon analysis of human gut microbiota using MinIONTM nanopore sequencing confers species-level resolution. BMC Microbiol. 21 :35.33499799
27 Kinde  I, Wu  J, Papadopoulos  N, Kinzler  KW, Vogelstein  B. 2011. Detection and quantification of rare mutations with massively parallel sequencing. Proc Natl Acad Sci U S A.  108 :9530–9535.21586637
28 Regan  EA, et al  2010. Genetic epidemiology of COPD (COPDGene) study design. COPD. 7 :32–43.20214461
29 LaFave  MC, Burgess  SM. 2014. sam2pairwise version 1.0.0. 10.5281/zenodo.11377
30 Li  H, et al  2009. The sequence alignment/map format and SAMtools. Bioinformatics. 25 :2078–2079.19505943
31 Li  H . 2018. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 34 :3094–3100.29750242
32 Harris  CR, et al  2020. Array programming with NumPy. Nature. 585 :357–362.32939066
33 Virtanen  P, et al  2020. Scipy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods.  17 :261–272.32015543
34 Hunter  JD . 2007. Matplotlib: a 2D graphics environment. Comput Sci Eng.  9 :90–95.
