Introduction

Epitranscriptomics describes the occurrence of chemical modifications on RNA that regulate stability, localization, splicing, translation or decay (Gilbert and Nachtergaele 2023). To date, over 170 RNA modifications are known to occur across all classes of RNA including eukaryotic mRNA (Cappannini et al. 2023). Of these, at least thirteen different chemical modifications have been characterized in mRNAs, comprising both cap-proximal and internal modifications (Anreiter et al. 2021). By contrast, current evidence indicates that IncRNAs harbor a more limited repertoire of modifications, primarily W⁶- methyladenosine (m⁶A), 5-methylcytosine (m⁵C), and N7-methylguanosine (m⁷G) (Yang et al. 2024a). One of the best studied and most abundant internal mRNA modifications is m⁶A. Functional studies of m⁶A methylases, demethylases and m⁶A binding proteins, as well as sequencing-based mapping methods have increased our understanding of m^sA as a reversible and dynamic modification that fine-tunes nearly all aspects of mRNA metabolism (Zaccara et al. 2019; He and He 2021; Murakami and Jaffrey 2022; Gilbert and Nachtergaele 2023; Sendinc and Shi 2023).

In addition to m⁶A, other widespread and functionally important marks include m⁵C, which is less abundant but influences mRNA export and stability, pseudouridine (ψ), which contributes to mRNA stability, translation and splicing, and adenosine-to-inosine (A to I) editing, which is involved in splicing, stability and recoding events that contribute to protein diversity (Sun et al. 2023).

Techniques based on short-read (Illumina) sequencing have advanced in recent years to achieve single-nucleotide resolution and quantitative mapping of RNA modifications. The approaches that have emerged to map m⁶A can be broadly classified as antibody-based, enzyme-assisted or chemical conversion methods (Moshitch-Moshkovitz et al. 2024).

Early antibody-based techniques (Dominissini et al. 2012; Meyer et al. 2012; Linder et al. 2015) provided the first transcriptome-wide insights into m⁶A architecture, albeit at limited resolution. Enzyme-assisted methods (Meyer 2019; Xiao et al. 2023) have improved precision but can introduce enzymatic biases. These biases are avoided with GLORI-Sequencing (GLORI), which selectively deaminates unmethylated adenosines, enabling theoretically unbiased, quantitative, single-nucleotide mapping of m⁶A (Liu et al. 2023).

Despite advances, the broad utility of these approaches is limited by a lack of isoform specificity, as well as biases associated with amplification, specificity of antibodies or variable adenosine conversion rates (Diensthuber and Novoa 2025). Oxford Nanopore Technologies (ONT) direct RNA sequencing (DRS) offers an alternative to these methods by analyzing native RNA molecules that have retained their modifications (Garalde et al. 2018). As single RNAs pass through nanopores, signal changes are interpreted to determine canonical and modified residues at single-nucleotide resolution on individual transcripts (Jain et al. 2016; Viehweger et al. 2019). For the now discontinued DRS SQK-RNA002 chemistry used with R9 flowcells, multiple research groups developed algorithms for detecting RNA modifications. These methods can be divided into error rate-based and current signal-based (Abebe et al. 2022). While the recent release of the DRS SQK-RNA004 chemistry and RNA-specific flowcells reduced the utility of many of these tools, introduction of a modification-aware basecaller (Dorado) provides an all-in-one toolkit for basecalling, RNA modification detection, and poly(A) tail length estimation. Early versions of Dorado (e.g. v0.5.0 and v0.6.0) enabled m⁶A (DRACH-context only) detection, while version v0.7.0 enabled detection of ψ and all-context m⁶A. More recent versions (v0.8.0 and v0.9.0) allow both DRACH- and all-context m⁶A detection, together with ψ, inosine and m⁵C (Zou et al. 2025).

In a recent overview of this emerging technology, Cruciani and Novoa highlighted the importance of rigorous evaluation of the limitations of Dorado and advocated forcefully for the establishment of best practices to ensure consistency and accuracy going forward (Cruciani et al. 2025). To date, only a handful of studies of this kind have been published. A series of studies from Novoa and colleagues showed that Dorado performs robustly on synthetic modified and unmodified RNA molecules but also found that overlaps in the identification of specific modifications and false-positive calls remain a challenge but can be mitigated by the inclusion of in vitro transcribed (IVT) controls (Cruciani et al. 2025; Diensthuber and Novoa 2025). Independently, Esfahani and colleagues used poly(A) RNA from a B-lymphocyte cell line to compare SQK-RNA002 and SQK-RNAOO4-derived DRS data. While they concluded that Dorado outperformed other RNA modification detection tools, IVT controls were necessary to reduce false-positive calls (Esfahani et al. 2025). Recent benchmarking efforts comparing Dorado and m6Anet on SQK-RNA004 chemistry have shown that both tools achieve high recall and strong stoichiometry correlations at DRACH motifs but exhibit higher false discovery rates on IVT controls and lower recalls of low- stoichiometry genomic m⁶A sites. However, these evaluations relied on previously published GLORI and eTAM-Seq datasets, rather than matched input material (Diensthuber and Novoa 2025; Zou et al. 2025). By directly comparing m⁶A detection by DRS coupled with Dorado and GLORI on the same RNA inputs, our study overcomes this limitation. Specifically, we systemically benchmark multiple versions of Dorado using poly(A) RNA from primary (normal) human dermal fibroblasts (NHDFs) and functional human neurons produced through terminal differentiation of HD10.6 cells (Raymon et al. 1999), along with IVT RNA. We included treatments with STM2457, an enzymatic inhibitor of the primary m⁶A methyltransferase METTL3 (Yankova et al. 2021), to robustly assess false-positive m⁶A detection rates and also compared the DRS outputs with those obtained from GLORI (Yang et al. 2024b) using the same RNA as input. While our analysis is centered on the sensitivity and specificity of m⁶A detection in DRS datasets, we also benchmark the ability of Dorado to reliably detect m⁵C, ψ, and inosine. Together, our data provides a robust workflow for DRS-based modification detection that reliably captures transcriptome-wide, isoform-level, site-specific stoichiometries of RNA modifications.By expanding m⁶A mapping on primary human fibroblasts and functional human neurons, we uncover mostly conserved m⁶A profiles but reveal increased site-specific m⁶A stoichiometries in neuronal cells.

