A seven-long noncoding RNA signature predicts overall survival for patients with early stage non-small cell lung cancer

Identification of seven lncRNAs as prognostic biomarkers in discovery group

By subjecting the lncRNA profiling data of discovery group to univariable Cox proportional hazards regression analysis, we obtained 457 lncRNAs whose P-value were less than 0.05. Among them, seven lncRNAs were selected for further analysis, including APTR, DHRS4-AS1, ITGA9-AS1, LINC01137, LOC101927972, RPARP-AS1 and SH3BP5-AS1. Next, we fitted the 7 lncRNAs expression, the corresponding survival status and survival time into the multivariable Cox regression. Using the coefficients obtained from the multivariable Cox regression, a risk-score formula was constructed as following: risk score = - 0.06 x APTR - 0.20 x DHRS4-AS1 - 0.35 x ITGA9-AS1 - 0.23 x LINC01137 - 0.25 x LOC101927972 - 0.06 x RPARP-AS1 - 0.20 x SH3BP5-AS1. Notably, in the formula, every lncRNA expression value is weight by a negative coefficient (Table 1). As shown in Figure 1A and 1B, the heatmap revealed that the expression level of seven lncRNAs was decreased accompanying with the higher risk scores. Furthermore, we calculated the association between risk score and cancer-related death (Figure 1C). Our data showed that mortality rate in high risk group was significant higher than low risk group (Figure 1D), indicating the seven lncRNAs might play a protective role in NSCLC.

Figure 1. The seven-lncRNA signature-based risk score in prognosis of overall survival in the discovery group. (A) The seven lncRNA-based risk score distribution. (B) The heatmap of the seven lncRNA expression profiles. (C) The vital status of patients in high- and low-risk groups. (D) The mortality rate in low- and high-risk score groups.

The prognostic values of seven-lncRNA signature in discovery and validation groups

Considering the relationship between the seven-lncRNA signature and patients’ deaths, we next explored whether the seven-lncRNA signature could affect patients’ OS time. Our data in discovery group (Meta-GEO dataset, n=682) showed that the patients who had low risk-scores were supposed to have a longer OS time than higher risk group (Figure 2A). To validate above finding, we employed two validation groups, TCGA and GSE31210 dataset, containing 994 and 226 patients, respectively. As expected, patients in high risk group had a significant increased mortality risk than low risk group either in TCGA (HR=1.34, 95% CI= 1.10-1.63, P =0.004) or GSE31210 datasets (HR=3.06, 95% CI= 1.58-5.93, P =0.002) (Figure 2B and 2C). Moreover, prognostic meta-analysis based on above three groups (n=1902) confirmed that the seven-lncRNA signature was a risk factor for NSCLC patients (combined HR=1.76, 95% CI= 1.26-2.47, P=0.001) (Figure 2D).

Figure 2. The association between seven-lncRNA signature and overall survival in discovery and two validation groups. Kaplan-Meier survival curves were plotted to estimate the overall survival probabilities for the low- versus high-risk group in the discovery group (A), validation group-1 (B) and validation group-2 (C). (D) Meta-analysis was performed using the prognostic results of three groups.

The seven-lncRNA signature was associated with OS in stage I patients

To explore the impacts of clinical characteristics on the prognostic values of the seven-lncRNA signature, we performed a set of predefined stratified analysis. As shown in Table 2, we stratified the NSCLC patients by five clinical characteristics (including age, gender, smoking status, pathological subtypes and AJCC stage). According to the results of stratified analysis, the most significant result was found in stage I patients in both discovery and two validation groups (all P<0.001) (Table 2). To validate above findings, we further divided stage I patients into Ia and Ib according AJCC system, and calculated the prognostic values of the seven-lncRNA signature in each subgroups, respectively. Kaplan-Meier method was used to visualize the OS probabilities between high- and low-risk groups. As expected, the overall survival time of high risk group was all significantly shorter than low risk group in either stage Ia or Ib patients (Figure 3). The prognostic meta-analysis showed that the stage Ia patients with high risk scores had a two-fold increased mortality risk than those with low risk scores (combined HR=2.17, 95% CI= 1.15-4.09, P=0.017) (Figure 3D). Similarly, the association between the seven-lncRNA signature and OS also reached a statistical significance in stage Ib patients (combined HR=1.88, 95% CI= 1.30-2.72, P=0.001) (Figure 3H).

Figure 3. The association between seven-lncRNA signature and overall survival in patients with stage I. The Kaplan-Meier survival curves of discovery group (A), validation group-1 (B) and validation group-2 (C) were plotted, and meta-analysis (D) was conducted in patients with stage Ia. The similar results were obtained from the patients with stage Ib in discovery group (E), validation group-1 (F) and validation group-2 (G), and the prognostic meta-analysis (H) was also conducted.

The seven-lncRNA signature was associated with OS in patients carrying wild-type KRAS or EGFR genes

EGFR and KRAS were the most frequent mutant genes and associated with poor prognosis in NSCLC. Previous studies suggested that the mutation frequency of above two genes was significantly decreased in patients with early stage. Bearing this in mind, we also performed stratified analysis based on EGFR or KRAS mutation status. Our data showed that the higher risk score was associated with higher mortality risk in wild-type KRAS or EGFR subgroups (Figure 4). Above results were consistent in the two validation groups, suggesting the seven-lncRNA signature acted as a risk factor for patients carrying wild-type KRAS or EGFR genes.

