Comprehensive mapping of RNA modification dynamics and crosstalk via deep learning and nanopore direct RNA-sequencing

Domain: chemistry
Published: 2026-01-14
Editorial outcome: minor_revision

Abstract
Despite the extensive studies of individual RNA modifications, the lack of methods to detect multiple modification types simultaneously has left the global epitranscriptomic landscape and its underlying crosstalk largely unexplored. Here, we present ORCA (Omni-RNA modification Characterization and Annotation), a deep learning framework that enables comprehensive mapping of RNA modification landscape using nanopore direct RNA sequencing. ORCA employs domain adversarial learning to detect and quantify a wide range of modifications by leveraging mixed stoichiometry-driven signal and sequence variability between modified and unmodified nucleotides. It also incorporates a transfer learning module for accurate annotation of modification types with minimal prior knowledge. Applying ORCA to multiple human cell lines reveals widespread, isoform-specific modification patterns, as well as intricate cooperative and competitive interactions among neighboring modification sites. This approach substantially expands the repertoire of known RNA modification sites and elucidates their spatial organization, revealing the emerging roles of RNA modifications in splicing regulation. ORCA thus provides an unbiased and generalizable framework for decoding RNA modification dynamics and their regulatory complexity across diverse biological contexts.

Introduction
RNA modifications represent a complex and dynamic layer of post-transcriptional regulation, with over 170 distinct chemical marks that regulate RNA stability1, splicing2,3, translation4,5, and subcellular localization6. While the functions of individual modification such as N6-methyladenosine(m6A), pseudouridine(Ψ), and 5-methylcytosine (m5C) have been extensively studied7,8, the combinatorial effects and crosstalk among different RNA modifications remain largely unexplored. Recent studies have revealed the coordinated roles for m6A and Ψ in modulating translation9, and a synergistic co-occurrence of m6A and m5C in plants under salt stress10, highlighting the emerging role of interaction between different modifications. Deciphering these interactions is crucial for understanding the multilayered regulatory mechanisms governing RNA biogenesis and function. However, progress in this area has been hindered by the lack of transcriptome-wide tools capable of simultaneously detecting and analyzing diverse RNA modifications and their interactions, limiting our ability to decode the full regulatory potential of the epitranscriptome.

Recent Illumina-based approaches using immunoprecipitation11–13 or chemical treatment14–17 have enabled transcriptome-wide profiling of individual RNA modification types, but are unable to capture the global epitranscriptomic landscape simultaneously8. Nanopore direct RNA sequencing (DRS) overcomes this limitation by directly sequencing native RNA molecules and recording ionic current signals that reflect each nucleotide’s chemical structure18,19. These inherent signals produce distinct ionic signal profiles and basecalling differences between modified and unmodified bases, encoding rich information about various RNA modifications within single molecules20. However, most existing DRS-based tools are either trained on in vitro synthesized datasets that are restricted to a few well-characterized modifications such as m6A9,10,21–25, m5C10,23,26 and Ψ9,10,20,27, or rely on comparative analyses to detect condition-specific modification changes21,23,28,29. Both strategies are unable to resolve the full spectrum of endogenous RNA modifications or reveal their complex interactions. While a few attempts have been made to simultaneously identify multiple RNA modifications and their associations9,10, these models remain constrained by the narrow scope of in vitro synthesized modification types, which limits their generalizability to unseen or uncharacterized modifications. Thus, the systematic and unbiased characterization of the full epitranscriptomic landscape and its underlying regulatory crosstalk remains a fundamental challenge.

To address these limitations, we present ORCA (Omni-RNA modification Characterization and Annotation), a deep-learning framework for comprehensive profiling of RNA modifications and their interactions at isoform and single-molecule resolution. ORCA employs an adversarial learning strategy to capture both signal- and sequence-level variations arising from the mixed stoichiometry of modified and unmodified nucleotides, thereby overcoming the limited detection scope of current DRS-based approaches. Extensive benchmarking shows that ORCA serves as a powerful tool for unbiased detection and stoichiometric quantification of RNA modifications, even for modification types absent from the training data, demonstrating ORCA’s broad generalizability across diverse modification types. Applying ORCA to human cell lines, we expand the known repertoire of RNA modification sites and uncover the widespread interactions among different modifications across transcript isoforms. Notably, ORCA reveal intricate cooperative and competitive relationships between neighboring modification sites, suggesting the complex crosstalk between RNA modifications and splicing regulation. Collectively, ORCA provides a robust and versatile approach for mapping the full spectrum of RNA modifications, revealing the regulatory complexity and isoform-specific crosstalk in the eukaryotic epitranscriptome.

Results
Deep-learning-based detection and annotation of various RNA modifications from direct RNA-seq data
To enable generalized detection of diverse RNA modifications from nanopore direct RNA-seq data, we developed a deep-learning framework (ORCA) to systematically identify multiple RNA modification types. Briefly, ORCA first aggregates the raw current signals and basecalled sequences from all reads aligned to a given genomic region, focusing on a 9-nucleotide window centered on each candidate site (Fig. 1a and “Methods”). Since RNA modifications exhibit mixed stoichiometry17,30, where not all copies of a given base are modified, modified positions should be characterized by elevated skewness in signal intensity distributions and increased basecalling errors (Supplementary Fig. 1). Thus, ORCA employs these polymorphic features to detect the presence of RNA modifications across the transcriptome. Afterward, ORCA integrates prior knowledge from established RNA modification databases for effective annotation of a wider range of modification types.

To accurately predict the presence of RNA modifications based on signal- and base-level features, we first constructed a robust and diverse training set comprising six types of RNA modifications from the in vitro synthetic ELIGOS sequences31 (Fig. 1b). Synthetic transcripts containing one of the six modified bases (m6A, m5C, Ψ, m1A, hm5C, and 5fC) or four canonical bases were randomly sampled and combined to simulate varying stoichiometries and sequencing depths (“Methods”). To mitigate k-mer bias from the limited sequence diversity of the ELIGOS sequences10, raw sequence or absolute current levels features were deliberately excluded to ensure that these features could represent generalized modification status without sequence preference. In total, over 7,000,000 sites were generated, with positions containing >10% modified transcripts designated as the positive set15. To develop a generalized model capable of accurately predicting diverse RNA modifications without being restricted to specific types, we implemented a domain adversarial learning framework32 (Fig. 1b). Here, a feature encoder comprising two LSTM layers that process the sequence in opposite directions was trained to capture contextual and sequential features and predict modification presence (modScore) and stoichiometry. Notably, a domain classifier was adversarially trained to minimize the models’ ability to discriminate between different modifications using the encoder’s output (“Methods”). This adversarial training strategy forced the encoder to learn generalized features that are shared across modifications, ensuring robust representation of modification status beyond training types while suppressing overfitting to the modification types used for training.

Considering occurrence of the same RNA modification across different transcriptomic positions often share conserved sequence contexts or signal patterns15,33, we implemented a transfer-learning strategy for modification type annotation (Fig. 1c). First, an autoencoder was trained to project all predicted modification sites into a low-dimensional embedding space using both signal- and base-level features, as well as k-mer frequency profiles that capture motif preference of modification sites. These modification sites were then provisionally annotated using public RNA modification databases (RMBase v3.034 and DirectRMDB35). Subsequently, the model was fine-tuned to predict the type of annotated modification sites, with unannotated sites randomly sampled as negative controls to reject low-confidence predictions and suppress false discoveries. Finally, the classifier’s predictions were transferred to all unannotated sites, enabling comprehensive and rigorous identification of unannotated RNA modification sites while minimizing dependence on pre-labeled training data. This framework ensures that ORCA can achieve stringent modification sites identification and annotation, with inherent flexibility to integrate new modification types with the emerging updates of RNA modification resources36,37.

Performance evaluation of RNA modification detection
To evaluate ORCA’s performance in detecting various RNA modifications, we first assessed its sensitivity and accuracy using the synthetic ELIGOS dataset31. Training datasets for six in vitro synthesized modifications were generated as described above, and 5-fold cross-validation was applied to evaluate prediction accuracy across modification types. As shown in Fig. 2a and Supplementary Fig. 2a, ORCA achieved high recall and precision across all six modifications, with an average area under the precision–recall curve (AUPRC) of 0.95 and average area under the receiver operating characteristic curve (AUROC) of 0.94. To quantify the overall performance, we further calculated the F1-score, which balances sensitivity and false discovery rate (FDR) (Fig. 2b). Across all modification types, ORCA consistently attained high F1-scores (0.971-0.976), reflecting its reliable and accurate detection capability. Given that modified sites represent only a small fraction of the transcriptome, we further estimated the false discovery rate of ORCA using an in vitro transcribed (IVT) human mRNA transcriptome31 devoid of endogenous modifications. As shown in Fig. 2c, ORCA exhibited a low false discovery rate of 2.25% using the default threshold (modScore > 0.9), demonstrating superior false discovery scores compared with multi-modification detection tools and comparable performance relative to several modification-specific methods (Supplementary Fig. 2b). In addition, ORCA demonstrated robust stoichiometry prediction (Supplementary Fig. 2c), showing its ability to accurately quantify diverse RNA modifications using the integrated signal- and base-level features, Together, these results suggested that ORCA provides accurate and unbiased prediction of multiple RNA modification types.

To assess the performance of ORCA in real-world transcriptomes, we benchmarked its ability to detect m6A modifications using DRS data from Mettl3 knockout (KO) and wild-type (WT) mouse embryonic stem cells31. We applied ORCA to predict modification sites in individual samples and analyzed site-specific differential modification levels after Mettl3 knockout (Methods). As shown in Fig. 2d, Mettl3-KO cells exhibited a significant global reduction in RNA modification, with 17.58% of modification sites showing >0.2 stoichiometric reduction. In contrast, only 5.78% of sites retained increased modification stoichiometries, consistent with Mettl3’s role as a primary m6A methyltransferase. For comparison, we evaluated both typical m6A-specific models (CHEUI-solo23, EpiNano-SVM31, m6Anet22, and TandemMod10) and comparative-based methods (CHEUI-diff23, EpiNano-Error21, Nanocompore29, and xPore28) for detecting differential m6A sites. The m6A sites identified by miCLIP212 and GLORI15 were collected as ground-truth benchmarks. To ensure a fair comparison between m6A-specific and comparative-based tools, the performance of m6A prediction was evaluated at both single-base level (only adenosines within DRACH motifs were considered modified) and 5-mer level (all nucleotides within DRACH motifs were treated as modified), respectively38. ORCA achieved an AUPRC of 0.42 at the single-base level (Supplementary Fig. 2d) and 0.43 at the 5-mer level (Fig. 2e), matching the performance of state-of-the-art m6A-specific and comparative-based algorithms. Among the top differentially modified sites (ranked by change of modScore), ORCA exhibited the highest proportion of sites overlapping DRACH motifs or modified 5-mers from m6A sequencing methods (Fig. 2f), confirming its high accuracy in identifying the biologically relevant m6A modifications.

To further validate the versatility of ORCA, we applied it to detect m5C modifications in HeLa cells following NSUN2 knockout23, and then evaluated the performance of ORCA and other tools against m5C sites reported by BS-seq39, bsRNA-seq40, and RNA-BisSeq41. For each tool, changes in predicted stoichiometry were then calculated, and a chi-squared test was employed to measure the reduction in m5C sites. Among the m5C sites that consistently detected across two biological replicates, ORCA identified a higher proportion of downregulated sites than most existing tools (Fig. 2g). While CHEUI-diff reported a marginally higher fraction of downregulated sites, ORCA uncovered a significantly larger absolute number of m5C sites with statistically significant downregulation (Fig. 2h and Supplementary Fig. 2e), indicating its strong ability to detect NSUN2-dependent m5C alterations.

Beyond m5C, we further tested ORCA’s ability to detect Ψ modifications in the ribosomal RNA dataset27. Compared to established tools (Tombo42, TandemMod10, and NanoPSU27), ORCA identified a comparable high number of validated Ψ sites (73/89) to TandemMod (75/89, Fig. 2i) and also showed strong orthogonal overlap with each approach (Supplementary Fig. 2f), suggesting its accuracy in Ψ detection. Collectively, these results demonstrated ORCA’s capacity as an effective framework for detecting a wider range of RNA modifications across different experimental conditions and modification types, including m6A, m5C and Ψ, and supported its utility for accurate and robust transcriptome-wide RNA modification discovery superior to canonical modification-specific and comparative-based tools.

ORCA enables zero-shot detection of unseen RNA modification types
A key limitation of current DRS-based RNA modification detection tools is their reliance on sophisticated training datasets, which typically derived from in vitro-synthesized transcripts10. However, many modification types are challenging to synthesize in vitro43, limiting the development of these modification-specific models. To determine whether our strategy is applicable to modifications not included in the training set, we comprehensively benchmarked the performance of ORCA using three complementary approaches: (1) zero-shot prediction of unseen modifications using the ELIGOS dataset, (2) prediction of modifications absent in training set with reference RNA modification sequencing data, and (3) evaluation of ribosomal RNA modifications from mass spectrometry (MS)-based databases (Fig. 3a).

To systematically evaluate the performance of ORCA in detecting unseen RNA modifications, we first assessed the precision and recall in zero-shot prediction using synthetic ELIGOS sequences31. For each target modification, we iteratively excluded it from the training set and trained ORCA on different combinations of other modifications to evaluate ORCA’s ability to generalize from arbitrary subsets to unseen targets. For zero-shot prediction of target modifications the absent in training set, ORCA achieved high prediction accuracy ~90% across all modification types and maintained considerable recall rates >30% for most modifications (Fig. 3b). Both accuracy and recall increased with the inclusion of additional training modification types (Fig. 3b and Supplementary Fig. 3a), indicating that ORCA can effectively extract generalized cross-modification features. In addition, we further evaluated ORCA’s ability to estimate stoichiometry for previously unseen modifications. Strikingly, our adversarial learning framework accurately encapsulated stoichiometry in the zero-shot prediction of all six modification types, achieving a strong linear correlation (average Pearson’s correlation coefficient=0.76) between predicted and ground-truth stoichiometries (Supplementary Fig. 3b). These results confirmed that ORCA can be effectively adapted to detect and quantify unseen RNA modifications without requiring prior training data for these modification types.