Results

Systematic benchmarking of ONT’s Dorado across human cell types

To evaluate the performance of ONT’s basecalling tool Dorado across different human cell types, we selected primary NHDFs and differentiated HD10.6 cells as a neuron-like model. Prior to isolation of total RNA, cells were treated for 48h with a small molecule inhibitor of METTL3 (STM2457) or with vehicle (DMSO). We then isolated poly(A) RNA and subjected the same input material to both DRS and GLORI. The same material was also used to generate IVT RNA (Fig. 1A). Dorado basecalling using version v.0.9.0 resulted in over 13M reads using poly(A) RNA from NHDFs treated with DMSO or STM2457, while we obtained over 5M and 10M reads from DMSO- and STM2457-treated HD10.6 cells, respectively. The rate of primary alignment to the human genome (hg38) was over 98% for all the sequenced samples. A detailed overview of read statistics generated by DRS can be found in Supplementary Table 1.

Figure 1.

To systematically benchmark Dorado, we performed basecalling using versions v0.5.0 through v0.9.0. Reads were aligned to the human genome and transcriptome to obtain comprehensive information on isoform-level modification detection. Further, we explored different aspects of detected modifications using the outputs generated by ONT’s Modkit tool. We assessed multiple filtering strategies using base and modification probabilities, as well as stoichiometry distributions, and evaluated them based on the IVT data and sensitivity to METTL3 inhibition. By filtering out fa Ise-positive calls using multiple control layers, we determined high-confidence m⁶A sites (Fig. IB). Finally, we compared these high-confidence predictions with sites identified by GLORI performed on the same input RNA. Using these results, we developed a benchmarked pipeline that filters m⁶A sites predicted by Dorado to minimize false-positive detection without impacting sensitivity (Fig. IC).

Stringent filtering of modification probabilities removes false positives

During modification-aware basecalling, each nucleotide is assigned a basecall probability score ranging from 0 to 1. Similarly, individual modified nucleotides are assigned a modification probability score in the same range. To establish a baseline for false-positive m⁶A predictions, we plotted basecall and modification probability distributions for an IVT poly(A) RNA dataset processed using Dorado v0.9.0 and aligned to human genome (Fig.2). We binned basecall and modification probability values for A and m⁶A (all-context and DRACH) and compared this data from DMSO- and STM2457-treated NHDF datasets.

Figure 2.

We observed that 75-90% of A basecalls (in all datasets) have probability scores above 0.8 (Fig. SIA). Since all A bases are considered during all-context m⁶A calling, the fraction of m⁶A sites is substantially smaller when compared to DRACH-context m^bA sites (Fig.2A,B, Fig.SIB,C). The modification probability distributions are effectively identical for DMSO-treated, STM2457-treated and IVT datasets below 0.98, while the profiles differ substantially at probabilities > 0.98. This suggests that sites detected below 0.98 are likely false positives, indicating that only high-probability m⁶A calls can be considered reliable. Here, STM2457 treatment shows a marked reduction in the number of reported m⁶A sites in both all-context and DRACH-only datasets when compared to the DMSO treatment (Fig. 2A,B).

In addition to m⁶A, we examined modification probability distributions for ψ, inosine and m⁵C, each of which can be basecalled with Dorado v0.9.0 (Fig. S2). For ψ and I, the modification probability distributions for the STM2457 and DMSO datasets only diverge from the IVT dataset at > 0.99 modification probability, thus providing a relevant threshold for filtering out false positives (Fig. S2A,B). By contrast, m⁵C modification probability distributions were similar between all three datasets and no filter could be adequately established (Fig. S2C). While this likely reflects a low rate of m⁵C installation, it also indicates that the m⁵C model utilized by Dorado requires further refinement.

