Precise Nanopore Signal Modeling Improves Unsupervised Single-Molecule Methylation Detection

Abstract

Base calling in nanopore sequencing is a dificult and computationally intensive problem, typically resulting in high error rates. In many applications of nanopore sequencing, direct analysis of raw signal is a viable alternative. Dynamic time warping (DTW) is an important building block for raw signal analysis. In this paper, we propose several improvements to DTW class of algorithms to better account for specifics of nanopore signal modeling. We have implemented these improvements in a new signal-to-reference alignment tool Nadavca. We demonstrate that Nadavca alignments improve unsupervised methylation detection over Tombo. We also demonstrate that by providing additional information about the discriminative power of positions in the signal, an otherwise unsupervised method can approach the accuracy of supervised models. Availability and implementation: Nadavca is available under MIT license at https://github.com/fmfi-compbio/nadavca. Nanopore sequencing data sets are available from ENA bioproject PRJEB64246. Jaminaea angkorensis reference genome assembly is available from Zenodo https://doi.org/10.5281/zenodo.8145315.

1. Background

DNA sequencing devices developed by Oxford Nanopore

Technologies (ONT), including the pocket-sized MinION, measure changes in electrical current as a DNA strand passes through a nanopore. The result is a sequence of signal measurements (also called a squiggle). In this paper, we provide an improved algorithm to align parts of the squiggle to a known reference DNA sequence and show that improved alignment accuracy benefits downstream analysis in de novo methylation identification.

Typically, the first step in the analysis of nanopore sequencing data is base calling, which translates a squiggle into a DNA sequence. The most successful base callers are based on complex machine learning models, such as recurrent neural networks [ 1, 2 ]. A typical median read accuracy achieved by a base caller for R9.4.1 nanopores employed in this study is approx. 97% [ 3 ], but when base called reads are piled up at a deep coverage, the consensus accuracy can be as high as 99.94% [ 4 ]. (ONT lists 98.3% accuracy for R9.4.1 and 99% for newer R10.4.1 nanopores; note that these numbers are based on accuracy in the mode of the read accuracy distribution and are generally higher than median accuracy.)

An alternative class of approaches avoids base calling and works directly with the raw signal sequence. The main motivation is to use the additional information present in the raw signal, which is potentially lost in base calling. Applications include mainly better modification calling and methylation detection. Base calling is also computationally intensive, and thus in some situations, working with the raw signal can be faster.

Basic properties of nanopore signal. The measured

electrical current depends mostly on the context of consecutive nucleotides passing through the nanopore.

A signal model predicts for each context of size the expected signal level. For the case of nanopores in R9.4.1 lfow-cells, which were used in this paper, a typical model uses the context of size = 6. The DNA moves through the pore at the speed of roughly 450 bases per second and the signal is read at a frequency of 4000 samples per second. This means that on average we have approx. 9 samples per one DNA context, but the actual number of samples per context can vary significantly.

Signal-to-reference alignment and dynamic time warping. One of the basic building blocks of nanopore signal analysis is the alignment of squiggles to a known reference sequence. The squiggle can be represented as a sequence of numbers 1 · · · . For each nucleotide of the reference sequence 12 · · · , we need to assign the corresponding interval +1 · · · of the squiggle which we call an event. Such squiggle-to-reference alignment facilitates further analysis, allowing to find significant diferences between the squiggle and the reference signal computed based on the signal model applied to the reference sequence. Identified diferences typically correspond to sequence variants or modifications [ 5 ]. Alignments also enable pooling together information from multiple reads aligned to a single reference location © 2023 Copyright for this paper by its authors. Use permitted under Creative Commons License and to visualize squiggles in the genomic context [ 6 ]. CPWErooUrckReshdoinpgs ISNc1e6u1r-3w-0s.o7r3g ACttEribUutRion W4.0oInrtekrnsahtioonpal (PCCroBYce4.0e).dings (CEUR-WS.org) The squiggle alignment is typically performed by a class of algorithms called dynamic time warping (DTW), Supervised and unsupervised modification detecoriginating in the field of speech processing [ 7 ]. The tion. It is known that modified nucleotides afect the alignment of a signal level to a position in the ref- signal levels when passing through nanopore [ 10 ]. For erence is evaluated by a scoring function (, , ), known modification patterns, such as CpG methylation, which typically takes into account the context around it is possible to prepare artificially methylated sequences position ; for example, the scoring function (, , ) = and either estimate a signal model for contexts containing ( − expectedsignal(− (− 1)/2 · · · +(− 1)/2))2 uses a methylated nucleotides [ 11, 12, 13 ] or to train a machine signal model based on the context of nucleotides (as- learning algorithm, such as a neural network, to recsuming is odd). Using a scoring function of this type, ognize signals from sequences that contain methylated we are looking for the alignment minimizing the score nucleotides [ 14, 15, 16 ]. We call this approach supera∑︀dy=n1a∑m︀ic=prog(ram, m, in)g, walhgiocrhitchamn sbiemciolamr ptoutseedquuesnincge lvaibseedlemd ettrhayinlaintigondadteateacrteiorne,qsuinirceedstuobsbtuainldtiasluacmhomuondteslos.f alignment. In contrast, unsupervised modification detection can be

