Predicting CpG methylation levels by integrating Infinium HumanMethylation450 BeadChip array data

    Predicting CpG methylation levels by integrating Infinium HumanMethylation450 BeadChip array data Shicai Fan, Kang Huang, Rizi Ai, Mengchi Wang, Wei Wang PII: DOI: Reference:

S0888-7543(16)30013-1 doi: 10.1016/j.ygeno.2016.02.005 YGENO 8802

To appear in:

Genomics

Received date: Revised date: Accepted date:

21 December 2015 19 February 2016 22 February 2016

Please cite this article as: Shicai Fan, Kang Huang, Rizi Ai, Mengchi Wang, Wei Wang, Predicting CpG methylation levels by integrating Infinium HumanMethylation450 BeadChip array data, Genomics (2016), doi: 10.1016/j.ygeno.2016.02.005

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Predicting CpG Methylation Levels by Integrating Infinium HumanMethylation450 BeadChip Array Data

SC R

IP

T

Shicai Fan1,2,*, Kang Huang1, Rizi Ai2, Mengchi Wang2, Wei Wang2,* 1 School of Automation Engineering, University of Electronic Science and Technology of China, Chengdu, China 2 Department of Chemistry and Biochemistry, University of California San Diego, La Jolla, CA, USA *Correspondence to: Wei Wang ([email protected]), Shicai Fan ([email protected])

CE P

TE

D

MA

NU

Abstract: The Infinium HumanMethylation450 BeadChip array, referred as 450K array hereinafter, has been widely adopted as an affordable technique to determine DNA methylation. Tens of thousands of data have been generated on diverse cell types and patient tissues, which has provided great insight into understanding the crucial roles of epigenetic modifications in many biological processes and diseases. The limitation of this technique is its coverage, which measures methylation levels of about 450,000 CpGs, accounting about 1.6% of all CpGs in the human genome. In the present study we developed and compared computational models to significantly expand the coverage of Illumina 450K (~11 folds). Using the whole genome bisulfite sequencing and Illumina 450K data in the human H1 embryonic stem cell, we showed that the predicted and measured methylation levels were well correlated. Our proposed model showed superior prediction accuracies compared to the existing methods on the same dataset. When applied to predict the DNA methylome on other cells, our proposed model achieved comparable performance in cross-validations, which indicates the generalizibility of the method. Our method would thus be invaluable to maximize the usage of the existing data. Keywords: DNA methylation, CpG loci, Prediction, 450K array data

AC

1. Introduction A key epigenetic modification in vertebrates is DNA methylation that plays crucial roles in many biological activities, such as X-chromosome inactivation, the establishment and maintenance of tissue-specific gene expression profiles, chromosomal domain silencing, stem cell differentiation, embryonic development and inflammation(Barlow, 2011; Bird, 2002; Heard, 2004; Jones, 2012; Laird, 2010). In human somatic cells, approximately 70–80% of CpG dinucleotides undergo methylation (Ehrlich et al., 1982). Aberrant methylation and unmethylation are related to various diseases(Baylin & Jones, 2011; Laird, 2010). Deciphering the role of DNA methylation requires accurate assessment of the genome-wide DNA methylation map of different human tissues, and various experimental detection methods have been developed. The most useful and popular approaches are bisulfite sequencing, bisulfite microarrays and enrichment-based methods (Brinkman et al., 2010; Down et al., 2008; Gu et al., 2010; Sandoval et al., 2011). Bisulfite sequencing has the most comprehensive genomic coverage, but the high sequencing depth makes it expensive. Although the costs of bisulfite microarrays and enrichment-based methods are relatively low, array-based platforms are constrained to priori-targeted regions, and enrichment-based methods have relatively

ACCEPTED MANUSCRIPT low resolution and high susceptibility to experimental biases(Bock, 2012).

SC R

IP

T

