Detection of DNA base modifications by deep recurrent neural network on Oxford Nanopore sequencing data — TEI-tagged view

Qian LiupersName, Li FangpersName, Guoliang YupersName, Depeng WangpersName, Chuan-Le XiaopersName, Kai WangpersName (2019date). “Detection of DNA base modifications by deep recurrent neural network on Oxford Nanopore sequencing data.” Nature Communicationstitle·j.

Object of study — DeepModterm·mth, a bidirectional LSTM RNNterm·mth over raw Oxford Nanopore electric signals

Research question — how to detect DNA modificationsterm·key where bisulfite and PacBio sequencing have inherent limitations

Contribution — DeepModterm·mth performs well for genome-scale detection of DNA modificationsterm·key across species

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Detection of DNA base modificationsterm·dom by deep recurrent neural networkterm·mth on Oxford Nanopore sequencingterm·key data Abstract DNA base modificationsterm·dom, such as C5-methylcytosineterm·dom (5mCterm·dom) and N6-methyldeoxyadenosineterm·dom (6mAterm·dom), are important types of epigeneticterm·key regulations. Short-read bisulfite sequencingterm·mth and long-read PacBio sequencing have inherent limitations to detect DNA modificationsterm·key. Here, using raw electric signals of Oxford Nanopore long-read sequencing data, we design DeepModterm·mth, a bidirectional recurrent neural networkterm·mth (RNNterm·mth) with long short-term memoryterm·mth (LSTMterm·mth) to detect DNA modificationsterm·key. We sequence a human genome HX1 and a Chlamydomonas reinhardtii genome using Nanopore sequencingterm·key, and then evaluate DeepModterm·mth on three types of genomes ( Escherichia coli , Chlamydomonas reinhardtii and human genomes). For 5mCterm·dom detection, DeepModterm·mth achieves average precisionterm·mth up to 0.99 for both synthetically introduced and naturally occurring modifications. For 6mAterm·dom detection, DeepModterm·mth achieves ~0.9 average precisionterm·mth on Escherichia coli data, and have improved performance than existing methods on Chlamydomonas reinhardtii data. In conclusion, DeepModterm·mth performs well for genome-scale detection of DNA modificationsterm·key and will facilitate epigeneticterm·key analysis on diverse species.

Introduction DNA base modificationsterm·dom are modified versions of standard nucleotides in a DNA molecule through the addition of chemical groups. For example, DNA methylations, including 5-methylcytosine (5mCterm·dom), 5-hydroxymethylcytosine (5hmC), and N6-methyldeoxyadenosineterm·dom (6mAterm·dom), are introduced into a DNA molecule by adding methyl or hydroxymethyl groups to nucleotides: 5mCterm·dom is introduced by adding methyl group at the 5th position of the pyrimidine ring of cytosines, whereas 6mAterm·dom is introduced by adding methyl group at the 6th position of the purine ring in adenines. There are also other different types of methylations, named according to the nucleotide type, the added molecule type and the modified position in nucleotides, such as 4mC (N4-methylcytosine), 5hmC (5-hydroxymethylcytosine), 5fC (5-formylcytosine), and 5caC (5-carboxylcytosine) introduced by methyltransferases 1 , and they generally do not disrupt base pairing. Some other modifications, such as 1mA (N1-methyladenine), 3mA (N3-methyladenine), 7mA (N7-methyladenine), 3mC (N3-methylcytosine), 2mG (N2-methylguanine), 6mG (O6-methylguanine), 7mG (N7-methylguanine), 3mT (N3-methylthymine), and 4mT (O4-methylthymine), may damage hydrogen bonds in base pairs 1 , 2 . DNA base modificationsterm·dom widely exist in different organisms 3 and are essential in various biological processes 4 , 5 , such as genomic imprinting, X-chromosome inactivation, genome stability, gene regulation, repression of transposable elements, aging, and carcinogenesis. For example, hyper-methylation of 5mCterm·dom in promoter regions usually represses gene transcription, and thus may regulate cellular processes such as cell differentiation and pluripotency. Local DNA 5mCterm·dom hyper-methylation and genome-wide hypo-methylation have been seen in cancer 6 , and the differential methylation of CpGterm·key islands can distinguish cancer cells from normal cells 7 or different tumor types 8 , and thus may be a potential cancer biomarker or a therapeutic target.

DNA modificationsterm·key, especially DNA methylations, could be detected by both short-read sequencing and long-read sequencing techniques. Bisulfite sequencingterm·mth with short-read techniques is widely used to call methylated cytosines 9 by converting un-methylated cytosines to uracil, but bisulfite sequencingterm·mth and its improved variations 10 were limited by the efficiency in conversion and the inability to assay repetitive genomic regions by short-read sequencing. Immunoprecipitation followed by short-read sequencing can detect DNA or RNA modifications 11 , 12 in genomic regions but it lacks single-base resolution. PacBio long-read sequencing can be used to directly detect DNA/RNA modifications 13 based on altered polymerase kinetics during sequencing 14 – 19 . However, this method has low signal-to-noise ratio for 5mCterm·dom modifications 20 , requires relatively high coverage for calling modifications, and is biased by incomplete and context-dependent enzymatic treatment of 5mCterm·dom detection using Tet1 21 .

Recently, several proof-of-concept studies on Oxford Nanopore sequencingterm·key techniques have demonstrated the feasibility to detect DNA modificationsterm·key based on the electric signal characteristics when a modified DNA molecule passes through nanopores 22 , 23 . DNA methylations at specific genomic positions can be measured with higher accuracy by comparing raw electric signals of methylated DNA copies with signals of the same un-methylated DNA copies 22 – 25 for known sequences without large prior training data set 3 . Furthermore, DNA methylation detection using Nanopore sequencingterm·key was used in a few studies 4 , 26 , 27 with the help of machine learning algorithms. Simpson et al. developed a HMMterm·mth (hidden Markov modelterm·mth) to measure long reads with complete methylations and long reads without methylations separately, and used the log-likelihood ratio to distinguish 5mCterm·dom from cytosine in Escherichia coli and Homo sapiens 4 . Similarly, Rand et al. used HMMterm·mth with a hierarchical Dirichlet process to identify three types of cytosine methylations (i.e., cytosine, 5mCterm·dom and 5hmC) and also 6mAterm·dom in E. coli with different phases 20 . Mclntyre et al. further improved the detection of 6mAterm·dom (implemented in mCallerterm·mth) in mouse, E. coli and Lambda phage DNA 26 using signal deviation of six 6-mer around positions of interest as input.

