Lung cancer datasets were downloaded from the GEO database. A total of 739 patients were utilized in the present study after filtering out samples without clinical survival information. It included 293 patients from GSE30219 (19), 226 patients from GSE31210 (20), 168 patients from GSE37745 (21), and 52 patients from GSE19188 (22). We selected these datasets that included >50 patients with survival status information. We followed the strategy of using the largest dataset (GSE30219) as training set. Three independent datasets (GSE31210, GSE37745 and GSE19188) were included in the present study as testing sets.

All the microarray raw data (CEl files) of four lung cancer cohorts were processed using the robust multichip average (RMA) algorithm for background adjustment (23). GATExplorer was used to process microarrays on a local computer for gene expression of lncRNAs (24). lncRNA mapper was obtained from GATExplorer, which included the probes that do not map to any coding region but that were mapped to a database for non-coding RNAs of human and mouse (derived from RNAdb (25). The coding potential analysis of the lncRNAs was carried out by CNCI to classify protein-coding or non-coding transcripts (26). Each lncRNA included at least a minimum of three probes mapping in the corresponding lncRNA entity. We created a risk-score formula according to the expression of these eight lncRNAs for survival prediction. Patients having higher risk scores were expected to have poorer survival outcomes.

The association between the lncRNA gene expression and patient survival was assessed by univariable Cox proportional hazards regression analysis along with a permutation test using BRB-ArrayTools (Biometric Research Branch) package (27) in the training set. We identified expression of several lncRNAs that were strongly correlated with survival. Considering that a smaller number of genes in the model would make the model more practical, we performed the random survival forests variable hunting (RSFVH) algorithm (28). Using a smaller number of genes selected fitted in a multivariable Cox regression model; we constructed a formula to predict survival in the training set. Each patient was assigned a risk score that is a linear combination of the expression levels of the significant lncRNAs weighted by their respective Cox regression coefficients (29). According to this risk score, patients in the training set were divided into low-risk and high-risk groups using the median risk score as the cut-off. The Kaplan-Meier method was used to estimate survival time, and other three independent testing groups were performed for validation. Differences in survival times between the low-risk and high-risk groups in each set were compared by the two-sided log-rank test, respectively.

GSEA was performed by the JAVA program (http://www.broadinstitute.org/gsea) using MSigDB C2 CP: canonical pathway gene set collection (1,320 gene sets available). Gene sets with a false discovery rate (FDR) value <0.05 after performing 1,000 permutations were considered to be significantly enriched (30). Cytoscape (version 2.8.2) and the Enrichment Map software were used to visualize the GSEA results (31). Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the co-expressed protein-coding genes with prognostic lncRNAs were performed to predict the biological function of prognostic lncRNAs using the DAVID Bioinformatics Tool (version 6.7) (32). Enrichment analysis was carried out using the functional annotation clustering options, and was limited to KEGG pathways in the 'Biological Process' categories.

As summarized in the workflow (Fig. 1), all analyses were performed in the training set (GSE30219) and validated in the testing set (GSE31210, GSE37745 and GSE19188). The training set (n=293) was analyzed for the detection of prognostic lncRNA genes. By subjecting the lncRNA expression data derived from the training set to univariable Cox proportional hazards regression analysis using the BBRB-ArrayTools, we identified a set of 36 lncRNAs that were strongly correlated with patient overall survival (p<0.001 and FDR <0.001) from a total of 5,635 lncRNAs. Based on the random survival forests model (see Materials and methods), the eight lncRNAs (Table I) were selected as predictors (Fig. 2). In Table I a list of these eight lncRNAs is shown with their obtained coefficient and variable importance values. Based on these results, BC030759 was the most relevant with overall survival in the training set (HR=3.513); the positive coefficients of the lncRNAs (AK021595, BC030759, AK000053, BC020384 and AK022024) indicated that their higher levels of expression were associated with shorter survival, and the negative coefficients of the other lncRNAs (AK124307, CR615992 and AF085995) indicated that their higher levels of expression were associated with longer survival. All of the eight lncRNAs have been verified in the ncRNA Expression Database (www.nred.matticklab.com) and these eight transcripts were classified as ncRNAs in this website (33). As coding potential analysis is commonly used to classify whether a transcript is of coding potential or not (34), we used CNCI to test those eight transcripts (26). This tool also suggested that all the eight transcripts were non-coding transcripts with no coding potential.

To investigate whether the eight-lncRNA signature could provide an accurate prediction of survival in lung cancer patients, we created a risk-score formula according to the expression of these eight lncRNAs for survival prediction in the training set GSE30219 (n=293), as follows: risk score, 0.212*AK021595+0.416*BC030759 +0.322*AK000053−0.165*AK124307+0.301*BC020384+0.42 3*AK022024−0.084*CR615992−0.459*AF085995. Then, we calculated the eight-lncRNA signature risk score for each patient in the training set. Patients were divided into a low-risk or high-risk group using the median risk score as cut-off value. Patients in the high-risk group had a shorter survival time than patients in the low-risk group (Fig. 3A). The association of the eight-lncRNA risk score and survival was also significant when it was evaluated as a continuous variable in the univariable Cox regression model (Table II).

To confirm our findings, we calculated the risk score in the testing sets including GSE31210 (n=226), GSE37745 (n=168) and GSE19188 (n=52). Similar to the training set findings, patients in the high-risk group had a shorter survival time than patients in the low-risk group (Fig. 3B–D). Meanwhile, patient survival throughout the follow-up in the low-risk group was better when compared to survival in the high-risk group. In the univariable Cox regression model, the risk score was similar with the high-risk group which had a shorter overall survival. The patient survival status (Fig. 4A) and lncRNA values (Fig. 4B) were analyzed independently in the training set. Some of the clinical information (stage and subtype) was not available for a substantial proportion of cases, thus we performed multivariate Cox regression analysis concerning age and gender. The result showed that the eight-lncRNA expression signature was independent of age and gender. Eight-lncRNA risk score, age (available in GSE31210 and GSE37745) and gender (available in GSE31210, GSE37745 and GSE19188) were defined as covariates. These results showed that risk score was an independent predictor of lung cancer patient survival (Table II).

GSEA was carried out to identify the associated biological processes and signaling pathways (30). We compared the gene expression profile of lung cancer patients in the low-risk and high-risk groups classified by the eight-lncRNA gene signature in the training set (GSE 30219). The gene sets with significantly different expression (FDR <0.01; p<0.005) were picked up, which implied that the signature may be involved in the cell cycle and DNA replication-related pathways (Fig. 5A and B), and it was visualized as an interaction network with Cytoscape (Fig. 5C). These related pathways were reported to affect cancer cell proliferation (35–37).