To further evaluate ORCA’s ability to predict unseen RNA modifications transcriptome-wide, we employed three different high-throughput RNA modification sequencing datasets to validate its applicability. First, we performed ONT direct RNA-seq and m6A-SAC-seq17 on the same mouse brain sample. To test ORCA’s capacity for zero-shot prediction, m6A was excluded from the training set to construct an m6A-absent model, which was then used for transcriptome-wide modification prediction. In total, 1000 m6A sites detected by m6A-SAC-seq were covered in the nanopore DRS data, of which 70.5% were confidently predicted as modified (modScore > 0.9) by the ORCA model trained without m6A-specific data. Moreover, these identified m6A sites showed significantly higher modification probabilities than randomly sampled DRACH motifs (P<1×10−308, Wilcoxon rank-sum test; Fig. 3c), confirming ORCA’s capacity to detect transcriptome-wide m6A without prior training on this modification. Similarly, we further evaluated ORCA’s ability to detect 2’-O-methylation (Nm) and inosine (I) modifications which were not included in the training set. For Nm detection, we analyzed public DRS data from mESC cells31 and benchmarked ORCA predictions against Nm sites identified by 2’-OMe-seq44. As shown in Fig. 3d, ORCA accurately predicted Nm modifications, with 74.4% of reference Nm sites successfully identified, and an overall significant enrichment of reference Nm and random control sites was also observed among ORCA’s prediction (P=1.30×10−24, Wilcoxon rank-sum test). For inosine prediction, we employed public DRS data from wild-type and FY-ADAR2 yeast strains engineered to express human ADAR2, which introduce A-to-I editing in yeast without an endogenous ADAR system45. A-to-I editing sites were identified as reference using Illumina RNA-Seq data from the same project (“Methods”), and differentially modified sites between WT and hADAR2-expressing yeast were identified using ORCA. As expected, reference inosine sites showed significantly higher modification scores in hADAR2 yeast, whereas no such patterns were observed for randomly sampled background adenosines (Fig. 3e). Taken together, these results indicated the ability of ORCA to detect previously unseen RNA modifications at the transcriptome scale, even in the absence of modification-specific training datasets.

Given that ribosomal RNA harbors a diverse array of RNA modifications—many of which are undetectable by either modification-specific or comparative-based tools (Fig. 3f), we employed human and yeast ribosomal RNA sequencing datasets to evaluate ORCA’s performance in predicting this broader spectrum of RNA modifications. Across human 28S and 18S rRNAs, ORCA successfully identified 78.1% of the 224 modification sites spanning 13 chemically distinct modification types (e.g., N4-acetylcytidine (ac4C) and m1acp3Ψ) supported by SILNAS mass spectrometry46. In parallel, a low false discovery rate (10.0%) was observed for unmodified bases, indicating high specificity in distinguishing modified from unmodified sites. Similarly, ORCA achieved comparable performance on yeast 18S and 25S rRNAs, accurately predicting 79.8% of orthogonally validated modifications, further supporting its robustness in identifying a wide range of RNA modification types across species (Supplementary Fig. 3c, d).

Furthermore, we evaluated ORCA’s ability to detect non-natural 4-thiouridine (4sU) modifications. Applying ORCA to K562 4sU pull down and DMSO control samples47, ORCA detected significantly elevated modification scores for U-containg 5-mers in 4sU pull down samples compared with non U-containing controls (P=1.38×10−20, one-sided Wilcoxon rank-sum test; Supplementary Fig. 4a). The predicted 4sU levels also strongly correlated with nascent RNA expression measured by orthogonal Illumina sequencing (P=1.67×10−16, Supplementary Fig. 4b), demonstrating ORCA’s sensitivity to 4sU incorporation. Taken together, these results demonstrate that ORCA can overcome the limitations of modification-specific models, enabling transcriptome-wide discovery of unseen RNA modifications through its generalized adversarial training strategy.

ORCA uncovers the transcriptome-wide landscape of a broad range of modifications
The comprehensive identification of diverse RNA modifications is essential for understanding the post-transcriptional processing of RNAs. However, current DRS-based tools typically rely on modification-specific models that target a limited subset of modifications or utilize comparative strategies to detect changes between different condition pairs. These limitations make it challenging to profile a broad spectrum of RNA modifications in individual samples. To further demonstrate ORCA’s ability to simultaneously detect multiple RNA modifications per sample, we applied it to the previously described Mettl3-KO mESCs dataset. All predicted modification sites were ranked by modScore and compared against known modification sites in the RMBase 3.0 database34, as well as sites predicted by modification-specific tools. Among the top-ranked predictions, over 40% were supported by either a public database or ONT-based RNA modification detection tools (Fig. 4a), indicating the high reliability of ORCA’s prediction. Specifically, a comprehensive catalog of well-characterized modifications including m6A, m5C, inosine (I), Ψ, m7G, m1A and 2′-O-methylation (Nm) were detected within the top 10,000 predicted sites (Fig. 4b). For example, 41.9% and 30.9% of predicted sites in mESCs and HeLa cells were supported by m6A sites, while 16.1% and 5.2% of sites were supported by Ψ modifications. This high concordance underscored ORCA’s ability to accurately resolve multiple RNA modifications in a single analysis.

For instance, ORCA identified adjacent modification sites in the 3’ UTR of the Lars2 transcript, where one Nm, two m5C, and five Ψ modification sites were consistently detected in both WT and Mettl3 KO mESCs. Notably, a previously unannotated site exhibited a dramatic reduction in modScore upon Mettl3 knockout, suggesting it may represent an unseen m6A site (Fig. 4c). Similarly, two m6A sites in the 3’ UTR of Ets2 transcript were robustly detected in WT mESCs but were absent in Mettl3 KO samples, consistent with substantial alterations in both sequence- and signal-level features following loss of the key m6A methyltransferase complex (Fig. 4d). Moreover, these predicted sites were also consistently supported by both modification-specific and comparative-based algorithms. We further extended this analysis to WT and NSUN2 KO HeLa cells. For exemplary modification sites in the 3’ UTR of RPL13A and in the small nuclear RNA RNA5-8SN1, ORCA revealed a selective reduction in m5C, but not in m6A or Ψ levels, consistent with the specificity of NSUN2 as an m5C methyltransferase41 (Supplementary Fig. 5). Taken together, these results demonstrated that ORCA enables the simultaneous detection of diverse RNA modifications while accurately resolves biologically relevant stoichiometry changes upon perturbation of specific modification writers.

Transfer learning of sequence and signal features enables accurate discovery of previously unannotated RNA modification sites
Based on accurate modification presence prediction, we further developed a transfer learning framework to annotate modification types using both sequence- and signal-level features of high-confidence sites curated from public databases (Fig. 5a and “Methods”). In addition to the features used for modification presence prediction, we incorporated k-mer occurrence frequencies to capture sequence similarity specific to each modification type. Given that chemical modifications affect current signals across a 5–6 nucleotide window as the strand passes through the nanopore48, a multi-task learning model was implemented to simultaneously predict both modification types and positional phase. The prediction results were then filtered based on the consistency between predicted modification types and the corrected nucleotide phases, and only predictions where the modification types matched the corresponding nucleotide position were retained to ensure accurate modification assignment. To avoid over-assignment of uncharacterized modifications to known categories, unannotated sites were also sampled as negative controls during training to ensure stringent and reliable identification of known modifications. Finally, the trained model was transferred to predict modification types for previously unannotated sites, enabling the discovery of RNA modification sites absent from existing databases. Overall, this approach balances the accuracy of modification type prediction with the sensitivity to uncover epitranscriptomic features.

To rigorously evaluate the performance of RNA modification annotation, ORCA was employed using a K562 DRS dataset from the SGNex project49. A total of 48,377 sites across 7 RNA modification types were then annotated with high confidence NGS-supported sites from RMBase 3.034 and DirectRMDB35, and the modification annotation model was trained as described previously (Fig. 5a). In cross-validation, ORCA achieved >90% precision and > 40% recall for most modification types (Fig. 5c, d), demonstrating robust accuracy and sensitivity. To further assess the specificity of modification type prediction, we iteratively masked each modification type during training and quantified the misclassification of masked sites into other categories. As shown in Fig. 5e and Supplementary Fig. 6a, ORCA maintained an average of >0.83 accuracy across all modification types. Moreover, ablation of the background negative-control class markedly increased false discovery (Supplementary Fig. 7), confirming that our strategy ensure accurate and stringent modification type prediction while minimize false positive assignments of unannotated modifications.

The trained model was subsequently applied to predict modification types across all unannotated modification sites. In total, 42,449 previously unannotated RNA modification sites were identified, whose gene body distributions closely consistent with those in the established databases (Fig. 5f). Notably, ORCA largely expanded the catalog of current modification sites, annotating 29% additional m6A sites and dramatically increasing the low-abundance modifications: >400% more m5C, Nm, Ψ and inosine, >1030% more m7G, and >178% more m1A sites compared to existing annotations (Fig. 5g). To validate these predictions, we then performed de novo motif analysis using XSTREME50 on previously unannotated m5C sites. Two canonical m5C motifs CUCC (88.3% of ORCA-annotated sites) and CGGG (8.8% of ORCA-annotated sites) were identified, aligning with known NSUN651 and NSUN2-dependent39 m5C sites. For example, ORCA predicted an m5C site in the CDT1 3’ UTR that was absent in both RMBase 3.0 and DirectRMDB, but was independently validated by UBS-seq16 in HeLa cells (Fig. 5i). Besides, ORCA also demonstrated high specificity for annotating m6A sites, with 26.9% sites supported by GLORI15, which was consistent with the validation rate of database-curated m6A sites (37.6%, Supplementary Fig. 6b). For other modifications, high validation rates were also observed for Ψ (BID-seq14), m7G (20% overlap with QKICLIP-seq11 peaks), and m1A (18.5% by m1A-seq13) (Supplementary Fig. 6c–e). Taken together, these results demonstrated the effectiveness of ORCA’s label-transfer learning strategy in discovering and annotating RNA modification sites that were previously unannotated in existing databases with high confidence.

Furthermore, we evaluated whether database composition introduces biases associated with modification detection technologies or cell line origins. Although cross-technology overlap for the four major modification types was generally limited (Supplementary Fig. 8a, b), ORCA’s annotation model trained on individual assays achieved high recall within orthogonal datasets, and performance further improved when complementary assays were incorporated while maintaining low FDR (Supplementary Fig. 8c, d). Consistent results were observed in cross-cell-line validation, where restricting training to a single cell line reduced annotation sensitivity but did not affect precision or FDR (Supplementary Fig. 9). Together, these results demonstrate that limited overlap across technologies or cell lines does not compromise annotation accuracy, and that integrating multiple public resources effectively mitigates technology-specific biases and improves sensitivity without sacrificing precision.

Characterization of RNA modification landscape and its regulatory crosstalk in human cell lines
To demonstrate the applicability of our method, we applied ORCA to characterize the RNA modification landscape in human cell lines DRS data from the SGNex project49. In summary, a total of 98,586 sites were detected across all samples, with 10,954 modification sites per cell line. Notably, 70.2% of these sites were consistently detected in at least two cell lines (Fig. 6a), which is consistent with the reported stable m6A modifications shared across human cell lines28. To further investigate the spatial associations between different modification types, we further calculated the genomic distances between adjacent modifications. Strikingly, a substantial proportion of modifications (33%) occurred within 20-nt of each other (Supplementary Fig. 10a). We therefore clustered proximal sites using a 20-nt window, yielding 13,633 clusters with an average of 2.85 modification per cluster (Fig. 6b, c). As m6A was the most abundantly detected modification, most clusters were m6A-enriched, while a high degree of association between m6A and other modifications, such as m5C and m1A, was also observed, highlighting the complex spatial organization and potential crosstalk among neighboring modification types (Fig. 6d).

To further investigate the regulatory interplay between these neighboring RNA modifications, we applied an expectation-maximization (EM)-based model to estimate the single-molecule co-occurrence patterns within modification clusters (Fig. 6e and “Methods”). Among 443,361 modification clusters, 7719 exhibited significant co-modification, while 39,906 were exclusively modified with competitive exclusion (Supplementary Fig. 6b). First, cooperative modified clusters were prioritized for downstream analysis. As shown in Fig. 6f, frequent co-occurrence of different modifications was substantially observed, with m5C and m6A emerging as the most prevalent combinatorial pattern. To further validate these predictions, we leveraged m6A-SAC-seq and m5C-TAC-seq52 datasets to extract short-read level co-modification evidence. For instance, a cooperative modification of a pair of m6A sites spaced 8 nucleotides apart in the 3’ UTR of DNAJB1 was identified in IM95 DRS data and independently confirmed in the HeLa m6A-SAC-seq data (Fig. 6g). Similarly, two co-occurring adjacent m5C sites in the 3’ UTR of HDGF were detected in Hct116 DRS data, which was also validated by HeLa m5C-TAC-seq (Supplementary Fig. 6c). These results demonstrated the cooperative modification between different types and also validated accuracy of ORCA in resolving spatial co-existence of neighboring RNA modifications at single-molecule resolution.

