Identification of a multidimensional transcriptome prognostic signature for lung adenocarcinoma.

BACKGROUND: Lung adenocarcinoma (LUAD) is one of the leading contributors to cancer-related deaths worldwide. The objective of the current study is to identify a multidimensional transcriptome prognostic signature by combining protein-coding gene (PCG) with long non-coding RNA (lncRNA) for patients with LUAD.
METHODS: We obtained LUAD PCG and lncRNA expression profile data from three datasets in the Gene Expression Omnibus database and conducted survival analyzes for these individuals.
RESULTS: We established a predictive model comprising the three PCGs (NHLRC2, PLIN5, GNAI3), and one lncRNA (AC087521.1). This model segregated patients with LUAD into low- and high-risk groups based on significant differences in survival in the training dataset (GSE31210, n = 226, log-rank test P < .001). Risk stratification of the model was subsequently validated in other two test datasets (GSE37745, n = 106, log-rank test P < .001; GSE30219, n = 85, log-rank test P = .006). Time-dependent receiver operating characteristic (timeROC) curve analysis demonstrated that the model correlated strongly with disease progression and outperformed pathological stage in terms of prognostic ability. Cox proportional hazards regression analysis revealed that the signature could serve as an independent predictor of clinical outcomes in patients with LUAD.
CONCLUSIONS: We describe a novel multidimensional transcriptome signature that can predict survival probabilities in patients with LUAD.

area under the ROC curve

Lung adenocarcinoma

overall survival

receiver operating characteristic

INTRODUCTION

Lung adenocarcinoma (LUAD or LAC), is a major histological subtype of lung cancer1, 2 and one of the most common malignant tumors with high incidence and mortality. Lack of typical symptoms and signs in the early stages, patients with LUAD often progress to advanced stages at the time of diagnosis.3 As higher stage tumors with higher rates of recurrence, there is a significant proportion of patients with LUAD less than 5‐year survival.4, 5, 6 Therefore, besides histological classification, it is urgently need to develop novel molecular prognostic signature for predicting the risk of disease recurrence and identifying high‐risk subgroup of patients with LUAD who might benefit from adjuvant treatment.

With the development of high‐throughput technology, gene expression profiles have been broadly used to identify more novel biomarkers. Protein‐coding genes (PCGs) are the most common biomarkers and involved in the many key biological processes which can be powerful predictors of survival in patients in different cancers.7, 8, 9, 10 Recently, long non‐coding RNAs (lncRNAs) are transcripts >200 nucleotides with little coding capacity. Long non‐coding RNA (lncRNA) becomes new participant in tumorigenesis due to their various functions in a variety of cancer gene regulatory mechanisms, and has important clinical implications in terms of prognosis.7, 8, 9, 10, 11, 12, 13, 14 Recent studies have constructed many lncRNA signature15, 16, 17 to predict the prognosis of patients. For instance, a 3‐lncRNA signature can be a new biomarker for the esophageal squamous cell carcinoma prognosis,18 an immune‐related 6‐lncRNA signature could improve prognosis prediction of glioblastoma multiforme19 and a potential signature of eight long non‐coding RNAs could predict survival in patients with non‐small cell lung cancer.20 The advantage of combining PCGs with lncRNAs as prognostic markers is to show the disorder alteration of patients with cancer in greater detail from multiple dimensions.14, 21, 22, 23

Here, we analyzed PCG and lncRNA expression profiles of LUAD from Gene Expression Omnibus and developed a multidimensional transcriptome prognostic signature to predict LUAD survival.

MATERIALS AND METHODS

Expression data of LUAD patients