Modification probability distributions have similar profiles fortranscriptome-aligned NHDF data (Fig. S3) as well as genome- and transcriptome-aligned data from HD10.6 cells (Fig. S4, Fig. S5). Of note, the fraction of high-probability inosine sites is increased in STM2457-treated cells, compared to DMSO-treated ones in genome-aligned data (Fig. S2B, Fig.S4D), which points to increased numbers of retained introns in the METTL3-inhibited sample. Accordingly, this effect is absent in transcriptome-aligned data (Fig. S3D, Fig. S5D).

While IVT transcriptomes serve as excellent negative controls, we also considered whether the RNA calibration strand (RCS) that is added during DRS library preparation could serve a similar function. The RCS is a synthetic Enolase II RNA derived by IVT¹³ that, similarly to the IVT transcriptomes, are devoid of RNA modifications. Despite the reduced diversity of sequence in the RCS data, the overall results were comparable to that from the IVT transcriptome data (Fig. S6A,B, Fig. 2A,B), demonstrating that the RCS can be used as a “built-in” negative control for modification detection with DRS. Finally, we determined that basecall and modification probability distributions remain generally consistent between Dorado versions although more m⁶A and ψ sites are detected with high modification probability in Dorado versions v0.8.0/v0.9.0 compared to v0.7.0 (Fig. S6C-E).

Filtering by stoichiometry further eliminates potential false positives

We next examined the stoichiometry distributions for all m⁶A sites in both unfiltered and filtered NHDF datasets (Fig. S7). Here, unfiltered datasets showed that the vast majority of m⁶A sites, 92% for all-context and 75% for DRACH-context, have stoichiometry values < 10% with less than 1% sites showing high (> 50%) stoichiometry. Application of our filtering strategy removed over 70% of the sites (Fig. 2E), regardless of the context, almost all of which showed stoichiometries of < 10%, without significantly changing the numbers of high stoichiometry sites detected (Fig. S7). A similar analysis of the IVT transcriptome revealed large numbers of false positive m⁶A sites with stoichiometry values below 10%, most of which were eliminated by application of our filters (Fig. S8A,B). However, the remaining presence of low stoichiometry false-positive sites indicates that an additional layer of filtering to exclude sites with < 10% stoichiometry is beneficial. An additional argument in support of our modification probability filtering strategy is the observation that treatment with the METTL3 specific inhibitor STM2457 dramatically reduced m^ŕA installation at sites with > 10% stoichiometry in untreated datasets in a context-independent manner (Fig. 2C,D, Supplementary Table 2). Similar results were observed for the transcriptome aligned NHDF and HD10.6 m⁶A stoichiometry distributions as well (Fig. S9, Fig. S10, Supplementary Table 2).

To determine whether our filtering strategy could be applied to other RNA modifications detected by Dorado, we generated stoichiometry distributions for ψ and inosine in unfiltered and filtered (> 0.99 modification probability) human genome-aligned DMSO- and STM2457-treated NHDFs and HDlO.6s. As with m⁶A (Fig. S7), 92% of detected ψ and 96% of inosine have less than 10% stoichiometry (Fig. S11, Fig. S12). A recent study that explores ψ in human mRNAs using 2-bromoacrylamide-assisted cyclization sequencing (BAGS) corroborates our results, showing that 82% of ψ modifications are detected at low levels (below 20%) (Xu et al. 2024).

Above 10% stoichiometry, in our data, there are more inosine sites in STM2457-treated cells than DMSO-treated cells, reflected also in the stoichiometry distributions (Fig. S11, Fig. S12, Supplementary Table 2). ψ and inosine stoichiometry distributions in our IVT data (Fig.S8C,D) further validates our filtering strategy of > 0.99 modification probability, as it substantially reduces the number of detected false-positive sites in the unfiltered data even at higher (> 50%) stoichiometries in case of inosine (Fig.S8D).

We have thus demonstrated that setting thresholds for basecall modification probabilities provides a robust method for identifying and filtering false-positive modification calls. For the purposes of this

study, we further evaluated the effect of placing the modification threshold on the modification data and restricted the downstream analyses to m⁶A sites with probability scores > 0.98.

RNA modification detection rates differ between Dorado versions

Dorado basecalling relies on models generated using Remora (github.com/nanoporetech/remora) that are frequently updated. This raises an important question as to whether these updates produce notably different outputs that could impact biological interpretation. To address this, we utilized multiple Dorado versions and Remora models to perform modification calling of our NHDF dataset. We restricted our analyses to m⁶A and ψ as these are the only modifications that could be detected in previous versions. Here, we observe that m⁶A and ψ stoichiometry distributions are generally consistent between different Dorado versions (Fig. S13), albeit with greater numbers of low stoichiometry (< 10%) sites detected in older Dorado versions. There is a large overlap of the detected > 1%-stoichiometry m⁶A (DRACH and all-context) and ψ sites between versions, although there are a lot of unique all-context m⁶A and ψ sites identified with just v0.7.0 and not v0.8.0-v0.9.0, that are reduced after applying our modification probability filter (Fig. S14A,B). ψ sites uniquely detected with different versions are found in slightly different sequence contexts (Fig. S14C,D). Thus, while RNA modification detection with Dorado is sensitive to changes in the underlying models/versions, this can be effectively controlled through the application of filters that only retain sites with high-confidence modification probability calls.

GLORI exhibits reduced specificity but converges with Dorado predictions at high stoichiometry sites