To investigate the interplay between RNA modifications and splicing regulation, we then focused on exclusively modified clusters that exhibited significant isoform-specific changes in K562 cells, and integrated ENCORE eCLIP-seq data53 of K562 cells to assess regulatory associations. The binding patterns of RNA-binding proteins, including splicing factors and RNA modification-associated proteins (writers, erasers, and readers, WERs) within these modification clusters were further calculated. Notably, many isoform-specific modification clusters, particularly those associated with m6A, showed significant enrichment of splicing factors and modification-associated WERs, suggesting the widespread coupling between m6A modifications and alternative splicing events (Fig. 6h). For example, two splicing factors ELAVL154 and FMR155 were significantly enriched in isoform-specific m6A/m5C modification clusters, consistent with the FMR1’s preference for binding m6A-modified RNAs56,57 (Fig. 6i, j). In one case, transcript-level analysis of RBIS revealed an isoform-specific exclusion pattern between neighboring m6A and m5C sites at exon 4 (Fig. 6k). m6A-modified reads were strongly associated with an upstream skipped exon, whereas the exon was consistently retained in m5C-modified reads. Furthermore, strong eCLIP-seq peaks for splicing factors MBNL1 and the RNA-binding protein FMR1 were also detected in the same region, consistent with their established roles in alternative splicing and m6A-mediated splicing regulation56. Together, these findings demonstrate that ORCA enables systematic characterization of the interactions between RNA modifications and splicing, offering a powerful platform for dissecting the multilayered regulation of the eukaryotic transcriptome.

Methods
Animal experiments
All mice used in this study were adult C57BL/6 mice and were purchased from SiPeiFu Biotechnology. Two adult mice were used for brain tissue dissection, with one mouse used for RNA 002 sequencing and the other for RNA 004 sequencing. Animals were maintained under conventional specific pathogen-free conditions and a 12-h light/12-h dark cycle at 25 °C and 40–60% humidity.

RNA isolation and nanopore direct RNA sequencing
Total RNA was extracted from two healthy adult mice brain using TRIzol (Invitrogen) according to the manufacturer’s instructions. RNA integrity and quality were assessed using the Agilent 5200 Fragment Analyzer System. Nanopore direct RNA-seq library was prepared using the Direct RNA Sequencing Kit (SQK-RNA002) from Oxford Nanopore Technologies following the manufacturer’s protocol and sequenced on an R9.4.1 flow cell (FLO-MIN106D) using a MinION Mk1B device for 72 h. An adult mouse brain direct RNA-seq library was also generated using the SQK-RNA004 sequencing kit and sequenced according to the manufacturer’s instructions on an FLO-MIN004RA flow cell for 72 h.

m6A-SAC-seq and data analysis
For m6A-SAC-seq experiments, 1 μg of total RNA was subjected to ribosomal RNA depletion using the RiboErase kit (human/mouse/rat, Kapa Biosystems). The rRNA-depleted total RNA was used directly for m6A-SAC-seq library preparation following the protocol described by He et al.17. Briefly, m6A modifications were selectively converted into allyl-labeled derivatives by MjDim1, followed by iodine-induced intramolecular cyclization. These modifications were subsequently converted into sequence mutations during by HIV-1 RT reverse transcription and detected via Illumina sequencing.

Sequencing reads were trimmed using Cutadapt66 (v2.10) and Fastp67 (v0.23.4). Reads aligning to rRNA sequences were removed using Bowtie268 (v2.3.4.3) and Samtools69 (v1.18). Cleaned reads were then mapped to the mm10 reference genome using STAR70 (v2.7.10b). PCR duplicates were collapsed using UMICollapse71 (v1.0.0), and deduplicated BAM files from biological replicates were merged with Samtools. Strand-specific BAM files were generated and processed with Samtools mpileup. Somatic variants were called using VarScan72 (2.3.9), and candidate m6A sites were identified based on mutation profiles and the presence of DRACH motifs.

Nanopore data preprocessing and feature extraction
The GRCh38 (human) and GRCm38 (mouse) reference transcriptome were obtained from the Ensembl database. For RNA002 data, raw nanopore fast5 files were basecalled using Guppy (v6.3.8) with rna_r9.4.1_70bps_hac model. For RNA004 data, raw nanopore pod5 files were basecalled using Dorado (v1.1.1+e72f1492) with rna004_130bps_hac@v5.2.0 model. POD5 files were transferred to BLOW5 format using bluecrab73 (v0.4.0) p2s and slow5tools74 (v0.8.0) merge commands. Basecalled reads were then aligned to the reference transcriptome using Minimap275 (v2.21, with the parameters “-ax splice -N 0 -uf -k14 --cs --secondary=no”). The alignment results were processed with samtools mpileup (v1.11) to generate per-base summary statistics. Ionic current signals were aligned to the reference sequence using the f5c76 eventalign (v1.11), an accelerated implementation of Nanopolish77, with the parameters ‘--min-mapq 0 --rna --signal-index --scale-events --secondary=no --collapse-events’ for RNA002 and ‘--pore RNA004 --min-mapq 0 --rna --signal-index --scale-events --secondary=no --collapse-events’ for RNA004 reads.

For modification presence prediction, both signal-level and sequence-level features were extracted within a ±2 k-mer window surrounding each candidate site. For signal-level features, raw electrical events from the eventalign output were standardized using the method defined in Nanopolish. Specifically, each event’s mean signal level was normalized by subtracting the expected reference mean and then dividing by the reference standard deviation. This normalization accounts for variation in signal intensity across different sequence contexts. The standardized signal values from all reads aligned to the same genomic position were then aggregated and interpolated into a fixed-length vector of 50 values to ensure consistent input dimensions for the model. For sequence-level features, rate of insertions, deletions, and mismatches, as well as statistical metrics including the mean, median, and standard deviation of sequence quality scores of all aligned at each position were extracted based on the ‘samtools mpileup’ result.

For modification-type annotation, three categories of features were included as input for the transfer-learning model. K-mer occurrence features were derived from the frequency of all 256 possible 4-mer motifs within an 11-nucleotide window centered on each modification site. Besides, signal-level features were computed using the event-level means and standard deviations across the window surrounding each modification site. Then a Gaussian mixture model was applied to partition each feature into divide into modified and unmodified clusters, and the mean, variance, and covariance of each component were extracted as model input. Finally, sequence-level features were obtained using the same strategy as described above.

ORCA model design
ORCA comprises two neural network models designed for predicting the presence of RNA modifications and inferring their types. For modification presence prediction, ORCA adopts a domain-adversarial neural network architecture composed of an encoder and two classifier branches. The encoder utilizes a bidirectional LSTM network to capture contextual and sequential dependencies from both sequence- and signal-level features within a 11-nucleotide window surrounding each candidate site. The encoded representations are then simultaneously passed into two parallel branches: (1) a modification predictor for predicting modification presence (modScore) along with an estimate of stoichiometry; and (2) a domain classifier with a gradient reversal layer aims to distinguish between different RNA modification types. The model is trained adversarially to optimize the encoder and modification predictor to accurately detect modification presence, while the encoder is simultaneously trained to learn representations that minimize the performance of the domain classifier, ensuring a generalizable representation of modification presence across a diverse range of modification types without introducing modification-specific bias.

For modification type inference, ORCA employs a transfer-learning framework consisting of an autoencoder and two classifier modules. During the pretraining phase, the autoencoder learns the global low-dimensional representation of all predicted modification sites. The encoded features are subsequently passed to multi-task prediction to produce probabilities corresponding to different modification types, and the phase represents the exact modification position in the input window. Then, only modification type predictions that correspond to the base at the inferred modification position are retained as valid outputs.

Model training
To train the presence prediction model, we generated a labeled dataset based on the public ELIGOS resource. Specifically, six in vitro transcribed RNA samples containing individual modifications (m1A, m6A, m5C, hm5C, 5fC, and Ψ) and a control sample composed of only canonical bases were obtained. For each modification type, raw nanopore reads were aligned to the reference sequence using Minimap2. Then, modified and unmodified reads were randomly sampled and combined to generate 3,000 training samples at varying sequencing depths (10–200) and modification rates (0–1) of individual modification sites. The sequence- and signal-level features were extracted as described above, resulting in a training set that simulates the realistic features of RNA modification types and stoichiometry levels. In total, over 7 million training samples were generated. Positions with >10% modified reads were labeled as positive samples. For each modification type, the dataset was randomly split into training and testing set at 4:1 ratio. Model optimization was conducted using the AdamW optimizer with a learning rate of 0.0005. The loss functions were defined as follows: negative log-likelihood loss (NLLLoss) for modification presence prediction, NLLLoss for domain label prediction, and mean squared error (MSELoss) for stoichiometry regression.

To train the modification type annotation model, we utilized 25 direct RNA-seq datasets generated from the MinION/GridION platform, spanning nine human cell lines from the SGNex project49. RNA modifications were first identified using the ORCA prediction module with a modScore threshold of 0.9 to ensure high-confidence predictions. Predicted modification sites were then annotated based on the presence of NGS-supported modifications from RMBase v3.034 and DirectRMDB35 using a distance threshold of ±2 nucleotides. The autoencoder was initially trained on all predicted modification sites using the AdamW optimizer with a learning rate of 0.002 and MSELoss to learn global feature representations. The model was then fine-tuned using the annotated subset to enable prediction of both modification type and phase. These downstream tasks employed cross-entropy loss and were optimized with AdamW at the same learning rate. For each modification type, the labeled dataset was split into training and testing set in a 1:4 ratio. To mitigate class imbalance, each mini-batch was constructed to contain similar numbers of samples from each modification type, preventing the loss function from being dominated by any single class. To avoid over-assignment of uncharacterized modifications to known categories, we also included 3,000 unannotated sites as negative controls during training to ensure a stringent and reliable identification of known modification types.

EM-based prediction of modification interactions
ORCA outputs read-level modification predictions for neighboring modification pairs through the EM-based model, which enables assessment of read-level linkage or mutual exclusivity between modification events. After performing RNA modification presence and type prediction and across all human cell line samples, only high-confidence sites detected in at least two samples were retained for downstream analysis. To define local modification clusters, transcriptomic distances between neighboring modification sites were calculated, and sites within 20 nucleotides were iteratively merged into the same cluster. To minimize potential signal interference caused by adjacent modifications, any neighboring sites located within 4 nucleotides of each other were excluded from clustering.

To evaluate potential co-modification between modification sites within each cluster, we employed a Local Outlier Factor (LOF) score-based strategy to infer the read-level modification states. For each site, raw signal features including event level mean, standard deviation, and dwell time were extracted from a ±5 k-mer (15 nucleotide) window. LOF scores were computed for each read and normalized to the range [0,1]. For every pair of neighboring modification sites, each read was represented as a two-dimensional coordinate (x, y), where x and y are the normalized LOF scores at the two positions, respectively.

Given the site-specific distributions of LOF scores, we employed an expectation-maximization (EM) algorithm to classify each read into one of four canonical modification states, including: d0 centered at (0,0) for unmodified reads; d1 centered at (1,1) for dual modification; and d2 and d3 centered at (1,0) and (0,1), respectively, for single-site modifications. Each state dj was assigned an initial mixing weight θj = 0.25. In the expectation step, the posterior probability that read ri (xi, yi) belongs to state dj was calculated as:
γij = θj * P(ri|dj) / Σk=0 to 3 θk * P(ri|dk)
where the likelihood P(ri|dj) = Euclidean(ri, dj) * s(ri, dj) was defined as the product of the Euclidean distance between ri and the center of dj, and a dispersion score s(ri, dj), such that:
s(ri, d0) = s(ri, d1) = 1 - |xi - yi|
s(ri, d2) = s(ri, d3) = 1 - |xi + yi - 1|
This formulation captures both the geometric proximity of the read to a canonical modification state and the consistency of LOF across the two sites.

In the maximization step, the mixing weights were updated as the mean of the posterior probability across all N reads:
θj = 1/N * Σi=1 to N γij
The EM process was iteratively repeated until convergence, defined by the change in posterior weights falling below a predefined threshold. After convergence, each read was assigned to the modification state with the highest posterior probability γij, enabling stratified analysis of co-modified, mutually exclusive, and unmodified read populations within each cluster.

To further quantify the interaction between each pair of modification sites, we defined a linkage score for each candidate pair as:
Linkage score = (θ1 * θ0) + (θ2 * θ3)
A modification site was considered cooperatively modified if θ1 > max (θ2, θ3), indicating an enrichment of simultaneously modified reads. Conversely, modification site pairs were considered as mutually exclusive if the linkage score was less than -0.2 and θ1 < min (θ2, θ3), suggesting that the two modifications tend not to co-occur in the same read.

False positive rate evaluation
To evaluate the false discovery rate of ORCA, the IVT human mRNA transcriptome was downloaded from the SRA database (accession number SRP166020). ORCA was employed to predict the modification presence in the IVT transcriptome. Among 7,576,597 sites with coverage great than 10, where only 2.25% (170,708) of sites were predicted as modified using the default modScore threshold of 0.9. TandemMod, EpiNano-SVM, NanoPSU and m6ANet were also applied to the same IVT transcriptome. For each tool, the cumulative distribution function of prediction scores was calculated, and false discovery performance was measured using the area under the cumulative distribution curve.

Model benchmarking
To evaluate ORCA’s performance, we benchmarked it against several representative tools for direct RNA modification detection using publicly available datasets. For m6A prediction, we included m6Anet22 (v2.0.1), TandemMod10 (v1.1.0), EpiNano21 (v1.2), CHEUI23 (v0.1), xPore28 (v2.1), Nanocompore29 (v1.0.4). For m5C detection CHEUI, TandemMod and xPore were also employed. For Ψ detection, Tombo (v1.5.1) and NanoPSU27 (v1.0) were used. Benchmarking was mainly performed on three datasets: Mettl3 knockout and wild-type (KO/WT) mouse embryonic stem cell (mESC) samples for m6A detection (SRP166020), NSUN2 KO/WT HeLa samples for m5C (SRP393373) and a mixed rRNA sample for Ψ detection (SRP329477). All sequencing reads were aligned using Minimap2. Most tools utilized using transcriptome-based alignment, except for EpiNano, which required genome-based alignment. For signal-to-reference alignment, most tools employed the eventalign module for f5c76 or Nanopolish77 with tool-specific parameter configurations. However, TandemMod used “tombo resquiggle” command to map raw ionic signals to the base called sequences. For xPore (v2.1), signal-level features were extracted using the “xpore dataprep” command, and differential modification analysis was performed with “xpore diffmod”. For ORCA, both wild-type and knockout samples were processed using the modification presence prediction model.