However, the sequential characteristic of Nanopore electric signals was not fully utilized in previous studies, and the performance of 5mCterm·dom/6mAterm·dom prediction may be improved using more sophisticated deep neural networks. In the current study, we developed a computational tool, DeepModterm·mth (Fig. 1 ), which takes a reference genome and long-read electric signals together with event information generated from Nanopore sequencingterm·key as input, and outputs modification summary for genomic positions of interest in a reference genome, together with modification prediction for bases of interest in a long read. The modification prediction model in DeepModterm·mth is a well-trained bidirectional recurrent neural networkterm·mth (RNNterm·mth) with long short-term memoryterm·mth (LSTMterm·mth) 28 units, which takes signal mean, standard deviation, and the number of signals of an event together with base information in the reference genome of an event and its neighbors as input, and makes modification prediction for the event. Then, after anchoring events with a reference genome based on the alignment of long reads, predicted modification summary for reference positions of interest can be generated in a BED format, indicating how many reads cover each genomic position and how many reads contain predicted modifications at genomic positions of interest. The prediction of DNA modificationterm·key by DeepModterm·mth is thus strand-sensitive and has single-base resolution. We note that one recent preprint 29 also used the RNNterm·mth framework (DeepSignal) but with limited evaluations. In the current study, we sequence a human genome HX1 and a Chlamydomonas reinhardtii genome using Nanopore sequencingterm·key techniques, and together with published Nanopore data for Escherichia coli and another human genome NA12878, we evaluate DeepModterm·mth on three types of genomes ( E. coli , C. reinhardtii, and human genomes) and show that it performs well on genome-scale detection of DNA modificationsterm·key. Fig. 1 The flowchart of DeepModterm·mth. RNNterm·mth: recurrent neutral network; LSTMterm·mth: long short-term memoryterm·mth. Several long reads (in yellow) were shown for demonstration purposes with their alignment and Nanopore signals (yellow lines). To view LSTM RNNterm·mth for modification prediction, nine LSTMterm·mth cells adjacent to the center position (with arrows in black) were shown Results A brief summary of the DeepModterm·mth algorithm We designed DeepModterm·mth to capture time-series characteristic of Nanopore signals for detecting DNA modificationsterm·key. DeepModterm·mth takes FAST5 files and a reference genome as inputs, where the FAST5 files are generated from Nanopore sequencer and must include event information after basecalling. The output files of DeepModterm·mth include BED files with coverage and methylation percentage information for genomic positions of interest. DeepModterm·mth works through several steps as illustrated in Fig. 1 .

First, a sequence of bases is extracted from events in each FAST5 file, and then aligned with a reference genome. Second, the raw signals of the events and the mapped reference base type of the events is obtained, and the summary of the raw signals of each of w events (by default, w = 21) with the reference base type of events (A or C or G or T) is used as input of a LSTM RNNterm·mth to predict whether the center event of the w events is generated from a modified base (the left part of Fig. 1 ). After repeating the prediction process by the LSTM RNNterm·mth for events of interest in a long read and then for all long reads, the sequence coverage and the methylation coverage are generated for genomic positions of interest in the reference genome. The whole-genome modification summary is formatted as BED files for performance evaluation and other analysis. In addition, a second neural network is also used for some types of modifications, such as 5mCterm·dom 30 to incorporate high correlation of methylation of a CpGterm·key site and its nearby CpGterm·key sites. The second neural network takes the predicted methylation percentage of a reference position of interest in the genome and the methylation percentage of its neighboring sites in both strands as input and generates a new methylation percentage of the position of interest. In this study, DeepModterm·mth was evaluated on three types of genomes with two types of modifications (5mCterm·dom and 6mAterm·dom), and the performance evaluations were described in detail below.

Performance of 5mCterm·dom detection on Nanopore sequencingterm·key data on E. coli DeepModterm·mth was first trained on a published Nanopore sequencingterm·key data on E. coli with positive and negative controls 4 (See the Methods section and Supplementary Notes ). The negative control (UMR for short) were amplified by polymerase chain reaction (PCR) and thus the reads contain no modifications; in comparison, reads in the positive control (CG_MSssI for short) were treated by the M.SssI methyltransferase after PCR amplification and thus the majority of CpGterm·key sites were methylated 2 . The summary of the two data sets was described in Table 1 . To select suitable hyper-parameters in DeepModterm·mth for the tradeoff between prediction performance and running time/consuming resources, two independent validation strategies were used on UMR and CG_MSssI. The first strategy is read-based-independent validationterm·mth, where 90% of the reads from the negative control and 90% of the reads from the positive control were randomly selected and used to train DeepModterm·mth, whereas the rest (10%) of the reads were used for testing. Under this validation strategy, reads in the training and testing groups might be aligned to the same genomic regions. In parallel, in the second region-based independent validation strategy, reads mapped to the genomic positions from 1,000,000 to 2,000,000 of E. coli reference genome were used for testing, whereas reads mapped to other genomic regions were used for training. We also evaluated two additional hyper-parameters in DeepModterm·mth: different w (i.e., the number of LSTMterm·mth units ranging from 7 to 51), 7-feature description of an event (base information, signal mean, standard deviation, and the number of signals associated with an event) and 57-feature description of an event (7-features description plus 50 discretized bins of signals). The performance was measured using accuracy, precisionterm·def, recall, and F1-scoreterm·mth (See Eqs. ( 3 – 6 ) in the Methods section for their definitions). The results of the two independent validation above were shown in Supplementary Table 1 for per-call performance (where the prediction of a base in each long read was considered separately without alignment). We found that (i) 57-feature description produced similar performance to 7-feature description, and (ii) larger w generated increasing accuracy, recall and F1-scoreterm·mth, but the memory usage and running time were much higher beyond w = 21. We thus used 7-feature description and w = 21 as the default setting in DeepModterm·mth. Cross-validation of region-based independent validation strategy was also provided in Supplementary Table 2 for DeepModterm·mth with 7-feature description and w = 21, and the performance values in five different regions ([0, 1000000], [1000000, 2000000], [2000000, 3000000], [3000000, 4000000], [4000000, 4700000]) were similar to each other. Table 1 Nanopore sequencingterm·key data sets used to evaluate DeepModterm·mth Genome Data set name Motif # reads Coverage Meth a Reference Escherichia coli UMR NA d 111,238 110X Neg b Simpson et al. 4 CG_MsssI C G c 69,899 67X 5mCterm·dom CG_SssI C G c 8679 19X 5mCterm·dom CG_MpeI C G c 23,593 39X 5mCterm·dom GCGC_HhaI G C GC c 18,180 50X 5mCterm·dom gaAttc_EcoRI GA A TTC c 16,661 27X 6mAterm·dom Stoiber et al. 3 gAtc_dam G A TC c 17,557 33X 6mAterm·dom tcgA_TaqI TCG A c 16,249 22X 6mAterm·dom Con1 NA d 23,762 34X Neg b Con2 NA d 34,170 40X Neg b Homo sapiens NA12878 C G c 30X 5mCterm·dom Jain et al. 31 HX1 C G c 4,827,155 30X 5mCterm·dom Current study Chlamydomonas reinhardtii C. reinhardtii NA d 772,817 126X 6mAterm·dom Current study a Methylation types b Negative control without any modifications c Underlined nucleotides in motifs were potential modified target. d No modifications or no motif information DeepModterm·mth’s cross-data set performance was further evaluated on another independent E. coli Nanopore data set 3 where there were two negative control samples without any modifications (denoted as Con1 and Con2 for short), and three positive control samples with 5mCterm·dom modifications. On the three positive control samples, three methylases (M.Hhal, M.Mpel, and M.Sssl) were used separately to synthetically introduce 5mCterm·dom for 70,180 G C GC motifs (the underlined bases represents modified cytosines by enzymes hereafter) and 693,586 C G motifs in E. coli , and thus the samples were denoted as GCGC_Hhal, CG_MpeI, and CG_SssI for short. DeepModterm·mth was used to make the modification predictions on the five Nanopore sequencingterm·key data sets (the three positive control data with 5mCterm·dom and two negative control without modifications). We mixed each of three positive control data with both Con1 and Con2 for performance evaluation of DeepModterm·mth.

