We first demonstrate the behaviour of both the existing methods and our proposed methods with toy examples of synthetic alignments. Subsequently, to demonstrate a major characteristic of scoring, i.e., the linearity, we merge multiple single-cell methylomes to test if the estimates of heterogeneity increase with the number of cells. In addition, we showed that comparing methylation heterogeneity can reveal differences between samples that might not be detectable by looking just at methylation levels. Finally, to test our methods on real data, we analyzed Arabidopsis methylome to profile non-CG methylation heterogeneity and human colorectal cancer data.

The behaviour of different scores in evaluating heterogeneity

We compiled a table of evaluation for several popular existing methods and our model-based methods (Table 1). The table includes specific features to be considered in methylation heterogeneity estimation and implementation. Overall, the main advantages of our model-based methods over existing methods lie in the possible extension to non-CG methylation sites, scoring linearity and the consideration of similarity between methylation patterns for unbiased and meaningful evaluation.

To assess the ability of our proposed method and other existing methods to precisely detect changes in methylation heterogeneity, we created toy examples with variable methylation patterns. To ensure a fair comparison between the methods, fully aligned reads are constructed to represent complete methylation patterns, and combinations of methylation patterns resembling alignments at a genomic region are simulated (see Fig. 2A–C, top panels). Firstly, we tested the hypothesis that the methylation heterogeneity would increase as new patterns occur. As a result, our three proposed models, as well as FDRP, qFDRP, MHL, ME, and EP showed such monotonic increasing trends (see dashed lines in Fig. 2A), while other methods showed differently. Next, we examined the importance of pattern similarity in the models, given that methylation patterns likely result from gradual changes associated with methylation maintenance. As shown in Fig. 2B and C, we would expect the ideal scores to increase when the patterns became more diverse (from left hand to right hand). We found that only PWS, MHL, and qFDRP were able to detect such changes in methylation patterns (see Fig. 2B and C). As a result, only two methods, PWS and qFDRP, aligned with both hypotheses. One specific concern for qFDRP is that the design of its score makes it easily influenced by the methylation level as described earlier (see Table 1). Our model-based method, specifically the PWS approach, demonstrated the ability to balance all these features effectively. Therefore, PWS was used for the following analyses.

Fig. 2

Evaluation of methylation heterogeneity methods. A–C Estimating methylation heterogeneity with synthetic datasets. Top panel lists combinations of methylation patterns at different loci. Circles are model-based methods, and triangles are existing methods. Dashed lines represent the methods with increasing trends. Four types of scores are used in the comparisons, model-based methods: AB, PWS and PHY; accordance-based methods: MC, PDR, FDRP; Entropy-based: ME, and probability-based: EP; Existing methods considering pattern similarity: qFDRP and MHL. D The methylation heterogeneity of merged mouse ESC and muscle single-cell methylome estimating by PWS and ME. E Genome-wide methylation heterogeneity ratios are plotted against different numbers of ESC single-cell methylomes. The black line represents the expected values given merged cells are all heterogeneous while the red represents linearity

Evaluating the scoring linearity using single-cell methylomes

When examining the real data, we expected the level of methylation heterogeneity to increase as new patterns are introduced. We processed a number of single-cell methylomes of mouse [10] from two different cell types, muscle cells and embryonic stem cells (ESC). We used the PWS method for the analysis as it is the only method that passed the previous evaluations in Figs. 2A–C. We also included ME and EP in the single-cell analysis as they share a similar data input format for such genome-wide analysis.

We would expect the methylation heterogeneity scores to increase from a pooled single-cell methylomes of one cell type to those of mixed two cell types. To this end, six single-cell methylomes of two cell types were pooled before the heterogeneity was estimated (i.e., the raw reads from multiple single-cell methylomes were pooled to form one merged methylome before they were aligned to the reference genome). We found that as expected methylation heterogeneity from PWS is increased in the mixed cell types (see Fig. 2D). In contrast, running the ME method on the same dataset, the mixed cell types showed a lower methylation heterogeneity score.