m6A detection in mESCs (Mettl3 KO/WT) samples
To benchmark m6A detection, we evaluated several tools using the Mettl3 knockout and wild-type mouse embryonic stem cell (mESC) sample. Experimentally detected m6A sites were obtained from the GEO database, using the union of GLORI15 (GSE210563) and miCLIP212 (GSE163491) datasets as ground-truth references.

For EpiNano (v1.2 SVM & Error mode), we first converted BAM files into TSV format using sam2tsv from jvarkit (v2023.09.07) and extracted basecalling error features with Epinano_Variants.py script. Then, Epinano_DiffErr.R and Epinano_Predict.py were used to identify differentiated modified sites across samples and to directly predict m6A modifications, respectively. For CHEUI (v0.1 diff & solo mode), feature extraction was performed using ‘CHEUI preprocess --m6A’, then followed by direct m6A prediction using CHEUI_predict_model1.py and CHEUI_predict_model2.py. In addition, signal-level differentiated analysis between samples was carried out using CHEUI_differentialRNAMod command. For m6Anet (v2.0.1), input data were preprocessed using “m6anet dataprep” command, and m6A sites were inferred with “m6anet inference”. For Nanocompore (v1.0.4), event-level signal alignments were first collapsed to the site level using “nanocompore eventalign_collapse” command, and sample comparisons were performed using “nanocompore samplcomp”. For TandemMod, raw signal features were extracted from Tombo-resquiggled fast5 files using extract_signal_from_fast5.py and extract_feature_from_signal.py, and modification predictions were performed using “TandemMod.py --run_mode predict”. For xPore (v2.1), signal-level features were extracted using the “xpore dataprep” command, and differential modification analysis was performed with “xpore diffmod”. For ORCA, both wild-type and knockout samples were processed using the modification presence prediction model.

To evaluate the m6A prediction performance, the difference in ORCA’s predicted modScore between paired samples was used to rank candidate m6A sites. The overlap between top differential modified m6A sites and the reference m6A dataset were further calculated. For tools that directly provide differential predictions, such as xPore and CHEUI-diff, the reported difference in modification rate was used as the ranking metric. For EpiNano-Error and Nanocompore, we used the delta_sum_err and P values, respectively, as provided in their outputs. To reduce tool-specific biases and ensure fair comparisons, transcriptome coordinates were converted to genomic coordinates, and only genomic sites reported in both samples were retained for downstream analysis. When multiple predictions are assigned to a single genomic site, the prediction with the highest score is selected to represent that position.

m5C detection in HeLa NSUN2 KO/WT samples
To evaluate the m5C detection performance, CHEUI, TandemMod, xPore, and ORCA were applied to NSUN2 knockout and wild-type HeLa cell samples using the same preprocessing and prediction workflow described above. Reference m5C sites identified from HeLa BS-seq39 (GSE122260), bsRNA-seq40 (GSE140995), and RNA-BisSeq41 (GSE93751) were directly downloaded from the GEO database and merged to construct a unified reference set. For comparison against WT and KO samples, the same strategy was applied here to ensure consistent genomic site-level comparisons, and only sites that were detected in both samples were retained for evaluation. The changes in predicted stoichiometry between WT and KO samples at m5C sites were calculated for comparison.

Ψ detection in mixed rRNA sequencing samples
For Ψ detection evaluation, we evaluated ORCA, Tombo, and NanoPSU using rRNA sequencing data. Reference Ψ sites supported by SILNAS-based mass spectrometry46 were used as the ground-truth. Base called reads were aligned to ribosomal RNAs from four species following the procedure described in NanoPSU27. For downstream comparison, only reads aligned to human 18S (NR_003286.4) and 28S (NR_003287.4) rRNA sequences were retained. For Tombo, Ψ sites were identified from resquiggled fast5 files using the ‘tombo detect_modifications’ command in de_novo detection mode. For NanoPSU, the recommended pipeline was applied, where alignment, remove_intron, extract_features, and prediction commands were subsequently performed for Ψ site prediction. For TandemMod, raw signal features were extracted from Tombo-resquiggled fast5 files using extract_signal_from_fast5.py and extract_feature_from_signal.py, and modification predictions were performed using ‘TandemMod.py --run_mode predict’ with the Ψ detection model.

Benchmarking of zero-shot modification detection
To evaluate ORCA’s zero-shot detection capability, training datasets were constructed by iteratively selecting different combinations of 2 to 5 RNA modification types from the full set of six RNA modifications (m6A, m5C, Ψ, m1A, hm5C, 5fC). For each combination, the selected modifications were used to train the ORCA modification presence prediction model as previously described. The trained models were then used to perform zero-shot prediction on the modification types excluded from training. Prediction accuracy and recall were calculated for each target modification type for evaluation. In addition, the Pearson correlation between predicted stoichiometry and simulated ground-truth were calculated to evaluate the generalizability of modification identification and quantification.

For zero-shot m6A prediction, an m6A-absent model was trained by excluding m6A from the training dataset. This model was then applied to direct RNA-seq data obtained from mouse brain tissue, and m6A sites detected by mouse brain m6A-SAC-Seq were used as a reference dataset. ORCA’s modScores were computed for all adenosine (A) sites located at the center of DRACH motifs across the transcriptome. The distribution of modScores between m6A-SAC-Seq-supported m6A sites and randomly selected background adenosines within DRACH motifs was compared using a two-sided Wilcoxon rank-sum test.

For Nm prediction, ORCA modScores were computed for sites supported by 2’-OMe-seq experiments44 downloaded from RMBase v3.034. All other transcriptome sites with read coverage greater than 10 were used as background controls. The distributions of modScores between reference Nm sites and background transcriptome sites were compared using a two-sided Wilcoxon rank-sum test.

For zero-shot prediction of inosine (A-to-I) RNA editing sites, a Schizosaccharomyces pombe direct RNA-seq dataset was downloaded from PRJEB46364. The ASM294v2 S. pombe genome and annotation were downloaded, and the recommended strategy described in DeepEdit45 was followed to establish a high-confidence set of inosines (A-to-I) RNA editing sites. Illumina RNA-seq reads from two hADAR2+ samples and two control samples were downloaded and individually aligned using HISAT278 and processed with bcftools69 mpileup for single-nucleotide variant calling. Candidate sites were defined as A-to-G substitutions supported by a minimum read coverage of 50 and a variant allele frequency exceeding 10%. Sites that consistently appeared in both hADAR2+ samples and were absent from any of the control samples were retained as the final ground-truth set for evaluation.

For the evaluation of rRNA modification predictions, high-confidence modification sites for human 18S/28S and yeast 18S/25S rRNAs were downloaded from a published SILNAS mass-spectrometry-based study46. The mixed-species rRNA direct RNA-seq library was then analyzed, and ORCA-predicted modScores were benchmarked against the reference modification sites to calculate the true-positive rate and false-discovery rate.

For the evaluation of non-natural modifications, publicly available 4sU-labeled K562 direct RNA-seq data47 and a matched DMSO control sample were downloaded and analyzed using ORCA. The change in modScore between the 4sU-labeled and control samples was compared between U-containing and non-U-containing 5-mers to measure the enrichment of modification signal at 4sU-incorporated sites. In addition, a transcript-level 4sU load was defined as the sum of ORCA-estimated modified counts across U-centered sites per transcript, and was compared with transcript abundance in matched short-read 4sU pulldown RNA-seq data using Spearman correlation to evaluate concordance between ORCA-derived 4sU signals and orthogonal measurements of 4sU incorporation.

Benchmarking of simultaneous prediction of multiple RNA modification types
For both mESCs and HeLa datasets, RNA modification sites were independently predicted from two WT replicates. Only sites commonly identified in both replicates were retained, and mean modScore across replicates was calculated and used for site ranking. To establish a reference modification set, predicted modification sites were compared against public databases, including DirectRMDB35 and RMBase v3.034, as well as orthogonal long-read-based prediction tools: m6Anet (m6A), EpiNano-SVM (m6A), CHEUI-solo (m6A & m5C), and NanoPSU (Ψ). To exclude SNV interference, predicted modification sites with unusually high mutation rates of insertion, deletion, or mismatch exceeding 0.5 were excluded from the prediction results. Then, the top 10,000 ranked sites were retained for downstream analysis.

Benchmarking of modification type annotation
To evaluate the performance of modification-type annotation, the K562 DRS dataset (K562_replicate6_run1) from the SGNex project was first employed for modification site prediction. A total of 48,377 high-confidence modification sites were annotated using NGS-supported sites from RMBase and DirectRMDB. The modification annotation model was trained as described above using a fivefold cross-validation. The annotation precision was measured by dividing the number of each correct modification annotation by the number assigned to it and other types. The recall rate was measured by dividing the number of correct modification annotations by the total number of modification sites in the validation site. The final trained model was applied to unannotated candidate sites for modification type inference, and prediction results were subsequently filtered based on the correspondence between the predicted modification type and the reference nucleotide. Specifically, an additional filtering step was applied for m6A prediction to retain only sites located within the consensus DRACH motif.

To evaluate the false discovery rate (FDR), each modification type was iteratively excluded from the training dataset. The trained model was then used to assess the fraction of excluded modification types incorrectly predicted as one of the included types. The number of such misclassified sites was used to quantify the FDR for each modification type.

For modification annotation validation, de novo motif analysis of ORCA-annotated m5C sites was performed using XSTREME50. For m7G validation, the QKI-CLIP peaks were downloaded from GSE193039, and intersected with the annotated m7G sites. For m1A validation, peaks identified by m1A-Seq (GSE70485) were downloaded, and each peak was expanded by ±150 nt around the center before intersecting with predicted m1A sites. For Ψ validation, single-nucleotide resolution Ψ sites from BID-seq (GSE179798) were downloaded and compared against predicted Ψ, non-Ψ sites, and unmodified controls. For m6A validation, m6A sites detected by GLORI were downloaded from GSE210563 and compared against public databases and ORCA-annotated m6A sites. To account for potential influence in neighboring nucleotides, coordinates from BID-seq and GLORI were converted to 5-mer regions prior to comparison.

Exclusive modification sites identification and RBP enrichment analysis
To quantify transcript-level variation in exclusively modified sites, each pair of unique genomic sites was treated as an individual unit of analysis. To focus on isoform-specific modifications, only modification pairs that were supported by reads spanning more than one transcript isoform were retained. For each modification pair, an n × 3 contingency table was constructed, where n represents the number of transcript isoforms, and the three columns correspond to the number of reads assigned to each distribution, including d1 (simultaneous modification at both sites), d2 (modification at the upstream site only), and d3 (modification at the downstream site only). A chi-square test was then applied to assess variation in modification patterns across different transcript isoforms. Resulting P values were used to rank the site pairs by isoform-specific modification heterogeneity.

To explore potential regulatory mechanisms underlying these isoform-specific modification patterns, the binding of RNA-binding proteins (RBPs) near these modification sites was analyzed. Genomic binding profiles from eCLIP-seq experiments for 139 RBPs in K562 cells were obtained from the ENCODE project in bedGraph format. All isoform-specific modification pairs were ranked by statistical significance, and enrichment of RBP binding around the top-ranked site pairs was assessed using a hypergeometric test. RBPs with P values ≤ 0.01 were considered significantly enriched. To further investigate the connection between isoform-specific modifications and splicing regulation, we performed an in-depth analysis using a curated list of 21 splicing factors from the SpliceAid F79 database and 7 known RNA modification regulators, including writers (RBM15), erasers (FTO), and readers (HNRNPA1, IGF2BP1, IFG2BP2, FMR1).

For single-read-level validation, raw reads from m5C-TAC-Seq (SRP459299) and m6A-SAC-Seq (SRP295164) raw reads were downloaded and aligned to the hg38 reference genome with STAR70 (v2.7.10b). Modification coordinates reported by m6A-SAC-Seq (GSE162356) and the m5C-TAC-Seq (provided in their supplementary data) were used as reference loci to quantify and visualize single-molecule co-occurrence between adjacent modification sites. The same m6A and m5C loci were independently visualized using ORCA predictions derived from IM95 and HCT116 direct RNA-seq samples, respectively.

Benchmarking on RNA004 chemistry data
For RNA004 chemistry, a prediction model was trained on curlcake IVT libraries65 containing seven synthetic RNA modification samples and an unmodified control. Training samples with defined modification fractions and read depths were generated using the same mixing strategy as for the RNA002 training sets, and the trained model was evaluated on held-out RNA004 test samples by ROC and precision–recall analysis for each modification type. For comparison with the vendor-provided state-of-the-art caller, Dorado RNA modification models (rna004_130bps_hac@v5.2.0 for m5C, Ψ, and inosine/m6A) were applied to the same curlcake reads. Read-level modification calls from Dorado were aligned to the reference using the Dorado aligner, and per-site modification statistics were obtained with modkit pileup. To place Dorado in the same simulated evaluation framework, synthetic sites spanning a range of sequencing depths and modification fractions were generated by random sampling, and the per-site modification fractions reported by Dorado were used as continuous prediction scores. Sites with a true simulated modification fraction of at least 0.1 were treated as positives, and ROC and precision–recall curves were computed using the Dorado scores.

To examine zero-shot detection on RNA004, leave-one-modification-out models were constructed by excluding each of the seven synthetic modification types in turn from the RNA004 curlcake training set and evaluating presence-prediction performance on the held-out type. De novo m6A detection across chemistries was further assessed on RNA002- and RNA004-based mouse brain direct RNA-seq libraries using models trained without m6A. For each library, ORCA modScores were calculated at adenosines located in DRACH motifs, and score distributions were compared between m6A-SAC-seq-supported sites and background DRACH positions. Overlaps of called m6A sites between the RNA002 and RNA004 datasets were summarized at the site level.