There were two types of evaluations on this independent E. coli data set. One type of evaluation was based on methylated cytosines within specific sequence motifs of interest. Under this type of evaluation, the number of the cytosines of interest in motifs from positive control is roughly equal to that from Con1 and Con2 as shown in Table 2 . In Table 2 , there are ~ 1.4 million CpGterm·key sites and ~ 0.14 million GCGC sites, where about half were 5mCterm·dom from positive control and the other half have no modification from the negative control. The performance of DeepModterm·mth was then evaluated using AUCterm·mth and APterm·mth (AUCterm·mth: area under curveterm·mth, and APterm·mth: average precisionterm·mth. See more details in Methods). The performance of DeepModterm·mth was shown in Fig. 2a and the confusion matrices were provided in Supplementary Table 3 . It can be seen from Fig. 2a that DeepModterm·mth achieved APterm·mth values of 0.990, 0.921 and 0.993, and AUCterm·mth values of 0.987, 0.906, and 0.988 for GCGC_Hhal, CG_MpeI, and CG_SssI, respectively, for all motif sites of interest (Fig. 2a ) for coverage ≥ 1. In particular, given a threshold of prediction methylation percentage ≥ 0.1 for a genomic position of interest, DeepModterm·mth gave precisionterm·def = 0.945 and recall = 0.97 for GCGC_Hhal, precisionterm·def = 0.835 and recall = 0.879 for CG_MpeI and precisionterm·def = 0.848 and recall = 0.986 for CG_SssI. As CG_SssI (19X) has much lower coverage than GCGC_Hhal (50X) and C G_MpeI (39X), if the threshold of prediction methylation percentage was set to a larger value for CG_SssI, i.e., ≥ 0.2, DeepModterm·mth achieved precisionterm·def = 0.96 and recall = 0.985 with significant improvement of precisionterm·def but a slight (0.001 point) decrease of recall. Meanwhile, a larger coverage threshold such as 5 can yield better performance, although the improvement is not significant as shown in Fig. 2a. Table 2 The number of modified and un-modified bases of interest used for evaluation on E. coli when coverage ≥ 1 Data set name Base of interest Modification # modified base in motif # non-modified base in motif # base not in motif CG_MpeI cytosine (C) 5mCterm·dom 693,518 693,427 3,326,971 CG_SssI 682,526 693,427 3,297,528 GCGC_HhaI 70,172 70,160 4,573,777 gaAttc_EcoRI adenine (A) 6mAterm·dom 277 280 1,003,885 gAtc_dam 7816 7831 989,363 tcgA_TaqI 6351 6403 989,375 Fig. 2 Evaluation of the performance of DeepModterm·mth on 5mCterm·dom prediction on E. coli , NA12878, and HX1. a APterm·mth (the outer) and AUCterm·mth (the inner) plots for 5mCterm·dom within sequence motifs in E. coli for three synthetically introduced 5mCterm·dom data sets by M.Mpel (CG_MpeI for CG motif), M.Sssl (CG_SssI for CG motif), and M.Hhal (GCGC_HhaI for GCGC motif), respectively. b , c AUCterm·mth and APterm·mth plots for 5mCterm·dom prediction of all cytosines in E. coli . d , e APterm·mth and AUCterm·mth of 5mCterm·dom prediction by DeepModterm·mth on NA12878. f , g APterm·mth and AUCterm·mth of 5mCterm·dom prediction by DeepModterm·mth on HX1. Cov: coverage. # base: total number of bases used in the evaluation Besides the evaluation on methylated cytosines within sequence motifs of interest above, the second type of evaluation on the independent E. coli data were based on all cytosines from both positive control and negative control, and the results were illustrated in Fig. 2b, c : DeepModterm·mth achieved AUCterm·mth of 0.985, 0.953, and 0.992 for all cytosines in the GCGC_Hhal data set, the CG_MpeI data set, and the CG_SssI data set, respectively, suggesting similar AUCterm·mth values when moving the classification of motif (GCGC or CG)-based cytosines to the classification of all cytosines. However, from Fig. 2c , the APterm·mth values of DeepModterm·mth to classify all cytosines are 0.771, 0.890, and 0.983 in the GCGC_Hhal data set, the CG_MpeI data set and the CG_SssI data set, respectively, which were lower than the performance for the prediction of cytosines in sequence motifs, indicating misclassifications for cytosines that are not present in motifs. When we checked which cytosines might be more likely to be misclassified, we found that cytosines closely adjacent to modified cytosines (i.e., one or two upstream/downstream cytosines) has much higher likelihood to be misclassified than cytosines in other positions, suggesting possible neighborhood effect caused by modifications. In fact, the misclassification of cytosines adjacent to modified cytosines may not be wrong but merely indicates the existence of modified cytosines that are adjacent to each other, and thus the prediction power of our method may be under-estimated in these scenarios. In summary, this analysis on several independent data sets clearly demonstrated the satisfactory performance of DeepModterm·mth in genome-wide, single-base resolution detection of 5mCterm·dom.

Performance of 5mCterm·dom prediction on independent Nanopore sequencingterm·key data on human genomes Cross-species 5mCterm·dom prediction performance of DeepModterm·mth was further evaluated on a human Nanopore sequencingterm·key data set. Under this evaluation, a DeepModterm·mth model, which was trained on E. coli data, were used to make 5mCterm·dom prediction on the NA12878 human cell line. Recently, NA12878 was sequenced by Jain et al. mainly using Nanopore R9.4 with ~ 30X coverage 31 . On this human Nanopore sequencingterm·key data, the well-trained DeepModterm·mth model on E. coli data were used to make 5mCterm·dom prediction for each Nanopore long read, and each cytosine in a long read was associated with a label indicating whether it was predicted to be 5mCterm·dom or not. Then, after all long reads were aligned against a reference genome, predicted methylation summary for the whole-genome was generated for each cytosine in the genome with a percentage to indicate how many modification predictions were made for that position, compared with the number of reads mapped to that position. Furthermore, to account for the high correlations of nearby 5mCterm·dom in CpGterm·key sites 30 , we designed another deep neural network, and the second neural network incorporates the predicted methylation percentages of a CpGterm·key site, the prediction on the opposite strand and the methylation prediction of its neighboring CpGterm·key sites, and outputs a final prediction percentage for a genomic position. This second deep neural network could improve APterm·mth by 1–3% points and AUCterm·mth by 3–5% for all chromosomes when coverage ≥ 1.