Next, we evaluated the methylation heterogeneity between different numbers of merged ESC single-cell methylomes for which we knew the compositions (see Fig. 2E). In a perfect setting, adding more cells of the same type would not increase heterogeneity. However, in this case of real data from ESC, each single cell methylome may not cover all expected patterns of ESC; based on the largely damaged DNA due to bisulfite treatments, the observed patterns may be very different between these single cell methylomes. Therefore, the methylation heterogeneity scores are likely to increase as the new patterns (from newly added single cells even from the same cell type) are added, we also expect to observe a gradual saturation of the methylation patterns, with the heterogeneity plateauing. First, different numbers of single-cell methylomes, i.e., 6, 8, 10, 12, 14, 16, and 18, were combined as merged methylomes to mimic bulk sequencing data. On those methylomes merged from many cells we would expect overall a higher methylation heterogeneity (indicative of cellular heterogeneity) than those from fewer cells.

We computed the genome-wide methylation heterogeneity ratio for each of the selected methods (for details of the procedure and calculations Additional file 1: Note S3). Overall, we observed a monotonic increase in the ratios as the number of single-cell methylomes increased with all methods. The extrapolated line (red) from 6 to 8 methylomes was drawn for PWS to demonstrate the expected linear increases per every 2 methylomes added. This line also revealed that ME and EP are likely to reach their plateaus quickly that are clearly deviated from being linear, suggesting that the two scores were less sensitive in detecting new patterns. The lower sensitivity revealed by such nonlinearity in real data application is less favoured (see Additional file 1: Fig. S1 for demonstration), particularly when different samples or regions were compared. Additionally, we found that none of these methods perfectly displayed the doubling property (see Fig. 2E black solid line). This could have occurred because in real data these single cells of ESC are typically not mutually exclusive groups. Still, we found the PWS heterogeneity is relatively linear compared to other methods. It also showed less deviation when the number of single-cell methylomes increased, making it a plausible scoring.

Comparing between methylation heterogeneity and methylation level

To determine the differences between methylation heterogeneity and the commonly used metric of methylation levels, we plotted the methylation heterogeneity estimated by the PWS method against the methylation levels of 3 replicates of samples from the human colorectal cancer (CRC) (Fig. 3A) and of Arabidopsis wild type methylome (Additional file 1: Fig. S2). As illustrated in Fig. 3A, the scatter plot indicated that the relationships between methylation heterogeneity and methylation varied across different cytosine contexts (i.e., CG, CHG and CHH, H = A, C or T). We observed a curve-shaped relationship between methylation heterogeneity and methylation level at regions of CpG methylation, that the regions with higher methylation heterogeneity have intermediate methylation levels that are found in both in human and in Arabidopsis. These regions are likely to reflect a dynamic process of epigenomic changes that are commonly observed in genic regions (Fig. 3C). In Arabidopsis we also profiled the methylation heterogeneity at non-CG sites. We found that the non-CH sites (i.e., CHG and CHH) showed very different relationship with methylation levels comparing to the CG sites (Additional file 1: Fig S2). While non-CG sites are low methylation, a fraction of them are highly methylated alone with a higher methylation heterogeneity. It is important to note that moderately methylated regions across all contexts retain a diverse range of heterogeneity, which could be easily overlooked when performing the evaluation using the methylation level alone. Furthermore, our PWS scores were able to detect the changes in methylation heterogeneity when the changes in methylation levels were not apparent. In brief, methylation heterogeneity could potentially complement the use of methylation levels for identifying minor changes that cannot be detected using methylation levels, and provides a different layer of biological information from the conventional methylation level.

Fig. 3