In our experience, the DTW-like algorithms applied in used when labeled training data is not available, or when the context of MinION squiggle alignment show a variety we are targeting previously uncharacterized DNA modiof artefacts. The scoring function allows squiggle to be ifcations. One option is to compare signal values in reads compressed or stretched arbitrarily in time, without any in the studied sample with signals in a control sample penalty. If several consecutive sequence contexts along devoid of modifications and look for statistically signifthe reference genome exhibit similar expected signal val- icant diferences [ 6, 17 ]. However, it is also possible to ues in the signal model (which is not uncommon), there ifnd putative modifications simply by comparing the acis a large uncertainty in how the corresponding observed tual signal with the expected signal predicted by a signal signal should be split between the events. Depending model trained on unmodified sequences (this approach is on the exact DTW implementation, this may lead to arti- implemented for example in Tombo [ 18 ]). The strength ifcially short or artificially long events. Small amounts of this method is that it requires only unmodified training of noise in the signal can easily lead to events split in- data to obtain the signal model and then can be applied correctly between multiple reference contexts. In fact, directly to biological samples from various organisms the real squiggles often contain transitional readouts be- without the need for modified training data or matched tween signal levels corresponding to consecutive signal unmodified dataset for each analysis. The crucial step in contexts, sometimes even alternating between two con- this approach is accurate alignment of signal to the refsecutive contexts for some time. If the squiggle does not erence to allow meaningful comparison with the model correspond exactly to the reference sequence (e.g. in case signal. Using Nadavca signal-to-reference alignments, of single nucleotide variants or base modifications), DTW we implemented a new tool nanometh for unsupervised tends to “hide” the diferences between the reference and modification detection. We demonstrate that Nadavca the squiggle within the neighboring events. This subse- alignments increase the accuracy of unsupervised modiquently afects the ability to detect such diferences. ifcation detection, and we also evaluate contribution of

In this paper, we propose several improvements to the individual features implemented in Nadavca. DTW framework to overcome these obstacles. Instead of simple maximization of the scoring function, we pro- 2. Methods pose to employ alternative methods similar to posterior decoding in probabilistic models [ 8, 9 ], to account for 2.1. Realignment of Signal to the uncertainty in the event boundaries. We enhance the signal model underlying the scoring function to take into Reference with Posterior Decoding account a wider context of the sequence and also specif- The core part of Nadavca aligns a portion of nanopore ically model transitions between the contexts. We also signal (values 1, · · · , ) to the corresponding part of employ a two-iteration approach to signal normalization the reference genome (bases 1, · · · , ). The objective to avoid artefacts resulting from incorrect shifting and is to improve the accuracy of an approximate alignment, scaling of the squiggle. We have implemented these im- resulting from aligning base-called reads to the reference. provements in an alignment tool called Nadavca, and we The output alignment assigns an interval of the signal have evaluated their efect on the accuracy of detection sequence , · · · , (also called event) to each base of base modifications in nanopore reads. To make the of the reference. For the basic version of the algorunning time practical, we use methods similar to banded rithm, we require that the intervals for adjacent bases and anchored sequence alignment. are adjacent in the signal sequence (i.e. +1 = + 1).