To evaluate DeepModterm·mth prediction performance on NA12878, we used 5mCterm·dom calls from bisulfite sequencingterm·mth of NA12878 as benchmark 32 . Owing to the heterogeneity of sequenced samples, we used completely methylated and completely un-methylated bases for the evaluation: a genomic position of a base was considered to be completely methylated if its methylation percentage ≥ 90% in both replicates of bisulfite sequencingterm·mth with coverage ≥ c ( c could be 1, 5, or 10), and to be completely un-methylated if its methylation is 0% in both replicates. The number of completely methylated and un-methylated bases for each chromosome in NA12878 were given in Supplementary Table 4 . We compared methylation percentage of DeepModterm·mth prediction against the completely methylated and un-methylated cytosines at reference positions of interest for each chromosome. The results on APterm·mth were shown in Fig. 2d , and the AUCterm·mth was given in Fig. 2e , where both APterm·mth and AUCterm·mth with coverage ≥ 1 were ~0.99, yet when coverage ≥ 5, the values of APterm·mth and AUCterm·mth were > 0.99. For example, on the chromosome 1, when coverage ≥ 5, DeepModterm·mth achieved the APterm·mth value of 0.998 and the AUCterm·mth value of 0.993, while on the chromosome 20, when coverage ≥ 5, DeepModterm·mth had the APterm·mth value of 0.998 and the AUCterm·mth value of 0.995, indicating accurate prediction on human genome, even when using a DeepModterm·mth model trained on E. coli data.

With NA12878 Nanopore data, DeepModterm·mth was also evaluated on those 5mCterm·dom sites where the methylation percentage has small difference (< 0.25) between the two replicates of bisulfite sequencingterm·mth to account for the heterogeneity of two bisulfite sequencingterm·mth batches. We found that the prediction methylation percentage has high correlation with the methylation percentage in both replicates (> 0.85 for the majority of chromosomes when coverage ≥ 5). We also noted that DeepSignal 29 trained on E. coli tried to make cross-species prediction of CpGterm·key methylation on NA12878, which achieved F1 = 0.949 (with precisionterm·def = 0.97 and recall = 0.93) on synthetically introduced 5mCterm·dom in human, whereas nanopolish achieved F1 = 0.89 (with precisionterm·def = 0.944 and recall = 0.844). Both methods required a CpGterm·key site with signals in both the template and complement strands and only output final prediction on template strand by summarizing predictions from both strands, and the prediction was based on synthetically introduced 5mCterm·dom (complete methylation and un-methylation). Our method on native 5mCterm·dom prediction in NA12878 achieved F1 = 0.983 (with precisionterm·def = 0.988 and recall = 0.979) and F1 values on each individual chromosome (calculated based on precisionterm·def and recall in Supplementary Table 4 ) were ~ 0.98. Although this is not a direct comparison, DeepModterm·mth still demonstrated its highly accurate prediction on NA12878.

DeepModterm·mth was further evaluated on another independent human data set HX1, by sequencing whole-blood sample 33 with ~ 30X coverage on the Nanopore platform and by two replicates of bisulfite sequencingterm·mth of the same sample. To make 5mCterm·dom prediction on HX1, we trained a model using completely methylated and completely un-methylated cytosines of chromosome 1–10 of NA12878, so that NA12878 and HX1 were basecalled with the same version of Albacore. Then, similar to the evaluation on NA12878, we compared 5mCterm·dom prediction by DeepModterm·mth on Nanopore sequencingterm·key data against completely methylated and un-methylated cytosines determined from two replicates of bisulfite sequencingterm·mth data for HX1 (see Supplementary Table 5 for details). After that, we calculated APterm·mth (Fig. 2f ) and AUCterm·mth (Fig. 2g ) and precisionterm·def and recall (Supplementary Table 5 ) for coverage ≥ 3, where the APterm·mth of DeepModterm·mth is > 0.99 and the AUCterm·mth is generally higher than 0.97 except on chromosomes X and Y. The performance was comparable to nanopolish, but both methods have already achieved excellent 5mCterm·dom methylation prediction on HX1. For example, on the chromosome 1, the APterm·mth value of DeepModterm·mth is 0.998 and the AUCterm·mth value is 0.979, demonstrating highly accurate predictions. We note that nanopolish was designed for those CpGterm·key sites which were completely methylated/un-methylated CpGterm·key sites within 10 bps and not semi-methylated, while DeepModterm·mth does not have this limitation.

Cross-species independent evaluation of 5mCterm·dom prediction on E. coli To further evaluate cross-species performance of DeepModterm·mth on the E. coli data, we trained a model using chromosome 1–10 of NA12878 (where both NA12878 and E. coli data were basecalled with the same version of Albacore). We then used DeepModterm·mth to make methylation prediction on several E. coli data sets, including Con1, Con2, GCGC_Hhal, CG_MpeI, and CG_SssI 3 . For each of three positive data set (GCGC_Hhal or CG_MpeI or CG_SssI), we mixed motif-based 5mCterm·dom with the corresponding motif sites in both Con1 and Con2 for performance evaluation of DeepModterm·mth, and the number of the cytosines of interest in motifs from positive control and from Con1 and Con2 were shown in Table 2 . With a coverage ≥ 1, DeepModterm·mth achieved APterm·mth = 0.95 on GCGC_Hhal, APterm·mth = 0.807 on CG_MpeI, and APterm·mth = 0.973 on CG_SssI. The AUCterm·mth values achieved by DeepModterm·mth are 0.942, 0.809, and 0.965 on GCGC_Hhal, CG_MpeI, and CG_SssI, respectively. On GCGC_Hhal and CG_SssI, the APterm·mth values were 2–4% points lower than the prediction by DeepModterm·mth trained on the E. coli data of CG_MSssI and UMR 4 , whereas the APterm·mth value (0.807) was 0.11 lower on CG_MpeI. In summary, the cross-species testing suggested that DeepModterm·mth trained on one species can make accurate prediction on the other species, when the same basecalling algorithm was used. We note that the data from NA12878 was generated from native DNA while the data from Stoiber et al. 3 were treated by PCR and enzymes: native DNAs may contain incompletely methylated cytosines due to cellular heterogeneity, which may result in slightly lower performance when training a model.