Genome-wide methylation heterogeneity profiles. A Mean methylation heterogeneity plotted with mean methylation levels in 3 replicates of adjacent normal samples of CRC. B Proportion of high (top 10%) and low (bottom 10%) methylation heterogeneity regions across different genomic features in Arabidopsis thaliana genome. C Metagene plot of A. thaliana CG methylation heterogeneity profile between highly and lowly expressed genes (top and bottom 25%). D Meta plots of A. thaliana CHG methylation heterogeneity between highly and lowly expressed TEs and their neighbouring regions

Methylation heterogeneity in plant at CG and non-CG sites

To reveal the genome-wide methylation heterogeneity in plants, we employed PWS to analyze an Arabidopsis wild-type methylome with a coverage of 58X [34]. We found that the regions of high methylation heterogeneity preferentially targeted transposable elements (TEs) in both CG and non-CG sites, which is different from regions with low methylation heterogeneity (Fig. 3B). In addition, the high methylation heterogeneity at CG sites are largely enriched at genebody comparing to non-CG sites; suggesting a differential preference between CG and non-CG. Subsequently we compare the methylation heterogeneity at genes and TEs of high and low expression (top and bottom 25%), see Fig. 3C and D. We observed a negative association between CG methylation heterogeneity and gene expression near transcription start sites (TSS), followed by a positive association toward the transcription end sites (TES); which indicates the dynamic epigenetic regulation of DNA methylation at promoter and genebody. We also found that lowly expressed TEs exhibited higher CHG methylation heterogeneity compared to highly expressed TEs (see Fig. 3D), suggesting that the methylation patterns at active TEs are highly variable in plant cells.

Our study produced the first map of methylation heterogeneity in the plant. High methylation heterogeneity regions were identified to be located at specific genomic features, which differed between CG and non-CG methylation heterogeneity. Methylation heterogeneity was demonstrated to be linked with transcriptional regulation. Our results illuminated the unique functions of CG and non-CG methylation heterogeneity in the Arabidopsis genome.

Strong association between genes with differential methylation heterogeneity and colorectal cancer-related diseases

Next, we wanted to demonstrate that the genomic regions with differential methylation heterogeneity may also be considered as biomarkers for phenotypes of interest. We downloaded and processed the human Reduced Representation Bisulfite Sequencing (RRBS) methylome data from CRC [35], which consisted of different stages, including stage III-IV CRC frozen tumours (tumour), normal-appearing mucosa as indicated by pathogens from the same patients (normal), and histologically confirmed matched normal samples collected from the margins on either side of the resected tumour (adjacent normal). The original study analyzed 10 samples per stage and found that the promoter methylation at specific cancer genes raised 40% to trigger the transcriptional changes at tumours, whereas at the adjacent normal the promoter methylation was only increased by 20% with no changes in expression, likely due to the lower changes in promoter methylation insufficient for triggering transcriptional changes.

As a demonstration of our method, we used 3 replicates from each normal, adjacent normal and tumour samples for methylation heterogeneity analyses using PWS method. The goal was to see if our PWS method is able to identifying putative biomarkers, as an alternative approach to the current approaches such as EWAS. A number of DNA methylation level studies have already shown that there existed methylation differences between say normal and normal-adjacent tissue [36] or between normal tissue and normal-tissue at risk of cancer development [37]. Therefore, we analyzed DNA methylation level in parallel, to assess the predictability between methylation heterogeneity and methylation level (Additional file 1: Fig. S4 for Venn diagram of differentially methylated regions, DMRs and differentially heterogeneous regions, DHRs).

In total 2,319 DHRs are identified between adjacent normal and normal (n = 911), and between tumour and normal samples (n = 1,558). These DHR are mostly found at genebody (Fig. 4A, left-panel). After normalising against RRBS genome we found the DHR are enriched at promoters, exons, 5’- and 3’ UTR but not from introns; suggesting a possible association with transcription.

Fig. 4