Benchmarking on ELIGOS and IVET training data
To assess the impact of training-set composition and sequencing chemistry on the presence-prediction model, additional models were trained on alternative IVT resources. One model was trained on IVET library10 from the TandemMod study using the same simulation pipeline and neural network architecture as the ELIGOS-derived RNA002 training set. Performance of the ELIGOS-trained and IVET-trained models was compared by computing ROC and precision–recall curves, in order to evaluate cross-dataset generalization and the influence of sequence diversity in the training data.

Benchmarking on the sequencing depth dependence of ORCA prediction model
The effect of read depth on presence-prediction performance was assessed by stratifying candidate sites into coverage bins and computing evaluation metrics within each bin. In the ELIGOS IVT datasets, sites with at least 10 supporting reads were grouped into four depth ranges. For each modification type and each depth bin, AUROC and AUPRC were calculated using the simulated modification labels as reference. A depth-stratified evaluation was also carried out for endogenous m6A detection in the mESC Mettl3 WT/KO datasets. Sites were grouped into the same four coverage bins, and precision–recall and ROC curves were computed separately for each bin.

Benchmarking on assay-specific and cell line-specific biases
Assay-specific and cell-line-specific biases were evaluated on K562 cell line by partitioning database-supported sites for m6A, m5C, m1A, and Nm according to the profiling technology or cell line. For each modification type in turn, annotation models with the same architecture and training procedure were trained in a multi-class setting in which training examples for the focal modification were restricted to sites from a single technology or cell line, whereas training examples for all other modification types and the background class always included all available database-supported sites. For each training configuration, recall and 1−FDR for the focal modification were quantified on sites detected by held-out technologies or cell lines.

Statistics and reproducibility
No statistical method was used to predetermine the sample size. The in vivo mouse experiment was performed once as a small proof-of-principle, using two adult C57BL/6 mice, and is not used for formal statistical comparisons or sex-specific analyses. This sample size was chosen to obtain one high-depth RNA002 and one high-depth RNA004 direct RNA-sequencing library in order to demonstrate technical feasibility rather than to estimate variability between animals. For all computational analyses, sample sizes were determined by the size of the available public datasets, and we used all reads or samples that passed predefined quality control criteria without additional subsampling. No data were excluded from the analyses.