Performance of 6mAterm·dom prediction on Nanopore sequencingterm·key data on E. coli To evaluate DeepModterm·mth on 6mAterm·dom modifications within specific sequence motifs, we used three positive control data sets by Nanopore sequencingterm·key of E. coli with synthetically introduced 6mAterm·dom 3 , where M.Taql, M.EcoRI and M.dam were used separately to methylase 30,914 TCG A motifs (denoted tcgA_TaqI for short), 1290 GA A TTC motifs (denoted gaAttc_EcoRI for short), and 38,240 G A TC motifs (denoted gAtc_dam for short), respectively. With region-based independent validation strategy, we evaluated DeepModterm·mth for the genomic positions from 1,000,000 to 2,000,000, whereas reads mapped outside of this region was used for training. After predicting methylation percentage for bases of interest in this region, we calculated APterm·mth and AUCterm·mth values by combining each of the three 6mAterm·dom positive control with Con1 and Con2 in Fig. 3a for those adenines in motifs, and the number of 6mAterm·dom in motifs and the corresponding non-methylated adenines were given in Table 2 . Figure 3a showed that DeepModterm·mth achieved APterm·mth = 0.884 and AUCterm·mth = 0.874 on gaAttc_EcoRI, APterm·mth = 0.858, and AUCterm·mth = 0.857 on tcgA_TaqI, and APterm·mth = 0.913 and AUCterm·mth = 0.903 on gAtc_dam with coverage ≥ 5, suggesting that DeepModterm·mth performed well for 6mAterm·dom modification detections. In particular, given a threshold of predicted methylation percentage ≥ 0.1 for a genomic position of interest, DeepModterm·mth achieved precisionterm·def = 1.0 and recall = 0.788 on gaAttc_EcoRI, precisionterm·def = 0.859 and recall = 0.841 on tcgA_TaqI, and precisionterm·def = 0.748 and recall = 0.929 on gAtc_dam. We also conducted cross-validation of the five different regions in E. coli reference genome (i.e., [0, 1000000], [1000000, 2000000], [2000000, 3000000], [3000000, 4000000], [4000000, 4700000]), and DeepModterm·mth achieved APterm·mth = 0.83 ± 0.037 on gaAttc_EcoRI, APterm·mth = 0.89 ± 0.073 on gAtc_dam, and APterm·mth = 0.81 ± 0.041 on tcgA_TaqI, showing varying but similar performance in different genomic regions. Fig. 3 Evaluation of the performance of DeepModterm·mth on 6mAterm·dom prediction on E. coli and C. reinhardtii. a APterm·mth (the outer) and AUCterm·mth (the inner) plots on E. coli for three synthetically introduced 6mAterm·dom data sets by EcoRI (gaAttc_EcoRI for GAATTC motif), TaqI (tcgA_TaqI for TCGA motif), and dam (gAtc_dam for GATC motif), respectively. b 6mAterm·dom prediction by DeepModterm·mth on C. reinhardtii (the outer plot is for genomic DNA digested in 0.5 h, whereas the inner plot is for genomic DNA digested in 12 h 34 ). Cov: coverage. # base: total number of bases used in the evaluation We next evaluated the performance of our methods on all methylated adenines regardless of whether adenines are within specific sequence motifs or not. When we calculated AUCterm·mth using the predictions for all adenines, the AUCterm·mth values increased from 0.874 to 0.898 for gaAttc_EcoRI, 0.857 to 0.92 for tcgA_TaqI, and 0.903 to 0.973 for gAtc_dam, indicating lower overall prediction errors. However, the APterm·mth values decreased. For example, for coverage ≥ 5, when we classified all adenines rather than gAtc-based adenines, the APterm·mth value decreased to 0.874, whereas the APterm·mth value decreased to 0.66 when we classified all adenines rather than tcgA-based adenines. This suggested a fraction of false positive predictions of 6mAterm·dom not in motifs of interest, compared with a small number of true positive modifications. For example, there was only 277 6mAterm·dom on gaAttc_EcoRI and 1,003,885 un-modified adenines, and even 1% misclassification of un-modified adenines is much larger than the number (277) of 6mAterm·dom, and will lower the overall APterm·mth values.

Cross-species independent evaluation of 6mAterm·dom prediction on C. reinhardtii As C. reinhardtii is among the first eukaryotes 19 , 34 , 35 with detailed characterization of 6mAterm·dom modifications, to further evaluate 6mAterm·dom prediction of DeepModterm·mth, we sequenced C. reinhardtii using Nanopore sequencingterm·key with ~ 126X coverage, and then used the model trained on E. coli to make 6mAterm·dom prediction for C. reinhardtii . The predictions were compared with the existing work of DA-6mA-seq results 34 , 35 where genomic DNA was digested in 0.5 h and 12 h. We plotted AUCterm·mth for the DeepModterm·mth prediction together with the prediction by mCallerterm·mth 26 with coverage ≥ 15 in Fig. 3b . mCallerterm·mth used current deviations of six consecutive 6-mers around a position of interest as input of a neural network for 6mAterm·dom prediction and it is known to be a highly accurate 6mAterm·dom predictor. In this evaluation, mCallerterm·mth achieved 0.74 AUCterm·mth values for both DA-6mA-seq results, and DeepModterm·mth achieved slightly higher AUCterm·mth value of 0.77 AUCterm·mth. When we used coverage > 50 for both tools (since this data has ~ 126X coverage), DeepModterm·mth achieved AUCterm·mth = ~ 0.8, whereas mCallerterm·mth achieved AUCterm·mth = ~ 0.75 for both DA-6mA-seq results. In summary, DeepModterm·mth can make accurate inferences on detection of 6mAterm·dom modifications in cross-species evaluation studies.

Discussion DeepModterm·mth developed in this study bridges the gap between the rapid growth of Nanopore sequencingterm·key data and the increasing need of detecting DNA modificationsterm·key at a genomic scale. DNA modificationsterm·key are not consistent through biological cell cycle and may vary between cell types. DNA modificationsterm·key, such as 5mCterm·dom, could be introduced by methyltransferases during cell development, and might also be de-methylated in some other biological processes such as cell differentiation. Specific enzymes are necessary to methylate cytosine or de-methylate 5mCterm·dom, which plays a critical role in epigeneticterm·key reprogramming, whereas aberrant DNA modificationsterm·key have been found to play important roles in human diseases such as various cancers 36 . Although bisulfite sequencingterm·mth and PacBio long-read techniques have achieved great progress to detect some types of DNA modificationsterm·key, especially 5mCterm·dom, the advent of Nanopore long-read sequencing has significant potential to detect major types of DNA modificationsterm·key from Nanopore raw electric signals. However, Nanopore sequencingterm·key data contains context-dependent electric signals which complicates the ability to pinpoint where modifications are located. Therefore, in this study, we propose a deep learning framework, DeepModterm·mth, to detect DNA modificationsterm·key from raw electric signals of events (together with mapped base types) from Nanopore sequencingterm·key in a context-dependent manner via the use of LSTMterm·mth. Evaluation on three types of species with two major type of DNA modificationsterm·key demonstrates that DeepModterm·mth performs accurately to identify both types of DNA modificationsterm·key and has great potential to be extended to other types of DNA modificationsterm·key. DeepModterm·mth thus can be used as a genomic tool to detect genomic modifications for the increasing volume of Nanopore sequencingterm·key data without additional sequencing cost.

There are several current limitations in DeepModterm·mth that we wish to discuss here. First, DeepModterm·mth is trained and evaluated on 5mCterm·dom and 6mAterm·dom, two common DNA modificationsterm·key, and it might not work well with other types of modifications or other different motifs before sufficient amounts of training data is available. Nevertheless, DeepModterm·mth provides a simple extension framework, and the model of other types of DNA modificationsterm·key could be easily trained. This limitation could be further addressed when unsupervised or one-class learning are used to build a prediction model in the future, with a tradeoff that the specific modification types might not be predicted. Second, DeepModterm·mth can only work on DNA data of Nanopore sequencingterm·key currently. RNA data of direct Nanopore sequencingterm·key contains natural RNA modifications, which are informative to interpret gene expression and gene regulation, but large-scale ground-truth RNA modification data are currently lacking for supervised machine learning. In addition, uracil in RNA may make Nanopore electric signals in RNA data different from that in DNA data, and thus additional work is needed to see whether the DeepModterm·mth framework needs modification to handle direct RNA Nanopore sequencingterm·key data. Third, DNA modificationsterm·key could affect Nanopore signals of several neighboring bases, and the impacts from nearby modifications within a few bases away may influence the prediction of a given base, yet this fact is not currently considered in DeepModterm·mth. Fourth, DeepModterm·mth relied on alignment tool to find correct reference positions of bases in long reads. Now, DeepModterm·mth supports two widely used aligners: BWA-MEM and minimap2 32 . If different aligners or improper parameters are used with poor alignment performance, DeepModterm·mth’s prediction will also be affected. Re-alignments with different tools or parameters might bring in new examined sites or exclude previously examined sites. This effect contributes much less to false modification prediction for those frequent modification types, but may become significant for those rarely occurring modifications. However, the majority of those limitations can be addressed by having more ground-truth modifications data sets for training the prediction models with proper parameters. When more and more modifications data sets become available, we believe that improved model in DeepModterm·mth could be learned to generate better prediction models.

