Meta-Analysis Based on Nonconvex Regularization.

The widespread applications of high-throughput sequencing technology have produced a large number of publicly available gene expression datasets. However, due to the gene expression datasets have the characteristics of small sample size, high dimensionality and high noise, the application of biostatistics and machine learning methods to analyze gene expression data is a challenging task, such as the low reproducibility of important biomarkers in different studies. Meta-analysis is an effective approach to deal with these problems, but the current methods have some limitations. In this paper, we propose the meta-analysis based on three nonconvex regularization methods, which are L1/2 regularization (meta-Half), Minimax Concave Penalty regularization (meta-MCP) and Smoothly Clipped Absolute Deviation regularization (meta-SCAD). The three nonconvex regularization methods are effective approaches for variable selection developed in recent years. Through the hierarchical decomposition of coefficients, our methods not only maintain the flexibility of variable selection and improve the efficiency of selecting important biomarkers, but also summarize and synthesize scientific evidence from multiple studies to consider the relationship between different datasets. We give the efficient algorithms and the theoretical property for our methods. Furthermore, we apply our methods to the simulation data and three publicly available lung cancer gene expression datasets, and compare the performance with state-of-the-art methods. Our methods have good performance in simulation studies, and the analysis results on the three publicly available lung cancer gene expression datasets are clinically meaningful. Our methods can also be extended to other areas where datasets are heterogeneous.

Introduction

With the rapid development of biotechnology and its wide applications, many database repositories of high-throughput gene expression data have been created and published. For example, Gene Expression Omnibus (GEO) currently has stored more than 2.76 million samples over 105,000 studies[1]. The gene expression datasets have been widely used in the prediction and diagnosis of diseases, and their application prospects are increasingly promising.

It is desirable to consider variable selection into the analysis of gene expression data due to its small sample size and high dimensionality. Variable selection not only enhances generalization by reducing overfitting, but also enhances interpretability by simplifying the model, i.e., identifying important biomarkers associated with the disease and helping to find the best solution for patients in the treatment process. For a single dataset, there exist many variable selection methods, such as Least Absolute Shrinkage and Selection Operator (LASSO)[2], L1/2 regularization[3-5], Minimax Concave Penalty (MCP)[6], Smoothly Clipped Absolute Deviation (SCAD)[7-9], Group LASSO[10], elastic net[11], Hard Ridge[12], SCAD-L2[13], Complex Harmonic Regularization (CHR) penalty[14] and so on. These methods are effective in discovering important biomarkers in a single dataset. However, it is well known that the analysis of gene expression data is still a challenging task due to high noise and low reproducibility of important biomarkers. There are two main reasons for this challenging task. One is that the decisive biomarkers that regulate the phenotypes are usually very sparse compared to the total number of biomarkers in the entire genome, and their effects are usually weak, therefore, the results of individual studies are not remarkable and difficult to reproduce. The other is that the different experimental datasets may come from inconsistent experimental conditions, sample preparation methods, measurement sensitivities or precision, and also from different study groups, biological sample selections. Therefore, the important genes in some studies may be not remarkable in other studies, which we call the data have the heterogeneity. The data heterogeneity reveals the complexity of gene expression data and significantly obstructs gene expression technology in clinical applications.

Since many genomic databases are publicly available, meta-analysis is an effective approach to address the heterogeneity among different datasets and make full use of different datasets. Meta-analysis is a significant technique for clinical diagnosis, which plays an important role in summarizing and synthesizing scientific evidence from multiple studies. Classic meta-analysis methods, which aggregate the summary statistics from individual datasets to obtain total scores and then evaluate them based on statistical significance of all studies, including p values[15], ranks[16,17], effect sizes[18-21]. Li and Tseng[22] apply Fisher’s method combining p values by summation of log-transformed p values, and the method increases the biological interpretation of meta-analysis results. Similar strategies can be applied to combine effect sizes of Random Effects Model (REM) or Fixed Effects Model (FEM) from individual studies. A comprehensive review of these methods is given in the researches[23-25]. These methods perform well in identifying differentially expressed genes, but they ignore the correlations between the covariates (genes). There are some approaches that attempt to model the preprocessed microarray datasets using latent variable-based models[26-30]. In general, latent variables are not observable in the data, but can be inferred from other observed variables. Huo et al.[31] use latent variable to quantify homogeneous and heterogeneous differentially expressed signals across studies to detect genes that are differentially expressed in only a subset of the combined studies. Rashid et al.[32] utilize a penalized Generalized Linear Mixed Model based on latent variable to select gene signatures and address between-study heterogeneity. These methods provide the potential to pool information across genes, making it possible to more clearly infer which genes are differentially expressed. Compared with the previous classical meta-analysis methods, these methods are more complex, which limit their application in practice. Recently, Zhang et al.[33] set different constant terms for multiple studies in the logistic regression model to measure the heterogeneity of the samples. This method assumes that the same variables in multiple studies should make the same contribution to their corresponding responses. In other words, this method conducts variable selection in an ‘all-in-or-all-out’ fashion. In this paper, we consider that some important genes in some studies are likely to be ineffective in other studies, and it is important to allow such flexibility.