Where statistical tests were applied, the specific test, the exact P values, and the definition of n are provided in the figure legends or Source Data. P values < 0.05 were considered statistically significant unless otherwise stated. All computational analyses were performed using Python Jupyter Notebooks with numpy, pandas, and scipy for numerical/statistical calculations and matplotlib/seaborn for plotting. All code to replicate the analysis are available as part of “Code availability”.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Dong et al presents a study of RNA modification detection method from Oxford Nanopre sequencing.\nThere are a lot of interetsing results and a few case studies. The figures and analysis are of rather high quality. I am\nimpressed by the overall quality of the manuscript.\nHowever, there are a few major flaws, making it not suitable for Nature Communications.\n1. All the results of the manuscript are based RNA002 chemistry, which is already discontinued, making the methodology in\npractice of very limited use. This is potentially the most significant flaw of the manuscirpt.\n2. There are quite a number widely occuring RNA modifications. The method accomodated only six of them, which is quite\nweird to me. It seems straightforward that the method should be trained on more RNA modification types.\n3. There are quite a few studies reporting the technical bias of different epitranscirptome profiling approaches. The\nconsistency of epitranscriptome datasets generated by two different technologies is usually less than 10%. The dataset\nstored in RMBase, DirectRMDB and etc all suffer from low-quality epitranscriptome profiling datasets. Methodology trained\non different datasets can behave very differently.", "reviewer_2": "A novel tool, ORCA, is presented for the identification of RNA modifications in data generated through the direct sequencing\nof native RNA. This is a very active area of research. Around 170 modifications exist, and tools have been developed that\nare able to identify only a small subset of these modifications. New tools are then needed for a more comprehensive\nmapping and characterization of the marks that make up the epitranscriptome.\nORCA represents a promising solution for the comprehensive mapping of the marks that are visible through the lens of\nNanopore sequencing. It leverages deep learning methods aimed at collectively identifying these marks – through\nadversarial training, i.e. forcing the learning of generalized features that are shared across modifications - and takes\nadvantage of large collections of RNA modification positions stored in public databases for annotating the identified hits. An\nextensive set of validations through orthogonal methods and datasets is included.\nAltogether, the data suggest that ORCA could represent an important new method in this rapidly evolving while competitive\nfield. The manuscript is well written, and results are well presented and convincing. Major and minor comments are listed\nbelow.\nMAJOR COMMENTS (1-13)\n1. TandemMod is one of the few tools currently available that aims at comprehensively identifying RNA modifications. I\nwould suggest a more thorough discussion of how ORCA’s features stand against it (including TandemMod model transfer\nabilities), and a more extensive comparison of the results, vs some independent reference “gold standard”. Such\ncomparison, vs TandemMod, is currently largely limited to m6A predictions.\n2. The concept of adversarial training, i.e. forcing the learning of generalized features that are shared across modifications is\ninteresting. It remains to be clarified what is the level of control that the Authors have on which/how many modifications are\nsuitable for this approach according to the collected features. Are the “shared features” necessarily shared among all\nmodifications under consideration, or they can be shared among sub-groups of a larger set of modifications under\nconsideration? Eventually, what is the search space, how many modifications are eventually targeted, considering the ability\nof the tool to proceed without training data for specific marks?\n3. What is the impact of the ELIGOS training set on ORCA's performance, in terms of experimental design (position,\ncomposition and prevalence of the modification) and quality of the Nanopore data?\n4. What are ORCA's CPU-GPU/memory requirements?\n5. The code for easily reproducing all analyses the analyses should be made available as an integral part of the study\nORCA employs a two-step strategy based on detection and annotation. The RNA modifications database (DBs) that are\nconsidered are an important resource. However this rises a number of points that should be clarified:\n6. These DBs are heavily populated by data about specific modifications (particularly m6A) and in specific cell lines (such as\nHEK and HeLa cells). How is the biased content of the considered modification DBs taken into account?\n7. Using modifications sites non-associated to DBs hits (unannotated sites) as negative controls, in the first stage, seems\nlike a problematic concept (especially for non-m6A modifications, see the previous point), even though these sites are\neventually re-considered, and potentially ex-novo annotated. The Authors should better discuss at which extent they need to\nrepresent negative controls, and what are the implications of this clearly over-simplified categorization.\n8. How is the DB’s composition impacting the results presented in Figure 6, about the clusters of adjacent modifications,\nwhich are heavily including m6A (being this over-represented in those DBs).\n9. More in general, how close ORCA's hits need to be to DB's records to be annotated by them?\nSequencing depth:\n10. How is ORCA impacted by the sequencing depth at a given genomic site?\n11. When analyzing various samples or conditions, the number of reads could greatly vary depending on technical reasons\n(e.g. the overall number of reads produced by flow cell) or for biological reasons (e.g. varying number of reads associated to\nthe same gene or genomic position as a consequence of differential gene expression between conditions/treatments). Is\nORCA in need to / able to get rid of this bias, for example by sub-sampling the reads of the larger dataset (at the sample\nand/or at the gene or site level) or by any other means?\nEvolution of the Nanopore platform:\n12. What is the compatibility of ORCA in terms of RNA002 and/or RNA004 Nanopore chemistry? Not proving that ORCA is\nable to work with RNA004 data would be a tremendous limitation.\n13. Regarding RNA004, a comprehensive comparison with Dorado's predictions of RNA modifications should be included.\nMINOR COMMENTS (14-24)\n14. Defining ORCA an \"universal framework for RNA modification detection\" is likely an overstatement, given the vast\nprevalence of untested (and currently largely untestable) modifications\n15. It remains unclear whether reads-level predictions are returned by ORCA\n16. How were the features selected (Fig. 1a)?\n17. ORCA focuses on a 9 nt window centered on each candidate site. What is the rationale, justification, and implications of\nthis choice?\n18. What about non-natural modifications, such as 4sU?\n19. It is unclear how the two-step strategy is connected to the mixed stoichiometry issue (lines 99-106)\n20. A 2.25% fraction of tested IVT sites identified as positive could be a very large number of modifications, in a dataset\nwhere none of these are expected. To which number of modifications is this percentage corresponding to?\n21. The chosen color palette often makes it very difficult to distinguish the various colors, e.g. in Fig. 2e\n22. The results about pseudouridine were not so convincing, due to the very little number of hits.\n23. When comparing ORCA's predictions with m6A-SAC-seq (line 244) could the author better described how many hits\nwere called in total by ORCA and how many of these could be associated with the DRACH motif\n24. Unclear statements: “increased sequence variability” (line 101); \".. we restricted base-level features to discrepancies\nbetween basecalled and reference sequences ..\" (lines 113); \"RNA modifications without type-specific bias\" (118); \"LSTM\nlayers\"; \"RNA modifications of the same type\" (line 128)."}, "author_response": "Response: We sincerely thank the reviewer for the thorough evaluation and the positive\nfeedback on our work.\nHowever, there are a few major flaws, making it not suitable for Nature Communications.\nResponse: Thank you for these constructive comments. We have substantially strengthened the\nmanuscript and expanded both methodological clarification and performance benchmarking to\ndirectly address the reviewer’s concerns. In the revised manuscript, we have conducted three\nmajor improvements:\n1. Updated support for the latest ONT RNA004 chemistry. We trained a new ORCA\nmodel on RNA004 datasets and show that ORCA achieves high accuracy on both IVT and\nmouse brain RNA004 samples, demonstrating that ORCA generalizes robustly across\nsequencing chemistries with minimal adaptation.\n2. Clarified the rationale for training dataset selection. We expanded the discussion on\nusing IVT derived training datasets for model training and detailed how domain-adversarial\nlearning enable ORCA to learn generalized representation of modification presence, thereby\nmitigating dependence on specific synthetic modification types.\n3. Expanded benchmarking of modification type annotation using diverse NGS-based\ntechnologies. We incorporated extensive benchmarking analyses of modification type\nannotation based on different combination of RNA modification profiling technologies. These\nresults show that ORCA achieves consistently high accuracy across technologies, and that\nintegrating information from complementary assays significantly improves sensitivity without\nintroducing technology-specific bias.\nTogether, these improvements directly address the reviewer’s concerns regarding platform\ngeneralizability, training dataset selection, and potential bias from different NGS-based datasets.\nWe believe that the expanded analyses substantially strength the manuscript and further\nhighlight ORCA as a broadly applicable and reliable framework for transcriptome-wide RNA\nmodification analysis.\n1. All the results of the manuscript are based RNA002 chemistry, which is already discontinued,\nmaking the methodology in practice of very limited use. This is potentially the most significant\nflaw of the manuscript.\nResponse: Thanks for the reviewer for this insightful comment. We agree that limited\ndependence on RNA002 chemistry would reduce practical utility, and we have therefore trained\nan additional ORCA model using RNA004 datasets to demonstrate ORCA’s robustness across\ndifferent sequencing chemistries.\nConsidering that ORCA identifies RNA modification by leveraging systematic shifts in\nsignal and sequence features arising from mixed stoichiometry at modified sites, we first\nexamined whether these polymorphic features are preserved across RNA002 and RNA004\ndatasets. To this end, we collected RNA004 curlcake IVT datasets containing seven types of\nsynthetic RNA modifications together with unmodified controls (PMID: 40001217), and\nevaluated the distribution of signal- and base-level features between modified and unmodified\nsites using our simulation pipeline. Consistent with our observations in RNA002, RNA004 data\nshowed the similar partial-stoichiometry dependent feature shifts (Rebuttal Figure 1),\nindicating the feasibility of generalizing ORCA’s strategy to the RNA004 platform.\nRebuttal Figure 1. Signal- and base-level feature shifts across RNA002 and RNA004\nsynthetic datasets. (a)Boxplots showing two representative signal-level features on RNA002\nplatform: the 2% (top) and 98% (bottom) quantiles of normalized current levels across all reads\ncovering same site. (b) Distribution of per-site basecalling quality for the same set of\nmodification types. (c) Distribution of per-site mutation rates across modification types. (d-f)\nCorresponding analyses using RNA004 curlcake datasets, demonstrating that the similar shifts\nof signal- and base-level features are preserved across sequencing chemistries.\nTherefore, we further evaluated ORCA’s performance on RNA004 curlcake dataset\n(Rebuttal Figure. 2a). We trained a new RNA004 model using these datasets, which achieved\nhigh performance (AUROC > 0.93) for all seven modification types(Rebuttal Figure. 2b-c),\nsuggesting ORCA’s architecture is inherently compatible with RNA004 chemistry. To test\nwhether ORCA truly leverage chemistry-agnostic polymorphic features, we trained a mixed-\nchemistry model using both RNA002 and RNA004 datasets (Rebuttal Figure. 2d). This joint\nmodel generalized effectively across both platforms and achieved strong performance on the\ncorresponding RNA002 and RNA004 test datasets (AUC > 0.8, Rebuttal Figure. 2e-f). These\nresults confirm that ORCA signal- and base-level polymorphisms features that are not restricted\nto a specific sequencing chemistry.\nRebuttal Figure 2. Evaluation of ORCA on RNA004 chemistry. (a) Schematic of the\nRNA004 specific model training workflow. (b-c) PR and ROC curve for RNA004 model\nevaluated on RNA004 test set. (a) Schematic of the cross-chemistry training workflow. (e-f)\nPR and ROC from cross-chemistry validation of the mixed-chemistry model, demonstrating\neffective generalization across both RNA002 and RNA004 chemistries.\nFurthermore, to evaluate ORCA’s de novo modification prediction performance on\nRNA004 chemistry, we conducted a zero-shot prediction analysis using the RNA004 curlcake\ndatasets. For each of the seven modification types, we trained an ORCA model with that\nmodification excluded from the training set and then evaluated its zero-shot prediction\nperformance. ORCA achieved robust de novo detection of unseen modifications, demonstrating\nstrong generalization to novel modifications under RNA004 chemistry (Rebuttal Figure. 3a-\nb). We further generated an additional mouse brain direct RNA-seq dataset using RNA004\nchemistry, and compared m6A sites detected on the RNA002 and RNA004 platforms using\ncorresponding non-m6A ORCA models (Rebuttal Figure. 3c). As shown in Rebuttal Figure.\n3d, ORCA exhibited consistent m6A detection performance across both chemistries, and the\npredicted m6A signals were highly reproducible between RNA002 and RNA004 datasets\n(Rebuttal Figure. 3e), confirming the accurate and robust de novo prediction performance on\nRNA004 datasets.\nRebuttal Figure 3. Evaluation of ORCA’s performance in de novo modification detection\nacross different platform. (a-b) PR and ROC curve for leave-one-out models evaluated on\ntest set of unseen modification type. (c) Experimental workflow of m6A detection evaluation\nusing non-m6A RNA002 and RNA004 models. (d) modScore of m6A sites detected by m6-\nSAC-seq and background DRACH motifs on RNA002 (left, P < 3*10-220, Wilcoxon rank-sum\ntest) and RNA004 datasets (right, P < 1*10-31, Wilcoxon rank-sum test). (e) Venn diagram\nshowing the overlap of ORCA-detected m6A sites between the RNA002 and RNA004 datasets.\nBesides, although RNA002 has been discontinued, the substantial archive of existing\nRNA002 datasets remains valuable for further integrative analysis. Public resources such as\nhuman poly(A) transcriptome maps (PMID: 31740818), Arabidopsis epitranscriptome profiles\n(PMID: 31931956), and the SG-NEx project (PMID: 40082608) provide extensive datasets that\ncan benefit from improved analytical frameworks. Thus, ORCA offers a powerful tool for\nintegrative analysis of these valuable resources. We have included discussion of these points in\nlines 99-105 and 529-541 of the revised manuscript.\n2. There are quite a number widely occuring RNA modifications. The method accomodated\nonly six of them, which is quite weird to me. It seems straightforward that the method should\nbe trained on more RNA modification types.\nResponse: Thanks for the comment. We agree with the reviewer that including a broader set\nof RNA modification types into the training dataset would further improve model performance.\nHowever, ORCA requires high-quality training sets with accurate ground-truth labels for both\nmodification presence and stoichiometry. Most endogenous direct RNA-seq datasets cannot\nprovide such stoichiometry labels, which cannot be reliably determined using current NGS-\nbased detection methods. Therefore, many widely occurring RNA modifications cannot be\ndirectly used as training data for ORCA.\nTo address this, we employed IVT-derived datasets in which modification types are\nexperimentally controlled and used a simulation pipeline to generate training samples with\ndefined stoichiometry by mixing modified and unmodified reads. This approach inherently\nlimits the modification types available for training, as only a limited set (m6A, m5C, m1A,\nhm5C, 5FC, Ψ in RNA002; m6A, m5C, m1U, hm5C, m1Ψ, Ψ, ac4C in RNA004) of\nmodifications is currently supported by high-quality IVT datasets.\nDespite this limitation, ORCA’s domain adversarial training framework enables the model\nto learn generalized features across different modification types. Thus, ORCA achieves robust\nperformance on unseen modification types, including m6A, m5C, m1A, 5fC, hm5C, Inosine,\nNm, m6Am, m7G, m1acp3Ψ, m62A, m3U, Ψm, ac4C, m1Ψ (Rebuttal Table 1), indicating\nthat the limited diversity in the training set does not substantially impair the model’s\ngeneralization ability. As additional IVT transcriptome resources become available, ORCA can\nalso be readily expanded to incorporate more modification types, further enhancing its\nperformance. We have added the discussion of these points in lines 99-105 and 529-541 of the\nrevised manuscript.\nRebuttal Table 1. RNA modification types used to evaluate ORCA’s de novo modification\ndetection performance.\nChemistry Modifications Sample\nm6A, m5C, m1A, 5fC, hm5C, Ψ, Nm ELIGOS IVT\nm1acp3Ψ, m62A, m3U, Ψm, m7G, Nm Ribosomal RNA\nRNA002\nInosine Yeast\nm6A Mouse brain\nRNA004 m6A, m5C, 5hmC, ac4C, m5U, Ψ, m1Ψ Curlcake IVT\nm6A Mouse brain\n3. There are quite a few studies reporting the technical bias of different epitranscirptome\nprofiling approaches. The consistency of epitranscriptome datasets generated by two different\ntechnologies is usually less than 10%. The dataset stored in RMBase, DirectRMDB and etc all\nsuffer from low-quality epitranscriptome profiling datasets. Methodology trained on different\ndatasets can behave very differently.\nResponse: Thanks for the thoughtful comment. We systematically evaluated potential batch\neffects across different RNA modification profiling approaches using four major modification\ntypes supported by multiple NGS-based assays in RMBase v3.0 and DirectRMDB. As the\nreviewer suggested, cross-technology overlaps for m6A, m5C, Nm and m1A were generally\nlimited (Rebuttal Figure. 4a-b). To assess the impact of this inconsistency, we trained\nORCA’s annotation model using sites derived from different combinations of technologies and\nevaluated its performance on sites detected by other assays.\nRebuttal Figure 4. Evaluation of cross-technology annotation across four RNA\nmodifications. (a) Heatmap showing pairwise intersections of different m6A technologies in\nRMBase v3.0. Values indicate the number of overlapping sites. (b) Intersections among\ntechnologies for m5C, m1A and Nm in RMBase v3.0. (c) Recall (left) and 1-FDR (right) of de\nnovo m6A discovery when training different combination of technologies (rows) and validating\non different dataset (columns). (d) Recall and 1-FDR of de novo discovery for m5C, and Nm.\n(e) Recall and 1-FDR of m1A discovery when training on specific technologies and evaluating\non held-out validation set.\nFor m6A, although overlap between methods was small, models trained on individual\nassays achieved consistently high recall within orthogonal methods such as MeRIP-seq, m6A-\nCLIP-seq, miCLIP and m6A-seq (Rebuttal Figure. 4c-d). Specifically, recall decreased when\nmodels were trained solely on antibody-based IP/CLIP datasets and evaluated on antibody-\nindependent assays (for example, m6A-seq on DART-seq and MAZTER-seq). However, recall\nwas largely improved once complementary antibody-independent datasets were incorporated\ninto the training set, highlighting the importance of combining multiple resources in our\nintegrative annotation framework. Importantly, FDR remained low across all models (Rebuttal\nFigure. 4d), indicating the high reliability of ORCA’s annotation strategy. Similar trends were\nobserved for m1A, while cross-technology prediction performance for Nm and m5C remained\nconsistently high (Rebuttal Figure. 4e).\nTogether, these results demonstrate that limited overlap across technologies does not imply\nunreliability, instead modification-associated features are reproducible across methods.\nIntegrating multiple assays from public databases effectively mitigates technology-specific\nbiases and improves sensitivity without sacrificing precision. Furthermore, ORCA can be\nfurther extended to incorporate emerging NGS-based RNA modification profiling approaches,\nunderscoring the flexibility and expandability of our annotation framework. We have included\ndiscussion of these results in lines 391-402 of the revised manuscript.\nReviewer #1 (Remarks on code availability):\nVery decent, with clear instruction\nResponse: Thanks for the reviewer’s comment.\nReviewer #2 (Remarks to the Author):\nA novel tool, ORCA, is presented for the identification of RNA modifications in data generated\nthrough the direct sequencing of native RNA. This is a very active area of research. Around\n170 modifications exist, and tools have been developed that are able to identify only a small\nsubset of these modifications. New tools are then needed for a more comprehensive mapping\nand characterization of the marks that make up the epitranscriptome.\nORCA represents a promising solution for the comprehensive mapping of the marks that are\nvisible through the lens of Nanopore sequencing. It leverages deep learning methods aimed at\ncollectively identifying these marks – through adversarial training, i.e. forcing the learning of\ngeneralized features that are shared across modifications - and takes advantage of large\ncollections of RNA modification positions stored in public databases for annotating the\nidentified hits. An extensive set of validations through orthogonal methods and datasets is\nincluded.\nAltogether, the data suggest that ORCA could represent an important new method in this\nrapidly evolving while competitive field. The manuscript is well written, and results are well\npresented and convincing. Major and minor comments are listed below.\nResponse: We sincerely thank the reviewer for the thorough evaluation and the positive\nfeedback on our work.\nMAJOR COMMENTS (1-13)\n1. TandemMod is one of the few tools currently available that aims at comprehensively\nidentifying RNA modifications. I would suggest a more thorough discussion of how ORCA's\nfeatures stand against it (including TandemMod model transfer abilities), and a more extensive\ncomparison of the results, vs some independent reference “gold standard”. Such comparison,\nvs TandemMod, is currently largely limited to m6A predictions.\nResponse: Thanks for the comment. We have added a comprehensive comparison between\nTandemMod and ORCA in the revised manuscript as follows (lines 515-527): “TandemMod\nemploys deep-learning models to identify multiple RNA modifications (including m6A, m1A\nand m5C) at the single-read level, and further incorporate transfer learning to predict\nadditional modification such as m7G, hm5C and pseudouridine using limited training examples.\nWhile this approach enables simultaneous detection of multiple RNA modification types, it still\nrelies on IVT-derived training sites, which restricts its ability to capture the full RNA\nmodification landscape. In contrast, ORCA leverages a domain-adversarial learning strategy\nto infer modification presence based on signal polymorphism, enabling the detection of a wide\nrange of RNA modification types without requiring corresponding IVT training sets. However,\nthis approach requires sufficient read depth and is less effective at low-coverage sites (<10\nreads). Further work integrating both strategies may enable robust de novo modification\ndetection at single molecule level.”\nAs the reviewer suggested, we have also expanded the benchmarking among TandemMod,\nORCA and other tools for m5C and Ψ using paired WT/KO experiments or orthogonally\nvalidated methods. For m5C prediction in HeLa NSUN2 knockout data, CHEUI shows the\nhighest precision in identifying downregulated m5C sites, while ORCA achieves the best\nsensitivity, recovering the largest number of downregulated m5C sites with a balanced overall\nperformance (Rebuttal Figure. 5a-b). For pseudouridine detection in ribosomal RNA,\nTandemMod identifies the largest number of known Ψ sites (75/89), whereas ORCA achieves\ncomparable performance (73/89) with substantial overlap, demonstrating that both tools show\nstrong predictive power for Ψ detection.\nFurthermore, we evaluated the false discovery rate using an IVT transcriptome. Five tools\nincluding TandemMod, ORCA, EpiNano-SVM, NanoPSU and m6Anet were applied to the\nmodification-free human IVT transcriptome, and the cumulative distributions of prediction\nscores were calculated to measure false-discovery performance via AUC. Compared with\nTandemMod (AUC = 0.771), ORCA showed lower false discovery signal (AUC = 0.821),\nwhich was surpassed only by the m6A-specific tool m6Anet. These results indicate that ORCA\nprovides reliable modification detection despite targeting a more general modification spectrum.\nWe have included discussion of these results in lines 158-162, 195-200, 202-205 of the revised\nmanuscript.\nRebuttal Figure 5. Evaluation of m5C and Ψ prediction performance. (a) Proportion of\ndownregulated m5C sites among top differentially modified sites after NSUN2 knockout (KO).\n(b) Number of significantly up- and down-regulated m5C sites (Chi-squared test, P < 0.05)\ndetected by each tool. Red: upregulated; Blue: downregulated. (c) Venn diagram comparing Ψ\nsites prediction of on human rRNA among ORCA, TandemMod, NanoPSU, and Tombo. (d)\nCumulative distribution of modification probability of five modification detection tools.\n2. The concept of adversarial training, i.e. forcing the learning of generalized features that are\nshared across modifications is interesting. It remains to be clarified what is the level of control\nthat the Authors have on which/how many modifications are suitable for this approach\naccording to the collected features. Are the “shared features” necessarily shared among all\nmodifications under consideration, or they can be shared among sub-groups of a larger set of\nmodifications under consideration? Eventually, what is the search space, how many\nmodifications are eventually targeted, considering the ability of the tool to proceed without\ntraining data for specific marks?\nResponse: Thanks for the comment. ORCA utilizes the systematic shifts in signal- and base-\nlevel features between modified and unmodified positions and thus can be theoretically\nextended to most RNA modifications that cause detectable signal perturbations without being\nrestricted to specific modification type subgroups or sequencing chemistries. We have\ndemonstrated that ORCA successfully captured modification-presence across 15 modification\ntypes using only 6 IVT-derived training modifications (Rebuttal Table. 1). Although\nadditional endogenous RNA modifications may also be detectable by ORCA, many of these\ncannot be reliably validated with current NGS-based methods and were therefore not included\nfor further evaluation.\nRebuttal Table 1. RNA modification types used to evaluate ORCA’s de novo modification\ndetection performance.\nChemistry Modifications Sample\nm6A, m5C, m1A, 5fC, hm5C, Ψ, Nm ELIGOS IVT\nm1acp3Ψ, m62A, m3U, Ψm, m7G, Nm Ribosomal RNA\nRNA002\nInosine Yeast\nm6A Mouse brain\nRNA004 m6A, m5C, 5hmC, ac4C, m5U, Ψ, m1Ψ Curlcake IVT\nm6A Mouse brain\nWe also examined the consistency of these polymorphism features across RNA002 and\nRNA004 chemistries using our constructed training datasets. As shown in Rebuttal Figure. 6,\nthese raw features exhibited similar perturbation between modified and unmodified sites.\nSignificant shifts in signal distribution feature (Rebuttal Figure. 6a), mean basecalling quality\n(Rebuttal Figure. 6b), and mutation rate (Rebuttal Figure. 6c) were observed across\nmodification types, indicating that these polymorphism patterns are broadly shared. After\napplying adversarial training, the model’s embedding space further reduced type-specific\nfeature separation while preserving the shared modification presence related variance\n(Rebuttal Figure. 6d). These results demonstrate that ORCA effectively maps different\nmodification types into a common latent space, enabling generalized modification detection\nand facilitating downstream identification and annotation analyses.\nRebuttal Figure 6. Shared polymorphism features across modifications and adversarially\nlearned embeddings. (a) Boxplots of two site-level signal features: the minimum (top) and\nthe maximum (bottom) of the normalized read-level current across all reads covering a site.\nShown for six modifications (m1A, m6A, m5C, hm5C, Ψ, 5fC) and unmodified sites. Boxes\nshow the median and interquartile range, whiskers extend to 1.5× the interquartile range,\noutliers omitted. Two-sided Wilcoxon tests. Significance levels: * P<0.05, ** P<0.01, ***\nP<0.001. Base quality (b) and mutation rate (c) distributions for the same sets of sites among\nsix modification types. (d) t-SNE of raw features (left) and of embeddings after adversarial\ntraining (right).\n3. What is the impact of the ELIGOS training set on ORCA's performance, in terms of\nexperimental design (position, composition and prevalence of the modification) and quality of\nthe Nanopore data?\nResponse: Thanks for the comment. For the ELIGOS-derived training set, we chose to use only\n5-mers with modification at the central position to minimize potential biases introduced by\nneighboring modified bases. Since ORCA does not rely on specific sequence context, the\ntrained model can also capture strong signal shifts arising from modifications located in\ndifferent positions of the 5-mer (Figure 4c). These signal shifts also result in reductions in\nbasecalling quality (PMID: 32710622), which is also incorporated in ORCA’s model. We\nfurther demonstrated that ORCA generalizes to a wide range of 15 modification types\n(Rebuttal Table. 1), indicating that the composition of the training set does not impose a strong\nbias toward any specific modification type.\nNevertheless, the sequence compositions of the ELIGOS dataset represents only a small\nsubset (10,144/262,144) of all possible 9-mer contexts. To assess whether broader sequence\ndiversity could improve model robustness, we performed additional benchmarking using the\nIVT cDNA library (IVET) from the TandemMod study. The model trained on IVET\nexhibited noticeably better cross-dataset validation performance that the model trained on\nELIGOS (Rebuttal Figure. 7), indicating that greater sequence diversity in the training set\nimproves prediction performance across diverse sequence contexts.\nRegarding modification prevalence, each ELIGOS sample contains only a single\nmodification type, resulting in that no two modification types co-occur within the same read\nand could impact the model’s performance to predict co-occurring modifications in very close\nproximity.\nWe have included a comprehensive discussion of these limitations in lines 542-556 of the\nrevised manuscript.\nRebuttal Figure 7. Comparison between ORCA models trained on IVET and ELIGOS.\n(a) Precision-recall and ROC curves of model trained on IVET dataset and evaluated on\nELIGOS. (b) Precision-recall and ROC curves of the model trained on the ELIGOS dataset and\nevaluated on IVET datasets.\n4. What are ORCA's CPU-GPU/memory requirements?\nResponse: All measurements were performed on a single NVIDIA RTX 3090 GPU. Training\na six-modification model with 300,000 samples per modification for 300 epochs took about 23\nhours. Prediction on a K562 sample containing ~4,000,000 candidate sites took about 8 minutes.\nTraining the annotation model on 470,000 sites with a batch size of 512 required roughly 10\nminutes, and applying the annotation to a full dataset completed within 1 minute.\n5. The code for easily reproducing all analyses the analyses should be made available as an\nintegral part of the study.\nResponse: Thanks for the comment. We have made all code required to reproduce our analyses\npublicly available on Zenodo (doi: 10.5281/zenodo.17785932) as suggested.\nORCA employs a two-step strategy based on detection and annotation. The RNA modifications\ndatabase (DBs) that are considered are an important resource. However this rises a number of\npoints that should be clarified:\n6. These DBs are heavily populated by data about specific modifications (particularly m6A)\nand in specific cell lines (such as HEK and HeLa cells). How is the biased content of the\nconsidered modification DBs taken into account?\nResponse: Thanks for the comment. To mitigate the impact of imbalance in dataset size across\ndifferent modification types, we applied a dataset-balancing strategy in which all modification\ntypes were sampled to equal counts in each batch during training the annotation model. This\nprocedure reduces the model being biased toward the majority m6A class, and the trained\nannotation model exhibit low false-positive rate towards abundant m6A modification.\nWe also extensively evaluated potential cell-type specific biases using three major\nmodification types (m6A, m5C and Ψ). As shown in Rebuttal Figure. 8, when the training set\nfor a given modification type was restricted to data from only a single cell line, the resulting\nmodel still achieved considerable performance for predicting sites detected in other cell lines.\nAlthough a clear reduction in sensitivity was observed, the model maintained high precision\nand low FDR across unseen cell lines, indicating that cell-line composition affects recall but\ndoes not compromise annotation accuracy, and that ORCA generalizes robustly across different\nbiological backgrounds (lines 397-399 of the revised manuscript).\nRebuttal Figure 8. De novo discovery evaluation among cell lines. Subset numbers barplot\n(top), de novo discovery heatmap (middle) and precision/1-FDR barplot (right) of m6A\n(a), m5C (b) and Ψ (c).\n7. Using modifications sites non-associated to DBs hits (unannotated sites) as negative controls,\nin the first stage, seems like a problematic concept (especially for non-m6A modifications,\nsee the previous point), even though these sites are eventually re-considered, and potentially ex-\nnovo annotated. The Authors should better discuss at which extent they need to represent\nnegative controls, and what are the implications of this clearly over-simplified categorization.\nResponse: Thanks for the comment. We agree that unannotated sites cannot be treated as true\nnegative. They represent a heterogenous background that may contain false-positives\npredictions, unlabeled instances of known modifications, and potentially uncharacterized\nmodification types. Accordingly, ORCA incorporates an open-set recognition strategy to\nmaximize the separation between these background sites and annotated modification classes,\nthereby reducing false discoveries when annotating novel modification sites. In practice, we\nrandomly down-sampled background sites to approximately 5% of the total number of\ndatabase-supported training sites.\nAblation of this background-class modeling demonstrates that removing this component\nmarkedly increases the false discovery rate across all seven modification types (Rebuttal\nFigure. 9). In contrast, including the background class results in substantially more stringent\nperformance and a consistently lowers FDR, indicating that the open-set recognition strategy\neffectively filters misannotated or ambiguous sites. We have included these results in lines 366-\n369 of the revised manuscript.\nRebuttal Figure 9. Ablation of the background open-set recognition strategy. Cross-\nvalidation heatmaps showing the annotation accuracy of models trained with (a) or without (b)\nthe background control.\n8. How is the DB's composition impacting the results presented in Figure 6, about the clusters\nof adjacent modifications, which are heavily including m6A (being this over-represented in\nthose DBs).\nResponse: Thanks for the comment. As suggested in previous points, we employed both dataset\nrebalancing and a background open-set recognition strategy to minimize potential bias towards\nthe over-represented m6A modification. These steps ensure that the annotation and clustering\nresult is not disproportionately influenced by the predominant abundance of m6A in existing\ndatasets. Moreover, the higher proportion of m6A observed in the clusters is consistent with\nLC-MS/MS studies reporting that m6A is by far the most abundant endogenous mRNA\nmodification, accounting for approximately 0.304% of all adenosines in adult mouse neurons\n(PMID: 30048615), which is substantially higher than any other modification type. Taken\ntogether, these observations indicated that the composition does not strongly impact the\nclustering results, and the predominance of m6A reflects its underlying biological prevalence.\n9. More in general, how close ORCA's hits need to be to DB's records to be annotated by them?\nResponse: Thanks for the comment. Considering that the sequencing signal can be altered by\nmodified bases at any position within 5-mer, ORCA adopts a five-base distance threshold for\nmodification annotation. Specifically, for each database record we constructed 5 flanking\ntraining examples that encode both modification type and phase (the position of modified base\nwithin 5-mer), and the annotation model was trained to infer these two labels simultaneously.\nDuring the annotation stage, we retained only predictions in which the inferred modification\ntype matches the base at the predicted phase position to reduce false discoveries. Finally,\npredictions from overlapping 5-mers corresponding to the same underlying modification site\nare merged to avoid redundant calls.\nSequencing depth:\n10. How is ORCA impacted by the sequencing depth at a given genomic site?\nResponse: Thanks for the comment. We evaluated ORCA across different sequencing depths\nusing both IVT ELIGOS and mESC Mettl3 datasets. As shown in Rebuttal Figure. 10,\nAUROC and AUPRC moderately increase with read depth, which is consistent with that\npolymorphism features are more accurately estimated at sufficient depth. However, with low\nread depth of 10 reads, ORCA could still achieve high performance (AUPRC ~ 0.9) across six\nmodification types in ELIGOS, and achieves ~80% AUPRC and ~90% AUROC compared to\nhigh coverage sites in mESC datasets. Taken together, these results indicates that ORCA could\nbe applied to DRS data with a wide range of sequencing depth. We have included discussion\nof this limitation in lines 542-548.\nRebuttal Figure 10. Evaluation of impact of sequencing depth on ORCA performance. (a)\nBar plots of AUPRC (top) and AUROC (bottom) for six modifications in the ELIGOS IVT\ndatasets, with sites grouped into read-depth bins as indicated in the legend. (b) Precision–\nrecall curves (top) and ROC curves (bottom) for m6A detection in mESC Mettl3 WT and KO\nsamples, with sites grouped into the same depth bins; the legend reports AUPRC and AUROC\nfor each bin.\n11. When analyzing various samples or conditions, the number of reads could greatly vary\ndepending on technical reasons (e.g. the overall number of reads produced by flow cell) or for\nbiological reasons (e.g. varying number of reads associated to the same gene or genomic\nposition as a consequence of differential gene expression between conditions/treatments). Is\nORCA in need to / able to get rid of this bias, for example by sub-sampling the reads of the\nlarger dataset (at the sample and/or at the gene or site level) or by any other means?\nResponse: Thanks for the comment. As demonstrated in the previous comment, ORCA exhibits\nrobust performance with a minimal read-depth requirement of 10 reads and achieves consistent\ncoverage-independent performance for high coverage sites (>30 reads). Thus, sites below the\nminimal read depth threshold are excluded, but no additional down-sampling of high-coverage\nsites is required to avoid reducing performance of modification prediction.\nEvolution of the Nanopore platform:\n12. What is the compatibility of ORCA in terms of RNA002 and/or RNA004 Nanopore\nchemistry? Not proving that ORCA is able to work with RNA004 data would be a tremendous\nlimitation.\nResponse: Thanks for the reviewer for this insightful comment. We agree that limited\ndependence on RNA002 chemistry would reduce practical utility, and we have therefore trained\nan additional ORCA model using RNA004 datasets to demonstrate ORCA’s robustness across\ndifferent sequencing chemistries (lines 529-541).\nConsidering that ORCA identifies RNA modification by leveraging systematic shifts in\nsignal and sequence features arising from mixed stoichiometry at modified sites, we first\nexamined whether these polymorphic features are preserved across RNA002 and RNA004\ndatasets. To this end, we collected RNA004 curlcake IVT datasets containing seven types of\nsynthetic RNA modifications together with unmodified controls (PMID: 40001217), and\nevaluated the distribution of signal- and base-level features between modified and unmodified\nsites using our simulation pipeline. Consistent with our observations in RNA002, RNA004 data\nshowed the similar partial-stoichiometry dependent feature shifts (Rebuttal Figure 1),\nindicating the feasibility of generalizing ORCA’s strategy to the RNA004 platform.\nRebuttal Figure 2. Signal- and base-level feature shifts across RNA002 and RNA004\nsynthetic datasets. (a)Boxplots showing two representative signal-level features on RNA002\nplatform: the 2% (top) and 98% (bottom) quantiles of normalized current levels across all reads\ncovering same site. (b) Distribution of per-site basecalling quality for the same set of\nmodification types. (c) Distribution of per-site mutation rates across modification types. (d-f)\nCorresponding analyses using RNA004 curlcake datasets, demonstrating that the similar shifts\nof signal- and base-level features are preserved across sequencing chemistries.\nTherefore, we further evaluated ORCA’s performance on RNA004 curlcake dataset\n(Rebuttal Figure. 2a). We trained a new RNA004 model using these datasets, which achieved\nhigh performance (AUROC > 0.93) for all seven modification types(Rebuttal Figure. 2b-c),\nsuggesting ORCA’s architecture is inherently compatible with RNA004 chemistry. To test\nwhether ORCA truly leverage chemistry-agnostic polymorphic features, we trained a mixed-\nchemistry model using both RNA002 and RNA004 datasets (Rebuttal Figure. 2d). This joint\nmodel generalized effectively across both platforms and achieved strong performance on the\ncorresponding RNA002 and RNA004 test datasets (AUC > 0.8, Rebuttal Figure. 2e-f). These\nresults confirm that ORCA signal- and base-level polymorphisms features that are not restricted\nto a specific sequencing chemistry.\nRebuttal Figure 2. Evaluation of ORCA on RNA004 chemistry. (a) Schematic of the\nRNA004 specific model training workflow. (b-c) PR and ROC curve for RNA004 model\nevaluated on RNA004 test set. (a) Schematic of the cross-chemistry training workflow. (e-f)\nPR and ROC from cross-chemistry validation of the mixed-chemistry model, demonstrating\neffective generalization across both RNA002 and RNA004 chemistries.\nFurthermore, to evaluate ORCA’s de novo modification prediction performance on\nRNA004 chemistry, we conducted a zero-shot prediction analysis using the RNA004 curlcake\ndatasets. For each of the seven modification types, we trained an ORCA model with that\nmodification excluded from the training set and then evaluated its zero-shot prediction\nperformance. ORCA achieved robust de novo detection of unseen modifications, demonstrating\nstrong generalization to novel modifications under RNA004 chemistry (Rebuttal Figure. 3a-\nb). We further generated an additional mouse brain direct RNA-seq dataset using RNA004\nchemistry, and compared m6A sites detected on the RNA002 and RNA004 platforms using\ncorresponding non-m6A ORCA models (Rebuttal Figure. 3c). As shown in Rebuttal Figure.\n3d, ORCA exhibited consistent m6A detection performance across both chemistries, and the\npredicted m6A signals were highly reproducible between RNA002 and RNA004 datasets\n(Rebuttal Figure. 3e), confirming the accurate and robust de novo prediction performance on\nRNA004 datasets.\nRebuttal Figure 3. Evaluation of ORCA’s performance in de novo modification detection\nacross different platform. (a-b) PR and ROC curve for leave-one-out models evaluated on\ntest set of unseen modification type. (c) Experimental workflow of m6A detection evaluation\nusing non-m6A RNA002 and RNA004 models. (d) modScore of m6A sites detected by m6-\nSAC-seq and background DRACH motifs on RNA002 (left, P < 3*10-220, Wilcoxon rank-sum\ntest) and RNA004 datasets (right, P < 1*10-31, Wilcoxon rank-sum test). (e) Venn diagram\nshowing the overlap of ORCA-detected m6A sites between the RNA002 and RNA004 datasets.\nBesides, although RNA002 has been discontinued, the substantial archive of existing\nRNA002 datasets remains valuable for further integrative analysis. Public resources such as\nhuman poly(A) transcriptome maps (PMID: 31740818), Arabidopsis epitranscriptome profiles\n(PMID: 31931956), and the SG-NEx project (PMID: 40082608) provide extensive datasets that\ncan benefit from improved analytical frameworks. Thus, ORCA offers a powerful tool for\nintegrative analysis of these valuable resources.\n13. Regarding RNA004, a comprehensive comparison with Dorado's predictions of RNA\nmodifications should be included.\nResponse: Thanks for the comment. We performed a systematic comparison between Dorado\nand ORCA using the RNA004 curlcake IVT datasets. Specifically, three dorado RNA\nmodification models (model1, model2 and model3) were used to basecall Ψ, m6A, and m5C\nmodification respectively. Read-level modification calls were then aggregated to site-level\nusing modkit, which is similar to ORCA’s site-level training set construction to ensure a fair\ncomparison. For benchmarking, the per-site modification fractions of Dorado were used as\ncontinuous scores, and sites with modification rate ≥ 0.1 were considered as positive\nreference. As shown in Rebuttal Figure. 14, Dorado achieved superior performance for m6A\n(AUPRC = 0.96), m5C (AUPRC = 0.94), whereas ORCA showed better performance for Ψ\n(AUPRC = 0.96). While Dorado have accurate read- and site-level predictions for three\nmodification types, ORCA has the advantage of being applicable to a broader set of\nmodification types including 5hmC, ac4C, m1Ψ, m5U which are not currently covered by\nDorado. We have added a discussion of this comparison in the revised manuscript (lines 529-\n534).\nRebuttal Figure 11. Comparison between ORCA and dorado. ROC curve of ORCA (a) and\ndorado (b) on modification prediction on RNA004 Chemistry data of IVT curlcake sample.\nMINOR COMMENTS (14-24)\n14. Defining ORCA an \"universal framework for RNA modification detection\" is likely an\noverstatement, given the vast prevalence of untested (and currently largely untestable)\nmodifications\nResponse: Thanks for the important comment. We’ve revised the statement to “effective\nframework for detecting a wide range of RNA modifications”, which more accurately reflect\nthe current scope of ORCA.\n15. It remains unclear whether reads-level predictions are returned by ORCA\nResponse: Thanks for the comment. We have added clarification that ORCA also outputs read-\nlevel modification predictions for neighboring modification pairs through the EM-based model,\nwhich enables assessment of read-level linkage or mutual exclusivity between modification\nevents (lines 703-705).\n16. How were the features selected (Fig. 1a)?\nResponse: Thanks for the comment. ORCA integrates both signal-level and sequence-level\nfeatures within a ±2 5-mer window surrounding each candidate site, as illustrated in Figure.\n1a. For the signal-level features, raw current events from the f5c eventalign output were\nstandardized at each position using the formula provided by f5c:\n𝑆𝑖𝑔𝑛𝑎𝑙_𝑚𝑒𝑎𝑛𝑙𝑒𝑣𝑒𝑙 −𝑆𝑖𝑔𝑛𝑎𝑙_𝑚𝑒𝑎𝑛𝑟𝑒𝑓\nFor every𝑆 g𝑡𝑎en𝑛o𝑑m𝑎𝑟ic𝑑 𝑖p𝑧o𝑒s𝑑it_i𝑙o𝑒n𝑣, 𝑒s𝑙t=andardized current levels from all aligned reads were\n𝑆𝑖𝑔𝑛𝑎𝑙_𝑠𝑡𝑑𝑣𝑟𝑒𝑓\naggregated and interpolated to a fixed-length vector of 50 values to ensure consistent\nrepresentation across sites.\nFor the sequence-level features, we extracted alignment and basecalling statistics from\nsamtools mpileup results. At each position, the insertion rate, deletion rate, mismatch rate, and\nthe mean, median and standard deviation of base quality across all aligned reads were calculated.\nFinally, for each candidate site, we concatenated signal-level and sequence-level features from\nthe central position and two upstream and downstream neighbours, yielding a 5×56 feature\nmatrix that serves as the input to the modification presence prediction model.\n17. ORCA focuses on a 9 nt window centered on each candidate site. What is the rationale,\njustification, and implications of this choice?\nResponse: Thanks for the comment. During nanopore direct RNA sequencing, the signal is\ndetermined by approximately 5 nucleotides passing through the pore (PMID: 31931956). A\nsingle RNA modification can therefore perturb the current profile across five consecutive 5-\nmers surrounding the modification site, which together span a 9-nt region. To capture this effect,\nORCA extract signal and alignment features from these 9-nt windows to each candidate site.\nFurthermore, Although the official pore model for RNA004 is parameterized as a 9-mer, several\nstudies have shown that the central five positions still dominate the current signal (PMID:\n40085000), so the 9-nt window used by ORCA remains a reasonable approximation for both\nRNA002 and RNA004 chemistries. In future work, the window could be extended, for example\nto 13 nt, to more fully exploit the extended 9-mer context of RNA004.\n18. What about non-natural modifications, such as 4sU?\nResponse: Thanks for the constructive comment. As the reviewer suggested, we added an\nevaluation of the non-natural 4sU modification using public direct RNA-seq data from 4sU-\nlabeled K562 cells (PMID: 31839405). In the comparison between 4sU-labeled and DMSO\ncontrol samples, a significant elevated modification level for U-containing 5-mers against non-\nU 5-mers was observed, suggesting that ORCA effectively detected 4sU incorporation at\nlabeled sites (Rebuttal Figure. 15a). Furthermore, the transcript-level 4sU load estimated by\nORCA in the 4sU pulldown DRS data exhibited a significant correlation (R=0.40, p<1.7×10-\n16) with transcript abundance measured by short-read 4sU pulldown RNA-seq (Rebuttal Fig.\n15b), further supporting that ORCA’s applicability to non-natural modifications. We have\nincluded these results in lines 283-290 of the revised manuscript.\nRebuttal Fig. 15. ORCA Prediction on non-natural modifications. (a) Distribution of\nΔmodScore in K562 across 5-mers with or without U. Boxes indicate median and interquartile\nrange; whiskers extend to 1.5× the interquartile range. (b) Scatter plot of transcript-level 4sU\nload in the K562 4sU pulldown nanopore DRS data, defined as the summed ORCA-estimated\n4sU counts per transcript, versus transcript abundance in the short-read 4sU pulldown data;\nboth axes are log10-transformed. The Spearman correlation coefficient and P value are shown\nin the lower tight.\n19. It is unclear how the two-step strategy is connected to the mixed stoichiometry issue (lines\n99-106)\nResponse: Thanks for the comment. We have revised the description to avoid further confusion:\n“Thus, ORCA employes these polymorphic features to detect the presence of RNA modifications\nacross the transcriptome. Afterwards, ORCA integrates prior knowledge from established RNA\nmodification databases for effective annotation of a wide range of modification types.” (lines\n102-105).\n20. A 2.25% fraction of tested IVT sites identified as positive could be a very large number of\nmodifications, in a dataset where none of these are expected. To which number of modifications\nis this percentage corresponding to?\nResponse: Thanks for the comment. In the IVT experiments, the 2.25% fraction corresponds\nto 170,708 candidate sites identified in the initial modification presence prediction step.\nHowever, these candidates were subsequently subjected to ORCA’s annotation model, and\nempirical analysis based on the K562 datasets indicates that fewer than 1% of them could be\nmatched to public resources or de novo annotated. Therefore, the final number of annotated\nfalse positive modifications in the IVT transcriptome should be fewer than 2,000, representing\nan extremely small proportion of all transcriptomic sites.\n21. The chosen color palette often makes it very difficult to distinguish the various colors, e.g.\nin Fig. 2e\nResponse: Thanks for the comment. We’ve revised the color palette in Figure. 2e and related\npanels to improve visual distinguishability and enhance clarity.\n22. The results about Ψ were not so convincing, due to the very little number of hits.\nResponse: Thanks for the comment. To ensure the reliability of evaluation results, only\nexperimentally validated Ψ modifications on rRNA were used for benchmarking the\nperformance of ORCA and other tools. Similar strategies were also employed in previous\nstudies (PMID: 34872593). Despite the limited number of sites, these high-confidence Ψ sites\ncould provide reliable benchmarking results.\n23. When comparing ORCA's predictions with m6A-SAC-seq (line 244) could the author better\ndescribed how many hits were called in total by ORCA and how many of these could be\nassociated with the DRACH motif\nResponse: Thanks for the comment. Here, a total of 1,489,442 sites were predicted as modified,\nwhere 162,259 (10.9%) were associated with DRACH motif. Among these candidate sites,\n109,187 overlapped with public m6A resources and 31,833 were newly annotated using\nORCA’s annotation pipeline.\n24. Unclear statements: “increased sequence variability” (line 101); \".. we restricted base-level\nfeatures to discrepancies between basecalled and reference sequences ..\" (lines 113); \"RNA\nmodifications without type-specific bias\" (118); \"LSTM layers\"; \"RNA modifications of the\nsame type\" (line 128).\nResponse: Thanks for the comment. We’ve revised these statements to “increased basecalling\nerrors”, “raw sequence or absolute current levels features were deliberately excluded to ensure\nthat these features could represent generalized modification status without sequence\npreference”, “To develop a generalized model capable of accurately predicting diverse RNA\nmodifications without being restricted to specific types”, “two LSTM layers that process the\nsequence in opposite directions”, “occurrence of the same RNA modification across\ndifferent transcriptomic positions” as suggested.\nReviewer #2 (Remarks on code availability):\nCode and documentation available on GitHub seem ok but not ideal.\nA more comprehensive documentation of the tools functionalities and how to chain them would\ncertainly benefit.\nResponse: Thanks for the comment. We’ve revised the documentation to better clarify\nORCA’s functionalities and provided detailed pipeline for main analyses in our manuscript.\nIt looks like the scripts used for performing the extensive evaluation and benchmarking\nanalyses presented in the manuscript is not available. I would consider this absolutely essential.\nResponse: Thanks for the comment. We’ve uploaded the codes for benchmarking and\ndownstream analysis to the Zenodo repository (doi: 10.5281/zenodo.17785932) as suggested.\nReviewer #1 (Remarks to the Author):\nGlad to confirm that all my comments are properly addressed. I am happy to recommend the\nmanuscript to be accepted at NC.\nResponse: We thank the reviewer for the positive feedback.\nReviewer #1 (Remarks on code availability):\nAll looks good with detailed instructions.\nResponse: Thanks for the positive and encouraging comments.\nReviewer #2 (Remarks to the Author):\nThe manuscript has been substantially revised. The rebuttal clearly outlines the novel analyses\nand results, convincingly addressing the raised points.\nResponse: We thank the reviewer for the positive comments.\nIn particular, the presented tool for RNA modifications detection, ORCA, has been\nconvincingly demonstrated to be compatible with both RNA002 and RNA004 Nanopore\ndRNA-seq data. This enables both the reanalysis of previously published data (or old in house\nunpublished data) and the analysis of novel datasets obtained with the latest RNA004 Nanopore\nchemistry.\nResponse: Thanks for the comments.\nAltogether this study represents an important step towards the ability to capture the whole\nepitranscriptome and ORCA is likely to position itself as an important tool in the field.\nResponse: Thanks for the comments.\nReviewer #2 (Remarks on code availability):\nI haven't further reviewed the code, compared to the first submission. However, I previously\nrequested the Authors to share scripts for enabling the repetition of all the extensive\nbioinformatic analyses presented in the study. A Zenodo link has been provided in the rebuttal\nand in the revised manuscript (doi: 10.5281/zenodo.17785932), but I failed to locate this\nfundamental resource.\nResponse: Thanks for carefully checking the code availability. We updated the Zenodo record\nat https://doi.org/10.5281/zenodo.17785932, where the scripts necessary to repeat the\nbioinformatic analyses are provided in the archive other_tools.zip."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Glad to confirm that all my comments are properly addressed. I am happy to recommend the manuscript to be accepted at\nNC.", "reviewer_2": "The manuscript has been substantially revised. The rebuttal clearly outlines the novel analyses and results, convincingly\naddressing the raised points.\nIn particular, the presented tool for RNA modifications detection, ORCA, has been convincingly demonstrated to be\ncompatible with both RNA002 and RNA004 Nanopore dRNA-seq data. This enables both the reanalysis of previously\npublished data (or old in house unpublished data) and the analysis of novel datasets obtained with the latest RNA004\nNanopore chemistry.\nAltogether this study represents an important step towards the ability to capture the whole epitranscriptome and ORCA is\nlikely to position itself as an important tool in the field."}, "author_response": "Response: We thank the reviewer for the positive comments.\nResponse: Thanks for the comments.\nResponse: Thanks for the comments."}]

Decision letter