We acquired the expression data and associated clinical information of patients with LUAD from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Then, we performed a probe re‐annotation pipeline to get both PCG and lncRNA expression data. Specifically, we downloaded GPL570 probe sequences from Affymetrix (http://www.affymetrix.com) website and aligned these probe sequences to the human lncRNA and PCG transcript sequences from GENCODE (http://www.gencodegenes.org/), using BLASTn by the followed steps:(a) only retained the probes that matched to one PCG or lncRNA transcript. (b) Removed the probes matched to more than one transcript. (c) Each transcript should be perfectly matched to more than three probes.24

Construction of a prognostic signature in the training dataset

Survival‐related PCGs and lncRNAs in training dataset were screened out by cox proportional hazards regression analysis (P < .05). In an effort to make the dataset manageable, we used the random survival forests‐variable hunting (RSFVH) algorithm to filter genes until nine PCGs and lncRNAs.18 Subsequently, in order to further identify the prognostic genes, multivariable cox regression analysis was performed and a model to estimate prognosis risk was constructed as follows17: N is the number of prognostic genes, ExpVluei is the expression value of lncRNAs, and β i is the estimated regression coefficient of lncRNAs in the Cox proportional hazards regression analysis. Each patient was assigned 511 risk scores, since nine genes form 29−1 = 511 combinations. We chose prognostic signatures with AUC > 0.7 and log‐rank P < .05 from all 511 combinations, which were calculated by ROC and Kaplan‐Meier (KM) analysis.

Statistical analysis and bioinformatics prediction analysis of the prognostic genes function

Utilizing the ROC and the timeROC analysis, we compared the predictive efficacy of pathological stage with that of the PCG‐lncRNA signature. Cox proportional hazards regression analysis was performed to test whether the signature was an independent prognostic indicator, with significance defined as P < .05. All analyzes were performed using R program (www.r-project.org), including timeROC, survival, and randomforestSRC (downloaded from Bio‐conductor).

The co‐expressed relationships between PCGs and lncRNAs of the selected signature and all other protein‐coding genes were computed using Pearson's test; values with P < .05 and an absolute value of the Pearson coefficient > 0.3 were selected. We used the R package clusterProfiler to make Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment.25

RESULTS

Patient characteristics and expression profiles

Expression profiles of 417 samples, along with corresponding clinical data of patients diagnosed with LUAD, were downloaded from GSE31210, GSE37745, and GSE30219. The median age of the LUAD patients was 61 years (30‐83 years), and all patients were categorized as stage I, II, III, or IV (Table 1). Then, GSE31210 (n = 226) and GSE37745 (n = 106) were served as training sets while GSE30219 (n = 85) dataset was validation set.

Table 1

Summary of patient demographics and clinical characteristics

Characteristic	GSE31210	GSE37745	GSE30219
Age (y)
>61	122	47	46
≤61	104	59	39
Sex
female	121	60	19
Male	105	46	66
Vital status
Living	191	29	40
Dead	35	77	45
Pathological stage
Stage Ⅰ	168	70
Stage Ⅱ	58	19
Stage Ⅲ		13
Stage Ⅳ		4

Identification of prognostic genes from the training dataset

Through probe reannotating the Affymetrix Human Genome U133 Plus 2.0 Array, we obtained the lncRNA and PCG expression profiles of the 417 LUAD patients. Then, we selected 1897 PCGs and 529 lncRNAs associated with survival of patients with LUAD via cox proportional hazards regression analysis (P < .05). Seven PCGs and two lncRNAs (ANGPT4, MESD, ZMYM5, MANF, NHLRC2, PLIN5, GNAI3, AC006128.1, AC087521.1) with a strong correlation to patient survival were found according to the importance score calculated by random survival forests‐variable hunting (RSFVH) (Figure 1A,B).

Figure 1

Screening steps of the prognostic PCG‐lncRNA signature in the training dataset. Random survival forests‐variable hunting analysis revealed the lowest error rate for the data as a function of trees (A), and the associated scores were used to filter genes (B). The AUC of all 511 signatures was calculated and the first nine AUC are shown in the plot. D, ROC analysis of the selected prognostic PCG‐lncRNA signature(C)

Construction of the prognostic multi‐gene signature in the training dataset

The seven PCGs and two lncRNAs could generate 29−1 = 511 signatures, and each signature corresponded to a risk score (); detailed in Methods). ROC analyzes were performed on all 511 signatures and compared their AUC. The PCG‐lncRNA combination composed of three PCG (NHLRC2, PLIN5, GNAI3), and one lncRNA (AC087521.1) with the largest AUC (0.76) and minimum number of genes was selected (Figure 1C,D, Table 2). The risk score equation was calculated as: Risk score = (1.32 × expression value of AC087521.1) + (2.56 × expression value of GNAI3) + (−1.00 × expression value of PLIN5) +  (−1.76 × expression value of NHLRC2). The hazard ratio of the selected signature in the training group was 20.84 (P < .05), indicating that the PCG‐lncRNA is a risk factor of LUAD.

