Workflow form A mapping in different human cell types using orthogonal sequencing approaches.

(A) Schematic of the experimental setup. RNA was isolated from NHDFs and differentiated HD10.6 cells treated with STM2457 or DMSO for 48h. Poly(A) RNA was subjected to DRS and GLORI from the same input material. IVT RNA was generated in parallel to assess potential false-positive modification calls. (B) Workflow for benchmarking Dorado basecalling. Reads were processed with Dorado versions v.0.5.0 through v0.9.0, aligned to the human genome or transcriptome, and analyzed with ONT’s Modkit tool (Rand 2024). Multiple filtering strategies were tested against IVT controls and METTL3 inhibition (STM2457). Resulting m⁶A sites were compared with GLORI data from the same input RNA. (C) Final benchmarked Dorado pipeline to detect high confidence m6A sites.

Modification probability of 0.98 and predicted m A/A proportion score of ≥ 10% are necessary cutoffs to identify m A sites accurately.

(A) (left) Dorado v0.9.0-generated modification probability distributions for all-context m6A from genome-aligned DMSO- and STM2457-treated NHDFs plotted against the distributions obtained from the IVT NHDF poly(A) RNA. DMSO probability score distribution is shown as a red solid line, STM2457 – red dashed line, IVT – grey solid line. (right) Bar plots showing the per-bin log2(sample/IVT) ratio of modification probability fractions for the DMSO- and STM2457-treated NHDFs conditions relative to the unmodified IVT reference. The horizontal navy dashed line marks the fold-enrichment threshold used to call enrichment. Bars are coloured by significance relative to the fold-enrichment threshold (red, log2 ≥ 1.5-fold; dark grey, not significant). The light grey ribbon indicates the 95% bootstrap confidence interval around the log2 ratio estimate at each bin. The vertical pink dashed line indicates the derived modification probability cutoff above which sites are considered genuinely modified. (B) Same as in (A), but for DRACH-context m6A sites. (C) Predicted m6A/A proportion (≥ 10%) distributions for all-context m6A from DMSO-(in red) or STM2457-treated (in pink) NHDF datasets (genome-aligned) are plotted using a 0.98 modification probability cutoff and a filter for ≥ 20 reads. (D) Same as in (C), but for DRACH-context m6A sites. (E) Overlap of m6A sites between all-context and DRACH models (A basecall accuracy ≥ 80%, coverage ≥ 20 reads, predicted m6A/A proportion ≥ 10%), comparing non-filtered to filtered (modification probability > 0.98) outputs.

Comparison of m A sites detected by Dorado and GLORI.

(A) Filtering strategy for Dorado vs GLORI to obtain a high-confidence list of m6A sites in NHDFs. (B) Predicted m6A/A proportion distributions are plotted from genome-aligned DMSO- or STM2457-treated NHDF datasets using the indicated Dorado models and the filtering strategy shown in (A). (C) Predicted m6A/A proportion distributions from GLORI performed on DMSO- or STM2457-treated NHDFs are plotted, using the filtering strategy shown in (A). In addition, m6A sites were filtered for DRACH motifs (right plot). (D) Overlap of high confidence m6A sites detected by Dorado vs GLORI. GLORI and Dorado datasets were filtered using increasing predicted m6A/A proportion cutoffs (from left to right). GLORI sites overlapping with Dorado-identified sites are shown in light grey, whereas sites unique to GLORI are shown in darker grey. (E) m6A sites in the SPEN gene detected by GLORI and Dorado. Boxplots (right) represent the predicted m6A/A proportions per m6A site shown on the left. (F) SPEN isoforms with m6A sites, determined by Dorado. Boxplots (right) show the distribution of predicted m6A/A proportions on the m6A sites on SPEN-201 and SPEN-202.

Comparative analysis of m⁶A modification patterns in NHDFs and HD10.6 cells using Dorado v0.9.0.