Some researchers propose the bi-level selection methods which consider the coefficients of each variable (gene) from all datasets as a group, and simultaneously shrink these groups and the variables within these groups by the penalty function to study the correlation between variables and identify important genes. Existing bi-level selection methods include composite MCP[34], group Bridge[35] and group exponential LASSO[36], meta-SVM[37] etc. These methods of the aforementioned references generally treat the coefficients of one gene from different datasets as a group, and conduct two levels selection. The first is to determine whether a particular gene is related to the response variable in all datasets, and the second is to determine which dataset contains the identified gene related to the response variable. These methods consider both the heterogeneity and the correlation between the datasets. However, for M independent datasets , each of which contains n samples and p variables, these methods consider to solve the problem which has the dimensional measurement matrix , the dimensional response and the Mp dimensional unknown coefficients , where the superscript T represents the transpose of the vector. Since the gene expression data has the characteristics of small sample size and high-dimensional, these methods greatly increase the variable dimension, so it may increase the difficulty of solving the problem.

Zhou and Zhu[38] propose a new group variable selection method “hierarchical LASSO” that can be used for gene-set selection. The hierarchical LASSO not only removes unimportant groups effectively, but also maintains the flexibility of selecting variables within the group. They also showed that the new method offers the potential for achieving the theoretical “oracle” property. Li et al.[39] propose meta-LASSO for variable selection with high-dimensional meta-analyzed data. The meta-LASSO not only improves the ability to identify important genes with the strength of multiple datasets, but also maintains the flexibility of selection between datasets to consider the data heterogeneity.

For many practical applications, LASSO often cannot find the most sparse solutions (this is extremely important for model selection), and it is inefficient when the errors in data have heavy tail distribution[2]. Zhao and Yu[40] give the Strong Irrepresentable Condition for the model selection consistency of LASSO, and show that to induce sparsity, LASSO shrinks the estimates for the nonzero coefficients too heavily. When Strong Irrepresentable Condition fails, the irrelevant covariates are correlated with the relevant covariates enough to be picked up by LASSO to compensate the over-shrinkage of the nonzero parameters. Therefore, to get universal consistency, some nonconvex regularization methods have been proposed in recent years, such as L1/2 penalty, Minimax Concave Penalty (MCP) and Smoothly Clipped Absolute Deviation (SCAD) penalty etc. These methods achieve both selection consistency and nearly unbiasedness, which make them widely applied in signal/image processing, statistics and machine learning, such as biological feature selection[14,41-44], compressed sensing and low rank matrix completion[8,45-48], sparse signals separation and image inpainting[49,50], and dictionary learning[51] etc.

In this paper, we propose the meta-analysis based on three nonconvex regularization methods (L1/2 regularization, MCP regularization and SCAD regularization), dubbed as meta-Half, meta-MCP and meta-SCAD respectively. Our methods combine the advantages of meta-analysis and the nonconvex regularization methods. We propose the efficient algorithms which apply the nonconvex iterative thresholding algorithms based on approximate message passing (Half-AMP, MCP-AMP and SCAD-AMP)[52,53] to solve our models. Furthermore, we apply our methods to the simulation data and three publicly available lung cancer gene expression datasets, and compare the performance of our methods with other four state-of-the-art methods, which are meta-LASSO, composite MCP, group Bridge and group exponential LASSO. The experiments results show that our methods have favorable performance.

Methodology

In this section, we study the meta-analysis based on the three nonconvex regularization methods (L1/2 regularization, MCP regularization and SCAD regularization).

Consider M independent datasets , each of which contains n samples. Denote and , where the superscript T represents the transpose of the vector, (i = 1, 2, ⋯, n) is ith sample in the mth dataset which contains p variables (genes), and y is the response variable, in this paper, we consider the response variable is a binary phenotype (for example, if the ith sample of the mth dataset is a disease patient, y is 1, and 0 otherwise). The p genes are assumed common in all datasets. We assume the conditional probability that y takes value 1 given the gene expression vector follows the logistic regression model where β is an intercept and is the unknown coefficients for the mth data. Due to heterogeneity between datasets, we allow β and β in (1) to vary with m. We hope to find the true nonzero components of β for each dataset.