Table 2

Characteristics of PCGs and lncRNA in the signature

Database IDa	Gene symbol	Gene name	Coefficientb	P b	Expression with poor prognosis	Chromosome location
						( GRCh38/hg38 )
ENSG00000196865	NHLRC2	NHL repeat containing 2	−1.76	.00	low	10:113854661‐113917194:1
ENSG00000214456	PLIN5	Perilipin 5	−1.00	.01	low	19:4522531‐4535224:‐1
ENSG00000065135	GNAI3	G protein subunit alpha i3	2.56	.00	high	1:109548611‐109618321:1
ENSG00000244953	AC087521.1		1.32	.00	high	11:43943787‐43947206:−1

Ensembl database.

Derived from the univariable Cox regression analysis in the training set.

The selected signature for survival prediction in the training and test datasets

In the training dataset, each patient was assigned a risk score by the prognostic model based on the PCG‐lncRNA signature. As the median risk score as a cutoff point, patients from the training dataset were divided into a high‐risk group (n = 113) and a low‐risk group (n = 113).Then Kaplan‐Meier survival analyzes were performed and found patients from the high‐risk group had a significantly lower overall survival rate (OS) than those from the low‐risk group (log‐rank test P < .001; Figure 2A). When applied the median risk score to the GSE37745 and GSE30219 sets, patients from the two test sets were also divided into two groups, respectively, namely high‐risk groups(n = 53/42) and low‐risk groups (n = 53/43).Similarly, the survival of patients in the high‐risk groups was significantly shorter than those in the low‐risk groups (GSE37745 median 2.78, 95% CI: 1.46‐4.01 vs 5.94 years, 95% CI: 4.11‐7.22, log‐rank test P < .001, Figure 2B; GSE30219 median 4.58, 95% CI: 2.33‐12.5 vs 16.25 years, 95% CI: 8.58‐16.73, log‐rank test P = .0063, Figure 2C).

Figure 2

The PCG‐lncRNA signature predicted overall survival of patients with LUAD. Kaplan‐Meier survival curves classified patients into high‐ and low‐risk groups by the PCG‐lncRNA signature in the GSE31210 (A) and GSE37745, GSE30219 (B, C) datasets. P values were calculated by log‐rank test

According to the gene expression, risk score distribution and survival status of patients, Figure 3 illustrated the association of the gene expression with the survival. In the training dataset (Figure 3A), GSE37745 (Figure 3B), and GSE30219 (Figure 3C), patients with high expression of NHLRC2 and PLIN5 or low‐risk scores had a higher probability of survival, and patients with high‐risk scores or high‐expressed AC087521.1 and GNAI3 had shorter survival time.

Figure 3

Risk score distribution, survival status, and gene expression patterns of patients with LUAD in the GSE31210 (A) and GSE37745 (B), GSE30219 (C) dataset

The selected signature is an independent prognostic indicator