(A) Overlap of high-confidence m⁶A sites between NHDFs and HD10.6 cells. (B) Correlation of predicted m6A/A proportions at shared m6A sites (n=38,217). (C) Boxplots showing the predicted m6A/A proportion distributions of the 38,217 detected in both NHDFs and HD10.6 cells. (D,E) Motif analysis showing m6A (all-context and DRACH-context) motif logos, 5-mer distributions and predicted modification proportions in NHDFs. (F,G) Same as in (D,E) but for HD10.6 cells. (H) KEGG pathway analysis (DAVID Bioinformatics Resources) of shared m6A-containing genes. (I) Predicted m6A/A proportions of representative examples of individual genes from enriched pathways (H). ACTG1 serves as a housekeeping control.

Cell viability assays in STM2457-treated NHDFs and HD10.6 cells.

Cell viability was assessed in NHDF and HD10.6 cells after 48h of DMSO or 30 µM STM2457 treatment using CellTiter-Glo (Promega). Statistical significance was determined from 10 replicate wells per condition using an unpaired, two-tailed Student’s t-test.

Evaluation of Dorado by base and modification probability distributions.

(A) Dorado v0.9.0-generated basecall accuracy probability values for A (all-context and DRACH-context basecaller models) from genome-aligned DMSO- and STM2457-treated NHDF samples against IVT poly(A) NHDF RNA. A probability score distribution is shown as a black solid line for DMSO, STM2457 – black dashed line, IVT – grey solid line. (B) Dorado v0.9.0-generated modification probability distributions, visualized as histograms, for all-context m6A from genome-aligned DMSO- and STM2457-treated NHDFs and IVT NHDFs. IVT modification probability score distribution is shown in grey, DMSO – bright red, STM2457 – dark red. (C) Same as in (B), but for DRACH-context m6A sites.

IVT-based estimation of false-positive m A detection across modification probability thresholds and ranges in NHDFs.

(A) False-positive m6A ratios across increasing modification probability cutoffs. For each cutoff (0.5, 0.7, 0.9, 0.98, and 0.99), the number of m6A calls detected in the IVT control was divided by the corresponding number detected in DMSO- or STM2457-treated NHDF samples. (B) False-positive m6A ratios calculated for discrete modification probability ranges (0–0.5, 0.5–0.7, 0.7–0.9, 0.9–0.98, 0.98–0.99, and 0.99–1). In both panels, the upper row shows DRACH-context m6A calls and the lower row shows all-context m6A calls; left and right plots show ratios relative to DMSO- and STM2457-treated samples, respectively. Each point represents an independently generated subsample of 200,000 reads. Statistical significance was assessed using pairwise comparisons derived from linear mixed-effects models, with sample included as a random effect and Benjamini–Hochberg correction for multiple comparisons. ****P < 0.0001.

IVT-based estimation of false-positive m A detection across modification probability thresholds and ranges in HD10.6 cells.

(A) False-positive m6A ratios across increasing modification probability cutoffs. For each cutoff (0.5, 0.7, 0.9, 0.98, and 0.99), the number of m6A calls detected in the IVT control was divided by the corresponding number detected in DMSO- or STM2457-treated NHDF samples. (B) False-positive m6A ratios calculated for discrete modification probability ranges (0–0.5, 0.5–0.7, 0.7–0.9, 0.9–0.98, 0.98–0.99, and 0.99–1). In both panels, the upper row shows DRACH-context m6A calls and the lower row shows all-context m6A calls; left and right plots show ratios relative to DMSO- and STM2457-treated samples, respectively. Each point represents an independently generated subsample of 200,000 reads. Statistical significance was assessed using pairwise comparisons derived from linear mixed-effects models, with sample included as a random effect and Benjamini–Hochberg correction for multiple comparisons. ****P < 0.0001.

Modification probability cutoffs are necessary to identify Ψ and inosine sites more accurately.

(A) (left) Dorado v0.9.0-generated basecall accuracy probability values for U (left) and modification probability distributions for Ψ (right) from genome-aligned DMSO- and STM2457-treated NHDF samples are plotted against IVT poly(A) RNA. Ψ probabilities are visualized in blue solid and dashed lines for DMSO and STM2457, respectively and in grey solid line for IVT. (right) Bar plots showing the per-bin log2(sample/IVT) ratio of modification probability fractions for the DMSO- and STM2457-treated NHDFs conditions relative to the unmodified IVT reference. The horizontal navy dashed line marks the fold-enrichment threshold used to call enrichment. Bars are coloured by significance relative to the fold-enrichment threshold (blue, log2 ≥ 1.5-fold; dark grey, not significant). The light grey ribbon indicates the 95% bootstrap confidence interval around the log2 ratio estimate at each bin. The vertical pink dashed line indicates the derived modification probability cutoff above which sites are considered genuinely modified. (B) Same as in (A) but showing A base accuracy probabilities and inosine sites. Inosine probabilities are visualized in purple solid and dashed lines for DMSO and STM2457, respectively and in grey solid line for IVT. (C) Same as in (A) but showing C base accuracy probabilities and m5C sites. m5C probabilities are visualized in green solid and dashed lines for DMSO and STM2457, respectively and in grey solid line for IVT.

