Construction of a prognostic immune-related lncRNA model and identification of the immune microenvironment in middle- or advanced-stage lung squamous carcinoma patients.

Background: Globally, non-small-cell lung cancer (NSCLC) has a high incidence, and NSCLC patients have poor prognoses. Lung squamous carcinoma (LUSC) is a major pathological type of NSCLC. LncRNAs play important roles in tumor progression and immune system functions. The aim of this study was to construct a predictive model with immune-related lncRNAs and to assess the immune microenvironment in middle- or advanced-stage LUSC patients.
Methods: RNA sequencing data and corresponding clinical LUSC data were downloaded from The Cancer Genome Atlas. Immune genes were obtained from the Molecular Signatures Database. Immune-related lncRNAs were identified by Pearson correlation analysis in R. The model was constructed using univariate and multivariate Cox regression analyses. Finally, we validated the prognostic immune-related lncRNA model in a cohort from the Fudan University Shanghai Cancer Center.
Results: Our risk model included four immune-related lncRNAs (LINC00944, AL034550.2, AC020907.1 and AC027682.6). Survival analysis revealed that overall and disease-free survival were shorter in the high-risk group than in the low-risk group. Independent prognostic analysis showed that our model could be used as an independent prognostic predictor. The high-risk group was positively associated with CD8+ T cells, B cells, myeloid dendritic cells, macrophages, regulatory T cells (Tregs) and cancer-associated fibroblasts and high expression of PD1 and CTLA4. Additionally, a low-risk score was correlated with lower half maximal inhibitory concentrations (IC50s) of cisplatin, docetaxel, vinorelbine and paclitaxel and a higher IC50 of gemcitabine. Gene set enrichment analysis suggested that these lncRNAs may participate in tumor progression and immune processes. Validation with the clinical cancer cohort demonstrated that higher risk scores were associated with a higher, but not statistically significant, likelihood of recurrence.
Conclusion: We established a risk score model including four immune-related lncRNAs. The model accurately predicts the prognosis of middle- or advanced-stage LUSC patients and provides an important reference for individualized treatment.

Introduction

Globally, lung cancer is still one of the most commonly diagnosed cancers (11.4%), second only to breast cancer in women (11.7%), and lung cancer remained the leading cause of cancer-related death (18%) in 2020 [1]. There are mainly two pathological categories of lung cancer: non-small-cell lung cancer (NSCLC, approximately 85%) and small cell lung cancer (SCLC, approximately 15%). The World Health Organization classifies NSCLC into three main types: adenocarcinoma, squamous carcinoma and large cell carcinoma, of which squamous carcinoma accounts for approximately 20–30% [2, 3]. Most patients are already at an advanced stage when discovered due to the dearth of classic signs and symptoms in the early stage, and the 5-year relative survival rate for patients with stage IV disease is only 4% [4]. This outcome shows that the prognoses of middle- or advanced-stage lung squamous carcinoma (LUSC) patients are poor.

In the past decennary, the importance of tumor microenvironment (TME) during the periods of initiation and development of lung cancer has been recognized [5, 6, 7]. The presence of innate immune cells and adaptive immune cells in the TME is conducive to the tumor growth and metastasis process [8]. The introduction of immune checkpoint inhibitors (ICIs), including antibodies against programmed death-1/death-ligand 1 (PD-1/PD-L1) and monoclonal antibodies against cytotoxic T-lymphocyte antigen-4 (CTLA4), ushered in a new era of immunotherapy for lung cancer patients. Many patients with NSCLC have been able to benefit from this therapy [9, 10]. However, in NSCLC, the predictive power of PD-L1, which is currently regarded as the most robust biomarker, is also insufficient in many situations, and the need for further, more complex biomarkers is immense [11]. Therefore, we hoped to select biomarkers appropriate for constructing a prognostic model for LUSC patients with middle- or advanced-stage disease and to assess their immune microenvironment.