The standard dynamic time warping (DTW) algorithm is based on dynamic programming. We define subproblem , , which represents the score of the best alignment of the squiggle segment 1 · · · to the reference segment 1 · · · , assuming that the signal is aligned to the reference base . , can be computed gradually using a simple recurrence , = (, , ) + min(− 1,′ , − 1,− 1), where (, , ) is the score of aligning signal value to the reference sequence at position . After finishing the computation, the fi- Figure 1: Extended signal model nal alignment (intervals [, ]) can be easily computed using standard methods for tracing back the dynamic programming solution. The alignment can thus be visualized supplied by Oxford Nanopore, which predicts the exas a path in this dynamic programming matrix. pected signal based on the sequence context of length

Our alignment algorithm is inspired by the posterior 6. Instead, we use a (restricted) 10-mer model which we decoding approach to sequence alignment [ 8, 9 ]. While compute in the following way. the standard DTW seeks the alignment with the highest First, we align the training dataset (described below) to overall score, our goal is to consider also a contribution the reference genome using Nadavca with the default 6of sub-optimal alignments. Many of these alignments mer model. Computing the mean signal for each 10-mer can have scores very close to the optimum, representing would lead to a huge model, requiring a gigantic amount uncertainty in the true alignment. Posterior decoding of data in order to cover each 10-mer with enough samalgorithms consider this uncertainty at each position of ples. Therefore, we restrict our 10-mer model by expressthe alignment. ing it as the sum of three separate sub-models: a base 6

To apply this principle to DTW, we first modify the mer model for the central part of the context and two scoring scheme and instead of the distance score (, , ), 4-mer models 1 and 2 that fine-tune the result based we use probabilistic score (, , ), which gives a prob- on more distant bases. For a context = 12 . . . 10, the ability of signal value being produced from the ref- expected level () of the signal is computed as follows: erence at position . Details of this probabilistic scor- () = 1(1 . . . 4) + (3 . . . 8) + 2(7 . . . 10) ing scheme are described below. Now the score of the (see Figure 1 for illustration). Models 1, 2, and alignment becomes the product of probabilities along the are estimated simultaneously for all 4-mers and 6-mers realignment path ∏︀=1 ∏︀= (, , ). Higher scores spectively, by fitting the observed signal levels using the indicate higher alignment quality. least-squares method. The probabilistic score (, , )

The key concept in posterior decoding is the posterior is now computed based on ∼ ( (− 4 . . . +5), ) probability , which for each and sums the scores of for constant = 0.3. all alignments that align signal value to the reference base . All values of , can be computed by Forward- 2.3. Transitional States Backward algorithm similarly as in the context of pair hidden Markov models [8, section 4.4]. Finally, we com- The simplest models [ 19 ] of the nanopore signal assume pute the resulting alignment as a path through the matrix that the signal can be partitioned into events, each event maximizing the product of individual posterior scores corresponding to one position in the reference and thus , (as in [ 9 ]), using a dynamic programming recurrence having a roughly uniform signal level influenced by the similar to the standard DTW. local sequence context. However, transitions between

In the real data applications, we impose additional re- events are not instantaneous, and as a result, the meastrictions. The minimum length of each event − + 1 sured signal may include intermediate values between is constrained by the parameter min_event_length. The the expected levels typical for two adjacent events or it very start and the very end of the signal may remain un- may even jump back and forth between these levels (see aligned. The length of unaligned portions is at most Figure 2). 2bandwidth + 1, where bandwidth is a parameter de- In Nadavca, we handle these phenomena by inserting a scribed below. Corresponding modifications of the al- special "transitional state” between each pair of adjacent gorithms are straightforward. reference bases and +1. The primary score for aligning signal to this transitional state is computed from a 2.2. Enhanced Signal Model special distribution as opposed to the -mer distribution used normally. While this distribution can depend on Accurate models of signal levels expected in diferent the sequence context, for simplicity we currently use a sequence contexts are important for precise sequence constant score 0.01. Thus, any signal aligned to these alignment. Tombo [ 18 ] by default uses the 6-mer model positions is penalized, regardless of the signal level and the sequence context. This simple heuristic works well in practice.