In summary, DeepModterm·mth together with Oxford Nanopore sequencingterm·key provide a useful tool to identify DNA modificationsterm·key at a genomic scale. We expect that DeepModterm·mth can be used in whole-genome epigeneticsterm·key analysis with more types of modifications in future, with appropriately trained prediction models on additional gold standard data sets generated by the community.

Methods We developed DeepModterm·mth for the detection of DNA modificationsterm·key from Nanopore sequencingterm·key data. DeepModterm·mth is a supervised deep RNNterm·mth with Nanopore sequencingterm·key data as input. In this section, several Nanopore sequencingterm·key data sets used in this study were described below together with how to build, train and evaluate DeepModterm·mth on Nanopore sequencingterm·key data generated for E. coli , C. reinhardtii, and two human genomes (NA12878 and HX1).

Existing Nanopore sequencingterm·key data Several Nanopore sequencingterm·key data sets were used in this study, including seven published positive and three negative control data sets for E. coli sequenced by Simpson et al. 4 and by Stoiber et al. 3 , one published Nanopore sequencingterm·key data set for a human genome NA12878 sequenced by Jain et al. 31 , a new Nanopore sequencingterm·key data set on HX1 and a new sequencing data set on C. reinhardtii both generated in this study.

The positive control samples for E. coli contain four data sets for 5mCterm·dom and three data sets for 6mAterm·dom where methylations were synthetically introduced by specific enzymes, whereas three negative control samples were amplified by PCR and thus contain no modified bases. A positive control of 5mCterm·dom data set (with the short name CG_MsssI) and a negative control data set (with a short name of UMR) were sequenced by Simpson et al. 4 , whereas the three 5mCterm·dom data sets (with short names of GCGC_Hhal, CG_MpeI, and CG_SssI) and three 6mAterm·dom data sets (with short names of tcgA_TaqI, gAtc_dam, and gaAttc_EcoRI) and two negative control data sets (with short names of Con1 and Con2) were sequenced by Stoiber et al. 3 . The detail description of these data sets can be found in Supplementary Notes and 3 , 4 . The sequencing data sets on two human genomes and on C. reinhardtii were based on native DNA molecules, and thus contain native modifications. Below, we will introduce Nanopore sequencingterm·key process for HX1 and C. reinhardtii , whereas the Nanopore sequencingterm·key data of NA12878 can be downloaded online.

Nanopore sequencingterm·key data on HX1 HX1 is a well-studied genome from a Chinese individual who has been sequenced by Illumina short-read sequencing and PacBio long-read sequencing 33 . Native human genomic DNA was extracted from fresh blood, and size selection was performed using Blue Pippin (Cassette kit: BUF7510; size range: 30–40 kb). DNA quality was assessed by running 1 µl on a genomic ScreenTape on the TapeStation 2200 (Agilent) to ensure a DNA Integrity Number > 7. Concentration of DNA was assessed using the dsDNA HS assay on a Qubit fluorometer (Thermo Fisher).

For library preparation, 2.0 µg size-selected (> 20 kb) genomic DNA was used as the input DNA of each library. End repair and dA-tailing was performed using NEBNext Ultra II End Repair/dA-tailing Module (catalog No. E7546). In all, 7 µl Ultra II End-Prep buffer, 3 µl Ultra II End-Prep enzyme mix were added to the input DNA. The total volume was adjusted to 60 µl by adding nuclease-free water (NFW). The mixture was incubated at 20 °C for 5 min and 65 °C for 5 min. A 1 × volume (60 µl) AMPure XP clean-up was performed and the DNA was eluted in 31 µl NFW. One microliter of the eluted dA-tailed DNA was quantified using the Qubit fluorometer. A total of ≥ 1.0 µg DNA should be retained if the process is successful.

Adaptor ligation was performed using the following steps. Twenty microliter Adaptor Mix (ONT, SQK-LSK108 Ligation Sequencing Kit) and 50 µl NEB Blunt/TA Master Mix (NEB, catalog No. M0367) were added to the 30 µl dA-tailed DNA. The mixture was incubated at room temperature for 10 min. The adaptor-ligated DNA was cleaned up using 40 µl of AMPure XP beads. The mixture of DNA and AMPure XP beads was incubated for 5 min at room temperature and the pellet was washed twice by 140 µl ABB (SQK-LSK108). The purified-ligated DNA was resuspended in 15.5 µl ELB (SQK-LSK108). A 1-µl aliquot was quantified by fluorometry (Qubit) to ensure ≥ 500 ng DNA was retained. The final library was prepared by mixing 35.0 µl RBF (SQK-LSK108), 25.5 µl LBB (SQK-LSK108), and 14.5 µl purified-ligated DNA. The library was loaded to R9.4 flow cells (FLO-MIN106, ONT) according to the manufacturer’s guidelines. GridION sequencing was performed using default settings for the R9.4 flow cell and SQK-LSK108 library preparation kit. The sequencing was controlled and monitored using the MinKNOW software developed by ONT. Nanopore sequencingterm·key generated 4,827,155 FAST5 files in total after using Albacore v2.3.1 basecalling, and there were ~ 91G bases in total.

Nanopore sequencingterm·key data on C. reinhardtii C. reinhardtii has been used to study 6mAterm·dom modification using DA-6mA-seq techinques 34 where a region containing 6mAterm·dom modifications were detected. In the current study, we sequenced the same strain of C. reinhardtii as previously published, using Nanopore sequencingterm·key techniques as described above. We generated 772,817 FAST5 files after using Albacore v2.31 for basecalling, there were 15G bases in total, corresponding to ~ 126X genome-wide coverage.

DeepModterm·mth framework DeepModterm·mth is a deep learning tool of bidirectional RNNterm·mth with long short-term memoryterm·mth (LSTMterm·mth) units. LSTM RNNterm·mth is a class of artificial neural network for modeling sequential behaviors with LSTMterm·mth to preclude vanishing gradient problem, and has achieved superior performance in handwriting recognition 37 , speech recognition 38 , and computational biology 39 . LSTM RNNterm·mth would also be a better method to model series of raw signals generated in Nanopore sequencingterm·key.