IVT-based estimation of false-positive Ψ, inosine, and m C detection across modification probability thresholds and ranges in NHDFs.

(A) False-positive modification ratios across increasing modification probability cutoffs (0.5, 0.7, 0.9, 0.98, and 0.99) for Ψ, inosine, and m5C. For each cutoff, the number of modification calls detected in the IVT control was divided by the corresponding number detected in DMSO-treated NHDF samples. (B) False-positive modification ratios calculated for discrete modification probability ranges (0–0.5, 0.5–0.7, 0.7–0.9, 0.9–0.98, 0.98–0.99, and 0.99–1). From left to right, plots show results for Ψ, inosine, and m5C. Each point represents an independently generated subsample of 200,000 reads. Statistical significance was assessed using pairwise comparisons derived from linear mixed-effects models, with sample included as a random effect and Benjamini–Hochberg correction for multiple comparisons. ****P < 0.0001; ns, not significant.

IVT-based estimation of false-positive Ψ, inosine, and m C detection across modification probability thresholds and ranges in HD10.6 cells.

(A) False-positive modification ratios across increasing modification probability cutoffs (0.5, 0.7, 0.9, 0.98, and 0.99) for Ψ, inosine, and m5C. For each cutoff, the number of modification calls detected in the IVT control was divided by the corresponding number detected in DMSO-treated HD10.6 samples. (B) False-positive modification ratios calculated for discrete modification probability ranges (0–0.5, 0.5–0.7, 0.7–0.9, 0.9–0.98, 0.98–0.99, and 0.99–1). From left to right, plots show results for Ψ, inosine, and m5C. Each point represents an independently generated subsample of 200,000 reads. Statistical significance was assessed using pairwise comparisons derived from linear mixed-effects models, with sample included as a random effect and Benjamini–Hochberg correction for multiple comparisons. ****P < 0.0001; ns, not significant.

Modification probability cutoffs are necessary to identify modification sites accurately, shown for transcriptome-aligned data.

Dorado v0.9.0-generated basecall accuracy probability values for A, U and C (left, A-E) and modification probability distributions for m6A (all-context), m6A (DRACH-context), Ψ, inosine and m5C (middle, A-E) in transcriptome-aligned DMSO-and STM2457-treated NHDFs, plotted against the values from IVT NHDF poly(A) RNA. (right) Bar plots showing the per-bin log2(sample/IVT) ratio of modification probability fractions for the DMSO- and STM2457-treated NHDFs conditions relative to the unmodified IVT reference. The horizontal navy dashed line marks the fold-enrichment threshold used to call enrichment. Bars are coloured by significance relative to the fold-enrichment threshold (blue, log2 ≥ 1.5 fold; dark grey, not significant). The light grey ribbon indicates the 95% bootstrap confidence interval around the log2 ratio estimate at each bin. The vertical pink dashed line indicates the derived modification probability cutoff above which sites are considered genuinely modified.

Modification probability distribution profiles in genome-aligned DMSO-and STM2457-treated HD10.6 cells.

Dorado v0.9.0-generated basecall accuracy probability values for A, U and C (left, A-E) and modification probability distributions for m6A (all-context), m6A (DRACH-context), Ψ, inosine and m5C (middle, A-E) in human genome-aligned DMSO-and STM2457-treated HD10.6 cells, plotted against the values from IVT HD10.6 poly(A) RNA. (right) Bar plots showing the per-bin log2(sample/IVT) ratio of modification probability fractions for the DMSO- and STM2457-treated HD10.6 conditions relative to the unmodified IVT reference. The horizontal navy dashed line marks the fold-enrichment threshold used to call enrichment. Bars are coloured by significance relative to the fold-enrichment threshold (blue, log2 ≥ 1.5-fold; dark grey, not significant). The light grey ribbon indicates the 95% bootstrap confidence interval around the log2 ratio estimate at each bin. The vertical pink dashed line indicates the derived modification probability cutoff above which sites are considered genuinely modified.

