A four-gene prognostic signature for predicting the overall survival of patients with lung adenocarcinoma

Data sources

We downloaded gene expression data (FPKM format), clinical information, and survival data, of patients with LUAD from The Cancer Genome Atlas (TCGA) from the University of California Santa Cruz Xenabrowser (UCSC Xena, https://xenabrowser.net/datapages/P) (Goldman et al., 2020) to act as a training cohort. We also acquired gene expression data and corresponding clinical information for the validation cohort under accession number GSE50081 from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). This validation dataset included 127 cases of lung adenocarcinoma.

Patients were included in our analyses if (1) clinical information (gender, age, T-stage, N-stage, M-stage), survival data (follow-up time and survival status), and gene expression information (mRNA expression levels) were complete and (2) the TNM stage was T1-4N0-2M0. We excluded cases if the histological type was not LUAD. After screening, a total of 317 patients were included in the training set and 127 patients were included in the validation set.

Prognostic gene signature screening

The set of training samples were randomly divided into three groups. The levels of gene expression in each group were then divided into high and low expression groups using the median value as the cutoff point. Next, we used a univariate Cox proportional hazard regression model to determine the association between gene expression and OS in each group; genes with a p-value < 0.05 (following the log rank test) were defined as being prognosis-related. We then identified the intersection between the three groups to identify candidate genes. Lasso regression and a multivariate Cox regression model were both used to carry out further screening; gene expression was used as continuous variable parameter. This additional analysis identified a prognostic gene signature consisting of four genes: CENPH, MYLIP, PITX3, and TRAF3IP3.

In the training set, the gene expression levels (z-score) of these four genes were used as covariates in a multivariate Cox regression model were applied. We then calculated risk scores based on gene expression levels and risk coefficients. Then, the entire cohort was divided into a high-risk group and a low-risk group using the median value as the threshold. Next, we compared survival data, clinical characteristics, and the immune microenvironment between the two different risk-groups and then validated the model in an external dataset. The immune microenvironment was estimated by TIMER2.0 (http://timer.comp-genomics.org/) (Li et al., 2020b), a public database, and six immune infiltration lymphocyte estimations were obtained.

Development and validation of the prognostic model

A multivariate Cox regression model was used to develop the prognostic mode. External validation was performed in an independent dataset. The degree of differentiation exhibited by the models was compared by comparing the area under the receiver operating characteristic curve (AUC-ROC). Calibration was evaluated using calibration curves.

Statistical analysis

All statistical analyses were performed with R software version 4.0.3. OS data were calculated using Kaplan–Meier curves, and the statistical difference between different groups was determined by the log-rank test. The influence of different parameters on OS was evaluated by univariate and multivariate Cox proportional hazard regression models. Hazard ratios (HRs) and the 95% confidence intervals (CIs) were generated using Cox proportional hazards models. In addition, receiver operating characteristic (ROC) analysis was carried out to compare the predictive accuracy of models. Pearson’s Chi-squared test was performed to compared the population demographics. A P-value < 0.05 was determined to be statistically significant.

Demographics of the study population

According to the inclusion and exclusion criteria, 317 cases were included in the TCGA cohort, while 127 cases were included in GSE50081. The median age of the two cohorts were 65-years and 70-years, respectively. The Stage of 5 patients in the TCGA cohort were not recorded; we defined these according to the 7th edition of the Union for International Cancer Control (UICC) TNM stage. The frequencies of T stage (p = 0.002), N stage (p < 0.001), and Stage (p < 0.001) in the two cohorts were different. Patient demographics are shown in Table 1.

Identification of a four-gene prognostic signature

The OS-related genes were identified from the training set using the Cox log-rank test and cut-off threshold of P < 0.05. The three randomized groups (A, B, C) contained 2061, 2742, and 699 OS-related genes, respectively. By identifying the intersection between the three gene sets, we identified 17 genes as candidate prognostic genes. Additional analysis, using Lasso and Cox stepwise regression, identified four prognostic genes: including two high risk genes (HR > 1; CENPH and PITX3) and two protective genes (HR < 1; MYLIP and TRAF3IP3). Figure 2 shows the gene screening process. Further details of the Cox regression model used to identify the four-genes are shown in Table 2. The Kaplan–Meier curves for the four genes in TCGA are shown in Fig. 3A.

Next, we validated the effect of the four prognostic genes on OS using the Kaplan–Meier Plotter (https://kmplot.com/analysis/) and obtained the same conclusion, as shown in Fig. 3B. Online analysis of the Gene Expression Profiling Interactive Analysis (GEPIA2, http://gepia2.cancer-pku.cn/) public database (Tang et al., 2019) showed that the expression of CENPH was up-regulated in tumors when compared to normal tissue and that expression levels increased with increasing tumor stage. In contrast, the other three genes were down-regulated but without statistical significance (—Log2FC— Cutoff = 1, p-value Cutoff = 0.01), as shown in Fig. 3C.

The risk score for the four-gene signature predicted the OS of patients with LUAD

Considering that the expression levels of the four identified genes exhibited a correlation with prognosis, we next explored the combined prognostic effect of these genes by using the ggrisk package in the R environment. Each patient was given an individual risk score by using a prognostic Cox regression model according to the expression level and its corresponding coefficients in the training set. Patients were then divided into a high risk group and a low risk group using the median risk score (0.02066664) as a cutoff. The same risk score method was performed for the validation set but using the median risk score (−0.00797257) as the cutoff point. Figure 4 shows the distribution of the risk scores, gene expression levels, and the survival status of patients in the training set. Two protective genes (MYLIP, TRAF3IP3) were highly expressed in the low-risk group, while risky genes (CENPH, PITX3) were highly expressed in the high-risk group.

Equations (1) and (2) show the formulae for calculating the risk score for the training and validation sets, respectively. These formulae used the same coefficients; these originated from the Cox model for gene risk as fitted to the training set (Table 2). (1)${\displaystyle \text{Risk score}\left(\text{Training set}\right)=TRAF3IP3\ast \left(-0.28593\right)+PITX3\ast \left(0.21193\right)+MYLIP\ast \left(-0.47968\right)+CENPH\ast \left(0.26692\right)}$ (2)${\displaystyle \text{Risk score}\left(\text{validation set}\right)=TRAF3IP3\ast \left(-0.28593\right)+PITX3\ast \left(0.21193\right)+MYLIP\ast \left(-0.47968\right)+CENPH\ast \left(0.26692\right)}$

Next, we performed univariate Cox regression analysis. Kaplan–Meier curves (Fig. 5) showed that patients in the high-risk group had a significantly shorter OS than those in the low-risk group (Training set: HR = 2.73, 95% CI [1.87–3.98], P < 0.001; Test set: HR = 2.72, 95% CI [1.52∼−4.86], P <  0.001). Subgroup analysis were further performed. In the training set (Fig. 6), the high-risk group had a shorter survival time than the low-risk group in almost every subgroup (Stage I, p = 0.003; Stage II, p = 0.007; Stage III, p = 0.047), although the p-value did not reach statistical significance in patients younger than 65 years (p = 0.058). Except for the Stage II subgroup (p = 0.075), all other subgroups showed significant survival differences in the test set (Fig. 7).

Risk-score was an independent prognostic factor for OS

Univariate Cox regression analysis for TCGA cohort suggested T stage (p < 0.001), N stage (p < 0.001), Stage (p < 0.001) and Risk group (p < 0.001) were significantly associated with OS (Table 3). The same results were obtained in the GES50081 cohort (Table 4). Next, we used the risk-group and other clinical factors (including age, gender, T-stage, N-stage and Stage) as covariates and created a multivariate Cox regression model for the training set and the validation set (Fig. 8). We found that only the risk-group was an independent prognostic factor in the training set, while risk-group and T3+4 were independent prognostic factors in the validation set. Stage and N stage in the validation set almost coincided exactly. Patients with Stage I were all N0, and 35 patients were Stage II; only 2 samples were N0. Consequently, multivariate analysis for Stage did not obtained result for Stage. After adjusting for other clinical factors, data showed that risk-group was an independent factor for the OS of patients with LUAD (training set: HR = 2.34, 95% CI [1.57–3.5], P  <  0.001; validation set: HR = 2.1, 95% CI [1.15–3.9], P = 0.017).

The relationship between risk-group and clinicopathological factors

Table 5 shows the proportions of patients in each risk-group by age, gender, T-stage, N-stage and Stage. It can be found that the high-risk proportion was significantly larger associated with a higher T-stage, N-stage, and Stage. This positive correlation also demonstrated that the risk score increased with tumor progression.

The immune microenvironment differed between risk groups

To determine if there were any differences with regards to tumor-infiltrating lymphocytes (TILs) relative abundance between high and low risk groups, we applied the online tool TIMER2.0 using gene expression data (Fig. 9). In the training set, the content of B-cells (p < 0.001), CD4+ T-cells (p < 0.001) was significantly lower in the high-risk group. In the validation set, the number of CD4+ T cells (p < 0.001) was significantly lower in the high-risk group.

We performed correlation analysis for the expression levels of the four genes, risk scores, and the relative abundance of TILs, in the two cohorts (Fig. 9). In the training set, the relative abundance of B-cells exhibited a significantly negative correlation with risk score (correlation coefficient = −0.38, p < 0.05). There was a negative correlation between CD4+T cells and risk score in both the training set and the validation set (correlation coefficient: −0.23, −0.49, p < 0.05). The expression levels of TAF3IP3 and all tumor-infiltrating lymphocytes were positively correlated in both datasets (correlation coefficient: 0.21–0.7, p < 0.05).

Development and validation of a composite prognostic model

Data derived from the training were used to develop the final prognostic model. As described previously, T-stage, N-stage, Stage and risk-group were shown to be prognostic factors via stepwise regression. We chose T-stage, N-stage and risk-group as parameters to establish the final model and referred to this final model as the composite prognostic model (named M5) which featured a combination of genetic and clinical factors. In addition, we performed model diagnostics for M5. Figure 10 shows that each covariate satisfies the proportional hazards risk hypothesis. None of the outliers affected the model estimates.

The performance of the composite model (AUC-ROC = 0.731) was significantly better than other parameters, including risk-group (M1, AUC-ROC = 0.643, p < 0.001), T-stage (M2, AUC-ROC = 0.642, p < 0.001), N-stage (M3, AUC-ROC = 0.648, p < 0.001), and Stage (M4, AUC-ROC = 0.646, p < 0.001). In the validation cohort, the AUC-ROC for M5 was 0.689, higher than other covariates; however, the difference between M2 and M5 was not significant (p = 0.11). The AUC values for M5 at 1, 3 and 5 years were 0.750, 0.737 and 0.719, in the training set, and 0.645, 0.766, and 0.725, in the validation set, respectively. Further details are given in Fig. 11.

Finally, we plotted a nomogram and calibration curves for the composite prognostic model (Fig. 12). The calibration curves showed a good match between the predicted probabilities and the actual probabilities (Fig. 13).