To assess the agreement of DRS Dorado-derived predictions with an orthogonal m⁶A detection method, we conducted GLORI on the same NHDF RNA preparations. To ensure a fair comparison of m⁶A sites across the two methods, we applied the following filtering strategies: For Dorado (DRACH model), this included thresholds on modification probability (A > 0.8, m⁶A > 0.98) and read coverage (> 20). For GLORI, we used filters based on coverage (> 20) and statistical significance (p-value < 0.05). Next, we extracted m⁶A sites that were present in both DMSO and STM2457 treated samples and removed all sites with < 10% stoichiometry in DMSO (Fig. 3A). This stepwise strategy reduced false positives and enabled a fair comparison of the high confidence m⁶A sites identified by both approaches.

As shown in the violin plots, both the all-context and DRACH-context Dorado models identify similar numbers of m⁶A sites (81,012 (all-context model) and 72,372 (DRACH model)), the stoichiometries of which were strongly reduced upon STM2457 treatment (Fig. 3B). When we analyzed the motifs of the m⁶A sites detected with the all-context model, we observed nearly 9,000 sites outside of DRACH motifs.

Importantly, the total number of m^sA sites within DRACH motifs was nearly identical to the number of m⁶A sites reported using the DRACH model (Fig.3B). Unexpectedly, GLORI detects a substantially higher number of sites when using all motifs (140,934) compared to DRACH motifs alone (27,926) (Fig. 3C). The fact that most GLORI-identified sites fall outside DRACH motifs suggests these are false positives. This interpretation is reinforced by the observation that the overlap between Dorado and GLORI datasets remains virtually unchanged when comparing GLORI “all motif” and “DRACH motifs” (n=l3,6l4) site sets, underscoring that the additional non-DRACH sites in GLORI (n= 113,008) are largely spurious (Fig. S15A).

We obtained similar results when comparing m⁶A sites obtained with Dorado and GLORI from HD10.6 cells. Overall, we detected a similar number of m⁶A sites in HD10.6 cells using the all-context and DRACH-specific Dorado models (77,728 and 67,222, respectively) (Fig. S15B). By contrast, we detected 174,282 vs 27,965 m⁶A sites using GLORI without or with DRACH filtering, respectively (Fig. S15C).

To evaluate how stoichiometry thresholds affect the agreement between GLORI and Dorado, we compared the overlap of high-confidence GLORI-identified sites with Dorado m⁶A (DRACH context) calls across increasing m⁶A stoichiometry cutoffs (Fig. 3D). At the lowest stoichiometry threshold applied (> 10%), 13,614 sites overlapped, representing less than half (48.8%) of GLORI sites. However, as the stringency increased, the relative overlap increased sharply: at > 20%, nearly 70% of GLORI sites overlapped (11,536 sites), at > 50% almost 90% overlapped (6,307 sites), and at the highest cutoff (> 80%), over 93% of GLORI sites (3,119 sites) were also detected by Dorado. These results demonstrate that while GLORI detects a broad set of potential m⁶A sites, only the sites with high stoichiometries overlap with Dorado, indicating high signal to noise ratios in the GLORI datasets.

Figure 3.

GLORI relies on the high-efficiency conversion of adenosine to inosine with overall gene A to I conversion rates needing to exceed 95-98% (Xu et al. 2024; Xie et al. 2025). Similarly, we observed the conversion rates for our NHDF and HD10.6 datasets to be 95% (Supplementary Table 3), but we found that small variations in gene conversion rates can strongly influence the sensitivity and specificity of GLORI, as only high-stoichiometry m⁶A sites overlapped with DRS-derived predictions (Fig. 3D).

DRS coupled with Dorado outperforms GLORI in isoform-specific m⁶A detection

Next, we wanted to highlight the advantage of long-read sequencing in terms of mapping reads with isoform-precision to the human transcriptome, which is not possible using GLORI due to its reliance on short-read sequencing. We first confirmed that mapping Dorado-basecalled sites to the human transcriptome results in the detection of similar numbers of m⁶A sites in both NHDFs and HD10.6 cells (Fig. S15D,E). To assess if DRS coupled with Dorado identifies distinct m⁶A profiles on different transcript isoforms, we analyzed the SPEN gene locus which was previously shown to contain more than 100 m⁶A sites (Liu et al. 2023). On a genomic level, we detected 42 m⁶A sites using DRS, while GLORI only detected 5 m⁶A sites (Fig. 3E). Out of the 7 total isoforms of SPEN, SPEN-201 and SPEN-202 were detectable in DRS-derived data, considering a minimal coverage of 20 (Supplementary Table 4). When plotting the m⁶A sites detected by transcriptome-aligned DRS data, we found that the SPEN-201 isoform exhibited more m⁶A sites, albeit with a lower median stoichiometry compared to the SPEN-202 isoform (Fig. 3F). This example highlights that DRS is superior compared with GLORI for isoform-level m^ŕA quantification.

Dorado profiling reveals cell type differences in m⁶A distribution and stoichiometry

