Penalized logistic regression based on L1/2 penalty for high-dimensional DNA methylation data

2.1Network-regularization

In this research, n samples were used, D={(X1,y1),(X2,y2),…,(Xn,yn) where Xi=(xi⁢1,xi⁢2,…,xi⁢p) is the methylation β-value of the i-th sample and p represents the total number of CpG sites, the dependent variable yi is a binary variable where 0 implies controls and 1 implies cases. The logistic regression is:

where φ is the regression coefficients. The logistic log-likelihood is defined as:

φ was obtained by minimizing the log-likelihood. In high dimensional application, it is not appropriate to solve the logistic model directly and may result in overfitting. Hence, the regularization approaches are employed to aim at the overfitting problem. The sparse logistic regression can be laid out as Eq. (2) when a regularization term is added:

where P⁢(φ) is the penalty function.

Lasso (L1) and L2, a well-recognized regularization approach was used in previous methods. The L2 does not have sparsity and L1 has a sparsity less than Lq (0 <q< 1). Nonetheless, when q lies closer to zero, results show a sparser Lq and subsequently more challenging to converge. Therefore, some researchers [12] investigated the properties of Lq (0 <q< 1) regularization and demonstrated the L1/2 regularization is particularly essential and crucial. The performance between Lq penalty and L1/2 has no significant diversity whereas the L1/2 regularization is much more facile to solve. Accordingly, the L1/2 regularization can be laid out as Lq (0 <q< 1) regularization which exhibits unbiasedness as well as oracle properties [12, 13, 14]. In high-dimensional DNA methylation data, disease-related CpG sites are very limited and therefore, in practice, the L1/2 penalty methodology would be more significant than the L0, L1 and L2 approaches. Consequently, the L1/2 penalty was favored in our logistic regression model.

Some methods have been provided in order to tackle highly correlated variables. Elastic net penalty (L1+L2) and HLR (L1/2+L2) emphasizes a grouping effect and tend to smooth the coefficient profiles. However, the pathway information was neglected in these methods. To merge CpG sites deduction into the analysis of high-demensional methylation data, we extended a network-based regularization technique designed for the L1/2 penalty.

The methylation data displays a strong group effect and thus previous research used a fully connected network (Fc.net) to describe the correlated CpG sites group patterns within a gene. In methylation data, the group effect of CpG sites is not only present within one gene and present within one CpG island. There are overlapping parts between groups and these overlapping parts correlate with both parts respectively. With the different previous network, we set the overlapping part as the central node and connect it with other correlated parts (Fig. 1). The network not only has the genome information or CpG island information, but also the two aspects of information integrated into the network. It can better reflect the relevance of CpG site.

Figure 1.

a. Previous fully connected network. b. Central node fully connected network.

The network information is represented in a graphed structure with p-dimensional Laplacian matrix L={la⁢b}. It is defined as:

where da is the total number of connections at vertex a in graph.

The penalty function in Eq. (3) is:

where ||⋅||1/2 is a L1/2 norm and a∼b illustrate the variables which are linked to the a-th predictor. The sparsity and smoothness are controlled by the parameters λ1 and λ2.

The effectiveness of the penalty function reduced significantly when two negatively correlated predictors are interacted; the signs of coefficients are thus predicted and added to the Laplacian matrix to overcome problem:

The adaptive net function can be written as:

Based on |βa|≈𝑠𝑔𝑛⁢(βa*)⁢βa for βa≈βa*, the adaptive penalty function can be written as:

2.2The coordinate descent algorithm

To solve regularization models, the coordinate descent algorithm adopted as a competent tool. Regarding the coordinate descent algorithm, we referred to previous research [15, 16, 11] and Eq. (2) can be linearized by Taylor series expansion at current estimates φ*:

where zi=xiT⁢φ*+(yi-f*⁢(xi))/wi, wi=f*⁢(xi)⁢(1-f*⁢(xi)), f*⁢(xi)=𝑒𝑥𝑝⁢(xiT⁢φ*)/(1+𝑒𝑥𝑝⁢(xiT⁢φ*)).

Next, the estimator:

where z~i(a)=∑j≠ixi⁢j⁢φj*, g⁢(a)=∑a∼b|φb*|da⁢db and s⁢(σ,γ) is an enhanced L1/2 thresholding operator for the coordinate descent algorithm [12, 13, 14].

where ϕγ⁢(σ)=𝑎𝑐𝑟𝑜𝑠𝑠⁢((γ/8)⁢(|σ|/3)2/3).

3.1Analyses of simulated data

The performance of the proposed simulation study quoting the simulation from Teschendorff et al. [17] and Su and Wang [11] was analyzed and evaluated. There were 600 groups, which were divided into 100 groups, 150 groups and 7 sets of 50 groups in accordance to their number of CpG sites. Each group comprised of at least 1 CpG site up to 9 CpG sites reciprocally. In total, there were 2500 CpG sites.

First, we simulated variables with the group effect ranging between 0 and 1. So we performed an inverse logit transformation on a multivariate normal distribution variable to represent the β-values of the i-th CpG site in the g-th group.

where sg is the size of group, i.e. 1⩽sg⩽9. In this simulation model, we set μ=(-1,…,-1)T, xi,g ranging between 0 (unmethylated) to 1 (completely methylated). The relationship of CpG sites within group is shown by σ. The covariance matrix σ is presented as follows:

The first condition is autoregressive (AR) model, and the second condition is compound symmetric correlation model. We set three different correlation coefficients ρ= 0.2, 0.5 and 0.7 for all conditions [18, 19].

Second, given the regression coefficients φ based on previous research, φg=(φ1,g,…,φpg,g)T is the coefficient of CpG sites within the g-th group. After that, one group from each of the 9 different groups was selected to set the regression coefficients. At this step, there were 45 CpG sites which have been assigned the regression coefficients value. The regression coefficients φk,g are presented as:

when δ is the strength of the true signals. The other sets of regression coefficients were set to 0.

In the simulation models, there were 45 pathogenic CpG sites in a total of 2500 CpG sites. Lastly, the yi is given by Bernoulli distribution. For each simulation set, there were 200 cases and 200 controls. There were nine simulation conditions based on different parameters, for instance the strength of the true signals.

Figure 2.

The ROC curve of every model.

We repeated simulations 100 times for each condition. We then used the 10-fold cross-validation (CV) approach in the training set in order to tune the optimal regularization parameters of the Lasso, Elastic-Net (Enet), L1+Fc.net, L1/2, HLR (L1/2+L2), L1/2+Fc.net. Note that, the Enet, L1+Fc.net, HLR (L1/2+L2) and L1/2+Fc.net methods have two-dimensional parameter surfaces in the 10-CV approach. Afterwards, the logistic regressions with the estimated tuning parameters were employed to build different classifiers. Lastly, the attained classifiers were adopted to the test set for further classification and prediction.

Figure 3.

The histogram of correlation between CpG sites.

Figure 2 shows the receiver operating characteristic curve (ROC curve) for every model. The green solid line (L1/2+Fc.net) is closer to the upper left corner in the system than other line. So the effect of L1/2+Fc.net is at optimal for the other algorithm in each model. Table 1 shows the total area under the averaged ROC curves and the MSE of every model respectively. From Table 1 it can be seen that L1/2+Fc.net also has a very good performance within all models. In general, our proposed enhanced L1/2 model achieved preponderant accuracy rates in all models in comparison to the other methods (the Lasso, Enet, L1+Fc.net, L1/2 and HLR (L1/2+L2)).

Figure 4.

The boxplot of correlation between CpG sites.

3.2Analyses of real data

To further evaluate the effectiveness of our proposed method, in this section, we examined the DNA methylation (ovarian cancer) data generated from Illumina Infiniumm HumanMethylation 27K Beadchip [20]. The data is accessible from NCBI (http://www.ncbi.nlm.nih.gov/).

The data was generated by llumina Infiniumm HumanMethylation 27K Beadchip that contains 22727 CpG sites. We first removed samples which were low in BS conversion efficiency or low in CpG coverage. After that, a total of 207 genes contained more than 3 CpG sites and 295 CpG islands contained more than 3 CpG sites in the data; samples with error were removed. Lastly, there were 156 controls case samples (Healthy sample), 120 pre-treatment case samples and 122 post-treatment case samples. For these three cases, we calculated the maximum correlation of CpG sites in each group (gene and CpG island).

Figure 3a–c shows the histogram of maximum sample correlation between CpG sites within genes in control, pre-treatment and post-treatment case where Fig. 3d–f shows the histogram of maximum sample correlation between CpG sites within CpG islands in control, pre-treatment and post-treatment cases. Figure 4 shows the boxplot of maximum sample correlation between CpG sites in gene or CpG islands. Based on Figs 3 and 4, the results show that most CpG sites within the same group have high correlation in pre-treatment case samples and post-treatment case samples whereas the control case samples only show a significant correlation.

Table 2 shows the AUC for each method from real data analysis. In real data, the enhanced L1/2 model also achieved higher accuracy rates. Tables 3 and 4 show the top 20 selected CpG sites for all methods. We further validated the chosen genes from the GeneCards Database (http://www.genecards.org). In Table 3, when comparing pre-treatment cases with controls, the algorithm L1/2+ central node Fc.net (gene and CpG island) found various genes (TIAM1 [21, 22], CST7 [25], TCEAL7 [23, 24] and RNF11 [26]) that were not found by L1+Fc.net (gene) and Enet. Likewise, HLR (L⁢1/2+L⁢2) was unable to find these genes (CST7, TCEAL7 and RNF11) where these genes (TIAM1, CST7, TCEAL7 and RNF11) were found to be correlated with cancer in previous research. On one hand, all methods were able to find genes (MPO [27, 28], PTPRCAP [29]); on the other hand, network penalty methods were able to find genes (CD248 [31] and HOXB5 [32]). In Table 4, our algorithm L1/2+ central node Fc.net (gene and CpG island) also found genes (CD200 [30], SRD5A2L [34], ENO1 [33]) which have not been found in L1+Fc.net(gene) and Enet. Additionally, gene CD200 and gene SRD5A2L also proved to be related to cancer. The L1+ Fc.net and L1/2+ central node Fc.net (gene and CpG island) algorithm, which has gene and island network information, also found gene (TNFAIP8 [35]) which was not found by Enet/HLR.