Computational methods to predict DNA methylation have been developed, including the methods predicting methylation status of CpG loci (Bhasin et al., 2005; Stevens et al., ; Yang et al., 2012) and CpG islands (Bock et al., 2006; Fan et al., 2008; Fang et al., 2006; Feltus et al., 2003; Feng et al., 2014; Previti et al., 2009; Zheng et al., 2013). DNA sequence features are the most commonly used in these models; inclusion of epigenetic features such as histone modifications could enhance the prediction performance(Bock et al., 2007; Ernst & Kellis, 2015; Fan et al., 2008). Furthermore, the whole genome methylome was predicted from DNA sequence features and methylation values measured by DNA immunoprecipitation (MeDIP-seq) and methylation-sensitive restriction enzyme (MRE-seq) sequencing (Stevens et al., 2014).

TE

D

MA

NU

A gap in predicting DNA methylome is to use Illumina HumanMethylation450 BeadChip array data (hereinafter referred as 450K array). The 450K array has been widely used in DNA methylation detection for human tissues, particularly in thousands of patient samples. Extension of the coverage of the 450K data would thus be invaluable to maximize the usage of the existing data. We presented here the first attempt to tackle this task. Based on the data of human H1 embryonic stem cells (ESC), we trained and compared three prediction models (logistic regression, random forest and support vector regression) by integrating the 450K array data and sequence features. We showed our models achieved satisfactory performance that is superior to an existing method EpiGRAPH(Bock et al., 2009). Furthermore, we demonstrated the generality of the trained model in predicting methylome in new cell types.

AC

CE P

2. Materials and Methods 2.1. DNA methylation data For model training, we downloaded the Whole Genome Bisulfite Sequencing (WGBS) data (accession number GSM432685) of H1 ESC from the GEO database as the gold standard for DNA methylation levels (dependent variable). The methylation level of each CpG locus was represented as a methylation ratio. The corresponding 450K array data of H1 ESC were also retrieved from the GEO database (accession number GSM853420) . The methylation value of each CpG locus was represented as a beta-value ranging from 0 (nonmethylated) to 1 (complete methylation). For the model test on independent data, the WGBS and 450K array data of the H9 ESC were also downloaded from the GEO database with accession numbers of GSM706059 and GSM853421, respectively. Both the WGBS and 450K array test data were quantile normalized to the corresponding WGBS and 450K training data, respectively. 2.2 Prediction features The key point of our model is to integrate the existing 450K array data with sequence features as the prediction features (independent variables). For each CpG locus, we collected the 450K array data of its 1000 bp flanking regions. The average methylation value of the 450K array data was regarded as one feature. Furthermore, we considered the occurring frequency of 1-mer and 2-mer DNA base pairs (20 features), NpN ratios (N=A/C/G/T) (16 features) in the 500bp flanking regions of a CpG. In total, 37 features were used as the input features for the model construction.

ACCEPTED MANUSCRIPT

T

2.3 Training and prediction In the model construction for the H1 ESC, we trained a model for each of the 24 chromosomes separately. For each of the 24 models, 3-fold cross-validation was conducted for 20 times and the averages were used to assess the model performance.

SC R

IP

In the field of bioinformatics, logistic regression (LR), random forest (RF) and support vector regression (SVR) are three widely used machine learning algorithms, and their performance varies on different datasets. In this study, we compared their performance on predicting DNA methylation levels based on 450K array and sequence features. The three models were constructed based on the glm function, the randomForest package and the e1071 package in R.

MA

NU

2.4 Performance metrics The performance of the models was evaluated using correlation coefficients (Pearson and Spearman), concordance, specificity (SP), sensitivity (SE), accuracy (ACC) and Matthews correlation coefficient (MCC). E[( X   X )(Y  Y )] Pearson correlation coefficient was calculated as  pearson  ,

 XY

CE P

TE

D