To assess the reproducibility of Dorado across DRS datasets derived from different cell types, we compared the m⁶A methylomes of NHDFs and HD10.6 cells. To ensure a fair comparison, we extracted genes that were expressed in both cell types. Next, we compared the total numbers of high confidence m⁶A sites in NHDFs and HD10.6 cells within the genes expressed in both cell types and found that 87% of HDlO.6-derived m⁶A sites overlap with those detected in NHDFs (Fig. 4A). Site-specific m⁶A stoichiometries from the shared gene set correlated strongly across cell types (Pearson r = 0.93), although we observed a global shift in the HD10.6 dataset, where m⁶A stoichiometries were generally higher than in NHDFs (Fig. 4B). Plotting the m⁶A stoichiometry distributions of the 38,217 m⁶A sites shared between both cell types showed a strong increase in the median m⁶A stoichiometry of HD10.6 cells (70.2%), compared to NHDFs (54%) (Fig. 4C).

Figure 4.

We next aimed to conduct a comparative motif analysis of m⁶A sites detected in NHDFs and HD10.6 cells. Sequence logos confirmed strong enrichment of the canonical DRACH consensus motif for both cell types in modification probability-filtered data, with the GGACU pentamer being the most prominent (Fig. 4D,E, Fig.Sl6). Quantification of 5-mer distributions revealed highly similar motif preferences across datasets. While the global motif architecture of m⁶A deposition is conserved between NHDFs and HD10.6 cells, we observed a trend toward higher stoichiometries in the HD10.6 dataset, particularly within the most abundant DRACH motifs, suggesting enhanced modification efficiency in the neuron-like context (Fig. 4F, G).

Metagene analysis showed, as known previously, that m^sA is enriched near the stop codons and in the 3’ UTRs (Fig. S17A). Several studies have examined the relationship between exon architecture and m⁶A distribution on mRNAs^1,2,3, and so we also examined m⁶A sites in monoexonic transcripts separately and found that monoexonic transcripts are more broadly m⁶A-modified across entire CDS (Fig. S17A). m⁶A (both all-context and DRACH-context) distribution profiles look less typical for all-stoichiometry or< 10% stoichiometry sites in both NHDF and HD10.6 datasets (Fig. S17B+C). For all-context m⁶A (Fig. S17B), modification filter also proves to be deciding for the distribution profile, together with > 10% stoichiometry.

To connect shared m⁶A sites with their functional context and illustrate their distribution on individual transcripts, we conducted a KEGG pathway analysis with the genes that contain shared m⁶A sites between NHDFs and HD10.6 cells. We observed a significant enrichment in pathways related to neurodegeneration, infection, and intracellular signaling, including axon guidance and mTOR signaling (Fig. 4H). To further illustrate these findings, we examined representative genes contained in one or more of these pathways. WNTB5, MAPK1, AKT1, NCK1 and NRP1 all showed robust methylation in both NHDFs and HD10.6 cells, with consistently higher stoichiometries observed in the neuronal-like context. In contrast, ACTG1, included as a housekeeping control not represented in the enriched pathways, displayed stable and comparatively uniform methylation across both cell types (Fig. 4I).

Together, these results demonstrate that while the global architecture of m⁶A deposition is largely conserved between primary fibroblasts and neuronal-like cells, but unexpectedly the neuronal context is characterized by a global shift toward higher m^bA stoichiometries. This difference is readily detectable in single genes. These findings suggest that cell type-specific regulation of m⁶A may fine-tune core cellular pathways relevant to that cell type.

Discussion

High rates of false positives remain a major problem in epitranscriptomic mapping tools, even with the advent of modification-aware basecalling models (Cruciani et al. 2025). Comparing the outputs of orthogonal methods performed on the same input material, combined with use of IVT RNA and/or specific methyltransferase inhibition, offers an effective mitigation strategy (Cruciani et al. 2025; Diensthuber and Novoa 2025). The results presented here establish a robust framework that minimizes false positives and enables isoform-specific and quantitative m⁶A mapping. Through a direct comparison with GLORI, we showed that combining DRS datasets with Dorado modification-aware basecalling provides greater scalability, isoform resolution, and reproducibility across different cell types. Further, the simplicity of DRS library preparation enhances consistency whereas GLORI is constrained by a laborious and technically challenging protocol that is extremely sensitive to adenosine conversion rates (Shen et al. 2024). Despite these limitations, we determined that the majority of m⁶A sites in our filtered DRS datasets could be confirmed by GLORI. By expanding our analysis across multiple human cell types, we further determine that the sites of m⁶A installation are conserved for the most part but that the stoichiometries show cell-type specific differences.

Crucially, the inclusion of transcriptome-wide IVT datasets and the use of a selective METTL3 inhibitor (STM2457) provides a robust approach to developing filters that dramatically decrease false-positive m⁶A detection. The impact of both approaches was readily visualized by plotting modification probability distributions, enabling our determination that retaining only m⁶A calls with confidence scores of 0.98 or greater retained almost all STM2457-sensitive m⁶A sites while eliminating the large numbers of false-positive m⁶A sites reported for the IVT datasets. This approach is supported by a recent study which highlighted that Dorado introduces large numbers of false positives outside canonical DRACH contexts (Xie et al. 2025). Our data confirms that in the unfiltered datasets, the number of false positives is higher in the all-context model compared to the DRACH model. This may reflect the degenerate nature of the canonical DRACH consensus and understudied differences in sequence preference between different cell types and organisms. Retaining only sites at which m⁶A stoichiometries exceeded 10% further reduced false-positive rates while also focusing the results on identification of sites that are more likely to be biologically significant. The combined effect of these filters was to provide a strong agreement between the m⁶A all-context and m⁶A DRACH model in the detection of high confidence m⁶A sites, as observed by filtering for DRACH motifs from m⁶A sites detected with the all-context model (Fig.3B).