Modification probability distribution profiles in transcriptome-aligned DMSO-and STM2457-treated HD10.6 cells.

Dorado v0.9.0-generated basecall accuracy probability values for A, U and C (left, A-E) and modification probability distributions for m6A (all-context), m6A (DRACH-context), Ψ, inosine and m5C (middle, A-E) in human transcriptome-aligned DMSO-and STM2457-treated HD10.6 cells, plotted against the values from IVT HD10.6 poly(A) RNA. (right) Bar plots showing the per-bin log2(sample/IVT) ratio of modification probability fractions for the DMSO- and STM2457-treated HD10.6 conditions relative to the unmodified IVT reference. The horizontal navy dashed line marks the fold-enrichment threshold used to call enrichment. Bars are coloured by significance relative to the fold-enrichment threshold (blue, log2 ≥ 1.5-fold; dark grey, not significant). The light grey ribbon indicates the 95% bootstrap confidence interval around the log2 ratio estimate at each bin. The vertical pink dashed line indicates the derived modification probability cutoff above which sites are considered genuinely modified.

Modification site probabilities detected with different versions and models of Dorado and compared to RCS.

(A,B) DRACH-context and all-context m6A modification probability distributions from genome-aligned DMSO (solid red line)- and STM2457 (dashed red line)-treated NHDFs are plotted against RCS negative control (grey). (C) Modification probability distributions for all-context m6A (Dorado v0.7.0 (model v5.0.0), v0.8.0 (model v5.1.0) and v0.9.0 (model v5.1.0)), the bottom panel showing the distributions in 0.95-1.00 modification range. (D) Same as in (C), but for DRACH-context m6A (Dorado v0.5.0 (model v3.0.1), v0.6.0 (model v3.0.1), v0.8.0 (model v5.1.0) and v0.9.0 (model v5.1.0)). (E) Same as in (C), but for Ψ (Dorado v0.7.0 (model v5.0.0), v0.8.0 (model v5.1.0), v0.9.0 (model v5.1.0)). All the plots in C-E are generated from genome-aligned DMSO-treated NHDF datasets and are plotted against RCS negative controls.

Evaluation of Dorado by the distributions of predicted m A/A proportions.

(A) Distributions of predicted m6A/A proportions for all-context m6A from DMSO-(in red) or STM2457-treated (in pink) NHDF datasets (genome-aligned and filtered for ≥ 20 reads) are plotted, comparing 0.98 modification probability-unfiltered (left) and filtered (right) sites. (B) Same as in (A), but for DRACH-context m6A sites.

Distributions of predicted m A/A, Ψ/U and inosine/A proportions in IVT datasets.

(A) All-context predicted m6A/A proportion distributions (all proportions (top) and ≥ 10% proportions (bottom)) in (genome-aligned) IVT RNA from DMSO-treated NHDFs, comparing modification probability-unfiltered and filtered data. The bottom panel shows distributions only for ≥ 10% proportion sites. (B) As in (A), but for DRACH-context m6A. (C) As in (A), but for Ψ. (D) As in (A), but for inosine.

Evaluation of Dorado by predicted modification proportion distributions in transcriptome-aligned NHDFs.

(A) Distributions of predicted m6A/A proportions for all-context m6A from DMSO-(in red) or STM2457-treated (in pink) NHDF datasets (transcriptome-aligned and filtered for ≥ 20 reads) are plotted, assessing the 0.98 modification probability cutoff. The bottom panel shows distributions only for ≥ 10%-proportion sites. (B) Same as in (A), but for DRACH-context m6A sites.

Evaluation of Dorado by distributions of predicted m A/A proportions in transcriptome-aligned HD10.6 cells.

(A) Distributions of predicted m6A/A proportions for all-context m6A from DMSO-(in red) or STM2457-treated (in pink) HD10.6 datasets (transcriptome-aligned and filtered for ≥ 20 reads) are plotted, assessing the 0.98 modification probability cutoff. The bottom panel shows distributions only for ≥ 10%-proportion sites. (B) Same as in (A), but for DRACH-context m6A sites.