Figure 4. Kaplan-Meier estimates of the overall survival of patients carrying wild-type KRAS or EGFR genes. (A) Kaplan-Meier survival curves were plotted to estimate the overall survival for patients carrying wild-type KRAS gene in two validation groups. (B) Kaplan-Meier survival curves were plotted to estimate the overall survival for patients carrying wild-type EGFR gene in two validation groups.

Cell cycle-related pathways were associated with the seven-lncRNA signature

Finally, to explore the potential mechanisms of the seven-lncRNA signature, we used the WGCNA method to cluster genes that highly correlated with the risk scores based on the profiling data of TCGA dataset [21]. We identified a total of 11 modules and found that pink module was the only positively correlated with the risk-score (P =2x10^-105) (Figure 5A and 5B). Pathway enrichment analysis was then performed on an online-based web tool “Metascape” (http://metascape.org/) using the genes in pink module. As shown in Figure 5C, genes were significantly enriched in cell cycle-related pathways, including “cell cycle”, “cell cycle phase transition” and “regulation of cell cycle process” pathways, suggesting the activation of cell cycle-related pathways might contribute to higher mortality risk in patients with high risk scores.

Figure 5. WGCNA predicted biological pathways associated with the seven-lncRNA signature. (A) The gene clusters obtained by WGCNA method. (B) The relationship between pink module membership and gene significance for risk score. (C) Significantly enriched pathways of the co-expressed genes in pink module.

Sample sources and study design

The raw data of gene expression and corresponding clinical information of NSCLC patients were downloaded from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) websites. Four GEO datasets (GSE30219, GSE31546, GSE37745 and GSE50081) were combined and used as the discovery group. After removal the five samples without survival data or enough clinical information, 682 patients from four GEO datasets were finally used for further analysis. Meanwhile, 994 patients from TCGA dataset were employed as validation group. Moreover, we also employed 226 patients from another independent GEO dataset (GSE31210) as the second validation group. All patients in GSE31210 belonged to the stage I and II of the American Joint Committee on Cancer (AJCC) system. The RNA expression profiling data of five GEO datasets were all performed on Affymetrix U133 Plus 2.0 microarray platform and the TCGA data was performed on Illumina sequencing platform. The clinical characteristics of discovery and validation groups were shown in Table 3.

Normalization and lncRNA annotation of GEO data and TCGA data

The raw data of five GEO datasets (GSE30219, GSE31210, GSE31546, GSE37745 and GSE50081) were downloaded as probe-level CEL files. Then, all raw data were quantile normalized using Robust Multi-array Average (RMA) method. With the lncRNA-specific probes presented on Affymetrix U133 Plus 2.0 downloaded from Affymetrix website (http://www.affymetrix.com) [15], we finally identified 2986 lncRNA transcripts with RefSeq transcript IDs. As for TCGA data, we downloaded the sequencing data with the type of fragments per kilobase of exon per million fragments (FPKM) from the TCGA website (https://tcga-data.nci.nih.gov/).

Construction of risk-score formula

After normalization and annotation of GEO raw data, we subjected 2986 lncRNAs to univariable Cox regression proportional hazards regression analysis to select lncRNAs which were associated with OS of NSCLC patients. Those lncRNAs with a statistical significance in univariable Cox regression were then selected into multivariable Cox regression to obtain the coefficients. By linearly combining the expression value of selected lncRNAs weighted by their coefficients, a risk-score formula was constructed as following [20]:

$\text{risk score }=\mathrm{ }\sum _{\mathrm{i}=1}^{\mathrm{N}}\left({Exp}_i\mathrm{*}{Coe}_i\right)$

In our formula, the $\mathrm{N}$ is the number of 7, the ${Exp}_i$ are the expression value of every 7 lncRNAs and the$\mathrm{ }{Coe}_i$ are their corresponding coefficients from the multivariable Cox regression analysis [27].

Survival analysis

According to above formula, the risk scores of NSCLC patients were calculated, by which patients were divided into high- and low-risk group with the cutoff of the median. The Kaplan-Meier method was used to assess the difference in survival time of high- and low- risk NSCLC patients. Difference of survival times between the two groups was considered as significant according to the P value of two-sided log-rank test less than 0.05.

Prognostic meta-analysis

To investigate the combined prognostic values of the seven-lncRNA signature of discovery and two validation groups, we performed a prognostic meta-analysis. All meta-analysis was performed on STATA software, version 12.0 (StataCorp LP, College Station, TX, United States). Pooled HR value was calculated using the random-effects model.

Weighted Correlation Network analysis (WGCNA)

To find the gene modules associated with our risk scores, we construct a co-expression network using the R package “WGCNA” according to our previous reports [28]. The soft thresholding power was selected to 9 to produce a weighted network. The enrolled genes were hierarchically clustered into 11 modules except the gray module.

Pathway enrichment analysis

The pink modules, the most significant modules being associated with risk score, were picked out to perform the pathway enrichment analysis. Pathway enrichment analysis of genes in pink module was performed on an online-based web tool “Metascape” (http://metascape.org/).