Main<\/h2>\n\nTransfer RNAs (tRNAs) are abundant small non-coding RNAs that play a pivotal role in decoding genetic information1<\/a>,2<\/a>,3<\/a><\/sup>. Based on their aminoacylation identity, tRNAs are subdivided into 20 accepting groups (alloacceptors), each comprising several tRNAs that translate synonymous codons with the same amino acid (isoacceptors). To fulfill their function as adapter molecules between the RNA and protein codes, tRNAs are extensively modified, containing on average 13 modifications per tRNA molecule4<\/a><\/sup>. Although some tRNA modifications are thought to be structural and static, others are dynamic and even reversible5<\/a>,6<\/a>,7<\/a>,8<\/a><\/sup>, affecting the fate and function of individual tRNA molecules2<\/a>,9<\/a>,10<\/a>,11<\/a>,12<\/a>,13<\/a><\/sup>. Notably, mutations in multiple tRNA modification enzymes have been associated with a wide variety of human diseases14<\/a>,15<\/a>,16<\/a>,17<\/a><\/sup>, highlighting their importance in proper cellular functioning.<\/p>\ntRNA modifications are present in all domains of life18<\/a>,19<\/a><\/sup> and are synthesized by dedicated tRNA-modifying enzymes that alter specific tRNAs in a site-specific fashion20<\/a><\/sup>. The chemical nature of these modifications is highly diverse and includes methylations, acetylations, isomerizations, deaminations and conjugation to amino acids, among others21<\/a>,22<\/a><\/sup>. Certain tRNA modifications are found only in a single tRNA species, whereas others are found in multiple tRNA species20<\/a>,23<\/a><\/sup>. For example, 2-lysidine (k2<\/sup>C) tRNA modifications are placed at position 34 of the anticodon of tRNAIle<\/sup>(AUA)24<\/a><\/sup>, whereas pseudouridine (\u03a8) can be placed at diverse positions of the tRNA molecule and in multiple tRNA isoacceptors25<\/a>,26<\/a>,27<\/a><\/sup>. Regarding their function, tRNA modifications can sometimes act as identity elements recognized by aminoacyl-tRNA synthetases28<\/a>,29<\/a>,30<\/a><\/sup>, and, without modifications, many tRNAs have poor aminoacylation capability31<\/a><\/sup>. On the other hand, tRNA modifications can affect the decoding preferences of tRNA molecules, especially those found at position 34 of the anticodon16<\/a>,32<\/a>,33<\/a>,34<\/a><\/sup>, restricting or increasing the wobbling capacity of the tRNAs and, consequently, changing the set of \u2018preferred\u2019 or \u2018optimal\u2019 codons that will be translated35<\/a>,36<\/a>,37<\/a>,38<\/a>,39<\/a>,40<\/a>,41<\/a>,42<\/a>,43<\/a>,44<\/a><\/sup>.<\/p>\nIn the last few years, it has been shown that some tRNA modifications are reversible8<\/a>,45<\/a>,46<\/a>,47<\/a><\/sup> and can be dynamically regulated upon environmental exposures48<\/a>,49<\/a>,50<\/a><\/sup>, across cell cycle stages51<\/a><\/sup> and upon tumorigenesis52<\/a>,53<\/a>,54<\/a>,55<\/a>,56<\/a><\/sup>. Similarly, tRNA abundances are also dysregulated upon environmental exposures, such as oxidative stress48<\/a>,57<\/a>,58<\/a><\/sup>, as well as in diverse types of cancer52<\/a>,59<\/a>,60<\/a>,61<\/a><\/sup>. Modulation of tRNA abundances and\/or tRNA modifications is generally regarded as a molecular strategy that allows cells to adapt to distinct physiological states or conditions, leading to increased expression of subsets of proteins that otherwise would remain poorly translated under \u2018normal\u2019 tRNA abundances60<\/a>,62<\/a>,63<\/a><\/sup>.<\/p>\nDespite the pivotal function that tRNAs play in cellular processes and their involvement in numerous human diseases, we currently lack a simple and cost-effective method to accurately quantify both tRNA abundances and their modifications systematically. On the one hand, tRNA modifications are typically identified and quantified with high accuracy using liquid chromatography coupled to mass spectrometry (LC\u2013MS) methodologies48<\/a>,64<\/a>,65<\/a>,66<\/a>,67<\/a>,68<\/a>,69<\/a><\/sup>. In these methods, RNA molecules are fragmented into oligomers or monomers, and their abundance is measured via UV absorption or MS\/MS. LC\u2013MS\/MS techniques using triple quadrupole-based detection are among the most sensitive, allowing limits of quantification in the low femtomole range48<\/a>,64<\/a>,70<\/a>,71<\/a>,72<\/a><\/sup>, but they typically cannot identify the tRNA isoacceptor that contained each detected modification. On the other hand, tRNA abundances can be determined using tRNA microarrays36<\/a>,59<\/a>,73<\/a><\/sup> and, more recently, by employing next-generation sequencing (NGS) technologies74<\/a>,75<\/a>,76<\/a>,77<\/a>,78<\/a>,79<\/a>,80<\/a>,81<\/a>,82<\/a><\/sup>, which require an initial conversion of the tRNA molecules into cDNA. Consequently, NGS-based methods are blind to most tRNA modifications, as these are typically erased during the reverse transcription step. Moreover, tRNA modifications that disrupt the Watson\u2013Crick base pairing, which are abundant in tRNAs, will interfere with the reverse transcriptase enzyme, causing it to drop off, producing truncated reads, in addition to misincorporations72<\/a>,83<\/a>,84<\/a>,85<\/a><\/sup> (Fig. 1a<\/a>). To overcome these limitations, a wide variety of improved tRNA sequencing protocols have been developed in recent years, which often include the use of highly processive reverse transcriptase enzymes75<\/a>,77<\/a><\/sup> and\/or cocktails of demethylases74<\/a>,75<\/a>,79<\/a><\/sup>. However, despite these improvements, NGS-based methods still suffer from the following caveats: (1) they introduce significant biases during the library preparation, caused by incomplete reverse transcriptions84<\/a>,86<\/a><\/sup>, incomplete demethylations87<\/a><\/sup> and polymerase chain reaction (PCR) amplification88<\/a><\/sup>, resulting in skewed representations of existing tRNA populations; and (2) they cannot detect most tRNA modifications, as these are typically \u2018erased\u2019 during the conversion of RNA to cDNA. Therefore, a simple, robust and efficient tRNA sequencing method is still needed.<\/p>\n\nFig. 1: Nano-tRNAseq can efficiently sequence both IVT and native tRNA populations.<\/b><\/figcaption>\n![\"Science](\"http:\/\/media.springernature.com\/lw685\/springer-static\/image\/art%3A10.1038%2Fs41587-023-01743-6\/MediaObjects\/41587_2023_1743_Fig1_HTML.png\")
<\/picture><\/a><\/div>\na<\/b>, Schematic of the modifications found in S. cerevisia<\/i>e cytoplasmic tRNA, shown in its usual secondary structure form with circles representing nucleotides and lines representing base pairs. Gray circles represent unmodified nucleotides; pink circles represent possible modification sites; and those with a black outline indicate modifications that cause errors during reverse transcription. Possible RNA modifications occurring at each position are listed in the surrounding boxes; modifications that cause misincorporation during reverse transcription are in green; and those that cause reverse transcription truncation are in blue. b<\/b>, Schematic overview illustrating the steps required for tRNA library preparation using Nano-tRNAseq (see Extended Data Fig. 1<\/a> for more details). c<\/b>, IGV snapshots of Nano-tRNAseq mapped reads from synthetic IVT tRNAs (upper panels) or biological tRNAs (lower panel). Positions with a mismatch frequency greater than 0.2 are colored, whereas those showing mismatch frequencies lower than 0.2 are shown in gray. d<\/b>, Scatter plot of tRNA abundances showing the replicability of Nano-tRNAseq when WT S. cerevisia<\/i>e tRNA biological replicates are sequenced. The correlation strength is indicated by Spearman\u2019s correlation coefficient (\u03c1). RT, reverse transcription.<\/p>\n<\/div>\nFull size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\nA promising alternative to the use of NGS-based technologies to characterize the tRNAome is the direct RNA sequencing (DRS) platform developed by Oxford Nanopore Technologies (ONT). This technology allows direct sequencing of native RNA molecules and, as such, can, in principle, detect and measure both tRNA modifications and tRNA abundances without the need for reverse transcription or PCR. Previous works have demonstrated that nanopores can capture tRNAs using solid-state or biological (ONT) nanopores89<\/a>,90<\/a>,91<\/a>,92<\/a><\/sup>. For example, sequencing of five distinct tRNAs was achieved using solid-state nanopores89<\/a><\/sup>, and tRNAs were shown to be distinguishable from other short RNAs using the MspA pore90<\/a><\/sup>. Later studies showed that, by lengthening the tRNA molecules with ligated adapter extensions, tRNAs could be sequenced, basecalled and mapped using biological nanopores92<\/a><\/sup>. However, in these studies, the proposed approach led to low sequencing yields of tRNA molecules (~20\u201340\u00d7 lower than expected for a DRS run) and did not report whether extant in vivo tRNA abundances and\/or tRNA modifications were recapitulated using this method.<\/p>\nHere we present Nano-tRNAseq, a nanopore-based approach that allows accurate and direct measurement of native tRNA molecules. The library preparation benefits from the 3\u2032 CCA overhang typically present in the mature tRNAs to incorporate a double ligation of RNA adapters at both the 5\u2032 and 3\u2032 ends of the tRNA molecules, which leads to an improved proportion of basecalled and mapped tRNA molecules. Moreover, we show that MinKNOW, the ONT proprietary software required to run nanopore sequencing experiments, erroneously discards the majority of tRNA reads, misinterpreting them as \u2018adapters\u2019, and also causes biases in the estimated tRNA abundances due to preferential capture of longer tRNAs (for example, tRNALeu<\/sup>, tRNAArg<\/sup> and tRNASer<\/sup>). To overcome these limitations, here we provide a computational framework that allows us to capture ~10\u00d7 more tRNA reads and accurately recapitulates tRNA abundances.<\/p>\nAltogether, our work provides a simple, cost-effective, high-throughput and reproducible method to accurately quantify tRNA abundances and capture tRNA modification changes simultaneously using native tRNA nanopore sequencing, providing a framework to study the tRNAome at single-molecule resolution. We envision that Nano-tRNAseq will contribute to the study of the biological function of tRNA modifications in a wide variety of contexts, such as cancer, stress exposures or viral infection, and opens the possibility of exploiting these molecules as biomarkers of human health and disease.<\/p>\n<\/div>\n<\/div>\n
\nResults<\/h2>\n\nStandard nanopore DRS results in low tRNA sequencing yields<\/h3>\n
Nanopore DRS is a well-established long-read sequencing technology to study RNA molecules, typically polyadenylated mRNAs93<\/a>,94<\/a>,95<\/a>,96<\/a>,97<\/a><\/sup>. Although several works have shown that this technology can also be used to study short RNA molecules, such as snoRNAs and snRNAs96<\/a>,98<\/a><\/sup>, DRS is inefficient at capturing RNA molecules shorter than 200\u2009nucleotides (nt) and is generally considered unable to capture sequences shorter than ~100\u2009nt96<\/a>,97<\/a><\/sup>, limiting its applicability to study short RNA populations, such as tRNAs. In addition, the first ~15\u2009nt at the 5\u2032 end of RNA molecules are typically lost in DRS runs93<\/a>,99<\/a><\/sup>, as this portion cannot be adequately basecalled due to the increase in the RNA translocation speed when the 5\u2032 end of the molecule exits the helicase99<\/a><\/sup>. To overcome these limitations, we reasoned that extension of the 5\u2032 and 3\u2032 ends of the tRNA would lead to improved sequencing of tRNA molecules, as these would now be beyond the ~100-nt threshold, in addition to capturing the sequence and modification information of 5\u2032 tRNA ends.<\/p>\nWe first attempted a modified tRNA DRS approach in which a 5\u2032 RNA adapter, complementary to the 3\u2032 CCA overhang present in mature tRNA molecules, was ligated to the 5\u2032 end of tRNAs that had been previously in vitro polyadenylated (Strategy A; Extended Data Fig. 1<\/a> and Methods<\/a>). A set of nine synthetic in vitro transcribed (IVT) tRNAs of various lengths and sequences (Supplementary Table 1<\/a>) were sequenced using this strategy. However, this approach produced poor sequencing yields (56,002 reads; Supplementary Table 2<\/a>) compared to a standard DRS run (~1\u20132 million reads). Moreover, only 7.5% of reads mapped uniquely to tRNAs using minimap2 (ref. 100<\/a><\/sup>) with recommended parameters (-ax map -ont -k15) (Supplementary Table 2<\/a>). Relaxation of the mapping parameters (-ax map-ont -k5), which had previously been shown to improve the mappability of highly modified RNAs98<\/a><\/sup>, did not significantly increase the number of mapped tRNA reads (Supplementary Table 2<\/a>).<\/p>\nNext, we altered our library preparation protocol to replace the poly(A) tail with a 3\u2032 DNA adapter complementary to the 5\u2032 RNA adapter (Strategy B), such that the two oligonucleotides could be pre-annealed and ligated to the tRNA (Extended Data Fig. 1<\/a>). However, this strategy also yielded a low number of sequenced reads (63,502 reads) and a low percentage of uniquely mapped reads (6.5%) (Supplementary Table 2<\/a>). We speculated that the low number of reads could be due to steric interference of the poly(A) preventing 5\u2032 ligation (Strategy A) or that the 3\u2032 DNA adapter is not basecalled (Strategy B). These scenarios would lead to low coverage of the 5\u2032 or 3\u2032 ends, respectively, decreasing the mappability of reads resulting from these two strategies.<\/p>\nExtending tRNA ends with RNA adapters improves basecalling<\/h3>\n
Based on the results of Strategy A and Strategy B, we rationalized that padding the 5\u2032 and 3\u2032 tRNA ends with RNA adapters, which can be accurately basecalled and mapped, would enable us to capture the entirety of the tRNA sequence. This approach, which we termed Nano-tRNAseq (Fig. 1b<\/a>), was the most successful at sequencing, basecalling and mapping both in vitro and native tRNA molecules using nanopore DRS. We should note that a recent work also proposed a similar solution to facilitate native tRNA nanopore sequencing92<\/a><\/sup>.<\/p>\nIn the first step, a 5\u2032 RNA adapter (orange) complementary to the CCA overhang of mature tRNAs is pre-annealed to a 3\u2032 RNA adapter (red) containing three DNA bases (pink) at the 3\u2032 end (Fig. 1b<\/a>). In preliminary Nano-tRNAseq runs, we used an RNA-only 3\u2032 adapter but observed that RNA-only 3\u2032 adapters led to increased self-ligation (Supplementary Table 2<\/a>), an issue that we mitigated by adding DNA bases to the end of the adapter (Extended Data Fig. 1<\/a>). Next, the pre-annealed 5\u2032 RNA and 3\u2032 RNA:DNA splint adapters were ligated to deacylated tRNAs. Knowing that an RNA:RNA ligation with an RNA bridge has a low efficiency101<\/a><\/sup>, various ligation times and the addition of a molecular crowding agent were tested to ensure that conditions that maximized ligation efficiency were chosen (Extended Data Fig. 2a,b<\/a>). Subsequently, ONT RTA oligoA and oligoB were pre-annealed and ligated to the tRNA molecule using T4 DNA Ligase (see Extended Data Fig. 3<\/a> for validation of each ligation step). This approach resulted in >200,000 basecalled reads (Supplementary Table 2<\/a>), thus significantly increasing sequencing output by up to fourfold relative to the previous strategies, and also with improved 5\u2032 and 3\u2032 coverage of both synthetic and biological tRNAs (Fig. 1c<\/a>). Additionally, we found Nano-tRNAseq to be highly replicable when sequencing native tRNAs (\u03c1\u2009=\u20090.984) (Fig. 1d<\/a>).<\/p>\nMapping parameters significantly affect tRNA read mappability<\/h3>\n
The alignment of native tRNA reads is challenging due to their short and highly modified nature. Indeed, native tRNAs contain a large proportion of mismatched bases, often originating from inaccurate basecalling of modified bases in DRS datasets94<\/a>,98<\/a>,102<\/a><\/sup>. As a consequence of these miscalled bases, the commonly used long-read mapper minimap2 with recommended settings (-ax map-ont -k15) aligned only a fraction (2.56%) of the reads (Fig. 2a\u2013c<\/a> and Supplementary Table 2<\/a>). We further tested a range of minimap2 parameters and observed only incremental improvements in mapping and an increase in false alignments (antisense mapped reads served as a proxy of mismapping) (Supplementary Table 3<\/a>).<\/p>\n\nFig. 2: The choice of mapping software and parameters markedly affects the number of mapped tRNA reads.<\/b><\/figcaption>\n![\"Science](\"http:\/\/media.springernature.com\/lw685\/springer-static\/image\/art%3A10.1038%2Fs41587-023-01743-6\/MediaObjects\/41587_2023_1743_Fig2_HTML.png\")
<\/picture><\/a><\/div>\na<\/b>,b<\/b>, IGV snapshots of reads mapped to IVT D. melanogaster<\/i> mitochondrial tRNAAla(UGC)<\/sup> (a<\/b>) and S. cerevisiae<\/i> tRNAPhe<\/sup> (b<\/b>), mapped using different mapping algorithms (minimap2, BWA-MEM and BWA-SW) and parameters. The 5\u2032 and 3\u2032 RNA adapter regions, ligated to the ends of the tRNA molecule, were included in the mapping references and are represented by an orange bar and a red bar, respectively. Positions with a mismatch frequency greater than 0.2 are colored, whereas those showing mismatch frequencies lower than 0.2 are shown in gray. c<\/b>, Bar plot depicting the effect of algorithm and parameter choice on the relative proportion of uniquely mapped reads (green) and mismapped reads (purple; reads mapping to antisense strands were used as a proxy to assess mismapping) (Supplementary Table 4<\/a>). The proportion of mapped reads (d<\/b>) and alignment identity (e<\/b>) for each template from the bar plot in c<\/b>, using either minimap2 or bwa mem -W13 -k6- T20. We should note that minimap2 alignment identity in S. cerevisiae<\/i> tRNAPhe<\/sup> was not computed because no reads were mapped to this tRNA using minimap2 with -ax map-ont -k15 parameters (Supplementary Table 5<\/a>). f<\/b>, Bar plot showing the effect of trimming the length of the 5\u2032 RNA adapter (reds) and 3\u2032 RNA adapter (blues) on tRNA read mappability (Supplementary Table 6<\/a>). The conditions used by Nano-tRNAseq are gray, whereas the effect of not using RNA adapters is black.<\/p>\n<\/div>\nFull size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\nTo improve the mappability of Nano-tRNAseq reads, we next tested the mapping algorithm BWA, a short-read mapping algorithm commonly used to map Illumina reads103<\/a>,104<\/a><\/sup>. Using sequencing data from a Nano-tRNAseq run that contained three different tRNA constructs (IVT Drosophila melanogaster<\/i> mitochondrial tRNAAla(UGC)<\/sup>, IVT Streptococcus pneumoniae<\/i> tRNASer(UGA)<\/sup> and native S. cerevisiae<\/i> tRNAPhe<\/sup>; see Supplementary Table 2<\/a> for a summary of the sequencing runs in this work), we found that the BWA-MEM aligner with recommended parameters outperformed minimap2 in terms of proportion of mapped reads (Fig. 2c<\/a>), in agreement with recent works92<\/a><\/sup>. Although more relaxed configurations of BWA-MEM aligned more reads, this also came at the expense of increased false alignments (Fig. 2c<\/a> and Supplementary Table 4<\/a>). An optimal balance between increased mapped reads and false alignments was found when using bwa mem with parameters -W13 -k6 -xont2d -T20, which mapped 54.63% of the reads, with very few false alignments (0.19%). When comparing the performance of the mapping algorithms in native tRNA molecules, the contrast was even more stark; although minimap2 mapped IVT tRNAs, it failed to map a single biological tRNA read (Fig. 2d<\/a> and Supplementary Table 5<\/a>). The alignment identity was similar to minimap2 for reads that mapped to IVT tRNAs but was slightly lower than the typical identity obtained in nanopore DRS runs, suggesting that short reads, even without modifications, cause a drop in the basecalling accuracy (Fig. 2e<\/a>). Notably, the alignment identity of S. cerevisiae<\/i> tRNAPhe<\/sup> was lower (~74.5%) than in synthetic tRNAs (81.8%), presumably due to the presence of base modifications present on endogenously modified S. cerevisiae<\/i> tRNAPhe<\/sup> (Fig. 2e<\/a> and Supplementary Table 5<\/a>).<\/p>\nWe then assessed whether the mappability of Nano-tRNAseq reads might be affected by the length of the 5\u2032 and 3\u2032 RNA adapters. To simulate different RNA adapter lengths, we trimmed one or both adapters from the reference sequences. We found that the absence of both RNA adapters had only a modest effect on the mappability of reads originating from IVT tRNAs (decrease of 6\u201311%) (Fig. 2f<\/a> and Supplementary Table 6<\/a>), whereas their absence had a major effect in the mappability of native S. cerevisiae<\/i> tRNAPhe<\/sup> (55% loss). Thus, we concluded that short and unmodified sequences can be aligned efficiently even without RNA adapters in the reference sequences, whereas short and modified reads greatly benefit from the extension with adapters, demonstrating that extending molecules with RNA adapters is essential for guiding the correct alignment of short reads enriched in \u2018mismatches\u2019, such as those derived from native tRNAs. For native S. cerevisiae<\/i> tRNAPhe<\/sup>, the read mappability plateaued at a 3\u2032 adapter length of 25\u2009nt (Fig. 2f<\/a>), suggesting that the 30-nt RNA portion of the 3\u2032 RNA:DNA oligo used in Nano-tRNAseq is more than sufficient to achieve optimal read mappability.<\/p>\nCustom MinKNOW improves yield and tRNA abundance estimates<\/h3>\n
A surprising feature of our initial tRNA sequencing runs was the low amount of sequenced reads. Although pore clogging caused by tRNA structure might partially explain the low sequencing yield92<\/a><\/sup>, we also noticed that the MinKNOW software classified a high proportion of reads as \u2018adapter-only\u2019 reads in real time. Hence, we hypothesized that a considerable fraction of tRNA reads might be discarded by the MinKNOW software due to their short signal lengths, as they resemble \u2018adapter-only\u2019 reads.<\/p>\nThe MinKNOW software is responsible for analyzing the continuous electrical current (signal intensity) measured at each pore, reporting the signal regions that correspond to \u2018reads\u2019 into FAST5 files, which are then basecalled to generate a FASTQ file (Fig. 3a<\/a>). We noted that MinKNOW, by default, reports reads that last at least 2\u2009seconds, roughly corresponding to RNA molecules of 140\u2009nt (assuming constant helicase processivity of 70\u2009nt per second in DRS) (Extended Data Fig. 4<\/a>). Considering that canonical tRNA molecules are ~73\u2009nt, this would imply that even after double RNA adapter ligation (where 24 and 30 RNA nucleotides are added to the 5\u2032 and 3\u2032 ends of the tRNA molecule, respectively), the size of the ligated tRNA molecule would still be below the threshold, possibly leading to misassignment of tRNA reads as \u2018adapter-only\u2019 reads. To alleviate this issue, we tested whether alternative MinKNOW configurations would improve the classification of tRNA reads and boost sequencing yields. To this end, the bulk \u2018raw\u2019 dump files were saved during the sequencing and were reprocessed using alternative MinKNOW configurations (Supplementary Table 7<\/a>).<\/p>\n\nFig. 3: Adjustment of MinKNOW parameters increases the number of sequenced and mapped tRNA reads.<\/b><\/figcaption>\n![\"Science](\"http:\/\/media.springernature.com\/lw685\/springer-static\/image\/art%3A10.1038%2Fs41587-023-01743-6\/MediaObjects\/41587_2023_1743_Fig3_HTML.png\")
<\/picture><\/a><\/div>\na<\/b>, MinKNOW software classifies continuous current passing through pores as open pore, adapter or strand (actual reads) and outputs fragments classified as strand to a FAST5 file, which are then basecalled to generate a FASTQ file. b<\/b>, Diagram showing the conceptual difference between default and custom MinKNOW read classification (Extended Data Fig. 4<\/a>). c<\/b>, Bar plot of sequencing yield in terms of basecalled and uniquely mapped reads obtained with default and custom configurations (Supplementary Table 7<\/a>). d<\/b>, Scatter plot of the relative fold change of uniquely mapped reads with respect to tRNA length (Supplementary Table 8<\/a>). e<\/b>, Histogram of read count and alignment length of IVT tRNA reads captured with default and custom configurations. f<\/b>, Bar plot of the relative proportion of IVT tRNA molecules D. melanogaste<\/i>r mitochondrial tRNAAla(UGC)<\/sup> and S. pneumoniae<\/i> tRNASer(UGA)<\/sup> and native S. cerevisiae<\/i> tRNAPhe<\/sup> reads recovered with default and custom settings (Supplementary Table 9<\/a>), where the dotted line indicates the expected proportion. g<\/b>, Expected versus observed log read counts of nine IVT and one native tRNA molecules captured using the custom MinKNOW configuration (Supplementary Table 10<\/a>). Spearman correlation (\u03c1) is shown.<\/p>\n<\/div>\nFull size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\nBy lowering the MinKNOW strand minimum duration to 1\u2009second and the adapter maximum duration to 2\u2009seconds (see Extended Data Fig. 4<\/a> for schematic), a configuration that we refer to as custom<\/i>, we captured ~12-fold more basecalled and ~4.5-fold more uniquely mapped tRNA reads compared to the default MinKNOW configuration (Fig. 3b,c<\/a> and Supplementary Table 7<\/a>). Notably, we found that the default MinKNOW configuration not only led to low sequencing outputs but also caused significant biases in the relative abundances of tRNA molecules. Specifically, we found a greater representation of shorter tRNAs in our custom configuration (Fig. 3d,e<\/a> and Supplementary Table 8<\/a>), suggesting that the default MinKNOW configuration is discarding shorter tRNA molecules and preferentially capturing longer ones, such as tRNA molecules with variable arms (for example, tRNALeu<\/sup>, tRNAArg<\/sup> and tRNASer<\/sup>). Moreover, the relative proportion of tRNA reads was better recapitulated using the custom configuration than default settings (Fig. 3f<\/a>), and the reported tRNA abundances using custom settings correlated well to the expected values (\u03c1\u2009=\u20090.93) (Fig. 3g<\/a> and Supplementary Tables 9<\/a> and 10<\/a>).<\/p>\nReverse transcription of tRNAs increases sequencing yield<\/h3>\n
We next questioned whether removing the tRNA structure, which can be achieved via linearization by reverse transcription, would further improve our sequencing yield. We should highlight that, in the case of DRS, the native RNA molecule is sequenced, whereas the cDNA strand is not (see Fig. 4a<\/a> schematic). A linear tRNA molecule may (1) reduce the clogging of pores, allowing more reads to be sequenced, and maintain the integrity of the flowcell longer and\/or (2) stabilize the tRNA translocation speed through the pore, improving the accuracy of basecalling algorithms. Notably, tRNAs are notoriously difficult to fully and accurately reverse transcribe due to their compact secondary and tertiary structures as well as their abundance of modifications that disrupt the Watson\u2013Crick base pairing74<\/a>,80<\/a>,86<\/a><\/sup> (Fig. 1a<\/a>). To examine whether tRNA linearization might improve sequencing yield, we tested a range of commercial reverse transcriptases and incubation conditions on both IVT and native tRNAs and examined their cDNA outputs using TapeStation (Extended Data Fig. 5a,b<\/a> and Methods<\/a>). We found that both Maxima and SuperScript IV at 60\u2009\u00b0C offered the best performance in the production of full-length cDNA products, and we opted to use Maxima at 60\u2009\u00b0C in our subsequent tRNA sequencing experiments (Extended Data Fig. 5b<\/a>).<\/p>\n\nFig. 4: Nano-tRNAseq can quantify tRNA abundance and RNA modifications as well as capture modification interdependencies.<\/b><\/figcaption>
Transfer RNAs (tRNAs) are abundant small non-coding RNAs that play a pivotal role in decoding genetic information1<\/a>,2<\/a>,3<\/a><\/sup>. Based on their aminoacylation identity, tRNAs are subdivided into 20 accepting groups (alloacceptors), each comprising several tRNAs that translate synonymous codons with the same amino acid (isoacceptors). To fulfill their function as adapter molecules between the RNA and protein codes, tRNAs are extensively modified, containing on average 13 modifications per tRNA molecule4<\/a><\/sup>. Although some tRNA modifications are thought to be structural and static, others are dynamic and even reversible5<\/a>,6<\/a>,7<\/a>,8<\/a><\/sup>, affecting the fate and function of individual tRNA molecules2<\/a>,9<\/a>,10<\/a>,11<\/a>,12<\/a>,13<\/a><\/sup>. Notably, mutations in multiple tRNA modification enzymes have been associated with a wide variety of human diseases14<\/a>,15<\/a>,16<\/a>,17<\/a><\/sup>, highlighting their importance in proper cellular functioning.<\/p>\n tRNA modifications are present in all domains of life18<\/a>,19<\/a><\/sup> and are synthesized by dedicated tRNA-modifying enzymes that alter specific tRNAs in a site-specific fashion20<\/a><\/sup>. The chemical nature of these modifications is highly diverse and includes methylations, acetylations, isomerizations, deaminations and conjugation to amino acids, among others21<\/a>,22<\/a><\/sup>. Certain tRNA modifications are found only in a single tRNA species, whereas others are found in multiple tRNA species20<\/a>,23<\/a><\/sup>. For example, 2-lysidine (k2<\/sup>C) tRNA modifications are placed at position 34 of the anticodon of tRNAIle<\/sup>(AUA)24<\/a><\/sup>, whereas pseudouridine (\u03a8) can be placed at diverse positions of the tRNA molecule and in multiple tRNA isoacceptors25<\/a>,26<\/a>,27<\/a><\/sup>. Regarding their function, tRNA modifications can sometimes act as identity elements recognized by aminoacyl-tRNA synthetases28<\/a>,29<\/a>,30<\/a><\/sup>, and, without modifications, many tRNAs have poor aminoacylation capability31<\/a><\/sup>. On the other hand, tRNA modifications can affect the decoding preferences of tRNA molecules, especially those found at position 34 of the anticodon16<\/a>,32<\/a>,33<\/a>,34<\/a><\/sup>, restricting or increasing the wobbling capacity of the tRNAs and, consequently, changing the set of \u2018preferred\u2019 or \u2018optimal\u2019 codons that will be translated35<\/a>,36<\/a>,37<\/a>,38<\/a>,39<\/a>,40<\/a>,41<\/a>,42<\/a>,43<\/a>,44<\/a><\/sup>.<\/p>\n In the last few years, it has been shown that some tRNA modifications are reversible8<\/a>,45<\/a>,46<\/a>,47<\/a><\/sup> and can be dynamically regulated upon environmental exposures48<\/a>,49<\/a>,50<\/a><\/sup>, across cell cycle stages51<\/a><\/sup> and upon tumorigenesis52<\/a>,53<\/a>,54<\/a>,55<\/a>,56<\/a><\/sup>. Similarly, tRNA abundances are also dysregulated upon environmental exposures, such as oxidative stress48<\/a>,57<\/a>,58<\/a><\/sup>, as well as in diverse types of cancer52<\/a>,59<\/a>,60<\/a>,61<\/a><\/sup>. Modulation of tRNA abundances and\/or tRNA modifications is generally regarded as a molecular strategy that allows cells to adapt to distinct physiological states or conditions, leading to increased expression of subsets of proteins that otherwise would remain poorly translated under \u2018normal\u2019 tRNA abundances60<\/a>,62<\/a>,63<\/a><\/sup>.<\/p>\n Despite the pivotal function that tRNAs play in cellular processes and their involvement in numerous human diseases, we currently lack a simple and cost-effective method to accurately quantify both tRNA abundances and their modifications systematically. On the one hand, tRNA modifications are typically identified and quantified with high accuracy using liquid chromatography coupled to mass spectrometry (LC\u2013MS) methodologies48<\/a>,64<\/a>,65<\/a>,66<\/a>,67<\/a>,68<\/a>,69<\/a><\/sup>. In these methods, RNA molecules are fragmented into oligomers or monomers, and their abundance is measured via UV absorption or MS\/MS. LC\u2013MS\/MS techniques using triple quadrupole-based detection are among the most sensitive, allowing limits of quantification in the low femtomole range48<\/a>,64<\/a>,70<\/a>,71<\/a>,72<\/a><\/sup>, but they typically cannot identify the tRNA isoacceptor that contained each detected modification. On the other hand, tRNA abundances can be determined using tRNA microarrays36<\/a>,59<\/a>,73<\/a><\/sup> and, more recently, by employing next-generation sequencing (NGS) technologies74<\/a>,75<\/a>,76<\/a>,77<\/a>,78<\/a>,79<\/a>,80<\/a>,81<\/a>,82<\/a><\/sup>, which require an initial conversion of the tRNA molecules into cDNA. Consequently, NGS-based methods are blind to most tRNA modifications, as these are typically erased during the reverse transcription step. Moreover, tRNA modifications that disrupt the Watson\u2013Crick base pairing, which are abundant in tRNAs, will interfere with the reverse transcriptase enzyme, causing it to drop off, producing truncated reads, in addition to misincorporations72<\/a>,83<\/a>,84<\/a>,85<\/a><\/sup> (Fig. 1a<\/a>). To overcome these limitations, a wide variety of improved tRNA sequencing protocols have been developed in recent years, which often include the use of highly processive reverse transcriptase enzymes75<\/a>,77<\/a><\/sup> and\/or cocktails of demethylases74<\/a>,75<\/a>,79<\/a><\/sup>. However, despite these improvements, NGS-based methods still suffer from the following caveats: (1) they introduce significant biases during the library preparation, caused by incomplete reverse transcriptions84<\/a>,86<\/a><\/sup>, incomplete demethylations87<\/a><\/sup> and polymerase chain reaction (PCR) amplification88<\/a><\/sup>, resulting in skewed representations of existing tRNA populations; and (2) they cannot detect most tRNA modifications, as these are typically \u2018erased\u2019 during the conversion of RNA to cDNA. Therefore, a simple, robust and efficient tRNA sequencing method is still needed.<\/p>\n a<\/b>, Schematic of the modifications found in S. cerevisia<\/i>e cytoplasmic tRNA, shown in its usual secondary structure form with circles representing nucleotides and lines representing base pairs. Gray circles represent unmodified nucleotides; pink circles represent possible modification sites; and those with a black outline indicate modifications that cause errors during reverse transcription. Possible RNA modifications occurring at each position are listed in the surrounding boxes; modifications that cause misincorporation during reverse transcription are in green; and those that cause reverse transcription truncation are in blue. b<\/b>, Schematic overview illustrating the steps required for tRNA library preparation using Nano-tRNAseq (see Extended Data Fig. 1<\/a> for more details). c<\/b>, IGV snapshots of Nano-tRNAseq mapped reads from synthetic IVT tRNAs (upper panels) or biological tRNAs (lower panel). Positions with a mismatch frequency greater than 0.2 are colored, whereas those showing mismatch frequencies lower than 0.2 are shown in gray. d<\/b>, Scatter plot of tRNA abundances showing the replicability of Nano-tRNAseq when WT S. cerevisia<\/i>e tRNA biological replicates are sequenced. The correlation strength is indicated by Spearman\u2019s correlation coefficient (\u03c1). RT, reverse transcription.<\/p>\n<\/div>\n Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n A promising alternative to the use of NGS-based technologies to characterize the tRNAome is the direct RNA sequencing (DRS) platform developed by Oxford Nanopore Technologies (ONT). This technology allows direct sequencing of native RNA molecules and, as such, can, in principle, detect and measure both tRNA modifications and tRNA abundances without the need for reverse transcription or PCR. Previous works have demonstrated that nanopores can capture tRNAs using solid-state or biological (ONT) nanopores89<\/a>,90<\/a>,91<\/a>,92<\/a><\/sup>. For example, sequencing of five distinct tRNAs was achieved using solid-state nanopores89<\/a><\/sup>, and tRNAs were shown to be distinguishable from other short RNAs using the MspA pore90<\/a><\/sup>. Later studies showed that, by lengthening the tRNA molecules with ligated adapter extensions, tRNAs could be sequenced, basecalled and mapped using biological nanopores92<\/a><\/sup>. However, in these studies, the proposed approach led to low sequencing yields of tRNA molecules (~20\u201340\u00d7 lower than expected for a DRS run) and did not report whether extant in vivo tRNA abundances and\/or tRNA modifications were recapitulated using this method.<\/p>\n Here we present Nano-tRNAseq, a nanopore-based approach that allows accurate and direct measurement of native tRNA molecules. The library preparation benefits from the 3\u2032 CCA overhang typically present in the mature tRNAs to incorporate a double ligation of RNA adapters at both the 5\u2032 and 3\u2032 ends of the tRNA molecules, which leads to an improved proportion of basecalled and mapped tRNA molecules. Moreover, we show that MinKNOW, the ONT proprietary software required to run nanopore sequencing experiments, erroneously discards the majority of tRNA reads, misinterpreting them as \u2018adapters\u2019, and also causes biases in the estimated tRNA abundances due to preferential capture of longer tRNAs (for example, tRNALeu<\/sup>, tRNAArg<\/sup> and tRNASer<\/sup>). To overcome these limitations, here we provide a computational framework that allows us to capture ~10\u00d7 more tRNA reads and accurately recapitulates tRNA abundances.<\/p>\n Altogether, our work provides a simple, cost-effective, high-throughput and reproducible method to accurately quantify tRNA abundances and capture tRNA modification changes simultaneously using native tRNA nanopore sequencing, providing a framework to study the tRNAome at single-molecule resolution. We envision that Nano-tRNAseq will contribute to the study of the biological function of tRNA modifications in a wide variety of contexts, such as cancer, stress exposures or viral infection, and opens the possibility of exploiting these molecules as biomarkers of human health and disease.<\/p>\n<\/div>\n<\/div>\n Nanopore DRS is a well-established long-read sequencing technology to study RNA molecules, typically polyadenylated mRNAs93<\/a>,94<\/a>,95<\/a>,96<\/a>,97<\/a><\/sup>. Although several works have shown that this technology can also be used to study short RNA molecules, such as snoRNAs and snRNAs96<\/a>,98<\/a><\/sup>, DRS is inefficient at capturing RNA molecules shorter than 200\u2009nucleotides (nt) and is generally considered unable to capture sequences shorter than ~100\u2009nt96<\/a>,97<\/a><\/sup>, limiting its applicability to study short RNA populations, such as tRNAs. In addition, the first ~15\u2009nt at the 5\u2032 end of RNA molecules are typically lost in DRS runs93<\/a>,99<\/a><\/sup>, as this portion cannot be adequately basecalled due to the increase in the RNA translocation speed when the 5\u2032 end of the molecule exits the helicase99<\/a><\/sup>. To overcome these limitations, we reasoned that extension of the 5\u2032 and 3\u2032 ends of the tRNA would lead to improved sequencing of tRNA molecules, as these would now be beyond the ~100-nt threshold, in addition to capturing the sequence and modification information of 5\u2032 tRNA ends.<\/p>\n We first attempted a modified tRNA DRS approach in which a 5\u2032 RNA adapter, complementary to the 3\u2032 CCA overhang present in mature tRNA molecules, was ligated to the 5\u2032 end of tRNAs that had been previously in vitro polyadenylated (Strategy A; Extended Data Fig. 1<\/a> and Methods<\/a>). A set of nine synthetic in vitro transcribed (IVT) tRNAs of various lengths and sequences (Supplementary Table 1<\/a>) were sequenced using this strategy. However, this approach produced poor sequencing yields (56,002 reads; Supplementary Table 2<\/a>) compared to a standard DRS run (~1\u20132 million reads). Moreover, only 7.5% of reads mapped uniquely to tRNAs using minimap2 (ref. 100<\/a><\/sup>) with recommended parameters (-ax map -ont -k15) (Supplementary Table 2<\/a>). Relaxation of the mapping parameters (-ax map-ont -k5), which had previously been shown to improve the mappability of highly modified RNAs98<\/a><\/sup>, did not significantly increase the number of mapped tRNA reads (Supplementary Table 2<\/a>).<\/p>\n Next, we altered our library preparation protocol to replace the poly(A) tail with a 3\u2032 DNA adapter complementary to the 5\u2032 RNA adapter (Strategy B), such that the two oligonucleotides could be pre-annealed and ligated to the tRNA (Extended Data Fig. 1<\/a>). However, this strategy also yielded a low number of sequenced reads (63,502 reads) and a low percentage of uniquely mapped reads (6.5%) (Supplementary Table 2<\/a>). We speculated that the low number of reads could be due to steric interference of the poly(A) preventing 5\u2032 ligation (Strategy A) or that the 3\u2032 DNA adapter is not basecalled (Strategy B). These scenarios would lead to low coverage of the 5\u2032 or 3\u2032 ends, respectively, decreasing the mappability of reads resulting from these two strategies.<\/p>\n Based on the results of Strategy A and Strategy B, we rationalized that padding the 5\u2032 and 3\u2032 tRNA ends with RNA adapters, which can be accurately basecalled and mapped, would enable us to capture the entirety of the tRNA sequence. This approach, which we termed Nano-tRNAseq (Fig. 1b<\/a>), was the most successful at sequencing, basecalling and mapping both in vitro and native tRNA molecules using nanopore DRS. We should note that a recent work also proposed a similar solution to facilitate native tRNA nanopore sequencing92<\/a><\/sup>.<\/p>\n In the first step, a 5\u2032 RNA adapter (orange) complementary to the CCA overhang of mature tRNAs is pre-annealed to a 3\u2032 RNA adapter (red) containing three DNA bases (pink) at the 3\u2032 end (Fig. 1b<\/a>). In preliminary Nano-tRNAseq runs, we used an RNA-only 3\u2032 adapter but observed that RNA-only 3\u2032 adapters led to increased self-ligation (Supplementary Table 2<\/a>), an issue that we mitigated by adding DNA bases to the end of the adapter (Extended Data Fig. 1<\/a>). Next, the pre-annealed 5\u2032 RNA and 3\u2032 RNA:DNA splint adapters were ligated to deacylated tRNAs. Knowing that an RNA:RNA ligation with an RNA bridge has a low efficiency101<\/a><\/sup>, various ligation times and the addition of a molecular crowding agent were tested to ensure that conditions that maximized ligation efficiency were chosen (Extended Data Fig. 2a,b<\/a>). Subsequently, ONT RTA oligoA and oligoB were pre-annealed and ligated to the tRNA molecule using T4 DNA Ligase (see Extended Data Fig. 3<\/a> for validation of each ligation step). This approach resulted in >200,000 basecalled reads (Supplementary Table 2<\/a>), thus significantly increasing sequencing output by up to fourfold relative to the previous strategies, and also with improved 5\u2032 and 3\u2032 coverage of both synthetic and biological tRNAs (Fig. 1c<\/a>). Additionally, we found Nano-tRNAseq to be highly replicable when sequencing native tRNAs (\u03c1\u2009=\u20090.984) (Fig. 1d<\/a>).<\/p>\n The alignment of native tRNA reads is challenging due to their short and highly modified nature. Indeed, native tRNAs contain a large proportion of mismatched bases, often originating from inaccurate basecalling of modified bases in DRS datasets94<\/a>,98<\/a>,102<\/a><\/sup>. As a consequence of these miscalled bases, the commonly used long-read mapper minimap2 with recommended settings (-ax map-ont -k15) aligned only a fraction (2.56%) of the reads (Fig. 2a\u2013c<\/a> and Supplementary Table 2<\/a>). We further tested a range of minimap2 parameters and observed only incremental improvements in mapping and an increase in false alignments (antisense mapped reads served as a proxy of mismapping) (Supplementary Table 3<\/a>).<\/p>\n a<\/b>,b<\/b>, IGV snapshots of reads mapped to IVT D. melanogaster<\/i> mitochondrial tRNAAla(UGC)<\/sup> (a<\/b>) and S. cerevisiae<\/i> tRNAPhe<\/sup> (b<\/b>), mapped using different mapping algorithms (minimap2, BWA-MEM and BWA-SW) and parameters. The 5\u2032 and 3\u2032 RNA adapter regions, ligated to the ends of the tRNA molecule, were included in the mapping references and are represented by an orange bar and a red bar, respectively. Positions with a mismatch frequency greater than 0.2 are colored, whereas those showing mismatch frequencies lower than 0.2 are shown in gray. c<\/b>, Bar plot depicting the effect of algorithm and parameter choice on the relative proportion of uniquely mapped reads (green) and mismapped reads (purple; reads mapping to antisense strands were used as a proxy to assess mismapping) (Supplementary Table 4<\/a>). The proportion of mapped reads (d<\/b>) and alignment identity (e<\/b>) for each template from the bar plot in c<\/b>, using either minimap2 or bwa mem -W13 -k6- T20. We should note that minimap2 alignment identity in S. cerevisiae<\/i> tRNAPhe<\/sup> was not computed because no reads were mapped to this tRNA using minimap2 with -ax map-ont -k15 parameters (Supplementary Table 5<\/a>). f<\/b>, Bar plot showing the effect of trimming the length of the 5\u2032 RNA adapter (reds) and 3\u2032 RNA adapter (blues) on tRNA read mappability (Supplementary Table 6<\/a>). The conditions used by Nano-tRNAseq are gray, whereas the effect of not using RNA adapters is black.<\/p>\n<\/div>\n Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n To improve the mappability of Nano-tRNAseq reads, we next tested the mapping algorithm BWA, a short-read mapping algorithm commonly used to map Illumina reads103<\/a>,104<\/a><\/sup>. Using sequencing data from a Nano-tRNAseq run that contained three different tRNA constructs (IVT Drosophila melanogaster<\/i> mitochondrial tRNAAla(UGC)<\/sup>, IVT Streptococcus pneumoniae<\/i> tRNASer(UGA)<\/sup> and native S. cerevisiae<\/i> tRNAPhe<\/sup>; see Supplementary Table 2<\/a> for a summary of the sequencing runs in this work), we found that the BWA-MEM aligner with recommended parameters outperformed minimap2 in terms of proportion of mapped reads (Fig. 2c<\/a>), in agreement with recent works92<\/a><\/sup>. Although more relaxed configurations of BWA-MEM aligned more reads, this also came at the expense of increased false alignments (Fig. 2c<\/a> and Supplementary Table 4<\/a>). An optimal balance between increased mapped reads and false alignments was found when using bwa mem with parameters -W13 -k6 -xont2d -T20, which mapped 54.63% of the reads, with very few false alignments (0.19%). When comparing the performance of the mapping algorithms in native tRNA molecules, the contrast was even more stark; although minimap2 mapped IVT tRNAs, it failed to map a single biological tRNA read (Fig. 2d<\/a> and Supplementary Table 5<\/a>). The alignment identity was similar to minimap2 for reads that mapped to IVT tRNAs but was slightly lower than the typical identity obtained in nanopore DRS runs, suggesting that short reads, even without modifications, cause a drop in the basecalling accuracy (Fig. 2e<\/a>). Notably, the alignment identity of S. cerevisiae<\/i> tRNAPhe<\/sup> was lower (~74.5%) than in synthetic tRNAs (81.8%), presumably due to the presence of base modifications present on endogenously modified S. cerevisiae<\/i> tRNAPhe<\/sup> (Fig. 2e<\/a> and Supplementary Table 5<\/a>).<\/p>\n We then assessed whether the mappability of Nano-tRNAseq reads might be affected by the length of the 5\u2032 and 3\u2032 RNA adapters. To simulate different RNA adapter lengths, we trimmed one or both adapters from the reference sequences. We found that the absence of both RNA adapters had only a modest effect on the mappability of reads originating from IVT tRNAs (decrease of 6\u201311%) (Fig. 2f<\/a> and Supplementary Table 6<\/a>), whereas their absence had a major effect in the mappability of native S. cerevisiae<\/i> tRNAPhe<\/sup> (55% loss). Thus, we concluded that short and unmodified sequences can be aligned efficiently even without RNA adapters in the reference sequences, whereas short and modified reads greatly benefit from the extension with adapters, demonstrating that extending molecules with RNA adapters is essential for guiding the correct alignment of short reads enriched in \u2018mismatches\u2019, such as those derived from native tRNAs. For native S. cerevisiae<\/i> tRNAPhe<\/sup>, the read mappability plateaued at a 3\u2032 adapter length of 25\u2009nt (Fig. 2f<\/a>), suggesting that the 30-nt RNA portion of the 3\u2032 RNA:DNA oligo used in Nano-tRNAseq is more than sufficient to achieve optimal read mappability.<\/p>\n A surprising feature of our initial tRNA sequencing runs was the low amount of sequenced reads. Although pore clogging caused by tRNA structure might partially explain the low sequencing yield92<\/a><\/sup>, we also noticed that the MinKNOW software classified a high proportion of reads as \u2018adapter-only\u2019 reads in real time. Hence, we hypothesized that a considerable fraction of tRNA reads might be discarded by the MinKNOW software due to their short signal lengths, as they resemble \u2018adapter-only\u2019 reads.<\/p>\n The MinKNOW software is responsible for analyzing the continuous electrical current (signal intensity) measured at each pore, reporting the signal regions that correspond to \u2018reads\u2019 into FAST5 files, which are then basecalled to generate a FASTQ file (Fig. 3a<\/a>). We noted that MinKNOW, by default, reports reads that last at least 2\u2009seconds, roughly corresponding to RNA molecules of 140\u2009nt (assuming constant helicase processivity of 70\u2009nt per second in DRS) (Extended Data Fig. 4<\/a>). Considering that canonical tRNA molecules are ~73\u2009nt, this would imply that even after double RNA adapter ligation (where 24 and 30 RNA nucleotides are added to the 5\u2032 and 3\u2032 ends of the tRNA molecule, respectively), the size of the ligated tRNA molecule would still be below the threshold, possibly leading to misassignment of tRNA reads as \u2018adapter-only\u2019 reads. To alleviate this issue, we tested whether alternative MinKNOW configurations would improve the classification of tRNA reads and boost sequencing yields. To this end, the bulk \u2018raw\u2019 dump files were saved during the sequencing and were reprocessed using alternative MinKNOW configurations (Supplementary Table 7<\/a>).<\/p>\n a<\/b>, MinKNOW software classifies continuous current passing through pores as open pore, adapter or strand (actual reads) and outputs fragments classified as strand to a FAST5 file, which are then basecalled to generate a FASTQ file. b<\/b>, Diagram showing the conceptual difference between default and custom MinKNOW read classification (Extended Data Fig. 4<\/a>). c<\/b>, Bar plot of sequencing yield in terms of basecalled and uniquely mapped reads obtained with default and custom configurations (Supplementary Table 7<\/a>). d<\/b>, Scatter plot of the relative fold change of uniquely mapped reads with respect to tRNA length (Supplementary Table 8<\/a>). e<\/b>, Histogram of read count and alignment length of IVT tRNA reads captured with default and custom configurations. f<\/b>, Bar plot of the relative proportion of IVT tRNA molecules D. melanogaste<\/i>r mitochondrial tRNAAla(UGC)<\/sup> and S. pneumoniae<\/i> tRNASer(UGA)<\/sup> and native S. cerevisiae<\/i> tRNAPhe<\/sup> reads recovered with default and custom settings (Supplementary Table 9<\/a>), where the dotted line indicates the expected proportion. g<\/b>, Expected versus observed log read counts of nine IVT and one native tRNA molecules captured using the custom MinKNOW configuration (Supplementary Table 10<\/a>). Spearman correlation (\u03c1) is shown.<\/p>\n<\/div>\n Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n By lowering the MinKNOW strand minimum duration to 1\u2009second and the adapter maximum duration to 2\u2009seconds (see Extended Data Fig. 4<\/a> for schematic), a configuration that we refer to as custom<\/i>, we captured ~12-fold more basecalled and ~4.5-fold more uniquely mapped tRNA reads compared to the default MinKNOW configuration (Fig. 3b,c<\/a> and Supplementary Table 7<\/a>). Notably, we found that the default MinKNOW configuration not only led to low sequencing outputs but also caused significant biases in the relative abundances of tRNA molecules. Specifically, we found a greater representation of shorter tRNAs in our custom configuration (Fig. 3d,e<\/a> and Supplementary Table 8<\/a>), suggesting that the default MinKNOW configuration is discarding shorter tRNA molecules and preferentially capturing longer ones, such as tRNA molecules with variable arms (for example, tRNALeu<\/sup>, tRNAArg<\/sup> and tRNASer<\/sup>). Moreover, the relative proportion of tRNA reads was better recapitulated using the custom configuration than default settings (Fig. 3f<\/a>), and the reported tRNA abundances using custom settings correlated well to the expected values (\u03c1\u2009=\u20090.93) (Fig. 3g<\/a> and Supplementary Tables 9<\/a> and 10<\/a>).<\/p>\n We next questioned whether removing the tRNA structure, which can be achieved via linearization by reverse transcription, would further improve our sequencing yield. We should highlight that, in the case of DRS, the native RNA molecule is sequenced, whereas the cDNA strand is not (see Fig. 4a<\/a> schematic). A linear tRNA molecule may (1) reduce the clogging of pores, allowing more reads to be sequenced, and maintain the integrity of the flowcell longer and\/or (2) stabilize the tRNA translocation speed through the pore, improving the accuracy of basecalling algorithms. Notably, tRNAs are notoriously difficult to fully and accurately reverse transcribe due to their compact secondary and tertiary structures as well as their abundance of modifications that disrupt the Watson\u2013Crick base pairing74<\/a>,80<\/a>,86<\/a><\/sup> (Fig. 1a<\/a>). To examine whether tRNA linearization might improve sequencing yield, we tested a range of commercial reverse transcriptases and incubation conditions on both IVT and native tRNAs and examined their cDNA outputs using TapeStation (Extended Data Fig. 5a,b<\/a> and Methods<\/a>). We found that both Maxima and SuperScript IV at 60\u2009\u00b0C offered the best performance in the production of full-length cDNA products, and we opted to use Maxima at 60\u2009\u00b0C in our subsequent tRNA sequencing experiments (Extended Data Fig. 5b<\/a>).<\/p>\n<\/picture><\/a><\/div>\nResults<\/h2>\n
Standard nanopore DRS results in low tRNA sequencing yields<\/h3>\n
Extending tRNA ends with RNA adapters improves basecalling<\/h3>\n
Mapping parameters significantly affect tRNA read mappability<\/h3>\n
<\/picture><\/a><\/div>\nCustom MinKNOW improves yield and tRNA abundance estimates<\/h3>\n
<\/picture><\/a><\/div>\nReverse transcription of tRNAs increases sequencing yield<\/h3>\n
![\"Science](\"http:\/\/media.springernature.com\/lw685\/springer-static\/image\/art%3A10.1038%2Fs41587-023-01743-6\/MediaObjects\/41587_2023_1743_Fig4_HTML.png\")
<\/picture><\/a><\/div>\n