The posterior decoding algorithm aligns transitional states to intervals of signal values, but these intervals are not constrained to have a minimum required length; they can even be empty. The intervals aligned to transitional states are omitted from Nadavca’s output, and thus the output events corresponding to the reference bases may have gaps between them. 2.4. Fast Implementation With Banded

Alignment

In the next step, we consider all matched bases from the BWA-MEM alignment. For each base called base matched to a reference base , we locate the interval ... in the signal that corresponds to according to the Events table provided by Guppy. Interval ... will be anchored to the reference base and we call the set of all anchors the guide alignment.

We use the guide alignment to determine the portion of the signal that needs to be aligned to the reference. It will extend from the start of the interval in the first anchor to the the end of the interval in tha last anchor, extended further on both ends by the bandwidth parameter.

The guide alignment will also define an area of interest for dynamic programming. Recall that our posterior alignment uses three passes of dynamic programming (forward, backward, and posterior), and all three of them will consider only such alignment paths that are contained within the area of interest. All values outside of the area of interest are considered as zeroes. For each reference base , the area of interest contains the interval of signal positions between the start of the closest anchor at or before base and the end of the closest anchor at or after base ; this interval is extended by bandwidth positions on both ends (Figure 4). 2.5. Signal Renormalization

To avoid slow computation of the full dynamic program

ming matrix, we use the banded alignment heuristic Raw nanopore signal needs to be normalized by applying around an approximate guide alignment obtained by read-specific shift and scale parameters (denoted and aligning the base called sequence to the reference. respectively). Given a single raw signal level raw and

In particular, we first base call each read, and align the parameters and , the normalized signal level is the base called sequence to the target genomic sequence (raw − )/ . This normalization can influence subseusing BWA-MEM [ 20 ]. The BWA-MEM alignment will quent alignment [ 21 ] and other processing. In Nadavca, determine the window of the target genome to which we initially normalize the signal by using the median we need to align the signal. This window is padded with signal value in a read as and the median of |raw − | several additional bases from base call to introduce con- in a read as . This median-based normalization is, howtext for the signal model; please note that these bases ever, sensitive to various nanopore data artefacts. For may be diferent from the target due to the adapters or this reason, we improve the normalization by an iterative barcodes present at the beginning or at the end of the procedure [ 22, 21, 18 ]. read (see Figure 3). The resulting sequence will serve as In particular, after the initial normalization, the rethe reference for the signal alignment. alignment algorithms assigns signal values to particular reference locations. For each such aligned location, our mM EDTA pH 8.0, resuspended in 100 mM EDTA pH 8.0, enhanced signal model predicts a particular expected 2 % (v/v) 2-mercaptoethanol and incubated for 15 min at signal level. In the next iteration, we update normaliza- room temperature. Cells were then pelleted, resuspended tion parameters and to minimize the sum of squared in 0.15 M NaCl, 100 mM EDTA pH 8.0 and disrupted by errors between the expected and observed signal val- vortexing with glass beads (0.4-0.5 mm, Sigma) three ues in the read (omitting signal values aligned to tran- times for 30 s, followed by addition of sodium dodecyl sitional states). Subsequently, we repeat the process by sulfate (SDS) to 0.1 % (w/v). Nucleic acids were obtained re-computing the alignment and recomputing normaliza- by series of phenol and phenol : chloroform : isoamylaltion again. We observed, that two iterations are suficient cohol (25 : 24 : 1) extractions, precipitated with an equal to achieve a good accuracy. volume of 96 % (v/v) ethanol, washed with 70 % (v/v) ethanol, and air-dried. The precipitate was dissolved in 2.6. Unsupervised Detection of DNA TE bufer (10 mM Tris-HCl, 1 mM EDTA, pH 8.0) and RNA was removed by RNase A (150 µg/mL digestion) for Modification Through Anomalies 30 min at 37°C. DNA was extracted by phenol : chloroDNA modifications typically cause small shifts in sig- form : isoamylalcohol (25 : 24 : 1), precipitated with 0.1 nal values. Given an accurate signal-to-reference align- M NaCl and two volumes of 96 % (v/v) ethanol, washed ment, this property can be exploited to detect DNA mod- with 70 % (v/v) ethanol, air-dried, dissolved in TE bufer ifications by comparison of expected signal values pre- and further purified using a Genomic-tip 100/G (Qiagen). dicted from a signal model to the observed signal values HMW DNA of the yeast Magnusiomyces capitatus NRRL in individual events produced by the alignment. Here, Y-17686 (CBS 197.35; [ 25 ]) was prepared essentially as we mostly follow the approach used by Tombo de novo described above, except that prior the phenol extractions, modification detection module [ 18 ]. However, we sub- yeast cells were resuspended in 1 M Sorbitol, 10 mM stitute our improved signal-to-reference alignment and EDTA pH 8.0 and converted to spheroplasts using zyenhanced signal model to demonstrate their usefulness. molyase 20T (0.125 mg/mL; Seikagaku) treatment for

