Abstract
The Illumina Infinium HumanMethylation450 Beadchips have been widely utilized in epigenome-wide association studies (EWAS). However, the existing two types of probes (type I and type II), with the distribution of measurements of probes and dynamic range different, may bias downstream analyses. Here, we propose a method, MGMIN (M-values Gaussian-MIxture Normalization), to correct the probe designs based on M-values of DNA methylation. Our strategy includes fitting Gaussian mixture distributions to type I and type II probes separately, a transformation of M-values into quantiles and finally a dilation transformation based on M-values of DNA methylation to maintain the continuity of the data. Our method is validated on several public datasets on reducing probe design bias, reducing the technical variation and improving the ability to find biologically differential methylation signals. The results show that MGMIN achieves competitive performances compared to BMIQ which is a well-known normalization method on β-values of DNA methylation.
1. Introduction
DNA methylation, as a well-known epigenetic marker, plays an essential role in biological processes and complex genetic diseases like cancer and diabetes (Irizarry et al., ; Paul et al., ). The Illumina Infinium HumanMethylation450 (450K) BeadChip (Bibikova et al., ) provides measurements of the level of methylation at over 480K CpG sites and has been widely used in epigenome-wide association studies (EWAS) and large-scale projects, such as The Cancer Genome Atlas (TCGA). The probes in the Infinium 450K BeadChip come in two different designs, type I (n = 135,501) and type II (n = 350,076), in order to increase the genomic coverage of the assay. However, the methylation values (β-values or M-values) derived from the two types of designs exhibit different distributions. Particularly, the type I probes possess a larger range of measurement than the type II probes (Dedeurwaerder et al., ). The differences between the two types of probe designs may impact the downstream analyses.
Several approaches have been published to correct the probe design bias. A peak-based correction (PBC) method normalizes type II probes to render them comparable with type I probes (Dedeurwaerder et al., ). In fact, PBC gets poor performance when the density distribution of methylation values does not show well-defined peaks. SQN (Touleimat and Tost, 2012) and SWAN (Maksimovic et al., ) select subset of probes with similar biological category to adjust the probe design bias. Beta MIxture Quantile dilation (BMIQ) is a model-based normalization approach to correct β-values of type II probes according to the beta distribution of β-values of type I probes, which appears to outperform PBC, SQN, and SWAN (Teschendorff et al., ).
In this work, we propose a method to correct the probe design bias based on the Gaussian Mixture Model (GMM) of the M-values of DNA methylation, which is called M-value Gaussian-MIxture Normalization (MGMIN). The method includes three steps: (i) fit Gaussian-mixture distributions to type I and type II probes separately, (ii) utilize a transformation of M-values into quantiles, (iii) perform a dilation transformation based on M-values to maintain the continuity of the data. We evaluate MGMIN using several independent datasets in terms of reducing the replicate technical variance and correcting the type II bias. By comparison with BMIQ, the results show that MGMIN improves the overall performance of normalization.
2. Materials and Methods
2.1. Measure DNA Methylation With M-value
The β-value of DNA methylation for each probe is defined by the ratio of the methylated intensity (M) and the overall intensity (sum of methylated intensity and unmethylated intensity: M + U): where α is a constant offset (by default, α = 100) to regularize the β-value when the overall intensity is low. The β-value falls between 0 and 1 which follows a Beta distribution naturally. A β-value of 0 indicates the CpG site of the measured sample is fully unmethylated and a value of 1 indicates that the CpG site is completely methylated.
The M-value is calculated by the log2 ratio of the methylated intensity (M) vs. the unmethylated intensity (U): where α here is also an offset (by default, α = 1) to counteract the big changes caused by small intensity estimation errors. An M-value close to zero indicates that the measured CpG site is about hemimethylated. A positive M-value suggests that more copies of the measured CpG site are methylated than unmethylated and a negative M-value means more copies of the CpG site are unmethylated. The M-value has been widely used in two-color expression microarray analysis (Du et al., ).
Due to more than 95% CpG sites have intensities more than 1,000 in Illumina methylation data, the α in β-value and M-value has an insignificant effect on observed results. So the relationship between β-value and M-value is shown as (with α ignored): According to the conclusions in Du et al. (), the M-value is more statistically valid in an analysis by modeling the distribution of M-values because of it's homoscedastic. So we choose to adjust the M-values of type II probes into the distribution property of type I probes to correct the probe design bias.
2.2. MGMIN: M-value Gaussian-MIxture Normalization
Gaussian Mixture Model (GMM) has been widely applied as a clustering method in analyzing gene-expression microarray data (Yeung et al., ; Pan et al., ) and used to detect differential gene expression (McLachlan et al., ). In this paper, we apply GMM to distinguish different methylation states of CpG sites for further correction. The M-values of a single 450K microarray can be viewed as a finite Gaussian mixture model of several methylation states (hypomethylated-U, hemimethylated-H, hypermethylated-F). The probability density function of the M-value for a single CpG site (Mi) is defined as:
where p(Mi, θ) represents the model density for Mi with unknown parameter vector θ, K is the number of different methylation states (components), is the probability density function of the kth Gaussian component, and πk is the mixing proportions which satisfy the constraint that and 0 ≤ πk ≤ 1. The parameter vector θ consists of the mixing proportions πk, the mean value μk and the standard deviation σk, which can be estimated by the EM algorithm.
Next, we describe MGMIN in detail. First, M-values of type I and type II probes are modeled by GMM separately. Let and denote the mean value and standard deviation where S ∈ (U, H, F) and T ∈ (I, II). KI and KII are the numbers of components for type I and type II probes, which are both set as 3 by default.
Second, each probe is assigned to hypomethylated (UT), hemimethylated (HT), or hypermethylated (FT) states by using the maximum probability criterion. Let () denote the UT probes with M-values smaller (larger) than , and let () represent the FT probes with M-values smaller (larger) than where T ∈ (I, II). Then, we calculate the probabilities of probes, i.e., where P represents the cumulative distribution function of the Gaussian component. These probabilities are transformed back to quantiles (M-value) by using the parameters and of type I probes, i.e., where P−1 returns the value of the inverse cumulative density function given the probability p and q is the normalized M-values for . The similar operation is performed on probes.
Then, we merge the , HII, and probes into one set G on which a conformal (shift + dilation) transformation is performed. Some parameters are identified as minG = min{MG}, maxG = max{MG} and . Similarly, the minimum value of and the maximum value of are also identified, i.e., and . Two distance values can be calculated as The new normalized maximum and minimum values of G-probes are expected to satisfy the constraint that where and are new normalized values for and , respectively. So the new normalized range value of set G is . The normalized M-values of set G, , is calculated by where is the dilation factor. So, the normalized M-values for type II probes consist of q for , , and q for . Lastly, the normalized M-values are transformed to β-values.
There are some important points to notice: (i) the initial values for μ and σ in EM algorithm are set as (−4,0,4) and (1,1,1) and small perturbations to the initial μ and σ will not affect the final model because MGMIN captures the natural property of the M-value of DNA methylation, (ii) KI will be changed to 4 automatically when is smaller than in order to ensure that can always be larger than and avoid the presence of an unexpected peak in transformed M-values of hypermethylated type II probes, (iii) if KI = 4, the FI will be the set of probes belonging to the component with the largest μ, while the UI contains the probes belonging to the component with the smallest μ and the other two components are assigned to HI, (iv) no thresholds need to be set by default or estimated by manual to distinguish the three different states of DNA methylation.
2.3. Datasets
We selected several public 450K datasets as following:
Dataset 1: GSE29290 downloaded from GEO considered in Dedeurwaerder et al. (). We used the three replicates (GSM15136, GSM15137 and GSM15138) from the HCT116WT cell-line and matched bisulfite pyrosequencing (BPS) date for nine type II probes of sample GSM815138 (r3) (Table 1 in Dedeurwaerder et al., ) to evaluate the performance of different methods.
Dataset 2: GSE38268 downloaded from GEO considered in Lechner et al. () which consists of 6 fresh frozen HNC samples. We selected 5 samples as same as (Teschendorff et al., ), of which 2 were HPV+ and 3 HPV− (GSM937820 to GSM937824).
Dataset 3: GSE38266 downloaded from GEO considered in Lechner et al. () which contains 21 FFPE HPV+ HNSCC samples and 21 FFPE HPV− HNSCC samples. Note that the entire quality of the dataset GSE38266 is not high.
Dataset 4: GSE95036 downloaded from GEO considered in Degli Esposti et al. () which contains 6 HPV+ HNC samples and 5 HPV− HNC samples.
3. Results
3.1. MGMIN Needs No Default Initial Values of Parameters
Similar to the mixture model of BMIQ, MGMIN applies Gaussian mixture models for M-values instead of beta-mixture models for β-values. MGMIN also uses quantile information to correct the M-values of the type II probes into a distribution which is comparable with that of type I probes. MGMIN complies the inherent Gaussian mixture distributions for M-values of type I and type II probes to avoid setting any parameters manually, which is different from the default breakpoints in BMIQ. Thus, MGMIN needs less manual intervention than BMIQ. However, MGMIN is slightly inferior to BMIQ on some dataset (Table 1) due to the entire low quality of the dataset. Note that the PPV of BMIQ on Dataset 3 is lower than that of no normalization (RAW).
Table 1
| Metric | Raw | BMIQ | MGMIN |
|---|---|---|---|
| nDMP | 51 (51a) | 239 (252a) | 220 |
| nTP (GSE38266) | 16 (13a) | 55 (51a) | 37 |
| PPV (GSE38266) | 0.31 (0.25a) | 0.23 (0.20a) | 0.17 |
| nTP (GSE95036) | 3 | 13 | 27 |
| PPV (GSE95036) | 0.06 | 0.05 | 0.12 |
Comparison of MGMIN and BMIQ on detecting the differentially methylated probes (DMPs) associated with HPV status was performed by counting the number of DMPs (Dataset 2), the number of validated differentially methylated probes (nTPs) (Dataset 3: GSE38266 and Dataset 4: GSE95036) and corresponding estimates for the positive predictive value (PPV = nTP/nDMPs).
Values reported in Teschendorff et al. ().
3.2. MGMIN Reduces Technical Variation
MGMIN is applied to Dataset 1 to identify the ability to improve reproducibility. The standard deviation (SD) for each probe across the three replicates was computed using no normalization (RAW), BMIQ, and MGMIN separately. As can be seen in Figure 1, both MGMIN and BMIQ almost made the density curves for the three replicates coincide with each other and reduced the technical variation significantly compared to no normalization. Compared to BMIQ, the standard deviation for type II probes adjusted by MGMIN is smaller (Figure 2). MGMIN also provided significant reduction of average absolute difference in β-values of type II probes between two samples in each of the three pairs of the three replicates (Figure 3).
Figure 1
Figure 2
Figure 3
3.3. MGMIN Reduces Probe Design Bias
MGMIN reduces the probe design bias via correcting the M-values of the type II probes such that the distribution curves for the M-values of the type I and type II probes show similar dynamic ranges and peaks (Figure 4). In Dedeurwaerder et al. (), the β-values for nine probes of type II by bisulfite pyrosequencing technique for sample GSM815138 (r3) were provided, which can be used as a gold-standard to evaluate the performance of different correction methods. Hence, we compared the normalized results of the nine type II probes in 450K arrays by MGMIN and BMIQ. As shown in Figure 5, although MGMIN performed slightly worse than BMIQ at the maximum value of the absolute deviation from BPS data, MGMIN significantly reduced the type II bias than BMIQ and raw data in terms of mean and root mean square error (RMSE) of the absolute deviation from the matched BPS values.
Figure 4
Figure 5
3.4. MGMIN Robustly Finds Informative Differential Methylation Probes Associated With HPV Status
The goal of a bias correction approach is to reduce the technical variation and identify the biological informative signals at the same time. We used a strategy similar to Teschendorff et al. () to compare the result between MGMIN and BMIQ in identifying the differential methylation probes (DMPs) associated with HPV status. First, Dataset 2 consisting of two HPV+ and three HPV− fresh frozen HNC samples were used as the training set to obtain the DMPs associated with HPV status by the limma method (Smyth, ) and an FDR threshold 0.35 which was as same as (Teschendorff et al., ). Both Dataset 3 and Dataset 4 described in the methods section were used as test set. We reanalyzed Dataset 2 and got similar numbers of DMPs to those reported in Teschendorff et al. () with no normalization method (Raw) or BMIQ method (shown in Table 1). The results in Table 1 shows that the positive predictive value (PPV) of MGMIN is slightly less than BMIQ in terms of GSE38266 (Dataset 3) whereas MGMIN outperforms BMIQ in GSE95036 (Dataset 4). The reason for MGMIN slightly inferior to BMIQ in Dataset 3 may be the entire low quality of the dataset (see Figure 6) which is that the ratio of samples passing filters is <0.9 (r = 0.88) under the least restrictive condition. Let τp represent the p-value threshold for bad probes and τr represent the threshold for the ratio of bad probes in a sample. The maximum value of τr is set to 0.3 here in our opinion because a sample with more than 30% bad probes is vulnerable. We can get the same test dataset from GSE38266 with the one described in Teschendorff et al. () which consists of 18 HPV+ and 14 HPV− samples under the following conditions: (i) τp = 1e − 4 or 1e − 3 and τr = 0.2 or 0.25, (ii) τp = 1e − 2 and τr = 0.1 or 0.15. Overall, MGMIN identified more true positive features than BMIQ.
Figure 6
4. Discussions
In this paper, we have proposed a method called MGMIN for correcting the probe design bias of type II probes in Illumina Infinium 450K BeadChips, which can reduce the technical variation and improve the ability to find biologically differential methylation signals. We have shown that MGMIN outperforms BMIQ on multiple evaluation datasets in correcting the type II design bias and improving the data quality.
Similar to BMIQ, MGMIN uses quantile information to correct the M-values of type II probes while leaving the M-values of type I probes unchanged. The three-state beta-mixture distribution model in BMIQ sets two default breakpoints (0.2, 0.75) to divide the β-values into three classes: hypomethylated, hemimethylated, and hypermethylated, which works well for most cases. However, the result curves of BMIQ show obviously inconsistent in some samples with high heterogeneity. We set 3 or 4 classes for probes depending on the result of to ensure that the fitted hypermethylated component of type II probes can be located in the left of the hypermethylated component of type I probes, which can partly eliminate the effects of the heterogeneity of samples.
Based on the results of Dataset 3, we think the high quality of dataset is the base of normalization, in other words, there is no meaning to correct the samples with low quality. It should be pointed out that the parameter estimation of MGMIN is slower than that of BMIQ (about 1.5 times), which can be relieved by reducing the number of iterations.
MGMIN can be used in the 450K methylation data preprocessing with other methods to normalize the values of the two type probes and improve the data quality.
Statements
Data availability statement
The datasets for this study can be found in GEO: GSE29290, GSE38268, GSE38266, and GSE95036.
Author contributions
ZW performed the experiments and wrote the manuscript. All authors read and revised the final manuscript.
Funding
This work has been supported by the National Key Research and Development Program of China (Nos: 2017YFC1201201 and 2017YFC0907503).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
References
1
BibikovaM.BarnesB.TsanC.HoV.KlotzleB.LeJ. M.et al. (2011). High density DNA methylation array with single CPG site resolution. Genomics98, 288–295. 10.1016/j.ygeno.2011.07.007
2
DedeurwaerderS.DefranceM.CalonneE.DenisH.SotiriouC.FuksF. (2011). Evaluation of the infinium methylation 450K technology. Epigenomics3, 771–784. 10.2217/epi.11.105
3
Degli EspostiD.SkliasA.LimaS. C.Beghelli-de la Forest DivonneS.CahaisV.Fernandez-JimenezN.et al. (2017). Unique DNA methylation signature in HPV-positive head and neck squamous cell carcinomas. Genome Med. 9:33. 10.1186/s13073-017-0419-z
4
DuP.ZhangX.HuangC.-C.JafariN.KibbeW. A.HouL.et al. (2010). Comparison of beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics11:587. 10.1186/1471-2105-11-587
5
IrizarryR. A.Ladd-AcostaC.WenB.WuZ.MontanoC.OnyangoP.et al. (2009). The human colon cancer methylome shows similar hypo-and hypermethylation at conserved tissue-specific cpg island shores. Nat. Genet. 41, 178–186. 10.1038/ng.298
6
LechnerM.FentonT.WestJ.WilsonG.FeberA.HendersonS.et al. (2013). Identification and functional validation of HPV-mediated hypermethylation in head and neck squamous cell carcinoma. Genome Med. 5:15. 10.1186/gm419
7
MaksimovicJ.GordonL.OshlackA. (2012). Swan: subset-quantile within array normalization for illumina infinium humanmethylation450 beadchips. Genome Biol. 13:R44. 10.1186/gb-2012-13-6-r44
8
McLachlanG. J.BeanR.JonesL. B.-T. (2006). A simple implementation of a normal mixture approach to differential gene expression in multiclass microarrays. Bioinformatics22, 1608–1615. 10.1093/bioinformatics/btl148
9
PanW.LinJ.LeC. T. (2002). Model-based cluster analysis of microarray gene-expression data. Genome Biol. 3:research0009-1. 10.1186/gb-2002-3-2-research0009
10
PaulD. S.TeschendorffA. E.DangM. A.LoweR.HawaM. I.EckerS.et al. (2016). Increased DNA methylation variability in type 1 diabetes across three immune effector cell types. Nat. Commun. 7:13555. 10.1038/ncomms13555
11
SmythG. K. (2005). Limma: linear models for microarray data, in Bioinformatics sand Computational Biology Solutions Using R and Bioconductor, eds GentlemanR.CareyV.HuberW.IrizarryR.DudoitS. (New York, NY: Springer), 397–420. 10.1007/0-387-29362-0_23
12
TeschendorffA. E.MarabitaF.LechnerM.BartlettT.TegnerJ.Gomez-CabreroD.et al. (2012). A beta-mixture quantile normalization method for correcting probe design bias in illumina infinium 450K DNA methylation data. Bioinformatics29, 189–196. 10.1093/bioinformatics/bts680
13
TouleimatN.TostJ. (2012). Complete pipeline for infinium® human methylation 450K beadchip data processing using subset quantile normalization for accurate dna methylation estimation. Epigenomics4, 325–341. 10.2217/epi.12.21
14
YeungK. Y.FraleyC.MuruaA.RafteryA. E.RuzzoW. L. (2001). Model-based clustering and data transformations for gene expression data. Bioinformatics17, 977–987. 10.1093/bioinformatics/17.10.977
Summary
Keywords
DNA methylation, design bias, normalization, M-value, Gaussian mixture model, Illumina Infinium 450K
Citation
Wang Z, Liu Y and Wang Y (2020) MGMIN: A Normalization Method for Correcting Probe Design Bias in Illumina Infinium HumanMethylation450 BeadChips. Front. Genet. 11:538492. doi: 10.3389/fgene.2020.538492
Received
27 February 2020
Accepted
21 September 2020
Published
27 October 2020
Volume
11 - 2020
Edited by
Zhongyu Wei, Fudan University, China
Reviewed by
Teng Zhixia, Northeast Forestry University, China; Zhen Tian, Zhengzhou University, China
Updates
Copyright
© 2020 Wang, Liu and Wang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yadong Wang ydwang@hit.edu.cn
This article was submitted to Computational Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.