Compared with the variable selection of single dataset model, the variable selection of the M datasets models are distinguishing and peculiar. On the one hand, each variable has M coefficients, which belong to the same explanatory variable. Therefore, there is some correlation or similarity, which makes it impossible to make coefficient estimation and variable selection separately, otherwise this correlation will be ignored. On the other hand, the significance of variables is not identical, so we cannot simply synthesize estimation. The penalization methods with meta-analysis make full use of this particularity to study data differences. These methods conduct variable selection by maximizing, where ℓ(β, β) is the log-likelihood for the mth dataset and has the following form P is a penalty function and λ is the regularization parameter that controls the complexity of the machine.

In this paper, we focus on the three nonconvex regularization methods (L1/2 regularization, MCP regularization and SCAD regularization), and through the hierarchical decomposition of coefficients that maintain the flexibility of variable selection as well incorporate the relationship between different datasets. We consider the following hierarchical reparameterization:

The parameter h is the effect of the jth gene, and the different m for ξ reflects the different effects of the jth gene among M datasets. If h = 0, then , this indicates that the jth gene is not significant in all M datasets. If h ≠ 0, then whether the β is equal to 0 depends on whether ξ is equal to 0. Since the M datasets may have heterogeneity (the M datasets may come from inconsistent experimental conditions, sample preparation methods, measurement sensitivities or precision, and also from different study groups, biological sample selections.), then one gene is important in some datasets may be not remarkable in other datasets. Through ξ contral β to keep the selection flexibility among M datasets. If the M datasets have no heterogeneity, then h = β for m = 1, ⋯, M defined in (2) and ξ = 1 for all j and m. With reparameterization (3), we propose a meta-analysis method based on nonconvex regularization. Our method selects important genes by solving where ℓ(β, , ξ) is the likelihood function and has the following form , , , , and  ⋅ ξ means the element-wise product. P( ⋅ ) is a nonconvex penalty function. In this paper, considering the three nonconvex penalty function, the L1/2 penalty, the MCP penalty and the SCAD penalty. The L1/2 penalty function is . MCP penalty function has the following form where . The SCAD penalty function has the following form we call these three nonconvex penalties for the methods (4) as “meta-Half”, “meta-MCP” and “meta-SCAD”, respectively.

Algorithm

In this section, we give the efficient algorithms (Algorithm 1) to solve our models. Note that we can assume that the mean of the predictor variable is zero (through the location transformation). (4) can be decomposed into two nonconvex problems, each of which views or ξ as fixed. We propose to iteratively solve β0, , and ξ in (4). First, we fix β0 and ξ in (4) to maximize . We next fix β0 and to maximize ξ. Finally, we maximize over β0 by fixing and ξ. Iterate these steps until the algorithm converges. Since at each step, the value of the objective function (4) decreases, the solution is guaranteed to converge. Specifically, the algorithm is described as follows

The iterative optimization algorithm for solving our meta-analysis based on nonconvex regularization models.

Step 3 and step 4 are general nonconvex regularization problem. We[52,53] propose the nonconvex iterative thresholding algorithms based on approximate message passing (Half-AMP, MCP-AMP and SCAD-AMP) to solve the nonconvex regularization problem, and verified the effectiveness of the algorithms through theoretical analysis and experiment. In this paper, for the two problems in step 3 and step 4, we apply the Half-AMP algorithm, the MCP-AMP algorithm and the SCAD-AMP algorithm to solve the meta-Half, meta-MCP and meta-SCAD respectively.

In order to solve the above two problems in step 3 and step 4, we first consider the solution of the traditional logistic regression model. Here we omit the intercept term (in fact, just rewrite the input variable as ), the logistic regression can be expressed as the following optimization problem Differentiating ℓ(β) with respect to β, we can get where To find the optimal solution of equation (6), we let , and use the Newton-Raphson iteration algorithm which requires computing the second derivative, the Hessian matrix has the following form Hence, given the current estimated value βold of β, the new estimated value βnew is updated as following where the value of the derivative (and second derivative) is calculated at the point βold. The equation (9) can be expressed by matrix form. Let be a N × P matrix, where the i-th row is ; W is a diagonal matrix, and the elements on the diagonal Let , then the formulas (7) and (8) can be expressed as Therefore, Newton-Raphson iteration (9) can be expressed as where It can be seen that each Newton-Raphson iteration actually solves the weighted least squares problem as follows

Based on the solution process of the traditional logistic regression model, a similar iterative algorithm can be used to solve the logistic regression with nonconvex penalties problem, and only a slight deformation of the formula (14) is needed to obtain the iterative algorithm. where P( ⋅ ; ⋅ ) is the nonconvex penalty function. It is easy to see that the minimization problem (15) is equivalent to the maximization problem (2). The minimization problem (15) can be solved by the nonconvex iterative thresholding algorithms based on approximate message passing[52,53] (which are based on linear regression  = β,  ∈ R,  ∈ R). The algorithms are according to the following iteration: where represents a measure of indeterminacy of the measurement system, in this paper, considering the case δ is fixed for N → ∞. For a vector u = (u1, u2, . . . , u), /N, . η is the thresholding function, in this paper, η represents the Half thresholding function, the MCP thresholding function and the SCAD thresholding function, respectively. The Half thresholding function is where The MCP thresholding function is where where sign(u) is sign function, The SCAD thresholding function is