In particular, let [, ] be an event aligned to the refer- 60-90 min at 37°C. Spheroplasts were then lyzed in 0.15 ence context . Let is the event mean, i.e. the mean of M NaCl, 100 mM EDTA pH 8.0, 0.1 % (w/v) SDS. signal values . . . . As an indicator of DNA modifica- Nanopore sequencing was carried out in a MinION tion, we compute the p-value that value comes from a Mk1B device using a FLO-MIN106 (R9.4.1, revD) flow distribution ( (), ) defined by the enhanced signal cell (Oxford Nanopore Technologies). The sequencing model. libraries were constructed essentially as described in

This simple method, however, does not lead immedi- the manufacturer’s instructions. Datasets of control unately to accurate modification detection. On one hand, modified DNA and DNA methylated in vitro were prethe shifts in signal values are small and thus individual pared using the PCR barcoding kit (SQK-PBK004, Oxp-values assigned to events are rarely significant. On the ford Nanopore Technologies). Briefly, HMW DNA was other hand, a typical modification afects several events sheared in a g-Tube (Covaris) at 2700 ×g or 4300 ×g in near the modified base. Following the method outlined a MiniSpin Plus centrifuge (Eppendorf). Fragmented in Tombo, we combine p-values from three adjacent po- DNA (∼ 200 ng) was treated using a NEBNext Ultra II sitions using the Fisher’s combined probability test [ 23 ] End repair / dA-tailing module (New England Biolabs) to form a more reliable indicator. and ligated to the Barcode adapters (BCA) from the se

Finally, we output scores for positions where the ref- quencing kit (SQK-PBK004) using a Blunt/TA Ligase Maserence sequence matches a user-supplied pattern recog- ter Mix (New England Biolabs). DNA samples (∼ 10 ng) nized by the methyltransferase enzyme potentially active were then amplified using the Rapid Barcode primers in a given sample. For each position with the pattern, (LWB, SQK-PBK004) and LongAmp Taq DNA polymerase we consider a window of length 11 centered on the po- (New England Biolabs) and treated by exonuclease I (New tentially modified base and report the maximum score England Biolabs). DNA aliquots (∼ 1 µg) were modified within this window. using 160 µM S-adenosyl methionine and methyltransferases M.EcoRI (40 U), M.BamHI (12 U), or M.HhaI (25 2.7. Nanopore Sequencing U) (New England Biolabs) for 1 hour at 37°C. The methyltransferase reactions were then stopped at 65°C for 20 High-molecular weight (HMW) genomic DNA was pre- min. The modification by M. TaqI (10 U) was carried out pared from an overnight culture of Jaminaea angkorensis at 65°C for 1 hour. DNA was purified on AMPure XP C5b (CBS 10918; [ 24 ]) grown in yeast-peptone-dextrose beads and eluted into 10 mM Tris.Cl pH 8.0, 50 mM NaCl. (YPD) medium (1 % (w/v) yeast extract, 2 % (w/v) peptone, Methylated and control barcoded DNAs were pooled in 2 % (w/v) glucose) at 28°C with constant aeration. Yeast equal ratio. Rapid adapters (RAP, SQK-PBK004) were cells were harvested by centrifugation, washed with 20 attached to the ends of pooled DNAs (∼ 400 ng). The library preparation, flow cell priming, and loading were completed as described in the protocol (SQK-PBK004, version PBK_9073_v1_revA_23May2018). 2.8. Training and Testing Data Sets