Comparison between methylation heterogeneity and methylation levels and the evaluation of PWS heterogeneity. A The proportion (left) and the enrichment plot (right) of DHRs at different genomic features. B Venn diagram of DHGs and DMGs. C Methylation heterogeneity heatmap of DHGs in the 3 stages of CRC samples. D IGV illustration of methylation heterogeneity estimated using PWS around CPXM2, with the DHR shaded in orange. Each blue bar indicates the mean methylation heterogeneity in bins of 400 bp, and the exact values of the bars at the DHR are labelled. E Composition of methylation patterns in all samples within a specific CG window 4 within the DHR from D for identifying a potential disease pattern

After we associated the DHRs with the genes that are co-localized with, 953 differentially heterogeneous genes (DHGs) are identified (Fig. 4B), whereas only 14 differentially methylated genes (DMGs) can be detected (15% methylation change and p-value < 0.05), including 2 genes C9orf69 and RAPGEFL reported in the original study (see Methods for identification of DHGs and DMGs). There is only one DMG, namely FK506-binding protein 10 (FKBP10), found to be also a DHG. This may suggest that the methylation level- and heterogeneity-based analyses actually targeted different sets of genes.

To track the changes of heterogeneity between stages, we plotted a heatmap of methylation heterogeneity using both tumour DHGs and adjacent normal DHGs (Fig. 4C). The heatmap shows there are clear changes of methylation heterogeneity from normal, adjacent normal to tumours, where most of the genes increased their heterogeneity towards tumours. A similar heatmap on DMGs was not able to reveal the differences accurately between the sample groups, as one normal sample is classified within the cancer group (Additional file 1: Fig. S5). We analyzed the enriched functions of the non-overlapping DHGs specific to either adjacent normal (Additional file 1: Fig. S6) or tumour (Additional file 1: Fig. S7) via ingenuity pathway analysis [38]. The enriched diseases and functions clearly indicated that DHGs identified by comparing adjacent normal samples against normal samples were involved in colorectal cancer-associated diseases; suggesting the changes of methylation heterogeneity at these genes are highly associated with the cancer progression, and the adjacent normal DHGs are predictive of CRC tumours. In summary, the DHG analysis complements conventional DMG approaches in the selection of regions associated with phenotypes of interest.

Identifying specific methylation patterns associating with increased heterogeneity

In total we identified 162 genes overlapping significantly between the tumour DHGs and adjacent normal DHGs (test for overlapping; p < 0.00001) (Fig. 4B; Additional file 1: Fig. S8). These genes showed strong and persistent changes of heterogeneity in DNA methylation towards tumour formation. The change of methylation heterogeneity at these genes may be suggesting specific methylation patterns emerging with the changes of cell types due to cancer formation or cell differentiation for example.

As an example, we found CPXM2 from the 162 overlapping DHGs. CPXM2 is a protein-coding gene that has been reported to be associated with several human disorders such as developmental diseases [39], Alzheimer’s disease and schizophrenia [40], and to promote tumour aggressiveness when active [41]. As shown in the screenshot of methylation heterogeneity at CPXM2 (Fig. 4D), an overlapping DHRs was constantly found at the promoter for comparisons between adjacent normal, and normal samples and between tumour and normal samples. The compositions of the methylation patterns at this particular DHR (Fig. 4E) revealed a specific methylation pattern labelled ‘1111’ in orange colour (fully methylated cytosines in a row) that seemed to be a “disease” pattern. It was not present in normal samples, but it started to appear in adjacent normal samples and became stabilized in tumour samples. Moreover, the proportions of reads showing this pattern increased in the presence of either an increased fully methylated ‘1111’ pattern or other partially methylated patterns that closely resembled ‘1111’ such as ‘0111’ or ‘1011’…etc.; instead of patterns resembling unmethylated ‘0000’, which was observed for most patterns in normal samples, the reads began to become similar to pattern ‘1111’. This verifies the ability of our model to detect changes in methylation patterns, which may serve as biomarkers for the early detection of disease.