In DeepModterm·mth, the input is a reference genome and FAST5 files generated by Nanopore sequencers with raw signals and base calls, and the output includes modification prediction for bases of interest in a long read and modification summary for genomic positions of interest in a reference genome in BED format. DeepModterm·mth contains several steps as shown in Fig. 1 : (i) the alignment of long reads to a reference genome, (ii) feature extraction from inputs, (iii) a deep learning framework for modification prediction, and (iv) an optional second neural network for considering methylation cluster effect of 5mCterm·dom in CpGterm·key sites, and (v) the modification summary for outputs. Each of the five steps was described below. (i) The input to DeepModterm·mth includes a reference genome and FAST5 files containing raw signals and events, which were generated by Nanopore sequencers with base calls. Each event is associated with a k-mer (e.g., 5-mer). In DeepModterm·mth, stay events without move were merged with the adjacent previous non-stay event with move to generate a single event. A sequence of bases associated with the resultant events was aligned to a reference genome using BWA-MEM with Nanopore-specific parameters 40 or minimap2. Each non-insertion base in long reads can be anchored with a genomic position in a reference genome for further analysis. (ii) In a FAST5 file, each event is associated with several sequential raw signals, and raw signals for all aligned bases in a long read were normalized using the method proposed by 3 and the normalized range was rescaled from −5 to 5. Then, the signal mean, standard deviation, and the number of signals associated with an event were extracted, plus a 4-vector feature for mapped reference bases. In the 4-vector feature, 1 indicates the mapped reference base is a specific nucleotide type, whereas 0 means otherwise. Thus, 7 features were used to describe an event (7-feature description for short and denoted by x i = [ f m , f d , f l , f A , f C , f G , f T ]). In the training and testing process, events were also labeled by modification or non-modification inferred from synthetically introduced modifications for motifs of interest or from bisulfite sequencingterm·mth (Please refer to the Supplementary Methods for more detail). (iii) The deep learning framework in DeepModterm·mth is a bidirectional RNNterm·mth of LSTMterm·mth units as shown in Fig. 1 . RNNterm·mth is used to capture the Nanopore sequencingterm·key characteristic of that a signal would be affected by several neighborhood bases, and LSTMterm·mth can overcome the effect of vanishing gradient problem in RNNterm·mth training process 28 . LSTM RNNterm·mth has been widely used in those fields where sequential order need to be preserved in training process of a deep learning framework such as handwriting recognition 37 , speech recognition 38 , and protein sequence analysis 41 . In DeepModterm·mth, bidirectional RNNterm·mth with LSTMterm·mth units was used with three hidden layers, where bidirectional RNNterm·mth were used to take into consideration both forward and reverse data flow from neighborhood bases, and each LSTMterm·mth unit contains multiple hidden nodes. In detail, an event and its ⌊ w /2⌋ upstream and ⌊ w /2⌋ downstream events (i.e., x = { x i - ⌊ w ∕ 2 ⌋ , ⋯ x i , ⋯ x i + ⌊ w ∕ 2 ⌋ } , where w is an odd number and ⌊ w /2⌋ is the floor of w /2) were used as the inputs of w LSTMterm·mth units in RNNterm·mth according to the sequential order in a long read, and 7 features ( x i ) of each event was the input of a LSTMterm·mth unit. In a LSTMterm·mth unit with x j as input, i −⌊ w /2⌋≤ j ≤ i +⌊ w /2⌋, then, the output p j is 1 p j = tanh ( f i 1 ∘ f f i 2 + v j - 1 ∘ f f ) ∘ f p where f p = sigmod ( f ( a p b p c p ) ) is activation function for output gate, f f = sigmod ( f ( a f b f c f ) ) is activation function for forgot gate, v j −1 is inner delay state, f i 1 = tanh ( f ( a i 1 b i 1 c i 1 ) ) is activation function for input, f i 2 = sigmod ( f ( a i 2 b i 2 c i 2 ) ) is activation function for input gate, f a b c = a b c * x j p j - 1 1 T , where superscript T is the tranpose of a matrix, p j −1 is the previous output, and ∘ is the element-wise multiplication. Then, both forward and backward data flow were captured as shown in Fig. 1 , and 3-layer RNNterm·mth with full connections was used to capture complicated relationship between signals and prediction target of modifications. Given a long read with sequential events, an event and its neighborhood could be used as input of this neural network, and a prediction label was generated for this event. This process can be repeated one by one for events of interest in a long read and then for all long reads available. In the training process, the prediction labels y ^ is treated by the softmax function, and cross-entropy E would be minimized to tune the parameters in RNNterm·mth. 2 E = y * - log ȳ + 1 - y * - log ( 1 - ȳ ) where y is the true labels of events regarding modifications, and ȳ = e ŷ n e ŷ m + e ŷ n e ŷ m e ŷ m + e ŷ n , and ŷ n is non-modification component for a prediction and ŷ m is modification component. (iv) A second deep neural network was also designed for considering high correlation of 5mCterm·dom modifications. Existing works 30 have demonstrated that 5mCterm·dom modification in human genome are highly correlated with nearby CpGterm·key sites and with the corresponding sites at the opposite strand. We also investigated the methylation percentage of a CpGterm·key site and of its nearby CpGterm·key sites and the corresponding sites at the opposite strand in both replicates of bisulfite sequencingterm·mth of NA12878, and found higher Pearson correlations between the methylation percentages of a CpGterm·key site and of its nearby CpGterm·key sites. To consider this effect, we designed the second deep neural network with four layers including an input layer, two hidden layers, and an output layer. The input of this layer is a 14-value vector: the predicted methylation percentage of a position, that of corresponding position of CpGterm·key site at the opposite strand, the number of nearby CpGterm·key sites within 25 bp, and the 11 discretized bins ([0, 0.05], [0.05, 0.15], [0.15, 0.25],…. [0.95, 1.0]) with value of how much percentage of nearby CpGterm·key sites has corresponding predicted methylation percentage. The first hidden layer has 100 hidden nodes, whereas the second hidden layer has 20 hidden nodes. Different layers were connected by full network with dropout (dropout rate = 0.7), and the output layer has sigmoid activation for outputting final methylation percentage for a CpGterm·key site. This network was only trained on the chromosome 1 of NA12878 with 100 epochs, and tested on all other chromosomes of NA12878 and HX1. (v) In a real application, given a set of FAST5 files, DeepModterm·mth needs a reference genome to generate prediction modification for events of interest in a long read, and a sequence of bases from events were aligned with a reference genome. After the alignment, a genomic position of interest in a reference genome was aligned with a set of events each from a long read and with a prediction modification label. The modification information of genomic positions of interest was summarized in a BED format in a strand-specific manner with single-base resolution, and contained total coverage, modification coverage and modification percentage: if a position of interest was covered by a long read with forward strand alignment, its forward coverage was increased by 1, and further if the aligned base in the long read was predicted to be modified, its forward modification coverage was increased by 1; similarly, if a position of interest was covered by a long read with reverse strand alignment, its reverse coverage was increased by 1, and further if the aligned base in the long read was predicted to be modified, its reverse modification coverage was increased by 1. Modification percentage of genomic positions of interest was calculated using its modification coverage divided by its coverage in total. In a real-world application, there is no filter for prediction sites, and reads with mapped quality < 10 would not be used (only tens of reads among thousands of passed reads had poor mapped quality in our experiments.).