All nanopore reads were base called by Albacore version

2.3.1. Base called reads were aligned to their respective reference genomes by BWA-MEM [ 26 ]. The reads that did not align to the reference for at least 80% of their length were discarded. Table 1 shows basic characteristics of resulting data sets.

Data set control0 was used for the purpose of training extended signal model. Mixtures of reads from control1 and met-* data sets were used for the purpose of testing the application of our methods to DNA methylation detection. Note that control0 and control1 data sets were produced from phylogenetically distant species and thus there is no overlap between training and testing sets.

3. Results and Discussion

3.1. Realignment by Nadavca Eliminates

Apparent Alignment Artefacts

Alignment of the signal to the reference sequence seg

ments the signal into events, each event corresponding To evaluate the efect of Nadavca alignments on downto a diferent context read by the pore. If the DNA moved stream analysis, we have implemented a simple tool through the pore at a constant speed, each event would nanometh for detection of methylation and other DNA span approximately 9 values. modifications from the nanopore signals. We have used

We have used Tombo [ 18 ], Nanopolish [ 11 ], and Na- the following general framework for detection of modifidavca to realign the signal to the reference sequence, cations (see Methods for details): using data set control1. Figure 5 shows the comparison of event length distributions. Only Nanopolish allows 1. Use Nadavca to align signal values of a read to events of length 0 (skip of the context), Tombo appar- the reference sequence ently has a minimum event length 3 (even though rarely 2. For each event implied by the alignment, compare it also outputs events of length 1 and 2), and in Nadavca the mean signal value within the event to the we require each event to be of length at least 2. With expected signal level from our enhanced signal Tombo, there is a clear bias against events of length 4, and model, producing a p-value preference for events of lengths 3 and 6; which appears 3. P-values from several adjacent positions are comto be an artefact of the alignment process. Interestingly, bined to a single score

Nanopolish does not report events of length 1 or 2. In case of Nadavca, the only apparent artefact is a higher abundance of events of length 2, which is a consequence of the requirement of minimum length 2 for each event.

3.2. Unsupervised Single-Molecule Single-Site DNA Modification Detection

We compare our results with de novo mode of Tombo.

Note that unlike most other methods, Nadavca and Tombo do not require any training, besides training the sequence context model. To make the results comparable, we use similar formulas for computing the p-values and scores as Tombo, however, substituting diferent signal alignment and diferent signal context model. For Tombo we have used parameters supplied directly with the tool.

We use a newly produced testing set from the genome of fungus Jaminaea angkorensis [ 24 ]. Genomic DNA was amplified by PCR, producing DNA devoid of any native modifications, and then one part (control) was left unmodified, and each of four additional parts were modified by a diferent methyltransferase enzyme. Each enzyme methylates either adenine or cytosine at a fixed position within a specific sequence pattern (for example, M.HhaI methylates the first cytosine in the pattern GCGC). In each test, we use a mixture of reads from the control sample and from one of the modified samples. We consider only positions matching the pattern of the corresponding enzyme (in the reference genome) and classify them as positives or negatives based on the sample of origin. Note that the signal model was trained on a data set from a diferent organism to avoid overfitting.

We do not attempt to set a threshold for calling a site methylated, instead we compute AUC score which evaluates each method over all choices of the threshold. In the computation of accuracy, each pattern occurrence within each read is considered separately. Table 2 shows an overview of results and comparison to Tombo de novo mode. On all of our methylation data sets, nanometh significantly outperforms Tombo. 3.4. Contribution of Individual Algorithm Features 3.3. Comparison to Supervised

Methylation Detection

Unlike Nadavca/nanometh, most tools for methylation detection are supervised, requiring a separate training set

We were also interested, how individual features of Na

davca afect the accuracy. Therefore, we started with a baseline modification detection tool and added individual features one at a time: posterior alignment, renormalization, modeling transitional signal states, and using an extended k-mer model. The results are shown in Table 3. Replacing Tombo resquiggle with posterior alignment with renormalization already helps to increase the ac

4. Conclusions

Acknowledgements. The strain of J. angkorensis was

