Immune-related eight-lncRNA signature for improving prognosis prediction of lung adenocarcinoma.

BACKGROUND: Lung adenocarcinoma (LUAD) is the leading cause of cancer-related deaths worldwide. Therefore, the identification of a novel prediction signature for predicting the prognosis risk and survival outcomes is urgently demanded.
METHODS: We integrated a machine-learning frame by combing the Cox regression and Least Absolute Shrinkage and Selection Operator (LASSO) regression model to identify the LUAD-related long non-coding RNA (lncRNA) survival biomarkers. Subsequently, the Spearman correlation test was employed to interrogate the relationships between lncRNA signature and tumor immunity and constructed the competing endogenous RNA (ceRNA) network.
RESULTS: Herein, we identified an eight-lncRNA signature (PR-lncRNA signature, NPSR1-AS1, SATB2-AS1, LINC01090, FGF12-AS2, AC005256.1, MAFA-AS1, BFSP2-AS1, and CPC5-AS1), which contributes to predicting LUAD patient's prognosis risk and survival outcomes. The PR-lncRNA signature has also been confirmed as the robust signature in independent datasets. Further parsing of the LUAD tumor immune infiltration showed the PR-lncRNAs were closely associated with the abundance of multiple immune cells infiltration and the expression of MHC molecules. Furthermore, by constructing the PR-lncRNA-related ceRNA network, we interrogated more potential anti-cancer therapy targets.
CONCLUSION: lncRNAs, as emerging cancer biomarkers, play an important role in a variety of cancer processes. Identification of PR-lncRNA signatures allows us to better predict patient's survival outcomes and disease risk. Finally, the PR-lncRNA signatures could help us to develop novel LUAD anti-cancer therapeutic strategies.

INTRODUCTION

Lung cancer is the most common cause of cancer‐related death worldwide, with a 5‐year survival rate of about 16.6%, in which lung adenocarcinoma (LUAD) is one of the most common types of lung cancers. , For LUAD patients, early surgical resection is currently the standard treatment. After surgical intervention, additional chemotherapy will contribute to further improve patient's survival rate. However, nearly half of all LUAD patients experience a relapse and will die as a result of the disease recurrence. Therefore, identification of a novel prediction signature for predicting the prognosis and high‐risk group of LUAD patients is urgently demanded.

Currently, histopathology images serve as the golden standard for the diagnosis of lung cancer primarily evaluated by pathologists or doctors. Since, this process is time‐consuming and is influenced by the subjective judgment of the observer, identifying new biomarkers will facilitate the detection of lung cancer with improved patient's survival outcomes. Non‐coding RNAs (ncRNAs) are a class of RNA molecules without protein‐coding potential. Long non‐coding RNA (lncRNA) is defined as an ncRNA of at least 200 nucleotides (nt) in length that molecularly resemble mRNA. , Accumulating evidence indicates that lncRNAs exert important regulatory functions on gene expression at the transcriptional, post‐transcriptional, and epigenetic levels. , , Recently, lncRNAs as emerging cancer biomarkers have been used in the diagnosis and prediction of survival in multiple cancer types. , Fox example, lncRNA LNBC3 facilitates non‐small cell lung cancer progression by stabilizing BCL6; overexpressed lncRNA LINC00691 was a poor prognosis biomarkers in lung cancer. Hence, in this study, we developed a set of lncRNA signature to improve the patient's survival outcomes of LUAD.

Herein, we integrated the Cox regression and Least Absolute Shrinkage and Selection Operator (LASSO) regression model discovering an eight‐lncRNA signature (NPSR1‐AS1, SATB2‐AS1, LINC01090, FGF12‐AS2, AC005256.1, MAFA‐AS1, BFSP2‐AS1, and CPC5‐AS1) for improving prognosis prediction for LUAD patient. Next, the lncRNA signature has been confirmed as the robust prognosis‐related lncRNAs in independent datasets. Notably, we also found the lncRNA signature was associated with the abundance of multiple immune‐cell infiltration. Finally, we also constructed the lncRNA signature‐related competing endogenous RNA (ceRNA) network for identifying more potential therapy targets for LUAD.

MATERIALS AND METHODS

Data acquisition and pre‐processing