Long noncoding RNAs (lncRNAs) are a kind of RNA with lengths longer than 200 bases that cannot be translated into proteins [12]. It has been demonstrated that a comprehensive and complicated position was taken by lncRNAs in regulating cancer initiation and development, and lncRNAs are extremely expected to serve as new biological markers and therapeutic targets [13, 14, 15, 16]. In addition, lncRNAs are essential and crucial in regulating immune gene expression as well as closely associate with the TME [14, 17, 18, 19]. Accordingly, it is worthwhile to construct an immune-related lncRNA model to predict the prognosis of LUSC patients with middle- or advanced-stage disease and to assess their immune microenvironment, which can contribute to the development of individualized treatment plans.

In our study, we comprehensively analyzed immune-related lncRNA data from 175 stage IIB-IV LUSC patients downloaded from The Cancer Genome Atlas (TCGA) to establish a prognostic model and identify immune cell infiltration in tumor tissues.

Materials and methods

Data collection and preprocessing

The RNA sequencing (RNA-seq) data of LUSC patients and the corresponding clinical data were downloaded from The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga/) in April 2021. GTF files were downloaded from Ensembl (http://asia.ensembl.org) for gene annotation to classify mRNAs and lncRNAs. We used the “DESeq2” package in R (4.0.5) to screen for differentially expressed mRNAs and lncRNAs (fold change = 1.0, P value <0.05), and the data were normalized by variance-stabilizing transformation (VST) [20]. The list of immune genes was obtained from the Molecular Signatures Database (MSigDB) and cross-referenced with the list of differentially expressed mRNAs. Then, the immune-related lncRNAs were identified by an analysis of the Pearson correlation of immune genes and differentially expressed lncRNAs with R version 4.0.5 (cor-Filter > 0.4 and P value <0.001). Finally, we selected stage IIB-IV patients with survival data >30 days (n = 175) and integrated their clinicopathological characteristics and survival information into the RNA-seq data.

Development of the immune-related lncRNA risk score model

Immune-related lncRNAs associated with prognosis were screened by univariate Cox regression analysis (P value <0.01). The screened lncRNAs along with clinicopathological features were included in multivariate Cox regression analysis to minimize the Akaike information criterion (AIC) value to obtain the best-performing prognostic model. Cox regression analysis was performed with the package “survival” in R (4.0.5). Finally, we identified four immune-related lncRNAs. According to the expression level and regression correlation coefficients (β) of these four lncRNAs, risk scores were calculated for each LUSC patient using the following equation: risk score = - 0.136673 × ExpressionAC020907.1 + 0.293362 × ExpressionAC027682.6 + 0.301214 × ExpressionAL034550.2 + 0.120110 × ExpressionLINC00944.

Identification and evaluation of the prognostic immune-related lncRNA model

The median risk score was used to divide LUSC patients into high- and low-risk groups. We assessed the accuracy of the risk score model with survival analysis and receiver operating characteristic (ROC) curves in R (4.0.5), using the “survival” and “survivalROC” packages. Then, we determined whether the model performance was influenced by other clinical characteristics and whether the model could be considered as an independent prognostic indicator with the “survival” package in R (4.0.5). The scatter plot combined with the heatmap shows the survival status and the differential expression of these four lncRNAs in the two risk subgroups. Moreover, we further explored whether the clinical characteristics and the expression of immune checkpoint molecules differed between the two subgroups. The immune subtypes of patients were determined by methods detailed in Vesteinn Thorsson et al [21]. Statistical analyses were performed with the chi-square test in Stata16 and the Wilcoxon test in R (4.0.5). A P value <0.05 was considered statistically significant. Finally, we applied the “scatterplot3d” and “limma” packages in R (4.0.5) for principal component analysis (PCA) and visualization of the prognostic model results.

GSEA and immune infiltration

We used gene set enrichment analysis (GSEA 4.1.0) to identify biofunctional differences between the two risk subgroups. We considered pathways with normalized enrichment score (NES) absolute values >1, nominal (NOM) p-val < 0.05, and false discovery rate (FDR) q-val < 0.25 to be significantly enriched. Immune cell infiltration estimation data from LUSC patients assessed by several mainstream methods (TIMER [22], CIBERSORT [23], quanTIseq [24], xCell [25], MCP-counter [26] and EPIC [27]) were downloaded from TIMER 2.0 [28] (http://timer.cistrome.org/). We calculated the correlation between the risk scores and the presence of each type of immune cell identified in the aforementioned methods and displayed their differences across the two subgroups. The Wilcox test and Spearman correlation test were performed in R (4.0.5). A P value <0.05 was considered statistically significant.

Predicting clinical chemotherapeutic response

Responses to clinical chemotherapeutic agents by these 175 LUSC patients, whose data obtained from the TCGA database, were predicted using the pRRopheticPredict() function in the “pRRophetic” package [29, 30], and statistical analysis was performed by t test or Wilcox test with a P value less than 0.05 as the threshold using R (4.0.5). The chemotherapeutic agents were selected according to the National Comprehensive Cancer Network (NCCN) guidelines [31].

Validation of our prognostic model in the Fudan University Shanghai Cancer Center cohort

We collected data from 14 LUSC patients received surgical resection from 2017 to 2018 at the Fudan University Shanghai Cancer Center. All enrolled patients did not undergo any chemotherapy, radiotherapy, or other treatments prior to surgery. Histopathological examination revealed that all enrolled patients had LUSC with stage IIB-IV disease. We detected the expression level of the four identified lncRNAs by quantitative real-time PCR and plotted a disease-free survival (DFS) curve. All samples were obtained from the Department of Biological repositories of the Fudan University Shanghai Cancer Center. All enrolled patients signed informed consent forms for the use of tissue upon admission. The study was supported by the Ethics Committee of Fudan University Shanghai Cancer Center (ethics number 050432-4-2108∗).

RNA extraction and quantification

We extracted total RNA of these 14 patients with TRIzol Reagent (Ambion), then reverse transcribed it into cDNA with RT SuperMix (+gDNA wiper) kit (Vazyme, R323-01). Quantitative real-time PCR was performed with SYBR qPCR Master Mix (Vazyme, Q712) and run on an Applied Biosystems® 7300Plus Real-Time PCR Instrument. Finally, we selected β-actin as an endogenous reference genes to calculate the relative expression level of the four lncRNAs in each sample using the 2-△Ct method. All PCR primers were directly synthesized (Sangon Biotech, Shanghai), and the sequences are shown in Table 1.

Table 1

qRT–PCR primer sequences.

LncRNA	Sequence
AC020907.1	Forward: TGCTGAAGTAACGAACTAACCTGAA
	Reverse: GGTGGTGATGGTGACGGTAATG
AC027682.6	Forward: TGCCATTTCAGCCAGCCTCAG
	Reverse: TGCCACCACCACCTTCAGGAA
AL034550.2	Forward: TGCTTCAGATTCAATGGATGCTACT
	Reverse: GCCCGAGGATGGAAACAGACT
LINC00944	Forward: CCTCTTAATCCTCTGTCCTCCATC
	Reverse: CTCTCCAGTGTTATGAAGTTCAAGT

Results

Acquisition of immune-related lncRNAs

We analyzed the mRNA and lncRNA expression profiles of 502 LUSC samples and 49 normal samples based on set thresholds (fold change >1 and P value <0.05) and obtained 7148 and 4748 differentially expressed genes (DEGs), respectively (Figure 1 and Supplementary Data 1). Subsequently, we downloaded the “immune system process” and “immune response” gene sets from MSigDB and cross-referenced them with the list of differentially expressed mRNAs. Immune-related genes were identified. Finally, we performed Pearson correlation analysis of differentially expressed lncRNAs and immune-related genes (Pearson correction coefficient >0.4 and P value <0.001), and obtained 755 immune-related lncRNAs.

Figure 1

Volcano plots of differentially expressed mRNAs (A) and lncRNAs (B) in LUSC versus normal samples from the TCGA database. Red dots are upregulated genes, blue dots are downregulated genes, and gray dots are no significantly different genes.

Construction and evaluation of an immune-related lncRNA risk score model

We selected stage IIB-IV LUSC patients (n = 175). A total of 8 immune-related lncRNAs were significantly associated with patient prognosis by univariate Cox regression analysis (P value <0.01), and the results are presented in a forest plot (Figure 2A). By including these 8 lncRNAs and patient clinicopathological features in the multivariate Cox regression analysis, a risk model of 4 immune-related lncRNAs was eventually established (Table 2). We calculated the risk score of each LUSC patient according to the following formula: risk score = - 0.136673 × ExpressionAC020907.1 + 0.293362 × ExpressionAC027682.6 + 0.301214 × ExpressionAL034550.2 + 0.120110 × ExpressionLINC00944. We divided the samples into a high-risk group (n = 87) and a low-risk group (n = 88) based on the median risk score.

Figure 2

A Forest plot of 8 immune-related lncRNAs associated with prognosis derived from univariate Cox regression analysis. Red and green indicate risk factors and protective factors, respectively. B-D Survival analysis showed that the higher the risk score was, the worse the prognosis. (B) Kaplan–Meier overall survival (OS) curve (P = 0.00019). (C) Kaplan–Meier disease-free survival (DFS) curve (P = 0.00030). (D) Scatterplot of survival status and risk score (P = 0.00004). E The risk score model of these 4 immune-related lncRNAs in LUSC patients. Risk score (top); survival status (middle); heatmap (bottom). F, G Forest plots of univariate independent prognostic analysis (F) and multivariate independent prognostic analysis (G). Red and green indicate risk factors and protective factors, respectively. H, I ROC curve analysis. (H) The area under the ROC curve (AUC) indicated the highest accuracy of our risk model compared to those of other clinical characteristics. (I) Time-ROC curve analysis of our risk model at 250, 500, 750, 1000, 2000, and 5000 days.

Table 2

The coefficients, hazard ratios (HRs) and P values of four immune-related lncRNAs that contributed to the risk model.

LncRNAs	Coef	HR (95% CI)	P value
AC020907.1	-0.136673	0.872256 (0.779413–0.976158)	0.0173
AC027682.6	0.293362	1.340928 (1.000794–1.79666)	0.0494
AL034550.2	0.301214	1.351499 (1.028015–1.776772)	0.0309
LINC00944	0.120110	1.12762 (0.966004–1.316277)	0.1281

Coef, regression coefficient; HR, hazard ratio; CI, confidence interval.

Not surprisingly, the overall survival (OS) and DFS of patients in the high-risk group were significantly lower than those in the low-risk group (Figure 2B-D). Survival analysis revealed that the survival rate of patients in the high-risk group was 69% after 1 year, 41% after 3 years, and 26% after 5 years, whereas the survival rate of patients in the low-risk group was 88% after 1 year, 60% after 3 years, and 56% after 5 years. We ranked the risk scores of these 175 LUSC patients by combining the survival outcomes and expression of these four lncRNAs (Figure 2E). The findings demonstrated that the expression of AC027682.6, LINC00944 and AL034550.2 was upregulated in the high-risk group, while that of AC020907.1 was downregulated, and the mortality rate rose with the increasing risk score.

To determine whether the model was independent of other clinical factors, such as age, sex, T stage, N stage, M stage and American Joint Committee on Cancer (AJCC) stage, we performed univariate and multivariate independent prognostic analyses. The outcomes (Figure 2F, G) suggested that only the risk score was correlated with OS, with statistically significant differences (P value <0.001). Therefore, the signature of four immune-related lncRNAs we identified can be used as an independent prognostic predictive factor.

The area under the ROC curve (AUC) was used to evaluate the accuracy of prediction. As shown in Figure 2H and Figure 2I, the signature of four immune-related lncRNAs we identified performed best (AUC = 0.737) compared with the predictive performances of other clinical characteristics, such as age, sex, T stage, N stage, M stage and AJCC stage. Time-ROC curve analysis suggested a robust predictive effect of the signature we identified.

Analysis of clinical features between the different risk subgroups

A series of chi-square tests were carried out to explore the differences in clinical features between the different risk subgroups, and the results were presented in Figure 3. There were more females and higher mortality in the high-risk group than in the low-risk group, with statistically significant differences (P value <0.05). Regarding immune subtypes, the high-risk group mainly exhibited C1 (wound healing) and C2 (IFN-γ dominant) types, whereas the low-risk group was more likely to exhibit the C1 (wound healing) type (P value <0.05).

Figure 3

Assessment of clinical characteristics grouped based on risk score. ∗P < 0.05, ∗∗P < 0.01.

Identification of immune status among the different risk subgroups

First, we conducted PCA of the distribution of the 175 LUSC patients using total gene expression profiles and the expression profiles of the four immune-related lncRNAs of the model (Figure 4A, B). The 175 LUSC patients were better divided into two groups by our risk model. This result suggested that the immune status of LUSC patients with middle- or advanced-stage disease differed in the two risk groups.

Figure 4

A, B PCA based on the whole gene expression profiles (A) and the four immune-related lncRNA signatures (B). Green dots represent the low-risk group, and red dots represent the high-risk group. C Correlation of infiltrating immune cells with risk scores assessed by six methods. D-F Evaluation of the expression of immune checkpoint molecules between the high- and low-risk groups. NS. -P>0.05, ∗∗∗P < 0.001.

Subsequently, to explore the correlation between the risk score and the TME, we adopted several currently acknowledged methods to assess immune cell infiltration in the 175 samples. We discovered that the high-risk group was positively associated with the presence of CD8+ T cells, B cells, myeloid dendritic cells (DCs), macrophages, regulatory T cells (Tregs) and cancer-associated fibroblasts (CAFs), which remained consistent across several methods (Figure 4C).

Finally, since immunotherapy has become a hot spot in the treatment of NSCLC patients, we also evaluated the expression of ICI-related biomarkers in the high- and low-risk groups. The high-risk group highly expressed PD1 (P value <0.001, Figure 4D), CTLA4 (P value <0.001, Figure 4E) and PD-L1 (P value >0.05, Figure 4F); however, PD-L1 expression was not significantly different.

Exploring the correlation between chemotherapy and risk scores

Chemotherapy is an important clinical treatment for middle- to advanced-stage lung squamous carcinoma patients. We evaluated the correlation between risk scores and the sensitivity of the common chemotherapeutic agents used for lung squamous carcinoma patient treatment. The results are presented in Figure 5 and illustrate that a low-risk score was associated with lower half maximal inhibitory concentrations (IC50s) of cisplatin (P value <0.01), docetaxel (P value <0.001), vinorelbine (P value <0.001), and paclitaxel (P value <0.001), while it was associated with a higher IC50 of gemcitabine (P value <0.05). The IC50 of etoposide was lower in the low-risk group than in the high-risk group, but the difference was not statistically significant (P value >0.05). This finding suggested that our risk model may be of value in predicting the efficacy of clinical chemotherapy. Low-risk group is more likely to benefit from the treatment with cisplatin, docetaxel, vinorelbine and paclitaxel, while high-risk group is more likely to benefit from the treatment with gemcitabine.

Figure 5

Box plots of the IC50 values for several chemotherapeutic agents in the high- and low-risk groups.

Functional enrichment analysis of the risk score model

To identify the biological functions of the genes identified in our immune-related lncRNA risk score model, we used gene set enrichment analysis (GSEA). We displayed the most significant enrichment terms based on Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Hallmark, and WikiPathways gene sets from the high-risk group (Figure 6A-D). These terms mainly included the Kit receptor signaling pathway (Figure 6E), the KRAS signaling pathway (Figure 6F) and many immune-related pathways. In addition, the GSEA results also showed that the immune response (Figure 6G) and immune system processes (Figure 6H) were enriched in the high-risk group.

Figure 6

A-D The most significantly enriched terms based on GO (A), KEGG (B), hallmark (C) and WikiPathways (D) gene sets obtained by GSEA. E-H GSEA of gene sets for the Kit receptor signaling pathway (E), KRAS signaling pathway (F), immune response (G) and immune system processes (H). NES, normalized enrichment score; FDR, false discovery rate.

Validation of the immune-related lncRNA risk score model in a clinical cancer cohort

To further validate our prognostic immune-related lncRNA model, we examined the relative expression of the four identified lncRNAs in 14 stage IIB-IV LUSC patients from Fudan University Shanghai Cancer Center. We calculated the risk scores of the 14 patients and divided them into high-risk and low-risk groups. The DFS curve showed a statistically nonsignificant trend: the higher the risk score was, the higher the likelihood of recurrence (P value >0.05, Figure 7).

Figure 7

Kaplan–Meier disease-free survival (DFS) curve of our clinical cancer cohort (P = 0.329).

Discussion

Currently, the treatment of lung cancer has entered the era of individualized medicine. Most NSCLC patients with middle- to advanced-stage disease are treated with chemotherapy and/or radiation therapy, and several targeted and immunotherapy drugs are also available for some NSCLC patients [4]. Therefore, it is not sufficient to rely solely on the AJCC TNM staging system for the selection of treatments and the prediction of prognosis for patients. In the era of precision medicine, more and more robust predictive models are urgently needed to enable individualized risk stratification and treatment.

It has been affirmed that lncRNAs are critical regulators of cancer pathways and are closely associated with most cancer signatures, such as sustained proliferation, replicative immortality, evasion of growth inhibitors, induction of angiogenesis, resistance to cell death and development of metastasis [13, 16]. Recent studies have also indicated the potential of lncRNAs as new diagnostic and prognostic biomarkers. LncRNA PCA3 has been used as a noninvasive prediagnostic marker for prostate cancer in clinical practice [32]. In patients with metastatic renal cell carcinoma (RCC), OS was worse in those with high lncARSR expression, and lncARSR could promote resistance to sunitinib [33]. In patients with colon cancer, RAMS11 overexpression was also demonstrated to be related to worse prognosis [34].

Based on the above information, we hoped to build a predictive model for NSCLC patients using immune-related lncRNAs. Immune-related lncRNA risk score models for lung adenocarcinoma have been reported in many studies [35, 36, 37, 38, 39, 40, 41]. Our research is the first to construct an immune-related lncRNA risk score model in LUSC patients with middle- or advanced-stage disease. Based on the RNA expression profiles and the clinical data of LUSC patients publicly available from TCGA database, we constructed a model using four immune-related lncRNAs assessed by difference analysis, Pearson correlation analysis, and univariate and multivariate Cox regression analysis. The calculated risk scores were interpreted by the following formula: risk score = - 0.136673 × ExpressionAC020907.1 + 0.293362 × ExpressionAC027682.6 + 0.301214 × ExpressionAL034550.2 + 0.120110 × ExpressionLINC00944. Among the lncRNAs identified, in breast cancer, previously LINC00944 was confirmed to be in positive association with the presence of tumor-infiltrating T cells and proapoptotic markers, and low expression of LINC00944 showed correlation with poorer OS and relapse-free survival (RFS) [42]. However, an oncogenic role of LINC00944 in RCC has also been reported [43]. In a pancancer study of immune-related lncRNAs, LINC00944 was validated as a cancer-associated lncRNA whose expression was significantly perturbed in breast, kidney, lung and colorectal cancers and was influenced by the T-cell receptor (TCR) signaling pathway in 16 cancer categories [44]. It has been proven that the expression of AL034550.2 is remarkably higher in CD8+ T cells from patients with polymyositis (PM) than in patients with dermatomyositis (DM) [45]. In addition, AC020907.1 is considered to be one of the possible oncogenes of lung adenocarcinoma (LUAD) [46]. However, we could not retrieve any relevant reports for AC027682.6. In terms of survival analysis, both the OS and DFS of high-risk patients were worse than those of low-risk patients, suggesting a poorer prognosis. The subsequent independent prognostic analysis, PCA, and AUC data all showed that our model could serve as an independent predictor with good accuracy.

Platinum-based chemotherapy is an important treatment for middle- to advanced-stage LUSC patients. We evaluated the sensitivity to common chemotherapy drugs in high- and low-risk patients. Low-risk patients had higher sensitivity to cisplatin, docetaxel, vinorelbine and paclitaxel, while high-risk patients were more sensitive to gemcitabine. This finding implies that our model might be potentially helpful in the therapeutic decision-making process. Immunotherapy has also been proven to improve the outcome of NSCLC patients [47]. Our results suggested that high-risk patients with a high expression of PD1 and CTLA4 were probably more likely to benefit from immunotherapy, but this finding remains to be validated by further studies. The important role of the TME in the development of primary lung cancer has been recognized over the past decade, and the use of the TME indicators as predictive biomarkers has been extensively studied [48]. Therefore, we further assessed the relationship between the risk score and the immune microenvironment, revealing that the risk score was positively correlated with the presence of CD8+ T cells, B cells, myeloid DCs, macrophages, Tregs and CAFs. CAFs are considered as a factor of bad prognosis in NSCLC patients [49, 50, 51, 52, 53, 54]. Tregs have been demonstrated to be associated with worse RFS in NSCLC patients without lymph node metastases [55]. O'Callaghan DS et al. have also reported that the infiltration of tumor islet Foxp3 (+) T-cells indicates worse prognosis among NSCLC patients [56]. Recently, a large study revealed that LUAD patients with high expression of immune-infiltrating Treg-related genes had worse OS [57]. A high density of mature DCs in tertiary lymphoid structures suggests a better prognosis [58]. In contrast, the prognostic role of tumor-associated macrophages (TAMs) has been controversial [59, 60], and more studies are required to further our understanding of the differences in local TAM phenotypes [48].

The GSEA results suggested possible biological functions of our model, indicating that these lncRNAs may have critical and important effects in the progression of LUSC as well as in the diagnosis and treatment of LUSC patients. c-Kit is a tyrosine kinase receptor that is deregulated and then activates downstream signaling in diseases such as cancer [61]. Activated signaling cascades include the Ras/Raf/MEK/MAPK and PI3K/AKT/RPS6K pathways. KIT stimulation is as well known to activate the JAK/STAT and PLC/PKC signaling pathways [62]. c-Kit activation has been demonstrated to be associated with a variety of human cancers, including hematopoietic malignancies, malignant melanomas, gastrointestinal stromal tumors, and small cell lung cancers [61]. KRAS belongs to the rat sarcoma virus (Ras) gene family, and Ras genes are key players in human tumor pathogenesis [63]. KRAS is a commonly mutated gene in NSCLC patients, occurring in approximately 30% of LUAD [64] and 6% of LUSC patients [65]. Mutations in Ras genes induce sustained activation of downstream pathways, leading to tumor growth, proliferation and survival [66].

However, there are still some shortcomings of our study. For our clinical cancer cohort, we were unable to obtain OS data limited by the relatively short follow-up period. Additionally, the DFS curve was not statistically significant, which we speculated might be related to the insufficient sample size. Currently, we are unable to obtain additional fresh tissue samples from LUSC patients eligible for enrollment. Moreover, we failed to retrieve publicly available datasets containing the expression profiles of these four lncRNAs, the clinicopathological characteristics and the survival outcomes of LUSC patients with middle- or advanced-stage disease. In the future, we will collect more clinical samples, expand the sample size and extend the follow-up time to further validate and evaluate our risk model.

In conclusion, we first constructed a risk score model of immune-related lncRNAs in LUSC patients with middle- or advanced-stage disease. The results demonstrated that our four immune-related lncRNA model performed well in individualized risk stratification and prognosis prediction.

Declarations

Author contribution statement

Qianqian Xue & Yue Wang: Conceived and designed the experiment; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Qiang Zheng: Conceived and designed the experiment; Contributed reagents, materials, analysis tools or data.

Lijun Chen: Performed the experiments.

Yan Jin & Xuxia Shen: Contributed reagents, materials, analysis tools or data.

Yuan Li: Conceived and designed the experiment; Contributed reagents, materials, analysis tools or data.

Funding statement

Prof. Yuan Li was supported by National Natural Science Foundation of China [81972171].

Prof. Yuan Li was supported by CSCO-PILOT Cancer Research Foundation [Y-2019AZMS-0492].

Data availability statement

Data included in article/supplementary material/referenced in article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.