DNA methylation analysis in plants: review of computational tools and future perspectives

Chemistry of bisulfite conversion and sequencing

Bisulfite sequencing is generally done in three main steps, namely (i) denaturing, (ii) bisulfite treatment and (iii) polymerase chain reaction (PCR) amplification. In bisulfite conversion, DNA is denatured in a process that separates the forward and reverse strands. This is followed by treatment with sodium bisulfite, which converts unmethylated cytosine into uracil—which is then converted to thymine during PCR [43]. Quantification of the abundance of each cytosine can be achieved via Sanger sequencing [44] or NGS technologies [45]. The DNA strands cease to be complementary after bisulfite conversion. Treatment of genomic DNA with sodium bisulfite [46] enables us to distinguish between highly similar (and yet different) methylated cytosine, which has the same base-pairing features as unmethylated cytosine. Mapping read sequences to a reference genome enables the determination of positions with matching and mismatching bases. This process enables identification of methylated and unmethylated bases.

Bisulfite sequencing can be accomplished with different sequencing kits depending on whether whole-genome bisulfite sequencing (WGBS) [18] or reduced-representation bisulfite sequencing (RRBS) [47, 48] is performed. Currently, WGBS remains the most informative method for generating DNA methylation data. It provides a huge wealth of data and requires no prior targeting. Unlike WGBS, which is expensive, RRBS can be performed more economically because it is restricted to CpG-enriched regions that make up a smaller portion of the genome. The restriction enzyme Msp1 cleaves at 5’-C*CGG-3’ targets (base preceding * is methylated), thereby, mainly CpG-rich regions are targeted—which is advantageous for large genomes.

Typical workflow for processing bisulfite sequencing data