Prior studies have emphasized the need for isoform-resolved m⁶A quantification by demonstrating widespread heterogeneity in m⁶A levels between genes and their isoforms (Guo et al. 2025; Xie et al. 2025). Here, we showed that aligning DRS datasets to the human transcriptome reveals differential m^ŕA installation patterns between isoforms, underscoring the importance of accurate, isoform-specific m⁶A quantification at single-molecule resolution and highlighting the advantage of DRS against short-read sequencing approaches such as GLORI (Fig. 3E,F).

Recently, the ability of Dorado and m6Anet (Hendra et al. 2022) to identify m^sA in DRS datasets from 293T and HeLa cells were systematically evaluated using IVT RNA and comparisons to published datasets obtained by GLORI and eTAM-seq (Zou et al. 2025). In this case the authors found that Dorado (v0.7.2) outperforms m6Anet but also advocated for the importance of applying filters to reduce false positive calls. Likewise, we found that IVT RNA served as a valuable negative control to detect false positive calls. Furthermore, using the RNA calibration strand (RCS) as a negative control provided a convenient and simpler alternative to IVT RNA, yielding comparable probability distributions and thereby enabling a built-in control for routine experiments (Viehweger et al. 2019). While m⁶A detection was generally reproducible across Dorado versions, model-specific differences are present, i.e. v0.7.0 produces substantially more low-confidence all-context m⁶A and ψ calls than v0.8.0-0.9.0, highlighting that the models remain under active development and that benchmarking should be repeated when switching to a new model. Recent work further highlights that continued improvements in both nanopore sequencing chemistry and basecalling models significantly enhance RNA modification detection accuracy, emphasizing the importance of revisiting benchmarking and filtering strategies as models evolve (Hewel et al. 2025).

Focusing on other modifications that are detected by Dorado, we observe that inosine and ψ detection appears robust when applying similar filtering strategies as for m⁶A, the latter supported by other studies (Esfahani et al. 2025; Verstraten et al. 2025). For inosine, we observed an enrichment of high-probability (> 0.99) inosine sites in genome-aligned data from STM2457-treated samples. Further studies are needed to determine whether this reflects increased intron retention or misclassification of modifications by the current inosine model. While DRS is the only approach that allows for robust isoform-specific m⁶A mapping, comparisons with orthogonal approaches remain essential to validate findings. When comparing Dorado and GLORI using the same input RNA, we unexpectedly observed that only 50% of m⁶A sites with stoichiometries above 10% that were identified by Dorado were also present in the GLORI data although increasing the stoichiometries led to increased overlap (Fig.3D). We attribute this to a sub-optimal gene-wide adenosine conversion rate of 95% during GLORI, compared to the 98-99% rate reported in the original study (Shen et al. 2024). Here, unconverted adenosines result in high signal-to-noise ratios which result in lower sensitivity, explaining why the overlap with Dorado increases at higher-stoichiometry sites. This is also reflected in the relatively high amount of false positive sites detected outside of DRACH motifs in the GLORI datasets (Fig. 3C). The sensitivity of GLORI to small fluctuations thus limits its utility in comparison to DRS. Whether this can be resolved with improvements to the GLORI methodology (Sun et al. 2025) requires further study.

Prior analyses of m^sA stoichiometry distributions across immortalized human cell lines, including cancer cell lines, have revealed broadly similar m⁶A stoichiometry distributions (Guo et al. 2025). By contrast, our analysis of primary fibroblasts (NHDFs) and neuronal cells (differentiated HD10.6 cells) revealed an overall increase in m⁶A stoichiometries in neurons compared to fibroblasts, especially in genes associated with neuronal function. While further studies are warranted to establish the underlying mechanisms and functional consequences of these differences, this aligns with the importance of m⁶A in neuronal function (Yu et al. 2021). This includes m⁶A maps from synaptic and nuclear RNA isolated from mouse dorsal hippocampus which showed that synapse-specific RNAs showed generally higher methylation levels (Sun et al. 2025). Intriguingly, neuronal regulation pathways utilize differentially spliced isoforms (Su et al. 2018) and a recent study revealed isoform-specific differences in m⁶A sites and stoichiometries in the human brain that are linked to m⁶A machinery expression and cell type context (Gleeson et al. 2025).

In summary, using matched RNA samples from different experimental conditions and from different human cell types analyzed with both Dorado and GLORI orthogonal methodologies, we establish a robust pipeline for the precise and quantitative transcriptome-wide detection of m⁶A in DRS datasets. Future studies should aim at extending this approach using orthogonal sequencing approaches to benchmark Dorado-based detection of ψ, m⁵C and inosine, as well as the newly released models for 2’O-methylation.