Distributions of predicted Ψ/U and inosine/A proportions in NHDFs.

(A) Distributions of all (left) and ≥ 10% (right) predicted Ψ/U proportions in the 0.99 modification probability-unfiltered and filtered data from genome-aligned DMSO-(darker color) and STM2457-treated (lighter color) NHDFs (filtered for ≥ 20 reads). (B) Same as in (A), but for inosine.

Distributions of predicted Ψ/U and inosine/A proportions in HD10.6 cells.

(A) Distributions of all and ≥ 10% predicted Ψ/U proportions in the 0.99 modification probability-unfiltered and filtered data from genome-aligned DMSO-(darker color) and STM2457-treated (lighter color) HD10.6 cells (filtered for ≥ 20 reads). (B) Same as in (A), but for inosine.

Relationship between predicted RNA modification proportions and nucleotide basecalling errors in NHDF and HD10.6 datasets.

(A) NHDF and (B) HD10.6 datasets showing, from top to bottom, DRACH-context m6A, all-context m6A, Ψ, inosine, and m5C sites. Prior to analysis, sites were filtered using modification probability thresholds of ≥0.98 for m6A and ≥0.99 for Ψ, inosine, and m5C, together with a minimum valid-read coverage of 20. For each modification type, scatterplots show the predicted modification proportion generated by Dorado and the percentage of overlapping reads basecalled as a nucleotide other than the corresponding reference nucleotide (A for m6A and inosine, T for Ψ, and C for m5C). Each point represents an individual modification site, with point color indicating density. Adjacent violin plots characterize the identity of alternative basecalls at sites for which ≥25% of all overlapping reads were basecalled as a non-reference nucleotide. For each site, the proportion of mis-basecalled reads assigned to each of the three alternative nucleotides was calculated, with individual points representing individual modification sites. Violin distributions therefore show the relative contribution of each alternative nucleotide among mis-basecalled reads at sites exhibiting high levels of basecalling mismatch.

Predicted m A/A and Ψ/U proportion distributions in different versions of Dorado.

Distributions of predicted modification proportions (all – left, ≥ 10% – right) in a line graph format for (A) all-context m6A detected in genome-aligned DMSO-treated NHDF dataset, with Dorado versions 0.7.0 (model v5.0.0), 0.8.0 (model v5.1.0) and 0.9.0 (model v5.1.0) colored in gradients of red, darkness increasing with the later versions. (B) DRACH-context m6A detected with Dorado versions 0.5.0 (model v3.0.1), 0.6.0 (model v3.0.1), 0.8.0 (model v5.1.0) and 0.9.0 (model v5.1.0), colored in gradients of red, darkness increasing with the later versions. (C) Ψ detected with Dorado versions 0.7.0 (model v5.0.0), 0.8.0 (model v5.1.0) and 0.9.0 (model v5.1.0), colored in gradients of blue, darkness increasing with the later versions. All the datasets were filtered for ≥ 20 reads.

Comparison of modification sites detected with different versions of Dorado.

Upset plots showing unique and shared detected modifications (filtered for ≥ 20 reads, ≥ 10% predicted modification proportion) in modification probability-unfiltered and filtered genome-aligned DMSO-treated NHDF data. The versions are ordered on the plot by the number of detected sites, in ascending order (on the left). Results are shown for the following models and versions: (A) all-context m6A detected with v.0.7.0 (model v5.0.0), v0.8.0 (model v5.1.0) and v0.9.0 (model v5.1.0), (B) DRACH-context m6A detected with v.0.5.0 (model v3.0.1), v0.6.0 (model v3.0.1), v0.8.0 (model v5.1.0) and v0.9.0 (model v5.1.0), (C) Ψ detected with v.0.7.0 (model v5.0.0), v0.8.0 (model v5.1.0) and v0.9.0 (model v5.1.0). (D) Nucleotide frequencies in the 5-base motif of Ψ sites uniquely for model v5.0.0 (basecalled with Dorado v0.7.0), common for models v5.0.0 (basecalled with v0.7.0) and v5.1.0 (basecalled with Dorado v0.8.0 and v0.9.0), and uniquely for model v5.1.0-detected (basecalled with v0.8.0 and v0.9.0) Ψ sites. Position 0 marks the modification site.

