Raw signal segmentation for estimating RNA modification from Nanopore direct RNA sequencing data

eLife Assessment

This study presents SegPore, a valuable new method for processing direct RNA nanopore sequencing data, which improves the segmentation of raw signals into individual bases and boosts the accuracy of modified base detection. The evidence presented to benchmark SegPore is solid, and the authors provide a fully documented implementation of the method. SegPore will be of particular interest to researchers studying RNA modifications.

https://doi.org/10.7554/eLife.104618.4.sa0

Significance of the findings:

Valuable: Findings that have theoretical or practical implications for a subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Solid: Methods, data and analyses broadly support the claims with only minor weaknesses

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Materials and methods
Results
Discussion
Appendix 1
Appendix 2
Data availability
References
Article and author information
Metrics

Abstract

Estimating RNA modifications from Nanopore direct RNA sequencing data is a critical task for the RNA research community. However, current computational methods often fail to deliver satisfactory results due to inaccurate segmentation of the raw signal. We have developed a new method, SegPore, which leverages a molecular jiggling translocation hypothesis to improve raw signal segmentation. SegPore is a pure white-box model with enhanced interpretability, significantly reducing structured noise in the raw signal. We demonstrate that SegPore outperforms state-of-the-art methods, such as Nanopolish and Tombo, in raw signal segmentation across three large benchmark datasets. Moreover, the improved signal segmentation achieved by SegPore enables SegPore+m6Anet to deliver state-of-the-art performance in site-level m6A identification. Additionally, SegPore surpasses baseline methods like CHEUI in single-molecule level m6A identification.

Introduction

RNA modifications play important roles in different diseases, such as Acute Myeloid Leukemia (Yankova et al., 2021) and Fragile X Syndrome (Prieto et al., 2020), as well as fundamental biological processes like cell differentiation (Bellodi et al., 2013; Lee et al., 2019) and immune response (Quin et al., 2021). To date, researchers have identified over 150 different types of RNA modifications (Boccaletto et al., 2022; Zimna et al., 2023; Ohira and Suzuki, 2024; Chen et al., 2023), highlighting the complexity and diversity of RNA regulation. These modifications are essential for proper RNA function, including maintaining secondary structures and facilitating accurate protein synthesis, such as the role of tRNA modifications in decoding mRNA codons (Agris et al., 2017).

The practical applications of RNA modifications are far-reaching. For instance, N1-methylpseudouridine (m1Ψ) has been used to enhance the efficacy of COVID-19 mRNA vaccines, highlighting their therapeutic potential (Nance and Meier, 2021). However, identifying and characterizing RNA modifications presents significant challenges due to the limitations of existing experimental techniques.

Traditional methods for RNA modification detection, such as MeRIP-Seq (Meyer et al., 2012), miCLIP (Linder et al., 2015), and m6ACE-Seq (Koh et al., 2019), rely on immunoprecipitation techniques that use modification-specific antibodies. While effective, these methods have several drawbacks. First, each method requires a separate antibody for each RNA modification, making it difficult to study multiple modifications simultaneously. Additionally, these techniques can only infer modification locations from next-generation sequencing (NGS) data, rather than measuring the modifications directly. As a result, they struggle to provide single-molecule resolution, limiting their accuracy and scope.

Recent advancements in direct RNA sequencing (DRS) by Oxford Nanopore Technologies (ONT) offer a promising alternative. DRS allows for the direct measurement of electrical currents as RNA molecules translocate through a nanopore, providing a potential avenue for detecting RNA modifications at the single-molecule level. Two versions of the direct RNA sequencing (DRS) kits are available: RNA002 and RNA004. Unless otherwise specified, this study focuses on RNA002 data. In this technique, five nucleotides (5-mers) reside in the nanopore at a time, and each 5-mer generates a characteristic current signal based on its unique sequence and chemical properties (Loman et al., 2015). By analyzing these signals, it is possible to infer the original RNA sequence and detect modifications like N6-methyladenosine (m6A).

The general workflow of Nanopore direct RNA sequencing (DRS) data analysis is as follows. First, the raw electrical signal from a read is basecalled using tools such as Guppy or Dorado (blawrence-ont and malton-ont, 2026), which produce the nucleotide sequence of the RNA molecule. However, these base-called sequences do not include the precise start and end positions of each ribonucleotide (or k-mer) in the signal. Because basecalling errors are common, the sequences are typically mapped to a reference genome or transcriptome using minimap2 (Li, 2018) to recover the correct reference sequence. Next, tools such as Nanopolish (Loman et al., 2015; Simpson et al., 2017) and Tombo (Stoiber et al., 2017) align the raw signal to the reference sequence to determine which portion of the signal corresponds to each k-mer. We define this process as the segmentation and alignment task (abbreviated as the segmentation task), which is referred to as ‘eventalign’ in Nanopolish. Based on this alignment, Nanopolish extracts various features—such as the start and end positions, mean, and standard deviation of the signal segment corresponding to a k-mer. This signal segment or its derived features is referred to as an ‘event’ in Nanopolish. The resulting events serve as input for downstream RNA modification detection tools such as m6Anet (Hendra et al., 2022) and CHEUI (Acera Mateos et al., 2024).

However, significant computational challenges remain. Segmenting the raw current signal into distinct 5-mers and distinguishing between normal nucleotides and their modified counterparts is a complex task. Current methods, such as Nanopolish, employ change-point detection methods to segment the signal and use dynamic programming methods and Hidden Markov Models (HMM) to align the derived segments to the reference sequence, but they are prone to noise and inaccuracies, which degrade performance in downstream tasks like RNA modification prediction. The root of this issue lies in the fact that these methods do not accurately model the physical process of Nanopore sequencing, particularly the dynamics of the motor protein that drives RNA through the pore. As a result, the segmentation process is not well-aligned with the actual translocation mechanics, leading to signal noise that hinders precise modification detection.

SegPore is a novel tool for direct RNA sequencing (DRS) signal segmentation and alignment, designed to overcome key limitations of existing approaches. By explicitly modeling motor protein dynamics during RNA translocation with a Hierarchical Hidden Markov Model (HHMM), SegPore segments the raw signal into small, biologically meaningful fragments, each corresponding to a k-mer sub-state, which substantially reduces noise and improves segmentation accuracy. After segmentation, these fragments are aligned to the reference sequence and concatenated into larger events, analogous to Nanopolish’s ‘eventalign’ output, which serve as the foundation for downstream analyses. Moreover, the ‘eventalign’ results produced by SegPore enhance interpretability in RNA modification estimation. While deep learning–based tools such as m6Anet classify RNA modifications using complex, non-transparent features (see Figure 4—figure supplement 1), SegPore employs a simple Gaussian Mixture Model (GMM) to distinguish modified from unmodified nucleotides based on baseline current levels. This transparent modeling approach improves confidence in the predictions and makes SegPore particularly well-suited for biological applications where interpretability is essential.

By introducing SegPore, we bridge the gaps left by traditional tools. SegPore utilizes a hierarchical hidden Markov model (HHMM) for more precise segmentation and combines it with signal alignment and a GMM for RNA modification prediction, offering both greater accuracy and interpretability. This integrated approach enables robust m6A modification detection at the single-molecule level, facilitating reliable, transparent predictions for both site-specific and molecule-wide m6A modifications.

Materials and methods

SegPore workflow

Workflow overview

Request a detailed protocol

An overview of SegPore is illustrated in Figure 1A, which outlines its five-step process. The output of Step 3 is the ‘events’, which is analogous to the output generated by the Nanopolish (v0.14.0) ‘eventalign’ command and can be used as input for downstream models such as m6Anet. An ‘event’ refers to a segment of the raw signal that is aligned to a specific k-mer on a read, along with its associated features such as start and end positions, mean current, standard deviation, and other relevant statistics. Step 4 allows for direct modification prediction at both the site and single-molecule levels. Notably, a key feature of SegPore is the k-mer parameter table, which defines the mean and standard deviation for each k-mer in either an unmodified or modified state. During the training phase, Steps 3~5 are iterated multiple times to stabilize the parameters in the k-mer table, which are subsequently fixed and applied for modification prediction on the test data. A detailed description of the model is provided in Appendix 1. Unless otherwise noted, the following analysis focuses on RNA002 data, using 5-mers rather than k-mers.

Figure 1

Download asset Open asset

SegPore workflow.

(A) General workflow. The workflow consists of five steps: (1) First, raw current signals are basecalled and mapped using Guppy and Minimap2. The raw current signal fragments are paired with the corresponding reference RNA sequence fragments using Nanopolish. (2) Next, the raw current signal of each read is segmented using a hierarchical hidden Markov model (HHMM), which provides an estimated mean ( $μ_{i}$ ) for each segment. (3) These segments are then aligned with the 5-mer list of the reference sequence fragment using a full/partial alignment algorithm, based on a 5-mer parameter table. For example, $A_{j}$ denotes the base ‘ $A$ ‘ at the $j$ -th position on the reference. In this example, $A_{1}$ and $A_{2}$ refer to the first and second occurrences of ‘ $A$ ’ in the reference sequence, respectively. Accordingly, $μ_{1}$ and $μ_{2}$ are aligned to $A_{1}$ , while $μ_{3}$ is aligned to $A_{2}$ (4). All signals aligned to the same 5-mer across different genomic locations are pooled together, and a two-component Gaussian Mixture Model (GMM) is used to predict the modification at the site-level or single-molecule level. One component of the GMM represents the unmodified state, while the other represents the modified state. (5) GMM is used to re-estimate the 5-mer parameter table. (B) Hierarchical hidden Markov model (HHMM). The outer HMM segments the current signal into alternating base and transition blocks. The inner HMM approximates the emission probability of a base block by considering neighboring 5-mers. A linear model is used to approximate the emission probability of a transition block. (C) Full/partial alignment algorithms. Rows represent the estimated means of base blocks from the HHMM, and columns represent the 5-mers of the reference sequence. Each 5-mer can be aligned with multiple estimated means from the current signal. (D) Gaussian mixture model (GMM) for estimating modification states. The GMM consists of two components: the green component models the unmodified state of a 5-mer, and the blue component models the modified state. Each component is described by three parameters: mean (μ), standard deviation ( $σ$ ), and weight ( $ω$ ).

Preprocessing

Request a detailed protocol

We begin by performing basecalling on the input fast5 file using Guppy (v6.0.1), which converts the raw signal data into ribonucleotide sequences. Next, we align the basecalled sequences to the reference genome using Minimap2 (v2.28-r1209) (Li, 2018), generating a mapping between the reads and the reference sequences. Nanopolish provides two independent commands: ‘polya’ and ‘eventalign’. The ‘polya’ command identifies the adapter, poly(A) tail, and transcript region in the raw signal, which we refer to as the poly(A) detection results. The raw signal segment corresponding to the poly(A) tail is used to standardize the raw signal for each read. The ‘eventalign’ command aligns the raw signal to a reference sequence, assigning a signal segment to individual k-mers in the reference. It also computes summary statistics (e.g. mean, standard deviation) from the signal segment for each k-mer. Each k-mer together with its corresponding signal features is termed an event. These event features are then passed into downstream tools such as m6Anet and CHEUI for RNA modification detection. For full transcriptome analysis (Figure 3), we extract the aligned raw signal segment and reference sequence segment from Nanopolish’s events for each read by using the first and last events as start and end points. For in vitro transcription (IVT) data with a known reference sequence (Figure 4), we extract the raw signal segment corresponding to the transcript region for each input read based on Nanopolish’s poly(A) detection results.

Due to inherent variability between nanopores in the sequencing device, the baseline levels and standard deviations of k-mer signals can differ across reads, even for the same transcript. To standardize the signal for downstream analyses, we extract the raw current signal segments corresponding to the poly(A) tail of each read. Since the poly(A) tail provides a stable reference, we standardize the raw current signals for each read, ensuring that the mean and standard deviation are consistent across the poly(A) tail region. This step is crucial for reducing variability between different reads and ensuring more accurate signal segmentation and modification prediction in subsequent steps. See Section 3 of Appendix 1 for more details.

Signal segmentation via hierarchical Hidden Markov model

Request a detailed protocol

The RNA translocation hypothesis (see details in the first section of Results) naturally leads to the use of a hierarchical Hidden Markov Model (HHMM) to segment the raw current signal. As shown in Figure 1B, the HHMM consists of two layers. The outer HMM divides the raw signal into alternating base and transition blocks, represented by hidden states ‘B’ (base) and ‘T’ (transition). Within each base block, the inner HMM models the current signal at a more granular level.

The inner HMM includes four hidden states: ‘prev’, ‘next’, ‘curr’, and ‘noise’. These correspond to the previous, next, and current 5-mer in the pore, while ‘noise’ refers to random fluctuations. Each raw current measurement is emitted from one of these hidden states, providing a detailed model of the signal within the base blocks. A linear model with a large absolute slope is used to represent sharp changes in the transition blocks.

To segment the signal, we first model the likelihood of the HHMM. Given the raw current signal $y$ of a read, the hidden states of the outer hidden HMM are denoted by $g$ . $y$ and $g$ are divided into $2 K + 1$ blocks $c$ , where $y^{(k)}$ , $g^{(k)}$ correspond to kth block and $c = (c_{1}, c_{2}, \dots, c_{k}, \dots, c_{2 K + 1})$ , $c_{k} \in {B, T}$ . Blocks with odd indices $(k = 1, 3, 5, \dots, 2 K + 1)$ are base blocks, while those with even indices $(k = 2, 4, \dots, 2 K)$ are transition blocks. The likelihood of the HHMM is given by

p (y, g) = p (y | g) p (g)

= p (y | g) {π_{g_{1}}^{o u t e r} \prod_{i = 2}^{N} T_{g_{i - 1} g_{i}}^{o u t e r}}

= {\prod_{i = 1}^{N} p (y_{i} | g_{i})} {π_{g_{1}}^{o u t e r} \prod_{i = 2}^{N} T_{g_{i - 1} g_{i}}^{o u t e r}}

= {\prod_{k = 0}^{K} p (y^{(2 k + 1)} | c_{2 k + 1} = B)} {\prod_{k = 1}^{K} p (y^{(2 k)} | c_{2 k} = T)} {π_{g_{1}}^{o u t e r} \prod_{i = 2}^{N} T_{g_{i - 1} g_{i}}^{o u t e r}}

where $T_{g_{i - 1} g_{i}}^{o u t e r}$ is the transition matrix of the outer HMM and $π_{g_{1}}^{o u t e r}$ is the probability for the first hidden state. The first part of Equation 2 is emission probabilities, and the second part is the transition probabilities. It is not possible to directly compute the emission probabilities of the outer HMM (Equation 3) since there exist dependencies for the current signal measurements within a base or transition block. Therefore, we use the inner HMM and linear model (Equation 4) to handle the dependencies and approximate emission probabilities.

The inner HMM models transition between the ‘prev’, ‘next’, ‘curr’, and ‘noise’ states. For the ‘prev’, ‘next’, and ‘curr’ states, a Gaussian distribution is used to model the emission probabilities, while a uniform distribution models the ‘noise’ state. The forward-backward algorithm is used to compute the marginal likelihood of the inner HMM, which approximates the emission probabilities for base blocks. For transition blocks, a standard linear model computes the emission probabilities.

For any given $g$ , we can calculate the joint likelihood (Equation 1). We enumerate different configurations and select the one with the highest likelihood. This process segments the raw current signal into alternating base and transition blocks, where one or more base blocks may correspond to a single 5-mer. Each base block is characterized by its mean and standard deviation, which are used as input for downstream alignment tasks.

Due to the large size of fast5 files (which can reach terabytes), parameter inference for this model is computationally intensive. To address this, we developed a GPU-based inference algorithm that significantly accelerates the process. More details on this algorithm can be found in Appendix 2.

In summary, the HHMM allows us to accurately segment the raw current signal, providing estimates of the mean and variance for each base block, which are crucial for downstream analyses such as RNA modification prediction.

5-mer and 9-mer parameter table

Request a detailed protocol

We downloaded the k-mer models ‘r9.4_180 mv_70bps_5-mer_RNA’ from the ONT GitHub repository (here; Brennen, 2023) for RNA002 data. The columns labeled ‘level_mean’ and ‘level_stdv’ in these models were used as the mean and standard deviation values for the unmodified 5-mers. These values serve as the initial parameters in the 5-mer parameter table for SegPore, which we refer to as the ONT 5-mer parameter table. In the RNA004 data analysis, we obtained the 9-mer parameter table from the source code of f5c (version 1.5). Specifically, we used the array named ‘rna004_130bps_u_to_t_rna_9mer_template_model_builtin_data’ from the file here (accessed on 17 October 2025).

The initialization of the k-mer parameter table is a critical step in SegPore’s workflow. By leveraging ONT’s established k-mer models, we ensure that the initial estimates for unmodified k-mers are grounded in empirical data. These initial estimates are refined through SegPore’s iterative parameter estimation process, enabling the model to accurately differentiate between modified and unmodified k-mers during segmentation and modification prediction tasks. The refined k-mer parameters also provide a foundation for downstream analysis, such as alignment and m6A identification.

Alignment of raw signal segment with reference sequence

Request a detailed protocol

After segmenting the raw current signal of a read into base and transition blocks using HHMM, we align the means of base blocks with the 5-mer list of the reference sequence.

The alignment process involves three main cases:

Base block matching with a 5-mer: In this case, a base block aligns with a 5-mer, which is modeled by a Gaussian distribution. The 5-mer may exist in either an unmodified or modified state, and the corresponding Gaussian parameters are retrieved from the 5-mer parameter table.
Base block matching with an insertion: Here, the base block aligns with an indel (‘-’), indicating an inserted nucleotide in the read.
Deletion in the read: In this case, an indel (0.0) matches with a 5-mer, representing a deleted nucleotide in the read.

The alignment score function models these matching cases as follows. For the first case (a base block matching a 5-mer), we calculate the probability that the base block mean is sampled from either the unmodified or modified 5-mer’s Gaussian distribution. The larger of the two probabilities is used as the match score. For the second and third cases, the alignment is treated as noise, and a fixed uniform distribution is used to calculate the match score.

An important distinction from classical global alignment algorithms (Needleman–Wunsch algorithm) is that one or multiple base blocks may align with a single 5-mer. Given the base block means $μ = (μ_{1}, μ_{2}, \dots, μ_{i}, \dots, μ_{m})$ and 5-mer list $s = (s_{1}, s_{2}, \dots, s_{j}, \dots, s_{n})$ , we define $(m + 1) \times (n + 1)$ the score matrix as $M$ (Figure 1C). The first row and column of the matrix are reserved for indels ("0.0" and "-"), representing insertions or deletions in the base blocks or 5-mers, respectively.

We denote the score function by $f$ . The recursion formula of the dynamic programming alignment algorithm is given by