In this work, we introduced Nadavca, a nanopore sig- kindly provided by Matthias Sipiczki (University of Denal aligner that incorporates several enhancements to brecen, Hungary). The strain of M. capitatus was kindly the DTW algorithm. Compared to existing aligners, Na- provided by Cletus P. Kurtzman and James Swezey (Agridavca’s output exhibits improved accuracy by eliminating cultural Research Service, Peoria, IL, USA) length distribution artifacts and eliminating the need for event segmentation as a preliminary step. We demon- References strated the eficacy of our alignment method by achieving state-of-the-art results in unsupervised methylation detection.

Most of the current tools for methylation detection in nanopore signal are supervised and based on deep learning (see e.g. Guppy, DeepSignal [ 14 ], DeepMod [ 28 ], Megalodon, or Remora). A disadvantage of these tools is that they require a substantial amount of labeled training data. Such data sets are available for the most common methylation types (e.g. CpG methylation) for which those tools were trained, but it is impossible to use these tools to discover less common or novel DNA modifications, where an unsupervised approach is needed. While unsupervised approaches are at a clear disadvantage when compared to supervised approaches in detection of common methylation types, by incorporating additional information such as afected signal positions due to methylation, our approach yields comparable results to Nanopolish (non-deeep-learning supervised approach) that has recently been shown to perform very well among the supervised methylation detection methods [ 27 ].

In the future, we aim to further refine our alignment technique by addressing alignment artifacts, particularly in regions where the expected signals for two events are similar, resulting in ambiguous alignments. Another task is to adapt our methods to the new generation R10 nanopores.

References

[1] H. Teng , M. D. Cao , M. B. Hall , T. Duarte , S. Wang , L. J. Coin , Chiron: translating nanopore raw signal directly into nucleotide sequence using deep learning , GigaScience 7 ( 2018 ) giy037 .

[2] V. Boža , B. Brejová , T. Vinař, DeepNano: deep recurrent neural networks for base calling in MinION nanopore reads , PloS One 12 ( 2017 ) e0178751 .

[3] V. Boža , P. Perešíni , B. Brejová , T. Vinař, Dynamic pooling improves nanopore base calling accuracy , IEEE/ACM Transactions on Computational Biology and Bioinformatics 19 ( 2021 ) 3416 - 3424 .

[4] R. R. Wick , L. M. Judd , K. E. Holt , Performance of neural network basecalling tools for oxford nanopore sequencing , Genome Biology 20 ( 2019 ) 1 - 10 .

[5] N. De Souza , Protein nanopores to detect dna methylation , Nature Methods 11 ( 2014 ) 8 - 8 .

[6] M. Stoiber , J. Quick , R. Egan , J. Eun Lee , S. Celniker , R. K. Neely , N. Loman , L. A. Pennacchio , J. Brown, De novo identification of DNA modifications enabled by genome-guided nanopore signal processing , bioRxiv ( 2017 ). URL: https://www.biorxiv.org/content/ early/2017/04/10/094672. doi: 10 .1101/094672.

[7] J. B. Kruskal , An overview of sequence comparison: Time warps, string edits, and macromolecules , SIAM Review 25 ( 1983 ) 201 - 237 .

[8] R. Durbin , S. R. Eddy , A. Krogh , G. Mitchison, Biological sequence analysis: probabilistic models of proteins and nucleic acids , Cambridge University Press, 1998 .

[9] G. Lunter , A. Rocco , N. Mimouni , A. Heger , A. Caldeira , J. Hein, Uncertainty in homology inferences: assessing and improving genomic sequence alignment , Genome research 18 ( 2008 ) 298 - 309 .

[10] L. Xu , M. Seki , Recent advances in the detection of base modifications using the Nanopore sequencer , Journal of Human Genetics 65 ( 2020 ) 25 - 33 .

[11] J. T. Simpson , R. E. Workman , P. Zuzarte , M. David , L. Dursi , W. Timp, Detecting DNA cytosine methylation using nanopore sequencing , Nature Methods 14 ( 2017 ) 407 .

[12] A. C. Rand , M. Jain , J. M. Eizenga , A. Musselman-Brown , H. E.