Methods

Cell culturing and STM2457 treatment

NHDFs (Lonza) were cultured in Dulbecco’s Modified Eagle Medium (DMEM, Gibco), supplemented with 5% heat-inactivated FBS (HI-FBS, Gibco) and 1% penicillin-streptomycin (PS, Lonza) under standard culture conditions. The HD10.6 cell line was originally derived from fetal human dorsal root ganglia (DRG) and immortalized by overexpression of the v-myc oncogene under tetracycline regulation. Upon addition of tetracycline and neurotrophins, these cells can be differentiated into neuron-like cells that show typical neuronal morphology and nociceptive properties(Raymon et al. 1999). For our experiments, we used HD10.6 cells that have been matured for 10 days. For maintenance of undifferentiated HD10.6 cells (a gift from Anna Cliffe, University of Virginia School of Medicine), T-75 flasks were coated with 17 µg/mL fibronectin in PBS for 15 min at 37°C before adding the cell suspension. Cells were grown in Advanced DMEM/F12 (Gibco), supplemented with lx Neurocult SMI (STEMCELLTechnologies), 2 mM L-glutamine (Gibco), lx Primocin (Invivogen) and 10 ng/mL prostaglandin El (Sigma). 0.5 ng/mL bFGF (PeproTech) were added freshly to each media change. For maturation, HD10.6 cells were seeded at a density of 25.000 cells/cm² in 6-well plates which were pre-coated overnight with 50 µg/mL poly-L-ornithine hydrobromide (Sigma) in 0.5M borate buffer pH 8.5 (Boston Bioproducts), and 1 µg/mL fibronectin (overnight). After cell attachment overnight, media was changed to Neurobasal media (Gibco), supplemented with lx Neurocult SMI, 2mM L-glutamine, lx Primocin, 1 µg/mL doxycycline, 50 ng/mL NGF (Alomone Labs), and 25 ng/mL CNTF, GDNF and NT-3 (PeproTech). Media was changed every three days until HD10.6 cells showed complete neuronal morphology after 12 days. NHDFs and matured HD10.6 cells were treated with 30 µM STM2457 (Sigma) or DMSO for 48h, harvested in TRIzol (Invitrogen), and stored at -80°C until RNA isolation.

RNA isolation

RNA was extracted by the TRIzol method according to the manufacturer’s instructions. Briefly, 0.2 mL chloroform per 1 mL TRIzol were added and mixed by vortexing. After centrifugation for 30 min at l2,OOOxg at 4°C, the aqueous phase was transferred to isopropanol (0.5 mL per 1 mLTRIzol), mixed and incubated for 10 min at room temperature. Subsequently, the samples were centrifuged for 10 min at l2,OOOxg at 4°C and the RNA pellets were washed twice in 75% fresh ethanol. After air-drying for 5-10 min, RNA was dissolved in nuclease-free H₂O and stored at -80°C until further processing. One third of each sample was used for DRS and the remaining same-input material was processed for GLORI. For both DRS and GLORI, poly(A) RNA was selected using the Dynabeads™ mRNA Purification Kit (Invitrogen), according to the manufacturer’s instructions.

GLORI-Sequencing

GLORI reactions were largely performed following the published protocol for RNA protection, deamination and deprotection with several minor adjustments (Shen et al. 2024). Briefly, poly(A) RNA was fragmented at 94°C for 2 min (NEB E6150S) and subjected to glyoxal protection, sodium-nitrite-mediated deamination (750 mM NaNO₂ at 16°C for 8h and 4°C for 6h), and deprotection in triethylammonium acetate-formamide buffer at 95°C for 10 min. Library construction was conducted using a ligation-based strategy as described by Liu and colleagues incorporating 5’adapters with 11-nt UMIs(Liu etal. 2025).

Indexed libraries were amplified with NEBNext Dual Primers Set 2 (E778OS), purified using AMPure XP beads, and validated on an Agilent TapeStation (D1000 High Sensitivity). Sequencing was performed on an Illumina NovaSeq 6000 platform (150 bp paired-end mode; >80 million reads per sample).

For data processing, we employed the same computational pipeline described in Liu et al (Liu et al. 2025), but implemented within a Snakemake-based framework that automates the complete analysis starting from raw fastq files (Koster and Rahmann 2012). Details on software versions, workflow availability, and repository links are provided in the Data and code availability section below.

Generation of a transcriptome-wide IVT RNA library

We sequenced the IVT RNA from both NHDFs and HD10.6 cells along with the native poly(A) RNA to evaluate potential false-positive predictions by Dorado. This was achieved by converting the poly(A) RNA to cDNA using a poly(A) primer and a T7-incorporating strand switch primer and using the cDNA as a template for IVT synthesis.

Nanopore Sequencing of poly(A) RNA

250-300 ng of isolated poly(A) RNA and 25 ng RNA calibration strand (RCS) were used as input for the standard DRS SQK-RNA004 protocol. Resulting libraries were loaded onto individual PromethlON RNA flow cells before sequencing on the PromethlON 2 Solo instrument for 48h.

Nanopore basecalling