To better understand the clinical significance of the PCG‐lncRNA signature in patients with LUAD, we also examined the association of the PCG‐lncRNA signature with a series of clinical parameters in the dataset. There was no association between the PCG‐lncRNA signature and clinicopathological parameters in the training and test datasets, except pathological stage in the training set (Chi‐square test, P < .05, Table 3). Therefore, we performed a cox proportional hazards regression analysis to assess predictive independence of the PCG‐lncRNA signature. The P values of the prognostic signature in the cox proportional hazards regression analysis from the training datasets were <.05, which showed that the PCG‐lncRNA signature risk score was an independent prognostic indicator for patients with LUAD and was not affected by clinical features including sex, age, and pathologic stage (high‐risk group vs low‐risk group, HR = 15.79, 95% CI 3.70‐67.33, P < .001, n = 226, Table 4). The independence of the PCG‐lncRNA signature was validated in two test sets (high‐risk group vs low‐risk group, HR = 2.27, 95% CI 1.42‐3.63, P < .001, n = 106/HR = 2.39, 95% CI 1.28‐4.48, P = .01; Table 4).

Table 3

Association of the PCG‐lncRNA signature with clinicopathological characteristics in patients with LUAD

Variables	Train group		P	Test group 1		P	Test group 2		P
	Low riska	High riska		Low riska	High riska		Low riska	High riska
Age			.69			.43			1
≤61	63	59		26	21		23	23
>61	50	54		27	32		20	19
Sex			.02			1.00			.27
Female	70	51		30	30		7	12
Male	43	62		23	23		36	30
Pathological stage			.00			.49
Stage Ⅰ	100	68		37	33
Stage Ⅱ	13	45		8	11
Stage Ⅲ				5	8
Stage Ⅳ				3	1

Low risk ≤ median of risk score, high risk > median of risk score; The Chi‐squared test; P value < .05 was considered significant.

Table 4

Univariable and multivariable Cox regression analysis of the signature with LUAD survival

Variables		Univariable analysis				Multivariable analysis
		HR	95% CI of HR		P	HR	95% CI of HR		P
			lower	upper			lower	upper
GSE31210 dataset(n = 226)
Age	>61 vs ≤61	1.43	0.73	2.78	.29	1.32		2.10	.23
Sex	Male vs Female	1.52	0.78	2.96	.22	1.26	0.80	1.97	.32
Pathological stage	II vs I,	4.23	2.17	8.24	.00	1.32	1.03	1.68	.03
PCG‐signature	High risk vs low risk	20.84	5.00	86.93	.00	2.32	1.46	3.69	.00
GSE37745 set (n = 106)
Age	>61 vs ≤61	1.34	0.68	2.61	.40	1.14	0.71	1.84	.58
Sex	Male vs Female	1.11	0.57	2.19	.75	1.43	0.90	2.26	.13
Pathological stage	III, IV vs I, II	2.30	1.16	4.55	.02	1.27	0.99	1.64	.06
PCG‐signature	High risk vs low risk	15.79	3.70	67.33	.00	2.27	1.42	3.63	.00
GSE30219 set (n = 85)
Age	>61 vs ≤61	1.88	1.03	3.41	.04	1.85	1.01	3.37	.04
Sex	Male vs Female	1.02	0.49	2.13	.95	1.37	0.65	2.90	.41
PCG‐signature	High risk vs low risk	2.29	1.24	4.22	.01	2.39	1.28	4.48	.01

Comparison of the survival prediction efficiency of the PCG‐lncRNA signature with pathologic stage

Since GSE30219 without pathological stage information, we performed ROC analysis in two datasets (GSE31210/GSE37745, n = 226/106) to compare the survival prediction efficiency of pathological stage and the PCG‐lncRNA signature. The AUC of the PCG‐lncRNA signature was bigger than AUC of the pathological stage (Signature‐AUC = 0.76/0.68 vs Stage‐AUC = 0.65/0.62, Figure 4A,B). The high predictive efficacy demonstrated the PCG‐lncRNA signature has important clinical significance.

Figure 4

Comparison of the survival predictive power of the signature and that of pathological stage by ROC analysis in the GSE31210, and GSE37745 datasets (A, B). Survival predictive power of the signature (C, E) and pathological stage (D, F) at 3, 5, 8 years in the GSE31210 and GSE37745 dataset was analyzed by ROC analysis

