Construction of an immune-related lncRNA signature as a novel prognosis biomarker for LUAD

Data acquisition

We obtained the LUAD transcriptome profiles and the corresponding clinical information of LUAD patients and normal samples from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) database in December 2020. We excluded samples with repeated data and those with a follow-up time of 0 days. We downloaded gene transfer format (GTF) files from Ensembl (http://asia.ensembl.org) to differentiate lncRNAs and mRNAs. We also downloaded several immune-related genes from The Immunology Database and Analysis Portal (ImmPort, https://www.immport.org/). We calculated the correlation coefficient between the immune-related genes and lncRNAs to identify irlncRNAs. The thresholds were set as |cor| > 0.4 with p < 0.001. To identify the differentially expressed irlncRNAs (DEirlncRNAs), we performed the R-x64-4.0.3 limma package for irlncRNA differential expression analysis (log fold change [FC] > 1, false discovery rate [FDR] < 0.05).

Matching DEirlncRNA pairs

The DEirlncRNAs were singly cyclically matched and a 0- or -1 matrix was constructed supposing A was a DEirlncRNA pair composed of lncRNA X and lncRNA Y; A was defined as 0 if the expression level of lncRNA X was lower than that of lncRNA Y; otherwise, A was defined as 1. Next, the matrix was further filtered. Regardless of whether its expression was 0 or 1, it could not directly predict the patient’s prognosis. When the number of DEirlncRNAs pairs with 0 or 1 expression accounted for <20% of the total pairs, they were excluded.

Construction of DEirlncRNA pairs and prognostic signature

Univariate Cox analysis was performed to screen the DEirlncRNA pairs whose expression levels were significantly related to the overall survival (OS) of the LUAD patients (p < 0.05). The modified least absolute shrinkage and selection operator (LASSO) regression was used to identify the DEirlncRNA pairs most relevant to OS by running the R-x64-4.0.3 glmnet package. The modified Lasso regression was run with 10-fold cross-validation (p < 0.05) and was set 1000 cycles with 1000 random stimulation for each cycle. Subsequently, the occurrence frequency of every pair in the Lasso regression was counted, and pairs with >100 frequency were chosen for subsequent Cox proportional hazard regression analysis. Then, the receiver operator characteristic curve (ROC) was plotted and the area under the curve (AUC) was counted. The model’s 1-, 2-, and 3-year ROC were plotted as well. The risk score with the constructed prognosis signature was calculated for every clinical sample with the formula as follows:

$$\begin{array}{c}\text{n}\\ \text{Risk score }=\sum \text{Coef }\left(\text{i}\right)\times \text{E}\left(\text{i}\right),\\ \text{i }=\text{1},\end{array}$$

Coef (i) and E(i) represent the regression coefficient of the multivariate Cox analysis for the DEirlncRNA pairs and the value of each DEirlncRNA pair expression, respectively. The maximum inflection point of the 1-year ROC curve was identified by determining the maximum Akaike information criterion (AIC) values, and it was recognized as the cut-off point for dividing each sample into different risk groups. The value of Youden index (Sensitivity+Specificity−1) is the largest at the maximum inflection point.

Validation of the risk prognosis model

Kaplan–Meier log-rank test was performed for comparing the survival between different risk groups to evaluate the predictive value of the prognostic signature. The specific risk scores of every patient were calculated by R-x64-4.0.3 packages survival, glmnet, survivalROC, and survminer. Several chi-square tests were performed to explore the correlation between the risk score and the clinicopathological characteristics. The figures was labeled as follows for visualization: ***<0.001, **<0.01, and *<0.05. The risk score differences in clinicopathological characteristics among different groups were calculated performing the Wilcoxon signed-rank test. Univariate and multivariate Cox regression analyses were used to validate whether the constructed model was an independent prognostic indicator for LUAD patients. The results were demonstrated with a forest map. The R-x64-4.0.3 packages utilized here were survival, pHeatmap, and ggupbr.

Correlation analysis of tumor immune microenvironment

Considering that the six DEirlncRNA pairs were associated with tumor immunity, we calculated the correlation coefficient between the risk prognosis model and the LUAD-infiltrating immune cells based on TIMER, xCELL, quanTIseq, MCPcounter, EPIC, CIBERSORT−ABS, and CIBERSORT software. The correlation coefficients are shown in a lollipop diagram. The differences in tumor-infiltrating immune cell abundance between different risk groups were explored by Wilcoxon signed-rank test. Spearman correlation analysis was used to calculate the correlation coefficients between the risk scores and the tumor-infiltrating immune cells using the R-x64-4.0.3 ggplot2 package. To further explore the TME of patients with LUAD, we calculated the patients’ StromalScore, ImmuneScore, and ESTIMATEScore (StromalScore + ImmuneScore) using the R-x64-4.0.3 limma and estimate packages. Next, the differences in the stromal cell and immune cell content around the tumor between different risk groups were calculated by the Wilcoxon signed-rank tests. Column diagrams were drawn with GraphPad Prism 8.0.2 for visualization.

Gene set enrichment analysis (GSEA)

We used GSEA 4.0.1 (https://www.gsea-msigdb.org/gsea/index.jsp) [14] to screen functional pathways in different risk groups. The two risk groups were performed on Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets. The thresholds were set as p < 0.05 and FDR < 0.25, a random sample permutation was set for 1000. All other settings were set according to their default values.

Expression analysis of immune checkpoint inhibitor (ICIs)-related immunosuppressive molecules

We used the R-x64-4.0.3 packages ggpubr and limma to analyze the differences in ICI-related immunosuppressive molecules’ expression level between different risk groups of the constructed model, and visualized the result with a violin plot.

The significance of the signature for guiding medication

National Comprehensive Cancer Network guidelines recommend anti-tumor drugs such as cisplatin, paclitaxel, gemcitabine, gefitinib, and erlotinib for treating lung cancer. To better guide our clinical medication, we calculated the median inhibitory concentration (IC₅₀) of common anti-tumor drugs used for treating LUAD. Wilcoxon signed-rank test was performed to analyze the difference in the IC₅₀ between different risk groups. We used the R-x64-4.0.3 pRRophetic and ggplot2 packages and visualized the results with several box drawings.

Consent for publication

Written informed consent for publication was obtained from all participants.

Availability of data and materials

All data analysed during the current study are accessible from the TCGA database (https://portal.gdc.cancer.gov/).

Identification of DEirlncRNAs

We identified DEirlncRNAs using a multi-step approach (Figure 1). We downloaded LUAD transcriptome data and corresponding clinical data from TCGA-LUAD cohort (54 normal samples and 497 LUAD samples). Then, we annotated the mRNAs and lncRNAs using Ensembl GTF files, and downloaded the immune-related genes. Next, we identified the irlncRNAs from the list of immune-related genes using co-expression analysis. Finally, 517 DEirlncRNAs were identified using R-x64-4.0.3 language (Figure 2A): 72 were downregulated and 445 were upregulated (Figure 2B).

Figure 1. The work flowchart for the study.

Figure 2. Identification of DEirlncRNAs. (A, B) The differential co-expression analysis identified 517 DEirlncRNAs. (C) Forest map shows the six DEirlncRNA pairs screened by Cox proportional hazard regression.

Construction of the prognostic signature

Using an iteration loop and a 0-or-1 matrix screening among 517 DEirlncRNAs, we identified 93,144 valid DEirlncRNA pairs. Univariate Cox analysis was first performed (Figure 2C). Subsequently, modified Lasso regression analysis (Figure 3A, 3B) yielded eight DEirlncRNA pairs; six were included for constructing the risk assessment model. To assess the model’s accuracy, we plotted an ROC for validation (Figure 3C). The curve indicated that the AUC for this model was 0.762. The 1-, 2-, and 3-year AUC values were 0.762, 0.734, and 0.746, respectively (Figure 3E). Compared with other clinical characteristics such as age, sex, and stage, the 1-year AUC value was the maximum (Figure 4A). According to the AIC values, the maximum inflection point on the 1-year ROC curve was recognized as the cut-off point (Figure 3D). Based on the cut-off point, the LUAD patients were distinguished into different risk groups.

Figure 3. Validation of the 6-DEirlncRNA pair signature. (A, B) The modified LASSO penalized regression analysis identified eight DEirlncRNA pairs most related to prognostics. (C, D) The ROC showed the optimal cut-off point based on AIC value. (E) The 1-, 2-, 3-year AUC values were >0.7. (F) High-risk patients exhibited worse OS than the low-risk patients. (G, H) Univariate (G) and multivariate (H) Cox regression analysis of risk score, age, sex, and tumor-node-metastasis (TNM) stage indicating that risk score (HR = 1.467 [1.343–1.603], p < 0.001) and stage (HR = 1.422, [CI] = [1.222–1.654], p < 0.001) were related to survival.

Figure 4. Validation of the 6-DEirlncRNA pair signature. (A) Calculation of the AUC for risk score, age, sex, and TNM stage shows that risk score is the maximum. (B, C) The risk curve (B) and scatter plot (C) indicates that patients with high risk score had worse survival. (D–G) Scatter diagrams and strip chart (G) show that N stage (D), clinical stage (E), and T stage (F) were significantly related to the risk score.

The prognostic signature is a powerful LUAD prognostic indicator

Based on the above division, there were 286 low-risk patients and 178 high-risk patients. Survival analysis (Figure 3F) indicated that the patients with high scores had a worse prognosis than patients in low-risk group (p < 0.001). We investigated the relationship between risk scores and clinicopathological characteristics with chi-square tests and generated a heat map. To assess the meaning of the prognosis model in clinical application, we compared the constructed model with clinicopathological characteristics using univariate (Figure 3G) and multivariate (Figure 3H) Cox analyses and ROC analysis. The results suggested that two factors, i.e., the risk score (hazard ratio [HR] = 1.467, confidence interval [CI] = 1.343–1.603, p < 0.001) and stage (HR = 1.422, [CI] = 1.222–1.654, p < 0.001), correlated with the survival. This model could also predict the survival status (AUC = 0.762; 1-year AUC = 0.762; 2-year AUC = 0.734; 3-year AUC = 0.746). The above results show that the prognostic signature can serve as a potential independent prognostic predictor.

We generated a risk curve (Figure 4B) and a scatter plot (Figure 4C) to visualize the survival status and risk score of each patient with LUAD, which indicated that patients in the high-risk group had higher risk coefficient and mortality compared to low-risk patients. Then, we investigated the correlation between the constructed model and clinicopathological characteristics using chi-square testing. The scatter diagrams and consequent strip chart (Figure 4G) indicated that N stage (Figure 4D), clinical stage (Figure 4E), and T stage (Figure 4F) were associated with the risk score significantly.

Correlation between LUAD prognostic signature and the tumor immune microenvironment

Figure 5A showed that the B cell, CD4+ T cell, M2 macrophage, neutrophil, and mast cell correlation values with the risk score were -0.368, -0.230, -0.333, -0.238, and -0.282, respectively, indicating that the abundance of these tumor-infiltrating immune cells correlated negatively with the LUAD risk score (Figure 5B–5D). We were also able to conclude that the StromalScore (Figure 6A), ImmuneScore (Figure 6B), and ESTIMATEScore (Figure 6C) associated negatively with the LUAD risk score and are all statistically significant.

Figure 5. Correlation between the LUAD 6-irlncRNA pair signature and immune cell subtype infiltration. (A) Estimation of the infiltrating immune cells by TIMER, CIBERSORT, xCELL, quanTIseq, MCPcounter, EPIC, and CIBERSORT-ABS. (B–D) High-risk patients exhibited higher content of B cells (B), CD4+ T cells (C), and M2 macrophages (D) compared to the low-risk patients.

Figure 6. Assessment of TME and GSEA of high- and low-risk patients. (A–C) Column diagrams show that the high-risk group had significantly lower StromalScore (A), ImmuneScore (B), and ESTIMATEScore (C) than the low-risk group. (D) GSEA results reveal that intestinal immune network for IgA production and the T and B cell receptor signaling pathway were enriched in the low-risk groups. In the high-risk groups, cell cycle, P53 signaling, DNA replication, adherens junction, actin cytoskeleton regulation, pathways in cancer, and TGF-β signaling pathways were active relatively.

Enrichment of the functional phenotypes of the irlncRNA signature

GSEA of the prognosis model (Figure 6D) revealed that several pathways were enriched in low-risk patients: the IgA production in intestinal immune network and the T and B cell receptor signaling pathway. However, cell cycle, P53 signaling, DNA replication, adherens junction, actin cytoskeleton regulation, pathways in cancer, and TGF signaling were enriched in high-risk patients.

The risk assessment model is a potential predictor for guiding clinical medication

In recent years, ICIs have been used in LUAD treatment. We determined whether there were statistical differences in ICI-related biomarkers expression between different risk groups, and discovered that a high risk score correlated negatively with high CTLA-4 (p < 0.001, Figure 7A), HAVCR2 (p < 0.05, Figure 7B), and PD-1 (p > 0.05, Figure 7D) expression, and correlated positively with high PD-L1 expression (p > 0.05, Figure 7C); There was no statistical difference in PD-1 and PD-L1 between the high- and low-risk groups. We also calculated the IC₅₀ of common anti-tumor drugs in different risk groups and discovered that patients with low risk scores were correlated with a higher IC₅₀ of chemotherapeutics such as cisplatin (p < 0.001, Figure 7E), paclitaxel (p < 0.001, Figure 7F), and gemcitabine (p < 0.001, Figure 7G) and targeted therapy such as erlotinib (p < 0.001, Figure 7H) and gefitinib (p < 0.001, Figure 7I), indicating that the prognostic signature acted as a promising predictor of anti-tumor drug sensitivity.

Figure 7. Clinical application of the 6-DEirlncRNA pair signature. (A–D) High risk score correlated negatively with high HAVCR2 (B) and PD-1 (D) expression, and correlated positively with high CTLA-4 (A) and PD-L1 (C) expression; There was no statistical difference in PD-1 and PD-L1 between different risk groups. (E–I) The signature could serve as a promising biomarker of anti-tumor drug sensitivity: low risk scores were associated with higher IC₅₀ for cisplatin (E), paclitaxel (F), gemcitabine (G), erlotinib (H), and gefitinib (I).