Basecalling was performed using Dorado v0.5.0, v0.6.0, v0.7.0, v0.8.0 and v0.9.0 (ONT) in super accuracy mode with modification detection enabled for m⁶A (DRACH context) in v.0.5.0, 0.6.0, 0.8.0 and 0.9.0, m⁶A (all context) in v0.7.0, 0.8.0 and 0.9.0, ψ in v0.7.0, v0.8.0 and v0.9.0, and inosine and m⁵C in v.0.9.0. Each dataset was processed with adapter trimming. Resulting unaligned BAM files were converted to FASTQ files containing modification tags for every read in the FASTQ read headers, using SAMtools vl.18 (Li et al. 2009) [samtools fastq -TMM,ML],

Genome and transcriptome level alignments

For genome-level alignment, reads were aligned against a combined human genome reference (GRCh38.pl4), and the RCS gene sequence, using minimap2 v2.26 [minimap2 - ax splice -kl4 -y -uf — secondary=no]. For transcriptome-level alignment, the reads were aligned against the human transcriptome (v47), obtained via GENCODE (Mudge et al. 2025) [minimap2 -t 8 -ax map-ont -y -L -p 0.99], Using SAMtools vl.18, resulting genome- and transcriptome-aligned SAM files were converted to BAM files retaining only primary alignments [samtools view -b -F2308 and samtools view -b -F2324, respectively]. BAM files were sorted [samtools sort] and indexed [samtools index].

Processing of RNA modification calls using Modkit

Probability distributions of the detected sites in human genome- and transcriptome-aligned reads were generated using sample-probs subcommand from Modkit v0.4.1 [sample-probs --hist --only-mapped --percentiles 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.85,0.9 -f 0.1]. For RCS reads -f value was increased to 0.25. The resulting .tsv files were imported into RStudio (v. 2024.12.1+563) and plotted using a custom R(v. 4.4.2) script (Dorado_mod_probability_distributions.R), with the following packages: dplyr (vl.1.4), ggplot2 (v3.5.1), reshape2 (vl.4.4), patchwork (vl.3.0).

To generate extended bedMethyl files with m⁶A-modified and unmodified base counts from the reads mapping against the potential modification sites, we used Modkit pileup subcommand with the thresholds [-filter-threshold A:0.8 -filter-threshold T:0.8 -filter-threshold C:0.8 -filter-threshold G:0.8 - motif A 0 —ref$REF] with and without modification threshold filter [-mod-thresholds a:0.98] to generate our “filtered” and “unfiltered” datasets, respectively. Canonical nucleotide calls (A, C, G, U) were retained with > 80% probability scores for all modification datasets. Modification thresholds were set at 98% for m⁶A and 99% for all the other modifications, —motif was adjusted to the base for which modification status was examined, which further restricted the analysis to only the positions that matched canonical nucleotides in the reference genome. The resulting sites were further filtered according to the read number per site (> 20 reads) for all the downstream analysis.

Gene Ontology analysis

GO analysis was performed using DAVID (Huang da et al. 2009). Genes harboring m⁶A sites were filtered for genes that were generally present in both datasets. Next, they were divided into NHDF-specific, HDlO.6-specific, and shared categories, and enrichment was assessed for biological processes. GO terms were considered significant at FDR < 0.05 with a minimum of 20 genes per GO term. For plotting, the top 20 GO terms that contained the highest gene numbers were extracted from the analysis of the shared gene sets.

Analysis and visualization of modification stoichiometries and distributions

Stoichiometry distribution histograms and metaplots were generated using Rstudio and R v4.4.2 with the following libraries: data.table, patchwork, scales, tidyr, dplyr, ggplot2⁴. Access to scripts used for analysis and visualization is detailed below.

Data availability

Raw fastq and pod5 datasets generated in this study are available via the European Nucleotide Archive (ENA) under the accession number PRJEB101380. R scripts used in the analysis and visualization of data presented herein are available in the https://github.com/DepledgeLab/HOLDEN repository. The Snakemake workflow implementing the GLORI analysis pipeline is available in the https://github.com/gp-micro/glori repository.

Supplementary figures

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Figure S8.

Figure S9.

Figure S10.

Figure S11.

Figure S12.

Figure S13.

Figure S14.

Figure S15.

Figure S16.

Figure S17.

Acknowledgements

DPD is supported by a German Centre for Infection Research (DZIF) Associate Professorship and the NIAID grants R01-AI170583 and R01-AI152543. DPD also receives funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC 2155, project number 390874280. EL is supported by the NIAID grant R01-AI170583 as well as the Hannover Biomedical Research School (HBRS) and the Center for Infection Biology (ZIB). DO is supported by the Walter Benjamin Programme of the DFG, project 554758329. ACW is supported by NIAID grants R01AI176335 and R01-AI170583.

Additional files

Supplementary Table S1.

Supplementary Table S2.

Supplementary Table S3.

Supplementary Table S4.

Additional information

Funding

German Center for Infection Research

NIAID (R01-AI170583)

NIAID (R01-AI152543)

NIAID (R01AI176335)

Deutsche Forschungsgemeinschaft (390874280)

Deutsche Forschungsgemeinschaft (554758329)