Network-Based Logistic Classification with an Enhanced L 1/2 Solver Reveals Biomarker and Subnetwork Signatures for Diagnosing Lung Cancer.

Identifying biomarker and signaling pathway is a critical step in genomic studies, in which the regularization method is a widely used feature extraction approach. However, most of the regularizers are based on L 1-norm and their results are not good enough for sparsity and interpretation and are asymptotically biased, especially in genomic research. Recently, we gained a large amount of molecular interaction information about the disease-related biological processes and gathered them through various databases, which focused on many aspects of biological systems. In this paper, we use an enhanced L 1/2 penalized solver to penalize network-constrained logistic regression model called an enhanced L 1/2 net, where the predictors are based on gene-expression data with biologic network knowledge. Extensive simulation studies showed that our proposed approach outperforms L 1 regularization, the old L 1/2 penalized solver, and the Elastic net approaches in terms of classification accuracy and stability. Furthermore, we applied our method for lung cancer data analysis and found that our method achieves higher predictive accuracy than L 1 regularization, the old L 1/2 penalized solver, and the Elastic net approaches, while fewer but informative biomarkers and pathways are selected.

1. Introduction

Identifying molecular biomarker or signaling pathway involved in a phenotype is a particularly important problem in genomic studies. Logistic regression is a powerful discriminating method and has an explicit statistical interpretation which can obtain probabilities of classification regarding the class label information.

A key challenge in identifying diagnosis or prognosis biomarkers using the logistic regression model is that the number of observations is much smaller than the size of measured biomarkers in most of the genomic studies. Such limitation causes instability in the algorithms used to select gene marker. Regularization methods have been widely used in order to deal with this problem of high dimensionality. For example, Shevade and Keerthi proposed the sparse logistic regression based on the Lasso regularization [1, 2]. Meier et al. investigated logistic regression with group Lasso [3]. Usually, the Lasso type procedures are often called L 1-norm type regularization methods. However, L 1 regularization may yield inconsistent selections when applied to variable selection in some situations [4] and often introduces the extra bias in the estimation [5]. In many genomic studies, we need a sparser solution for interpretation and accurate outcomes, but L 1 regularization has a gap to meet these requirements. Thus, a further improvement of regularization is urgently required. L (0 < q < 1) regularization can assuredly generate more sparse and precise solutions than L 1 regularization. Moreover, L 1/2 penalty can be taken as a representative of L (0 < q < 1) penalty and has demonstrated many attractive properties which do not appear in some L 1 regularization approaches, such as unbiasedness, sparsity, and oracle properties [6-8].

So far, we observed dense molecular interaction information about the disease-related biological processes and gathered it through databases focused on many aspects of biological systems. For example, BioGRID records collected various biological interactions from more than 43,468 publications [9]. These regulatory relationships are usually represented by a network. Combining these pieces of graphic information extracted from the biological process with an analysis of the gene-expression data had provided useful prior information to detective noise and removes confounding factors from biological data for several classification and regression models [10-14].

Inspired by the aforementioned methods and ideas, here, we define a network-constrained logistic regression model with L 1/2 penalty following the framework established by [11], where the predictors are based on the gene-expression data with biologic network knowledge. The proposed model is aimed at identifying some biomarkers and subnetworks regarding diseases. In order to achieve a better prediction, we use an enhanced half thresholding algorithm for L 1/2  regularization, which is more efficient than the old half thresholding approach in the literature [6, 15, 16].

The rest of the paper is organized as follows. In Section 2, we proposed a new version of the network-constrained logistic regression model with L 1/2 regularization. In Section 3, we presented an enhanced half thresholding method for L 1/2 regularization and the corresponding coordinate descent algorithm. In Section 4, we evaluated the performance of our proposed approach on the simulated data and presented the applications of the proposed methods to an analysis of lung cancer data. We concluded the paper with Section 5.

2. L 1/2 Penalized Network-Constrained Logistic Regression Model

Generally, assuming that dataset D has n samples, D = {(X 1, y 1), (X 2, y 2),…, (X , y )}, where X = (x , x ,…, x ) is ith sample with p genes and y is the corresponding variable that takes a value of 0 or 1. Define a classifier f(x) = e /(1 + e ) and the logistic regression is defined aswhere β = (β 1,…, β ) are the coefficients to be estimated. We can obtain β by minimizing the log-likelihood function of the logistic regression. Following [11], to combine biological network with an analysis of the gene microarray data, we used a Laplacian constraint approach here. Consider a graph G = (V, E), where V is the set of genes that meet p explanatory variables and E is the set of edges. If gene u and gene v are connected, then there is an edge between gene u and gene v, which is denoted by E = 1; else E = 0. w denotes the weight of edge E . The normalized Laplacian matrix L for G is defined bywhere d and d are the degrees of genes u and v, respectively. The degrees of gene u (or v) describe the number of the edges that connected with u (or v). For λ ≥ 0, the network-constrained logistic regression model is presented aswhere the first term in (3) is the log-likelihood function of the logistic model and the second term is a network constraint based on the Laplacian matrix, which induces a smooth solution of β on the graph.