Theoretical Properties

In this section, we study the theoretical properties of the meta-Half method. The meta-Half has the following uniform form there are two tuning parameters λ and λ in (19), we first show that the two tuning parameters can be simplified into one. Specifically, let λ = λλ, we can show that (19) is equivalent to

Lemma 1.

If is a local maximizer of (19). Then there exists a local maximizer of (20) such that and . Vice versa.

The proof is in the Supplementary. This lemma indicates that although (19) and (20) may provide different h and ξ, the final fitted models from them are the same. Therefore, we only need to tune one parameter λ = λλ other than tune λ and λ separately in practice.

We then show that (20) can also be written in an equivalent form using the original regression coefficients β.

Lemma 2.

Suppose is a local maximizer of (20), for j = 1, 2, ⋯, p, let , and ,

(a) If , then ;

(b) If , then and , .

The proof is in the Supplementary.

Theorem 1.

If is a local maximizer of (20), then with , is a local maximizer ofwhere . On the other hand, if is a solution of (21), then is a solution of (20), where , , where , , and .

The proof is in the Supplementary. If we regard one gene’s effects among all datasets as a “group”, then (21) imposes an L1 penalty on each group and a square root penalty on individual elements within a group. The following theorem shows the theoretical properties of the meta-Half.

Theorem 2.

The meta-Half method possesses sparsity, unbiasedness and oracle properties.

The proof is in the Supplementary.

Experiments

