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): β-value=MM+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): M-value=log2(M+α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., 2010).

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): β=2M2M+1;M=log2(β1-β) According to the conclusions in Du et al. (2010), 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.

Gaussian Mixture Model (GMM) has been widely applied as a clustering method in analyzing gene-expression microarray data (Yeung et al., 2001; Pan et al., 2002) and used to detect differential gene expression (McLachlan et al., 2006). 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 (M_i) is defined as: (1)p(Mi;θ)=∑k=1KπkN(Mi|μk,σk2)

where p(M_i, θ) represents the model density for M_i with unknown parameter vector θ, K is the number of different methylation states (components), N(Mi|μk,σk2) is the probability density function of the kth Gaussian component, and π_k is the mixing proportions which satisfy the constraint that ∑k=1Kπk=1 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 μTS and σTS denote the mean value and standard deviation where S ∈ (U, H, F) and T ∈ (I, II). K_I and K_II 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 (U_T), hemimethylated (H_T), or hypermethylated (F_T) states by using the maximum probability criterion. Let UTL (UTR) denote the U_T probes with M-values smaller (larger) than μTU, and let FTL (FTR) represent the F_T probes with M-values smaller (larger) than μTF where T ∈ (I, II). Then, we calculate the probabilities of UIIL probes, i.e., (2)p=P(MUIIL|μIIU,(σIIU)2) where P represents the cumulative distribution function of the Gaussian component. These probabilities are transformed back to quantiles (M-value) by using the parameters μIU and σIU of type I probes, i.e., (3)q=P-1(p|μIU,(σIU)2) where P⁻¹ returns the value of the inverse cumulative density function given the probability p and q is the normalized M-values for UIIL. The similar operation is performed on FIIR probes.

Then, we merge the UIIR, H_II, and FIIL probes into one set G on which a conformal (shift + dilation) transformation is performed. Some parameters are identified as minG = min{M_G}, maxG = max{M_G} and ΔGM=maxG-minG. Similarly, the minimum value of FIIR and the maximum value of UIIL are also identified, i.e., minF=min{FIIR} and maxU=max{UIIL}. Two distance values can be calculated as ΔUG=minG-maxU ΔGF=minF-maxG The new normalized maximum and minimum values of G-probes are expected to satisfy the constraint that maxG′=min{FIIR′}-ΔGF minG′=max{UIIL′}+ΔUG where FIIR′ and UIIL′ are new normalized values for FIIR and UIIL, respectively. So the new normalized range value of set G is ΔGM′=maxG′-minG′. The normalized M-values of set G, MGII′, is calculated by (4)MGII′=minG′+df(MGII-minG) where df=ΔGM′/ΔGM is the dilation factor. So, the normalized M-values for type II probes consist of q for UIIL, MGII′, and q for FIIR. MII′=(qUIIL, MGII′,qFIIR) 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) K_I will be changed to 4 automatically when μIF-σIF is smaller than μIIF-σIIF in order to ensure that μIF can always be larger than μIIF and avoid the presence of an unexpected peak in transformed M-values of hypermethylated type II probes, (iii) if K_I = 4, the F_I will be the set of probes belonging to the component with the largest μ, while the U_I contains the probes belonging to the component with the smallest μ and the other two components are assigned to H_I, (iv) no thresholds need to be set by default or estimated by manual to distinguish the three different states of DNA methylation.

We selected several public 450K datasets as following:

Dataset 1: GSE29290 downloaded from GEO considered in Dedeurwaerder et al. (2011). 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., 2011) to evaluate the performance of different methods.

Dataset 2: GSE38268 downloaded from GEO considered in Lechner et al. (2013) which consists of 6 fresh frozen HNC samples. We selected 5 samples as same as (Teschendorff et al., 2012), of which 2 were HPV+ and 3 HPV− (GSM937820 to GSM937824).

Dataset 3: GSE38266 downloaded from GEO considered in Lechner et al. (2013) 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. (2017) which contains 6 HPV+ HNC samples and 5 HPV− HNC samples.

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).

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).

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. (2011), 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.

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. (2012) 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, 2005) and an FDR threshold 0.35 which was as same as (Teschendorff et al., 2012). 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. (2012) 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. (2012) 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.