Systematic benchmarking of tools for CpG methylation detection from nanopore sequencing

Abstract

DNA methylation plays a fundamental role in the control of gene expression and genome integrity. Although there are multiple tools that enable its detection from Nanopore sequencing, their accuracy remains largely unknown. Here, we present a systematic benchmarking of tools for the detection of CpG methylation from Nanopore sequencing using individual reads, control mixtures of methylated and unmethylated reads, and bisulfite sequencing. We found that tools have a tradeoff between false positives and false negatives and present a high dispersion with respect to the expected methylation frequency values. We described various strategies to improve the accuracy of these tools, including a consensus approach, METEORE ( https://github.com/comprna/METEORE ), based on the combination of the predictions from two or more tools that shows improved accuracy over individual tools. Snakemake pipelines are also provided for reproducibility and to enable the systematic application of our analyses to other datasets.

Conflict of interest statement

The authors declare no competing interests.

Figures

Fig. 1. Analysis pipeline for 5mC detection at CpG sites from nanopore sequencing.

The diagram describes the approach used to test the tools Nanopolish, Tombo, DeepSignal, Megalodon, Guppy, and DeepMod. Snakemake pipelines and command lines used are available at https://github.com/comprna/METEORE. A Snakemake pipeline was not developed for Megalodon and DeepMod as they can be run with a single command. Excluding DeepMod, all tools produce predictions per individual read and per CG site. In addition, all tools predict the methylation frequency at each genome site from fast5 input files. Methods that currently only accept single-read fast5 format are indicated.

Fig. 2. Accuracy analysis per CpG site on control mixture dataset 1.

a Violin plots showing the predicted methylation frequencies (y axis) for each control mixture set with a given proportion of methylated reads (x axis). The Pearson’s correlation (r), coefficient of determination (r2), and root mean square error (RMSE) are given for each tool. b For each method, we indicate the proportion of sites predicted outside a 10% window around the expected methylation proportion, i.e., each predicted site in the m% dataset was classified as “outside” if its predicted percentage methylation was outside the interval [(m − 5)%, (m + 5)%] for intermediate methylation values, or outside the intervals [0,5%] or [95%, 100%] for the fully unmethylated or fully methylated sets, respectively. The percentage is indicated on top of each bar, except for 100%. c Empirical cumulative distribution function (ECDF) plot showing the number of true negatives (y axis) for each tool according to different thresholds for the predicted methylation frequency below which a site was called unmethylated (x axis), using the dataset of 100 fully unmethylated sites. d ECDF plot showing the number of true positives (y axis) for each tool according to different thresholds for the predicted methylation frequency above which a site was called fully methylated (x axis), using the dataset of 100 fully methylated sites.

Fig. 3. Model accuracy at the individual read level and per-site accuracy analysis on control mixture dataset 2.

a Receiver operating characteristic (ROC) curves showing the false-positive rate (x axis) and true positive rate (y axis) for the predictions at individual read levels for the five methods tested, using reads from 0 and 100% methylated sets. b Precision–recall (PR) curves showing the recall (x axis) and precision (y axis) for the predictions at individual read levels for the five methods tested, using reads from 0 and 100% methylated sets. c ROC curves for METEORE for the random forest (RF) model (parameters: max_depth = 3 and n_estimator = 10) combining two methods, as well as combining the five methods. The curves were built from the average of a tenfold cross-validation with mixture dataset 1. Similar plots for another RF model using default parameters and a regression (REG) model are shown in Supplementary Fig. 3. d PR curves for the same models as in c. e Violin plots showing the predicted methylation frequencies (y axis) for each control mixture set with a given proportion of methylated reads (x axis) from the mixture dataset 2 for the five tested tools plus METEORE combining Megalodon and DeepSignal using a random forest (RF) or a regression (REG) model. The Pearson’s correlation (r) and coefficient of determination (r2) are given for each tool. f We indicate the proportion of sites predicted outside a window around the expected methylation proportion, i.e., a site in the m% dataset was “outside” if the predicted percentage methylation was outside the interval [(m − 5)%, (m + 5)%] for intermediate methylation sets, or outside the intervals [0,5%] or [95%,100%] for the fully unmethylated or fully methylated sets, respectively. The percentage is indicated on top of each bar, except for 100%.

Fig. 4. Comparison of CpG methylation predictions from nanopore with whole-genome bisulfite sequencing (WGBS).

a Distribution of methylation calls (n = 1743) from nanopore (y axis) across three WGBS methylation bins: no or low methylation (0.0–0.3) (n = 793), intermediate methylation (0.3–0.7) (n = 264), and high or full methylation (0.7–1.0) (n = 686). In the boxplots, the lower and upper boundaries of the box are the first and third quartiles of the data, respectively, with the median indicated by a thick black line. The lower and upper whiskers extend to 1.5 times the interquartile range. The outliers are represented by the black dots. b Pearson’s correlation (r) (y axis) between methylation frequencies calculated from nanopore reads and WGBS at sites by each of the tested tools (combining predictions from both strands) at each level of minimal input coverage, i.e., minimum number of nanopore reads considered per site as reported from the BAM file (x axis). c Mean reported coverage (using the coverage reported by each tool for all sites) at each value of minimum input coverage in b. METEORE (RF) is the combination of Megalodon and DeepSignal using a random forest (parameters: max_depth = 3 and n_estimator = 10), and METEORE (REG) is the combination of Megalodon and DeepSignal using a regression model.

Fig. 5. Comparison of CpG methylation predictions from nanopore with whole-genome bisulfite sequencing (WGBS) along Cas9-targeted regions.

Locally estimated scatterplot smoothing (LOESS) line plots of methylation calls frequency predictions (left y axes) from WGBS Illumina and from nanopore data using seven tools: Nanopolish, DeepSignal, Megalodon, Tombo, Guppy, DeepMod, and METEORE random forest (RF) and regression (REG) models. The plots include the nanopore read coverage (right y axes), shown as a light gray area. The panels below show the GC content of the region, using a window size of 50 bases. a–c Three of our ten target regions and d one of the eight target regions from Gilpatrick et al.. The depicted regions are a chr6:392228-401463, which covers the first and second introns of gene IRF4; b chr1:159199780-159212236, which covers the genes ACKR1 and CADM3;c chr2:1480363-1494141, which covers the gene TPO; and d chr3:49352525-49366169, which covers the genes GPX1. A zoom in of the LOESS line plots with the individual methylation calls and logos for the sequence context around the CpG site in four CpG islands (CGIs) of our target regions are shown in Supplementary Fig. 15.