Directly computing (3) performs poorly for both prediction and biomarker selection purposes when the gene number p ≫ the sample size n. Therefore, the regularization approach is vitally needed. When adding a regularization term to (3), the sparse network-constrained logistic regression can be written aswhere λ 1 > 0 is a regularization parameter. In Zhang et al. [13], the authors used Lasso (L 1) which has the regularization term P(β) = ∑ |β | to penalize (4). However, the result of the Lasso type (L 1) regularization is not good enough for interpretation, especially in genomic research. Besides this, L 1 regularization is asymptotically biased [17, 18]. To improve the solution's sparsity and its predictive accuracy, we need to think beyond L 1 regularization to L penalties. In mathematics, L   (0 < q < 1) type regularization |β| = ∑|β| with the lower value of q would lead to better solutions with more sparsity and gives asymptotically unbiased estimates [17]. Moreover, L 1/2 penalty can be taken as a representative of L   (0 < q < 1) penalty and has permitted an analytically expressive thresholding representation [6, 7]. Therefore, we proposed a novel L 1/2 net approach based on L 1/2 regularization to penalize the network-constrained logistic regression model, as shown in where |β|1/2 = ∑ |β |1/2.

3. A Coordinate Descent Algorithm for the Network-Constrained Logistic Model with the Enhanced L 1/2 Thresholding Operator

L 1/2 penalty function is nonconvex, which raises numerical challenges in fitting the models. Recently, the coordinate descent algorithms [19] for solving nonconvex regularization models (SCAD [20], MCP [21]) have shown significant efficiency and convergence [22]. Since the computational burden increases only linearly with the feature number p, the coordinate descent algorithm can be a powerful tool for solving high-dimensional problems. Its standard procedure can be demonstrated as follows: for every coefficient β   (j = 1,2,…, p), to partially optimize the target function with respect to β , and fix the remaining elements β (k = 1,2,…, p  and  k ≠ j) at their most recently updated values. The specific form of updating β depends on the thresholding operator of the penalty.

In this paper, we present an enhanced L 1/2 thresholding operator for the coordinate descent algorithm:where φ (ω) = arccos⁡((λ/8)(|ω|/3)−3/2), π = 3.14, , and as the partial residual for fitting β .

Remark. This enhanced L 1/2 thresholding operator outperforms the old L 1/2 thresholding (3/4)(λ)2/3 introduced in [6, 15, 16]. We know that the quantity of the regularization solutions depends seriously on the value of the regularization parameter λ. Based on this enhanced L 1/2 thresholding operator, when λ is chosen by some efficient strategies for the parameter tuning, such as cross validation, the convergence of algorithm (6) is proved [7].

The Laplacian matrix L is nonnegative definite; thus, it can be written as L = SS by Cholesky decomposition. Following C. Li and H. Li [11] approach, (4) can be expressed aswhere , , , and γ is the regularization parameter and can be expressed as .

One-term Taylor series expansion for (7) can be written aswhere is the estimated response and is the weight for the estimated response. is the evaluated value under the current parameters. Thus, we can redefine the partial residual for fitting current as and . The procedure of the coordinate descent algorithm for L 1/2 penalized network-constrained logistic model is described as follows.

Algorithm 1 (the coordinate descent algorithm for L 1/2 penalized network-constrained logistic model).

We consider the following.

Step 1. Initialize all β (m) ← 0 (j = 1,2,…, p) and y , X , and set m ← 0, γ chosen by cross validation.

Step 2. Calculate Z(m) and W(m) and approximate the loss function (8) based on the current β(m).

Step 3. Update each β (m) and cycle over j = 1,…, p, until β (m) does not change.