where  X and Y are the means of X and Y ,  X and  Y are the standard deviations of X and Y , and X and Y represent the fitted methylation value and the real WGBS data, respectively. Spearman’s correlation coefficient was calculated 6 di 2 as:  spearman  1  , where n is the sample size, and the n raw methylation n(n 2  1) values X i and Yi are converted to ranks xi and yi , with di  xi  yi .

AC

For the concordance measurement, we used the same definition as in the Stevens et al. study(Stevens et al., 2014), which is the percent of CpGs with a predicted methylation value difference of less than 0.25 with the real WGBS value. For SE, SP, ACC and MCC, we defined the methylation status as +1 when the methylation value was larger than 0.5 and as -1 when the methylation value was not larger than 0.5. The four metrics are defined as below: TN SP  TN  FP TP SE  TP  FN TP  TN ACC  TP  FN  TN  FP TP  TN  FP  FN MCC  (TP  FN )  (TP  FP)  (TN  FP)  (TN  FN ) where TN, TP, FN and FP represent the number of true-negatives, true-positives, false-negatives and false-positives, respectively. 3. Results 3.1. Correlation-related measurements We observed the correlations between WGBS data and all the 37 features, and found

ACCEPTED MANUSCRIPT

T

the 450K array feature was with the highest correlation coefficient (0.7602), indicating that 450K array data would be the most correlated feature. By using the 3-fold cross-validation test 20 times, we obtained the average performance on each chromosome. The average correlations and concordances with LR, RF and SVR are shown in Figure 1-3.

SC R

IP

The performances of all three algorithms are comparable. The Pearson correlation coefficients between the predicted and measured methylation values are all about 0.85 (the Spearman correlation coefficients are approximately 0.7), and the concordance measurements are mostly larger than 80%. Overall, the predicted methylation values are highly correlated with the real WGBS data.

MA

NU

The prediction results showed that the RF model performs slightly better. The average Pearson correlation coefficient of the RF model is approximately 0.02 larger than LR and 0.04 larger than SVR. In the RF model, the prediction results on chromosome 18 are the best: the correlation coefficient reaches 0.91(0.82), and the concordance reaches 88%; the results on chromosome 19 are the worst: the correlation coefficient is 0.84(0.73), and the concordance is 82%. 3.2. Prediction measurements after binarization

TE

D

After binary processing of the predicted methylation values with the criteria mentioned in the Methods section, the SP, SE, ACC and MCC of the three models above were investigated. The box plots of the performance on the 24 chromosomes using the three prediction models are shown in Figure 4.

AC

CE P