TimeROC analysis was performed in the training dataset (n = 226), and we found that the AUC of the PCG‐lncRNA signature was greater than the AUC of the pathological stage (Signature‐AUC = 0.73/0.78/0.84 at 3/5/8 years vs Stage‐AUC = 0.75/0.64/0.73 at 3/5/8 years, Figure 4C,D). We also observed the same results in the GSE37745 dataset (Signature‐AUC = 0.64/0.63/0.62 at 3/5/8 years vs Stage‐AUC = 0.58/0.55/0.57 at 3/5/8 years, Figure 4E,F).

Gene oncology and KEGG pathway enrichment analysis

To characterize the molecular function of lncRNAs and PCGs in the PCG‐lncRNA signature, firstly, we screened out their co‐expressed protein‐coding genes from the GSE31210 and GSE37745 datasets and computed pearson correlation coefficients. Of these, 2654 protein‐coding genes were highly correlated with at least one of the selected genes, in the GSE31210 and GSE37745 datasets (Pearson correlation coefficient > 0.3/<−0.3, P < .05). Gene oncology and KEGG pathway enrichment analysis of the 2654 protein‐coding genes demonstrated that they were enriched in 38 gene oncology terms (GO terms) and KEGG pathways, including ncRNA transcription, response to insulin and snRNA transcription by RNA polymerase II checkpoint (P < .05, Figure 5).

Figure 5

Gene oncology and KEGG pathway enrichment analysis of the PCGs and lncRNAs in the signature

DISCUSSION

In recent years, the development of high‐tech sequencing makes novel PCGs or lncRNA signatures become a hot topic in cancer prognostic research. Although pathological stage is a commonly used prognostic method in clinical practice, its accuracy and effectiveness are insufficient for patients with LUAD.

In this study, we examined the clinical information and gene expression data of GSE31210 and identified a PCG‐lncRNA signature which could predict the survival of patients with LUAD. The PCG‐lncRNA signature was closely correlated with the overall survival rate of patients with LUAD in two test sets, indicating it could be a reliable indicator of survival. Cox proportional hazards regression analysis was performed to assess the independence of the selected PCG‐lncRNA signature in predicting the overall survival of patients with LUAD in the training and the test dataset. The PCG‐lncRNA signature maintained its correlation with the overall survival rate when coupled with age, gender, and pathological stage as covariables. This suggests that the predictive power of the PCG‐lncRNA signature is independent of these other clinical features. ROC analysis co‐founds that the prognostic ability of the signature is stronger than pathological stage, indicating that the signature could be an additional biomarker of the pathological stage.

We found high expression of GNAI3, AC087521.1 was associated with a short survival time (HR > 1, P < .05) and NHLRC2, PLIN5 was associated with a long survival time (HR < 1, P < .05). There was a study demonstrated that expression of GNAI3 shared a tight relationship with the prognosis of patients with hepatocellular carcinoma,26 but few research reported the function of AC087521.1, NHLRC2, and PLIN5 in cancer. While our study explored the function of these four prognostic genes by bioinformatic analysis, the biological role of them in LUAD tumorigenesis is still not clear and warrants further study. Additionally, experimental studies on these genes are needed to deepen understanding of the prognostic mechanisms behind the PCG and lncRNA, and enhance our understanding of their functional roles. In 417 LUAD samples, we confirmed that the signature is an effective marker for LUAD patients' prognosis, but this conclusion that selected PCG‐lncRNA signature may complement the pathological stage in a clinical setting would benefit from additional study.

In conclusion, using bioinformatics analysis, we identified a PCG‐lncRNA signature composed of AC087521.1, GNAI3, NHLRC2, and PLIN5 that accurately predicted the overall survival of patients with LUAD based on three LUAD independent datasets. However, additional large‐scale study is needed before the current results can be applied in clinical settings.

CONFLICT OF INTEREST

No potential conflicts of interest were disclosed.

AUTHORS' CONTRIBUTIONS

Jing Ye: data analysis, interpretation, and drafting, Hui Liu, Zhi‐Li Xu, Ling Zheng: data collection; Rong‐Yu Liu: study design, study supervision, and final approval of the manuscript. All authors read and approved the final manuscript.