The Cancer Genome Atlas (TCGA) level 3 gene/lncRNA/miRNA expression data (n = 585, 526 cancer samples, 59 normal samples), and clinical data (527 patients) of LUAD were obtained from the Genomic Data Commons (GDC, available at https://www.cancer.gov/tcga). The GSE31210 independent LUAD datasets were obtained from publicly available Gene Expression Omnibus (GEO, available at https://www.ncbi.nlm.nih.gov/geo/)(15). For the gene/lncRNA/miRNA expression data, we removed the genes/lncRNAs that were not expressed over 70% of the samples.

Differential expression analysis of lncRNAs

We calculated the fold change (FC) of each lncRNA between cancer and normal samples, that is, , where and represent the mean expression of an lncRNA in cancer or normal samples, respectively. Then, we employed the Mann–Whitney U test to test whether the lncRNA expression was differentially expressed in cancer (H: there was no difference between cancer and normal samples). The lncRNAs with P‐value <0.05 and |FC| >1.5 were considered as the differentially expressed lncRNAs.

Identification of prognosis‐related lncRNA signature

We first randomly sub‐grouped the TCGA LUAD cancer samples into the training set and test set (70% and 30% of cancer samples, respectively). The univariate Cox regression that was operated in the R package Survival (v3.1–12) was employed on the training set to evaluate the correlation between the lncRNA expression levels and the patient's overall survival (OS). Then, we constructed the Least Absolute Shrinkage and Selection Operator (LASSO) regression model by R package glmnet (v4.0–2) to further screen the prognosis‐related lncRNA signature (PR‐lncRNA). , The features selected by the LASSO model were used to fit for a multivariate Cox risk regression model, where OS is the dependent variable and other clinical information is the covariate. Finally, the lncRNA with P‐value <0.01 as the PR‐lncRNA signature, which is significantly contributed to the LUAD patient's survival outcomes. The whole pipeline above was also performed in the test set and an independent dataset (GSE31210 ) to confirm the PR‐lncRNA signature was robust.

Functional enrichment analysis

Pearson correlation analysis was used to evaluate the co‐expression relationship between lncRNA signature and mRNAs. The significant co‐related mRNAs were used to annotate the biological functions of lncRNAs (p‐value <0.05). The functional enrichment analysis was performed by the R package clusterProfiler (3.16.1).

Survival analysis

Kaplan‐Meier survival curve that was operated in the R package survminer (v0.4.8) was used to prove the difference of overall survival (OS) between the high‐risk group and low‐risk group, and the bilateral logarithmic rank test was used to evaluate whether it was statistically significant. The risk score for each patient was calculated according to the linear combination of expression values weighted by the coefficient from the multivariate (univariate) Cox regression analysis: where n denotes the number of the PR‐lncRNA signature (n = 8), represents the coefficient of univariate Cox regression analysis when n = 1, was the coefficient of multivariate Cox regression analysis when n > 1, was the expression level of kth PR‐lncRNA expression of patient i. The median risk score was used as the cut‐off to divide patients into high‐ and low‐risk groups.

Tumor immunoactivity analysis

We first characterized the immune‐cell infiltration abundance using the TIMER. Next, we calculated Spearman's correlation between the immune‐cell abundance and the expression levels of PR‐lncRNA used R package stats (v4.0.2). The immune‐cell pathway enrichment was performed by the Immlnc (http://bio‐bigdata.hrbmu.edu.cn/ImmLnc/).

Construction of competing endogenous RNA regulatory network

We first identified the differentially expressed miRNAs and mRNAs in LUAD using the Mann–Whitney U test (p < 0.05). Next, we calculated the coupled Spearman's correlation coefficient (rho) among the expression levels of miRNAs, PR‐lncRNAs, and mRNAs. The raw P‐values () were adjusted by multiple hypotheses using a permutation method. For each lncRNA, the expression value was held consistent, and 10,000 random miRNAs and mRNAs were used to perform the same Spearman's correlation test, generating sets of 10,000 permutation P‐values (). Then, the empirical P‐value () was corrected, that is, which introduces a pseudo‐count of 1. The significantly positively correlated lncRNA‐mRNA pairs ( and ), and negatively correlated miRNA‐lncRNA and miRNA‐mRNA pairs ( and ) were integrated as the miRNA‐lncRNA‐mRNA regulatory axes. The ceRNA network was visualized using the Cytoscape.

Statistical analysis

All analyses were conducted using R (v3.6.3) software.

RESULTS

An eight‐lncRNA signature shows association with the prognosis of LUAD

Long non‐coding RNA (lncRNA) is an emerging biomarker, which played an important role in cancer development. To explore the role of lncRNA in LUAD patient survival, we first identified the differentially expressed lncRNAs between LUAD cancer and normal samples (See the Methods section). A total of 309 up‐regulated and 119 down‐regulated lncRNAs were recognized (Figure 1A). Next, we developed a computational model by combing the LASSO regression and Cox regression models to identify the prognosis‐related lncRNA (PR‐lncRNA) signature (Figure 1B‐C, see the ‘Methods’ section). The PR‐lncRNA signature consisting of eight lncRNAs (NPSR1‐AS1, SATB2‐AS1, LINC01090, FGF12‐AS2, AC005256.1, MAFA‐AS1, BFSP2‐AS1, and CPC5‐AS1), contributed to the LUAD patient survival outcomes. Of them, NPSR1‐AS1, SATB2‐AS1, FGF12‐AS2, and AC005256.1 could be acted as the poor prognosis markers, while LINC01090, MAFA‐AS1, BFSP2‐AS1, and CPC5‐AS1 were treated as the protective factors for survival in LUAD patients. In the previous studies, NPSR1‐AS1 had been reported as the poor prognosis markers for LUAD ; SATB2‐AS1 could be broadly involved in multiple cancer development , , ; FGF12‐AS2 may be served as a potential target for the treatment of non‐small‐cell lung carcinoma. The evidence further highlighted the PR‐lncRNA signature might exert an essential contribution for LUAD processes and have the potential clinical utility.

FIGURE 1

The PR‐lncRNAs contributes to patient survival outcomes in LUAD. (A) Volcano plot showing the differentially expressed lncRNAs in LUAD. x‐axis represents the log2(Fold Change). Y‐axis shows the log10(p‐value). (B‐C) Lasso regression model screened characteristic lncRNA. Variation curve of regression coefficient and β value. (D) Forest plot of multivariate Cox regression results, which include p‐value and confidence interval of hazard ratios

The PR‐lncRNA signature contributes to the patient survival outcome prediction

We next used the nomogram method to build a more intuitive prediction model of one‐ and three‐year survival time probability based on the PR‐lncRNA expression. Under the same gene expression fluctuations, the scores of the PR‐lncRNAs that are more important to the patient's survival will rise (Figure 2A). The calibration curve also showed that there was a good prediction capability for the one‐ and three‐year survival time probability (Figure 2B). Besides, the Receiver Operating Characteristic (ROC) curve reflected that the sensitivity and specificity of the PR‐lncRNA predictive model. We observed the promising potential for predicting patient's survival outcomes (Figure 2C). These results suggested that the PR‐lncRNA signature was a reliable biomarker contributing to the LUAD patient's survival outcomes.

FIGURE 2

PR‐lncRNA signature could act as the promising survival biomarker of LUAD patient. (A) Nomogram showing the prediction model of survival probability. (B) Calibration curve of nomogram. (C) ROC curve reflecting the performance of the Cox regression model in predicting the probability of survival of patients at different time nodes

PR‐lncRNA signature prompts LUAD patients' risk score

We next constructed the risk scoring model based on the multivariate Cox regression analysis and patient's overall survival (OS) time (see the ‘Methods’ section). , The risk scoring model could be expressed as Risk score = (0.070 * NPSR1‐AS1 + 0.067 * SATB2‐AS1‐0.057 * LINC01090 + 0.076 * FGF12‐AS2 + 0.498 * AC005256.1 – 0.040 * MAFA‐AS1 ‐ 0.044 * BFSP2‐AS1 ‐ 0.040 * GPC5‐AS1). Based on the model above, we calculated the risk score for each patient in the TCGA training set. The samples were subsequently divided into the high‐risk group (n = 179) and low‐risk group (n = 180) by considering the median values of all patients' risk scores as the cut‐off point. We found there was an obvious survival difference between the high‐ and low‐risk group (Figure 3A) and the risk score was also a poor prognosis marker in LUAD (p < 0.0001) (Figure 3B). The same results were also repeated in the TCGA test set, implying the PR‐lncRNA signature was a robust survival‐related biomarker (Figure 3CandD).

FIGURE 3

PR‐lncRNA signature is associated with the patient risk score in LUAD. A‐B. Risk scores, survival status, the expression of PR‐lncRNA signature (A), and Kaplan‐Meier curves (B) for two groups in the training set. C‐D. Risk scores, survival status, the expression of PR‐lncRNA signature (C), and Kaplan‐Meier curves (D) for two groups in the test set

Evaluation of the PR‐lncRNA signature robustness for LUAD patient survival risk

To further confirm the PR‐lncRNA signature is a robust biomarker in LUAD, we employed the risk scoring model in an independent LUAD dataset (GSE31210, n = 246). Consistent with the results in the TCGA cohort, the risk score was also identified as a poor prognosis marker (p = 0.027) (Figure 4A‐B). Moreover, the overexpressed PR‐lncRNAs NPSR1‐AS1 (p = 0.0084) and STB2‐AS1 (p = 0.0015) were also closely associated with the shorter survival times (Figure 4C‐D). Furthermore, by exploring the relationships between the disease stages and patient's risk score, we found the more advanced disease stages usually presented the higher the PR‐lncRNA risk score (Figure 4E). These results further supported that the PR‐lncRNA signature was a robust prognosis indicator in LUAD.

FIGURE 4

PR‐lncRNA signature is a robust prognosis indicator in LUAD. A‐B. Risk scores, survival status, the expression of PR‐lncRNA signature (A), and Kaplan‐Meier curves (B) for two groups in the GSE31210. C‐D. The Kaplan‐Meier curves of the high‐/low‐ expression group of NPSR1‐AS1 (C) and STB2‐AS1 (D) in the GSE31210 cohort. Patients were divided by the median expression value of lncRNA. E. Boxplot showing the relationship between LUAD disease stages and patients’ PR‐lncRNA risk score. The p‐values were calculated by the Mann–Whitney U test

PR‐lncRNA signature shows crosstalk with the tumor immunoactivity

To investigate the role of PR‐lncRNA signature in the tumor immune microenvironment, we systemically analyze the correlation between PR‐lncRNA expression and immune‐cell infiltration abundance. Notably, the PR‐lncRNAs generally showed significant correlations with multiple immune‐cell infiltration abundance. Of them, PR‐lncRNA BFSP2‐AS1 exhibited significant associations with multiple immune cells, that is, B cell, CD4+ T cell, CD8+ T cell, Dendritic, Macrophage, and Neutrophil, implying BFSP2‐AS1 may be an immune‐related lncRNA in LUAD (Figure 4A). Besides, we also performed the immune pathway enrichment analysis using the Immlnc and constructed the PR‐lncRNAs and immune pathway network (Figure 5B). In the network, we observed the PR‐lncRNAs AC005256.1 and SATB2‐AS1 were engaged in the greatest number of immune pathways. Both of them could be involved in the Antimicrobials, Cytokines, and Antigen Processing and Presentation pathways. Additionally, AC005256.1 was involved in the Cytokine Receptors pathway, SATB2‐AS1 was specifically enriched in the Chemokines pathways (Figure 5B). The previous study has been demonstrated that SATB2‐AS1 could inhibit tumor metastasis and affects the tumor immune‐cell microenvironment in colorectal cancer by regulating SATB2 expression. Furthermore, the PR‐lncRNAs showed correlations with the expression levels of Major Histocompatibility Complex (MHC). MHC class Ⅰ and MHC class Ⅱ are mainly responsible for the antigen presentation required for tumor cell recognition and subsequent T cell‐mediated killing. We found the PR‐lncRNA signature was also closely associated with the expression levels of MHC class Ⅰ (Figure 5C) and MHC class Ⅱ (Figure 5D), implying the PR‐lncRNAs may be involved in the tumor immunoactivity alterations. In summary, these results indicated that PR‐lncRNAs were immune‐related lncRNAs and showed crosstalk with the LUAD tumor immunoactivity.

FIGURE 5

PR‐lncRNAs are LUAD tumor immunoactivity‐related lncRNAs. A. Bubble plot showing the correlation between tumor immune‐cell abundance and the expression of PR‐lncRNA. B. Network diagram showing the immune‐related pathway and PR‐lncRNA associations. C‐D. Heatmaps showing the correlation between the expression levels of MHC and PR‐lncRNAs. (C) MHC‐I, (D) MHC‐II. The p‐values were calculated by Spearman's correlation test

PR‐lncRNA related competing endogenous RNA identifies potential therapy targets

To better understand the biological functions of PR‐lncRNAs, we explored the PR‐lncRNA‐related competing endogenous RNA (ceRNA) regulatory pattern. We identified the miRNA‐lncRNA‐mRNA regulatory axis to constructed the PR‐lncRNA‐related ceRNA network (Figure 6A, see the ‘Methods’ section). The network contented two miRNAs, one lncRNA, and eight mRNAs comprising 14 ceRNA axes. Notably, only PR‐lncRNA NPSR1‐AS1 was recognized as the partner of the ceRNA axis, implying its importance for understanding the LUAD development at molecular levels. Subsequently, we systemically annotated the biological pathways involved in the mRNAs in the ceRNA network (Figure 6B). Remarkably, mRNA KRAS was engaged in multiple KEGG pathways, such as the mTOR signaling pathway, PI3K‐Akt pathway, VEGF signaling pathway, Longevity regulating pathway, etc. Previous studies have also demonstrated KRAS was an oncogene in various cancer types. , Therefore, we next focused on the KRAS‐related ceRNA axis (has‐miR‐365a‐3p ‐ NPSR1‐AS1 ‐ KRAS). The KRAS expression was significantly up‐regulated in LUAD (Figure 6C) and positively correlated with the NPSR1‐AS1 expression that was also significantly up‐regulated in LUAD (Figure 6D and E), indicating they were co‐expressed in LUAD. Besides, the has‐miR‐365a‐3p expression was negatively associated with the expression levels of NPSR1‐AS1 (Figure 6F) and KRAS (Figure 6G), which also significantly elevated (Figure 6H). All these implying that there was a potential competing relationship between NPSR1‐AS1 and KRAS for has‐miR‐365a‐3p. Finally, overexpressed KRAS was also identified as a poor prognosis biomarker in the TCGA LUAD cohort and GSE50081 cohort (Figure 6I and J). These results highlighted PR‐lncRNA related ceRNA axis was a promising potential therapy target for LUAD treatment.

FIGURE 6

Integrative PR‐lncRNA related ceRNA analysis. A. ceRNA network showing the competing endogenous RNA relationships among PR‐lncRNAs, microRNAs, and mRNAs. B. Circos plot showing the key mRNA‐associated KEGG pathways. C‐H. The hsa‐miR‐365a‐3p‐NPSR1‐AS1‐KRAS axis in LUAD. The p‐values were calculated by the Mann–Whitney U test (C), (E), and (H) and Spearman's correlation test (D), (F), and (G). I‐J. The KM curve showing the patient survival outcome of high‐/low‐ expression of KRAS in TCGA LUAD (I) and GSE50081 (J)

DISCUSSION

Accumulating evidence suggested that the lncRNAs were promising biomarkers for predicting patients’ survival outcomes in multiple cancers. , Identification of novel lncRNA signatures would lead to facilitate improving the tumor detection and prognosis treatment. LUAD, a common poor prognosis type of lung cancer, has a patient 5‐year survival rate of only 16.6%. Herein, we integrated Cox regression and LASSO regression models discovering an eight‐lncRNA signature (PR‐lncRNA signature) that was closely related to patient's prognosis outcomes. The PR‐lncRNA signature as a robust prognostic marker for LUAD demonstrated a favorable ability to predict patient's survival outcome and prognostic risk.

Previous studies have demonstrated tumor immune‐cell infiltration plays an important role in the development and progression of cancer. , The survival outcome of LUAD patients is closely related to the level of tumor immune‐cell infiltration. Notably, we also found the PR‐lncRNAs showed crosstalk with multiple immune‐cell infiltration abundances, such as B cell, CD4+ T cell, CD8+ T cell, Dendritic, Macrophage, and Neutrophil. Further parsing of immune pathway enrichment analysis suggested that PR‐lncRNAs could be involved in various immune‐related pathways. Finally, we also observed the expression levels of PR‐lncRNAs were associated with the Major Histocompatibility Complex (MHC), implying they might participate in antigen presentation required for tumor cell recognition and subsequent T cell‐mediated killing.

LncRNAs have sophisticated functions through various pathways, and the ceRNA hypothesis complicates the relationship between miRNAs, lncRNAs, and mRNAs. In this study, we also constructed a ceRNA network for the PR‐lncRNA signature. The analysis highlighted the PR‐lncRNA NPSR1‐AS1 was an essential ceRNA biomarker for LUAD survival outcomes. We found NPSR1‐AS1‐related ceRNA regulatory axis (has‐miR‐365a‐3p ‐ NPSR1‐AS1 ‐ KRAS) could be involved in multiple cancer‐related pathways. And the oncogene KRAS was also identified in various cancer types. , These results indicated that the PR‐lncRNA related ceRNA network will be contributed to finding more anticancer therapy targets.

In summary, our study identified a set of eight‐lncRNA signature, which contributes to predicting the LUAD patients' prognosis risk and survival outcomes. The lncRNA signature has been characterized to be closely associated with the abundance of multiple immune cells infiltration and the expression of MHC molecules. Finally, by constructing a ceRNA network, we found more PR‐lncRNA‐related anti‐cancer therapy targets.

CONFLICTS OF INTEREST

The authors declared no potential conflicts of interest in terms of the research, authorship, and/or publication of this article.

AUTHOR CONTRIBUTIONS

Y.C. and X.Z. performed analysis. J.L. collected the dataset. M.Z. designed the study. All authors read and approved the final manuscript.

CONSENT FOR PUBLICATION

Not applicable.