METHOD FOR CDNA LIBRARY CONSTRUCTION AND ANALYSIS FROM TRANSFER RNA
20260117220 · 2026-04-30
Inventors
Cpc classification
C12N15/1068
CHEMISTRY; METALLURGY
International classification
Abstract
The present invention relates to a method for the generation of a cDNA library from transfer RNA (tRNA) comprising (a) optionally ligating at least one DNA adapter to 3-end of tRNA, wherein 3-end of the DNA adapter is preferably a chain terminator dideoxycytidine, preferably under one or more of the following conditions: (i) crowding reagent at 5% to 35%, preferably 15% to 30%, and most preferably about 25%, (ii) MgCl.sub.2 concentration of 1 mM to 15 mM, preferably 3 to 12 mM and most preferably about 10 mM (iii) a temperature of 12 C. to 37 C., preferably of about 25 C. (iv) a pH of 6.0 to 9.0, preferably 6.5 to 8.0 and most preferably about 7.0, (v) reducing agent concentration of 0.1 mM to 10 mM, preferably 0.5 to 5 mM and most preferably about 1 mM, and (vi) a reaction time of at least 30 min, preferably at least 1.5 h and most preferably at least 3 h; and (b) reverse transcription in a primer-dependent reaction, in case step (a) is present, or in a template-switching reaction, in case step (a) is absent, of tRNAs into cDNA by a group II intron reverse transcriptase, under the following conditions: (i) KCl or NaCl at a concentration of 20 mM to 250 mM, preferably 50 to 100 mM and most preferably about 75 mM, (ii) MgCl.sub.2 at a concentration of 0.5 mM to 15 mM, preferably 1 to 5 mM and most preferably about 3 mM, (iii) a temperature of 30 C. to 65 C., preferably of about 42 C., and (iv) a reaction time of at least 2 h, preferably at least 8 h and most preferably at least 15 h, and preferably (v) a pH of 6.5 to 9.5, preferably 7.0 to 8.5 and most preferably about 8.0, and/or (vi) a reducing agent (DTT) at a concentration of 1 mM to 12.5 mM, preferably 3 to 8 mM and most preferably about 5 mM.
Claims
1. A method for the generation of a cDNA library from tRNA, preferably a pool of tRNAs, comprising (a) optionally ligating at least one DNA adapter to the 3-end of tRNA, preferably the tRNAs in the pool, wherein the 3-end of the DNA adapter is preferably a chain terminator dideoxycytidine, preferably under one or more of the following conditions: (i) crowding reagent at 5% to 35%, preferably 15% to 30%, and most preferably about 25%, (ii) MgCl.sub.2 concentration of 1 mM to 15 mM, preferably 3 to 12 mM and most preferably about 10 mM, (iii) a temperature of 12 C. to 37 C., preferably of about 25 C., (iv) a pH of 6.0 to 9.0, preferably 6.5 to 8.0 and most preferably about 7.0, (v) reducing agent concentration of 0.1 mM to 10 mM, preferably 0.5 to 5 mM and most preferably about 1 mM, (vi) a reaction time of at least 30 min, preferably at least 1.5 h and most preferably at least 3 h; and (b) reverse transcription in a primer-dependent reaction, in case step (a) is present, or in a template-switching reaction, in case step (a) is absent, of the tRNA, preferably of the tRNAs in the pool into cDNA(s) by a group II intron reverse transcriptase, under the following conditions: (i) KCl or NaCl at a concentration of 20 mM to 250 mM, preferably 50 to 100 mM and most preferably about 75 mM, (ii) MgCl.sub.2 at a concentration of 0.5 mM to 15 mM, preferably 1 to 5 mM and most preferably about 3 mM, (iii) a temperature of 30 C. to 65 C., preferably of about 42 C., and (iv) a reaction time of at least 2 h, preferably at least 8 h and most preferably at least 15 h, and preferably (v) a pH of 6.5 to 9.5, preferably 7.0 to 8.5 and most preferably about 8.0, and/or (vi) a reducing agent (DTT) at a concentration of 1 mM to 12.5 mM, preferably 3 to 8 mM and most preferably about 5 mM.
2. The method of claim 1 further comprising prior to step (a), (a) the purification of tRNAs from an RNA preparation.
3. The method of claim 2 further comprising prior to step (a), (a) the isolation of an RNA preparation from eukaryotic cells.
4. The method of claim 1, further comprising (c) the construction of a sequencing library by PCR amplification of the cDNA(s) as obtained in step (b), wherein the cDNA(s) is/are optionally circularized prior to the amplification.
5. The method of claim 1, wherein the DNA adapter comprises at least two, three or four DNA adapters that are distinguished from each other by a barcode sequence.
6. The method of claim 4, further comprising (d) sequencing the sequencing library as obtained in step (c).
7. The method of claim 6, further comprising (e) aligning the sequencing reads as obtained in step (d) to known tRNA reference sequences.
8. The method of claim 6, further comprising (f) aligning the sequencing reads as obtained in step (d) to a reference of known tRNA transcript sequences, wherein the reference comprises information on the identity and location of known modified ribonucleosides in tRNAs.
9. The method of claim 8, wherein the reads are aligned using a short read alignment algorithm, preferably by using the Genomic Short-read Nucleotide Alignment Program to the reference generated in (f), wherein the reference is preferably clustered into clusters of tRNA genes sharing an anticodon.
10. The method of claim 8, wherein reads are aligned to RNA reference sequences in single nucleotide polymorphism (SNP)-tolerant mode, using known modified ribonucleotides as potential sites of mismatch to the reference sequence.
11. The method of claim 8, further comprising (g) subjecting the reads aligned to clusters to a deconvolution algorithm that assigns the aligned reads to unique tRNA species.
12. The method of claim 8, further comprising after read deconvolution (h) the analysis of one or more of read coverage, 3-CCA, differential tRNA abundance and modification profiling.
13. The method of claim 1, wherein the group II intron reverse transcriptase is a thermostable group II intron reverse transcriptase.
14. The method of claim 13, the thermostable group II intron reverse transcriptase is a TGIRT RT or the MarathonRT.
15. The method of claim 14, wherein TGIRT has the sequence of any one of SEQ ID NOs: 1 to 5 or a sequence being at least 80% identical thereto and/or wherein the MarathonRT has the sequence of SEQ ID NO: 6 or a sequence being at least 80% identical thereto.
Description
[0100] The figures show.
[0101]
[0104]
[0110]
[0115]
[0122]
[0128]
[0134] The examples illustrate the invention.
Example 1Design
Efficient Sequencing Library Generation from Native Eukaryotic tRNA Pools
[0135] To develop a method for high-resolution tRNA quantitation, it was focused on improving the efficiency of full-length cDNA synthesis from endogenously modified tRNAs by TGIRT. This enzyme can attach adapter sequences to RNA by template switching (Mohr et al., 2013), which circumvents potential hindrances to 3 adapter ligation and RT posed by tRNA structure (Katibah et al., 2014; Qin et al., 2016; Zheng et al., 2015). TGIRT can also read through a subset of Watson-Crick face modifications more efficiently than other commercial RTs (Li et al., 2017), albeit with reduced fidelity (Katibah et al., 2014; Qin et al., 2016; Zheng et al., 2015). Despite these advantages, RT stops at modified sites in tRNA are still pervasive in TGIRT-mediated reactions (Clark et al., 2016; Zheng et al., 2015), and cDNA yield is extremely low (Zhao et al., 2018; Zheng et al., 2015).
[0136] As TGIRT is active in a wide range of conditions (Mohr et al., 2013), it was asked whether its efficiency on tRNA templates can be further improved. To test this, tRNA pools were first purified from S. cerevisiae and human K562 cells by gel size selection of 60-100 nt RNAs from total RNA. These were then used, along with a synthetic unmodified E. coli tRNA-Lys-UUU, in template-switching TGIRT reactions. The cDNA yield from all templates was minimal under conditions previously used for tRNA sequencing (450 mM salt, 60 C.; Katibah et al., 2014; Qin et al., 2016; Zheng et al., 2015) but dramatically improved at lower temperatures and salt concentration (
a Comprehensive Computational Framework for tRNA Sequencing Data Analysis
[0137] We reasoned that the increase in full-length cDNA reads would reduce alignment ambiguity. However, given TGIRT's low fidelity at modified sites, it was expected many tRNA-derived reads to contain multiple mismatches to the reference genome. Another source of mismatches are non-templated nucleotides added to 3 cDNA ends by TGIRT and other RTs (Chen and Patton, 2001; Mohr et al., 2013). Such read extensions are penalized by most algorithms but can be recognized and dynamically processed (soft-clipped) by some aligners. It was therefore asked how two short-read aligners commonly used for tRNA analysis-Bowtie (Langmead et al., 2009) and Bowtie 2 (Langmead and Salzberg, 2012)would perform on a tRNA sequencing dataset from human HEK293T cells obtained with our improved library construction protocol (
[0138] We first generated a non-redundant reference of 420 mature tRNA transcripts from 596 curated nuclear- and mitochondrial-encoded tRNA genes retrieved from GtRNAdb and mitotRNAdb (Chan and Lowe, 2016; Juhling et al., 2009;
[0139] In contrast, Bowtie 2's lack of mismatch restrictions and ability to soft-clip read ends make it seem more suited for tRNA read mapping. High mismatch tolerance, however, compounds the problem of misalignment: while Bowtie 2 increased alignment rates of our HEK293T-derived dataset to 82%, most mapped reads (85%) could not be assigned to a single reference (
[0140] Given these limitations, it was reasoned that an accurate tRNA read analysis workflow requires solutions to two main challenges: alignment bias against reads with modification-induced misincorporations and multi-mapping of reads from nearly identical tRNAs. To tackle the first issue, it was taken advantage of the comprehensive annotation of tRNA modifications in MODOMICS (Boccaletto et al., 2018) and used these data to enable position-specific mismatch tolerance during alignment (
[0141] On the basis of these premises, a new computational workflow was developed to suit the intricacies of tRNA sequencing data (
[0142] This computational workflow dramatically improved both the efficiency and accuracy of tRNA read alignment. Both clustering and SNP tolerance at modified sites prevented data loss for defined tRNA subsets. A cluster ID of 0.95 maximized unique transcript resolution and minimized multi-mapping for human tRNAs, yielding 86% uniquely mapped and only 2.5% ambiguously aligned reads (
Example 2Results
The mim-tRNAseq Workflow Alleviates tRNA Sequencing Bias
[0143] To benchmark our workflow, mim-tRNAseq was used to analyze HEK293T tRNAs and compared our results with those published for the same cell type with DM-tRNAseq (Zheng et al., 2015) and from the closely related HEK293 T-Rex Flp-IN line (Lin et al., 2014) obtained with hydro-tRNAseq (Gogakos et al., 2017) or QuantM-tRNAseq (Pinkard et al., 2020). To distinguish experimental from computational differences, the published datasets was also re-analyzed using the computational pipeline as described herein (
mim-tRNAseq Improves tRNA coverage and abundance estimates
[0144] We extended our analysis to single-cell and multicellular eukaryotes by preparing mim-tRNAseq libraries from exponentially growing S. cerevisiae and S. pombe, as well as D. melanogaster BG3-c2 cells and human induced pluripotent stem cells (hiPSCs) with a normal karyotype. The optimal cluster ID threshold was determined as 0.90 for budding yeast and 0.95 for fission yeast, Drosophila, and human tRNA pools. These settings yielded between 85% and 91% of uniquely mapped reads (
[0145] Most reads in mim-tRNAseq datasets mapped to cytosolic tRNA, with mitochondrial tRNA fractions ranging from 0.5% in budding yeast to 3% in hiPSCs. Importantly, nearly all mapped reads (>96%) spanned the post-transcriptionally added 3-CCA stretch, indicating that they originate from mature tRNA. This was not due to bias toward A-ending RNA species, as our workflow accurately captured the 3:1 ratio of two synthetic E. coli tRNA-Lys-UUU tRNAs with either 3-CCA or 3-CC spiked in prior to library construction. cDNA circularization also did not introduce appreciable length bias, as tRNA coverage after alignment mirrored initial cDNA size (
[0146] We asked whether these experimental and computational advances would enable more accurate tRNA quantitation. It was first sought to compare our measurements of absolute tRNA abundance with data obtained with an orthogonal, hybridization-based approach. Absolute RNA quantitation by, for example, Northern blotting or arrays requires highly specific probes and careful comparisons of signal in serial sample dilutions to calibration curves with known target amounts. The design of specific probes for tRNAs, however, is extremely challenging: even with full-length probes, a difference of at least 8 nt is required to avoid cross-hybridization (Dittmar et al., 2004, 2006). Probe design is particularly problematic for human tRNA pools, which can contain >400 tRNA species from 57 anticodon families. As the major tRNA transcript for each anticodon family can differ among cell types (Ishimura et al., 2014), probe selection can unduly influence measurement accuracy. In contrast, the 41 anticodon families of S. cerevisiae consist of only 54 tRNA species, and most major anticodon variants differ sufficiently in sequence to be distinguished by hybridization. Therefore fluorescence intensity measurements for 39 of the 41 budding yeast anticodon families obtained by direct hybridization to a tRNA microarray (Tuller et al., 2010) were compared to the fraction of reads mapping to those anticodon families in mim-tRNAseq datasets. This comparison yielded Pearson's r=0.75 (p=3.810.sup.8), corroborating the quantitative nature of mim-tRNAseq.
[0147] The main regulatory elements for tRNA transcription are intrinsic and overlap with conserved structural regions of mature tRNAs, and it remains unclear how selective tRNA gene expression is achieved in metazoans (Ishimura et al., 2014; Kutter et al., 2011; Schmitt et al., 2014). In rapidly growing yeast cells, however, nearly all tRNA loci are transcribed (Harismendy et al., 2003; Roberts et al., 2003). tRNA gene copy number thus positively correlates with the abundance of tRNA anticodon families during exponential growth measured by hybridization (R.sup.2=0.47 in microscale thermophoresis and R.sup.2=0.60 in tRNA microarray; Jacob et al., 2019; Tuller et al., 2010). The superior resolution of mim-tRNAseq were leveraged to probe this relationship at the level of individual tRNA transcripts (
[0148] The correlation between gene copy number and tRNA abundance was also lower in Drosophila BG3-c2 cells (adjusted R.sup.2=0.79) and hiPSCs (adjusted R.sup.2=0.62). The values were similar regardless of whether copy numbers for all predicted human tRNA genes or only the high-confidence tRNA gene set were used. These findings are consistent with differential tRNA gene use in distinct cell types (Dittmar et al., 2006; Ishimura et al., 2014; Kutter et al., 2011; Schmitt et al., 2014) and highlight that mechanisms beyond gene copy number shape metazoan tRNA pools.
mim-tRNAseq Captures Differences in tRNA Abundance and Aminoacylation
[0149] To establish whether mim-tRNAseq can accurately detect differences in tRNA abundance, first the tRNA pools of karyotypically normal hiPSCs were compared with those in two aneuploid human cell lines (K562 and HEK293T). Of the 368 cytosolic tRNA species resolved quantitatively by mim-tRNAseq, 205 were undetectable in one or more cell lines (0.005% of tRNA-mapped reads). Remarkably, more than half of the detectable tRNAs were differentially expressed, some by up to 3 orders of magnitude (adjusted p<0.05;
[0150] We validated the changes in relative abundance by Northern blotting for two tRNA species: tRNA-Arg-UCU-4 and tRNA-Gly-CCC-2, which differ sufficiently from their isodecoders to avoid probe cross-hybridization and represent tRNAs with low and high abundance. tRNA-Arg-UCU-4 and its mouse ortholog are highly expressed in the central nervous system and are also present at low levels in HEK293T cells (Ishimura et al., 2014; Torres et al., 2019). mim-tRNAseq detected 6- to 8-fold lower levels of tRNA-Arg-UCU-4 in K562 and hiPSCs versus HEK293T, and a similar 5- to 10-fold decrease was observed by Northern blotting (
[0151] We then confirmed the ability of mim-tRNAseq to accurately measure tRNA aminoacylation. Charged tRNAs have periodate-resistant 3 ends and can be quantified as a fraction of tRNAs with 3-CCA versus 3-CC following oxidation and -elimination (Evans et al., 2017). mim-tRNAseq data from oxidized tRNA of wild-type (WT) yeast and a trm.sup.7 strain were compared, which has a tRNA-Phe-GAA charging defect (Han et al., 2018). This defect was evident by a 2.5-fold decrease in 3-CCA proportions for both tRNA-Phe-GAA isodecoders in tRNA pools from trm.sup.7 cells in the absence of other changes in aminoacylation status (
Improved Readthrough Facilitates the Discovery and Annotation of Watson-Crick Face tRNA modifications
[0152] Mismatches to reference and/or premature RT stop signatures are frequently used to detect Watson-Crick face RNA modifications (Clark et al., 2016; Ebhardt et al., 2009; Hauenschild et al., 2015; Katibah et al., 2014; Li et al., 2017; Motorin et al., 2007; Qin et al., 2016; Ryvkin et al., 2013; Safra et al., 2017; Zheng et al., 2015), but their analysis is prone to both experimental and computational artifacts (Sas-Chen and Schwartz, 2019). As tRNA-derived reads are particularly misalignment-prone with standard algorithms (
[0153] In contrast, mim-tRNAseq abrogated nearly all RT stops and yielded reproducibly high levels of mismatches coinciding with frequently modified tRNA positions (
[0154] The only RT blocks remaining in mim-tRNAseq datasets were at rare hypermodified positions. These include 2-methylthio-derivatives of A37 (ms.sup.2t.sup.6A/ms.sup.2i.sup.6A in human cytosolic tRNA-Lys-UUU and 3-4 mitochondrial tRNAs in Drosophila and human cells) and rare stretches of two modified sites (m.sup.2.sub.2G26/27 and 20/20a N.sup.3-[3-amino-3-carboxypropyl]-uridines [acp.sup.3U];
[0155] We then examined whether different modifications are marked by specific signatures of nucleotide misincorporation. This can depend on the processivity and fidelity of an RT, the reaction conditions, and the sequence context of the modified site (Hauenschild et al., 2015; Katibah et al., 2014; Li et al., 2017; Qin et al., 2016; Safra et al., 2017). Signature analysis is especially challenging when RT stops are pervasive, as mismatches at read ends stemming from non-templated nucleotide addition during RT may manifest as misincorporation and lead to spurious modification calls (Sas-Chen and Schwartz, 2019). As mim-tRNAseq enables near-complete modification readthrough (
[0156] To validate the specificity of these signatures, misincorporation patterns at G37 in tRNA-Phe-GAA from WT and trm.sup.7 yeast were compared (
[0157] This remarkable consistency enables the use of misincorporation signatures not only for mapping RNA modifications but also for predicting their identity. The datasets from S. cerevisiae, S. pombe, and Drosophila BG3-c2 cells and human cells were therefore probed for misincorporation-inducing modifications not annotated in MODOMICS. Such sites were identified by a mismatch frequency of >10% and the presence of a distinct misincorporation signature to limit spurious modification calls due to genomic misannotation or SNPs. Modification type was then predicted by combining information on the canonical tRNA position, nucleotide type, and misincorporation signature in comparison with known sites (
Accurate Quantitation of RNA Modification Stoichiometry on the Basis of Misincorporation Rates
[0158] Proportions of RT stops and/or misincorporations are widely used to estimate RNA modification levels (Arimbasseri et al., 2015; Clark et al., 2016; Gogakos et al., 2017; Ryvkin et al., 2013), but whether such measurements are quantitative is unknown. Misincorporation rates at individual modified positions in mim-tRNAseq datasets varied remarkably across tRNA species (
[0159] These findings enabled us to profile modified tRNA fractions with single-transcript resolution in cells from four eukaryotic species. Misincorporation rates were 100% at all instances of 134 and of wyosine derivatives at position 37, suggesting these modifications are present in stoichiometric levels. A similar trend was observed for m.sup.2.sub.2G26, with a clear separation between a majority of fully modified tRNAs and a very small number of transcripts with 10%-30% misincorporation. In contrast, the modified fractions of m.sup.1G, m.sup.3C, and m.sup.1A varied substantially among individual tRNAs independently of sequence context. Instances of very high misincorporation were detectable for all three modifications (m.sup.1A, 100%; m.sup.3C, 94%; m.sup.1G, 88%), indicating that mim-tRNAseq can capture high stoichiometry at these sites if it is present (
[0160] In contrast to the regulatory roles of anticodon loop modifications, m.sup.1A58 is important for the maturation and stability of initiator tRNA-Met in yeast (Anderson et al., 1998) and may play a similar role in other eukaryotic tRNA species. A sequence comparison of budding yeast tRNAs with high or low m.sup.1A58 levels revealed no notable differences, however, indicating that sequence alone is unlikely to be a major determinant of modification stoichiometry at this position.
[0161] To examine whether the stoichiometry of misincorporation-inducing tRNA modifications differs in distinct cell types or states, log odds ratios of misincorporation proportions were calculated across all tRNA positions (see STAR methods). There were very few statistically significant changes when comparing mim-tRNAseq datasets from hiPSCs and HEK293T or K562 cells, suggesting that most tRNAs are modified to a similar extent in these cell lines. A comparison of datasets from WT and trm.sup.104 or trm.sup.14 yeast, however, revealed the striking precision of our approach in detecting transcripts with large reductions in m.sup.1G9 or m.sup.2.sub.2G26. Unexpectedly, in trm.sup.14 yeast cells that lack m.sup.2.sub.2G26, there were also differences in modification levels at other tRNA sites. These included a 3- to 6.5-fold increase in m.sup.1G9 levels in four tRNAs (tRNA-Lys-CUU-1, tRNA-Thr-AGU-1, tRNA-Arg-CCU-1, and tRNA-Asn-GUU-1) and a 2-fold decrease in m.sup.3C32 of tRNA-Ser-UGA-1. m.sup.1G9 levels in tRNA-Lys-CUU-1 and tRNA-Thr-AGU-1 also increase upon Trm.sup.10 overexpression in yeast (Swinehart et al., 2013). Sequence comparisons between tRNAs with increased versus unchanged m.sup.1G9 levels in trm.sup.14 cells indicate that a U7: A66 pair rather than G7: C66 pair may be linked to m.sup.1G9 hypermethylation in the absence of m.sup.2.sub.2G26. These findings reveal an interdependence between Watson-Crick face modifications at distinct tRNA sites and suggest that their stoichiometry is determined by structural features.
Example 3Discussion
[0162] The abundance, charging, and modification status of individual tRNA species can differ in distinct cellular environments. Measuring these properties on a global scale, however, has not been feasible because of technical limitations. No library construction method so far has allowed the efficient reverse transcription of these highly modified RNAs, while the lack of computational tools suited to the complexity of tRNA sequencing data has been another major methodological gap.
[0163] We describe conditions that permit near-complete tRNA modification readthrough by TGIRT, dramatically improving cDNA yield and the fraction of full-length products from tRNA templates. All but one rare tRNA modification roadblock are resolved by mim-tRNAseq, which alleviates the bias of existing tRNA quantitation methods toward low-modified tRNAs species. Our library construction protocol circumvents the need to purify enzymes for modification removal (Zheng et al., 2015) or RT (Zhou et al., 2019), which can introduce unwanted variation. Also described are multiple conceptual advances in the analysis of tRNA sequencing data, including the use of modification annotations, which permits position-specific mismatch tolerance during read alignment. Collectively, these advances enable the efficient and accurate mapping and analysis of tRNA-derived reads with single-transcript resolution.
[0164] One poignant example of the substantial improvements in our computational workflow concerns IRNAs with 134, which is essential for wobble pairing during decoding. Inosines are interpreted as cytosines during RT, resulting in the stoichiometric presence of G in sequencing reads. When using Bowtie or Bowtie 2 to align tRNA datasets from human cells, it was found that reads with G34 were frequently mapped to nearly identical tRNA isoacceptors with U34. Such misalignment can have wide-ranging implications, as it would not only skew abundance estimates but can also lead to spurious conclusions about tRNA modification status and stoichiometry. These findings highlight the importance of both sensitivity and accuracy of read alignment in the context of analyzing tRNA transcriptomes.
[0165] The robust misincorporation signatures deposited by TGIRT reveal the location, type, and stoichiometry of Watson-Crick face base modifications in tRNA. Calibration measurements of observed versus expected modified fractions in existing approaches for sequencing-based modification analysis are either lacking (Ryvkin et al., 2013; Zheng et al., 2015) or display a non-linear relationship (Zhou et al., 2019), likely because of persistent RT stops. In contrast, mim-tRNAseq enables efficient readthrough of almost all tRNA modifications, while modification ID is also discernible by highly specific misincorporation patterns. Improved readthrough permits accurate measurements of modification stoichiometry from misincorporation rates alone, evident from calibration curves with near-perfect linear regression for m.sup.1G and m.sup.2.sub.2G (R.sup.2=0.97). Performing this calibration with mixtures of endogenously modified tRNA pools shows that our entire workflow is free of bias toward low-modified tRNAs.
[0166] Remarkably, although some tRNA positions are almost always fully modified (e.g., m.sup.2.sub.2G26 and 134), others are sub-stoichiometric in some tRNA species. This is in line with a model in which some modifications are deposited because of overlapping substrate specificities in RNA modification enzymes (Phizicky and Alfonzo, 2010). Indeed, methylation at G9 in some yeast tRNAs is enhanced when they lack m.sup.2.sub.2G26, while methylation of C32 is decreased, suggesting that a tRNA conformational change upon m.sup.2.sub.2G26 loss (Steinberg and Cedergren, 1995) might change the affinity of other modification enzymes for individual tRNAs.
[0167] In summary, mim-tRNAseq is a sensitive and accurate start-to-finish technique for quantitation of tRNA abundance and charging, which also reports on the presence and stoichiometry of misincorporation-inducing RNA modifications. The robust library construction workflow and the easy-to-use and freely available computational toolkit make mim-tRNAseq broadly applicable for studying key aspects of tRNA biology in a range of organisms and cell types. Our experimental workflow can also be implemented for the discovery and quantitation of modified sites in other RNA species.
Example 4Material and Methods
TABLE-US-00001 REAGENT or SOURCE IDENTIFIER RESOURCE Commercial reagents mTeSR1 STEMCELL Cat# 85850 Technologies Micro Bio-Spin P30 Bio Rad Cat# 7326251 columns, RNase-free Glycogen Ambion Cat# AM9510 T4 Polynucleotide Kinase New England Cat# M0201L Biolabs T4 RNA ligase 2 New England Cat# M0373L (truncated KQ) Biolabs SUPERase In Ambion Cat# AM2694 TGIRT InGex Cat# TGIRT50 Superscript III Invitrogen Cat# 18080044 AMV RT Promega Cat# M9004 CircLigase ssDNA ligase Lucigen Cat# CL4115K
TABLE-US-00002 REAGENT or RESOURCE SOURCE IDENTIFIER KAPA HiFi DNA Roche Cat# KK2102 Polymerase DNA Zymo Cat# D4013 Clean&Concentrator-5 Research PCR purification kit Immobilon NY+ Millipore Cat# INYC00010 Deposited data Raw and analyzed This paper GEO: GSE152621 sequencing data DM-tRNAseq raw data Zheng et al., GEO: GSE66550 for H. sapiens HEK293T 2015 Hydro-tRNAseq raw data Gogakos GEO: GSE95683 for H. sapiens HEK293 T- et al., 2017 Rex Flp-IN QuantM-tRNAseq raw Pinkard GEO: GSE141436 data et al., 2020 for H. sapiens HEK293 T- Rex Flp-IN
Experimental Models: Cell Lines
TABLE-US-00003 D. melanogaster BG3-c2 P. Becker, N/A cells LMU HEK293T cells O. N/A Griesbeck, MPIN HPSI0214i-kucg_2 cells Kilpinen Cat# 77650065 et al., 2017; ECACC
Experimental Models: Organisms/Strains
TABLE-US-00004 S. cerevisiae: strain Euroscarf N/A BY4741 S. cerevisiae: strain Euroscarf N/A BY4741 trm1::kanMX S. cerevisiae: strain Euroscarf N/A BY4741 trm7::kanMX S. cerevisiae: strain Euroscarf N/A BY4741 trm10::kanMX S. pombe: strain ED668 S. Braun, N/A h+ LMU
Oligonucleotides
TABLE-US-00005 RNA sequences, primers This paper N/A for library construction, and probes for primer
TABLE-US-00006 REAGENT or RESOURCE extension and Northern IDENTIFIER blotting, SOURCE Software and algorithms mim-tRNAseq v0.2.5.6 This paper https://github.com/nedialkova-lab/mim-tRNAseq Bowtie v1.2.2 Langmead http://bowtie-bio.sourceforge.net/index.shtml et al., 2009 Bowtie2 v2.3.3.1 Langmead http://bowtie-bio.sourceforge.net/bowtie2/index.shtml and Salzberg, 2012 GSNAP v2019-02-26 Wu and http://research-pub.gene.com/gmap/ Nacu, 2010 Samtools v1.11 Li et al., 2009 http://samtools.sourceforge.net/ Bedtools v2.29.2 Quinlan and https://bedtools.readthedocs.io/en/latest/ Hall, 2010 BLAST+ v2.9.0 Camacho https://blast.ncbi.nlm.nih.gov et al., 2009 Infernal v1.1.2 Nawrocki and http://eddylab.org/infernal/ Eddy, 2013 usearch Edgar, 2010 https://www.drive5.com/usearch/ v10.0.240_i86linux32 R/DESeq2 v1.26.0 Love et al., https://bioconductor.org/packages/release/bioc/html/DESeq2.html 2014 R/ComplexHeatmap Gu et al., https://www.bioconductor.org/packages/release/bioc/html/ComplexHeatmap.html v2.2.0 2016 Python/Biopython v1.70 Cock et al., https://biopython.org/ 2009
Data and Code Availability
[0168] The accession number for the sequencing data reported in this paper is GEO: GSE152621. The mim-tRNAseq computational pipeline is available under a GNU public License v3 at https://github.com/nedialkova-lab/mim-tRNAseq. A package description and installation guide are available at https://mim-trnaseq.readthedocs.io/en/latest/.
Experimental Model and Subject Details
Cell Lines and Strains
[0169] S. cerevisiae cells (BY4741 wild-type, trm.sup.7, trm.sup.14 and trm 104) were grown in yeast extract-peptone-dextrose (YPD) medium. S. pombe cells (ED668 h+, ade6-M216 ura4-D18 leu1-32) were cultured in yeast extract with supplements (YES). Overnight cultures were diluted to an optical density 600 (OD.sub.600) of 0.05, grown at 30 C. at 250 revolutions per minute, and harvested at OD.sub.600=0.5 by rapid filtration and snap-freezing in liquid nitrogen. D. melanogaster BG3-c2 cells were cultured at 26 C. in Schneider's Drosophila Medium (GIBCO) supplemented with 10% fetal calf serum, 1% penicillin/streptomycin, and 10 g/ml human insulin. HEK293T cells were grown at 37 C. and 5% CO.sub.2 in DMEM supplemented with 10% fetal bovine serum (Sigma Aldrich). The HPSI0214i-kucg_2 human induced pluripotent stem cell line (obtained from HipSci; Kilpinen et al., 2017) was cultured at 37 C. and 5% CO.sub.2 in mTeSR1 (STEMCELL Technologies). K562 cells were grown at 37 C. and 5% CO.sub.2 in RPMI 1640 supplemented with 10% fetal calf serum and 2 mM L-Glutamine.
RNA Isolation
[0170] RNA from Drosophila BG3-c2, HEK293T, and human iPS cells was isolated with Trizol (Sigma Aldrich) according to the manufacturer's instructions. For total RNA isolation from yeast, frozen cells were resuspended in 100 mM sodium acetate pH=4.5, 10 mM EDTA pH=8.0, 1% SDS (1 mL per 50 OD.sub.600 units). An equal volume of hot acid phenol (pH=4.3) was added, and the cell suspension was vortexed vigorously followed by incubation at 65 C. for 5 min (S. cerevisiae) or 45 min (S. pombe) with intermittent mixing. After addition of 1/10 volume 1-Bromo-3-chloropropane (BCP, Sigma Aldrich), samples were centrifuged at 10,000g for 5 min and the aqueous phase was transferred to a new tube. Following an additional round of hot acid phenol/BCP and a round of BCP only extraction, RNA was precipitated from the aqueous phase by the addition of 3 volumes of 100% ethanol. Pellets were washed in 80% ethanol, briefly air-dried, and resuspended in RNase-free water. For RNA isolation from yeast under conditions that preserve tRNA charging, frozen cells were resuspended in ice-cold 100 mM sodium acetate pH=4.5, 10 mM EDTA pH=8.0. One volume of cold acid phenol (pH=4.3) was added and cells were lysed with 500 m-diameter glass beads by three rounds of vortexing for 45 s with a 1-min incubation on ice in between. One-tenth volume of BCP was then added and the samples were centrifuged at 10,000g/4 C. for 5 min, followed by a second round of cold phenol-BCP and one round of BCP-only extraction. RNA was ethanol-precipitated from the aqueous phase and pellets were washed in 80% ethanol containing 50 mM sodium acetate, pH=4.5, briefly air-dried, and resuspended in 50 mM sodium acetate pH=4.5, 1 mM EDTA PH=8.0. RNA concentration was determined with NanoDrop and samples were frozen at 80 C. in single-use aliquots.
RNA Oxidation and -Elimination
[0171] To measure tRNA charging levels, RNA oxidation and -elimination were performed as described (Evans et al., 2017) with minor modifications. 25 g of total RNA were resuspended in 10 mM sodium acetate pH 4.5 and oxidized by the addition of freshly prepared NalO.sub.4 to a final concentration of 50 mM in a 58 L volume for 30 min at 22 C. The reaction was quenched by addition of 6 L 1 M glucose for 5 min at 22 C. RNA was purified with Micro Bio Spin P30 columns (BioRad) followed by two rounds of ethanol precipitation in the presence of 0.3M sodium acetate pH=4.5. Pellets were resuspended in 20 L RNase-free water and -elimination was performed by addition of 30 l 100 mM sodium borate pH=9.5 (freshly prepared) for 90 min at 45 C. RNA was recovered with Micro Bio Spin P30 columns followed by ethanol precipitation, resuspended in RNase-free water, quantified on a NanoDrop, and stored at 80 C. in single-use aliquots.
tRNA Purification by Gel Size Selection
[0172] Two synthetic RNA standards corresponding to E. coli tRNA-Lys-UUU with intact 3-CCA (5-GGGUCGUUAGCUCAGUUGGUAGAGCAGUUGACUUUUAAUCAAUUGGUCGCAGGUU CGAAUCCUGCACGACCCACCA-3) or a 3-CC (5-GGGUCGUUAGCUCAGUUGGUAGAGCAGUUGACUUUUAAUCAAUUGGUCGCAGGUU CGAAUCCUGCACGACCCACC-3) were added to total RNA in a 3:1 molar ratio at 0.06 pmol/g, followed by incubation at 37 C. in 50 mM Tris-HCl PH=9.0 to deacylate tRNAs. Deacylation was omitted for samples subjected to oxidation and -elimination. Total RNA was subsequently dephosphorylated with 10 U of T4 PNK (NEB) at 37 C. for 30 min and purified by ethanol precipitation in 0.3M sodium acetate pH=4.5 with 25 g glycogen (Ambion) as a carrier. RNA was resolved on a denaturing 10% polyacrylamide/7M urea/1TBE gel alongside Low Range ssRNA marker (NEB) and visualized with SYBR Gold. Species migrating at the size range of mature tRNAs (60-100 nt) were excised and gel slices were crushed with disposable pestles. Low-retention tubes and tips (Biotix, Axygen) were used for all subsequent steps of sequencing library construction to maximize nucleic acid recovery. Following addition of 400 l gel elution buffer (0.3M sodium acetate pH=4.5, 0.25% SDS, 1 mM EDTA pH=8.0), the gel slurry was incubated at 65 C. for 10 min, snap-frozen on dry ice, and thawed at 65 C. for 5 min. RNA was eluted overnight at room temperature with continuous mixing. Gel pieces were removed with Costar Spin-X centrifuge tube filters and RNA was recovered from the flow-through by ethanol precipitation in the presence of 25 g of glycogen. This protocol typically recovers 5%-10% of total RNA in the 60-100 nt fraction, consistent with estimates of tRNA proportions in cells (Warner, 1999).
3 Adapter Ligation
[0173] 50 to 200 ng of gel-purified tRNA was ligated to one of four adapters with distinct barcodes: [0174] I1: 5-pGATATCGTCAAGATCGGAAGAGCACACGTCTGAA/ddC/-3; [0175] I2: 5-pGATAGCTACAAGATCGGAAGAGCACACGTCTGAA/ddC/-3; [0176] I3: 5-pGATGCATACAAGATCGGAAGAGCACACGTCTGAA/ddC/-3;
[0177] I4: 5-pGATTCTAGCAAGATCGGAAGAGCACACGTCTGAA/ddC/-3 (barcodes italicised; underlined sequence complementary to RT primer). The adapters are blocked by the 3 chain terminator dideoxycytidine to prevent concatemer formation, and 5-phosphorylated to enable pre-adenylation by Mth RNA ligase prior to ligation (McGlincy and Ingolia, 2017). Ligation was performed for 3 hours at 25 C. in a 20-l reaction volume containing pre-adenylated adapter and RNA substrate in a 4:1 molar ratio, 1T4 RNA Ligase Reaction Buffer, 200 U of T4 RNA ligase 2 (truncated KQ; NEB), 25% PEG 8000, and 10 U SUPERase In (Ambion). Ligation products were separated from excess adapter on denaturing 10% polyacrylamide/7M urea/1TBE gels. Bands migrating at 95-125 nt were excised and ligation products were recovered from crushed gel slices.
Reverse Transcription
[0178] All reactions contained 125 nM primer, 125 nM template and 500 nM TGIRT (InGex) or 200 U Superscript III (Invitrogen). To prime reverse transcription in template-switching reactions, a synthetic RNA/DNA duplex with a single-nucleotide 3 overhang was generated by annealing an RNA oligonucleotide (5-GAGCACACGUCUGAACUCCACUCUUUCCCUACACGACGCUCUUCCGAUCU-3) to a DNA oligonucleotide (5-PRAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGGAGTTCAGACGTGTGCTCN-3). The DNA oligonucleotide contained a phosphorylated A/G at its 5 end, which is a preferred substrate for CircLigase used in subsequent cDNA circularization (Heyer et al., 2015; McGlincy and Ingolia, 2017). For primer-dependent reverse transcription reactions, adapter-ligated RNA and RT primer (5-pRNAGATCGGAAGAGCGTCGTGTAGGGAAAGAG/iSp18/GTGACTGGAGTTCAGACGTG TGCTC-3; underlined sequence complementary to 3 adapter, 5-RN to ameliorate potential biases during circularization) were mixed in MAXYMum Recovery PCR Tubes (Axygen), denatured at 82 C. for 2 min and annealed at 25 C. for 5 min in a Thermocycler. TGIRT reactions were assembled in a 20-l final volume by combining template and primer with 10 U SUPERase In, 5 mM DTT (from a freshly made 100 mM stock) and manufacturer-recommended TGIRT buffer (20 mM Tris-HCl PH=7.6, 450 mM NaCl, 5 mM MgCl.sub.2) or low salt buffer (50 mM Tris-HCl PH=8.3, 75 mM KCl, 3 mM MgCl.sub.2). After TGIRT addition, samples were pre-incubated at reaction temperature for 10 min (primer-dependent reactions) or 22 C. for 30 min (template-switching reactions), initiated by addition of dNTPs to a final concentration of 1.25 mM, and incubated in a Thermocycler for 1 hour or 16 hours. For Superscript III RT, template and primer were denatured at 75 C. for 5 min and chilled on ice, and reverse transcription was performed in the presence of 1 First-Strand Buffer, 5 mM DTT, 0.5 mM dNTPs, 10 U SUPERase In, and 200 U Superscript III (Invitrogen) at 57 C. for 60 min.
[0179] Template RNA was subsequently hydrolyzed by the addition of 1 l 5M NaOH and incubation at 95 C. for 3 min and reaction products were separated from unextended primer on denaturing 10% polyacrylamide/7M urea/1TBE gels. Gels were stained with SYBR Gold, the region between 60 and 150 nt was excised and cDNA was eluted from crushed gel slices in 400 l 10 mM Tris-HCl PH=8.0, 1 mM EDTA at 70 C./2000 rpm for 1 hour in a Thermoblock, followed by ethanol precipitation in 0.3M sodium acetate pH=5.5 in the presence of 25 g glycogen.
cDNA Circularization and Library Construction PCR
[0180] Purified cDNA was circularized with CircLigase ssDNA ligase (Lucigen) in 1 reaction buffer supplemented with 1 mM ATP, 50 mM MgCl.sub.2, and 1M betaine for 3 hours at 60 C., followed by enzyme inactivation for 10 min at 80 C. One-fifth of circularized cDNA was directly used for library construction PCR with a common forward (5 AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCT.Math.C-3) and (5-unique indexed reverse primers CAAGCAGAAGACGGCATACGAGATNNNNNNGTGACTGGAGTTCAGACGTGT.Math.G-3, asterisks denote a phosphorothioate bond and NNNNNN corresponds to the reverse complement of an Illumina index sequence). Amplification was performed with KAPA HiFi DNA Polymerase (Roche) in 1 GC buffer with initial denaturation at 95 C. for 3 min, followed by five to six cycles of 98 C. for 20 s, 62 C. for 30 s, 72 C. for 30 s at a ramp rate of 3 C./sec. PCR products were purified with DNA Clean&Concentrator 5 (Zymo Research) and resolved on 8% polyacrylamide/1TBE gels alongside pBR322 DNA-Mspl Digest (NEB). The 130-220 bp region of each lane was excised and DNA was eluted from crushed gel slices in 400 l water with continuous mixing at room temperature overnight. After ethanol precipitation in 0.3M sodium acetate pH=5.5 and 25 g glycogen, libraries were dissolved in 10 l 10 mM Tris-HCl pH=8.0, quantified with the Qubit dsDNA HS kit, and sequenced for 150 cycles on an Illumina NextSeq platform.
Northern Blotting
[0181] Two micrograms of total RNA were resolved on denaturing 10% polyacrylamide/7M urea/1TBE gels. RNA was transferred to Immobilon NY+membranes (Millipore) in 1TBE for 40 min at 4 mA/cm.sup.2 on a TransBlot Turbo semi-dry blotting apparatus (Bio-Rad) and crosslinked at 0.04 J in a Stratalinker UV crosslinker. Membranes were incubated at 80 C. for one hour and pre-hybridized in 20 mM Na.sub.2HPO.sub.4 PH=7.2, 5SSC, 7% SDS, 2Denhardt, 40 g/ml sheared salmon sperm DNA at 55 C. for 4 hours. The buffer was exchanged 10 pmol and 5-end .sup.32P-labeled probe (Arg-UCU-4: 5-CGGAACCTCTGGATTAGAAGTCCAGCGCGCTCGTCC-3; Gly-CCC-2: 5-CGGGTCGCAAGAATGGGAATCTTGCATGATAC-3) was added, followed by hybridization at 55 C. overnight. Membranes were washed three times in 25 mM Na.sub.2HPO.sub.4 PH=7.5, 3SSC, 5% SDS, 10Denhardt, once in 1SSC, 10% SDS, and exposed to PhosphorImager screens, which were subsequently scanned on a Typhoon FLA 9000 (GE Healthcare). Band intensity was quantified with ImageQuant (GE Healthcare).
Primer Extension Analysis of m.sup.1G37
[0182] The extent of RT arrest at m.sup.1G37 in tRNA-Leu-UAA and tRNA-Pro-UGG from S. cerevisiae was quantified via primer extension with AMV RT, an enzyme with low processivity at this modification (Werner et al., 2020). The primers were designed to enable a 4-nucleotide extension to m.sup.1G37 (tRNA-Leu-UAA: 5-CGCGGACAACCGTCCAAC-3; tRNA-Pro-UGG: 5-TGAACCCAGGGCCTCT-3) and 5-end-labeled with -32P-ATP. 3 g of total RNA from exponentially growing yeast cells was mixed with 1 pmol end-labeled primer and incubated at 95 C. for 3 min followed by slow cooling to 37 C. RT reactions were assembled by adding 15 U AMV RT (Promega), 0.5 mM dNTPs, 20 U SUPERase In (Ambion) and 1AMV RT buffer in a 5-l volume. Following incubation at 37 C. for 45 min, reactions were stopped by addition of 5 l 2RNA loading dye (47.5% Formamide, 0.01% SDS, 0.01% bromophenol blue, 0.005% Xylene Cyanol, 0.5 mM EDTA), boiled at 95 C. for 5 min, and resolved on a denaturing 15% PAA/7M urea/1TBE gel. The gel was exposed at 80 C. to a PhosphorImager screen, which was scanned on a Typhoon FLA 9000 (GE Healthcare). Band intensity was quantified with ImageQuant (GE Healthcare).
Read Preprocessing
[0183] Sequencing libraries were demultiplexed using cutadapt v2.5 (Martin, 2011) and a fasta file (barcodes.fa) of the first 10 nt for the four different 3 adapters (see 3 adapter ligation above). Indels in the alignment to the adapter sequence were disabled with--no-indels. Following demultiplexing, reads were further trimmed to remove the two 5-RN nucleotides introduced by circularization from the RT primer with-u 2. In both processing steps, reads shorter than 10 nt were discarded using-m 10. Example commands for demultiplexing and 5 nucleotide trimming: [0184] cutadapt--no-indels-a file: barcodes.fa-m 10-0 mix1_{name}_trim.fastq.gz mix.fastq.gz [0185] cutadapt-j 40-m 10-u 2-0 mix1_barcode1_trimFinal.fastq.gz mix1_barcode1_trim.fastq.gz
Modification Indexing and Clustering mim-tRNAseq uses modification data from MODOMICS (Boccaletto et al., 2018) to guide accurate alignment of short reads from tRNAs. A prepackaged set of data is available for S. cerevisiae, S. pombe, C. elegans, D. melanogaster, M. musculus, H. sapiens and E. coli, and can be specified with the-species parameter. For other organisms, mim-tRNAseq requires a fasta file of predicted genomic tRNA sequences (-t) and a tRNAscan-SE out file containing information about tRNA introns (-o), both of which should be obtained from GtRNAdb (Chan and Lowe, 2016) or from running tRNAscan-SE (Lowe and Chan, 2016) on the genome of interest. Lastly, a user-generated sample input file is required which contains two tab-separated columns specifying the path to trimmed tRNAseq reads in fastq format, and the experimental condition of each fastq file. Additionally, a mitochondrial tRNA fasta reference is supplied with the prepackaged data inputs listed above, or may be supplied (-m) for custom genomes as a fasta file obtained from mitotRNAdb (Juhling et al., 2009). mim-tRNAseq automatically removes nuclear-encoded mitochondrial tRNAs (nmt-tRNAs) and tRNA species with undetermined anticodons (where applicable), generates mature, processed tRNA sequences (with appended 3-CCA if necessary, 5-G for tRNA-His, and spliced introns), and fetches species-matched MODOMICS entries accordingly. Transcript sequences are then matched to MODOMICS entries using BLAST in order to index all known instances of residues modified at the Watson-Crick face within each tRNA. An additional modifications file for modifications reported in the literature but not yet added to MODOMICS may be supplied and is automatically processed by the pipeline (e.g., 134 annotation; Arimbasseri et al., 2015; Torres et al., 2015). tRNA clustering is enabled with the--cluster parameter, which utilizes the usearch--cluster_fast algorithm (Edgar, 2010) to cluster tRNA sequences by a user-defined sequence identity threshold (customizable with--cluster-id). Regardless of the chosen threshold, only tRNAs sharing an anticodon are clustered to maintain isoacceptor resolution in cases where tRNA transcripts differ by a single nucleotide in the anticodon. The clusters are re-centered based on the number of identical sequences, and this is used to re-cluster and improve the selection of a representative centroid/parent sequence for each cluster (https://www.drive5.com/usearch/manual7/recenter.html). Polymorphisms between cluster members are recorded, and mismatches at these sites during alignment are tolerated, but they are not included in misincorporation analysis for modified sites. Since inosine is interpreted as a G during reverse transcription, annotated inosines are changed to G in tRNA reference sequences.
Read Alignment and Modification Discovery
[0186] After clustering, reads are aligned using GSNAP to the representative centroid cluster sequences of mature tRNA transcripts. By enabling SNP-tolerant alignment with--snp-tolerance, the indexed modified sites are treated as pseudo-SNPs to allow modification-induced mismatches at these sites in a sequence- and position-specific manner. Soft-clipping during alignment in combination with the GSNAP parameter--ignore-trim-in-filtering=1 ensures that non-templated nucleotide extensions are not counted as mismatches during alignment. Mismatch tolerance outside of indexed SNPs is controlled using the--max-mismatches parameter, where an integer of allowed mismatches per read can be provided, or a relative mismatch fraction of read length between 0.0 and 0.1 can be supplied (default 0.1). If--remap is specified, then misincorporation analysis is performed and new, unannotated modifications are called where--misinc-thresh (total misincorporation proportion at a residue; default is 0.1 or 10%) and--min-cov (minimum total coverage for a cluster) regulate the calling of new modifications, which exclude mismatch sites between cluster members appearing as misincorporations in this analysis. The existing SNP index is then updated with these new sites, and realignment of all reads is performed with a mismatch tolerance set using--remap-mismatches. New potential inosine sites are classified for position 34 where a reference A nucleotide is misincorporated with a G in 95% or more total misincorporation events. Both--remap and--max-mismatches are extremely useful for detecting unknown modifications in poorly annotated tRNAs, subsequently allowing more accurate and efficient read alignment, which improves the results of all downstream analyses. Users should consider a low mismatch tolerance during remap to avoid inaccuracy resulting from lenient alignment parameters. a relative mismatch fraction of 0.075 during remapping (--remap-mismatches 0.075) is recommended. Only uniquely mapped reads are retained for post-alignment analyses.
Read Deconvolution
[0187] This process aims to recapitulate the single-transcript resolution of--cluster-id 1 (see above), but with the alignment accuracy and decreased multi-mapping achieved at lower--cluster-id values. The deconvolution algorithm first searches each cluster of tRNA reference sequences for single-nucleotide differences that distinguish among cluster members. For this, each nucleotide in a reference sequence is assessed for uniqueness at that position when compared to all other reference sequences in the cluster. If a nucleotide is unique in position and identity for a specific tRNA reference in the cluster, it is catalogued. Then, after alignment, each read is assessed for mismatches to the cluster parent to which it was aligned. These are then scanned individually to find potential matches to the previously catalogued set for the cluster which can distinguish unique tRNA references. Based on the presence and identity of a unique distinguishing mismatch, a read is then be assigned to a specific tRNA reference within a cluster. Depending on the organism and/or cluster ID threshold, unique distinguishing mismatches may not always be present for all tRNA references in a cluster. Reads without distinguishing mismatches remain assigned to the cluster parent, which is then marked as not fully deconvolved. Using this algorithm, uniquely aligned reads are assigned to individual tRNA sequences in the reference (where possible) before any of the downstream analyses detailed below. For differential expression analyses of reads summed per tRNA anticodon, read deconvolution is not necessary and therefore not performed.
Modification, RT Stop, Readthrough and 3 CCA End Analyses
[0188] Following read deconvolution, all other mismatched positions for the read are extracted from alignment records in bam files, and converted into positions relative to the unique transcript to which the read was assigned (or the cluster parent if definitive assignment is not possible). The identity of the misincoporated nucleotide is recorded to enable signature analysis, and the counts of mismatches for each of the four nucleotides for all reads with the misincorporation are normalized relative to total read coverage at that position. Stops during reverse transcription are extracted from the alignment start position of each read relative to the reference (5 read ends correspond to cDNA 3 ends during RT) and normalized to total read counts for the unique tRNA. Similarly, readthrough for each position is calculated as the fraction of reads that stop at a position relative to read coverage at each position (as opposed to stop proportions which are normalized to total tRNA read coverage). This value is then subtracted from one to estimate the proportion of reads per position that extend beyond that site, and the minimum value in a 3-nucleotide window centered around the modification is recorded. Using a 3-nucleotide window ensures that potential variance in the position at which the RT stalls due to the modification is accounted for. Taking the minimum value of readthrough for these 3 nucleotides reduces the likelihood of readthrough overestimation. Misincorporation, stop data, and readthrough per unique tRNA sequence, per position are output as tab-separated files, and global heatmaps showing misincorporation and stop proportions across all unique tRNA sequences are plotted per experimental condition. Misincorporation signatures are also plotted for well-known conserved modified tRNA sites (9, 20, 26, 32, 34, 37 and 58) separated by upstream and downstream sequence context to assess potential factors influencing misincorporation signatures. Lastly, the dinucleotide at the 3 ends of reads is quantified, so long as the read aligns to the conserved 3-CCA tail of the reference. Proportions of transcripts with absent 3 tails, 3-C, 3-CC and 3-CCA are calculated per unique tRNA sequence and plotted pairwise between conditions for quantitation and comparison of functional tRNA pools, or tRNA charging fractions in periodate oxidation experiments.
Post-Alignment Analyses
[0189] The cluster deconvolution algorithm allows coverage analysis, novel modification discovery and read counting for tRNA quantitation to be done at the level of unique tRNA sequences. Coverage is calculated as the depth of reads at all positions across a tRNA sequence and plotted using custom R scripts. Cytosolic tRNAs with low read coverage can be filtered at the coverage analysis step by supplying a minimum coverage threshold to--min-cov. Unique tRNA sequences filtered out here are excluded from all downstream analyses, except differential expression analysis by DESeq2 (Love et al., 2014) where all unique tRNA sequences are included. Normalized coverage (read fraction relative to library size) is plotted per sample in 25 bins across gene length in a metagene analysis. Normalized coverage is also scaled relative to the second last bin to account for potential differences in 3 CCA intactness. Read counts per unique tRNA sequence are summed to calculate read counts per isoacceptor family (all tRNAs sharing an anticodon). These counts are subsequently used by a DESeq2 pipeline for count transformations, sample distance analysis using distance matrix heatmaps, PCA plots, and differential expression analysis at the level of isoacceptor families and unique tRNA transcripts (only for completely resolved clusters). In the case that only one experimental condition is supplied, or if there are no replicates for one or more conditions, differential expression analysis is not performed on these samples, but a normalized counts table is still produced for investigations into tRNA abundance.
Data Analysis with the Mim-tRNAseq Package
[0190] The following parameters were used for the analysis of mim-tRNAseq generated sequencing datasets (see mimseq--help or https://mim-trnaseq.readthedocs.io/en/latest/intro.html for full explanations of parameters; package version v0.2.5.6): [0191] S. cerevisiae:--cluster--cluster-id 0.90--snp-tolerance--min-cov 2000--max-mismatches 0.1--control-condition Exp--cca-analysis--remap--remap-mismatches 0.075 [0192] S. pombe:--cluster--cluster-id 0.95--snp-tolerance--min-cov 2000--max-mismatches 0.1--control-condition Exp--cca-analysis--remap--remap-mismatches 0.075 [0193] D. melanogaster:--cluster--cluster-id 0.95--snp-tolerance--min-cov 2000--max-mismatches 0.1--control-condition bg3--cca-analysis--remap--remap-mismatches 0.075 [0194] H. sapiens:--snp-tolerance--cluster--cluster-id 0.95--min-cov 2000--max-mismatches 0.1--control-condition kiPS--cca-analysis--remap--remap-mismatches 0.075
tRNA Read Alignment with Bowtie and Bowtie 2
[0195] To test previously used alignment strategies as in DM-tRNAseq (Zheng et al., 2015) or ARM-seq (Cozen et al., 2015), a non-redundant set of reference human tRNA transcripts was created by fetching the full set of 610 predicted tRNA genes for human genome hg19 from GtRNAdb (Chan and Lowe, 2016) and the 22 mitochondrially encoded human tRNA genes from mitotRNAdb (Juhling et al., 2009). Following intron removal and addition of 3 CCA (for nuclear-encoded transcripts) and 5-G (for tRNA-His), a curated set of 596 genes (excluding anticodon/isotype mismatch and nuclear-encoded mitochondrial tRNAs) were collapsed into 420 unique sequences. Corresponding Bowtie and Bowtie 2 indices were built from this set of references. Bowtie alignment was performed with a maximum of 3 allowed mismatches per read (-v 3), filtering for uniquely aligning reads (-m 1) and ensuring the best alignment from the best stratum (i.e., reads with the least number of mismatches) were reported (--best--strata). Bowtie 2 alignments were performed in very sensitive local mode (--very-sensitive--local) and up to 100 alignments per read were allowed (-k 100). Read quality scores were ignored for alignment score and mismatch penalty calculation (--ignore-quals) with increased penalties for ambiguous characters (N) in reference or read (--np 5). Output alignments in SAM format were reordered to match read order in input fastq file (--reorder). The alignment commands for both algorithms are given below: [0196] bowtie-v 3-m 1--best--strata--threads 40-S [0197] bowtie2--local-x-k 100--very-sensitive--ignore-quals--np 5--reorder-p 40-U [0198] QuantM-tRNAseq data for HEK293 T-Rex Flp-IN cells downloaded from the NCBI Gene Expression Omnibus repository was adapter-trimmed and analyzed with Bowtie 2 as described in (Pinkard et al., 2020): [0199] bowtie2--local--score-min G, 1,8-D 20-R 3-N 1-L 10-i S, 1, 0.5
Sequence Logo Analysis
[0200] Alignment files for uniquely aligned reads from human HEK293T and S. cerevisiae cells were utilized to generate frequency plots of untemplated nucleotide additions by TGIRT, and 5 sequence logos in each sample. Briefly, CIGAR strings for each unique alignment were assessed for GSNAP soft-clipped nucleotides representing untemplated additions. The number of additions per read were recorded and plotted as frequency histograms. Since a total of 3 additions or less were present in >90% of reads analyzed, sequence logos were generated using the Python package Logomaker (Tareen and Kinney, 2020) for these reads using soft-clipped residues and the first 10 nucleotides after them. For the logo representing all cataloged tRNA genes, mature tRNA transcript sequences was used from each genome present in GtRNAdb, and generated a multiple sequence alignment of these using Infernal (Nawrocki and Eddy, 2013). A sequence logo was then generated from the first 11 nucleotides of each aligned tRNA transcript (in order to include G-1 for tRNA-His, plus 10 additional nucleotides as in the uniquely aligned read logo above).
Differential Modification Analysis
[0201] To test for global differential modification between two conditions, first, misincorporation proportion and coverage data generated by mim-tRNAseq were used to calculate absolute counts of modified and unmodified bases per position for each resolved tRNA transcript. Then, log odds ratios (log OR) were calculated for each position, x, as follows:
where M.sub.a and M.sub.b are the counts of modified nucleotides at position x in condition a and b, and U.sub.a and U.sub.b are the counts of unmodified nucleotides at position x in condition a and b, respectively. Significance for each log OR was determined with chi-square tests using the respective modified and unmodified nucleotide counts for each condition in a two-dimensional contingency table for the Pearson's chi-square test. Correction for multiple testing was performed with the FDR method. Following significance tests, log OR values were filtered for FDR-adjusted p values 0.01, absolute log 2 fold-changes 1, and total misincorporation at the given position of 10% or more in at least one of the conditions to ensure only sites with high-confidence misincorporation levels are kept. The resulting log OR were used in generating heatmaps for individual contrasts between cell types or experimental conditions.
REFERENCES
[0202] Anderson J., Phan L., Cuesta R., Carlson B. A., Pak M., Asano K., Bjrk G. R., Tamame M., Hinnebusch A. G. The essential Gcd10p-Gcd14p nuclear complex is required for 1-methyladenosine modification and maturation of initiator methionyl-tRNA. Genes Dev. 1998; 12:3650-3662. [0203] Arimbasseri A. G., Blewett N. H., Iben J. R., Lamichhane T. N., Cherkasova V., Hafner M., Maraia R. J. RNA Polymerase Ill output is functionally linked to tRNA dimethyl-G26 modification. PLOS Genet. 2015; 11:e1005671. [0204] Arimbasseri A. G., Iben J., Wei F. Y., Rijal K., Tomizawa K., Hafner M., Maraia R. J. Evolving specificity of tRNA 3-methyl-cytidine-32 (m.sup.3C32) modification: a subset of tRNAsSer requires N6-isopentenylation of A37. RNA. 2016; 22:1400-1410. [0205] Boccaletto P., Machnicka M. A., Purta E., Piatkowski P., Baginski B., Wirecki T. K., de Crcy-Lagard V., Ross R., Limbach P. A., Kotter A. MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. 2018; 46 (D1): D303-D307. [0206] Camacho C., Coulouris G., Avagyan V., Ma N., Papadopoulos J., Bealer K., Madden T. L. BLAST+: architecture and applications. BMC Bioinformatics. 2009; 10:421. [0207] Chan P. P., Lowe T. M. GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes. Nucleic Acids Res. 2016; 44 (D1): D184-D189. [0208] Chen D., Patton J. T. Reverse transcriptase adds nontemplated nucleotides to cDNAs during 5-RACE and primer extension. Biotechniques. 2001; 30:574-580, 582. [0209] Clark W. C., Evans M. E., Dominissini D., Zheng G., Pan T. tRNA base methylation identification and quantification via high-throughput sequencing. RNA. 2016; 22:1771-1784. [0210] Cock P. J. A., Antao T., Chang J. T., Chapman B. A., Cox C. J., Dalke A., Friedberg I., Hamelryck T., Kauff F., Wilczynski B., de Hoon M. J. L. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics. 2009; 25:1422-1423. [0211] Cozen A. E., Quartley E., Holmes A. D., Hrabeta-Robinson E., Phizicky E. M., Lowe T. M. ARM-seq: AlkB-facilitated RNA methylation sequencing reveals a complex landscape of modified tRNA fragments. Nat. Methods. 2015; 12:879-884. [0212] Dittmar K. A., Mobley E. M., Radek A. J., Pan T. Exploring the regulation of tRNA distribution on the genomic scale. J. Mol. Biol. 2004; 337:31-47. [0213] Dittmar K. A., Goodenbour J. M., Pan T. Tissue-specific differences in human transfer RNA expression. PLOS Genet. 2006; 2: e221-e229. [0214] Ebhardt H. A., Tsang H. H., Dai D. C., Liu Y., Bostan B., Fahlman R. P. Meta-analysis of small RNA-sequencing errors reveals ubiquitous post-transcriptional RNA modifications. Nucleic Acids Res. 2009; 37:2461-2470. [0215] Edgar R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010; 26:2460-2461. [0216] Ellis S. R., Morales M. J., Li J. M., Hopper A. K., Martin N. C. Isolation and characterization of the TRM1 locus, a gene essential for the N2, N2-dimethylguanosine modification of both mitochondrial and cytoplasmic tRNA in Saccharomyces cerevisiae. J. Biol. Chem. 1986; 261:9703-9709. [0217] Evans M. E., Clark W. C., Zheng G., Pan T. Determination of tRNA aminoacylation levels by high-throughput sequencing. Nucleic Acids Res. 2017; 45: e133. [0218] Gamper H. B., Masuda I., Frenkel-Morgenstern M., Hou Y.-M. Maintenance of protein synthesis reading frame by EF-P and m()G37-tRNA. Nat. Commun. 2015; 6:7226. [0219] Gerber A. P., Keller W. An adenosine deaminase that generates inosine at the wobble position of tRNAs. Science. 1999; 286:1146-1149. [0220] Gogakos T., Brown M., Garzia A., Meyer C., Hafner M., Tuschl T. Characterizing expression and processing of precursor and mature human tRNAs by hydro-tRNAseq and PAR-CLIP. Cell Rep. 2017; 20:1463-1475. [0221] Goodenbour J. M., Pan T. Diversity of tRNA genes in eukaryotes. Nucleic Acids Res. 2006; 34:6137-6146. [0222] Gu Z., Eils R., Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016; 32:2847-2849. [0223] Guy M. P., Podyma B. M., Preston M. A., Shaheen H. H., Krivos K. L., Limbach P. A., Hopper A. K., Phizicky E. M. Yeast Trm.sup.7 interacts with distinct proteins for critical modifications of the tRNAPhe anticodon loop. RNA. 2012; 18:1921-1933. [ [0224] Han L., Guy M. P., Kon Y., Phizicky E. M. Lack of 2-O-methylation in the tRNA anticodon loop of two phylogenetically distant yeast species activates the general amino acid control pathway. PLOS Genet. 2018; 14: e1007288. [0225] Harismendy O., Gendrel C.-G., Soularue P., Gidrol X., Sentenac A., Werner M., Lefebvre O. Genome-wide location of yeast RNA polymerase III transcription machinery. EMBO J. 2003; 22:4738-4747. [0226] Hauenschild R., Tserovski L., Schmid K., Thuring K., Winz M.-L., Sharma S., Entian K.-D., Wacheul L., Lafontaine D. L. J., Anderson J. The reverse transcription signature of N-1-methyladenosine in RNA-Seq is sequence dependent. Nucleic Acids Res. 2015; 43:9950-9964. [0227] Heyer E. E., Ozadam H., Ricci E. P., Cenik C., Moore M. J. An optimized kit-free method for making strand-specific deep sequencing libraries from RNA fragments. Nucleic Acids Res. 2015; 43: e2. [0228] Hoffmann A., Fallmann J., Vilardo E., Morl M., Stadler P. F., Amman F. Accurate mapping of tRNA reads. Bioinformatics. 2018; 34:1116-1124 [0229] Ishimura R., Nagy G., Dotu I., Zhou H., Yang X.-L., Schimmel P., Senju S., Nishimura Y., Chuang J. H., Ackerman S. L. RNA function. Ribosome stalling induced by mutation of a CNS-specific tRNA causes neurodegeneration. Science. 2014; 345:455-459. [0230] Jackman J. E., Montagne R. K., Malik H. S., Phizicky E. M. Identification of the yeast gene encoding the tRNA m.sup.1G methyltransferase responsible for modification at position 9. RNA. 2003; 9:574-585. [0231] Jacob D., Thuring K., Galliot A., Marchand V., Galvanin A., Ciftci A., Scharmann K., Stock M., Roignant J.-Y., Leidel S. A. Absolute quantification of noncoding RNA by microscale thermophoresis. Angew. Chem. Int. Ed. Engl. 2019; 58:9565-9569. [0232] M., Stadler P. F., Ptz J. tRNAdb 2009: compilation of tRNA sequences and tRNA genes. Nucleic Acids Res. 2009; 37: D159-D162. [0233] Karaca E., Weitzer S., Pehlivan D., Shiraishi H., Gogakos T., Hanada T., Jhangiani S. N., Wiszniewski W., Withers M., Campbell I. M., Baylor Hopkins Center for Mendelian Genomics Human CLP1 mutations alter tRNA biogenesis, affecting both peripheral and central nervous system function. Cell. 2014; 157:636-650. [0234] Katibah G. E., Qin Y., Sidote D. J., Yao J., Lambowitz A. M., Collins K. Broad and adaptable RNA structure recognition by the human interferon-induced tetratricopeptide repeat protein IFIT5. Proc. Natl. Acad. Sci. USA. 2014; 111:12025-12030. [0235] Kilpinen H., Gonalves A., Leha A., Afzal V., Alasoo K., Ashford S., Bala S., Bensaddek D., Casale F. P., Culley O. J. Common genetic variation drives molecular heterogeneity in human iPSCs. Nature. 2017; 546:370-375. [0236] Kirchner S., Ignatova Z. Emerging roles of tRNA in adaptive translation, signalling dynamics and disease. Nat. Rev. Genet. 2015; 16:98-112. [0237] Kutter C., Brown G. D., Gonalves A., Wilson M. D., Watt S., Brazma A., White R. J., Odom D. T. Pol Ill binding in six mammals shows conservation among amino acid isotypes despite divergence among tRNA genes. Nat. Genet. 2011; 43:948-955. [0238] Langmead B., Salzberg S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods. 2012; 9:357-359. [0239] Langmead B., Trapnell C., Pop M., Salzberg S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10: R25. [0240] Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R., 1000 Genome Project Data Processing Subgroup The sequence alignment/map format and SAMtools. Bioinformatics. 2009; 25:2078-2079. [0241] Li X., Xiong X., Zhang M., Wang K., Chen Y., Zhou J., Mao Y., Lv J., Yi D., Chen X.-W. Base-resolution mapping reveals distinct m.sup.1A methylome in nuclear- and mitochondrial-encoded transcripts. Mol. Cell. 2017; 68:993-1005.e9. [0242] Lin Y.-C., Boone M., Meuris L., Lemmens I., Van Roy N., Soete A., Reumers J., Moisse M., Plaisance S., Drmanac R. Genome dynamics of the human embryonic kidney 293 lineage in response to cell biology manipulations. Nat. Commun. 2014; 5:4767. [0243] Love M. I., Huber W., Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15:550. [0244] Lowe T. M., Chan P. P. tRNAscan-SE On-line: integrating search and context for analysis of transfer RNA genes. Nucleic Acids Res. 2016; 44 (W1): W54-W57. [0245] Maehigashi T., Dunkle J. A., Miles S. J., Dunham C. M. Structural insights into +1 frameshifting promoted by expanded or modification-deficient anticodon stem loops. Proc. Natl. Acad. Sci. USA. 2014; 111:12740-12745.] [0246] Martin M. M. E. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 2011; 17:10-12. [0247] Masuda I., Matsubara R., Christian T., Rojas E. R., Yadavalli S. S., Zhang L., Goulian M., Foster L. J., Huang K. C., Hou Y.-M. tRNA methylation is a global determinant of bacterial multi-drug resistance. Cell Syst. 2019; 8:302-314.e8. [0248] McGlincy N. J., Ingolia N. T. Transcriptome-wide measurement of translation by ribosome profiling. Methods. 2017; 126:112-129. [0249] Mohr S., Ghanem E., Smith W., Sheeter D., Qin Y., King O., Polioudakis D., lyer V. R., Hunicke-Smith S., Swamy S. Thermostable group II intron reverse transcriptase fusion proteins and their use in cDNA synthesis and next-generation RNA sequencing. RNA. 2013; 19:958-970. [0250] Motorin Y., Helm M. Methods for RNA modification mapping using deep sequencing: established and new emerging technologies. Genes (Basel) 2019; 10:35. [0251] Motorin Y., Muller S., Behm-Ansmant I., Branlant C. Identification of modified residues in RNAs by reverse transcription-based methods. Methods Enzymol. 2007; 425:21-53. [0252] Nawrocki E. P., Eddy S. R. Infernal 1.1:100-fold faster RNA homology searches. Bioinformatics. 2013; 29:2933-2935. [0253] Pernod K., Schaeffer L., Chicher J., Hok E., Rick C., Geslain R., Eriani G., Westhof E., Ryckelynck M., Martin F. The nature of the purine at position 34 in tRNAs of 4-codon boxes is correlated with nucleotides at positions 32 and 38 to maintain decoding fidelity. Nucleic Acids Res. 2020; 48:6170-6183. [0254] Phizicky E. M., Alfonzo J. D. Do all modifications benefit all tRNAs? FEBS Lett. 2010; 584:265-271. [0255] Pinkard O., McFarland S., Sweet T., Coller J. Quantitative tRNA-sequencing uncovers metazoan tissue-specific tRNA regulation. Nat. Commun. 2020; 11:4104-4115. [0256] Qin Y., Yao J., Wu D. C., Nottingham R. M., Mohr S., Hunicke-Smith S., Lambowitz A. M. High-throughput sequencing of human plasma RNA by using thermostable group II intron reverse transcriptases. RNA. 2016; 22:111-128. [0257] Quail M. A., Otto T. D., Gu Y., Harris S. R., Skelly T. F., McQuillan J. A., Swerdlow H. P., Oyola S. O. Optimal enzymes for amplifying sequencing libraries. Nat. Methods. 2011; 9:10-11. Quinlan A. R., Hall I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010; 26:841-842. [0258] Roberts D. N., Stewart A. J., Huff J. T., Cairns B. R. The RNA polymerase III transcriptome revealed by genome-wide localization and activity-occupancy relationships. Proc. Natl. Acad. Sci. USA. 2003; 100:14695-14700. [0259] Ryvkin P., Leung Y. Y., Silverman I. M., Childress M., Valladares O., Dragomir I., Gregory B. D., Wang L. S. HAMR: high-throughput annotation of modified ribonucleotides. RNA. 2013; 19:1684-1692. [0260] Safra M., Sas-Chen A., Nir R., Winkler R., Nachshon A., Bar-Yaacov D., Erlacher M., Rossmanith W., Stern-Ginossar N., Schwartz S. The m.sup.1A landscape on cytosolic and mitochondrial mRNA at single-base resolution. Nature. 2017; 551:251-255. [0261] Sas-Chen A., Schwartz S. Misincorporation signatures for detecting modifications in mRNA: Not as simple as it sounds. Methods. 2019; 156:53-59. [0262] Schmitt B. M., Rudolph K. L. M., Karagianni P., Fonseca N. A., White R. J., Talianidis I., Odom D. T., Marioni J. C., Kutter C. High-resolution mapping of transcriptional dynamics across tissue development reveals a stable mRNA-tRNA interface. Genome Res. 2014; 24:1797-1807. [0263] Steinberg S., Cedergren R. A correlation between N.sup.2-dimethylguanosine presence and alternate tRNA conformers. RNA. 1995; 1:886-891. [0264] Swinehart W. E., Henderson J. C., Jackman J. E. Unexpected expansion of tRNA substrate recognition by the yeast m.sup.1G9 methyltransferase Trm.sup.10. RNA. 2013; 19:1137-1146. [0265] Tareen A., Kinney J. B. Logomaker: beautiful sequence logos in Python. Bioinformatics. 2020; 36:2272-2274. [0266] Torres A. G., Pineyro D., Rodrguez-Escriba M., Camacho N., Reina O., Saint-Lger A., Filonava L., Batlle E., Ribas de Pouplana L. Inosine modifications in human tRNAs are incorporated at the precursor tRNA level. Nucleic Acids Res. 2015; 43:5145-5157. [0267] Torres A. G., Reina O., Stephan-Otto Attolini C., Ribas de Pouplana L. Differential expression of human tRNA genes drives the abundance of tRNA-derived fragments. Proc. Natl. Acad. Sci. USA. 2019; 116:8451-8456. [0268] Tuller T., Carmi A., Vestsigian K., Navon S., Dorfan Y., Zaborske J., Pan T., Dahan O., Furman I., Pilpel Y. An evolutionarily conserved mechanism for controlling the efficiency of protein translation. Cell. 2010; 141:344-354. [0269] Warner J. R. The economics of ribosome biosynthesis in yeast. Trends Biochem. Sci. 1999; 24:437-440. [0270] Werner S., Schmidt L., Marchand V., Kemmer T., Falschlunger C., Sednev M. V., Bec G., [0271] Ennifar E., Hobartner C., Micura R. Machine learning of reverse transcription signatures of variegated polymerases allows mapping and discrimination of methylated purines in limited transcriptomes. Nucleic Acids Res. 2020; 48:3734-3746. [0272] Wu T. D., Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics. 2010; 26:873-881. [0273] Xu H., Yao J., Wu D. C., Lambowitz A. M. Improved TGIRT-seq methods for comprehensive transcriptome profiling with decreased adapter dimer formation and bias correction. Sci. Rep. 2019; 9:7953. [0274] Zhao C., Liu F., Pyle A. M. An ultraprocessive, accurate reverse transcriptase encoded by a metazoan group II intron. RNA. 2018; 24:183-195. [0275] Zheng G., Qin Y., Clark W. C., Dai Q., Yi C., He C., Lambowitz A. M., Pan T. Efficient and quantitative high-throughput tRNA sequencing. Nat. Methods. 2015; 12:835-837. [0276] Zhou H., Rauch S., Dai Q., Cui X., Zhang Z., Nachtergaele S., Sepich C., He C., Dickinson B. C. Evolution of a reverse transcriptase to map N1-methyladenosine in human messenger RNA. Nat. Methods. 2019; 16:1281-1288. [0277] Zhuang F., Fuchs R. T., Sun Z., Zheng Y., Robb G. B. Structural bias in T4 RNA ligase-mediated 3-adapter ligation. Nucleic Acids Res. 2012; 40: e54.