Comparison of m A sites detected by Dorado and GLORI-Sequencing in NHDFs and HD10.6s.

(A) Upset plot showing the number of m6A sites detected by GLORI and Dorado in NHDFs and their overlaps. Dorado sites represent high-confidence m6A sites obtained after filtering DRACH-model detected data. Similarly, only high confidence m6A sites detected by GLORI are shown but were additionally filtered for DRACH sites. (B) Distributions of predicted m6A/A proportions are plotted from genome-aligned DMSO- or STM2457-treated HD10.6 datasets using the indicated Dorado models. (C) Distributions of predicted m6A/A proportions from GLORI performed on DMSO- or STM2457-treated NHDFs are plotted. In addition, m6A sites were filtered for DRACH motifs (right plot). (D) Distributions of predicted m6A/A proportions are plotted from transcriptome-aligned DMSO- or STM2457-treated NHDF datasets using the indicated Dorado models. (E) Distributions of predicted m6A/A proportions are plotted from transcriptome-aligned DMSO- or STM2457-treated HD10.6 datasets using the indicated Dorado models.

Isoform-level m A prediction using the Dorado DRACH model in NHDF data.

(A) m6A sites in the KRTAP1-5 gene detected by GLORI and Dorado. Boxplots (right) represent the predicted m6A/A proportions per m6A site shown on the left. (B) m6A sites in the HSPA8 gene detected by GLORI and Dorado. Boxplots (right) represent the predicted m6A/A proportions per m6A site shown on the left. Below, HSPA8 isoforms with m6A sites, determined by Dorado. Boxplots (right) show the distribution of predicted m6A/A proportions on the m6A sites on HSPA8-201 and HSPA8-202. (C) Same as in B, but for CITED2.

Non-DRACH m A calls predicted by the Dorado DRACH model in NHDF data.

(A) Top 15 most frequently detected non-DRACH motifs, ranked by occurrence (left); bars are colored by the central dinucleotide, and percentages indicate the proportion of all non-DRACH calls. Bold labels mark the 4 most abundant motifs discussed in the text. Mean predicted m A/A proportions for each motif in DMSO- and STM2457-treated cells (right). Motifs conforming to DRAC but not DRACH are indicated with asterisks. (B) Representative genome browser view of an ACACU call located one to two nucleotides from a canonical DRACH motif, with the called and adjacent positions marked. (C) Read-level basecalls at a GGACG site, where the terminal position is called as G in 51% and A in 49% of reads, rendering the context non-DRACH and DRACH in approximately equal read fractions. (D) Representative AAAAA/AAAAC calls within a homopolymeric region. Neighboring DRACH motifs are indicated in blue.

Non-DRACH m6A calls predicted by the Dorado all-context model in NHDF data.

(A) Same as in Fig. S23A, but for m A sites predicted with the all-context Dorado model. (B) Representative locus showing an STM2457-sensitive GGAUU call without evidence of positional shift, sequence variation, or homopolymer context. (C) STM2457-insensitive ACAGA calls within MAT2A, shown relative to the METTL16 consensus (UACAGAGAA).

Metagene plots showing m A (all-context and DRACH-context) distribution profiles in NHDF and HD10.6 data filtered at different predicted m A/A proportions.

(A) Metagene plots showing the density of m6A (all-context – dark red, DRACH-context – bright red) sites across all annotated poly(A) transcripts (upper plot) and monoexonic poly(A) transcripts (lower plot) from genome-aligned DMSO-treated NHDF datasets filtered with 0.98 modification probability, ≥ 20 coverage, and ≥ 10% predicted m6A/A proportion. (B) The density of all-context m6A sites (0.98 modification probability-unfiltered – left, filtered – right) is plotted across all annotated poly(A) transcripts of NHDF (dark red) and HD10.6 (bright red) datasets, filtered with different predicted m6A/A proportions (descending: all proportions, ≥ 50%, ≥ 10% and < 10%. All datasets were filtered additionally for ≥ 20 coverage. (C) Same as (A), but for DRACH-context m6A sites.