Funding. The research presented in this paper has been supported by funding from the Slovak Research and Development Agency grant APVV-18-0239 (JN), and by Olsen , M. Akeson , B. Paten , Mapping DNA methylation with high-throughput nanopore sequencing , Nature Methods 14 ( 2017 ) 411 - 413 .

[13] I. Lee , R. Razaghi , T. Gilpatrick , M. Molnar , A. Gershman , N. Sadowski , F. J. Sedlazeck , K. D. Hansen , J. T. Simpson , W. Timp, Simultaneous profiling of chromatin accessibility and methylation on human cell lines with nanopore sequencing , Nature Methods 17 ( 2020 ) 1191 - 1199 .

[14] P. Ni , N. Huang , Z. Zhang , D.-P. Wang , F. Liang , Y. Miao , C.-L. Xiao , F. Luo , J. Wang, DeepSignal: detecting DNA methylation state from Nanopore sequencing reads using deep-learning , Bioinformatics 35 ( 2019 ) 4586 - 4595 .

[15] A. B. McIntyre , N. Alexander , K. Grigorev , D. Bezdan , H. Sichtig , C. Y. Chiu , C. E. Mason , Single-molecule sequencing detection of N 6-methyladenine in microbial reference materials , Nature Communications 10 ( 2019 ) 579 .

[16] Q. Liu , L. Fang , G. Yu , D. Wang , C.-L. Xiao , K. Wang , Detection of DNA base modifications by deep recurrent neural network on Oxford Nanopore sequencing data , Nature Communications 10 ( 2019 ) 2449 .

[17] Q. Liu , D. C. Georgieva , D. Egli , K. Wang , NanoMod: a computational tool to detect DNA modifications using Nanopore long-read sequencing data , BMC Genomics 20 ( 2019 ) 31 - 42 .

[18] Oxford Nanopore Technologies, Tombo documentation, 2018 . https://nanoporetech.github.io/tombo/.

[19] N. J. Loman , J. Quick , J. T. Simpson , A complete bacterial genome assembled de novo using only nanopore sequencing data , Nature Methods 12 ( 2015 ) 733 - 735 .

[20] H. Li , Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM , arXiv preprint arXiv:1303.3997 ( 2013 ).

[21] V. Boža , B. Brejová , T. Vinař, Improving nanopore reads raw signal alignment , arXiv preprint arXiv:1705.01620 ( 2017 ).

[22] M. David , L. J. Dursi , D. Yao , P. C. Boutros , J. T. Simpson , Nanocall: an open source basecaller for Oxford Nanopore sequencing data , Bioinformatics 33 ( 2016 ) 49 - 55 .

[23] R. A . Fisher, Statistical methods for research workers , Oliver and Boyd (Edinburgh) , 1925 .

[24] M. Sipiczki , E. Kajdacsi, Jaminaea angkorensis gen . nov., sp. nov., a novel anamorphic fungus containing an S943 nuclear smallsubunit rRNA group IB intron represents a basal branch of microstromatales , International Journal of Systematic and Evolutionary Microbiology 59 ( 2009 ) 914 - 920 .

[25] B. Brejová , H. Lichancová , F. Brázdovič , E. Hegedűsová, M. F. Jakubková , V. Hodorová , V. Džugasová , A. Baláž , L. Zeiselová , A. Cillingová , M. Neboháčová , V. Raclavský , Ľubomír Tomáška , B. F. Lang , T. Vinař , J. Nosek , Genome sequence of the opportunistic human pathogen Magnusiomyces capitatus , Current Genetics 65 ( 2019 ) 539 - 560 .

[26] H. Li , Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM , arXiv preprint arXiv:1303.3997 ( 2013 ).

[27] Z. W. Yuen , A. Srivastava , R. Daniel , D. McNevin , C. Jack , E. Eyras, Systematic benchmarking of tools for CpG methylation detection from nanopore sequencing , Nature Communications 12 ( 2021 ) 3438 .

[28] Q. Liu , L. Fang , G. Yu , D. Wang , C. L. Xiao , K. Wang , Detection of DNA base modifications by deep recurrent neural network on Oxford Nanopore sequencing data , Nature Communications 10 ( 2019 ) 2449 .