Performance measurements For M completely modified and N completely un-modified bases (or motifs of interest), assume that DeepModterm·mth generated P modified predictions and Q un-modified predictions, accuracy, precisionterm·def, recall, APterm·mth, and AUCterm·mth are used to evaluate the performance. 3 accuracy = P ∩ M + Q ∩ N M + N 4 precisionterm·def = P ∩ M P 5 recall = P ∩ M M 6 F 1 - score = 2 * precison * recall precisionterm·def + recall where |*| is the size of *. Thus, accuracy is the fraction of correct predictions over all modified/un-modified cases in a data set, precisionterm·def is a percentage of correct predictions of completely modified cases over all modified predictions, while recall is the number of correct prediction of completely modified cases divided by the number of completely modified cases in a data set. F1−score is a tradeoff metric of precisionterm·def and recall. To evaluate binary classification of completely methylated positions and completely un-methylated positions in a reference genome, APterm·mth (which has been widely used in information retrieval) is weighted mean of precisionsterm·def achieved at each threshold of predicted methylation percentage. As predicted methylation percentage decreases from 1 to 0, recall usually increases, whereas precisionterm·def might decrease, and APterm·mth is in practice calculated by summing all precisionterm·def with the recall difference at two adjacent thresholds in an ordered list as weights by scikit-learn. APterm·mth could evaluate how a classifier performs for the predictions of modifications (The prediction of un-modifications is not fully considered in this measurement). The range of the five measurements above is from 0.0 to 1.0. AUCterm·mth is area under receiver operating characteristic curve, usually ranging from 0.5 to 1.0, and AUCterm·mth can evaluate how a classifier performs for all predictions (considering both modifications and un-modifications). For all the measurements, the larger the value, the better the classification is.

Training DeepModterm·mth We trained a 5mCterm·dom model of DeepModterm·mth on the positive control CG_MSssI and the negative control UMR of E. coli 4 . To validate the deep learning model and select optimal hyper-parameters, two independent validation strategies were used on CG_MSssI and UMR. One strategy is read-based, where all long reads were divided into two groups: one group was used to train DeepModterm·mth with 90% positive-control long reads and 90% negative-control long reads, and the other group was used for testing DeepModterm·mth with 10% positive-control long reads and 10% negative-control long reads. Since the data set has thousands of bases in each of thousands of long reads, each time, a small portion of training data of positive control and another small portion of training data of negative control were loaded to train DeepModterm·mth, and this process was repeated until negative control data were loaded multiple times. Then, DeepModterm·mth was tested on the testing group. In this strategy, long reads in testing groups might align with same reference sequences as some reads in training data. Thus, we also employed a second strategy, which is a region-based independent validation, where all long reads or bases of long reads mapped to the genomic positions from 1,000,000 to 2,000,000 of E. coli were used for testing, no matter whether the long reads were from positive control or negative control; the rest of long reads were used for training DeepModterm·mth. The loading process for training DeepModterm·mth was similar to read-based independent validation strategy did above.

As no prior knowledge can guarantee which hyper-parameters would be optimal, to build a well-trained model for modification prediction, different hyper-parameters were tested during this training and validation process, such as window size w from 7, to 11, to 15, to 21, to 31, and then to 51 with a step of 10. The detail performance of this validation was shown in Supplementary Table 1 . According to the detail performance, the final well-trained model was selected for making methylation prediction on other data sets for 5mCterm·dom prediction. A 6mAterm·dom model of DeepModterm·mth was trained in a similar way on Con1, Con2, tcgA_TaqI, gAtc_dam, and gaAttc_EcoRI.

Testing 5mCterm·dom model of DeepModterm·mth on an independent E. coli data set 5mCterm·dom and negative control Nanopore sequencingterm·key data sets from another research group 3 were used for cross-data set performance evaluation of DeepModterm·mth: the three 5mCterm·dom data sets are GCGC_Hhal, CG_MpeI, and CG_SssI, whereas the negative control data sets are Con1 and Con2 (Please refer to the Supplementary Notes for their detail description). FAST5 files for each testing group (5mCterm·dom data or negative control) together with a reference genome ( E. coli strand K-12 sub-strand MG1655 42 ) were taken as input of DeepModterm·mth, and each base of interest in long reads was predicted with a label to indicate whether it was methylated. Then, methylation prediction was summarized for each reference position of interest associated with a coverage, a methylation coverage and a methylation percentage. As in a positive control, a nucleotide of a specific motif were all methylated and negative control has no modification for all bases, we mixed each positive control with the negative control (Con1 and Con2), and calculated APterm·mth and AUCterm·mth to evaluate the performance of DeepModterm·mth.

Testing 5mCterm·dom model of DeepModterm·mth on human data sets of NA12878 and HX1 All long reads of NA12878 in FAST5 files together with a reference genome GRCh38 were used as input of DeepModterm·mth, and DeepModterm·mth made modification prediction for each base of interest in long reads with a label to indicate whether it was methylated. Then, predicted methylation was summarized for each genomic position of interest in GRCh38 associated with a total coverage, a methylation coverage and a methylation percentage. As native DNA sequencing in NA12878 contain naturally occurring modifications, and heterogeneity of sequenced DNA molecules generally exists, the criteria below was used to obtain high-quality methylation and non-methylation labels from bisulfite sequencingterm·mth 32 of NA12878 for bases in long reads: for a cytosine of a position in GRCh38 with both > 90% of methylations and coverage ≥ c ( c would be 1, 3, 5, or 10) in two replicates of bisulfite sequencingterm·mth, a cytosine in a long reads aligned with that position was considered to be modified and the cytosine at that position in the reference genome was called complete methylation; and if a cytosine of a position in GRCh38 has 0% methylations in both replicates of bisulfite sequencingterm·mth, a cytosine in a long read aligned with that position was considered to be un-modified and the cytosine at that position in the reference genome was called complete un-methylation. Bases in long reads aligned with all other positions were not considered, and other positions of interest in GRCh38 without complete (un)-methylations were not used for binary classification performance. Then, the DeepModterm·mth prediction percentage of methylation at reference positions of interest (for example, CpGterm·key sites for 5mCterm·dom) was compared against completely methylated or completely un-methylated reference positions, and APterm·mth and AUCterm·mth were calculated to evaluate binary classification performance of DeepModterm·mth. The testing of DeepModterm·mth on HX1 data had the similar process to what we did for NA12878.

Testing 6mAterm·dom model of DeepModterm·mth on C. reinhardtii To test 6mAterm·dom prediction on C. reinhardtii , we downloaded version 4.0 of C. reinhardtii 43 from JGI (Joint Genome Institute) as the reference genome, and input all FAST5 files of C. reinhardtii Nanopore sequencingterm·key data to DeepModterm·mth to obtain methylation prediction for each genomic position of interest in C. reinhardtii . The prediction results were compared with DA-6mA-seq results 34 , 35 for all GATC, TCGA. and GAATTCC motifs. In DA-6mA-seq, an adenine in motifs, which has an absolute score larger than 5 was considered as 6mAterm·dom modifications, whereas other adenine in motifs is considered to not be modified. mCallerterm·mth was run with a coverage ≥ 15 and motifs of interest.

Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Supplementary information Supplementary Information Peer Review File Reporting Summary