Before reads are mapped to a reference genome, the sequencing quality of reads can be checked with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) or NGS QC Toolkit [49], followed by removing low-quality bases and adapters with, among others, Trim Galore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore), cutadapt [50] or Trimmomatic [51]. However, some WGBS data processing tools integrate various analytic steps—enabling data preprocessing, read alignment, a more robust statistical analysis that output statistics such as read coverage, the percentage of uniquely aligned reads and statistics on the three methylation contexts (CpG/CHG/CHH). One such tool is gemBS [52], which is a recently published pipeline for processing and analysis of WGBS data. The pipeline integrates data preprocessing and analysis steps from adaptor trimming through downstream statistical analysis of mapping results. gemBS uses the high-performance read aligner GEM3 [53] as a dependency and BScall (embedded in samtools, bcftools; http://samtools.sourceforge.net), which is a variant caller for bisulfite sequencing data. Both GEM3 and BScall support single and paired-end reads. Further reading on the generic workflow of analyzing WGBS is found in the work of Liang et al. [54] and Wrecyzcka et al. [55].

Non-bisulfite-based methods and related bioinformatics tools

While bisulfite sequencing methods represent the most popular approaches for analyzing epigenomic data, there are other approaches within the field of DNA modification-based methods. These approaches include methylated DNA immunoprecipitation (MeDIP)-seq and MethylCap-seq (a robust procedure for genome-wide profiling of DNA methylation) in MeDIP analyses [56] where the genomic DNA is randomly sheared, sonicated and immunoprecipitated with an antibody recognizing 5-methylcytidine. Precipitated DNA can either be sequenced or hybridized to microarrays. MethylCap-seq uses the methyl-CpG-binding domain of MeCP2 [57] while Oxidative bisulfite sequencing (oxBS) [58] is used to specifically detect 5-methylcysteine and 5-hydroxymethylcytosine (5hmC) that can be also done with ‘Tet’-assisted bisulfite sequencing [59]. CAB and fCAB are used for the recognition of 5caC [60]. Notably, the presence/absence of 5hmC in plants remains contentious. Some scholars claim that 5hmC is present in plants [61, 62] while others claim it’s absent [63]. A comprehensive overview of the various tools is given at https://omictools.com/medip-seq-category.

Tools for analyzing epigenomics datasets

Bismark [64] and BSMap [65], as one of the 1st published tools for quantifying epigenomic datasets, had to address the challenge of attaining high-read mapping efficiency to enable a sensitive sequence search. Bowtie [66], Merman [67], SNAP (http://snap.cs.berkeley.edu) and Bowtie2 [68] have been used as dependencies in epigenomics tools, for instance, BS-Seeker [69], BS-Seeker2 [70], BS-Seeker3 [71], BRAT-nova [72], WALT [73] and Bismark, which are currently among the most commonly applied tools for mapping bisulfite methylation data. We outlined the most common tools for mapping bisulfite sequencing data along with tools that allow for the detection and analysis of differentially methylated regions (DMRs). The program parameters as well as input and output data formats are specified in Table S1. This table provides an overview of the main tools for mapping and analysis of epigenomic data—particularly for bisulfite sequencing data. Additionally, we also categorized the tools into three major classes, namely (i) mapping, (ii) statistical analysis and (iii) complete pipelines (Table S1). The defining features for each tool, such as their ability to handle single or double-stranded sequence data as well as their ability to process data and perform downstream statistical analysis, are also provided. Reviews by Adusulalli et al. [74], Shafi et al. [75] and Wrecyzcka et al. [55] complement our overview Table S1. The most frequently applied computational epigenetics methods were applied and tested using DNA methylation data, particularly with data acquired from bisulfite sequencing experiments. Therefore, there are many statistical procedures available for analyzing methylome data—categorized into the parametric and non-parametric approach. Both approaches are widely used in the literature [76]. For instance, MethylMix [77] is an excellent example of a parametric approach that uses Bayesian mixture models to identify DNA methylation states of genes as either hypo- or hypermethylated. The method entails fitting a distribution function onto the frequencies of DNA methylation counts. The advantage of using non-parametric models is that no prior knowledge of the data distribution is required. However, when such knowledge is available, then parametric models are the preferred choice for modeling such data. MethylMix quantifies the effect of DNA methylation on genes, which is interesting for integrative studies that aim at establishing the association between the methylation states of the individual genes and their expression profiles. Investigating such associations unravels any hidden variations within and between samples (or tissues) as illustrated in [78–80]. Lea et al. [81] discussed the applications of mixed models on DNA methylation in plant epigenetics. They specifically focused on the binomial mixed model with the sampling-based algorithm (MACAU, mixed model association for count data via data augmentation) for the approximation of parameters and computation of P-values.Other modeling frameworks are based on algorithms that integrate various analytical steps resulting in the detection of DMRs across the entire genome, for instance, (i) the weighted optimization algorithm proposed in [82] (which is an extension of MethylKit [83]) and (ii) ChAMP.DMR [84], which applies the Bumphunter [85] or ProbeLasso Algorithm [86]. An example of a non-parametric model is the Bayesian approach based on the Dirichlet process beta mixture model—which is used for clustering methylation profiles [76]. The model considers the DNA methylation expressions consisting of an infinite number of beta mixturex distributions [87, 88].

DNA methylation: plant physiology and pathophysiology

Investigating the dynamics of DNA methylation in plant growth and development requires the analysis of samples from different plant tissues (e.g. [34]). To our knowledge, no existing software has been developed specifically for the analysis of plant physiology and pathophysiology. However, there are many studies analyzing bisulfite data using samples from different plant developmental stages (from seedlings to mature plants). For instance, Bismark—in leaf tissues from bread wheat seedlings [89], BSMap—for various datasets from different tissues in A. thaliana [90] and BS-Seeker2—for young Zea mays leaves [91]. With rapid advancements in the development of software/tools for analysis of epigenomes, we are optimistic such tools will soon be available to the public.

DMRs and their significance

Genomic regions (or bases) with different methylation profiles between samples are known as DMRs. This is also referred to as differentially methylated CpG sites since the CpG-methylated sites occur in much larger numbers compared to the non-CpG contexts (CHG and CHH) [92, 93]. Peak detection enables the identification of CpG islands—which are essential for differentiating methylation profiles between samples (typically between controls and test samples). CpG islands are not randomly distributed in the genome but are instead grouped close together [94]. Long stretches of non-dense CpG sites, known as CpG shores, can also be detected. Combining the methylation profiles of both CpG-islands and CpG-shores enables more efficient comparative analysis of DNA methylation profiles between samples.

Various statistical algorithms have been proposed for identifying DMRs—the most popular ones being methylKit [83], metilene [95], DMRcaller [96] and Bumphunter [85]. For elaborate discussions on the DMR detection methods and a discussion on choosing the right method for DMR detection, see [97, 98]. The tools are written and compiled in different programming languages (e.g. R, Python, Perl, Java, C and C++; Table S1). Essentially, such tools are used to identify DMRs from either targeted regions of the genome or from the whole genome. Critical considerations have to be made, e.g. the choice of experimental designs for experiments and statistical methods for data analysis [99]. DMRs are intricately linked to transcription and the abundance of CpG sites (CpG islands). A high concentration of CpG sites is often found within the promoter regions of genes—so it is essential to accurately identify such sites. Methylation of promoter regions influences the level of transcription—heavy methylation disrupts transcription, and de-methylation leads to transcription reactivation [100–102].

Peak identification and normalization are crucial initial steps in analyzing DNA methylation data and visualization and can be useful for comparing datasets and judging the performance and agreement between tools. Post-processing and visualization of (differentially) methylated sites enable high-resolution exploration and comparison of regions in the genome for variations in methylation profiles. Therefore, tools like BiQ [103] and BSeQC [104] have aided quality control and visualization of methylation data, thereby enabling researchers to explore data attributes and perform data quality control before analysis. There are many methods for clustering methylation marks such as the dynamic genome warping [105] approach that uses hierarchical clustering and the combination of different epigenomics analytic platforms and data integrative modules. Dynamic genome warping has been demonstrated to be a reliable way to get more meaningful results from datasets (for instance, [106]). To utilize this method, Liang et al. [54] developed a webserver to analyze WGBS data and their platform includes major steps for detection and mapping of the conversion rate, detection of DMRs and their association with gene expression. Wreczycka et al. [55] discussed data requirements and computational attributes for specific software and assess bisulfite sequencing data analysis methods, alignment and data processing, detection of differential methylation and assess strategies for handling large epigenetic datasets. In contrast, our work highlights existing asymmetries between mapping tools and contrasts their computational capabilities.

Another important aspect in plant epigenetics is how hypomethylation and hypermethylation affects gene expression. The concept of hypomethylation and hypermethylation is not limited to plants as they have also been extensively studied in cancer progression in humans [107], coronary heart disease [108] and eukaryotes in general [109]. The division of DMRs into hypo- and hypermethylated enables investigations into the influence of both types of methylation on gene expression. Many computational tools have integrated modules that enable the extraction and quantification of the extent of hypo- and hypermethylation in genes. One such tool is MethylMix, which requires that changes in a gene’s methylation state must also agree with its expression profile. Additionally, it requires a treatment and control sample (for agricultural studies) or healthy and disease conditions (for clinical studies).

Downstream analyses of bisulfite methylome data

After data processing and calling of methylation sites, down-stream analysis can be performed—including the functional annotation of DMRs and analysis of the associated pathways influenced by the targeted genes. Such analysis enables the assignment of functions and gene annotation as seen in the overviews of Bioinformatics omicX tools (https://omictools.com/epigenomics-category). Examples of tools for performing down-stream analysis are given in Table 1.

Technical challenges: conversion rate, repetitive regions and DMRs

The main challenges in the analysis of DNA methylation data include incomplete methylation patterns and overdispersion of read mappings [110–112]. Here, overdispersion means the presence of variability in the reads compared to the expected read distributions based on a given model structure. When epigenomics marks coincide with repetitive regions in the genome, mapping tools need to keep reads that map to multiple genomic locations—making these tools slower and computationally memory-intensive. This problem can be partly circumvented through parallel computing using multiple threads, especially for larger repetitive plant genomes.

Conversion rates

As a method for studying DNA methylation, bisulfite conversion involves the conversion of cytosine to uracil (while 5-methylcytosine, 5-mC remains unchanged). Bisulfite sequence conversion rates vary for different datasets. It is essential for conversion rates to be determined accurately to ensure the reliability of downstream data analysis. Reliable results can be obtained from datasets with bisulfite conversion rates higher than ~0.999 (see, e.g. [113]—demonstrated using their tool MethQA). However, they urge caution for datasets with lower conversion rates. Modern commercially available bisulfite sequence conversion kits generally indicate conversion efficiencies of 90–100% [114]. An elaborate discussion on methods for estimating conversion rate from bisulfite DNA methylation data is provided in [115, 116].

Description of experiment: benchmarking selected tools

We aimed to determine how the well-established computational epigenomics methods perform on a small genome such as A. thaliana with ~130 Mbp (TAIR10) compared to a genome with a high repeat content and much larger genome size such as bread wheat—taking chromosome 1A (Chr1A) for demonstration purpose [117]. We used bisulfite sequencing data from two studies (with accession numbers SRR429549 [118, 119] for A. thaliana and ERR1141918 [89] for Triticum aestivum; data from NCBI) and applied four methods: BSMap [65], Bismark [64], BS-Seeker3 and segemehl [120]. Our analysis focused on the speed and agreement of common methylated sites between the tools. BS-Seeker3 was the fastest, followed by BSMap, while Bismark and segemehl were the slowest irrespective of genome size—especially for multiple threads (Figure 1A and B). When using a single thread, segemehl (keeping reads that mapped a maximum of three times) performed slowest compared with the other methods. Overall, the computation time required for the T. aestivum (Chr1A) dataset is significantly longer than those from A. thaliana (Figure 1A and B). When comparing the reported sites, we found that, for A. thaliana, 562 051 sites are shared among all four tools. While most sites were overlapping between BSMap, BS-Seeker3 and Bismark, likely because they use the same mapping software, segemehl reported only ~10% of these sites. However, for T. aestivum, ~101 944 sites were reported with most of them being reported in segemehl (Figure 1C and D). The existence of such asymmetries requires more attention and is certainly worth taking into consideration when using the different computational tools. Other studies on comparisons of the performance of epigenetics analysis tools, specifically focusing on mapping short reads for bisulfite sequencing data, can be found in the work of Tran et al. [121]. Several studies have also compared runtime and memory consumption of different epigenomics tools, such as Tran et al. [121] who compared the five bisulfite short read-mapping tools BSMap, Bismark, BS-Seeker, Bisulfite Sequencing Scorer (BiSS) and BRAT-BW and Bismark performed best on real data, followed by BiSS, BSMap and BRAT-BW and BS-Seeker. Recently, Huang et al. [71] proposed BS-Seeker3—a fast mapping tool for bisulfite data and compared it performance for runtime and sensitivity to sister tools like Bismark, BRAT-nova and BSMap. Additional to being accurate and versatile, Huang et al. concluded that BS-Seeker3 is an ultra-fast pipeline to process bisulfite-converted reads. The tool also aids visualization of methylation data, hence justifying its comparability to the other three tools (Bismark, BRAT-nova and BSMap).

We simulated reads from A. thaliana and bread wheat using the tool by Sherman (https://www.bioinformatics.babraham.ac.uk/projects/sherman) to test the performances of the four tools by comparing the precision and sensitivity along all chromosomes (Figure 2). The sensitivity, also sometimes referred to as recall, is defined as TP/(TP + FN). The precision is defined as TP/(TP + FP), where TP—true positive, FN—false negative and FP—false positive. We observed best performances for the Bismark, followed by BSMap and segemehl, while BS-Seeker3 seemed to have a lower sensitivity in A. thaliana compared to the other tools. For bread wheat a similar order to performances of tools was observed when reads where simulated for each subgenomes of chromosome 1 with the three genome copies. All scripts were provided in GitHub (https://github.com/jomony/EPItools/blob/master/README.md).

Feature comparison between the tools and related literature benchmarking

To further benchmark the performance of the tools, we used bisulfite sequencing data from five plant genomes. These genomes consist of the dicots: A. thaliana (genome size, ~0.13 Gb; SRR4295494), Arabidopsis lyrata (~0.21 Gb; SRR3880297) and Glycine max (~1.2 Gb, SRR5079790) and also the monocots: T. aestivum (chromosome 1A; size, ~0.67 Gb; ERR1141918) and Oryza sativa (~0.43Gb; SRR7265433). Figure 3A shows the results of a comparative analysis of the memory footprint analysis of the performance of the four tools benchmarked using data from five genomes. These results come from mapping the bisulfite reads data to their respective reference genomes. Association analysis was performed for each of the four tools as seen in the linear regression model fits (Figure 3B–E). The results show that the genome sizes for each of the five genomes are significantly correlated to the memory footprint analysis (P < 0.05).

The key attributes and parameters for the four tools are summarized in Table S2. This table presents a summary of the supported features in the four tools (BSMap, BS-Seeker3, Bismark and segemehl). Such features are essential for deciding on which tool to use for mapping reads and data analysis. Examples of such features can also be found in the work of Guo et al. [70] and Tran et al. [121]. Lee et al. [122] evaluated the mapping accuracy and mapping rates for Bismark, BSMap and BS-Seeker2 as a function of the error rates. Using WGBS data, they assessed the influence of the error rates on the mapping rates and mapping accuracy and observed that at low error rates (<4%), BSMap had a higher mapping rate than Bismark and BS-Seeker2. On the contrary, BSMap had a lower mapping accuracy than Bismark and BS-Seeker2. They also showed that mapping accuracy is independent of the methylation level.

A discussion on benchmarking approaches with a focus on short sequence mapping tools is found in the work of Hatem et al. [123]. They assess the performance of various aligners for the read mapping tools and benchmark them using criteria such as mapping percentage, running time and memory footprint. Variations in parameters such as seed length, base quality and single- or paired-end reads on the mapping quality are also evaluated. Benchmarking of tools by comparing the performance of each tool based on multiple attributes can be achieved in various ways, for instance, by assessing (i) the effect of the read length and sequencing error, (ii) the effect of data processing and (iii) the effect of varying parameters in the tools. These are some of the approaches discussed by Tran et al. [121]. They compared the performance of epigenomic mapping tools such as BSMap, Bismark, BS-Seeker, BRAT-BW [124] and the BiSS [125]. Tran et al. primarily benchmarked the performance of the tools basing on mapping efficiency (as the percentage of reads that map uniquely to the genome) and the central processing unit (CPU) time.

Outlook

In the near future, there is a need for more comparative analyses to explore the epigenomes of diverse plants in different development stages together with various stress factors. This would enable the discovery of exclusive and common epigenetic regulatory mechanisms. Uncovering and exploiting such mechanisms could potentially promote adaptation to changing environmental conditions. Moreover, a large number of methylomes are required to study the effect of the environment and stress conditions on the epigenomic state of a single plant [126, 127]. Resources like the 1001 Epigenomes Project (https://1001genomes.org) in A. thaliana are exciting datasets to aid in our understanding of the role of the epigenome. However, it remains unclear whether the observations in these studies are directly applicable to crops.

Computational tools are instrumental for bridging the gap between mapping of sequenced reads, the accurate prediction of methylated sites and their statistical analysis However, this effort is hampered by variations in the size of epigenomic marks and the complexity associated with normalizing peaks. The need to increase crop yield on the same amount, and in some cases dwindling, of arable land is another important aspect that requires advancements in epigenomics studies. Several studies have shown that during seed and grain development, the plant epigenome changes and leads to gene silencing. Therefore, a change in the epigenetic state of a plant would result in an increase in its likelihood of adapting from one geographical location to another or to different environmental conditions.

Lämke and Bäur [128] argued that such modifications have the potential to provide a mechanistic basis for stress memory in plants. This enables plants to respond more efficiently to recurring stress from the environment, for instance, drought and salinity stress [129], a topic that was reviewed by Golldack et al. [130] (and more recently by Yang and Guo [131] and Abhinandan et al. [132]). This might enable plants to prepare their offspring for future attacks from stressors and to improve their adaptation to specific stress factors [130]. Plant adaptation to stress might enable us to explore new ways to improve yield, for instance, by shortening or prolonging the time for grain development, by finding ways to regulate the expression of the three homeologs in wheat or by interfering with fruit ripening (as seen in tomatoes [133–135] and other fruits like peach, apples and strawberries [136]). A more intriguing discussion on the epigenetic mechanisms of plant stress response and adaptation to different environmental conditions was reviewed in [137–139].

In this review, we have discussed the use of bioinformatics tools to study DNA methylation data in plants. Notably, several studies in humans and mouse were successfully performed using popular tools like BSMap, BS-Seeker/BS-Seeker2/BS-Seeker3, Bismark in mouse and segemehl in human cancer cell lines. For the analysis of bisulfite sequence data, most of the fundamentals of the chemical background and methylation principles are the same; however, the major difference between the use of such tools in plants and animals (specifically, in humans and mouse) is the genome structure organization and the presence of predominantly more CHG/CHH methylation contexts in plants. The most predominant context of DNA methylation in mammals is the symmetric CG—estimated to be at ~70−80% of CG dinucleotides genome-wide [140]. The mechanisms of regulation and function of DNA methylation vary in animals and plants [141, 142]. These variations in regulation and function mechanism, coupled with genome structure differences and complexity levels, is a motivating factor for integrating small subtle differences in mapping and analysis tools for epigenome data. Another important difference of plants and animals is how they are able to demethylate their genome. So far, enzymes removing directly the methyl group from cytosines have not been identified in plants, but they are important components of mammalian DNA methylation homeostasis. Plants use either passive mechanisms (not maintaining methylation during DNA replication) or base excision and subsequent repair for direct removal of methylated cytosines. Unlike with the human genome, the CHG/CHH contexts that are more abundant in plants [143] need to be integrated into the mapping and analysis of methylome data. Many plants have large and repetitive genomes compared to that of humans. Such large genomes are a limiting factor in the analysis since they require a lot of computational resources. The sequence mapping to references and statistical computational time for large genomes such that of bread wheat (~17 Gb) and barley (~5.3 Gb) is likely to scale linearly.

Concluding remarks

In the last decade, there has been tremendous progress in the development of tools for analyzing epigenomic data; however, numerous challenges remain. For instance, the visualization capacity of many tools remains either inadequate or lacks essential modules for handling and displaying statistical outcomes from the resulting analysis. Additionally, the ability of these tools to scale-up and to handle large genomes remains an issue for further exploration. Technically, most computational tools for analyzing epigenomic data perform well for datasets from organisms with a genome size that is smaller than the human genome (~3 Gb). For much larger and complex genomes, more computational resources are required, and the genome structure (whether diploid, hexaploidy or tetraploid) and repetitive nature of the genome have to be taken into consideration during mapping to a reference genome. This is demonstrated in our example where we compared the mapping efficiency for Arabidopsis and a wheat chromosome; however, the complexity in genome structure, the presence of TEs and the lack of consistent gene annotations for some plants remain a major obstacle to advancing epigenetic research.

In the next decade, there is likely to be a steady improvement in sequencing methods and performance of already existing computational algorithms. Recently, it was shown that even well-established sequencing methods might be prone to errors, leading to misleading results, e.g. DNA immunoprecipitation sequencing [144]. Discovering and amending such errors can lead to new findings from the previous studies and limit these errors’ damage to future studies. This will aid further epigenetic research not only in plants but also in life sciences in general. Additionally, a few tools have the capability to effectively get more information out of low-coverage data. Developing new tools or improving on existing ones to attain optimal results using low coverage data and fewer replicates would save experiment and sequencing costs. A high sequence coverage allows for good data quality and enables robust statistical analysis [145]. Achieving high sequence coverage can be quite expensive and the minimum desired coverage can depend on the research objectives at hand. Typically, a coverage value of 5–10 × is sufficient for many comparative studies and for achieving reliable methylation calls [145]. However, studies have demonstrated that coverage values as low as 2× is sufficient [146]. Accurate identification of DMRs in large samples, especially between multiple conditions, remains a challenge—despite tremendous progress already made in this area.