In this section, we analyze the performance of our methods (meta-Half, meta-MCP and meta-SCAD) by simulation and real-data analysis. We compare these three methods with other four methods, which are meta-LASSO, composite MCP, group Bridge and group exponential LASSO. The codes of our methods are available at GitHub (https://github.com/zhhui019/meta-nonconvex). The meta-LASSO is implemented by Li et al.[39]. The composite MCP, the group Bridge and the group exponential LASSO are implemented by Patrick Breheny and Yaohui Zeng’s R package “grpreg”.

Simulations

Simulation studies are performed to compare the performance of the proposed meta-Half, meta-MCP and meta SCAD with the meta-LASSO, composite MCP, group Bridge and group exponential LASSO.

Generate simulated data

In this simulation, we use the normal distribution to generate the gene expression (m = 1, 2, ⋯, M; i = 1, 2, ⋯, n) with M = 10 datasets, each dataset contains n = 50 samples, and each sample contains p = 1, 000 genes. The response y is generated from a logistic model where and we suppose that the intercept term We let simulate possible data heterogeneity, for m = 1, 2 ⋯, M; j = 1, 2 ⋯, 10, α are generated from N(3, 0. 52) and θ are generated from Bernoulli(π0), for m = 1, 2 ⋯, M; j = 11, 12 ⋯, 1000, let . This means that the first 10 genes of each dataset are important to the response with probability π0. The value α demonstrates whether the jth gene is important in the mth dataset, and the value θ demonstrates different levels of heterogeneity among different datasets, in this simulation, considering π0 = 0.9, 0.5, 0.2 to represent the low, medium and high heterogeneity. We run 30 replicates and report the average measurement.

For the all methods, the tuning parameters are selected by minimizing the BIC: where is the estimated coefficients in the mth dataset, λ is the tuning parameter, S is the number of non-zero elements of ,  is the log-likelihood for the mth dataset and has the form (5).

Analysis of simulation

The variable selection performance of the seven methods is evaluated using the selection sensitivity, specificity and accuracy of coefficient β. The sensitivity is the proportion of non-zero ’s that are correctly estimated as non-zero, the specificity is the proportion of zero ’s that are correctly estimated as zero and the accuracy is the proportion of ’s that are correctly estimated.

The simulation results are summarized in Table 1 (The variable selection performance of the seven methods are evaluated using the selection sensitivity, specificity and accuracy of coefficient β). Table 1 shows that the specificity and accuracy of the coefficients β of all seven methods are similar. The sensitivity trend of coefficient β for all seven methods with the varying levels of heterogeneity is shown in Fig. 1. Figure 1 shows that the sensitivity (the proportion of non-zero ’s that are correctly estimated as non-zero) of the composite MCP, the group Bridge and the group exponential LASSO dramatically decreases as π0 increases, while the sensitivity of meta-Half, meta-MCP, meta-SCAD and meta-LASSO remains above 0.9 for π0 = 0.2, 0.5, 0.9. When π = 0.2, the sensitivity of meta-Half, meta-MCP and meta-SCAD are 0.9693, 0.9651 and 0.9738, respectively, which are significantly higher than other methods. This result shows that our proposed meta-Half, meta-MCP and meta-SCAD have the superior performance when data heterogeneity is strong (π0 is small). With the weakening of data heterogeneity(π = 0.5, 0.9), the performance of the four meta methods (meta-Half, meta-MCP, meta-SCAD and meta-LASSO) tends to be comparable. The specificity and accuracy of the coefficients for all seven methods are similar.

Table 1

The sensitivity, specificity and accuracy of coefficient β of the seven methods: presented values are the mean (standard error).

		π = 0.2	π = 0.5	π = 0.9
meta-Half	Sensitivity	0.9693 (1.70E − 03)	0.9215 (4.97E − 03)	0.9229 (1.30E − 03)
	Specificity	0.9862 (2.26E − 04)	0.9903 (6.36E − 05)	0.9837 (1.40E − 03)
	Accuracy	0.9861 (2.24E − 04)	0.9901 (6.50E − 05)	0.9835 (1.41E − 03)
meta-MCP	Sensitivity	0.9651 (1.90E − 03)	0.9362 (3.90E − 02)	0.9205 (2.70E − 03)
	Specificity	0.9884 (4.18E − 05)	0.9840 (2.02E − 05)	0.9846 (1.15E − 02)
	Accuracy	0.9883 (4.16E − 05)	0.9838 (2.05E − 05)	0.9840 (1.13E − 02)
meta-SCAD	Sensitivity	0.9738 (1.60E − 03)	0.9306 (2.07E − 03)	0.9392 (2.50E − 03)
	Specificity	0.9903 (1.44E − 05)	0.9853 (1.20E − 04)	0.9505 (4.20E − 03)
	Accuracy	0.9903 (1.42E − 05)	0.9850 (1.44E − 05)	0.9504 (4.10E − 03)
meta-LASSO	Sensitivity	0.9065 (8.42E − 02)	0.9217 (6.60E − 02)	0.9425 (6.48E − 02)
	Specificity	0.9710 (2.72E − 03)	0.9869 (3.07E − 03)	0.9940 (1.71E − 03)
	Accuracy	0.9708 (2.79E − 03)	0.9866 (3.02E − 03)	0.9935 (1.69E − 03)
composite MCP	Sensitivity	0.8454 (1.46E − 01)	0.5428 (1.23E − 01)	0.3167 (7.60E − 02)
	Specificity	0.9988 (6.02E − 04)	0.9992 (4.85E − 04)	0.9984 (7.02E − 04)
	Accuracy	0.9985 (6.12E − 04)	0.9969 (1.02E − 03)	0.9922 (1.18E − 03)
group Bridge	Sensitivity	0.8734 (7.77E − 02)	0.6856 (1.11E − 01)	0.2842 (6.27E − 02)
	Specificity	0.9997 (2.84E − 04)	0.9999 (1.05E − 04)	0.9999 (3.49E − 05)
	Accuracy	0.9994 (3.42E − 04)	0.9983 (6.08E − 04)	0.9934 (6.74E − 04)
group exponential LASSO	Sensitivity	0.8809 (8.70E − 02)	0.7315 (1.67E − 01)	0.4661 (2.29E − 01)
	Specificity	0.9984 (1.10E − 03)	0.9981 (9.33E − 04)	0.9976 (1.33E − 03)
	Accuracy	0.9981 (1.09E − 03)	0.9967 (8.01E − 04)	0.9928 (1.37E − 03)

Figure 1

The sensitivity trend of coefficient β for all seven methods with the varying levels of heterogeneity.

Real-Data analysis

In this section, we apply our methods (meta-Half, meta-MCP and meta-SCAD) to three publicly available lung cancer gene expression datasets, and compare our three methods with other four methods including meta-LASSO, composite MCP, group Bridge and group exponential LASSO.

Lung cancer datasets

The three publicly available lung cancer microarray datasets come from disparate platforms and can be download from GEO (https://www.ncbi.nlm.nih.gov/gds/). The three datasets are described as follows:

GSE10072 dataset. The dataset is gene expression signature of cigarette smoking, it contains 107 final expression samples from 58 tumors and 49 non-tumor tissues from 20 never smokers, 26 former smokers, and 28 current smokers, each sample has 22283 genes. The original gene expression data is provided by Landi et al.[54].

GSE19188 dataset. The dataset is expression data for early stage non-small-cell lung cancer (NSCLC), it contains 156 samples from 91 tumor tissues and 65 adjacent normal lung tissue samples, each sample has 54675 genes. The more information can be found in Hou et al.[55].

GSE19804 dataset. The dataset is non-smoking female lung cancer in Taiwan, it contains 120 samples from 60 tumors and 60 normal tumor tissues, each sample has 54675 genes. The more information can be found in Lu et al.[56].

Each dataset is divided into two parts, about 70 percent of the datasets as training samples and the other 30 percent as testing samples. Table 2 lists the details of the three datasets.

Table 2

The description of three publicly available lung cancer gene expression datasets.

Dataset	No. of Probs	Classes (Class 0/Class 1)	No. of samples (Class 0/Class 1)	Affymetrix Platform
GSE10072	22284	Normal/ Lung Cancer	107 (49/58)	U133A
GSE19188	54675	Normal/ Lung Cancer	156 (65/91)	U133 Plus 2.0
GSE19804	54676	Normal/ Lung Cancer	120 (60/60)	U133 Plus 2.0

The original Affymetrix data was first normalized and log-transformed by a robust multi-array average (RMA) method[57]. After that, downloading and installing the appropriate custom chip definition files (CDFs) packages according to the type of microarray platform. The CDF package is necessary for probe annotation for Affymetrix data. The probes of the normalized data can be successfully mapped to Entrez Gene IDs by annotation packages in Bioconductor[58]. If multiple probes match a single Entrez ID, we calculated the median of values of those probes as the expression value for this gene.

We extract common genes from the three gene expression datasets as the merged set of genes. There are 13515 common genes in three datasets and our analysis is based on those 13515 genes. We use a random partition in three lung cancer datasets, and apply aforementioned seven methods to select important genes, with the optimal tuning parameters chosen by the BIC as discussed above. We repeat this procedure 30 times and report the average measurement and standard error.

Evaluating the classification performance

Table 3 demonstrates the prediction performance of the seven methods in three lung cancer datasets. The sensitivity, specificity and accuracy of training and testing predictions for all seven methods are shown in Fig. 2.

Table 3

Performance comparisons of different methods in three lung cancer datasets. Presented values are the average (standard error).

Methods	Training data			Testing data
	Accuracy	Sensitivity	Specificity	Accuracy	Sensitivity	Specificity
meta-Half	0.9766 (2.69E-02)	0.9673 (9.10E-03)	0.9903 (1.54E-05)	0.9449 (2.16E-02)	0.9464 (7.93E-03)	0.9437 (1.66E-05)
meta-MCP	0.9768 (5.37E-04)	0.9677 (1.70E-03)	0.9903 (8.56E-05)	0.9452 (1.99E-03)	0.9466 (6.94E-03)	0.9439 (1.35E-04)
meta-SCAD	0.9727 (3.13E-03)	0.9608 (1.77E-03)	0.9903 (2.19E-02)	0.9528 (4.64E-03)	0.9464 (4.80E-03)	0.9577 (2.50E-02)
meta-LASSO	0.9309 (1.01E-02)	0.8722 (1.87E-02)	0.9994 (2.15E-03)	0.8953 (2.11E-02)	0.8291 (3.58E-02)	0.9792 (1.25E-02)
composite MCP	0.9353 (1.45E-02)	0.9221 (2.05E-02)	0.9519 (1.85E-02)	0.8656 (1.94E-02)	0.8283 (2.23E-02)	0.9060 (2.85E-02)
group Bridge	0.6410 (1.88E-02)	0.4039 (2.59E-02)	0.9508 (2.23E-02)	0.6255 (2.06E-02)	0.3240 (3.15E-02)	0.9317 (2.98E-02)
group exponential Lasso	0.9385 (9.78E-03)	0.9155 (1.64E-02)	0.9655 (1.48E-02)	0.8942 (1.85E-02)	0.8432 (2.47E-02)	0.9589 (2.05E-02)

Figure 2

Training and testing prediction performance of different methods on lung cancer datasets. (a) Training. (b) Testing.

As shown in Table 3 and Fig. 2, for the training dataset and testing dataset, the sensitivity and accuracy of meta-Half, meta-MCP, meta-SCAD are consistently higher than the other four methods, and the specificity of all methods are similar. This result shows that our three methods are more effectively distinguish whether an individual is a disease patient compared to the other four methods. Therefore, our three methods have superior performance than the other four methods in the prediction and diagnosis of diseases.

Analysis of the selected genes

Table 4 gives the names of genes selected in each dataset. We focus on the gene WIF1 which is bolded in the Table 4. WIF1, a secreted Wnt antagonist, is a downstream gene of the Wnt/β-catenin pathway, which exerts inhibition through direct binding to Wnt proteins[59]. WIF1 was found to be silenced by methylation in various human carcinomas including lung[60], oral[61], nasopharyngeal[62], esophageal[63], breast[64] and colon cancer[65] etc.

Table 4

Gene selections of seven methods in three lung cancer datasets.

	GSE10072			GSE19188			GSE19804
meta-Half	CXCL13	MMP12	COL11A1	CXCL13	MMP12	COL11A1	CXCL13	SPINK1	AGER
	TOX3	SPINK1	FCN3	TOX3	SPINK1	FCN3	SPP1	COL10A1	GPM6A
	SPP1	COL10A1	WIF1	SPP1	COL10A1	WIF1	MMP12	TOX3
	GPM6A	AGER		GPM6A	AGER
meta-MCP	CXCL13	SPP1	COL11A1	CXCL13	SPP1	COL11A1	CXCL13	SPP1
	AGER	MMP12	FCN3	AGER	MMP12	FCN3	AGER	MMP12
	PPAP2C	COL10A1	WIF1	PPAP2C	COL10A1	WIF1	PPAP2C	COL10A1
	GPM6A	TOX3		GPM6A	TOX3		GPM6A	TOX3
	TOP2A	TMEM100		TOP2A	TMEM100		TOP2A	TMEM100
meta-SCAD	CXCL13	MMP12	COL11A1	CXCL13	MMP12	COL11A1	CXCL13	MMP12	COL11A1
	CYP4B1	SPINK1	FCN3	CYP4B1	SPINK1	FCN3	CYP4B1	SPINK1	FCN3
	SPP1	AGER	WIF1	SPP1	AGER	WIF1	SPP1	AGER
meta-LASSO	PPBP	SFTPC	SPP1	PPBP	SFTPC	SPP1	PPBP	SFTPC	SPP1
	CLDN10	AKR1B10	SFTPD	CLDN10	AKR1B10	SFTPD	CLDN10	AKR1B10	SFTPD
	UPK3B	APOLD1	XIST	UPK3B	APOLD1	XIST	UPK3B	APOLD1	XIST
	SPINK1	COL10A1	WIF1	SPINK1	COL10A1	WIF1	SPINK1	COL10A1	WIF1
	HLA-DQA1 /// LOC100509457			HLA-DQA1 /// LOC100509457			HLA-DQA1 /// LOC100509457
composite MCP	SOSTDC1	COL11A1		SYNE1			SOSTDC1	COL11A1
group Bridge	P2RY14			P2RY14	ATP1A2		P2RY14
group	GDF10	FABP4	COL11A1	GDF10	FABP4	COL11A1	GDF10	FABP4	COL11A1
exponential
LASSO

As shown in Table 4, our three methods (meta-Half, meta-MCP and meta-SCAD) all select gene WIF1 on both datasets GSE10072 and GSE19188, but not select gene WIF1 on dataset GSE19804. As we known that GSE10072 dataset is gene expression signature of cigarette smoking, GSE19188 dataset is expression data for early stage non-small-cell lung cancer (NSCLC) and GSE19804 dataset is non-smoking female lung cancer in Taiwan. Huang et al. showed that WIF1 is significantly associated with the smoking behavior in NSCLC patients[66]. It shows that our three methods can more realistically identify the important biomarkers from different datasets which have heterogeneity. The meta-LASSO selects gene WIF1 in all three lung cancer datasets, and the genes selected by meta-LASSO in the three lung cancer datasets are the same. The other three methods (composite MCP, group Bridge and group exponential LASSO) cannot select gene WIF1 in all three lung cancer datasets. Therefore, our three methods are superior to the other four methods when applied in the heterogeneity datasets.

The number of genes selected by meta-Half, meta-MCP, meta-SCAD and meta-LASSO are 11, 13, 9 and 13 respectively. Figure 3 shows the overlap of commonly selected genes across the four different methods (meta-Half, meta-MCP, meta-SCAD, meta-LASSO) in three lung cancer datasets. The other three methods (composite MCP, group Bridge and group exponential LASSO) select fewer genes, so we don’t show the genes they selected in Fig. 3. As shown in Fig. 3(a), for the datasets GSE10072 and GSE19188, seven common genes are selected by meta-Half, meta-MCP and meta-SCAD, which are CXCL13, COL11A1, SPP1, MMP12, AGER, WIF1 and FCN3. Two common genes are selected by meta-Half, meta-MCP, meta-SCAD and meta-LASSO, which are SPP1 and WIF1. Figure 3(b) shows that for dataset GSE19804, four common genes selected by meta-Half, meta-MCP and meta-SCAD are CXCL13, SPP1, MMP12 and AGER. One common genes are selected by meta-Half, meta-MCP, meta-SCAD and meta-LASSO, which is SPP1. More unique non-overlapping sets of genes are selected by our three methods and meta-LASSO. In addition, some of the aforementioned genes have been reported in the literature. COL11A1 is collagen type XI alpha 1 chain. The over-expression of COL11A1 reportedly correlates with lymph node metastasis and poor prognosis in non-small cell lung cancer and ovarian cancer[67]. Zhang et al. suggest that SPP1 and AGER are risk factors for lung adenocarcinoma, and these two genes may be utilized in the prognostic evaluation of patients with lung adenocarcinoma[68]. The advanced glycosylation end-product specific receptor (AGER) belongs to the immunoglobulin superfamily, whose abnormal expression has been detected in lung cancer[69]. MMP12 is matrix metallopeptidase 12 and may play a role in aneurysm formation and mutations in this gene are associated with lung function and chronic obstructive pulmonary disease (COPD)[70]. WIF1 was found to be silenced by methylation in lung[60]. Lea et al. shows that the Ficolin-3, encoded by the FCN3 gene and expressed in the lung and liver, is a recognition molecule in the lectin pathway of the complement system[71]. The aforementioned genes CXCL13, MMP12, AGER and FCN3 are only selected by our three methods, and the gene COL11A1 is selected by our three methods and group exponential LASSO.

Figure 3

Overlap of commonly selected genes across the different methods in lung cancer datasets. (a) GSE10072 and GSE19188. (b) GSE19804.

To make it easier to demonstrate the interplay between the selected genes from the different methods, we construct a network of interactions among the genes using the cBioPortal[72,73]. Figures 4, 5 and 6 show the interactive network of the genes selected by our three methods in three lung cancer datasets. Most of the genes selected by our three methods are linked to the frequently altered neighbor genes from the TCGA lung adenocarcinoma dataset. The expression of SPP1 is controlled by TP53. TP53 is tumor protein p53, this gene encodes a tumor suppressor protein containing transcriptional activation, DNA binding, and oligomerization domains. Mutations in this gene are associated with a variety of human cancers[74]. MMP12 and TOP2A are targeted by certain cancer drugs, and are only selected by our three methods.

Figure 4

Network view of the genes selected from meta-Half in lung cancer datasets. The genes corresponding to the selected variables are highlighted by a thicker black outline. The rest of the nodes correspond to the genes that are frequently altered and are known to interact with the highlighted genes (based on publicly available interaction data). The nodes are gradient color-coded according to the alteration frequency based on microarray data derived from the TCGA lung cancer dataset via cBioPortal. (a) GSE10072 and GSE19188. (b) GSE19804.

Figure 5

Network view of the genes selected from meta-MCP in lung cancer datasets. (a) GSE10072 and GSE19188. (b) GSE19804.

Figure 6

Network view of the genes selected from meta-SCAD in lung cancer datasets. (a) GSE10072 and GSE19188. (b) GSE19804.

In this part, we analyze the genes selected by the four methods (meta-Half, meta-MCP, meta-SCAD and meta-LASSO) in three lung cancer datasets. According to the network of interactions between genes, among the genes selected by our three methods, we find that some genes are connected to other frequently altered genes in publicly available datasets, and some genes are targeted by certain cancer drugs. Some functions may also need to be verified in the future. Results demonstrate that our three methods have good performance in the high-dimensionality gene expression data with heterogeneity.

Conclusion

With the rapid development of biotechnology and its wide applications, a large number of publicly available gene expression datasets have been produced. However, due to the gene expression datasets have the characteristics of small sample size, high dimensionality and high noise, the application of biostatistics and machine learning methods to analyze gene expression data is a challenging task, such as the low reproducibility of important biomarkers in different studies. The low reproducibility of important biomarkers is mainly caused by the heterogeneity of the different datasets. These problems reveal the complexity of gene expression data and significantly obstruct biotechnology in clinical applications. Meta-analysis is an effective approach to deal with these problems. It plays an important role in summarizing and synthesizing scientific evidence from multiple studies, and provides a more comprehensive understanding of the biological systems, but the current methods have some limitations. The nonconvex regularization method is an effective approach for variable selection developed in recent years. In this paper, we combine the advantages of meta-analysis and the nonconvex regularization method, and propose three novel methods, dubbed as meta-Half, meta-MCP and meta-SCAD, respectively. Through the hierarchical decomposition of coefficients, our methods not only consider the data heterogeneity to maintain the flexibility in selecting variables on different datasets, but also consider the correlation between multiple datasets to improve the ability of identifying important biomarkers. We give the efficient algorithms which apply the nonconvex iterative thresholding algorithms based on approximate message passing (Half-AMP, MCP-AMP and SCAD-AMP) to solve our models and study the theoretical property of meta-Half. The theoretical property analysis of MCP-AMP and SCAD-AMP are the future work. We prove meta-Half possesses sparsity, unbiasedness and oracle properties. Furthermore, we apply our methods to the simulation data and three publicly available lung cancer gene expression datasets, and compare the performance of our methods with other four methods, which are meta-LASSO, composite MCP, group Bridge and group exponential LASSO. Simulation studies demonstrate our methods have the superior performance when data heterogeneity is strong. In the three publicly available lung cancer gene expression datasets, the analysis results show that our three methods have good performance in the gene expression data of small sample size and high dimensionality from different sources (heterogeneity), and the selected important biomarkers have clinical significance. Our methods can also be extended to other areas where datasets are heterogeneous.

LaTeX Supplementary File.