Large-scale gene expression analysis reveals robust gene signatures for prognosis prediction in lung adenocarcinoma

Data collection and pre-processing

Gene expression datasets were obtained from the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo) and ArrayExpress (http://www.ebi.ac.uk/arrayexpress/) databases. Prior to downstream analysis, we firstly mapped the array probes to the respective gene symbol by using the array annotations. To calculate the gene expression more conveniently, we used the average expression values of genes matching multiple probes.

Binarization of gene expression levels from multiple datasets

The first seven datasets used in this study was merged by Lim’s merging method to remove batch effect as they were produced by the same microarray platform (Table 1). For each gene of the merged dataset and the 3 additional datasets, the expression values were binarized as high or low expression when the expression values higher or lower than its corresponding median, respectively. Based on the binarized gene expression pattern for each gene and each sample, we then merged the seven datasets and three addition datasets.

Overrepresentation enrichment analysis (ORA)

Overrepresentation enrichment analysis, which used hypergeometric test, was also implemented at WEB-based Gene Set Analysis Toolkit (WebGestalt) (Wang et al., 2017). The Reactome pathways were selected as the functional database (Fabregat et al., 2018). We chose 0.05 as the threshold of the p-value for significant pathways.

Gene set enrichment analysis

The gene set enrichment analysis was implemented in javaGSEA (Subramanian et al., 2005) (version 3.0). The database with GMT files was customized by NCI-PID pathways (Schaefer et al., 2009) selected from all canonical pathways. The genes were pre-ranked based on the Z statistic in Cox model. 10,000 permutations were used to calculate the enrichment significance.

Cox-regression based survival analysis

Cox-regression model was used to evaluate the differences of overall survival between patients from two conditions. This analysis was implemented in R programming software (R Core Team, 2018) with the survdiff function. To visualize the overall survival for each group, we used Kaplan–Meier curves to estimate the survival probability. The hazard.ratio function in survcomp package (Haibe-Kains et al., 2008) was used to calculate the hazard ratios and corresponding p-values. The risk score for each patient was predicted by the Cox model with “linear predictor” type based on the 25 genes selected by MMPC algorithm (Brown, Tsamardinos & Aliferis, 2004), which was implemented in predict.coxph function.

Summary of enrolled datasets for discovery

A total of 10 non-small cell lung cancer (NSCLC) gene expression datasets were collected from Gene Expression Omnibus (GEO) or ArrayExpress database. Tumor samples should be characterized by primary lung adenocarcinoma histology, and with overall survival. Notably, 309 tumor samples from seven datasets, including GSE10245 (Kuner et al., 2009), GSE10445 (Broet et al., 2009), GSE19188 (Hou et al., 2010), GSE28571 (Micke et al., 2011), GSE31210 (Okayama et al., 2012), GSE33356 (Lu et al., 2011), and GSE50081 (Der et al., 2014), were produced by the same microarray platform (Affymetrix Human Genome U133 Plus 2.0 Array), which were merged and normalized by Lim et al. (2018). In addition, another three datasets, GSE68465 (Director’s Challenge Consortium for the Molecular Classification of Lung et al., 2008), GSE67639 (Roepman et al., 2009), and GSE13213 (Tomida et al., 2009), were also incorporated in the present study. Finally, a total of 1,308 LUAD cases were collected for further analysis, 741 (56.65%) of whom were dead (Table 1).

Identification of prognostic genes by meta-analysis-based Cox regression model

To robustly identify the prognostic genes associated with overall survival of lung adenocarcinoma, we integrated the ten gene expression datasets, and discretized the normalized expression value for each gene as high and low expression status within each dataset, which could avoid the batch effect by different platforms. Cox proportional hazard regression analysis was then performed on the discretized expression status for each gene. Given a stringent threshold at BH-adjusted p-value < 0.01, we successfully identified 42 genes significantly associated with LUAD overall survival, including 21 positively and 21 reversely correlated genes (Fig. 1A).

To further investigate functional roles of the prognostic genes, we performed overrepresentation enrichment analysis (ORA) on these genes. We identified seven pathways significantly enriched by the prognostic genes (Fig. 1B, p-value < 0.05). Remarkably, the cell cycle genes, such as CHEK1, PCNA, RRM2, BUB1B, and CDC6, were reversely correlated with patients’ overall survival, which were significantly enriched in pathways, such as G1/S-specific transcription, G1/S transition, cell cycle checkpoints, and cell cycle. Moreover, GALNT3 and GALNT12, also reversely correlated with overall survival, were involved in O-linked glycosylation of mucins, indicating that O-linked glycosylation of mucins played key roles in LUAD progression. In addition, we also identified two prognostic genes, ABAT and STX1A, which participated in GABA synthesis, release, reuptake and degradation. Notably, cell cycle (Sithanandam et al., 2003), O-linked glycosylation of mucins (Hauselmann & Borsig, 2014), and GABA synthesis, release, reuptake, and degradation (Al-Wadei et al., 2012) have been reported to be involved in tumorigenesis or tumor progression. The results based on the enrichment analysis indicated that pathways such as cell cycle, O-linked glycosylation of mucins, and GABA synthesis, release, reuptake and degradation were the hallmarks of tumor progression and short overall survival.

Identification of prognostic pathways by GSEA

To identify the prognostic pathways for LUAD patients, we ranked the genes based on their significance levels by Cox regression-based meta-analysis. The gene set enrichment analysis was then performed on the ranked gene set. Given the stringent thresholds (FDR < 0.05 for pathways, and log-rank test p-value < 0.05 for pathway genes in core enrichment), we identified the Aurora-A pathway, the Aurora-B pathway, and the FOXM1 transcription factor network as prognostic pathways in LUAD (Fig. 2).

The Aurora-A and Aurora-B pathways were responsible for G2-M transition in cell cycle (Sithanandam et al., 2003), and the expression levels of two key kinases, Aurora-A and Aurora-B, were significantly higher in high-risk group than low-risk group (Figs. 2A–2C). As FOXM1 is a transcription factor, which was a famous oncogene (Gartel, 2017), its target genes, such as CCNA2, CCNB1, CCNB2, CCNE1, TGFA, BIRC5, CDK2, CENPF, CENPA, and AURKB, were closely associated with overall survival of LUAD patients. Particularly, the transcription factor FOXM1, overexpression of which could significantly shorten the overall survival of LUAD patients, was also involved in G2-M transition (Fig. 2C). The result suggested that hyperactive G2-M transition in cell cycle was an indicator of poor prognosis in LUAD.

Development of a gene expression signature-based prognostic model in LUAD

As we described above, the univariate Cox proportional hazard regression analysis successfully identified 42 prognostic genes. To further select signatures used for multivariable Cox regression model, we employed MMPC algorithm, which is a constraint based feature selection algorithm (Brown, Tsamardinos & Aliferis, 2004). We then selected 25 genes from the 42 prognostic genes, including ABAT, BCAR3, CTSF, DEAF1, ENC1, ETV5, FAM117A, FZD2, GALNT12, GALNT3, GJB3, KDM6A, KYNU, PCNA, PFKP, PLEK2, RASGRP2, SERPIND1, SGSH, TLE1, TMEM38B, TMEM57, TRIM45, USP47, and VWA1, at the threshold of p-value < 0.1 for MMPC algorithm. Finally, we built a multivariable Cox regression model on the 25 genes for overall survival prediction (Table 2). Based on the multivariable Cox regression model, risk score for each patient in the training set was calculated, and the 1,308 patients were classified into high- and low-risk groups. Kaplan–Meier curves showed that patients in the high-risk group had significantly shorter overall survival than those in the low-risk group (log-rank test P < 0.0001) (Fig. 3A).

In addition, we also investigated the performance of our stratification in specific stage, gender, and age group of LUAD in the training set. As no samples were stratified into TNM stage IV group in the training set, we only focused on the performance of the model in another three stages (I, II, and III). The overall survival difference between high- and low-risk groups in training set was observed in TNM stages I, II, and III, male/female, and old/young groups (Figs. 3B–3H, log-rank test, P < 0.005), in accordance with the performance in all samples. These results indicated that our stratification in training set was independent on TNM stages, gender, and age.

Evaluation of the gene expression signature-based prognostic model in the validation sets

To evaluate the performance of the prognostic model in independent datasets, we collected two LUAD gene expression datasets, TCGA-LUAD (The Cancer Genome Atlas-lung adenocarcinoma, n = 502) (Cancer Genome Atlas Research N, 2014) and GSE37745 (n = 106) (Botling et al., 2013). The risk scores for the patients from validation sets based on the prognostic model were calculated. By using the same model and criteria, patients in the validation sets were classified into high-risk and low-risk groups. Similar with that in training set, the overall survival of the patients in high-risk group was significantly worse than that of low-risk group patients in the two validation sets (P < 0.001) (Figs. 4A–4B). Notably, the stratification still showed significant predictive ability in overall survival by adjusting the cofactors including age, gender, smoking status, tumor stage in TCGA cohort (P < 0.0001, Table 3). The distribution of the risk score, overall survival status along with the corresponding expression profiles of the 25 prognostic genes from two validation sets were showed in Figs. 4C–4D, which were ranked according to the risk score value. The 25 prognostic genes were significantly differentially expressed between the two risk groups (P < 0.05). The results indicated that the 25-gene signature based prognostic model showed high and robust performance in both training and the two validation sets.

Evaluating the performance of gene expression signature-based prognostic model within TNM stages, gender, and age groups

With high performance of the gene expression signature-based prognostic model in all LUAD patients from both training and validation sets, it was also necessary to investigate its performance in specific stage, gender, and age group of LUAD. As no samples were stratified into TNM stage IV group in the training set, we only focused on the performance of the model in another three stages (I, II, and III). For validation of the prognostic prediction value within TNM stages, gender, and age groups, Cox regression coefficients and dichotomization cut-off threshold generated from the training set were directly applied to the two validation sets. Similarly, significant overall survival difference was observed between high- and low-risk groups with each TNM stage, male/female, and old/young age groups in both of the validation datasets (Fig. 5, P < 0.05), except samples in male and old group of GSE37745, which may be resulted from its small sample size. These findings further validate the robustness of the gene expression-based signatures in predicting survival in lung adenocarcinoma.

Comparing signatures of 25 genes with known prognostic signatures in predicting LUAD prognosis

To demonstrate the robustness of the signatures of 25 genes in predicting LUAD prognosis, we built three more Cox models based on three signature gene sets found by previous studies (Der et al., 2014; Guo et al., 2006; Zhao, Li & Tian, 2018), which were selected from single dataset, and predicted the stratification of the two validation sets. We found that the three models showed worse ability in predicting the prognosis of patients in GSE37745 (Figs. 6B, 6D, and 6F), as compared with our signatures of 25 genes based Cox model (Fig. 4B), which may be caused by small sample size (n = 106). Although they had improved performance in TCGA-LUAD cohort (n = 502) (Figs. 6A, 6C and 6E), the significance levels of the three models were still worse than our model (Fig. 4A). These results suggested that the signatures of 25 genes were more robust than those selected by only one dataset.