Step 3.1. Compute Z ((m) ← ∑ W (m)(Z (m) − ∑ x β (m)) and .

Step 3.2. Update β (m) ← Enhanced_Half(ω (m), γ).

Step 4. Let m ← m + 1, β(m + 1) ← β(m).

If β(m) dose not converge, then repeat Steps 2 and 3.

4. Simulation and Application

4.1. Analyses of Simulated Data

We evaluate the performance of four methods: the network-constrained logistic regression models with L 1 regularization (L 1 net), L 1/2 regularization with old thresholding value (3/4)(λ)2/3 (L 1/2 net) and with the enhanced thresholding value (enhanced L 1/2 net), and the Elastic net regularization approach (Elastic net). We first simulated the graph structure to mimic gene regulatory network: assuming that the graph consists of 200 independent transcription factors (TFs) and each TF regulates 10 unlike genes, so there are a total of 2200 variables, X = (x 1, x 2,…, x ), p = 2200. The training and the independent test data sets include the sample sizes of 100, respectively. Each TF x and its regulated genes x were generated by the normal distribution N(0,1). We set the correlation rate between x and its regulated gene x as 0.75, x = (1 − 0.75) × x + (0.75) × x . The binary responder y (1 ≤ i ≤ 100), which is associated with the matrix X of TFs and their regulated genes, is calculated based on the following formula and rule:where , for Model 1, and ε ~ N(0, σ 2).

Model 2 was defined similar to Model 1, except that we considered the case when the TF can have positive and negative effects on its regulated genes at the same time:

In these two models, the 10-fold cross validation approach was conducted on the training datasets to tune the regularization parameters of the enhanced L 1/2 net, L 1/2 net, and L 1 net. Both penalized parameters for L 1 and ridge regularization in the Elastic net were tuned by the 10-fold cross validation on the two-dimensional parameter surfaces. We repeated the simulations over 100 times and then computed the misclassification error, the sensitivity, and the specificity averagely for each net model on the test datasets.

Table 1 summarizes the simulation results from each regularization net model. In general, our proposed enhanced L 1/2 net model achieved the smallest misclassification errors in Models 1 (9.22%) and 2 (10.76%) compared with the other regularization methods including the old L 1/2 thresholding method (9.85% for Model 1 and 10.83% for Model 2), L 1 net (11.81% for Model 1 and 13.21% for Model 2), and the Elastic net (13.12% for Model 1 and 14.14% for Model 2). Meanwhile, the enhanced L 1/2 net resulted in the highest sensitivity in Model 1 (98.5%) compared with the other methods. Moreover, the enhanced L 1/2 net obtained the best specificity in Model 2 (98.7%) amongst the other approaches. To sum up, the enhanced L 1/2 net outperforms the other three algorithms in terms of prediction accuracy, sensitivity, and specificity.

Table 1

Simulation results of the enhanced L 1/2 net, L 1/2 net, L 1 net, and Elastic net, respectively.

Model	Misclassification errors (%)				Sensitivity (%)				Specificity (%)
	Eh_L 1/2	L 1/2	L 1	Elastic	Eh_L 1/2	L 1/2	L 1	Elastic	Eh_L 1/2	L 1/2	L 1	Elastic
1	9.22	9.85	11.81	13.12	0.985	0.971	0.968	0.873	0.969	0.970	0.962	0.981
	(0.36)	(0.31)	(0.41)	(0.12)	(0.00)	(0.00)	(0.02)	(0.00)	(0.00)	(0.01)	(0.01)	(0.00)
2	10.76	10.83	13.21	14.14	0.939	0.939	0.943	0.835	0.987	0.981	0.987	0.980
	(0.33)	(0.36)	(0.24)	(0.23)	(0.00)	(0.00)	(0.01)	(0.00)	(0.02)	(0.01)	(0.01)	(0.00)

Simulation results (averaged over 100 runs) for comparison of misclassification errors, sensitivity, and specificity used the enhanced L 1/2 net, L 1/2 net, L 1 net, and the Elastic net, respectively. The standard errors are given in parentheses.

4.2. Analysis of Lung Cancer

In this section, we merged the protein-protein interaction (PPI) network (see http://thebiogrid.org/) with a lung cancer (LC) gene-expression dataset [23] to demonstrate the performance of our proposed enhanced L 1/2 net method. The gene-expression dataset contains the expression profiles of 22284 genes for 107 patients, in which 58 had lung cancer. To test the generalization ability of the proposed method, we divided the dataset into the training set (sample size n = 70; 38 LC, 32 non-LC) which covered 2/3 samples of the dataset and the test set (sample size n = 37; 20 LC, 17 non-LC) which covered the other 1/3 specimens of the dataset. The 10-fold cross validation approach was conducted on the training dataset to tune the regularization parameters. By combining the gene-expression data with the PPI network, the final PPI network includes 8619 genes and 28293 edges.

Figures 1 –4 display the solution paths of the four regularization net methods for the LC dataset in one sample run. Here, x-axis displays the values of the running lambda (the running lambda of L 1 penalty in the Elastic net approach), and x-axis at the top (degrees of freedom) means the number of nonzero coefficients of beta. y-axis is the values of the coefficients beta which measure the gene importance. The predictive model builds from the training set and then tests its predictive performance on the test set. The detailed results were represented in Table 2.

Figure 1

The solution paths of the enhanced L 1/2  net for the lung cancer dataset in one sample run.

Figure 2

The solution paths of L 1/2 net for the lung cancer dataset in one sample run.

Figure 3

The solution paths of L 1 net for the lung cancer dataset in one sample run.

Figure 4

The solution paths of the Elastic net for the lung cancer dataset in one sample run.

Table 2

The results of the enhanced L 1/2 net, L 1/2 net, L 1 net, and Elastic net on LC dataset, respectively.

	Selected genes	Connected genes	Connected edges	Cross validation error	Test error
Eh_L 1/2 net	171	54	41	6/70	5/37
L 1/2 net	193	61	47	6/70	6/37
L 1 net	500	150	121	7/70	6/37
Elastic	636	337	510	6/70	6/37

Results of analysis of LC gene expression dataset by four procedures, including the number of genes selected, the number of linked PPI network genes, the number of linked PPI network edges, the CV error, and test errors.

As shown in Table 2, the enhanced L 1/2 net selected the fewest number of genes and edges compared to L 1/2 net, L 1 net, and the Elastic net. Meanwhile, the predictive performance of the enhanced L 1/2 net outperforms the other three regularization net algorithms.

To further evaluate the performance of the enhanced L 1/2 net procedure, we report its capacity of identifying the biomarkers related to lung cancer. NK2 homeobox 1 (Nkx2-1) protein regulates transcription of genes specific for lung. It is used as a biomarker to determine lung cancer in anatomic pathology. It also has a critical role in maintaining lung tumor cells [24, 25]. Epidermal growth factor receptor (EGFR) is known to play a key role in cell proliferation and apoptosis. EGFR overexpression and activity could result in tumor growth and progression [26] and somatic mutations within the tyrosine kinase domain of EGFR, which have been identified in a subset of lung adenocarcinoma [27, 28]. The enhanced L 1/2 net (Figure 5) and L 1/2 net successfully identified these two important biomarkers for LC. However, neither L 1 net nor the Elastic net selected them both.

Figure 5

Subnetworks identified by the enhanced L 1/2 net for lung cancer datasets (only those genes that are linked on the PPI network are plotted).

Except to identify these two significant biomarkers (EGFR and Nkx2-1), the enhanced L 1/2 net also selected several pathways that were associated with lung cancer. For example, one of the subnetworks includes genes involving molecular proliferation (e.g., genes ARF4, EGFR, DCN, BRCA1, and ITIH5). As these gene express significantly and continuously, it promotes lung cancer progression. On the other hand, this group is linked to ENO1. We are unable to get a clear testimony to sustain this relationship by looking at PPI database. However, a recent report [29] has demonstrated that ENO1 is the promising biomarker that may provide more diagnostic efficacy for lung cancer. This link implies a functional relationship and suggests the important role of ENO1 in lung cancer.

All these results reveal that the enhanced L 1/2 net is more reliable than L 1/2 net, L 1 net, and the Elastic net approaches for selecting key markers from high-dimensional genomic data. Another advantage of our proposed method is that it has the ability to recognize novel and potential relationships with biologic significance. It is mentionable that our proposed method is inclined to identify fewer but more informative genes (or edges) than L 1 net and the Elastic net approaches in genomic data and that means the proposed method has allowed the researcher to more easily concentrate on the key targets for functional studies or downstream applications.

5. Conclusions

In biological molecular research, especially for cancer, the analysis of combining biological pathway information with gene-expression data may play an important role to search for new targets for drug design. In this paper, we use the enhanced L 1/2 solver to penalized network-constrained logistic regression model to integrate lung cancer gene-expression with protein-to-protein interaction network. We develop the corresponding coordinate descent algorithm as a novel biomarker selection approach. This algorithm is extremely fast and easy to implement. Both simulation and real genomic data studies showed that the enhanced L 1/2 net is a ranking procedure compared with L 1/2 net (using the old thresholding operator), L 1 net, and the Elastic net in the selection of biomarker and subnetwork.

We successfully identified several important clinical biomarkers and subnetwork that are driving lung cancer. The proposed method has provided new information to investigators in biological studies and can be the efficient tool for identifying cancer related biomarker and subnetwork.

“Sub-networks identified by the L_1 net and the Elastic net for lung cancer datasets (only those genes that are linked on the PPI network are plotted). Nodes colored based on higher (red) to lower (green) coefficients in the model.”