M (i, j) = m a x {\begin{cases} M (i - 1, j - 1) + f (μ_{i}, s_{j}) \\ M (i, j - 1) + f (0.0, s_{j}) \\ M (i - 1, j) + f (μ_{i}, -) \\ M (i - 1, j) + f (μ_{i}, s_{j}) \end{cases},

where $f$ is the score function. It can be seen that we can still align $μ_{i}$ with $s_{j}$ after we have aligned $μ_{i - 1}$ with $s_{j}$ , which fulfills the special consideration that one or multiple base blocks might be aligned with one 5-mer.

We implement two types of alignment algorithms (Figure 1C) based on the score matrix $M$ :

Full Alignment Algorithm: This algorithm aligns the full list of base block means with the full list of 5-mers, similar to classical global alignment. It traces back from the $(m + 1, n + 1)$ position in the score matrix.
Partial Alignment Algorithm: This aligns the full list of base blocks with the initial part of the 5-mer list, with no indels allowed in the base block means or the 5-mer list. The trace back starts from the maximum value in the last row of the score matrix.

A detailed description of both alignment algorithms is provided in Appendix 1. The output of the alignment algorithm is an alignment that pairs the base blocks with the 5-mers from the reference sequence for each read (Figure 1C). Base blocks aligned to the same 5-mer are concatenated into a single raw signal segment (referred to as an ‘event’), from which various features—such as start and end positions, mean current, and standard deviation—are extracted. A detailed derivation of the mean and standard deviation is provided in Section 5.3 in Appendix 1. In the remainder of this paper, we refer to these resulting events as the output of eventalign analysis, which also represents the final output of the segmentation and alignment task.

Modification prediction

Request a detailed protocol

After obtaining the eventalign results, we estimate the modification state of each motif using the 5-mer parameter table. Specifically, for each 5-mer, we compare the probability that its mean is sampled from either the modified or unmodified 5-mer’s Gaussian distribution. If the probability under the modified 5-mer distribution is higher, the 5-mer is predicted to be in the modification state, and vice versa for the unmodified state.

To estimate the overall modification rate at a specific genomic location, we pool all reads that map to that location on the reference sequence. The modification rate is calculated as the proportion of reads that are predicted to be in the modification state at that location. A detailed description of the modification state prediction process can be found in Appendix 1.

GMM for 5-mer parameter table re-estimation

Request a detailed protocol

To improve alignment accuracy and enhance modification predictions, we use a Gaussian Mixture Model (GMM) to iteratively fine-tune the 5-mer parameter table. As illustrated in Figure 1A, the rows of the 5-mer parameter table represent the 5-mers, while the columns provide the mean and standard deviation for both the unmodified and modified states. We denote a 5-mer by $s$ , with its relevant parameters as $μ_{s, u n}$ , $δ_{s, u n}$ , $μ_{s, m o d}$ , $δ_{s, m o d}$ .

Using the alignment results from all reads, we collect the base block means aligned to the same 5-mer (denoted by $s$ ) across different reads and genomic locations, particularly those with high modification rates. A two-component GMM is then fit to these base blocks. In this process, the mean of the first component is fixed to $μ_{s, u n}$ . From the GMM, we update the parameters $δ_{s, u n}$ , $ω_{s, u n}$ , $μ_{s, m o d}$ , $δ_{s, m o d}$ , and $ω_{s, m o d}$ , where $ω_{s, u n}$ , $ω_{s, m o d}$ represent the weights for the unmodified and modified components, respectively. Afterward, we manually adjust the 5-mer parameter table using heuristics to ensure that the modified 5-mer distribution is significantly distinct from the unmodified distribution (see details in Section 7 of Appendix 1).

This re-estimation process is only performed on the training data. The initial 5-mer parameter table is based on the table provided by ONT. After each iteration of updating the table, we rerun the SegPore workflow from the alignment step onward. Typically, the process is repeated three to five times until the 5-mer parameter table stabilizes (when the average change of mean values of all 5-mers is less than 5e-3). Once a stabilized 5-mer parameter table is estimated from the training data, it is used for RNA modification estimation in the test data without further updates. A more detailed description of the GMM re-estimation process is provided in Section 6 of Appendix 1.

m6A site level benchmark

Request a detailed protocol

The HEK293T wild-type (WT) samples were downloaded from the ENA database (accession number PRJEB40872), while the HCT116 samples were obtained from ENA under accession PRJEB44348. The reference sequence (Homo_sapiens.GRCh38.cdna.ncrna_wtChrIs_modified.fa) was downloaded from https://doi.org/10.5281/zenodo.4587661. Ground truth data were sourced from Supplementary Data 1 of Pratanwanich et al., 2021. Fast5 files for the test dataset (mESC WT samples, mESCs_Mettl3_WT_fast5.tar.gz; Jenjaroenpun et al., 2021) were retrieved from the NCBI Sequence Read Archive (SRA) under accession SRP166020.

During training, we initialized the 5-mer parameter table using ONT’s data. The standard SegPore workflow was executed on the training data (HEK293T WT samples), using the full alignment algorithm. The 5-mer parameter table estimation was iterated five times. For mapping, reads were first aligned to the cDNA and subsequently converted to genomic locations using Ensembl GTF file (GRCh38, v9). The same 5-mer across different genomic locations was pooled together. A 5-mer was considered significantly modified if its read coverage exceeded 1500 and the distance between the means of the two Gaussian components in the GMM was greater than 5 picoamperes (pA). As a result, modification parameters were specified for ten significant 5-mers, as illustrated in Figure 3—figure supplement 1A.

With the estimated 5-mer parameter table from the training data, we then ran the SegPore workflow on the test data. The transcript sequences from GENCODE release version M18 were used as the reference sequence for mapping, with the corresponding GTF file (gencode.vM18.chr_patch_hapl_scaff.annotation.gtf) downloaded from GENCODE to convert transcript locations into genomic coordinates. It is important to note that the 5-mer parameter table was not re-estimated for the test data. Instead, modification states for each read were directly estimated using the fixed 5-mer parameter table. Due to the differences between human (Figure 3—figure supplement 1A) and mouse (Figure 3—figure supplement 1B), only six 5-mers were found to have m6A annotations in the test data’s ground truth (Figure 3—figure supplement 1C). For a genomic location to be identified as a true m6A modification site, it had to correspond to one of these six common 5-mers and have a read coverage of greater than 20. SegPore derived the ROC and PR curves for benchmarking based on the modification rate at each genomic location.

In the SegPore+m6Anet analysis, we fine-tuned the m6Anet model using SegPore’s eventalign results to demonstrate improved m6A identification. We started with the pre-trained m6Anet model (available at here; Hendra, 2025, model version: HCT116_RNA002) and fine-tuned it using the eventalign results from HCT116 samples. SegPore’s eventalign output provided the pairing between each genomic location and its corresponding raw signal segment, which allowed us to extract normalized features such as the normalized mean $μ_{i}$ , standard deviation $σ_{i}$ , dwell time $l_{i}$ (number of data points in the event). For genomic location $i$ , m6Anet extracts a feature vector $x_{i} = {μ_{i - 1}, σ_{i - 1}, l_{i - 1}, μ_{i}, σ_{i}, l_{i}, μ_{i + 1}, σ_{i + 1}, l_{i + 1}}$ , which was used as input for m6Anet. Feature vectors from a randomly selected 80% of the genomic locations were used for training, while the remaining 20% were set aside for validation. We ran 100 epochs during fine-tuning, selecting the model that performed best on the validation set.

Ground truth data and the performance of other methods (Tombo v1.5.1, Nanom6A v2.0, m6Anet v1.0, and Epinano v1.2.0) on the mESC dataset were provided through personal communications with Prof. Luo Guanzheng, the corresponding author of the benchmark study referenced (Zhong et al., 2023).

m6A single molecule level benchmark

Request a detailed protocol

The benchmark IVT data for single molecule m6A identification was downloaded from NCBI-SRA with accession number SRP166020. CHEUI was used for the benchmark. For the ground truth, every A in every read of the IVT_m6A sample is treated as a m6A modification and every A in every read of the IVT_normalA is treated as a normal A. Detailed methods for calculating modification probability at single molecule level were provided in Appendix 1, ‘Modification state estimation on single molecule level’ and ‘Modification probability estimation on single molecule level’.

Results

RNA translocation hypothesis

Accurate segmentation of raw current signals in direct RNA sequencing remains a major challenge, largely because the precise dynamics of RNA translocation through the pore are not fully understood. To address this, we propose a hypothesis that better reflects the physical movement of the RNA molecule. In traditional basecalling algorithms such as Guppy and Albacore, we implicitly assume that the RNA molecule is translocated through the pore by the motor protein in a monotonic fashion, that is the RNA is pulled through the pore unidirectionally. In the DNN training process of Guppy and Albacore, we try to align the current signal with the reference RNA sequence. The alignment is unidirectional, which is the source of the implicit monotonic translocating assumption.

Contrary to the conventional assumption of monotonic translocation, the raw current data suggests that the motor protein drives RNA both forward and backward during sequencing. Figure 2B illustrates this with several example fragments of DRS raw current signal (Pratanwanich et al., 2021), where each fragment roughly corresponds to three neighboring 5-mers. The raw current signals, as shown in Figure 2B, strongly support this hypothesis, with several instances of measured current intensities matching both the previous and next 5-mer’s baseline. These repeated patterns, observed across multiple DRS datasets, provide empirical evidence of the jiggling translocation. Similar patterns are widely observed across the whole data. This suggests that the RNA molecule may move forward and backward while passing through the pore. This observation is also supported by previous reports (Craig et al., 2017; Caldwell and Spies, 2017), in which the helicase (the motor protein) translocates the DNA strand through the nanopore in a back-and-forth manner. Depending on ATP or ADP binding, the motor protein may translocate the DNA/RNA forward or backward by 0.5~1 nucleotides.

Figure 2 with 1 supplement see all

Download asset Open asset

RNA translocation hypothesis.

(A) Jiggling RNA translocation hypothesis. The top panel shows the raw current signal of Nanopore direct RNA sequencing, with gray areas representing SegPore-estimated transition blocks. We focus on three neighboring 5-mers, considering the central 5-mer (CTACG) as the current 5-mer. The RNA molecule may briefly move forward or backward during the translocation of the current 5-mer. If the RNA molecule is pulled backward, the previous 5-mer is placed in the pore, and the current signal (“prev” state, red dots) resembles the previous 5-mer’s baseline (mean and standard deviation highlighted by red lines and shades). If the RNA is pushed forward, the current signal (“next” state, blue dots) is similar to the next 5-mer’s baseline. (B) Example raw current signals supporting the jiggling hypothesis. The dashed rectangles highlight base blocks, with red and blue points representing measurements corresponding to the previous and next 5-mer, respectively. Red points align closely with the previous 5-mer’s baseline, and blue points match the next 5-mer’s baseline, reinforcing the hypothesis that the RNA molecule jiggles between neighboring 5-mers. The raw current signals were extracted from mESC WT samples of the training data in the m6A benchmark experiment.

Based on the reported kinetic model (Craig et al., 2017), we hypothesize that the RNA is translocated through the pore in a jiggling manner. This “jiggling” hypothesis presents a significant departure from the traditionally accepted view of unidirectional RNA translocation and aligns better with recent evidence from studies on DNA translocation (Craig et al., 2017; Caldwell and Spies, 2017). Incorporating this hypothesis into segmentation algorithms could lead to more accurate predictions of RNA modifications by accounting for the dynamic nature of translocation. On average, the motor protein sequentially translocates 5-mers on the RNA strand forward, and each 5-mer resides in the pore for a short time. During the short period of a single 5-mer, the motor protein may swiftly drive the RNA molecule forward and backward by 0.5~1 nucleotide in the translocation process of the current 5-mer (Figure 2A), which makes the measured current intensity occasionally similar to the previous or the next 5-mer. When the motor protein does not move the RNA molecule, we hypothesize that the RNA molecule undergoes slight thermal fluctuations, causing it to oscillate slightly within the pore and produce a stable current close to its baseline. In contrast, sharp changes in current intensity between consecutive 5-mers define transition blocks, where one 5-mer is replaced by the next.

We also assume that the raw current signal of a read can be segmented into a series of alternating base and transition blocks. In the ideal case, a base block corresponds to the base state where the 5-mer resides in the pore and jiggles between neighboring 5-mers, that is the current 5-mer can transiently jump to the previous or the next 5-mer. A transition block corresponds to the transition state between two consecutive base states where one 5-mer translocates to the next 5-mer in the pore. The current signal should be relatively flat in the base blocks, while a sharp change is expected in the transition blocks. One challenge we encountered was the overestimation of transition blocks. This may be due to a base block actually corresponding to a sub-state of a single 5-mer, rather than each base block corresponding to a full 5-mer, leading to inflated transition counts. To address this issue, SegPore’s alignment algorithm was refined to merge multiple base blocks (which may represent sub-states of the same 5-mer) into a single 5-mer, thereby facilitating further analysis.

SegPore Workflow

The SegPore workflow (Figure 1) consists of five key steps: (1) Preprocess fast5 files to pair raw current signal segments with corresponding RNA sequence fragments; (2) Segment each raw current signal using a hierarchical hidden Markov model (HHMM) into base and transition blocks; (3) Align the derived base blocks with the paired RNA sequence; (4) Fit a two-component GMM to estimate the modification state at the single-molecule level or the modification rate at the site level; (5) Use the results from Step (4) to update relevant parameters. Steps (3) to (5) are iterative and continue until the estimated parameters stabilize.

The final outputs of SegPore are the events and modification state predictions. SegPore’s events are similar to the outputs of Nanopolish’s ‘eventalign’ command, in that they pair raw current signal segments with the corresponding RNA reference 5-mers. Each 5-mer is associated with various features, such as start and end positions, mean current, and standard deviation, derived from the paired signal segment. For 5-mers that exhibit one clearly unmodified component and one clearly modified component, SegPore reports the modification rate at each site, as well as the modification state of that site on individual reads.

A key component of SegPore is the 5-mer parameter table, which specifies the mean and standard deviation for each 5-mer in both modified and unmodified states (Figure 1A). Since the peaks (representing modified and unmodified states) are separable for only a subset of 5-mers, SegPore can provide modification parameters for these specific 5-mers. For other 5-mers, modification state predictions are unavailable.

Segmentation benchmark

To evaluate SegPore’s performance in raw signal segmentation, we compared SegPore with Nanopolish (v0.14.0) and Tombo (v1.5.1) using three Nanopore direct RNA sequencing (DRS) datasets: two HEK293T datasets (wild type and Mettl3 knockout) (Pratanwanich et al., 2021) and the HCT116 dataset (Chen et al., 2025). Nanopolish and SegPore employed the ‘eventalign’ method to align 5-mers on each read with their corresponding raw signals, producing the mean and standard deviation (std) of the aligned signal segments. Tombo used the ‘resquiggle’ method to segment the raw signals, but the resulting signals are not reported on the absolute pA scale. To ensure a fair comparison with SegPore, we standardized the segments using the poly(A) tail in the same way as SegPore (See preprocessing section in Materials and methods).

To benchmark segmentation performance, we used two key metrics (details provided in Appendix 1, Section 8): (1) the log-likelihood of the segment mean, which measures how closely the segment matches ONT’s 5-mer parameter table (used as ground truth), and (2) the standard deviation (std) of the segment, where a lower std indicates reduced noise and better segmentation quality. If the raw signal segment aligns correctly with the corresponding 5-mer, its mean should closely match ONT’s reference, yielding a high log-likelihood. A lower std of the segment reflects less noise and better performance overall.

As shown in Table 1, SegPore consistently achieved the best performance averaged on all 5-mers across all datasets, with the highest log-likelihood and the lowest std values. These results suggest that SegPore provides a more accurate segmentation of the raw signal with significantly reduced noise compared to Nanopolish and Tombo. It is worth noting that the data points corresponding to the transition state between two consecutive 5-mers are not included in the calculation of the standard deviation in SegPore’s results in Table 1. However, their exclusion does not affect the overall conclusion. Because SegPore contains on average ~6 points per event within the transition state, we similarly removed the first and last three data points of each event for Nanopolish and Tombo prior to recalculating the metrics. The updated results are presented in Table 2.

Table 1

Segmentation benchmark on RNA002 data.

TestDataset	Avg. std ( $\hat{σ}$ ) ↓			Avg. log p ( $\hat{L}$ ) ↑
TestDataset	Nanopolish	Tombo	SegPore	Nanopolish	Tombo	SegPore
HEK293T(WT)	3.073	4.187	2.736	–2.871	–3.749	–2.778
HEK293T(KO)	2.948	4.204	2.670	–2.856	–3.749	–2.745
HCT116	3.167	4.076	2.872	–2.872	–3.704	–2.746

Table 2

Segmentation benchmark on the RNA002 data excluding boundaries.

TestDataset	Avg. std ( $\hat{σ}$ ) ↓			Avg. log p ( $\hat{L}$ ) ↑
TestDataset	Nanopolish	Tombo	SegPore	Nanopolish	Tombo	SegPore
HEK293T(WT)	2.862	4.090	2.736	–2.838	–3.772	–2.778
HEK293T(KO)	2.856	4.120	2.670	–2.878	–3.729	–2.745
HCT116	3.058	4.053	2.872	–2.869	–3.708	–2.746

To provide a more intuitive comparison, the segmentation results for two example raw signal clips are illustrated in Figure 3—figure supplement 2. These examples demonstrate the clearer, more precise segmentation achieved by SegPore compared to Nanopolish.

To evaluate SegPore’s performance on RNA004 data, we compared it with f5c (v1.5) (Gamaarachchi et al., 2020) and Uncalled4 (Kovaka et al., 2025) (v4.1.0) using three public DRS datasets: the S. cerevisiae dataset (Watson et al., 2025), the curlcake IVT and m6A datasets (Cruciani et al., 2025). The RNA002 data provides reference current levels for 5-mers, whereas RNA004 provides reference values for 9-mers, with Uncalled4 normalizing them to approximately zero mean and unit variance. As there are currently no established poly(A) detection methods available for RNA004, we used f5c to standardize the raw signals of each read prior to segmentation. SegPore was then applied to perform segmentation on the standardized signals. We computed the same benchmarking metrics—average log-likelihood and standard deviation—using the standardized raw signals and the segmentation results from f5c, Uncalled4, and SegPore. The 9-mer parameter table in pA scale for RNA004 data provided by f5c (see Materials and methods) was used in the analysis. As shown in Table 3, SegPore achieved the best overall performance across all datasets, indicating its robustness and suitability for RNA004 data. Moreover, we find that the jiggling hypothesis remains valid for RNA004, as illustrated in Figure 2—figure supplement 1.

Table 3

Segmentation benchmark on RNA004 data.

TestDataset	Avg. std ( $\hat{σ}$ ) ↓			Avg. log p ( $\hat{L}$ ) ↑
TestDataset	f5c	Uncalled4	SegPore	f5c	Uncalled4	SegPore
S. cerevisiae	3.129	3.970	3.125	–2.541	–3.835	–2.489
curlcake(IVT)	3.424	4.729	3.359	–2.716	–2.904	–2.515
curlcake(m6A)	3.520	4.971	3.392	–3.118	–3.524	–2.599

m6A identification at the site level

We evaluated SegPore’s performance in raw signal segmentation and m6A identification using independent public datasets as both training and test data. Since m6A modifications typically occur at DRACH motifs (where D denotes A, G, or U, and H denotes A, C, or U) (Linder et al., 2015), this study focuses on estimating m6A modifications on these motifs.

To begin, we estimated the 5-mer parameter table for m6A modifications using Nanopore direct RNA sequencing (DRS) data from three wild-type human HEK293T cell samples (Pratanwanich et al., 2021). The fast5 files from all samples were concatenated, and the full SegPore workflow was run to obtain the 5-mer parameter table. The parameter estimation process was iterated five times to ensure stabilization, refining the modification parameters for ten 5-mers where the modification state distribution significantly differs from the unmodified state.

Next, we applied SegPore’s segmentation and m6A identification to test data from wild-type mouse embryonic stem cells (mESCs; Zhong et al., 2023). Given the comparable methods and input data requirements, we benchmarked SegPore against several baseline tools, including Tombo, MINES (Lorenz et al., 2020), Nanom6A (Gao et al., 2021), m6Anet, Epinano (Liu et al., 2019), and CHEUI (Acera Mateos et al., 2024). By default, MINES and Nanom6A use eventalign results generated by Tombo, while m6Anet, Epinano, and CHEUI rely on eventalign results produced by Nanopolish. In Figure 3C, ‘Nanopolish +m6Anet” refers to the default m6Anet pipeline, whereas ‘SegPore +m6Anet’ denotes a configuration in which Nanopolish’s eventalign results are replaced with those from SegPore. Based on the output of SegPore eventalign, we fine-tune m6Anet using the HCT116 data at all DRACH motifs, aiming to demonstrate that the performance of SegPore benefits downstream models. Additionally, due to the differences in the availability of ground truth labels for specific k-mer motifs between human and mouse (Figure 3—figure supplement 1), six shared 5-mers were selected to demonstrate SegPore’s performance in modification prediction directly. By utilizing the 5-mer parameter table derived from the training data, SegPore employs a two-component GMM to calculate the modification rates at the selected m6A sites.

Figure 3 with 3 supplements see all

Download asset Open asset

m6A identification at the site level.

(A) Histogram of the estimated mean from current signals mapped to an example m6A-modified genomic location (chr10:128548315, GGACT) across all reads in the training data, comparing Nanopolish (left) and SegPore (right). The x-axis represents current in picoamperes (pA). (B) Histogram of the estimated mean from current signals mapped to the GGACT motif at all annotated m6A-modified genomic locations in the training data, again comparing Nanopolish (left) and SegPore (right). The x-axis represents current in picoamperes (pA). (C) Site-level benchmark results for m6A identification across all DRACH motifs, showing performance comparisons between SegPore+m6Anet and other methods. (D) Benchmark results for m6A identification on six selected motifs at the site level, comparing SegPore and other baseline methods.

Figure 3—source data 1 site_mod_rate_mES_WT_all_motifs.: https://cdn.elifesciences.org/articles/104618/elife-104618-fig3-data1-v1.csv
Download elife-104618-fig3-data1-v1.csv
Figure 3—source data 2 site_mod_rate_mES_WT_selected_motifs.: https://cdn.elifesciences.org/articles/104618/elife-104618-fig3-data2-v1.csv
Download elife-104618-fig3-data2-v1.csv

SegPore demonstrated improved segmentation compared to Nanopolish. Figure 3A shows the eventalign results at an example genomic location with m6A modifications. SegPore’s results show a more pronounced bimodal distribution in the raw signal segment mean, indicating clearer separation of modified and unmodified signals. Furthermore, when pooling all reads mapped to m6A-modified locations at the GGACT motif, SegPore exhibited more distinct peaks (Figure 3B), indicating reduced noise and potentially enabling more reliable modification detection.

We evaluated m6A predictions using two approaches: (1) SegPore’s segmentation results were fed into m6Anet, referred to as SegPore+m6Anet, which works for all DRACH motifs and (2) direct m6A predictions from SegPore’s Gaussian Mixture Model (GMM), which is limited to the six selected 5-mers shown in Figure 3—figure supplement 1C that exhibit clearly separable modified and unmodified components in the GMM (see Materials and methods for details).

In terms of m6A identification, SegPore performed strongly on the test data. Using miCLIP2 (Körtel et al., 2021) data as the ground truth, we calculated the area under the curve (AUC) for both the receiver operating characteristic (ROC) and precision-recall (PR) curves. SegPore+m6Anet achieved the best performance with an ROC AUC of 89.0% and a PR AUC of 43.8% (Figure 3C). For six selected m6A motifs, SegPore achieved an ROC AUC of 82.7% and a PR AUC of 38.7%, earning the third best performance compared with deep learning methods m6Anet and CHEUI (Figure 3D). It is noteworthy that SegPore’s GMM for m6A estimation is a very simple model, utilizing far fewer parameters than DNN-based methods. Achieving a decent performance with such a simple model is a significant accomplishment. These results highlight SegPore’s robust performance in m6A identification. For practical applications, we recommend taking the intersection of m6A sites predicted by SegPore and m6Anet to obtain high-confidence modification sites, while still benefiting from the interpretability provided by SegPore’s predictions.

m6A identification at the single molecule level

SegPore naturally identifies m6A modifications at the single-molecule level, which is crucial for understanding the heterogeneity of RNA modifications across individual transcripts. We benchmarked SegPore against CHEUI using an in vitro transcription (IVT) dataset containing two samples—one transcribed with m6A and the other with adenine (Jenjaroenpun et al., 2021). This dataset provides clear ground truth for m6A modifications at the single-molecule level, with all adenosine positions replaced by m6A in the ivt_m6A sample, and all adenosines unmodified in the ivt_normalA sample. We used 60% of the data for training and 40% for testing, with both SegPore and CHEUI estimating the m6A modification probability at each adenosine site on each read. Based on these probabilities, we calculated the ROC-AUC and PR-AUC by varying the modification probability threshold.

As shown in Figure 4A, SegPore outperformed CHEUI on this benchmark dataset, achieving better performance in both PR-AUC (94.7% vs 90.3%) and ROC-AUC (93% vs 89.9%). These results clearly demonstrate SegPore’s accuracy and robustness in detecting single-molecule m6A modifications.

Figure 4 with 2 supplements see all

Download asset Open asset

m6A identification at the single-molecule level.

(A) Benchmark results for single-molecule m6A identification on IVT data. SegPore shows better performance compared to CHEUI in both PR-AUC and ROC-AUC. (B) Comparison of ‘eventalign’ results from SegPore and Nanopolish for five consecutive k-mers. Note that DRS is sequenced from 3’ to 5’, so the k-mers enter the pore from right to left. A total of 100 reads were randomly sampled from transcript locations A1 (positions 711–719) in both the IVT_normalA and IVT_m6A samples (SRA: SRP166020). Each line represents an individual read, and the y-axis shows the raw signal intensity in picoampere (pA). Pink lines represent the IVT_m6A sample, and gray lines represent the IVT_normalA sample. The k-mers ‘GCGGA’, ‘CGGAC’, ‘GGACT’, ‘GACTT’, and ‘ACTTT’ all contain N6-Methyladenosine (m6A) in the IVT_m6A sample. SegPore’s results show clearer separation between m6A and adenosine, especially for ‘CGGAC’ and ‘GGACT’, compared to Nanopolish. (C) The upper panel shows the modification rate for selected genomic locations in the example gene ENSMUSG00000003153. The lower panel displays the modification states of all reads mapped to this gene. The black borders in the heatmap highlight the biclustering results, showing distinct modification patterns across different read clusters labeled C1 through C6.

Figure 4—source data 1 ENSMUSG00000003153.10.info.: https://cdn.elifesciences.org/articles/104618/elife-104618-fig4-data1-v1.csv
Download elife-104618-fig4-data1-v1.csv
Figure 4—source data 2 ENSMUSG00000003153.10.reads.: https://cdn.elifesciences.org/articles/104618/elife-104618-fig4-data2-v1.csv
Download elife-104618-fig4-data2-v1.csv

Next, we demonstrate SegPore’s interpretability through an example comparison of raw signal clips from both ivt_m6A and ivt_normalA samples (Figure 4B). The raw signal segments (means) are aligned to a randomly selected m6A site as well as its neighboring sites in the reference sequence. Each line represents an individual read, with pink lines from the ivt_m6A sample and gray lines from the ivt_normalA sample. SegPore’s segmentation clearly distinguishes between m6A and adenosine at the single-molecule level, while Nanopolish’s segmentation shows less distinction.

Interestingly, we observed that the position of m6A within a 5-mer can affect the signal intensity. For instance, a clear difference between m6A and adenosine is evident when m6A occupies the fourth position in the 5-mer (e.g. ‘CGGAC’), but this difference is less pronounced when m6A is in the second position (e.g. ‘GACTT’).

To illustrate the benefits of single-molecule m6A estimation, we present an example gene (ENSMUSG00000003153) from the mESC dataset used in site-level m6A identification. Figure 4C shows a heatmap of highly modified genomic locations (modification rate >0.1), where rows represent reads and columns represent genomic locations. Biclustering reveals six clusters of reads and three clusters of genomic locations.

The heatmap in Figure 4C suggests heterogeneity in m6A modification patterns across different reads of the same gene. Biclustering reveals that modifications at g6 are specific to cluster C4, g7 to cluster C5, and g8 to cluster C6, while the first five genomic locations (g1 to g5) show similar modification patterns across all reads. Additionally, high modification rates are observed at the 3rd and 5th positions across the majority of reads. These results suggest that m6A modification patterns can vary significantly even within a single gene. This observation highlights the complexity of RNA modification regulation and underscores the importance of single-molecule resolution for understanding RNA function at a finer scale.

Discussion

One of the main computational challenges in direct RNA sequencing (DRS) lies in accurately segmenting the raw current signal. We developed a segmentation algorithm that models the jiggling property in the physical process of DRS, resulting in cleaner current signals for m6A identification at both the site and single-molecule levels. Our results demonstrate that SegPore’s segmentation enables clear differentiation between m6A-modified and unmodified adenosines. We believe that the de-noised current signals will be beneficial for other downstream tasks, such as the estimation of m5C, pseudouridine, and other RNA modifications. Nevertheless, several open questions remain for future research. In SegPore, we assume a drastic change between two consecutive 5-mers, which may hold for 5-mers with large difference in their current baselines but may not hold for those with small difference. As with other RNA modification estimation methods, SegPore can be affected by misalignment errors, particularly when the baseline signals of adjacent k-mers are similar. These cases may lead to spurious bimodal signal distributions and require careful interpretation. Another key question concerns the physical interpretation of the derived base blocks. Ideally, one base block would correspond to a single 5-mer, but in practice, multiple base blocks often align with one 5-mer. We hypothesize that the HHMM may segment a 5-mer’s current signal into multiple base blocks, where the 5-mer oscillates between different sub-states, each characterized by distinct baselines.

Currently, SegPore models only the modification state of the central nucleotide within the 5-mer. However, modifications at other positions may also affect the signal, as shown in Figure 4B. Therefore, introducing multiple states for a 5-mer to account for modifications at multiple positions within the same 5-mer could help to improve the performance of the model. This approach, however, would significantly increase model complexity—introducing two states per location results in $2^{5} = 32$ modification states per 5-mer, a challenge for simple GMMs. It remains to be investigated how to model multiple modification states properly to improve the performance in RNA modification estimation. While the current two-state assumption simplifies the model, it has proven effective, as demonstrated by improved performance on m6A benchmarks at both site and single-molecule levels.

A key advantage of SegPore is its interpretability, which sets it apart from DNN-based methods. SegPore provides a clearer understanding of RNA modification predictions. A limitation of DNN-based approaches is that they often struggle to differentiate between m6A and adenosine signals, as their intensity levels are quite similar (Figure 4B, Figure 3—figure supplement 3). This poses a challenge when applying DNN-based methods to new datasets without short read sequencing-based ground truth. In such cases, it is difficult for researchers to confidently determine whether a predicted m6A modification is genuine. SegPore encodes current intensity levels for different 5-mers in a parameter table, where unmodified and modified 5-mers are modeled using two Gaussian distributions. One can generally observe a clear difference in the intensity levels between 5-mers with an m6A and those with adenosine, which makes it easier for a researcher to interpret whether a predicted m6A site is genuine (see Figure 4—figure supplement 1). A challenge is the relatively small number of 5-mers that show significant changes in their modification states. To improve accuracy, larger training datasets and expanding the scope to 7-mers or 9-mers could help capture more context, potentially revealing more significant baseline changes.

The limited improvement of SegPore combined with m6Anet over Nanopolish+m6Anet in bulk in vivo analysis (Figure 3) may be explained by several factors: potential alignment inaccuracies due to pseudogenes or transcript isoforms, the complexity of in vivo datasets containing additional RNA modifications (e.g. m5C, m7G) affecting signal baselines, and the fact that m6Anet is specifically trained on events produced by Nanopolish rather than SegPore. Additionally, the lack of a modification-free control (in vitro transcribed) sample in the benchmark dataset makes it difficult to establish true baselines for each k-mer. Despite these limitations, SegPore demonstrates clear improvement in single-molecule m6A identification in IVT data (Figure 4), suggesting it is particularly well suited for in vitro transcription data analysis.

Although SegPore provides clear interpretability, there is potential to explore DNN-based models that can directly leverage SegPore’s segmentation results. Currently, m6Anet computes features (e.g., mean and standard deviation) from raw signal segments, which are then fed into a neural network for m6A prediction. However, a more direct approach—where raw signal segments are used as input to a DNN—could allow for the extraction of more intricate features that may exist within the signal but are currently underexplored. These high-order features might capture subtle aspects of the raw signal, leading to improved m6A estimation. Since SegPore currently models only the mean and standard deviation, further work involving advanced DNNs could extend beyond this, uncovering finer patterns in the signal that traditional statistical models might miss.

Computation speed is also a concern when processing fast5 files. We addressed this by implementing a GPU-accelerated inference algorithm in SegPore, resulting in a 10- to 20-fold speedup compared to the CPU-based version. We believe that the GPU implementation will unlock the full potential of SegPore for a wider range of downstream tasks and larger datasets. More details are provided in Appendix 2. SegPore’s running times on datasets of varying sizes, using a single NVIDIA DGX-1 V100 GPU and one CPU core, are provided in Figure 4—figure supplement 2.

In summary, we developed a novel software SegPore that considered the conformation changes of motor protein to segment raw current signal of Nanopore direct RNA sequencing. SegPore effectively masks out noise in the raw signal, leading to improved m6A identification at both site and single-molecule levels.

Appendix 1

Notation

$\| \cdot \|$	Take the length of the inside vector.
$Y$	All reads in the input multi-fast5 file. $Y = (y_{1}, y_{2}, \dots, y_{i}, \dots)$ . $y_{i}$ are the raw current measurements of ith read. Assuming there are $N$ reads, i.e. $\| Y \| = N$ .
$y_{i}$	The current signal of ith read. $y_{i} = (y_{i 1}, y_{i 2}, \dots y_{i j}, \dots y_{i l_{i}})$ where $l_{i}$ is the length of $y_{i}$ i.e. $\| y_{i} \| = l_{i}$ .
$y_{i}^{(k)}$	The kth mapped current signal fragment of ith read $y_{i}$ , given by eventalign using Nanopolish.
$R$	Reference sequences for each read. $R = (r_{1}, r_{2}, \dots r_{i}, \dots)$ .
$r_{i}$	The reference sequence corresponding to ith read $y_{i}$ , given by basecalling and minimap2. $r_{i} = (r_{i 1}, r_{i 2}, \dots r_{i j}, \dots, r_{i m_{i}})$ , where $r_{i j} \in {A, C, G, T}$ and $m_{i}$ is the length of $r_{i}$ i.e. $\| r_{i} \| = m_{i}$ .
$r_{i}^{(k)}$	The kth reference sequence fragment of $r_{i}$ corresponding to ${y_{i}}^{(k)}$ . Note that $\| {y_{i}}^{(k)} \|$ is the larger than $\| {r_{i}}^{(k)} \|$ in general.
$s_{i}$	The 5mer list corresponds to $r_{i}$ . We have $s_{i} = (s_{i 1}, s_{i 2}, \dots, s_{i j}, \dots, s_{i m_{i}})$ , where $s_{i j} = r_{i, j - 2} r_{i, j - 1} r_{i j} r_{i, j + 1} r_{i, j + 2}$ and $s_{i j} \in {A A A A A, A A A A C, A A A A G, A A A A T, \dots, T T T T T}$ . Note that $\| s_{i} \| = \| r_{i} \| = m_{i}$ .
${s_{i}}^{(k)}$	The k-mer list fragment corresponds to the kth reference sequence fragment ${r_{i}}^{(k)}$ . $\| s_{i}^{(k)} \| = \| r_{i}^{(k)} \|$ .
$Z$	The estimated states of all reads given by the Hierarchical Hidden Markov Model (HHMM). $Z = (z_{1}, z_{2}, \dots, z_{i}, \dots)$ .
$z_{i}$	The hidden states for each current measure of ith read. $z_{i} = (z_{i 1}, z_{i 2}, \dots, z_{i j}, \dots, z_{i l_{i}})$ , where $z_{i j} \in {0, 1, 2, 3, 4}$ corresponds to $y_{i j}$ and $\| y_{i} \| = \| z_{i} \| = l_{i}$ , where $z_{i}^{(k)} = 0$ means the transition state and $z_{i j} = 1, 2, 3, 4$ means base states.
${z_{i}}^{(k)}$	The kth fragment in $z_{i}$ , which corresponds to ${y_{i}}^{(k)}$ . Note that $\| z_{i}^{(k)} \| = \| y_{i}^{(k)} \|$ .
$μ_{i}$	The mean values of base segments given by HHMM, which are estimated from $y_{i}$ and $z_{i}$ . $μ_{i} = (μ_{i 1}, μ_{i 2}, \dots, μ_{i j}, \dots, μ_{i n_{i}})$ , where $j = 1 \dots n_{i}$ indexes the base segments of ith read. In general, we have $l_{i} > n_{i} > m_{i}$ .
${μ_{i}}^{(k)}$	The kth fragment of $μ_{i}$ which corresponds to ${y_{i}}^{(k)}$ and ${z_{i}}^{(k)}$ . $μ_{i}^{(k)} = (μ_{i 1}^{(k)}, μ_{i 2}^{(k)}, \dots, μ_{i j}^{(k)}, \dots)$ . In general we have $\| y_{i}^{(k)} \| > \| μ_{i}^{(k)} \| > \| s_{i}^{(k)} \|$ .
${u_{i}}^{(k)}$	The modification state of each 5mer in ${s_{i}}^{(k)}$ . $u_{i}^{(k)} = (u_{i 1}^{(k)}, u_{i 2}^{(k)}, \dots, u_{i j}^{(k)}, \dots)$ , where $u_{i j}^{(k)} \in {u n, m o d}$ indicates if ${s_{i j}}^{(k)}$ is in unmodified or modified state. $\| u_{i}^{(k)} \| = \| s_{i}^{(k)} \|$ .
$T^{r e f}$	The reference 5mer parameter table provided by Oxford Nanopore Technology (ONT). $T^{r e f} = {(μ_{s, u}, σ_{s, u}), \forall s \in {A A A A A, A A A A C, \dots, T T T T T}, u = ‘ ‘ u n "}$ .
$T^{k m e r}$	The parameter table for all 5mers. $T^{k m e r} = {(μ_{s, u}, σ_{s, u}), \forall s \in {A A A A A, A A A A C, \dots, T T T T T}, \forall u \in {‘ ‘ u n ", ‘ ‘ m o d "}}$ .

Introduction

Given the parameter table $T^{k m e r}$ for all $k$ -mers, the overall workflow of SegPore consists of the following major steps: (1) Base-calling and mapping, (2) preprocessing of nanopore raw signal, (3) segmenting nanopore raw signal using hierarchical hidden Markov model (HHMM), and (4) aligning raw signal segments with corresponding $k$ -mer list. The input of Segpore is the raw current signal of Nanopore direct RNA sequencing. The output of the Segpore workflow is the pairing of raw signal segments with corresponding $k$ -mers with inferred state (unmodified state or modified state). Unless otherwise noted, the following analysis focuses on RNA002 data, using 5-mers rather than $k$ -mers.

As illustrated in Appendix 1—figure 1, we denote the input FAST5/POD5 file as a set of raw current signal reads $Y = (y_{1}, y_{2}, \dots, y_{i}, \dots)$ . We use $| Y |$ to denote the number of reads and $| y_{i} |$ to denote the number of measurements / data points of $y_{i}$ . We then perform basecalling using Guppy and map basecalled sequence to the reference sequences (‘Basecalling, mapping and preprocessing’), yielding a pairing between the ith read $y_{i}$ and its corresponding reference sequence $r_{i}$ . We assume there are $k$ matched fragments $[(y_{i}^{(1)}, r_{i}^{(1)}), (y_{i}^{(2)}, r_{i}^{(2)}), \dots, (y_{i}^{(k)}, r_{i}^{(k)})]$ . For the jth fragment $y_{i}^{(j)}$ , we use HHMM (‘Hierarchical hidden Markov model’) to estimate the hidden states for each measurement $z_{i}^{(j)} = (z_{i 1}^{(j)}, z_{i 2}^{(j)}, \dots)$ , where $z_{i m}^{(j)} \in {0, 1, 2, 3, 4}$ and $| y_{i}^{(j)} | = | z_{i}^{(j)} |$ . Based on $y_{i}^{(j)}$ and $z_{i}^{(j)}$ , we can convert it into a list of low variation segments, with their means denoted by $μ_{i}^{(j)} = (μ_{i 1}^{(j)}, μ_{i 2}^{(j)}, \dots)$ . As a result, the number of measurements in $y_{i}^{(j)}$ is larger than the number of segments in $μ_{i}^{(j)}$ , that is $| y_{i}^{(j)} | = | z_{i}^{(j)} | > | μ_{i}^{(j)} |$ . After that, we perform full or partial alignment between the raw signal segments $μ_{i}^{(j)}$ and reference sequence $r_{i}^{(j)}$ , which is converted into a list of 5-mers $s_{i}^{(j)} = (s_{i 1}^{(j)}, s_{i 2}^{(j)}, \dots)$ and $∣ \cdot ∣$ . Based on $μ_{i}^{(j)}$ and $s_{i}^{(j)}$ , we perform full or partial alignment (‘Alignment of signal segments with reference sequence’) to match each 5-mer with one or multiple current signal segments, as well as estimating its modification state $u_{i}^{(j)} = (u_{i 1}^{(j)}, u_{i 2}^{(j)}, \dots)$ , where $u_{i m}^{(j)} \in {‘ ‘ u n ", ‘ ‘ m o d "}$ .

Appendix 1—figure 1

Download asset Open asset

SegPore workflow with a fixed 5-mer parameter table $T^{k m e r}$ .

The input is the raw current signal $y_{i}$ and reference sequence $r_{i}$ for ith read. The output is the paired current signal segments µ_i, 5-mer list $s_{i}$ and modification states $u_{i}$ .

However, $T^{k m e r}$ is not fully known beforehand (Appendix 1—table 1), especially the parameters for the modification states of 5-mers. We have to estimate the parameters from the training data through an iterative process (‘Estimation of k-mer parameters table’). After obtaining the estimated $T^{k m e r}$ , we derive the modification states from new test data using SegPore workflow with fixed $T^{k m e r}$ .

Appendix 1—table 1

Example of T^kmer (1024×4).

‘-’ represents ‘NA’.

modle_kmer	un_mu	un_sigma	mod_mu	mod_sigma
AAAAA	108.9	2.68	-	-
AAAAC	107.75	2.68	-	-
AAAAG	101.72	2.68	-	-
AAAAT	112.77	2.68	-	-
AAACA	99.38	2.41	95.28	2.44
AAACC	100	1.32	96.82	2.42
AAACG	101.01	3.45	-	-
AAACT	106.91	2.18	101.98	2.29
AAAGA	110.54	4.06	-	-
AAAGC	107.69	4.06	-	-
AAAGG	108.29	4.06	-	-
AAAGT	108.73	4.06	-	-
AAATA	114.11	3.11	-	-
...	...	...	...	...
TTTTT	80.78	1.97	-	-

Basecalling, mapping, and preprocessing

This section describes the details about the basecalling and mapping process, as well as raw data preprocessing for standardizing the raw signal of all reads.

First, we obtain the matching between raw current signals and reference sequences. The raw multi-fast5 reads (input file) were basecalled using Guppy (v6.0.1), which results in a FASTQ file. Then we align the yielded FASTQ file to the reference sequences (we use cDNA reference) using minimap2 (v2.24-r1122)(Li, 2021), which provides a bam file that matches the basecalled reads with the reference sequence. After that, we feed the raw fast5 file, basecalled FASTQ file, bam file, and reference sequence to Nanopolish (v0.13.2) to get (1) matching fragments between the raw current signal $y_{i}^{(k)}$ and the reference sequence $r_{i}^{(k)}$ (as well as $s_{i}^{(k)}$ ) and (2) the current signal fragment for poly(A) part of the mRNA molecule.

Then, we standardize the raw current signal based on the poly(A) tail. Due to technical reasons, there exist variations in the current signal of the same sequence fragment across different reads. It is necessary to perform a standardization. Since all reads have a poly(A) tail in Nanopore direct RNA sequencing, we could standardize the reads such that the mean and std of poly(A) part is the same for all reads. The standardization is given by

\begin{aligned} s t a n d a r d i z e d_s i g n a l = \frac{r a w_s i g n a l - μ_{p o l y A}}{σ_{p o l y A}} * σ_{s t a n d P o l y A} + μ_{s t a n d P o l y A} \end{aligned}

where $μ_{s t a n d P o l y A} = 108.9$ and $σ_{s t a n d P o l y A} = 1.67$ (RNA002). Note that $y_{i}$ refers to the standardized signal throughout this document.

Hierarchical hidden Markov model

This section describes the hierarchical hidden Markov model (HHMM) to segment the raw current signal of a read. The input of HHMM is $y_{i}$ and the output is $z_{i}$ . For notation simplicity, we use $y$ and $z$ to represent $y_{i}$ and $z_{i}$ in this section, respectively.

Model assumptions

Here we have several assumptions regarding the raw signal. We first assume there is a sharp change of the baseline levels of the current signal between two consecutive 5-mers. This means there is a large change of current signal within a short period of time. If we fit a line to the corresponding signals, we expect to see a large slope. The second assumption is that the 5-mer resides in the Nanopore, oscillates forward and backward during a 5-mer event, while staying in the pore most of the time. Based on these two assumptions, we can divide the raw current signal into two states (hidden states of the outer HMM): the base states for 5-mer event and the transition state for the transition between two consecutive 5-mer events. We assume that the current signal always starts with a base state, then switches to a transition state, after that switches back to a base state or transition state alternately, and in the end, the signal should end in a base state. For the base state, we use an inner HMM to model the oscillations, which has four hidden states: ‘curr’, ‘prev’, ‘next’, and ‘noise’. For the transition state, we use a linear regression model to model the abrupt change of the current signal.

Outer HMM

We define the following notations for the HHMM. As illustrated in Appendix 1—figure 2, we assume there are $N$ time points of the current signal, that is $y = (y_{1}, y_{2}, \dots, y_{i}, \dots, y_{N})$ . We use $g = (g_{1}, g_{2}, \dots, g_{i}, \dots, g_{N})$ , $g_{i} \in {‘ ‘ B ", ‘ ‘ T "}$ to indicate the hidden states of the outer HMM, where $B$ and $T$ stand for ‘base’ and ‘transition’, respectively. Similarly, we use $h = (h_{1}, h_{2}, \dots, h_{i}, \dots, h_{N})$ , $h_{i} \in {‘ ‘ c u r r ", ‘ ‘ p r e v ", ‘ ‘ n e x t ", ‘ ‘ n o i s e "}$ , respectively. We then integrate $g$ and $h$ into a single vector $z = (z_{1}, z_{2}, \dots, z_{i}, \dots, z_{N})$ , where $z_{i} = {\tilde{g}}_{i} \cdot {\tilde{h}}_{i}$ . Here ${\tilde{g}}_{i}$ $(‘ ‘ B "\to 1, ‘ ‘ T "\to 0)$ and ${\tilde{h}}_{i}$ (‘curr’→1, ‘prev’→2, ‘next’→3, ‘noise’→4) are the numerical mappings of $g_{i}$ and $h_{i}$ . Based on $g$ , we can derive that $y$ is divided into $2 K + 1$ blocks, which is defined by $c = (c_{1}, c_{2}, \dots, c_{i}, \dots, c_{2 K + 1})$ , $c_{i} \in {‘ ‘ B ", ‘ ‘ T "}$ . Note that $c_{i} = ‘ ‘ B " \forall i = 1, 3, \dots, 2 K + 1$ and $c_{i} = ‘ ‘ T " \forall i = 2, 4, \dots, 2 K$ . As a result, we could also represent $y = (y^{(1)}, y^{(2)}, \dots, y^{(i)}, \dots, y^{(2 K + 1)})$ , where $y^{(i)}$ represent the ith block. The base blocks are given by $(y^{(b_{1})}, y^{(b_{2})}, \dots, y^{(b_{i})}, \dots, y^{(b_{K + 1})})$ , where $b_{i} = 2 i - 1$ and $i \in {1, 2, \dots, K + 1}$ . The transition blocks are given by $(y^{(t_{1})}, y^{(t_{2})}, \dots, y^{(t_{i})}, \dots, y^{(t_{K})})$ , where $t_{i} = 2 i$ and $i \in {1, 2, \dots, K}$ . Similarly, we have $g = (g^{(1)}, g^{(2)}, \dots, g^{(i)}, \dots, g^{(2 K + 1)})$ and $h = (h^{(1)}, h^{(2)}, \dots, h^{(i)}, \dots, h^{(2 K + 1)})$ .

Appendix 1—figure 2

Download asset Open asset

Example of HHMM where $g_{i} \in {‘ ‘ B ", ‘ ‘ T "}$ , $h_{i} \in {‘ ‘ c u r r ", ‘ ‘ p r e v ", ‘ ‘ n e x t ", ‘ ‘ n o i s e "}$ and $c_{i} \in {‘ ‘ B ", ‘ ‘ T "}$ .

In this figure, ‘c’, ‘p’, ‘n’, and ‘o’ represent ‘curr’, ‘prev’, ‘next’, and ‘noise’ respectively in $h$ . $| y | = | g | = | h | a n d | c | = 11$ .

We use $T^{o u t e r} =$ ${T_{B B}^{o u t e r}, T_{B T}^{o u t e r}, T_{T B}^{o u t e r}, T_{T T}^{o u t e r}}$ to denote the transition probabilities between the base state and the transition state in the outer HMM. We use $T^{i n n e r} =$ ${T_{h_{i} h_{j}}^{i n n e r} \forall h_{i}, h_{j} \in {‘ ‘ c u r r ", ‘ ‘ p r e v ", ‘ ‘ n e x t ", ‘ ‘ n o i s e "}}$ to denote the transition probabilities between four hidden states of the inner HMM. For base blocks, we have inner HMM parameters $ϕ = (ϕ^{(b_{1})}, ϕ^{(b_{2})}, \dots, ϕ^{(b_{K + 1})})$ and $σ = (σ^{(b_{1})}, σ^{(b_{2})}, \dots, σ^{(b_{K + 1})})$ . For transition blocks, we have linear coefficients $β_{0} = (β_{0}^{(t_{1})}, β_{0}^{(t_{2})}, \dots, β_{0}^{(t_{K})})$ , $β_{1} = (β_{1}^{(t_{1})}, β_{1}^{(t_{2})}, \dots, β_{1}^{(t_{K})})$ and noise std $σ_{ϵ} = (σ_{ϵ}^{(t_{1})}, σ_{ϵ}^{(t_{2})}, \dots, σ_{ϵ}^{(t_{K})})$ .

The joint likelihood of the HHMM is given by

\begin{array}{ll} p (y, g) & = p (y ∣ g) p (g) \\ = p (y ∣ g) π_{g_{1}}^{o u t e r} \prod_{i = 2}^{N} T_{g_{i - 1} g_{i}}^{o u t e r} \\ = \prod_{i = 1}^{N} p (y_{i} ∣ g_{i}) π_{g_{1}}^{o u t e r} \prod_{i = 2}^{N} T_{g_{i - 1} g_{i}}^{o u t e r} \end{array},

where $π_{g_{1}}^{o u t e r} = 1$ is the initial probability of hidden state of the $y_{1}$ (always ‘B’ state). However, it is not possible to directly derive $p (y_{i}, h_{i} | g_{i})$ , since $y_{i}$ are not independent given $g_{i}$ . Here we use different models for each block $c_{k}$ to calculate the joint probability $p (y^{(k)} | g^{(k)})$ to approximate $\prod_{i = 1}^{N} p (y_{i} | g_{i})$ . Therefore, we have

\begin{array}{ll} \prod_{i = 1}^{N} p (y_{i} ∣ g_{i}) = \prod_{j = 0}^{2 K + 1} p (y^{(j)} ∣ c_{j}) \\ = \prod_{j = 0}^{K} p (y^{(2 j + 1)} ∣ c_{2 j + 1} = ‘ ‘ B ") \prod_{j = 1}^{K} p (y^{(2 j)} ∣ c_{2 j} = ‘ ‘ T ") \end{array}

The first part is the marginal likelihood for base blocks, while the second part is the marginal likelihood for transition blocks. The specific likelihood formulas are given in the following sections.

For the outer HMM, the transition probabilities $T^{o u t e r}$ are pre-specified and the emission probabilities are obtained from the inner HMM (‘Inner HMM for base state’) and linear regression model (‘Linear model for transition state’). The initial probability for the first hidden state $π_{g_{1}}^{o u t e r}$ is set to 1, since the first state should always be “B”.

Inner HMM for base state

This section provides the derivations for calculating the likelihood of the ith base block $b_{i}$ , where $i \in {1, 2, \dots, K + 1}$ . We assume the data $y^{(b_{i})} = (y_{1}^{(b_{i})}, y_{2}^{(b_{i})}, \dots, y_{N_{b_{i}}}^{(b_{i})})$ is generated by an HMM with four states: ‘curr’, ‘prev’, ‘next’, and ‘noise’. The descriptions of the states are given by Appendix 1—table 2, where parameters for the noise state are given and parameters for other states need to be estimated from the data.

Appendix 1—table 2

Hidden states of inner HMM for the $b_{i}$ th block.

State	Event	Emission distribution	Parameter estimation
‘curr’	current 5-mer resides in pore	$N (ϕ^{(b_{i})}, (σ^{(b_{i})})^{2})$	Yes
‘prev’	previous 5-mer resides in pore	$N (ϕ^{(b_{i - 1})}, (σ^{(b_{i - 1})})^{2})$	Yes
‘next’	next 5-mer resides in pore	$N (ϕ^{(b_{i + 1})}, (σ^{(b_{i + 1})})^{2})$	Yes
‘noise’	unknown noise	$U n i f (l b, u b)$	No

The likelihood for the kth block is given by

\begin{array}{ll} p (y^{(b_{i})} ∣ c_{b_{i}} = ``B``) = \sum_{h^{(b_{i})}} p (y^{(b_{i})}, h^{(b_{i})} ∣ c_{b_{i}} = ``B``) \\ = \sum_{h^{(b_{i})}} p (y^{(b_{i})} ∣ ϕ, σ, h^{(b_{i})}, c_{b_{i}} = ``B``) p (h^{(b_{i})} ∣ c_{b_{i}} = ``B``) \\ = \sum_{h^{(b_{i})}} [\prod_{j = 1}^{N_{b_{i}}} p (y_{j}^{(b_{i})} ∣ h_{j}^{(b_{i})}, ϕ, σ) π_{h_{1}^{(b_{i})}}^{inner} \prod_{j = 2}^{N_{b_{i}}} T_{h_{j - 1}^{(b_{i})} h_{j}^{(b_{i})}}^{inner}] \\ = \sum_{h^{(b_{i})}} [\prod_{j = 1}^{N_{b_{i}}} {N {(y_{j}^{(b_{i})} ∣ ϕ^{(b_{i})}, (σ^{(b_{i})})^{2})}^{I (h_{j}^{(b_{i})} = ``curr``)} N {(y_{j}^{(b_{i})} ∣ ϕ^{(b_{i} - 1)}, (σ^{(b_{i} - 1)})^{2})}^{I (h_{j}^{(b_{i})} = ``prev``)} \\ N {(y_{j}^{(b_{i})} ∣ ϕ^{(b_{i} + 1)}, (σ^{(b_{i} + 1)})^{2})}^{I (h_{j}^{(b_{i})} = ``next``)} Unif {(y_{j}^{(b_{i})} ∣ l b, u b)}^{I (h_{j}^{(b_{i})} = ``noise``)}} \\ π_{h_{1}^{(b_{i})}}^{inner} \prod_{j = 2}^{N_{b_{i}}} T_{h_{j - 1}^{(b_{i})} h_{j}^{(b_{i})}}^{inner}] \end{array}

where $π_{h_{1}^{(b_{i})}}^{i n n e r} = 0.25$ is the probability of the first hidden state and $I (\cdot)$ is the indicator function

For the inner HMM, the transition probabilities $T^{i n n e r}$ are pre-specified and the initial probability for the first hidden state $π_{h_{1}^{(b_{i})}}^{i n n e r} = 0.25$ is set to 0.25. Other parameters (see Appendix 1—table 2) need to be estimated.

By comparing the emission probabilities of the outer HMM (Equation 3) and the marginal likelihood of the inner HMM (Equation 4), we have

p (y^{(b_{i})} = c_{b_{i}} = ‘ ‘ B ") = \prod_{j = 1}^{| y^{(b_{i})} |} p (y_{j}^{(b_{i})} | g_{j}^{(b_{i})})

To approximate $p (y_{j}^{(b_{i})} | g_{j}^{(b_{i})})$ , we assume all values of $p (y_{j}^{(b_{i})} | g_{j}^{(b_{i})})$ are the same for any given $j$ , we have

p (y^{(b_{i})} | c_{b_{i}} = ‘ ‘ B ") = \prod_{j = 1}^{| y^{(b_{i})} |} p (y_{j}^{(b_{i})} | g_{j}^{(b_{i})}) = [\hat{p} (y_{j}^{(b_{i})} | g_{j}^{(b_{i})})]^{| y^{(b_{i})} |}

The approximated emission probability of each data point in the ith base block is given by

\hat{p} (y_{j}^{(b_{i})} ∣ g_{j}^{(b_{i})}) = \sqrt[∣ y^{(b_{i})} ∣]{\hat{p} (y^{(b_{i})} ∣ c_{b_{i}} = ‘ ‘ B ")},

where $\hat{p} (y^{(b_{i})} | c_{b_{i}} = ‘ ‘ B ")$ can be obtained by our inference algorithm.

Note that the hidden states of the inner HMM $h$ can be provided as a by-product of our inference algorithm. We will use it to derive the final composite hidden state $z$ of the HHMM (‘Parameter inference’), which is utilized by downstream alignment tasks (‘Post-processing of alignment path’).

Linear model for transition state

This section provides the derivations for calculating the likelihood of the ith transition block $t_{i}$ , where $i \in {1, 2, \dots, K}$ . We denote the input data by $y^{(t_{i})} = (y_{1}^{(t_{i})}, y_{2}^{(t_{i})}, \dots, y_{N_{t_{i}}}^{(t_{i})})$ and its corresponding index by $x^{(t_{i})} = (1, 2, \dots, N_{t_{i}})$ . We assume a linear model for data points in the transition state, given by

y^{(t_{i})} = β_{1}^{(t_{i})} x^{(t_{i})} + β_{0}^{(t_{i})} + ϵ,

where $β_{1}^{(t_{i})}$ is the slope, $β_{0}^{(t_{i})}$ is the intercept and $ϵ \sim N (0, (σ_{ϵ}^{(t_{i})})^{2})$ is the Gaussian noise.

The likelihood of the linear model is given by

\begin{array}{ll} p (y^{(t_{i})} ∣ c_{t_{i}} = ‘ ‘ T ") & = p (y^{(t_{i})} ∣ x^{(t_{i})}, β_{0}^{(t_{i})}, β_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})}) \\ = \prod_{j = 1}^{N_{t_{i}}} p (y_{j}^{(t_{i})} ∣ x_{j}^{(t_{i})}, β_{0}^{(t_{i})}, β_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})}) \\ = \prod_{j = 1}^{N_{t_{i}}} N (y_{j}^{(t_{i})} - β_{1}^{(t_{i})} x_{j}^{(t_{i})} - β_{0}^{(t_{i})} ∣ 0, (σ_{ϵ}^{(t_{i})})^{2}) \end{array}

Similar to the inner HMM, the approximated emission probability for each data point in ith transition block is given by

\hat{p} (y_{j}^{(t_{i})} | g_{j}^{(t_{i})}) = | y^{(t_{i})} | \sqrt{\hat{p} (y^{(t_{i})} | c_{t_{i}} = ‘ ‘ T ")},

\hat{p} (y^{(b_{i})} = c_{b_{i}} = ‘ ‘ T ")

Parameter Inference

Due to the complexity of GPU-accelerated inference algorithm, the details are provided in Appendix 2.

The inference algorithm will infer the following parameters: (1) the hidden states $g$ of the outer HMM; (2) the hidden states of the inner HMM $h$ ; (3) relevant parameters for emissions (Appendix 1—table 2) of each current signal fragment of the $2 k + 1$ th block (base block) $y^{(2 k + 1)}$ , where $k \in (0, 1, \dots, K)$ and $c_{2 k + 1} = ‘ ‘ B "$ ; and (4) the linear coefficients and noise std of the $2 k t h$ block (transition block) $y^{(2 k)}$ , where $k \in (1, 2, \dots, K)$ and $c_{2 k} = ‘ ‘ T "$ .

We then integrate $g$ and $h$ into a single vector $z = (z_{1}, z_{2}, \dots, z_{i}, \dots, z_{N})$ , where $z_{i} = {\tilde{g}}_{i} \cdot {\tilde{h}}_{i}$ and ${\tilde{g}}_{i}, {\tilde{h}}_{i}$ are the numerical mappings of $g_{i}, h_{i}$ . $z$ is the final output of HHMM.

Alignment of signal segments with reference sequence

This section describes the alignment algorithms for aligning an input fragment of signal segments $μ_{i}^{(k)}$ and the corresponding 5-mer list $s_{i}^{(k)}$ . After post-processing of the alignment path, we could get the modification states $u_{i}^{(k)}$ . For notation simplicity, we drop the indices and denote $μ_{i}^{(k)}$ , $s_{i}^{(k)}$ and $u_{i}^{(k)}$ by μ, $s$ and $u$ , respectively. In general, one 5-mer $s_{j}$ is aligned with $k$ signal segments $(μ_{i}, μ_{i + 1}, \dots, μ_{i + k - 1})$ , where $k \geq 1$ .

Here we discuss different scenarios of aligning $s_{j}$ with $μ_{i}$ . The first scenario (scenario 1) is that $μ_{i}$ is generated by unmodified 5-mer $s_{j}$ , for which we could compute the likelihood using the Normal distribution defined by $(μ_{s_{j}}^{u n}, σ_{s_{j}}^{u n}) = T_{s_{j}, u n}^{r e f}$ . The second scenario (scenario 2) is the $μ_{i}$ is generated by modified 5-mer $s_{j}$ , that is $N (μ_{s_{j}}^{m o d}, σ_{s_{j}}^{m o d})$ , where $(μ_{s_{j}}^{m o d}, σ_{s_{j}}^{m o d}) = T_{s_{j}, m o d}^{r e f}$ . The third scenario (scenario 3) is that $μ_{i}$ is generated by an unknown nucleotide insertion event, whose probability is given by a pre-defined uniform distribution $p_{u n i f}$ . The fourth scenario (scenario 4) is that the 5-mer $s_{j}$ represents a deletion event on the reference sequence, whose probability is given by the same uniform distribution $p_{u n i f}$ . We define the scoring function as

f (s_{j}, μ_{i}) = {\begin{cases} m a x {l o g (N (μ_{i} ∣ μ_{s_{j}}^{u n}, (σ_{s_{j}}^{u n})^{2})), l o g (N (μ_{i} ∣ μ_{s_{j}}^{m o d}, (σ_{s_{j}}^{m o d})^{2}))}, & (s c e n a r i o 1, 2), \\ l o g (p_{u n i f}) # s_{j} =" - " o r μ_{i} = 0.0 & (s c e n a r i o 3, 4) \end{cases}

We can modify the classical global alignment algorithm to align μ with $s$ . The log probability of aligning $s_{j}$ with $μ_{i}$ could be used as the matching score (Equation 11) in the standard global alignment algorithm. Since one 5-mer $s_{j}$ could be aligned with multiple signal segments $(μ_{i}, μ_{i + 1}, \dots, μ_{i + k - 1})$ , we need to modify the recursion such that $μ_{i + 1}$ could still be aligned with $s_{j}$ even if $μ_{i}$ is already aligned with $s_{j}$ . Therefore, we have the following recursions in our alignment algorithm (Appendix 1—figure 3a). Here we denote the scoring matrix by $M$ ( $m \times n$ ) and the traceback matrix by $T$ , where $| μ | = m - 1$ and $| s | = n - 1$ .

Appendix 1—figure 3

Download asset Open asset

Alignment of signal segments with reference sequence.

M (i, j) = m a x {\begin{cases} M (i - 1, j - 1) + f (μ_{i}, s_{j}), & s c e n a r i o 1, 2 T (i, j) = 1 \\ M (i, j - 1) + f (0, s_{j}), & s c e n a r i o 3 T (i, j) = 2 \\ M (i - 1, j) + f (μ_{i}, -), & s c e n a r i o 4 T (i, j) = 3 \\ M (i - 1, j) + f (μ_{i}, s_{j}), & s c e n a r i o 1, 2 T (i, j) = 4 \end{cases}

where both $μ_{i - 1}$ and $μ_{i}$ are aligned with $s_{j}$ when $T (i, j) = 4$ . This is the only difference between our recursion and the recursion of the standard global alignment algorithm.

Based on the proposed recursion, we developed two alignment algorithms designed for different problems: the full alignment algorithm (‘Full alignment algorithm’) allows insertions and deletions in the reference sequence, while the partial alignment algorithm (‘Partial alignment algorithm’) does not allow any insertions or deletions in the partial reference sequence. The full alignment algorithm tries to align all the current signal segments to the whole reference sequence, while the partial alignment algorithm tries to align all the current signal segments to a consecutive sub-sequence of the 5-mer list from the beginning.

The output of the alignment algorithms is a matching path between the current signal segments μ and the reference 5-mer list $s$ . We denote it by $γ = (μ^{p}, s^{q})$ , where $μ^{p} = (μ_{p_{1}}, μ_{p_{2}}, \dots, μ_{p_{N}})$ and $s^{q} = (s_{q_{1}}, s_{q_{2}}, \dots, s_{q_{N}})$ . The length of the alignment path is given by $| γ | = N$ . $p = (p_{1}, p_{2}, \dots, p_{N})$ are the indexes of the current signal segment μ in the alignment path $γ$ , while $q = (q_{1}, q_{2}, \dots, q_{N})$ are the indexes of the 5-mer list $s$ . Note that there are duplicated values in $p$ and $q$ to allow one 5-mer to match with multiple current signal segments. Appendix 1—figure 3b provides an illustration of the alignment result.

To obtain the modification states $u$ , we will first process the optimal alignment path $γ$ into $\hat{γ}$ (‘Post-processing of alignment path’), from which we can derive the modification states $u$ (‘Modification state estimation on single molecule level’).

Full alignment algorithm

The full alignment algorithm is used when the read may possess differences (insertions or deletions) compared with the reference sequence, which is the general case. The algorithm is described in Appendix 1–algorithm 1.

Appendix 1–algorithm 1. Full alignment algorithm.
Input: μ, s, T^kmer Output: Alignment path $γ$ // Add special symbol before μ and s μ ={0.0, μ}, s ={“–”, s} // Get the length of and μ and s m = \|μ\|, n = \|s\| // Initialize the score matrix M and traceback matrix T Initialize a m × n matrix M, with all elements set to 0.0 Initialize a m × n matrix T, with all elements set to 0 // Initialize the first row of M and T for i ∈ {1, 2, … , n-1} do M(0, i) ← M(0, i - 1) + f (0.0, s_i) T(0, i)←2 end // Intialize the first column of M and T for i ∈ {1, 2, · · · , m − 1} do M(i, 0) ← M(i − 1, 0) + f(µ_i , “ − ”) T(i, 0) ← 3 end // Loop for all matrix for i ∈ {1, 2,···,m − 1} do for j ∈ {1, 2,···, n − 1} do v₁ ← M(i − 1, j − 1) + f(μ_i, s_j) v₂ ← M(i, j − 1) + f(0.0, s_j) v₃ ← M(i − 1, j) + f(μ_i, “ − ”) v₄ ← M(i − 1, j) + f(μ_i, s_j) M(i, j) ← max{v₁, v₂, v₃, v₄} T(i, j) ← argmax{v₁, v₂, v₃, v₄} # index of maximum end end // Traceback Traceback from position (m − 1, n − 1) to position (0, 0) using T to generate the optimal alignment path γ.

Appendix 1–algorithm 1. Full alignment algorithm.

Input: μ, s, T^kmer
Output: Alignment path

γ

// Add special symbol before μ and s
μ ={0.0, μ}, s ={“–”, s}
// Get the length of and μ and s
m = |μ|, n = |s|
// Initialize the score matrix M and traceback matrix T
Initialize a m × n matrix M, with all elements set to 0.0
Initialize a m × n matrix T, with all elements set to 0
// Initialize the first row of M and T
for i ∈ {1, 2, … , n-1} do
M(0, i) ← M(0, i - 1) + f (0.0, s_i)
T(0, i)←2
end
// Intialize the first column of M and T
for i ∈ {1, 2, · · · , m − 1} do
M(i, 0) ← M(i − 1, 0) + f(µ_i , “ − ”)
T(i, 0) ← 3
end
// Loop for all matrix
for i ∈ {1, 2,···,m − 1} do
for j ∈ {1, 2,···, n − 1} do
v₁ ← M(i − 1, j − 1) + f(μ_i, s_j)
v₂ ← M(i, j − 1) + f(0.0, s_j)
v₃ ← M(i − 1, j) + f(μ_i, “ − ”)
v₄ ← M(i − 1, j) + f(μ_i, s_j)
M(i, j) ← max{v₁, v₂, v₃, v₄}
T(i, j) ← argmax{v₁, v₂, v₃, v₄} # index of maximum
end
end
// Traceback
Traceback from position (m − 1, n − 1) to position (0, 0) using T to generate the optimal
alignment path γ.

Partial alignment algorithm

The partial alignment algorithm is used when the read must agree with the reference sequence (no insertions or deletions). If the sequencing data is generated by in vitro transcription from a single template, all reads should agree with the reference sequence. However, the read may only match the reference sequence from the start to a middle point due to the break of the RNA molecule in the reverse transcription.

The partial alignment algorithm is a modification of the full alignment algorithm by considering this special case. Since the alignment may not be full, we traceback from the maximum value in the last row of the scoring matrix μ. This means all current signal segments must be aligned, while only the first part of the 5-mer list needs to be aligned. Appendix 1—figure 3c provides an illustration of the alignment result.

The partial alignment algorithm is described in Algorithm 1–algorithm 2. We can see that scenarios 3 and 4 have been removed from the algorithm compared with the full alignment algorithm, that is $T (i, j) = 2 o r 3$ . In other words, there should be no insertions or deletions in the 5-mer list $s$ .

Appendix 1–algorithm 2. Partial alignment algorithm.
Input: μ, s, T^kmer Output: Alignment path γ // Add special symbol before μ and s μ = {0.0, μ}, s = {“ − ”, s} // Get the length of μ and s m = \|μ\|, n = \|s\| // Initialize the score matrix M and traceback matrix T Initialize a m × n matrix M, with all elements set to 0.0 Initialize a m × n matrix T, with all elements set to 0 // Initialize the first row of M and T for i ∈ {1, 2,···, n − 1} do M(0, i) ← M(0, i − 1) + f(0.0, s_i) T(0, i) ← 2 end // Initialize the first column of M and T for i ∈ {1, 2,···,m − 1} do M(i, 0) ← M(i − 1, 0) + f(µ_i , “ − ”) T(i, 0) ← 3 end // Loop for all matrix for i ∈ {1, 2, · · · ,m − 1} do for j ∈ {1, 2, · · · , n − 1} do v₁ ← M(i − 1, j − 1) + f(μ_i, s_j) v₄ ← M(i − 1, j) + f(μ_i, s_j) M(i, j) ← max{v₁, v₄} if v₁ ≥ v₄ then T(i, j) ← 1; else T(i, j) ← 4; end end end // Traceback Find the maximum value M(m − 1, k) in the last row of the scoring matrix M. Traceback from position (m − 1, k) to position (0, 0) using T to generate the optimal alignment path γ.

Appendix 1–algorithm 2. Partial alignment algorithm.

Input: μ, s, T^kmer
Output: Alignment path γ
// Add special symbol before μ and s
μ = {0.0, μ}, s = {“ − ”, s}
// Get the length of μ and s
m = |μ|, n = |s|
// Initialize the score matrix M and traceback matrix T
Initialize a m × n matrix M, with all elements set to 0.0
Initialize a m × n matrix T, with all elements set to 0
// Initialize the first row of M and T
for i ∈ {1, 2,···, n − 1} do
M(0, i) ← M(0, i − 1) + f(0.0, s_i)
T(0, i) ← 2
end
// Initialize the first column of M and T
for i ∈ {1, 2,···,m − 1} do
M(i, 0) ← M(i − 1, 0) + f(µ_i , “ − ”)
T(i, 0) ← 3
end
// Loop for all matrix
for i ∈ {1, 2, · · · ,m − 1} do
for j ∈ {1, 2, · · · , n − 1} do
v₁ ← M(i − 1, j − 1) + f(μ_i, s_j)
v₄ ← M(i − 1, j) + f(μ_i, s_j)
M(i, j) ← max{v₁, v₄}
if v₁ ≥ v₄ then
T(i, j) ← 1;
else
T(i, j) ← 4;
end
end
end
// Traceback
Find the maximum value M(m − 1, k) in the last row of the scoring matrix M.
Traceback from position (m − 1, k) to position (0, 0) using T to generate the optimal
alignment path γ.

Post-processing of alignment path

In the alignment path $γ$ , one $k$ -mer might be aligned to noise (scenario 4) or one current signal segment may be aligned with ‘-’ (scenario 3). Firstly, we remove such pairs from the alignment path. Then, we re-compute a weighted mean for each remaining $k$ -mer, since one $k$ -mer may be aligned to multiple current signal segments. Here, we combine these segments and calculate the weighted mean over corresponding data points with hidden state 1, that is $z_{i} = 1$ .

Let us assume there are $m$ segments aligned to the same $k$ -mer $s_{i}$ , with $n_{j}$ denotes the number of data points assigned to ‘curr’ state ( $z_{i} = 1$ ) in the j-th segment and $μ_{j}$ denotes the mean of the j-th segment. The weighted mean, denoted by ${\hat{μ}}_{i}$ , corresponding to the $k$ -mer $s_{i}$ , is defined as follows:

{\hat{μ}}_{i} = \sum_{j = 1}^{m} \frac{n_{j}}{n_{1} + n_{2} + \dots + n_{m}} μ_{j} .

In the end, we obtain a one-to-one correspondence between each $k$ -mer and its newly derived weighted mean, denoted by $\hat{μ} = ({\hat{μ}}_{1}, {\hat{μ}}_{2}, \dots)$ . Note that the length of $\hat{μ}$ and $s$ is the same, that is $| \hat{μ} | = | s |$ . We denote the post-processed alignment path by $\hat{γ} = (\hat{μ}, s)$ .

Modification estimation

Modification state estimation on single molecule level

We are interested in estimating the modification state for 5-mer after alignment. $T^{k m e r}$ has provided the mean and std for each 5-mer, as well as its modifications. Note that $T^{k m e r}$ may only contain modification parameters for only a subset of 5-mers. For those 5-mers without modification parameters, their mean and std are set to ‘NA’. We could evaluate the probability density of a given 5-mer under the unmodified distribution or the modified distribution. By comparing the likelihoods, we could infer the modification state. Given the post-processed alignment path $\hat{γ} = (\hat{μ}, s)$ , we estimate the modification state of jth 5-mer $s_{j}$ by

u_{j} = {\begin{cases} ‘ ‘ u n " & i f μ_{s_{j}}^{u n} = ‘ ‘ N A " o r σ_{s_{j}}^{u n} = ‘ ‘ N A ", \\ ‘ ‘ u n " & i f N ({\hat{μ}}_{j} ∣ μ_{s_{j}}^{u n}, (σ_{s_{j}}^{u n})^{2}) \geq N ({\hat{μ}}_{j} ∣ μ_{s_{j}}^{m o d}, (σ_{s_{j}}^{m o d})^{2}) \\ ‘ ‘ m o d " & i f N ({\hat{μ}}_{j} ∣ μ_{s_{j}}^{u n}, (σ_{s_{j}}^{u n})^{2}) < N ({\hat{μ}}_{j} ∣ μ_{s_{j}}^{m o d}, (σ_{s_{j}}^{m o d})^{2}) \end{cases}

where $(μ_{s_{j}}^{u n}, σ_{s_{j}}^{u n}) = T_{s_{j}, u n}^{k m e r}$ and $(μ_{s_{j}}^{m o d}, σ_{s_{j}}^{m o d}) = T_{s_{j}, m o d}^{k m e r}$ .

In the end, we get the modification state vector $u = (u_{1}, u_{2}, \dots, u_{i}, \dots, u_{| \hat{γ} |})$ for each 5-mer in the post-processed alignment path $\hat{γ}$ .

Modification probability estimation on single molecule level

In addition to estimating the modification state, it is also crucial to determine the modification probability of each 5-mer on each read, which is essential for accurately benchmarking the performance at the single-molecule level.

Given the post-processed alignment path $\hat{γ} = (\hat{μ}, s)$ , we estimate the modification probability of jth 5-mer $s_{j}$ by

p_{j} = {\begin{cases} ‘ ‘ N A " & i f μ_{s_{j}}^{u n} = ‘ ‘ N A " o r σ_{s_{j}}^{u n} = ‘ ‘ N A " \\ \frac{p ({\hat{μ}}_{j} ∣ T_{s_{j}, m o d}^{k m e r})}{p ({\hat{μ}}_{j} ∣ T_{s_{j}, m o d}^{k m e r}) + p ({\hat{μ}}_{j} ∣ T_{s_{j}, u n}^{k m e r})}, & i f μ_{s_{j}}^{u n} \neq ‘ ‘ N A " and σ_{s_{j}}^{u n} \neq ‘ ‘ N A " \end{cases}

where $(μ_{s_{j}}^{u n}, σ_{s_{j}}^{u n}) = T_{s_{j}, u n}^{k m e r}$ and $(μ_{s_{j}}^{m o d}, σ_{s_{j}}^{m o d}) = T_{s_{j}, m o d}^{k m e r}$ . Indeed, the $p_{j}$ is the posterior probability given the distribution of modified state:

p_{j} = p (T_{s j, m o d}^{k m e r} ∣ {\hat{μ}}_{j}) = \frac{p ({\hat{μ}}_{j} ∣ T_{s j, m o d}^{k m e r}) \cdot p (T_{s j, m o d}^{k m e r})}{p ({\hat{μ}}_{j} ∣ T_{s j, m o d}^{k m e r}) \cdot p (T_{s j, m o d}^{k m e r}) + p ({\hat{μ}}_{j} ∣ T_{s j, u n}^{k m e r}) \cdot p (T_{s j, u n}^{k m e r})}

where we assume that the prior probabilities are equal, $p (T_{s j, m o d}^{k m e r}) = p (T_{s j, u n}^{k m e r})$ .

In the end, we get the modification probability vector $p = (p_{1}, p_{2}, \dots, p_{i}, \dots, p_{| \hat{γ} |})$ for each 5-mer in the post-processed alignment path $\hat{γ}$ .

Modification rate estimation on site level

In real data analysis, we need to decide the modification rate of a genomic location, which quantifies the probability of a read mapped to this location to be in the modified state. For a genomic location, we denote the 5-mer from reads mapped to this location as $(s_{1}, s_{2}, \dots, s_{N})$ . Note that all $s_{i}, i = 1, 2, \dots, N$ are identical 5-mers on different reads. The modification state of these 5-mers on different reads is denoted by $(u_{1}, u_{2}, \dots, u_{N})$ . The modification rate is given by

\begin{matrix} η = \frac{\sum_{i = 1}^{N} u_{i} = ‘ ‘ m o d "}{N} \end{matrix}

Estimation of k-mer parameters table

This section provides the inference algorithms for 5-mer parameter table $T^{k m e r}$ . This is necessary for training data, as we do not have accurate prior parameters for $T^{k m e r}$ .

Given input data $Y = (y_{1}, y_{2}, \dots, y_{i}, \dots, y_{N})$ , we first perform base calling and mapping to get $S = (s_{1}, s_{2}, \dots, s_{i}, \dots, s_{N})$ run HHMM for all reads such we get $M = (μ_{1}, μ_{2}, \dots, μ_{i}, \dots, μ_{N})$ . We assume there are $n_{i}$ fragments for the ith read $y_{i} = (y_{i}^{(1)}, y_{i}^{(2)}, \dots, y_{i}^{(n_{i})})$ . Similarly, we have $s_{i} = (s_{i}^{(1)}, s_{i}^{(2)}, \dots, s_{i}^{(n_{i})})$ and $μ_{i} = (μ_{i}^{(1)}, μ_{i}^{(2)}, \dots, μ_{i}^{(n_{i})})$ .

The inference algorithm (Appendix 1–algorithm 3) then estimate $T^{k m e r}$ given $S$ , $M$ and $T^{r e f}$ (https://github.com/nanoporetech/kmer_models; Brennen, 2023). As illustrated in Appendix 1—figure 4, the inference algorithm infers the parameters iteratively. With a given $T^{k m e r}$ , we can use it to obtain the pairing between raw current signal segments and 5-mers, as well as modification states, through the alignment algorithm (‘Alignment of signal segments with reference sequence’). Given this pairing, we extract all current signal segments for a specific 5-mer, which may correspond to different genomic locations and different reads. As there are two states for this 5-mer (unmodified and modified), we expect to observe two peaks in the density of the current signal segments. Therefore, we could fit a two-component Gaussian mixture model to the current signal segments to estimate relevant parameters of the given 5-mer in $T^{k m e r}$ .

Appendix 1—figure 4

Download asset Open asset

Estimation of 5-mer parameters.

$T^{k m e r}$ is initialized using $T^{r e f}$ . In training, our algorithm infers the parameters iteratively. After training, $T^{r e f}$ is fixed in testing.

Appendix 1–algorithm 3. Inference algorithm of k-mer parameter table.
Input: S, M, T^ref Output:T^kmer // Initialize T^kmer, ∀s ∈ {AAAAA,AAAAC, · · · , TTTTT}, ∀u ∈ {“un”, “mod”} for (μ_s,u, σ_s,u) ∈ T^kmer do (μ_s,u, σ_s,u) = T^ref_s,“un” end // iterate T times to update T^kmer for round = 1, 2 · · · T do // align all fragments of all reads Perform full alignment algorithm (Appendix 1–algorithm 1) for μ_i^(k) and s_i(k) to get the alignment path γ_i^(k) , where i = 1 · · · N and k = 1 · · · n_i. Post-process (Post-processing of alignment path) the alignment path γ_i^(k) into ${\hat{γ}}_{i}^{(k)}$ = $[({\hat{μ}}_{1}^{(k)}, s_{1}^{(k)}), ({\hat{μ}}_{2}^{(k)}, s_{2}^{(k)}), \dots]$ . Concatenate all ${\hat{γ}}_{i}^{(k)}$ and re-index into one vector $\hat{Γ} = [(\hat{μ}, s_{1}), (\hat{μ}, s_{2}), \dots]$ . // parameter estimation for s0 ∈ {AAAAA,AAAAC, · · · , TTTTT} do Collect all ${\hat{μ}}_{j}$ in $\hat{Γ}$ where sj = s0 into one vector $x = {{\hat{μ}}_{j}, \forall s_{j} = s_{0}, \forall ({\hat{μ}}_{j}, s_{j})} \in \hat{Γ}$ . Fit x using a two-component Gaussian Mixture Model with the first component mean fixed to $μ_{2_{0}}$ , where ( $μ_{2_{0}}$ , $σ_{2_{0}}$ ) = $T_{s_{0}, ‘ ‘ u n "}^{r e f}$ , and obtain the std of the first component σ₁, the mean and std of the second component μ₂, σ₂. Update the parameters of 5-mer s0 in $T^{k m e r}$ : $T_{s_{0}, ‘ ‘ u n "}^{k m e r} \leftarrow (μ_{s_{0}, σ_{1}})$ and $T_{s_{0}, ``mod"}^{k m e r} \leftarrow (μ_{2}, σ_{2})$ end Update T^kmer manually based on several heuristic criteria. end

Appendix 1–algorithm 3. Inference algorithm of k-mer parameter table.

Input: S, M, T^ref
Output:T^kmer
// Initialize T^kmer, ∀s ∈ {AAAAA,AAAAC, · · · , TTTTT}, ∀u ∈ {“un”, “mod”}
for (μ_s,u, σ_s,u) ∈ T^kmer do
(μ_s,u, σ_s,u) = T^ref_s,“un”
end
// iterate T times to update T^kmer
for round = 1, 2 · · · T do
// align all fragments of all reads
Perform full alignment algorithm (Appendix 1–algorithm 1) for μ_i^(k) and s_i(k) to get the alignment path γ_i^(k) ,
where i = 1 · · · N and k = 1 · · · n_i.
Post-process (Post-processing of alignment path) the alignment path γ_i^(k) into

{\hat{γ}}_{i}^{(k)}

[({\hat{μ}}_{1}^{(k)}, s_{1}^{(k)}), ({\hat{μ}}_{2}^{(k)}, s_{2}^{(k)}), \dots]

.
Concatenate all

{\hat{γ}}_{i}^{(k)}

and re-index into one vector

\hat{Γ} = [(\hat{μ}, s_{1}), (\hat{μ}, s_{2}), \dots]

.
// parameter estimation
for s0 ∈ {AAAAA,AAAAC, · · · , TTTTT} do
Collect all

{\hat{μ}}_{j}

\hat{Γ}

where sj = s0 into one vector

x = {{\hat{μ}}_{j}, \forall s_{j} = s_{0}, \forall ({\hat{μ}}_{j}, s_{j})} \in \hat{Γ}

.
Fit x using a two-component Gaussian Mixture Model with the first component
mean fixed to

μ_{2_{0}}

, where (

μ_{2_{0}}

σ_{2_{0}}

) =

T_{s_{0}, ‘ ‘ u n "}^{r e f}

, and obtain the std of the first
component σ₁, the mean and std of the second component μ₂, σ₂.
Update the parameters of 5-mer s0 in

T^{k m e r}

T_{s_{0}, ‘ ‘ u n "}^{k m e r} \leftarrow (μ_{s_{0}, σ_{1}})

and

T_{s_{0}, ``mod"}^{k m e r} \leftarrow (μ_{2}, σ_{2})

end
Update T^kmer manually based on several heuristic criteria.
end

In practice, the following heuristic criteria are used: (1) for a given 5-mer, perform GMM for each genomic location with the same 5-mer and only pool those sites with high modification rate (2) there must be two peaks in the density plot of a given 5-mer; (3) the mean of the two peaks are far apart; (4) the weight of second component (modified state) must be large, e.g. >0.1; (5) the weight of the first component (unmodified state) should be larger than the second component (modified state), as majority nucleotides are unmodified; (6) limit the size of std for both components in the GMM ( $σ_{1} < 5$ , $σ_{2} < 5$ ) to comply with the physical settings.

In the end, we will get an updated $T^{k m e r}$ from the training data, which contains the modification parameters for the strongest 5-mers. As illustrated in Appendix 1—table 1, we do have parameters for the modification state for certain $k$ -mers such as ‘AAACA’, ‘AAACT’ etc.

Evaluation of segmentation and event alignment

The output of SegPore eventalign is the alignment path γ = (μ, s), which represents a one-to-one correspondence between the mean list μ = (μ₁, μ₂,...) and the k-mer lists = (s₁, s₂, ...). In addition, the corresponding standard deviation list σ = (σ₁, σ₂, ...) is also obtained along the alignment path. Appendix 1—table 3 is an example of the SegPore eventalign output, which illustrates one alignment path. The column ref_kmer corresponds to s, the column mean corresponds to μ, and the column stdv corresponds to σ.

Appendix 1—table 3

SegPore eventalign output example.

read_idx	read_name	contig	pos	strand	ref_kmer	model_kmer	mean	stdv	start_idx	end_idx	event_len
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	57	+	GGTGT	GGTGT	79.339	2.414	47296	47321	18
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	58	+	GTGTC	GTGTC	89.809	3.579	47276	47292	14
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	59	+	TGTCT	TGTCT	117.017	6.287	47251	47272	18
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	60	+	GTCTT	GTCTT	76.448	1.038	47148	47247	87
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	61	+	TCTTA	TCTTA	76.036	1.044	47117	47144	25
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	62	+	CTTAG	CTTAG	78.886	1.127	47092	47113	21
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	63	+	TTAGT	TTAGT	92.971	2.307	47007	47088	73
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	64	+	TAGTG	TAGTG	121.129	4.2	46882	46983	90
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	65	+	AGTGT	AGTGT	90.22	4.554	46822	46878	33
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	67	+	TGTGC	TGTGC	104.566	2.563	46798	46818	8
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	68	+	GTGCT	GTGCT	91.513	2.51	46778	46794	15
80741	003655f4-ff77-409d-a29e-20dc778d1c47	A1	69	+	TGCTT	TGCTT	111.485	4.001	46754	46774	10

Given the eventalign results, we use two metrics to evaluate the performance of the segmentation and event alignment: (1) the average std $\hat{σ}$ , and (2) the average log-likelihood $\hat{L}$ . Assuming there are $N$ reads in total (each read has one path), and the nth read has $K_{n}$ events ( $| γ_{n} | = K_{n}$ ). The average std $\hat{σ}$ is defined as

\hat{σ} = \frac{1}{N} \sum_{n = 1}^{N} {\frac{1}{K_{n}} \sum_{k = 1}^{K_{n}} σ_{k}},

and the average log-likelihood $\hat{L}$ is defined as

\hat{L} = \frac{1}{N} \sum_{n = 1}^{N} {\frac{1}{K_{n}} \sum_{k = 1}^{K_{n}} \log N (μ_{k} | μ_{s_{k}}^{r e f}, σ_{s_{k}}^{r e f})}

where $(μ_{s_{k}}^{r e f}, σ_{s_{k}}^{r e f}) = T_{s_{k}}^{r e f}$ .

As shown in Appendix 1—figure 5, the red line represents the event mean $μ_{k}$ and the shaded area represents the std $σ_{s_{k}}^{e s t}$ for event $k$ . A poorly segmented raw signal corresponding to an event will exhibit a large standard deviation. Therefore, a smaller $\hat{σ}$ indicates lower variations within the raw signal segment, signifying better performance.

Appendix 1—figure 5

Download asset Open asset

Illustration of segmentation and event alignment.

If an event is aligned to the correct reference $k$ -mer, the mean $μ_{k}$ will be close to the reference $μ_{s_{k}}^{r e f}$ and the log-likelihood $\hat{L}$ will be large. So higher $\hat{L}$ means more similar results to ONT’s estimates and better performances.

Appendix 2

Introduction

We have introduced the hierarchical hidden Markov model (HHMM) in Appendix 1. The HHMM is a two-layer HMM, consisting of the outer HMM and the inner HMM. The outer HMM specifies the base blocks and transition blocks. The inner HMM specifies the base blocks. Linear models are used for the transition blocks. However, we have not provided the details for the parameter inference. This document aims to provide a comprehensive introduction to the inference algorithm.

As we know, the raw fast5 files are very large, which contain millions of reads. Each read has 20,000 measurements in the raw current signal. This poses a significant computational burden to the HHMM inference algorithm. Therefore, we used GPU to accelerate the inference algorithm.

We will first introduce the general inference algorithm, then describe the implementation details.

For simplicity, we assume the input is a raw current signal of a single read with $N$ measurements $y = (y_{1}, y_{2}, \dots, y_{N})$ . The output is the hidden states of the outer HMM $g = (g_{1}, g_{2}, \dots, g_{N})$ and the hidden states of the inner HMM $h = (h_{1}, h_{2}, \dots, h_{N})$ , where $g_{i} \in {‘ ‘ B ", ‘ ‘ T "}$ and $h_{i} \in {‘ ‘ c u r r ", ‘ ‘ p r e v ", ‘ ‘ n e x t ", ‘ ‘ n o i s e "}$ .

Given a hidden state configuration of the outer HMM $g$ , it partitions the raw current signal $y$ into alternating base blocks and transition blocks $c = (c_{1}, c_{2}, \dots, c_{i}, \dots, c_{2 K + 1})$ , where $c_{i} \in {‘ ‘ B ", ‘ ‘ T "}$ denotes the label of ith block. Note that $c_{i} = ‘ ‘ B " \forall i = 1, 3, \dots, 2 K + 1$ and $c_{i} = ‘ ‘ T " \forall i = 2, 4, \dots, 2 K$ . Now we can denote $y = (y^{(1)}, y^{(2)}, \dots, y^{(2 K + 1)})$ .

The general idea of the inference algorithm is that we will calculate the likelihood for all possible configurations of the hidden states of the outer HMM $g$ and choose the one with the largest likelihood as our final estimation. Due to the special constraints of $g$ , we can enumerate all possible configurations (‘Enumeration of hidden states of outer HMM’). Here we denote the full configuration set of $g$ by $G$ . The general inference algorithm is discussed in detail in ‘General inference algorithm’.

The general inference algorithm provides the analytic forms for calculating the exact likelihood. However, a direct implementation of the algorithm (Python and C++) requires huge amounts of computation resources and time to handle a typical fast5 file. Therefore, we implement a GPU version of the inference algorithm (‘Implementation details’), in which base blocks from different configurations are simultaneously inferred on different GPU cores.

Enumeration of hidden states of outer HMM

This section describes how to get the full configuration set $G = (g_{1}, g_{2}, \dots, g_{i}, \dots, g_{M})$ given the raw signal $y$ . Obviously, it is not possible to enumerate all hidden state $g_{i} = (g_{i 1}, g_{i 2}, \dots, g_{i j}, \dots, g_{i N})$ (where $g_{i j} \in {‘ ‘ B ", ‘ ‘ T "}$ ) because there are $2^{N}$ configurations and $N$ can be larger than 20,000. However, there exist special properties of the raw signal $y$ : (1) the signal should be in alternating base and transition blocks (2) the hidden state $g_{i}$ does not change within a block (3) the average length of base blocks is larger than that of the transition blocks (4) the variation of base blocks is smaller than that of transition blocks. Based on these data properties, we could partition the raw signal $y$ into atomic level base and transition blocks (‘Atomic block initialization’). Multiple consecutive atomic level base and transition blocks can be merged into a single larger base block. As shown in Appendix 2—figure 1a, we partition the raw signal $y$ into 15 atomic-level blocks $(B 1, T 1, B 2, T 2, \dots, B 7, T 7, B 8)$ . We can merge $(B 1, T 1, B 2)$ into a new base block or $(B 5, T 5, B 6, T 6, B 7)$ into a new base block (Appendix 2—figure 1b).

Given a pre-defined partition $g_{0}$ of raw signal $y$ , we can derive the $2 K + 1$ atomic blocks

$(B_{1}, T_{1}, B_{2}, T 2, \dots, B_{K}, T_{K}, B_{K + 1}))$ . To facilitate the description of the enumeration algorithm, we define $s e g (i, j) = (B_{i + 1}, T_{i + 1}, B_{i + 2}, \dots, T_{j - 1}, B_{j})$ to be the merged block of all blocks in $s e g (i, j)$ , where $0 \leq i < j \leq K + 1$ . A path is then a list of non-overlapping and consecutive blocks. For example, the path corresponding to $g_{i}$ in Appendix 2—figure 1b is given by $[s e g (0, 2), T_{2}, B_{3}, T_{3}, B_{4}, T_{4}, s e g (4, 7), T_{7}, B_{8}]$ . To reduce the number of paths, we set the maximum number of base blocks in any merged block to $N_{m a x}$ . The enumeration algorithm is given by Appendix 2–algorithm 1.

Appendix 2–algorithm 1. Enumeration of hidden states of outer HMM.
Input: g₀, N_max, K Output: $G$ Derive (B₁, T₁, B₂, T₂,··· ,B_K, T_K, B_K+1)) from g₀ $G \leftarrow {}$ for path in recur_enum(0) do Convert path to g Append g to $G$ end // Recursive function to enumerate all paths starting from position st Function recur_enum(st): res = [] if st = K then return [B_K+1] end // “[” is inclusive and “)” is exclusive for i ∈ [st, K + 1) do if i − st > N_max then return res end tmp_seg = seg(st, i) all_paths_from_i = recur_enum(i) for path in all_ paths_ from_i do tmp_path = concatenate tmp_seg with path Append tmp_path to res end end return res;

Appendix 2–algorithm 1. Enumeration of hidden states of outer HMM.

Input: g₀, N_max, K
Output:

G

Derive (B₁, T₁, B₂, T₂,··· ,B_K, T_K, B_K+1)) from g₀

G \leftarrow {}

for path in recur_enum(0) do
Convert path to g
Append g to

G

end
// Recursive function to enumerate all paths starting from position st
Function recur_enum(st):
res = []
if st = K then
return [B_K+1]
end
// “[” is inclusive and “)” is exclusive
for i ∈ [st, K + 1) do
if i − st > N_max then
return res
end
tmp_seg = seg(st, i)
all_paths_from_i = recur_enum(i)
for path in all_ paths_ from_i do
tmp_path = concatenate tmp_seg with path
Append tmp_path to res
end
end
return res;

General inference algorithm

Given all possible hidden state configurations of the outer HMM $G$ , we want to calculate the likelihood for all $g \in G$ . In order to calculate the likelihood, we need to specify the parameters for the emission distributions.

Given $g$ , we can derive that $y$ is divided into $2 K + 1$ alternating base and transition blocks. We denote all blocks by $y = (y^{(1)}, y^{(2)}, \dots, y^{(i)}, \dots, y^{(2 K + 1)})$ , where $y^{(i)}$ represent the ith block. The base blocks are given by $(y^{(b_{1})}, y^{(b_{2})}, \dots, y^{(b_{i})}, \dots, y^{(b_{K + 1})})$ , where $b_{i} = 2 i - 1$ and $i \in {1, 2, \dots, K + 1}$ . The transition blocks are given by $(y^{(t_{1})}, y^{(t_{2})}, \dots, y^{(t_{i})}, \dots, y^{(t_{K})})$ , where $t_{i} = 2 i$ and $i \in {1, 2, \dots, K}$ .

We pre-specify the transition probabilities $T^{o u t e r} = {T_{B B}^{o u t e r}, T_{B T}^{o u t e r}, T_{T B}^{o u t e r}, T_{T T}^{o u t e r}}$ between the base state and the transition state in the outer HMM, as shown in Appendix 2—table 1. The pre-specification has two considerations: (1) the probability of a transition from B to T after 10 consecutive B states is 0.1; (2) the probability of a transition from T to B after 10 consecutive T states is larger than 0.9; (3) to reduce computation load in the parameter inference. We pre-specify the initial probability of the first hidden state to 1, since $g$ always starts with a ‘B’ state.

Appendix 2—table 1

Transition parameters of the outer HMM ( $T^{o u t e r}$ ).

	‘B’	‘T’
‘B’	0.99	0.01
‘T’	0.10	0.90

For inner HMM, we pre-specify the transition probabilities $T^{i n n e r}$ shown in Appendix 2—table 2, which is based on the following consideration (1) large proportion of the data should come from the ‘curr’ state, that is it is easy to transit from other states to ‘curr’ state; (2) it is unlikely to transit from ‘prev’ to ‘next’, and vice versa; (3) it is unlikely to transit from other states to the ‘noise’ state. We pre-specify the initial probability of the first hidden state to 0.25, which gives equal probabilities for all four states. We pre-specify the parameters of the ‘noise’ state ( $U n i f (l b, u b)$ ) to $l b = 50$ and $u b = 130$ , which are estimated from the current signals across different datasets to represent a reasonable range.

Appendix 2—table 2

Transition parameters of the inner HMM ( $T^{i n n e r}$ ).

	‘curr’	‘prev’	‘next’	‘noise’
‘curr’	0.925	0.025	0.025	0.025
‘prev’	0.300	0.500	0.100	0.100
‘next’	0.300	0.100	0.500	0.100
‘noise’	0.300	0.100	0.100	0.500

The remaining parameters, mainly associated with emission distributions, need to be estimated. For base blocks, we have inner HMM parameters $ϕ = (ϕ^{(b_{1})}, ϕ^{(b_{2})}, \dots, ϕ^{(b_{K + 1})})$ and $σ = (σ^{(b_{1})}, σ^{(b_{2})}, \dots, σ^{(b_{K + 1})})$ . For transition blocks, we have linear coefficients $β_{0} = (β_{0}^{(t_{1})}, β_{0}^{(t_{2})}, \dots, β_{0}^{(t_{K})})$ , $β_{1} = (β_{1}^{(t_{1})}, β_{1}^{(t_{2})}, \dots, β_{1}^{(t_{K})})$ and noise std $σ_{ϵ} = (σ_{ϵ}^{(t_{1})}, σ_{ϵ}^{(t_{2})}, \dots, σ_{ϵ}^{(t_{K})})$ .

The general idea of the inference algorithm is to calculate the joint likelihood of the outer HMM given for any $g$ . However, it is not possible to calculate the emission probability for each measurement as there exist dependencies with a base or transition block. Here we first calculate the joint likelihood for a base block (inner HMM) or a transition block (linear model), then average the joint likelihood to each measurement to use it as substitute for the emission probability of each measurement in the outer HMM. Here we use $ω = (ω_{1}, ω_{2}, \dots, ω_{N})$ to denote the approximated emission probability of each measurement. Note that $ω_{i}$ is derived either based on the inner HMM or the linear model, and $| ω | = | y | = | g | = | h |$ .

The general inference algorithm is given by Appendix 2–algorithm 2.

Appendix 2–algorithm 2. General inference algorithm for HHMM.
Input: y, $G$ , T^outer, π^outer, T^inner, π^inner, Output: best path $\tilde{g}$ of outer HMM, best path $\tilde{h}$ of inner HMM, ϕ, σ, β₀, β₁, *σ_ϵ* // Initialization set all elements of $\tilde{ω}, \tilde{g} a n d \tilde{h}$ to 0, and $∣ \tilde{g} ∣=∣ \tilde{h} ∣=∣ y ∣$ $H \leftarrow {}$ Γ ← {} // Calculate likelihood for each g for g ∈ $G$ do $h \to \tilde{h}$ Divide y into base and transition blocks (y⁽¹⁾, y⁽²⁾, · · · , y^(2K+1)) based on g // transition state parameter estimation Estimate ${\hat{ω}}^{(t_{i})}, {\hat{β}}_{0}^{(t_{i})}, {\hat{β}}_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})}$ using linear model (Appendix 2–algorithm 3), where t_i = 2i and i ∈ {1, 2,···, K} Update relevant part of ω using ${\hat{ω}}^{(t_{i})}$ . // base state parameter estimation Initialize ϕ^(bi) and σ⁽^bi⁾ use the empirical mean and std of y^(bi), where b_i = 2i − 1 and i ∈ {1, 2,···, K + 1} ϕ^new ← ϕ σ^new ← σ for round = 1, 2,···, R do for i ∈ (1, 2,···, K + 1) do b_i ← 2i − 1 // inner HMM parameter estimation Estimate ${\hat{ω}}^{(b_{i})}, {\hat{h}}^{(b_{i})}, {\hat{ϕ}}^{(b_{i})}, {\hat{σ}}^{(b_{i})}$ using Forward-backward algorithm (Appendix 2–algorithm 4) given $y^{(b_{i})}, ϕ^{(b_{i})}, σ^{(b_{i})}, ϕ^{(b_{i - 1})}, σ^{(b_{i - 1})}, ϕ^{(b_{i + 1})}, σ^{(b_{i + 1})}, T^{inner}, π^{inner}$ Update corresponding elements of ϕ^new and σ^new using ${\hat{ϕ}}^{(b_{i})}, {\hat{σ}}^{(b_{i})}$ Update relevant part of h and ω using ${\hat{h}}^{(b_{i})}$ and ${\hat{ω}}^{(b_{i})}$ end ϕ ← ϕ^new σ ← σ^new end Calculate the joint likelihood of the outer HMM using $g, ω, T^{o u t e r}, π^{o u t e r}$ $γ = π_{g_{1}}^{outer} \prod_{i = 2}^{N} T_{g_{i - 1} g_{i}}^{outer} \prod_{j = 1}^{N} ω_{j}$ Append h and γ to $H$ and Γ, respectively. end Choose the largest likelihood from Γ and retrieve the corresponding g and h from $G$ and $H$ , respectively. $\tilde{g} \leftarrow g$ $\tilde{h} \leftarrow h$

Appendix 2–algorithm 2. General inference algorithm for HHMM.

Input: y,

G

, T^outer, π^outer, T^inner, π^inner,
Output: best path

\tilde{g}

of outer HMM, best path

\tilde{h}

of inner HMM, ϕ, σ, β₀, β₁, σ_ϵ
// Initialization
set all elements of

\tilde{ω}, \tilde{g} a n d \tilde{h}

to 0, and

∣ \tilde{g} ∣=∣ \tilde{h} ∣=∣ y ∣

H \leftarrow {}

Γ ← {}
// Calculate likelihood for each g
for g ∈

G

h \to \tilde{h}

Divide y into base and transition blocks (y⁽¹⁾, y⁽²⁾, · · · , y^(2K+1)) based on g
// transition state parameter estimation
Estimate

{\hat{ω}}^{(t_{i})}, {\hat{β}}_{0}^{(t_{i})}, {\hat{β}}_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})}

using linear model (Appendix 2–algorithm 3),
where t_i = 2i and i ∈ {1, 2,···, K}
Update relevant part of ω using

{\hat{ω}}^{(t_{i})}

.
// base state parameter estimation
Initialize ϕ^(bi) and σ⁽^bi⁾ use the empirical mean and std of y^(bi), where b_i = 2i − 1 and
i ∈ {1, 2,···, K + 1}
ϕ^new ← ϕ
σ^new ← σ
for round = 1, 2,···, R do
for i ∈ (1, 2,···, K + 1) do
b_i ← 2i − 1
// inner HMM parameter estimation
Estimate

{\hat{ω}}^{(b_{i})}, {\hat{h}}^{(b_{i})}, {\hat{ϕ}}^{(b_{i})}, {\hat{σ}}^{(b_{i})}

using Forward-backward
algorithm (Appendix 2–algorithm 4) given

y^{(b_{i})}, ϕ^{(b_{i})}, σ^{(b_{i})}, ϕ^{(b_{i - 1})}, σ^{(b_{i - 1})}, ϕ^{(b_{i + 1})}, σ^{(b_{i + 1})}, T^{inner}, π^{inner}

Update corresponding elements of ϕ^new and σ^new using

{\hat{ϕ}}^{(b_{i})}, {\hat{σ}}^{(b_{i})}

Update relevant part of h and ω using

{\hat{h}}^{(b_{i})}

and

{\hat{ω}}^{(b_{i})}

end
ϕ ← ϕ^new
σ ← σ^new
end
Calculate the joint likelihood of the outer HMM using

g, ω, T^{o u t e r}, π^{o u t e r}

γ = π_{g_{1}}^{outer} \prod_{i = 2}^{N} T_{g_{i - 1} g_{i}}^{outer} \prod_{j = 1}^{N} ω_{j}

Append h and γ to

H

and Γ, respectively.
end
Choose the largest likelihood from Γ and retrieve the corresponding g and h from

G

and

H

,
respectively.

\tilde{g} \leftarrow g

\tilde{h} \leftarrow h

Implementation details

Atomic block initialization

As mentioned in ‘Enumeration of hidden states of outer HMM’, we need to partition the raw signal $y$ into $2 K + 1$ atomic blocks

$(B_{1}, T_{1}, B_{2}, T 2, \dots, B_{K}, T_{K}, B_{K + 1})$ . If we can identify all atomic transition blocks $(T_{1}, T_{2}, \dots, T_{K})$ , then the atomic base blocks are automatically determined.

We assume atomic transition blocks $(T_{1}, T_{2}, \dots, T_{K})$ have higher variations than atomic base blocks. If we fit a line to a transition block, the absolute value of the slope of the line should be large. Based on this idea, we can obtain the atomic transition blocks as follows:

For each measurement $i$ , fit another line to $(y_{i - 1}, y_{i}, y_{i + 1})$ and obtain the slope $β_{1}$ , fit a line to $(y_{i - 2}, y_{i - 1}, y_{i}, y_{i + 1}, y_{i + 2})$ and obtain the slope $β_{2}$ , use $β_{i} = | 0.5 * (β_{1} + β_{2}) |$ as the absolute slope for $y_{i}$ .
Smooth the obtain slopes $β = (β_{1}, β_{2}, \dots, β_{N})$ by taking the mean of a sliding window of size. Now we get the smoothed $\bar{β}$
Find peaks from $\bar{β}$ using the find_peaks() function in the Python Scipy.signal package. Here, set the minimum distance of detected peaks to 10. Now we get a list of peak positions $θ = (θ_{1}, θ_{2}, \dots, θ_{K})$ .
Expand $θ_{k}$ to the atomic transition block $T_{k}$ . We fit a 3-segment spline to

$(y_{θ_{k} - 8}, \dots, y_{θ_{k} - l}, \dots, y_{θ_{k}}, \dots, y_{θ_{k} + r}, \dots, y_{θ_{k} + 8})$ , where $0 < l, r < 8$ . The slopes of the first segment $(y_{θ_{k} - 8}, \dots, y_{θ_{k} - l - 1})$ and the third segment $(y_{θ_{k} + r + 1}, \dots, y_{θ_{k} + 8})$ are set to 0, while a standard linear model is fit to the second segment $(y_{θ_{k} - l}, \dots, y_{θ_{k}}, \dots, y_{θ_{k} + r})$ . In the fitting, the slope of the second segment should be larger than 5.0. After the fitting of the spline, we take the residuals and calculate the joint likelihood assuming each residual follows a Gaussian distribution $N (0, 4)$ . By enumerating all $l, r$ , we pick the $\hat{l}, \hat{r}$ with the largest joint likelihood. Therefore, the final kth transition block $T_{k}$ corresponding to the interval $[θ_{k} - l, θ_{k} + r]$ .

Emission probabilities of linear model

For the ith transition block, the measurements are given by $y^{(t_{i})} = (y_{1}^{(t_{i})}, y_{2}^{(t_{i})}, \dots, y_{N_{t_{i}}}^{(t_{i})})$ , where $t_{i} = 2 i$ and $i \in 1, 2, \dots, K$ . We denote the corresponding index of $y^{(t_{i})}$ by $x^{(t_{i})} = (1, 2, \dots, N_{t_{i}})$ . Here we want to get the emission probabilities for each measurement ${\hat{ω}}^{(t_{i})}$ , which can be used by the outer HMM. Detailed steps are provided in the following Appendix 2–algorithm 3.

Appendix 2–algorithm 3. Emission probabilities of linear model.
Input: $y^{(t_{i})}$ , $x^{(t_{i})}$ Output: ${\hat{ω}}^{(t_{i})}$ , $β_{0}^{(t_{i})}$ , $β_{1}^{(t_{i})}$ , $σ_{ϵ}^{(t_{i})}$ initialize ${\hat{ω}}^{(t_{i})}$ to a 1 × N_ti all zero vector. Fit a linear model of y^(ti) versus x^(ti), obtain the maximum likelihood estimates $β_{0}^{(t_{i})}, β_{1}^{(t_{i})}$ and $σ_{ϵ}^{(t_{i})}$ Calculate the maximum log joint likelihood $\log p (y^{(t_{i})} \| x^{(t_{i})}, β_{0}^{(t_{i})}, β_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})})$ = $\log N (y^{(t_{i})} - β_{1}^{(t_{i})} x^{(t_{i})} - β_{0}^{(t_{i})} ∣ 0, (σ_{ϵ}^{(t_{i})})^{2})$ Set each element of ${\hat{ω}}^{t_{i}}$ to exp $(\frac{1}{N_{t_{i}}} \log p (y^{(t_{i})} ∣ x^{(t_{i})}, β_{0}^{(t_{i})}, β_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})}))$

Appendix 2–algorithm 3. Emission probabilities of linear model.

Input:

y^{(t_{i})}

x^{(t_{i})}

Output:

{\hat{ω}}^{(t_{i})}

β_{0}^{(t_{i})}

β_{1}^{(t_{i})}

σ_{ϵ}^{(t_{i})}

initialize

{\hat{ω}}^{(t_{i})}

to a 1 × N_ti all zero vector.
Fit a linear model of y^(ti) versus x^(ti), obtain the maximum likelihood estimates

β_{0}^{(t_{i})}, β_{1}^{(t_{i})}

and

σ_{ϵ}^{(t_{i})}

Calculate the maximum log joint likelihood

\log p (y^{(t_{i})} | x^{(t_{i})}, β_{0}^{(t_{i})}, β_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})})

\log N (y^{(t_{i})} - β_{1}^{(t_{i})} x^{(t_{i})} - β_{0}^{(t_{i})} ∣ 0, (σ_{ϵ}^{(t_{i})})^{2})

Set each element of

{\hat{ω}}^{t_{i}}

to exp

(\frac{1}{N_{t_{i}}} \log p (y^{(t_{i})} ∣ x^{(t_{i})}, β_{0}^{(t_{i})}, β_{1}^{(t_{i})}, σ_{ϵ}^{(t_{i})}))

Emission probabilities of inner HMM

This section describes how to obtain the approximated emission probabilities of the outer HMM for a the ith base block $y^{(b_{i})} = (y_{1}^{(b_{i})}, y_{2}^{(b_{i})}, \dots, y_{N_{b_{i}}}^{(b_{i})})$ , which is modeled using an inner HMM. Here $b_{i} = 2 i - 1$ and $i \in {1, 2, \dots, K}$ .

For the inner HMM, we list all the parameters in Appendix 2—table 3.

Appendix 2—table 3

Parameters of the inner HMM for the $b_{i}$ th block.

Parameter	Fixed parameter	Description
$T^{i n n e r}$	Yes	Transition probabilities between hidden states (‘curr’, ‘prev’, ‘next’, ‘noise’), specified in Appendix 2—table 2
$π_{k}^{i n n e r}$	Yes	probability of the first hidden state, $π_{k}^{i n n e r} = 0.25$ , where $k \in {1, 2, 3, 4}$
$l b$	Yes	lower bound of uniform distribution, $l b = 50$
$u b$	Yes	upper bound of uniform distribution, $u b = 130$
$ϕ^{(b_{i})}$	No	Gaussian mean ofith base block
$σ^{(b_{i})}$	No	Gaussian std ofi th base block
$ϕ^{(b_{i - 1})}$	No	Gaussian mean of $i - 1$ th base block
$σ^{(b_{i - 1})}$	No	Gaussian std of $i - 1$ th base block
$ϕ^{(b_{i + 1})}$	No	Gaussian mean of $i + 1$ th base block
$σ^{(b_{i + 1})}$	No	Gaussian std of $i + 1$ th base block

We want to obtain the approximated emission probabilities ${\hat{ω}}^{(b_{i})}$ , hidden states ${\hat{h}}^{(b_{i})}$ , estimated parameters ${\hat{ϕ}}^{(b_{i})}$ and ${\hat{σ}}^{(b_{i})}$ for the ith base block $y^{(b_{i})}$ . To achieve this, we follow the Forward-backward algorithm in Bishop and Nasrabadi, 2006 (Chapter 13.2). The general idea of the Forward-backward algorithm is to use a $Q$ -function to approximate the marginal likelihood, then use an expectation-maximization (EM) algorithm for the parameter estimation.

We introduce the following notations for the Forward-backward algorithm. We use $Z = (z_{1}, z_{2}, \dots, z_{N_{b_{i}}})$ to denote the hidden states, where $Z$ is a $N_{b_{i}} \times 4$ matrix, $z_{j} = (z_{j 1}, z_{j 2}, z_{j 3}, z_{j 4})$ and $z_{j k} \in {0, 1}$ . If $z_{j k}$ is 1, then it means $h_{j}$ corresponds to the kth hidden state. We use $α$ ( $N_{b_{i}} \times 4$ ) to denote the forward probabilities and $β$ ( $N_{b_{i}} \times 4$ ) to denote the backward probabilities. $γ$ ( $N_{b_{i}} \times 4$ ) denotes the posterior distribution calculated from $α$ and $β$ . $ξ$ ( $(N_{b_{i}} - 1) \times 4 \times 4$ ) denotes the joint distribution of two successive hidden states. We use $θ = (ϕ^{(b_{i})}, σ^{(b_{i})})$ to denote all parameters to be estimated.

\begin{array}{ll} Q (θ, θ^{old}) & = \sum_{Z} p (Z ∣ y^{(b_{i})}, θ^{old}) \log p (y^{(b_{i})}, Z ∣ θ) \\ = \sum_{k = 1}^{4} γ_{1 k} \log π_{k}^{i n n e r} + \sum_{n = 2}^{N_{b_{i}}} \sum_{j = 1}^{4} \sum_{k = 1}^{4} ξ_{n - 1, j, k} \log T_{j k}^{i n n e r} \\ + \sum_{n = 1}^{N_{b_{i}}} {γ_{n 1} \log p (y_{n}^{(b_{i})} ∣ ϕ^{(b_{i})}, σ^{(b_{i})}) + γ_{n 2} \log p (y_{n}^{(b_{i})} ∣ ϕ^{(b_{i - 1})}, σ^{(b_{i - 1})}) \\ + γ_{n 3} \log p (y_{n}^{(b_{i})} ∣ ϕ^{(b_{i + 1})}, σ^{(b_{i + 1})}) + γ_{n 4} \log p (y_{n}^{(b_{i})} ∣ l b, u b)} \end{array}

Here we only need to update $α$ , $β$ , $γ$ , $ξ$ , $ϕ^{(b_{i})}$ and $σ^{(b_{i})}$ in the inference process of inner HMM for ith base block. Our inference algorithm is a modification of the standard Forward-backward algorithm based on the defined $Q$ -function, as described in Appendix 2–algorithm 4.

Appendix 2–algorithm 4. Inference algorithm of inner HMM.
Input: $y^{(b i)}, ϕ^{(b i)}, σ^{(b i)}, ϕ^{(b i - 1)}, σ^{(b i - 1)}, ϕ^{(b i + 1)}, σ^{(b i + 1)}, T^{i n n e r}, π^{i n n e r}, l b, u b$ Output: ${\hat{ω}}^{(b i)}, {\hat{h}}^{(b i)}, {\hat{ϕ}}^{(b i)}, {\hat{σ}}^{(b i)}$ Initialize ${\hat{h}}^{(b i)}$ and ${\hat{ω}}^{(b i)}$ to the $1 \times N_{b_{i}}$ all zero vector. Initialize α, β, γ, ξ // EM algorithm for round = 1, 2, · · · , R do E-step: update α, β, γ, and ξ given $y^{(b i)}, ϕ^{(b i)}, σ^{(b i)}, ϕ^{(b i - 1)}, σ^{(b i - 1)}, ϕ^{(b i + 1)}, σ^{(b i + 1)}, T^{i n n e r}, π^{i n n e r}, l b, u b$ M-step: calculate ${\hat{ϕ}}^{(b i)}$ and ${\hat{σ}}^{(b i)}$ by Equation 2 and Equation 3 $ϕ^{(b i)} \leftarrow {\hat{ϕ}}^{(b i)}$ $σ^{(b i)} \leftarrow {\hat{σ}}^{(b i)}$ Calculate the Q-function (log marginal likelihood) $\hat{q}$ given all estimated parameters end Derive ${\hat{h}}^{(b_{i})}$ from $γ$ by taking the state with highest probability for each data point, i.e. $h_{n}^{(b_{i})} = a r g m a x ([γ_{n 1}, γ_{n 2}, γ_{n 3}, γ_{n 4}])$ Set each element of ${\hat{ω}}^{(b_{i})}$ to $e x p (\hat{q} / N_{b_{i}})$ .

Appendix 2–algorithm 4. Inference algorithm of inner HMM.

Input:

y^{(b i)}, ϕ^{(b i)}, σ^{(b i)}, ϕ^{(b i - 1)}, σ^{(b i - 1)}, ϕ^{(b i + 1)}, σ^{(b i + 1)}, T^{i n n e r}, π^{i n n e r}, l b, u b

Output:

{\hat{ω}}^{(b i)}, {\hat{h}}^{(b i)}, {\hat{ϕ}}^{(b i)}, {\hat{σ}}^{(b i)}

Initialize

{\hat{h}}^{(b i)}

and

{\hat{ω}}^{(b i)}

to the

1 \times N_{b_{i}}

all zero vector.
Initialize α, β, γ, ξ
// EM algorithm
for round = 1, 2, · · · , R do
E-step: update α, β, γ, and ξ given

y^{(b i)}, ϕ^{(b i)}, σ^{(b i)}, ϕ^{(b i - 1)}, σ^{(b i - 1)}, ϕ^{(b i + 1)}, σ^{(b i + 1)}, T^{i n n e r}, π^{i n n e r}, l b, u b

M-step: calculate

{\hat{ϕ}}^{(b i)}

and

{\hat{σ}}^{(b i)}

by Equation 2 and Equation 3

ϕ^{(b i)} \leftarrow {\hat{ϕ}}^{(b i)}

σ^{(b i)} \leftarrow {\hat{σ}}^{(b i)}

Calculate the Q-function (log marginal likelihood)

\hat{q}

given all estimated parameters
end
Derive

{\hat{h}}^{(b_{i})}

from

γ

by taking the state with highest probability for each data point, i.e.

h_{n}^{(b_{i})} = a r g m a x ([γ_{n 1}, γ_{n 2}, γ_{n 3}, γ_{n 4}])

Set each element of

{\hat{ω}}^{(b_{i})}

e x p (\hat{q} / N_{b_{i}})

ϕ^{(b_{i})} = \frac{\sum_{n = 1}^{N_{b_{i}}} γ_{n 1} y_{n}^{(b_{i})}}{\sum_{n = 1}^{N_{b_{i}}} γ_{n 1}}

σ^{(b_{i})} = \frac{\sum_{n = 1}^{N_{b_{i}}} γ_{n 1} (y_{n}^{(b_{i})} - ϕ^{(b_{i})})^{2}}{\sum_{n = 1}^{N_{b_{i}}} γ_{n 1}}

GPU-accelerated parameter inference

Given the general inference algorithm (‘General inference algorithm’) and algorithms for calculating the emission probabilities for the inner HMM (Appendix 2–algorithm 4) and linear model (Appendix 2–algorithm 3), we have the full algorithm for the parameter inference. The actual inference, however, can be really slow since the input fast5 data are huge. Here we use heuristics and GPU parallelization to accelerate the parameter inference.

We divide the raw current signal of each read $y$ into overlapping regions, which start at an atomic base block and also end with an atomic base block. As shown in Appendix 2—figure 2, we process the rth region $y^{(r)}$ independently to get the segmentation results $z^{(r)}$ . This process is based on the following heuristics: the segmentation result of one region will not affect another region when they are far away from each other. The inference is performed in two rounds. In the first round, we obtain the segmentation results for non-overlapping regions. The center part of the segmentation is not affected by the flank parts, which may be affected by neighboring regions. In the second round, we merge relevant parts of neighboring regions near the touch point as a new region based on segmentation results of the first round. After that, we perform the inference for the new regions to get the final segmentations for the whole read.

We use GPU to parallelize the inference algorithm of inner HMM (Alg. 4). For each region, we first enumerate all possible paths and calculate the likelihood for each candidate path (‘Enumeration of hidden states of outer HMM’). As shown in Appendix 2—figure 2, we need to perform inner HMM parameter inference for different paths of different regions of one read. We consider parallelizing the computation of different merged blocks $s e g (i, j)$ in different paths of different regions. We first collect all merged blocks and sort them by size and only keep unique blocks, which can be identified by the combinatorial key of region index, start atomic block index, end atomic block index. Then we run the inner HMM inference (Alg. 4) of each merged block on a GPU core. Note that all GPU cores perform the same computation step of Alg. 4 simultaneously, but for different merged blocks.

Appendix 2—figure 1

Download asset Open asset

Illustration of enumeration of hidden states of outer HMM.

(a) Raw signal and initialized to atomic segmentations. (b) Examples of the full configuration set $G = (g_{1}, g_{2}, \dots, g_{i}, \dots, g_{M})$ .

Appendix 2—figure 2

Download asset Open asset

Illustration of GPU-accelerated parameter inference.

Data availability

The data utilized in this study are obtained from publicly available repositories. Details regarding the accession number and data processing can be found in Methods. The resulting data and the source code are hosted on GitHub (https://github.com/guangzhaocs/SegPore copy archived at Cheng, 2025).

The following previously published data sets were used

1. Pratanwanich PN
(2021) NCBI BioProject
ID PRJEB40872. Differential RNA modifications from dRNA-seq.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB40872
(2021) NCBI BioProject
ID PRJEB44348. SGNEx: A systematic benchmark of Nanopore long read RNA sequencing for transcript level analysis in human cell lines.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB44348
1. Jenjaroenpun P
2. Wongsurawat T
3. Wadley TD
4. Wassenaar TM
5. Liu J
6. Dai Q
7. Wanchai V
8. Akel NS
9. Jamshidi-Parsian A
10. Franco AT
11. Boysen G
12. Jennings ML
13. Ussery DW
14. He C
15. Nookaew I
(2020) NCBI Sequence Read Archive
ID SRP166020. Decoding Epitranscriptional Landscapes form Native RNA Sequences.

https://trace.ncbi.nlm.nih.gov/Traces/?view=study&acc=SRP166020

References

1. Acera Mateos P
2. J Sethi A
3. Ravindran A
4. Srivastava A
5. Woodward K
6. Mahmud S
7. Kanchi M
8. Guarnacci M
9. Xu J
10. W S Yuen Z
11. Zhou Y
12. Sneddon A
13. Hamilton W
14. Gao J
15. M Starrs L
16. Hayashi R
17. Wickramasinghe V
18. Zarnack K
19. Preiss T
20. Burgio G
21. Dehorter N
22. E Shirokikh N
23. Eyras E
(2024) Prediction of m6A and m5C at single-molecule resolution reveals a transcriptome-wide co-occurrence of RNA modifications
Nature Communications 15:3899.

https://doi.org/10.1038/s41467-024-47953-7
- PubMed
- Google Scholar
(2017) The importance of being modified: The role of RNA modifications in translational fidelity
The Enzymes 41:1–50.

https://doi.org/10.1016/bs.enz.2017.03.005
- PubMed
- Google Scholar
1. Bellodi C
2. McMahon M
3. Contreras A
4. Juliano D
5. Kopmar N
6. Nakamura T
7. Maltby D
8. Burlingame A
9. Savage SA
10. Shimamura A
11. Ruggero D
(2013) H/ACA small RNA dysfunctions in disease reveal key roles for noncoding RNA modifications in hematopoietic stem cell differentiation
Cell Reports 3:1493–1502.

https://doi.org/10.1016/j.celrep.2013.04.030
- PubMed
- Google Scholar
Book
1. Bishop CM
2. Nasrabadi NM
(2006) Pattern Recognition and Machine Learning
Springer.

https://doi.org/10.1117/1.2819119
- Google Scholar
Software
1. blawrence-ont
2. malton-ont
(2026) Dorado, version ba44a01
GitHub.

https://github.com/nanoporetech/dorado
1. Boccaletto P
2. Stefaniak F
3. Ray A
4. Cappannini A
5. Mukherjee S
6. Purta E
7. Kurkowska M
8. Shirvanizadeh N
9. Destefanis E
10. Groza P
11. Avşar G
12. Romitelli A
13. Pir P
14. Dassi E
15. Conticello SG
16. Aguilo F
17. Bujnicki JM
(2022) MODOMICS: a database of RNA modification pathways. 2021 update
Nucleic Acids Research 50:D231–D235.

https://doi.org/10.1093/nar/gkab1083
- PubMed
- Google Scholar
Software
1. Brennen F
(2023) Kmer_models, version 4e56dae
GitHub.

https://github.com/nanoporetech/kmer_models
1. Caldwell CC
2. Spies M
(2017) Helicase SPRNTing through the nanopore
PNAS 114:11809–11811.

https://doi.org/10.1073/pnas.1716866114
- PubMed
- Google Scholar
1. Chen AY
2. Owens MC
3. Liu KF
(2023) Coordination of RNA modifications in the brain and beyond
Molecular Psychiatry 28:2737–2749.

https://doi.org/10.1038/s41380-023-02083-2
- PubMed
- Google Scholar
1. Chen Y
2. Davidson NM
3. Wan YK
4. Yao F
5. Su Y
6. Gamaarachchi H
7. Sim A
8. Patel H
9. Low HM
10. Hendra C
11. Wratten L
12. Hakkaart C
13. Sawyer C
14. Iakovleva V
15. Lee PL
16. Xin L
17. Ng HEV
18. Loo JM
19. Ong X
20. Ng HQA
21. Wang J
22. Koh WQC
23. Poon SYP
24. Stanojevic D
25. Tran H-D
26. Lim KHE
27. Toh SY
28. Ewels PA
29. Ng H-H
30. Iyer NG
31. Thiery A
32. Chng WJ
33. Chen L
34. DasGupta R
35. Sikic M
36. Chan Y-S
37. Tan BOP
38. Wan Y
39. Tam WL
40. Yu Q
41. Khor CC
42. Wüstefeld T
43. Lezhava A
44. Pratanwanich PN
45. Love MI
46. Goh WSS
47. Ng SB
48. Oshlack A
49. Göke J
50. SG-NEx consortium
(2025) A systematic benchmark of Nanopore long-read RNA sequencing for transcript-level analysis in human cell lines
Nature Methods 22:801–812.

https://doi.org/10.1038/s41592-025-02623-4
- PubMed
- Google Scholar
Software
1. Cheng G
(2025) SegPore, version swh:1:rev:336603875525555b46e2c577dcce1c0ac29eb9f2
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:9c796fb3d8a3b723800d3118855f23ca318af578;origin=https://github.com/guangzhaocs/SegPore;visit=swh:1:snp:6bc0911addbe516c2cb8ee192d2cd5fb38eb6326;anchor=swh:1:rev:336603875525555b46e2c577dcce1c0ac29eb9f2
1. Craig JM
2. Laszlo AH
3. Brinkerhoff H
4. Derrington IM
5. Noakes MT
6. Nova IC
7. Tickman BI
8. Doering K
9. de Leeuw NF
10. Gundlach JH
(2017) Revealing dynamics of helicase translocation on single-stranded DNA using high-resolution nanopore tweezers
PNAS 114:11932–11937.

https://doi.org/10.1073/pnas.1711282114
- PubMed
- Google Scholar
(2025) De novo basecalling of RNA modifications at single molecule and nucleotide resolution
Genome Biology 26:38.

https://doi.org/10.1186/s13059-025-03498-6
- PubMed
- Google Scholar
(2020) GPU accelerated adaptive banded event alignment for rapid comparative nanopore signal analysis
BMC Bioinformatics 21:343.

https://doi.org/10.1186/s12859-020-03697-x
- PubMed
- Google Scholar
1. Gao Y
2. Liu X
3. Wu B
4. Wang H
5. Xi F
6. Kohnen MV
7. Reddy ASN
8. Gu L
(2021) Quantitative profiling of N⁶-methyladenosine at single-base resolution in stem-differentiating xylem of Populus trichocarpa using Nanopore direct RNA sequencing
Genome Biology 22:22.

https://doi.org/10.1186/s13059-020-02241-7
- PubMed
- Google Scholar
1. Hendra C
2. Pratanwanich PN
3. Wan YK
4. Goh WSS
5. Thiery A
6. Göke J
(2022) Detection of m6A from direct RNA sequencing using a multiple instance learning framework
Nature Methods 19:1590–1598.

https://doi.org/10.1038/s41592-022-01666-1
- PubMed
- Google Scholar
Software
1. Hendra C
(2025) M6anet, version 590ec27
GitHub.

https://github.com/GoekeLab/m6anet
1. Jenjaroenpun P
2. Wongsurawat T
3. Wadley TD
4. Wassenaar TM
5. Liu J
6. Dai Q
7. Wanchai V
8. Akel NS
9. Jamshidi-Parsian A
10. Franco AT
11. Boysen G
12. Jennings ML
13. Ussery DW
14. He C
15. Nookaew I
(2021) Decoding the epitranscriptional landscape from native RNA sequences
Nucleic Acids Research 49:e7.

https://doi.org/10.1093/nar/gkaa620
- PubMed
- Google Scholar
1. Koh CWQ
2. Goh YT
3. Goh WSS
(2019) Atlas of quantitative single-base-resolution N⁶-methyl-adenine methylomes
Nature Communications 10:5636.

https://doi.org/10.1038/s41467-019-13561-z
- PubMed
- Google Scholar
1. Körtel N
2. Rücklé C
3. Zhou Y
4. Busch A
5. Hoch-Kraft P
6. Sutandy FXR
7. Haase J
8. Pradhan M
9. Musheev M
10. Ostareck D
11. Ostareck-Lederer A
12. Dieterich C
13. Hüttelmaier S
14. Niehrs C
15. Rausch O
16. Dominissini D
17. König J
18. Zarnack K
(2021) Deep and accurate detection of m6A RNA modifications using miCLIP2 and m6Aboost machine learning
Nucleic Acids Research 49:e92.

https://doi.org/10.1093/nar/gkab485
- PubMed
- Google Scholar
1. Kovaka S
2. Hook PW
3. Jenike KM
4. Shivakumar V
5. Morina LB
6. Razaghi R
7. Timp W
8. Schatz MC
(2025) Uncalled4 improves nanopore DNA and RNA modification detection via fast and accurate signal alignment
Nature Methods 22:681–691.

https://doi.org/10.1038/s41592-025-02631-4
- PubMed
- Google Scholar
1. Lee H
2. Bao S
3. Qian Y
4. Geula S
5. Leslie J
6. Zhang C
7. Hanna JH
8. Ding L
(2019) Stage-specific requirement for Mettl3-dependent m⁶A mRNA methylation during haematopoietic stem cell differentiation
Nature Cell Biology 21:700–709.

https://doi.org/10.1038/s41556-019-0318-1
- PubMed
- Google Scholar
1. Li H
(2018) Minimap2: pairwise alignment for nucleotide sequences
Bioinformatics 34:3094–3100.

https://doi.org/10.1093/bioinformatics/bty191
- PubMed
- Google Scholar
Software
1. Li H
(2021) Minimap2, version fe35e67
GitHub.

https://github.com/lh3/minimap2/releases/tag/v2.24
(2015) Single-nucleotide-resolution mapping of m6A and m6Am throughout the transcriptome
Nature Methods 12:767–772.

https://doi.org/10.1038/nmeth.3453
- PubMed
- Google Scholar
1. Liu H
2. Begik O
3. Lucas MC
4. Ramirez JM
5. Mason CE
6. Wiener D
7. Schwartz S
8. Mattick JS
9. Smith MA
10. Novoa EM
(2019) Accurate detection of m6A RNA modifications in native RNA sequences
Nature Communications 10:4079.

https://doi.org/10.1038/s41467-019-11713-9
- Google Scholar
(2015) A complete bacterial genome assembled de novo using only nanopore sequencing data
Nature Methods 12:733–735.

https://doi.org/10.1038/nmeth.3444
- PubMed
- Google Scholar
1. Lorenz DA
2. Sathe S
3. Einstein JM
4. Yeo GW
(2020) Direct RNA sequencing enables m⁶A detection in endogenous transcript isoforms at base-specific resolution
RNA 26:19–28.

https://doi.org/10.1261/rna.072785.119
- PubMed
- Google Scholar
1. Meyer KD
2. Saletore Y
3. Zumbo P
4. Elemento O
5. Mason CE
6. Jaffrey SR
(2012) Comprehensive analysis of mRNA methylation reveals enrichment in 3’ UTRs and near stop codons
Cell 149:1635–1646.

https://doi.org/10.1016/j.cell.2012.05.003
- PubMed
- Google Scholar
1. Nance KD
2. Meier JL
(2021) Modifications in an emergency: The role of N1-Methylpseudouridine in COVID-19 Vaccines
ACS Central Science 7:748–756.

https://doi.org/10.1021/acscentsci.1c00197
- PubMed
- Google Scholar
1. Ohira T
2. Suzuki T
(2024) Transfer RNA modifications and cellular thermotolerance
Molecular Cell 84:94–106.

https://doi.org/10.1016/j.molcel.2023.11.041
- PubMed
- Google Scholar
1. Pratanwanich PN
2. Yao F
3. Chen Y
4. Koh CWQ
5. Wan YK
6. Hendra C
7. Poon P
8. Goh YT
9. Yap PML
10. Chooi JY
11. Chng WJ
12. Ng SB
13. Thiery A
14. Goh WSS
15. Göke J
(2021) Identification of differential RNA modifications from nanopore direct RNA sequencing with xPore
Nature Biotechnology 39:1394–1402.

https://doi.org/10.1038/s41587-021-00949-w
- PubMed
- Google Scholar
(2020) Post-translational modifications of the fragile X mental retardation protein in neuronal function and dysfunction
Molecular Psychiatry 25:1688–1703.

https://doi.org/10.1038/s41380-019-0629-4
- PubMed
- Google Scholar
1. Quin J
2. Sedmík J
3. Vukić D
4. Khan A
5. Keegan LP
6. O’Connell MA
(2021) ADAR RNA modifications, the epitranscriptome and innate immunity
Trends in Biochemical Sciences 46:758–771.

https://doi.org/10.1016/j.tibs.2021.02.002
- PubMed
- Google Scholar
1. Simpson JT
2. Workman RE
3. Zuzarte PC
4. David M
5. Dursi LJ
6. Timp W
(2017) Detecting DNA cytosine methylation using nanopore sequencing
Nature Methods 14:407–410.

https://doi.org/10.1038/nmeth.4184
- PubMed
- Google Scholar
Preprint
1. Stoiber M
2. Quick J
3. Egan R
4. Eun Lee J
5. Celniker S
6. Neely RK
7. Loman N
8. Pennacchio LA
9. Brown J
(2017) De Novo identification of DNA modifications enabled by genome-guided nanopore signal processing
bioRxiv.

https://doi.org/10.1101/094672
- Google Scholar
(2025) Duplexed direct RNA sequencing protocol using polyadenylation and polyuridylation
Microbiology Resource Announcements 14:e0104124.

https://doi.org/10.1128/mra.01041-24
- PubMed
- Google Scholar
(2021) The N6-methyladenosine RNA modification in acute myeloid leukemia
Current Opinion in Hematology 28:80–85.

https://doi.org/10.1097/MOH.0000000000000636
- PubMed
- Google Scholar
1. Zhong Z-D
2. Xie Y-Y
3. Chen H-X
4. Lan Y-L
5. Liu X-H
6. Ji J-Y
7. Wu F
8. Jin L
9. Chen J
10. Mak DW
11. Zhang Z
12. Luo G-Z
(2023) Systematic comparison of tools used for m⁶A mapping from nanopore direct RNA sequencing
Nature Communications 14:1906.

https://doi.org/10.1038/s41467-023-37596-5
- PubMed
- Google Scholar
(2023) The expanding role of RNA modifications in plant RNA polymerase II transcripts: highlights and perspectives
Journal of Experimental Botany 74:3975–3986.

https://doi.org/10.1093/jxb/erad136
- PubMed
- Google Scholar

Article and author information

Author details

Guangzhao Cheng

Department of Computer Science, Aalto University, Espoo, Finland

Contribution
Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0008-6160-3637
Aki Vehtari

Department of Computer Science, Aalto University, Espoo, Finland

Contribution
Resources, Supervision, Writing – original draft

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2164-9469
Lu Cheng
1. Department of Computer Science, Aalto University, Espoo, Finland
2. Institute of Biomedicine, University of Eastern Finland, Kuopio, Finland
Contribution
Conceptualization, Resources, Software, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing

For correspondence
lu.cheng.ac@gmail.com

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-6391-2360

Funding

Research Council of Finland (335858)

Guangzhao Cheng
Lu Cheng

Research Council of Finland (358086)

Guangzhao Cheng
Lu Cheng

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Prof. Zhijie Tan from Wuhan University for a useful discussion about the molecule dynamics of Nanopore sequencing, Dr. Dan Zhang from Sichuan University for helpful tutorials about Nanopore analysis workflows. We also thank Prof. Luo Guanzheng for sharing the m6A benchmark baseline results. GC and LC acknowledge the computational resources provided by the Aalto Science-IT project. This work was supported by Research Council of Finland grants [335858, 358086 to GC and LC].

Version history

Preprint posted: October 10, 2024
Sent for peer review: December 4, 2024
Reviewed Preprint version 1: March 6, 2025
Reviewed Preprint version 2: October 15, 2025
Reviewed Preprint version 3: November 10, 2025
Version of Record published: March 2, 2026

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.104618. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.