Identification of immune-associated lncRNAs as a prognostic marker for lung adenocarcinoma.

BACKGROUND: Lung adenocarcinoma (LUAD) accounts for the largest proportion of lung cancer patients and has the highest morbidity and mortality worldwide. Accumulating evidence shows that immune-associated long non-coding RNAs (lncRNAs) play a role in LUAD, although their predictive value for immunotherapy treatment and cancer-related death remains poorly investigated.
METHODS: Gene expression profiles and clinical data were obtained from The Cancer Genome Atlas. We constructed a risk model by univariate and multivariate Cox regression and least absolute shrinkage and selection operator regression analysis and subsequently divided each sample into low- or high-risk category. Survival and receiver operating characteristic (ROC) analyses were applied to assess the prognostic value of the model. Additionally, immune and somatic mutation status were analysed between the two risk groups. Finally, the model was applied to pancreatic ductal adenocarcinoma (PDAC) samples to explore the applicability of the model in other cancers.
RESULTS: We obtained data from 499 LUAD patients and randomised the samples into a training set (N=351) and validation set (N=148) at a 7:3 ratio. We detected 7 immune-associated lncRNAs (AP000695.2, AC026355.2, LINC01843, ITGB1-DT, LINC01150, AL590226.1 and AC091185.1) that were applicable for establishing a risk signature. Survival analysis revealed that patients categorised in the high-risk group had shorter overall survival (OS) than those in the low-risk group. ROC analyses showed excellent AUC values in all data sets (>0.65 at 1, 3, and 5 years). Notably, ESTIMATE algorithm and analysis of PCA, (ss)GSEA, and somatic mutations revealed that the high-risk group had a stronger immunosuppressive status and a higher tumour mutation burden (TMB). Moreover, patients in the low-risk group responded better to immunotherapy due to higher levels of immune-checkpoint receptor genes and TLS-related genes. Our model using the 7 immune-associated lncRNAs showed similar applicability for PDAC patients.
CONCLUSIONS: We constructed a model for risk signatures based on 7 immune-associated lncRNAs and showed its prognostic value for identifying immune and somatic mutation characteristics in LUAD patients, which may assist clinical treatment plans and elucidate molecular mechanisms of LUAD immunity. 2021 Translational Cancer Research. All rights reserved.

Introduction

Lung cancer is one of the most prevalent malignancy with approximately 2.1 million new cases and 1.59 million deaths worldwide in 2018 (1). Lung adenocarcinoma (LUAD) is the most common pathological subtype of non-small cell lung cancer (NSCLC) and amounts to 40% of all lung cancers across the world (1,2). Previous studies have reported that despite the development of targeted molecular therapies and treatment strategies, the 5-year survival rates in LUAD patients remain low, especially in advanced patients (3-5). Cancer immunotherapy shows great promise and immune-checkpoint inhibitors for PD-L1, CTLA-4, and TIM-3 have demonstrated effectiveness in treating advanced lung cancer (6-8). However, the selection criteria for sensitive populations suitable for immunotherapy remains undefined (7,9). A growing number of studies have recognized that the tumour microenvironment (TME) plays a critical role in LUAD and the fraction of stromal cells in the TME may influence the patients’ prognosis (10). Recent studies have indicated that a high immune score is associated with improved survival in LUAD, which suggests that the immune microenvironment may be a prognostic factor (11).

LncRNAs have a length of over 200 nucleotides without coding functions and studies have demonstrated an array of biological functions (12,13). An increasing number of abnormal biological behaviours in tumours have been linked to lncRNAs including the unusual activity of the genetic material and changes in the TME, which are involved in processes promoting tumour progression such as immune escape, DNA damage, metabolic disorders, epithelial-mesenchymal transition, cell stemness, and chemical resistance (14). Previous studies have also investigated the regulation of tumour immunity by lncRNAs (13,15,16). For example, lncRNA LNMAT1 affects lymphatic metastasis of bladder cancer by inducing CCL2 in the TME via the recruitment of tumour-associated macrophages (17). In addition, Lnc-SNHG1 promotes immune escape of breast cancer by increasing the differentiation of regulatory T cells (18). The prognostic value of lncRNAs in NSCLC have also received an increasing amount of interest. For instance, a risk signature composed of 7 lncRNAs has been proposed by Lin et al. to predict overall survival (OS) for early-stage NSCLC patients (19) and Miao et al. showed similar findings for 8 lncRNAs in elderly NSCLC patients (20). However, the prognostic value of immune-associated lncRNAs has not been widely studied.