From Figure 4, the prediction accuracies on each chromosome with the three prediction models are all approximately 0.9, and the MCCs are close to 0.8 (most MCCs in the RF model are larger than 0.8), indicating that the prediction results are precise. Similar to the correlation coefficients results, the RF model (The SP, SE, ACC and MCC are 0.94  0.02 , 0.87  0.03 , 0.91  0.01 and 0.81  0.02 , respectively) outperforms the LR and SVR models. Therefore, the RF model was selected for the following analyses. Furthermore, we compared the importance of all the features in the RF model by calculating the decrease of node impurity (Gini index). The box plot of the top 5 most important features among the 24 chromosomes were shown in Figure 5, from which we could also see the key contribution of 450K array feature in prediction. 3.3. Performance comparison Among the developed models for DNA methylation prediction, we selected EpiGRAPH for comparison because of its popularity and availability of user-friendly web service (Bock et al., 2009)(http://epigraph.mpi-inf.mpg.de/). Using the same training and test datasets used in our model construction, we obtained all parameters from EpiGRAPH, such as chromosome organization, histone methylation features, DNA sequences, DNA structure, epigenome_and_chromatin structure, and evolutionary history. In the prediction methods setup, SVM, LR and RF were selected to predict the methylation state of the test data. For the reason of computational cost, we implemented the cross-validation only on

ACCEPTED MANUSCRIPT

T

the methylation data of chromosomes 18 and Y, which showed the best and worst performance in our RF model. The prediction performance on chromosome 18 for EpiGRAPH is shown in Table.1. Generally, the SVM with RBF kernel showed better results than the other two methods in EpiGRAPH but still significantly worse than our RF model, which indicates that the feature of the 450K array data is informative in the prediction.

SC R

IP

The prediction results for the methylation data of chromosome Y are shown in Table 2. Similar to the results for chromosome 18, our model integrating the 450K array data and DNA sequence features outperforms the three models of EpiGRAPH.

MA

NU

3.4 Prediction performance in other tissues To investigate the capacity of the model for generalization to other tissues, we predicted the methylome of H9 ESC using the 450K array data of H9 ESC and sequence features. The Pearson and Spearman correlation coefficients are shown in Figure 6(a), and the concordances are shown in Figure 6(b). Except for the chromosome X, the Pearson correlation coefficients are 0.80  0.01 , and the concordances are 88%  1% , which indicates that the prediction model trained on H1 is applicable to H9.

TE

D

After the binary process of the predicted methylation values, we calculated the SE, SP, ACC and MCC of the prediction results. The results are shown in Figure 6(c). Similar to the correlation coefficient performance, the prediction results on H9 are comparable to the performance of the cross-validation test.

AC

CE P

Figure 6 shows that the performance metrics, including the correlation coefficients, SE, SP, ACC and MCC, are satisfactory, which indicates that our proposed model can be applied to predict the DNA methylation levels of other tissues. The outlier with relatively poor performance in Figure 6 is chromosome X, indicating that the methylation patterns of sex chromosome in different tissues have a larger degree of diversity. It is worth of noting that H9 and H1 are female and male respectively. The DNA sequence features are the same in different tissues, our model is generalizable because the 450K array data reflect tissue specificity. Therefore, the DNA methylome with base pair resolution in other tissues can be estimated using our proposed model as long as the 450K array data are provided. 4. Discussion We have proposed a methylation level prediction model for CpG loci by integrating 450K array data and DNA sequence features. The performance evaluated based on cross-validation testing and the inter-tissue prediction results demonstrate the method’s effectiveness. Thousands of 450K array data from various normal and diseased tissues that are publicly available in GEO and The Cancer Genome Atlas (TCGA) make our model easy to implement. For samples with only the DNA methylation landscape derived from 450K BeadChip array, our model can be used to obtain their DNA methylome with greater coverage (11 folds with our selected parameters), and has the potential to make novel discoveries when applied to the existing disease data. In this study, 450K array data are informative features for tissue specificity, and

ACCEPTED MANUSCRIPT

SC R

IP

T

sequence features depict general sequence preference for DNA methylation status. Here, only 1-mer, 2-mer base pairs and NpN ratios (N=A/C/G/T) were used as DNA sequence features in the model construction. Some other sequence features reported to be associated with DNA methylation, such as transcription factor binding sites (TFBS), RepeatMasker annotation and evolutionary information, could be added to the input features (Bock et al., 2006; Gaidatzis et al., 2014). Furthermore, the combination of 450K array data, DNA sequence features, histone modification marks and chromatin states data(Zhang et al., 2015) could show better performances if the corresponding histone modification data are available. However, the additional data are not always available and imputing DNA methylation levels based on 450K array alone would be invaluable to re-analyze the existing thousands of 450K array data on precious disease samples.

NU

Competing interests

The authors declare no competing financial interests.

MA

Acknowledgment

D

This work is partially supported by the National Natural Science Foundation of China under No. 61503061 and CIRM RB5-07012. 5. References

AC

CE P

TE

Barlow D. P. (2011). Genomic Imprinting: A Mammalian Epigenetic Discovery Model. Annual Review Genetics, Vol 45 45: 379-403. Baylin S. B., and Jones P. A. (2011). A decade of exploring the cancer epigenome - biological and translational implications. Nature Reviews Cancer 11: 726-734. Bhasin M., Zhang H., Reinherz E. L., and Reche P. A. (2005). Prediction of methylated CpGs in DNA sequences using a support vector machine. Febs Letters 579: 4302-4308. Bird A. (2002). DNA methylation patterns and epigenetic memory. Genes & Development 16: 6-21. Bock C. (2012). Analysing and interpreting DNA methylation data. Nat Rev Genet 13: 705-19. Bock C., Halachev K., Buch J., and Lengauer T. (2009). EpiGRAPH: user-friendly software for statistical analysis and prediction of (epi)genomic data. Genome Biology 10. Bock C., Paulsen M., Tierling S., Mikeska T., Lengauer T., and Walter J. (2006). CpG island methylation in human lymphocytes is highly correlated with DNA sequence, repeats, and predicted DNA structure. Plos Genetics 2: 243-252. Bock C., Walter J., Paulsen M., and Lengauer T. (2007). CpG island mapping by epigenome prediction. Plos Computational Biology 3: 1055-1070. Brinkman A. B., Simmer F., Ma K., Kaan A., Zhu J., and Stunnenberg H. G. (2010). Whole-genome DNA methylation profiling using MethylCap-seq. Methods 52: 232-236. Down T. A., Rakyan V. K., Turner D. J., Flicek P., Li H., Kulesha E., Graf S., Johnson N., Herrero J., Tomazou E. M., Thorne N. P., Backdahl L., Herberth M., Howe K. L., Jackson D. K., Miretti M. M., Marioni J. C., Birney E., Hubbard T. J. P., Durbin R., Tavare S., and Beck S. (2008). A Bayesian deconvolution strategy for immunoprecipitation-based DNA methylome analysis. Nature Biotechnology 26: 779-785. Ehrlich M., Gama-Sosa M. A., Huang L. H., Midgett R. M., Kuo K. C., McCune R. A., and Gehrke C. (1982). Amount and distribution of 5-methylcytosine in human DNA from different types of tissues of cells. Nucleic Acids Res 10: 2709-21. Ernst J., and Kellis M. (2015). Large-scale imputation of epigenomic datasets for systematic annotation of diverse human tissues. Nat Biotechnol 33: 364-76. Fan S. C., Zhang M. Q., and Zhang X. G. (2008). Histone methylation marks play important roles in predicting the methylation status of CpG islands. Biochemical and Biophysical Research Communications 374: 559-564. Fang F., Fan S. C., Zhang X. G., and Zhang M. Q. (2006). Predicting methylation status of CpG islands in the human brain. Bioinformatics 22: 2204-2209. Feltus F. A., Lee E. K., Costello J. F., Plass C., and Vertino P. M. (2003). Predicting aberrant CpG

ACCEPTED MANUSCRIPT

AC

CE P

TE

D

MA

NU

SC R

IP

T

island methylation. Proceedings of the National Academy of Sciences of the United States of America 100: 12253-12258. Feng P., Chen W., and Lin H. (2014). Prediction of CpG island methylation status by integrating DNA physicochemical properties. Genomics. Gaidatzis D., Burger L., Murr R., Lerch A., Dessus-Babus S., Schubeler D., and Stadler M. B. (2014). DNA Sequence Explains Seemingly Disordered Methylation Levels in Partially Methylated Domains of Mammalian Genomes. Plos Genetics 10. Gu H., Bock C., Mikkelsen T. S., Jager N., Smith Z. D., Tomazou E., Gnirke A., Lander E. S., and Meissner A. (2010). Genome-scale DNA methylation mapping of clinical samples at single-nucleotide resolution. Nature Methods 7: 133-U69. Heard E. (2004). Recent advances in X-chromosome inactivation. Current Opinion in Cell Biology 16: 247-255. Jones P. A. (2012). Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nature Reviews Genetics 13: 484-492. Laird P. W. (2010). Principles and challenges of genome-wide DNA methylation analysis. Nature Reviews Genetics 11: 191-203. Previti C., Harari O., Zwir I., and del Val C. (2009). Profile analysis and prediction of tissue-specific CpG island methylation classes. Bmc Bioinformatics 10. Sandoval J., Heyn H. A., Moran S., Serra-Musach J., Pujana M. A., Bibikova M., and Esteller M. (2011). Validation of a DNA methylation microarray for 450,000 CpG sites in the human genome. Epigenetics 6: 692-702. Stevens M., Cheng J. B., Li D., Xie M., Hong C., Maire C. L., Ligon K. L., Hirst M., Marra M. A., Costello J. F., and Wang T. (2014). Estimating absolute methylation levels at single-CpG resolution from methylation enrichment and restriction enzyme sequencing methods. Genome Res 23: 1541-53. Yang Y., Nephew K., and Kim S. (2012). A novel k-mer mixture logistic regression for methylation susceptibility modeling of CpG dinucleotides in human gene promoters. Bmc Bioinformatics 13. Zhang W., Spector T. D., Deloukas P., Bell J. T., and Engelhardt B. E. (2015). Predicting genome-wide DNA methylation using methylation marks, genomic position, and DNA regulatory elements. Genome Biol 16: 14. Zheng H., Wu H. W., Li J. P., and Jiang S. W. (2013). CpGIMethPred: computational model for predicting methylation status of CpG islands in human genome. Bmc Medical Genomics 6.

SC R

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE

D

MA

NU

Figure 1. The correlation coefficients and concordance of logistic regression (LR). (a) Pearson (in green) and Spearman (in yellow) correlation coefficients; (b) prediction concordance.

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE

D

MA

Figure 2. The correlation coefficients and concordance of random forest (RF). (a) Pearson (in green) and Spearman (in yellow) correlation coefficients; (b) prediction concordance.

SC R

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE

D

MA

NU

Figure 3. The correlation coefficients and concordance of support vector regression (SVR). (a) Pearson (in green) and Spearman (in yellow) correlation coefficients; (b) prediction concordance.

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE

Figure 4. Prediction performance after binarization with the LR, RF and SVR models. The RF model has the best performance.

TE

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

Figure 5. The top 5 most important features in RF model.

AC

CE P

TE

D

MA

NU

SC R

IP

T

ACCEPTED MANUSCRIPT

Figure 6. The prediction performance on the H9 ESC data using the RF model. (a) Pearson (in green) and Spearman (in yellow) correlation coefficients; (b) prediction concordance; (c) Sensitivity, Specificity, Accuracy and MCC.

ACCEPTED MANUSCRIPT

Table 1. Performance comparison between models from EpiGRAPH and our model on chr18. The RF performance in our model is better than the three models from EpiGRAPH. SE

SP

ACC

SVM in EpiGRAPH RFEpiGRAPH in EpiGRAPH LR in EpiGRAPH RF in our model

0.90 0.90 0.85 0.94

0.89 0.82 0.82 0.92

0.89 0.86 0.83 0.93

AC

CE P

TE

D

MA

NU

SC R

IP

T

Chr18

MCC 0.79 0.73 0.67 0.86

ACCEPTED MANUSCRIPT Table 2. Performance comparisons between models from the EpiGRAPH and our model for chrY. The RF performance in our model is better than the three models from EpiGRAPH. SE

SP

ACC

MCC

SVM in EpiGRAPH RF in EpiGRAPH LR in EpiGRAPH RF in our model

0.81 0.96 0.67 0.92

0.64 0.26 0.55 0.73

0.77 0.82 0.65 0.85

0.40 0.33 0.19 0.74

AC

CE P

TE

D

MA

NU

SC R

IP

T

ChrY

ACCEPTED MANUSCRIPT

Hightlights: Significant expansion of the coverage of Infinium HumanMethylation450 BeadChip in determining DNA methylation;



Performance evaluation and comparison of different computational models;



Generalizable models to expand DNA methylation measurement in diverse tissues.

AC

CE P

TE

D

MA

NU

SC R

IP

T