In this study, we obtained 12 immune-associated and prognostic lncRNAs according to the co-expression network and univariate Cox proportional hazards regression analysis. Then we constructed a risk signature of 7 immune-associated lncRNAs that clustered LUAD samples into low- or high-risk based on this signature and we validated the signature’s prognostic and clinical value. Unlike previous studies, we analysed the immune status using 29 immune-associated gene sets, which represented diverse immune cell types, functions, and pathways based on previous literature (21-23), and employed the ESTIMATE algorithm to verify consistency of the results. Furthermore, we analysed the characteristics of genomic alterations including somatic mutations and copy number variations in multiple risk groups. Finally, we applied the risk model to pancreatic ductal adenocarcinoma (PDAC) samples to explore whether the model is valid for other tumours.

We present the following article in accordance with the MDAR checklist (available at http://dx.doi.org/10.21037/tcr-20-2827).

Methods

Data acquisition and processing

Gene expression profiles and clinicopathological features of LUAD and PDAC patients were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) with standardised gene expression in fragments per kilobase million. The following exclusion criteria were applied to ensure data integrity: (I) patients with survival time less than 30 days and (II) missing information on survival status, TNM, AJCC stage, age, and gender. We included a total of 499 LUAD patients and 170 PDAC patients in this study. The 499 LUAD patients were randomly assigned to a training set (N=351) or a validation set (N=148) in a 7:3 ratio (Table S1). The gene transfer format was downloaded from GENCODE (https://www.gencodegenes.org/) to convert the Entrez IDs to gene IDs (24,25) and the lncRNA profiles were extracted from mRNA expression data. We obtained 332 immune‐related genes categorised as immune system process (M13664) or immune response (M19817) through the molecular signatures database (http://www.broadinstitute.org/gsea/msigdb/index.jsp) (21,26). The expression of immune-related genes was further correlated with lncRNA expressions based on Pearson’s correlation analysis.

Risk score construction and survival analysis

Univariate and multivariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression analysis (27) were performed on the immune-associated lncRNAs based on lncRNAs expressions and survival data from patients in the training group. The risk scores were obtained according to the following algorithm as described previously (20,28).

Where Expi is the normalised expression level of each of the immune-associated lncRNAs and Coei is the weighted regression coefficient of each item (20). LUAD patients were subsequently divided into high- and low-risk group based on the cut-off value of the median risk score of the training group. Kaplan-Meier survival analysis and receiver operating characteristic (ROC) curve was used to differentiate the survival time between low- and high-risk LUAD patients.

Evaluation of immune status and somatic mutations status

We used 29 immune signatures in the single-sample gene-set enrichment analysis (ssGSEA) to analyse the immune status of the high- and low-risk samples as described previously (23). ESTIMATE algorithm was performed to evaluate the immune cell infiltration level (immune score), stromal content (stromal score), and tumour purity (29). Expression patterns in the different groups were investigated by principal component analysis (PCA) and gene set enrichment analysis (GSEA) was performed to study the different biological processes and pathways enriched in the high- or low-risk groups (http://www.broadinstitute.org/gsea/index.jsp). All gene somatic mutations of LUAD and PDAC samples were also downloaded from TCGA. We then analysed the variant classification and type, classification summary, single nucleotide variation (SNV) class, variants per sample, and top ten mutated genes.

Statistical analysis

R version 4.0.1 and SPSS Statistics version 23.0 were used for all statistical analyses. To ensure equal clinical baseline of the training and validation set, we applied the Chi-squared test, Fisher’s exact test, or Student’s t-test. Wilcoxon-Mann-Whitney or Kruskal-Wallis tests were used to compare different groups depending on the number of comparisons. All reported P values were two-tailed and P<0.05 was considered statistically significant.

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Model construction of 7 immune-associated lncRNAs in the training group

The expression of 332 immune-associated genes in LUAD samples are shown in https://cdn.amegroups.cn/static/public/TCR-20-2827-1.xlsx. A total of 788 immune-associated lncRNAs were identified by the Pearson’s correlation analysis with |R| >0.5 and P<0.001 (https://cdn.amegroups.cn/static/public/TCR-20-2827-2.xlsx). Univariate Cox proportional hazards regression analysis was then performed on these 788 lncRNAs in the training set, resulting in 12 lncRNAs with a P value of <0.01 (https://cdn.amegroups.cn/static/public/TCR-20-2827-3.xlsx). We further applied LASSO for the 12 lncRNAs, of which the following 9 lncRNAs significantly correlated with the survival of LUAD patients: AP000695.2, AC026355.2, LINC01843, ITGB1-DT, DPYD-AS1, LINC01150, AL590226.1, AC091185.1, and AC005034.3 (). Multivariate Cox regression was performed and the coefficients were obtained to construct a risk score composed of 7 immune-associated lncRNAs with the following formula.

Figure 1

Construction of immune-associated lncRNA risk signature and prognostic analysis of the training set in LUAD samples. (A) LASSO coefficient profiles of lncRNAs from univariate Cox proportional hazards regression and ten-fold cross-validation of the 9 prognostic genes. (B) The risk score and survival status of patients in the training set. (C) The expression level of the 7 prognostic lncRNAs in low- and high-risk patients.

Risk score =0.251 × AP000695.2 - 0.334 × AC026355.2 + 0.085 × LINC01843 + 0.123 × ITGB1-DT - 0.476 × LINC01150 - 0.583 × AL590226.1 + 0.524 × AC091185.1 (). According to the median risk score, the training set was divided into the high- or low-risk group. The relationship between risk score and cancer-related mortality showed that the mortality rate in the low-risk group was significantly lower than in the high-risk group (). Moreover, a heatmap of the expression level of the 7 immune-associated lncRNAs revealed that AP000695.2, LINC01843, ITGB1-DT, and AC091185.1 may play a harmful role in LUAD whereas AC026355.2, LINC01150, and AL590226.1 may exert protective effects ().

Table 1

Seven immune-related lncRNAs significantly associated with the OS of LUAD patients in the training group

lncRNA_symbol	Ensemble ID	Coefficient	Univariate Cox regression analysis
			HR	95% CI lower	95% CI higher	P value
AP000695.2	ENSG00000233818	0.251	1.411	1.153	1.726	<0.01
AC026355.2	ENSG00000236385	−0.334	0.771	0.647	0.918	<0.01
LINC01843	ENSG00000251169	0.085	1.111	1.042	1.184	<0.01
ITGB1-DT	ENSG00000229656	0.123	1.194	1.075	1.326	<0.01
LINC01150	ENSG00000229671	−0.476	0.527	0.328	0.848	<0.01
AL590226.1	ENSG00000228401	−0.583	0.397	0.217	0.727	<0.01
AC091185.1	ENSG00000253476	0.524	1.343	1.09	1.655	<0.01

Validation of the prognostic value of immune-associated lncRNA signature

Further examination of the validation and combination set performed with the same algorithm, risk-score formula, and cut-off value showed similar results as expected and confirmed the findings outlined above (). We next explored the relationship between the lncRNA risk signature and the OS of LUAD patients. Both the training and validation groups showed that patients in the low-risk group exhibited longer OS than the high-risk group (). Moreover, the lncRNA risk signature showed excellent AUC values (>0.65 at 1, 3, and 5 years) in the ROC analysis in all cohorts ( and Figure S1), indicating effective prediction of survival by our lncRNA risk signature.

Figure 2

Validation of the 7 immune-associated lncRNA risk signature and survival analysis. (A) Prognostic analysis of the lncRNA risk signature in the validation set and combination set. (B) Survival curve of low- and high-risk groups in the training, validation, and combination set. (C) 3-year time-dependent ROC analysis of the lncRNA risk signature in the training, validation, and combination set.

Clinical correlation and independent risk factor analysis

We subsequently explored the underlying mechanisms by comparing the correlation between the lncRNA risk signature and clinical pathological characteristics. We observed that TNM stage, gender, and AJCC stage were significantly associated with the expression of the 7 lncRNAs in the combination set ( and Figure S2), which was particularly remarkable in AP000695.2, ITGB1-DT, LINC01150, and AL590226.1. ITGB1-DT expression was positively correlated with T stage (P<0.01), N stage (P<0.05), and AJCC stage (P<0.01) and its expression level was higher in males than in females (P<0.01). Interestingly, the expression of the four lncRNAs significantly associated with OS independently (). Kaplan-Meier analysis was performed after risk stratification by AJCC stage, gender, age, and TNM stage ( and Figure S2). Patients in the low-risk group showed improved OS compared with patients with high-risk for stage I/II (P<0.001), stage III/IV (P=0.02799), age ≤65 (P<0.001), age >65 (P<0.001), male (P<0.001), and female (P<0.001). Multivariate Cox analysis was performed to define the independent risk factors in the training, validation, and combination sets. The results suggested that in the combination set, AJCC stage and the lncRNA risk signature were independent prognostic factors that were significantly associated with OS. However, only the lncRNA risk signature was an independent prognostic factor in the training and validation set (). These data further support that the 7 immune-associated lncRNAs may be an independent predictor for the prognosis of LUAD patients.

Figure 3

Relationship between the expressions of the 7 lncRNAs and clinicopathological parameters and Kaplan-Meier survival curves in groups stratified by different clinical parameters. (A) The relationship between the expression of the 7 lncRNAs and TMN stage, gender, and AJCC stage. (B) The survival curve was plotted according to the expression levels of AP000695.2, ITGB1-DT, LINC01150, and AL590226.1 in the combination set. (C) The difference in survival curves between the high- and low-risk group stratified by age, stage, and gender. *, P<0.05, **, P<0.01, ***, P<0.001.

Table 2

Multivariate Cox regression analysis of clinicopathologic factors and immune-related lncRNAs signature for OS in training, validation and combined sets

Variable	Training set			Validation set			Combination set
	HR (95% CI)	P value		HR (95% CI)	P value		HR (95% CI)	P value
Age, years
≤65	1			1			1
>65	1.39 (0.93–2.09)	0.11		0.64 (0.34–1.21)	0.17		1.16 (0.84–1.61)	0.37
Gender
Female	1			1			1
Male	1.04 (0.67–1.57)			1.00 (0.54–1.85)	0.99		1.00 (0.72–1.40)	0.98
T stage
T1	1	0.86		1			1
T2	1.10 (0.65–1.84)	0.73		0.90 (0.45–1.82)	0.77		1.04 (0.69–1.55)	0.86
T3-T4	1.43 (0.69–3.01)	0.34		1.68 (0.51–5.51)	0.4		1.41 (0.78–2.56)	0.26
N stage
N0	1			1			1
N1	1.04 (0.49–2.19)	0.92		1.91 (0.51–7.16)	0.34		1.19 (0.66–2.15)	0.56
N2-N3	0.76 (0.22–2.55)	0.65		1.31 (0.23–7.54)	0.76		0.89 (0.34–2.32)	0.81
NX	1.37 (0.18–10.24)	0.76		12.049 (1.25–116.04)	0.03		2.20 (0.52–9.16)	0.28
M stage
M0	1			1			1
M1	0.76 (0.23–2.51)	0.65		0.52 (0.06–4.47)	0.55		0.79 (0.30–2.07)	0.64
MX	0.90 (0.56–1.44)	0.67		1.23 (0.60–2.57)	0.57		0.93 (0.63–1.37)	0.71
AJCC stage
Stage I	1			1			1
Stage II	1.71 (0.78–3.76)	0.18		2.66 (0.64–11.02)	0.18		2.01 (1.07–3.79)	0.03
Stage III-IV	3.32 (0.89–12.42)	0.07		4.95 (0.64–38.41)	0.13		3.33 (1.17–9.49)	0.02
Risk score
Low risk	1			1			1
High risk	2.61 (1.67–4.08)	<0.001		3.48 (1.86–6.51)	<0.001		2.69 (1.89–3.81)	<0.001

Low risk score was associated with highly active immune status

To investigate the immune status in the different risk groups, 29 immune-associated gene sets were analysed in the combination set. ssGSEA was performed to visualise the enrichment levels, functions, or pathways of immune cells in LUAD patients. We found that the low-risk group showed a separate cluster from the high-risk group with a higher ssGSEA score (). Similar results were obtained when the low- or high-risk group was scored by the ESTIMATE algorithm. We observed that the immune and stromal scores were significantly elevated in the low-risk group compared with the high-risk group (P<0.001, ). In contrast, tumour purity was significantly lowered in the low-risk group (P<0.001, ). We used PCA to evaluate the different distribution patterns in the immune-associated gene expression set and the whole genome expression set between low- and high-risk LUAD patients. The status of the patients in the low- and high-risk group was well distinguished in the immune-associated gene expression set. While the separation was not as clear when using the whole genome expression set, the two groups were still distinguishable (). These results showed that the 7 immune-associated lncRNAs may represent the immune status of the patient, in which the low-risk group harboured the highest number of immune and stromal cells in the TME whereas the high-risk group showed the highest number of tumour cells.

Figure 4

The low- and high-risk groups showed differential immune status in the combination set. (A) Overall immune status and tumour purity of low- and high-risk groups calculated by ssGSEA and ESTIMATE algorithm. (B) Comparison of the immune and stromal scores and ESTIMATE scores between low- and high-risk groups. (C) Comparison of tumour purity between low- and high-risk groups. (D) PCA analysis of the distribution patterns in low- and high-risk groups from the immune-associated gene expression set and whole-genome expression set. (E) Comparison of the expression levels of TLS signature and TLS hallmark genes between low- and high-risk groups. (F) Comparison of the expression levels of PD-L1, CTLA-4, and TIM-3 between low- and high-risk groups. **, P<0.01, ***, P<0.001.

Recently, tertiary lymphatic structures (TLSs) have gained scientific interest as a pathologic marker for immunotherapy in a variety of tumours (30-32). To investigate whether the risk score was related to TLSs in the LUAD immune micro-environment, we compared the expression level of the TLS signatures (CD79B, CD1D, CCR6, LAT, SKAP1, CETP, EIF1AY, RBP5, and PTGDS) (30) and TLS-hallmark genes (CCL19, CCL21, CXCL13, CCR7, CXCR5, SELL, and LAMP3) (33) between the low- and high-risk group (). These data suggested that TLSs may be more active in the low-risk group and may explain why patients in the low-risk group had a higher immune status.

In addition, we compared the expression level of PD-L1, CTLA-4, and TIM-3 in the two groups, which may be associated with a patient’s immunotherapeutic responsiveness. We found higher PD-L1, CTLA-4, and TIM-3 expression levels in the low-risk group compared with the high-risk group (P<0.1, P<0.01, and P<0.001, respectively, ). This indicated that patients in the low-risk group may respond better to immunotherapy than those in the high-risk group.

Association of cancer-related pathways and tumour mutation burden (TMB) with immune-associated lncRNA signature

We employed GSEA to discover the biological processes and pathways in the low- and high-risk groups in the combination set. As expected, the low-risk group showed significant enrichment in the gene sets related to the immune system in LUAD such as “IMMUNE_RESPONSE” (P=0.02, NSE =−1.94), “IMMUNE_SYSTEM_PROCESS” (P=0.004, NSE =−1.95), “KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION” (P=0.039, NSE =−1.68), and “GO_CYTOKINE_ RECEPTOR_ ACTIVITY” (P=0.002, NSE =−1.99). These data further underlined that a low-risk score may represent a higher immune status in LUAD. The high-risk group had the following significantly enriched gene sets: “GO_MITOTIC_NUCLEAR_DIVISION” (P<0.001, NSE =2.52), “GO_DNA_REPLICATION” (P<0.001, NSE =2.49), “GO_REGULATION_OF_SIGNAL_TRANSDUCTION_BY_P53_CLASS_MEDIATOR” (P<0.001, NSE =2.43), “KEGG_CELL_CYCLE” (P<0.001, NSE =2.32), and “KEGG_P53_SIGNALING_PATHWAY” (P=0.002, NSE =1.98, ), which have been reported to determine cancer cell division and proliferation. We observed unusually active signalling in tumour suppressor protein p53 (TP53). TP53 is the most commonly mutated gene in more than 50% of all human cancers including NSCLC (34). Next, we identified the somatic mutations in patients with LUAD and investigated the relationship between the risk score and the TMB. illustrates the variant classification, variant type, and SNV class and shows the variants per sample, variant classification summary, and top ten mutated genes in LUAD patients. TTN (41%), MUC16 (40%), RYR2 (34%), CSMD3 (34%), TP53 (47%), LRP1B (29%), USH2A (27%), ZFHX4 (27%), XIRP2 (22%), and KRAS (25%) were the top 10 mutant genes in LUAD, which were comparable to previous reports (34,35). TP53 was the most frequently mutated gene in LUAD and the main types of mutations were missense mutations (). TMB was also significantly higher in the high-risk group (P<0.001, ), which may explain why patients in the high-risk group had reduced OS compared with low-risk patients.

Figure 5

The low- and high-risk groups were associated with differential pathways and somatic mutation status. (A) GSEA analysis revealed differential enrichment of pathways in low- and high-risk groups. (B) The variant classification and type and SNV class in LUAC. (C) Variants per sample, variant classification summary, and top ten mutated genes in LUAC. (D) The mutant spectrum of the top 30 mutated genes with the largest number in LUAC. (E) Comparison of TMB levels between low- and high-risk groups. ***, P<0.001.

Prognostic value of immune-associated lncRNA risk signature may be applicable in other tumours

In the GSEA analysis outlined above, we found that “KEGG_PANCREATIC_ CANCER” (P=0.01, NSE =1.77) was also significantly enriched in the high-risk group (). Thus, we examined whether the prognostic value of the 7 immune-associated lncRNAs is also applicable in PDAC. We downloaded the expression of the 7 lncRNAs in PDAC samples from TCGA and established survival curves for the low- and high-risk groups. PDAC patients with low-risk scores had a significantly longer OS than high-risk patients (P<0.001, ). ssGSEA and ESTIMATE analyses revealed that low- and high-risk groups were segregated (). The low-risk group showing higher ssGSEA score and proportion of immune and stromal cells, which was consistent with our results in LUAD samples (P<0.001, ). However, unlike in LUAD patients, the mutation variants per sample in PDAC were lower (33.5/Mb vs. 140/Mb) (), although similar results were obtained when comparing the TMB between the low- and high-risk groups (P<0.001, ). These data indicated that the 7 immune-associated lncRNAs may be able to predict patient survival in other tumours as well.

Figure 6

Validation of the 7 immune-related lncRNA risk signature in PDAC samples. (A) GSEA analysis revealed that high-risk patients were associated with “KEGG_PANC REATIC_ CANCER” (P=0.01, NSE =1.77). (B) Survival curve of low- and high-risk groups in PDAC patients. (C) Overall immune status and tumour purity of low- and high-risk groups analysed by ssGSEA and ESTIMATE algorithm in PDAC patients. (D) Comparison of the immune and stromal scores, ESTIMATE scores, and tumour purity between low- and high-risk groups. (E) The variant classification and type, SNV class, variants per sample, variant classification summary, and top ten mutated genes in PDAC. (F) TMB levels of low- and high-risk groups in PDAC. ***, P<0.001.

Discussion

NSCLC is one of the most common causes of cancer morbidity and mortality globally, in which LUAD accounts for a large proportion (1,2). Despite advances in treatment methods in recent years, the prognosis of patients remains poor (3-5). Immunotherapy has been regarded as the most promising treatment since its introduction to the clinic. However, determining the eligibility of a patient is still an unsolved issue (6-8). One of the causes is the lack of effective biomarkers or risk models to guide physicians (7,8). In the recent years, lncRNA has attracted a significant amount of attention due to its wide range of biological effects in carcinogenesis and tumour progression (13,15,16). The role of lncRNA in immune regulation has also been reported. For instance, LncNKILA enhances T cell sensitivity by NF-κB-induced apoptosis and knockdown lncNKILA effectively attenuates tumour-specific cytotoxic T cells and improves patient survival in lung cancer (36). Therefore, immune-associated lncRNAs can be used to distinguish different immune statuses and predict a patient’s prognostic risk. In this study, we constructed a novel model to predict prognosis and survival using immune-related lncRNA.

A total of 788 immune-associated lncRNAs were derived from the expression of 332 immune-related genes in LUAD samples from TCGA based on Pearson’s correlation analysis. Univariate and multivariate Cox regression and LASSO regression analysis was performed to construct an immune-associated lncRNA signature consisting of 7 lncRNAs, which was able to determine the prognosis of samples in the training, validation, and combination group. 1-, 3-, and 5-year ROC analyses also demonstrated the feasibility of our prognostic lncRNA risk signature. Through stratified analysis, we further found that the immune-associated lncRNA risk signature may be used to distinguish the survival outcomes of patients with different clinical variables. Moreover, multivariate Cox analysis revealed that the lncRNA risk signature was an independent prognostic factor. These 7 immune-associated lncRNAs were first reported in lung cancer and AP000695.2, ITGB1-DT, LINC01150, and AL590226.1 appeared to be closely related to clinical characteristics and independently predicted survival, although further investigations are required. These data revealed that immune-associated lncRNAs may be good prognostic indicators.

To investigate the TME in LUAD, ssGSEA and ESTIMATE analyses were performed to evaluate the immune and stromal status as well as tumour purity. We observed elevated immune activity in the low-risk group compared with the high-risk group, resulting in a significantly higher immune and stromal score while the high-risk group had higher tumour purity, suggesting that innate and adaptive immunity in the TME of low-risk LUAD patients were more aggressive. These findings were confirmed in subsequent GSEA analysis. Combining these results with clinical parameters, we postulate that the expression of immune checkpoint genes is a necessary indicator for screening patients before immunotherapy and that TLS is the most promising biomarker for guiding immunotherapy (28-31). We examined the differences in the expression of these genes including PD-L1, CTLA4, TIM3, and TLS signature and TLS hallmark genes (6-8) in the low- and high-risk groups. The expression of most of these genes were higher in the low-risk group, suggesting that the low-risk group may better respond to immunotherapy. The data above revealed that according to the risk model of the 7 immune-associated lncRNAs, a higher proportion of immune cells and stromal cells were present in the TME. ssGSEA and ESTIMATE analyses indicated a stronger immune-related response, which was reflected in the high expression of signals for immune-related pathways and TLS-related genes in tissues of low-risk patients, indicating that low-risk patients with LUAD may be more responsive to immunotherapy.

Somatic mutations have been confirmed to be a key aspect of carcinogenesis and cancer development in previous studies (37-39). Gao et al. demonstrated that the expression of some lncRNAs is affected by somatic mutations and were commonly downregulated while carrying low mutation frequencies and non-silent mutations in most cancer types, thus influencing the molecular pathogenesis (38). In our study, we found that the high-risk group was more likely to harbour enriched pathways related to cancer cell division and proliferation. In addition, the signalling pathway of the most commonly mutated gene TP53 was also significantly related to the high-risk group. Consequently, we also found that TMB was significantly higher in the high-risk group. We speculate that the immune-related lncRNA risk signature can identify patients with a highly active immune status and low TMB in LUAD and such patients have a better prognosis than other patients.

Extrapolation of the immune-related lncRNA risk signature from LUAD to other cancers may further confirm the role of immune-associated lncRNAs in tumour progression and prognosis. In this study, we investigated the applicability of our lncRNA risk signature model in PDAC samples and observed similar consistency as in LUAC samples in predicting prognosis, immune status, and TMB status in low- and high-risk groups. Our study is in agreement with the work of Ruess et al. (40), who collectively refer to some NSCLC and PDAC as mutant KRAS-driven tumours. In our study, KRAS was in the top ten mutated genes in both LUAD and PDAC samples and accounted for 25% and 77%, respectively. Further exploration of the relationship between KRAS and immune lncRNA is required. Taken together, we demonstrated that our immune-related lncRNA risk signature may be of prognostic value in certain types of tumour patients based on its ability to assess the immune status and TMB. Further studies are required to define this special type of tumours.

In conclusion, we constructed an immune-related lncRNA signature composed of 7 lncRNAs that indicated significant differences in immune microenvironment and somatic mutations in LUAD patients. This may play an important role in predicting the prognosis of patients and providing information related to immunotherapy effectiveness. We provided a new perspective on the role of lncRNA in the development and the molecular mechanism of immunity in LUAD. Through the preliminary exploration of the immune-related lncRNA risk signature in PDAC, we suggest the extension of the use of this risk signature to pancreatic cancer in the future. There are certain limitations in this study. Firstly, the results of the study were based on TCGA database, which lacked detailed information on resection extent, subtypes of LUAD, radiotherapy, and chemotherapy. Secondly, the 7 immune-associated lncRNAs have not been reported before and the mechanisms through which the prognostic immune-associated lncRNAs modulate the progression of LUAD require further investigation.