Introduction

Lung cancer is the leading cause of cancer-related deaths worldwide, accounting for nearly 1.5 million deaths annually. In the United States, lung cancer causes 158,000 fatalities, representing 28% of cancer deaths in men and 26% in women.1 Despite extensive research over the past fifty years, the 5-year survival rate for lung cancer remains dismal at 17.7%.2,3 This poor prognosis is largely due to the fact that most cases are diagnosed at an advanced stage, with 57% of lung cancer cases being diagnosed at stage IV, where the survival rate drops to less than 8%.4

Among the subtypes of lung cancer, lung squamous cell carcinoma (LUSC) represents approximately 40% of cases and is often associated with a poorer clinical outcome compared to lung adenocarcinoma.5 LUSC is more prevalent in smokers and is characterized by limited therapeutic options, especially due to the lack of specific biomarkers for its initiation and progression. Given its aggressive nature and poor response to standard therapies, the identification of novel biomarkers for LUSC is critical for improving prognosis and treatment strategies.

One promising area of investigation is immunogenic cell death (ICD), a form of regulated cell death that triggers an anti-tumor immune response by activating cytotoxic T cells and promoting inflammation.6,7 ICD is marked by the release of damage-associated molecular patterns (DAMPs) such as calreticulin, HMGB1, ATP, and annexin A1, which enhance the immune system’s recognition of tumors.8,9 Recent studies suggest that ICD plays a crucial role in cancer progression and response to immunotherapy, yet its impact on LUSC prognosis remains poorly understood.

The goal of this study is to explore the role of ICD-related genes in the survival of LUSC patients. By integrating transcriptomic data from 504 LUSC samples with clinical follow-up information, we identified distinct ICD-related gene expression patterns that correlate with significant differences in survival rates and treatment responses. This study aims to uncover potential therapeutic targets associated with ICD that could offer new avenues for LUSC treatment.

Materials and Methods Data Download and Processing

RNA sequencing data and associated clinical details, including treatment information, patient age, and tumor stage, for 530 cases of lung squamous cell carcinoma (LUSC) were extracted from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/, This study is exempt from approval in accordance with the national legislative guidelines, such as Article 32 (1) and (2) of the “Measures for the Ethical Review of Life Science and Medical Research Involving Human Subjects” issued by China on February 18, 2023).10 In addition, 34 genes related to immune cell dysfunction (ICDRGs) were identified in previous research and are listed in a Supplementary Table 1. The patients were categorized based on clinical characteristics such as tumor stage, grade, and age, which were reviewed statistically.

Survival data from the Gene Expression Omnibus (GEO) database were also utilized. Specifically, the GSE30219 dataset11 provided expression profiles and survival outcomes for 307 LUSC patients, while the GSE37745 dataset12 contributed data from an additional 196 patients. To identify differentially expressed ICDRGs within the TCGA dataset, we used the “limma” package, applying a P-value cutoff of 0.05. The expression correlation of these ICDRGs was subsequently analyzed.13

Importantly, the data from both RNA-seq (TCGA) and microarray (GEO) platforms were harmonized using appropriate normalization techniques to mitigate any batch effects, enabling integration without bias toward platform-specific differences. The goal of our analysis was to build a robust predictive model based on gene expression data, which remains consistent regardless of the sequencing technology used. This approach ensures the validity and generalizability of our model.

All clinical data used in the analysis were obtained from public repositories, so patient consent or ethical committee approval was not required.

Unsupervised Clustering Analysis Identified ICDRGs Expression Patterns

To dissect the immune dysfunction landscape in LUSC, Unsupervised clustering methodology was employed to identify distinct patterns of immune cell dysfunction (ICD) based on the expression of ICDRGs.14 We employed a differential expression-focused approach. Using the ConsensusClusterPlus package (v1.62.0), we performed unsupervised clustering of AML patients with parameters set to maxK = 6, clusterAlg = “hc” (hierarchical clustering), and distance = “Pearson”. This process determined the optimal number of immune subtypes (K) by tracking the cumulative distribution function and its relative change.This approach involved the analysis of 482 tumor samples using the established ICD model. The ConsensClusterPlus package was utilized to conduct this analysis,15 which was repeated 1000 times to ensure the stability of the clusters. The overall survival rate was computed utilizing the Kaplan-Meier survival curve.

Analysis of Immune Microenvironment Differences Between Subtypes

To explore immune-related differences between subtypes, the “GSVA” package was used to conduct GSVA score analysis on LUSC patients using differentially expressed ICDRGs, calculating the difference in GSVA scores between the two expression patterns of ICDRGs.16 CIBERSORT methodology, along with the R packages e1071 and preconditioningCore, as well as the LM22 signature from the CIBERSORT website, were employed to determine the fraction of 22 immune cell-infiltrating tumors.17 Gene set variation and enrichment data were analyzed to explore biological processes. The “ESTIMATE” package assessed tumor purity, immune scores, stromal scores, and ESTIMATE scores.18 Spearman correlation analysis confirmed the differential expression of immune checkpoint molecules. The Kruskal test investigated variations in immunological checkpoints, HLA-related genes, and immune subtype expression among different classifications.

Functional Enrichment Analysis of Genes

To investigate the biological functions associated with ICDRGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment studies were conducted using several R packages, such as “org.Hs.Eg.Db”, “ggplot2”, and “clusterProfiler”, to investigate the biological processes enriched by the genes in key modules.19 Any GO or KEGG terms that showed significant enrichment were examined further with a p-value threshold of 0.05.

Establishing a Prognosis Signature Related to ICDGs

To establish a robust prognostic model, univariate Cox regression analysis was conducted on all candidate ICDRGs to identify prognostic genes with significant predictive value (P<0.05) in the training set. Prior to multivariate Cox analysis to create a signature, the range of possible genes was further examined using the lowest absolute shrinkage and selection operator Cox regression model through the “Survival” R package.20 The TCGA samples (n = 541) were separated into a training group (n = 272) and a test group (n = 269) using a propensity score matching model. We utilized LASSO Cox analysis to identify genes most prognostically significant for overall survival (OS). To mitigate the risk of overfitting, we applied a 10-fold cross-validation using the glmnet package within the R statistical environment. Subsequently, we calculated a risk score for each patient based on their RNA expression levels, using the following formula:

(1)

where coefi represents the coefficients for each gene and Xi denotes the standardized counts for these genes. Utilizing clinical follow-up data from the lung cancer datasets GSE30219 and GSE37745, we computed risk scores for each sample within the TCGA training set and the GEO validation set. Samples were then categorized into high-risk (above the median risk score) and low-risk (at or below the median risk score) groups. The Kaplan-Meier curve method was employed to assess the association between these groups and the actual survival outcomes. The robustness of the risk model was further validated by examining risk scores from an external independent dataset.

Chemotherapy Response Analysis

To explore potential therapeutic implications, drug sensitivity analysis was performed using the Cancer Drug Sensitivity Genomics (GDSC) database.21 The pRRophetic package in R was employed to predict drug responses based on gene expression data, identifying differential drug sensitivity patterns between LUSC subtypes. This analysis revealed that high-risk patients were more resistant to certain chemotherapeutic agents while being more susceptible to immune checkpoint inhibitors. These findings provided critical insights into the potential for personalized treatment strategies based on ICDRG expression profiles.22

Quantitative Reverse Transcriptase PCR (qRT-PCR)

Predicting genes by risk models of protein and protein function interactions, we identified five core genes. These five genes may be the key markers of lung cancer prognosis, and further experimental studies are needed. A549 (Human Lung Squamous Cell) and BEAS-2B (non-small cell lung cancer cell line) were obtained (Mall North South Chuanglian Biotechnology Co., LTD). Control groups were created by combining A549 with BEAS-2B cells to compare the expression of AKR1B1, LOX, SERPINA1, SERPINA5 and GPC3 between normal and lung cancer cell lines. Total RNA was extracted using the Redzol kit (Beijing SBS Gene Technology Co., Ltd), and qRT-PCR was performed using the SYBR® Premix Ex Taq™ II Kit. The relative mRNA expression levels were calculated using the 2−ΔΔCt method, with GAPDH as the internal reference gene. The forward primer sequences were as follows: (1) AKR1B1: F-5′-GCGCAACCAATCAGAAGGC-3′, (2) GPC3: F-5′-AGGGCTAGACTTACAGATTGGC-3′, (3) LOX: F-5′-GTGGGCGAAGGTACAGCATA-3′, (4) SERPINA5: F-5′-AGCAATGCGGTCGTGAT-3′, (5) SERPINA1: F-5′-CAAGAAGGAGGAGGGGGTC-3′, (6) GAPDH: F-5′-TGAAGGTCGGAGTCAACGGATT-3′. The reverse primer sequences were as follows: (1) AKR1B1: R-5′-AACCACATTGCCCGACTCAT-3′, (2) GPC3: R-5′-ATGTAGCCAGGCAAAGCACT-3′, (3) LOX: R-5′-TGACAACTGTGCCATTCCCA-3′, (4) SERPINA5: R-5′-TCCGGTCCAGGAGGTAG-3′, (5) SERPINA1: R-5′-GAAGACGGCATTGTCCTGTG-3′, (6) GAPDH:R-5′-CCTGGAAGATGGTGATGGGATT-3′.

Cell Culture and Lentiviral Transduction

Human non-small cell lung cancer A549 cells (ATCC) were cultured in Dulbecco’s Modified Eagle Medium (DMEM; Gibco) supplemented with 10% fetal bovine serum (FBS; Gibco) and 1% penicillin/streptomycin (Sigma-Aldrich) at 37°C in a humidified 5% CO2 incubator. To establish stable knockdown cell lines, lentiviral particles expressing short hairpin RNA (shRNA) targeting AKR1B1, LOX, SERPINA1, SERPINA5, and GPC3 (sequences listed in Supplementary Table 2) were synthesized (GenePharma). Cells were transduced with lentivirus at a multiplicity of infection (MOI) of 10 in the presence of 8 μg/mL polybrene (Sigma-Aldrich). After 48 hours, puromycin (2 μg/mL; Invitrogen) was added to select stable clones. Non-targeting shRNA (scrambled sequence) was used as the negative control (NC).

CCK-8 Proliferation Assay

Cells were seeded in 96-well plates at a density of 5 × 10³ cells/well (100 μL medium/well) and incubated for 0, 24, 48, 72, 96, or 120 hours. At each time point, 10 μL of CCK-8 reagent (Dojindo Molecular Technologies) was added to each well, followed by 2 hours of incubation at 37°C. Absorbance was measured at 450 nm (reference wavelength: 630 nm) using a microplate reader (BioTek Synergy H1). Three independent biological replicates were performed for each group.

Data Analysis

Raw absorbance values were normalized to the 0-hour baseline for each group. Relative proliferation rates were calculated using the formula:

Statistical significance was determined by two-way ANOVA with Tukey’s post hoc test using GraphPad Prism 9.0. Data are presented as mean ± standard deviation (SD), with *p < 0.05 and **p < 0.01 considered statistically significant.

Statistical Analysis

Survival analysis was conducted using the R survival package, and the Log Rank test was utilized to assess the survival rates of each group. The Kruskal–Wallis test was employed to compare the differences between two or more sets of data, while the Wilcoxon test was used to compare the differences between the two groups. Kaplan-Meier technique was applied to generate survival curves for each subgroup within the dataset. The Spearman correlation analysis was used to determine the correlation coefficient. To account for the multiple hypothesis testing and control the false-positive rate, we applied the Benjamini-Hochberg correction for multiple comparisons where appropriate. Statistical significance was determined at P < 0.05 for all analyses. All statistical calculations were performed using R versions 4.1.0 and 4.0.0.23

Results Expression Levels of ICDRGs in LUSC

Transcriptomes and clinical data for 504 normal and LUSC tissue samples were obtained from the TCGA database. This study encompasses 452 LUSC cancer patients with both clinical information and gene expression profiles. The expression of 28 ARGs was assessed using the Wilcoxon test, resulting in the identification of 14 ICDRGs with high tumor expression (|log2(FC)|, P value < 0.05). The expression correlation of these ICDRGs was examined based on their tumor expression levels (Figure 1a). Utilizing an unsupervised clustering method, two distinct regulatory patterns were identified based on the expression levels of survival-related ICDRGs in the TCGA database (Figure 1b). A total of 258 cases were classified into ICDRGs group 1, while 272 cases were classified into ICDRGs group 2 (Figure 1c and d). Analyzing immune cell dysfunction-related genes, two types of immune cell dysfunction-related scores were assessed using the GSVA method. The immune cell dysfunction-related scores for group 1 were found to be higher than those for group 2 (Figure 1E). To evaluate the survival difference between the two groups, patients from different groups (TCGA) were analyzed across various ranges, and their complete survival information was examined (Figure 1F), Figure 1G shows the grouping of samples and the heatmap of the expression level of ICDRGs.

Figure 1 Consensus clustering of ICDRGs in LUSC. (a) Box plot showing the expression of ICDRGs in tumor and normal LUSC samples, indicating significant differential expression between the two groups. (b) Heatmap displaying the expression correlation of ICDRGs, providing insight into the relationships between these genes across the LUSC samples. (c) Consensus matrices were generated to assess the stability of clustering. The figure shows how the stability increases with the number of clusters, confirming the robustness of the clustering results. (d) The consensus cumulative distribution function (CDF) curve was used to explore the optimal number of clusters, and the sharp increase in the CDF indicates the appropriate cluster stability. (e) GSVA (Gene Set Variation Analysis) scores for ICDRGs across different patient groups, highlighting how distinct subtypes exhibit different levels of immune-related gene activation. (f) Prognostic analysis comparing survival outcomes for the two ICDRG-based subgroups (Group 1 and Group 2), stratified by LUSC tumor stages. Significant differences in survival between the groups suggest the clinical relevance of ICDRGs in patient prognosis. (g) Unsupervised clustering heatmap of all ICDRGs across the TCGA LUSC cohort, with annotations for ICDRG groups, tumor stage, age, grade, and survival status. This heatmap visually supports the hypothesis that ICDRG expression correlates with patient outcomes and clinical features. *Statistical significance indicated by *p <0.05, **p <0.01, ***p <0.001 and ****p <0.0001.

Immune Cell Infiltration Evaluation

Additionally, the 22 immune cell subtypes in the TCGA samples were examined within two subtype clusters using the CIBERSORT method. The analysis revealed that memory B cells, T cells CD4 memory resting, follicular helper cells, T cells gamma delta, macrophages M0, macrophages M2, resting mast cells, and eosinophils were infiltrated by the two subtype clusters, contributing significantly to the overall immune cell infiltration (Figure 2a). Furthermore, group 1, which had more diverse immune cell populations, exhibited higher immune, stromal, and estimated scores compared to other subtype groups, such as group 2, which contained more pure tumor cells (Figure 2b). Cancer has demonstrated a strong response to immune checkpoint inhibition, which blocks inhibitory signals of T cell activation. Figure 2c displays the expression of HLA family genes in the two clusters, with group 1 showing higher HLA gene expression. Immunological cells express immunological checkpoints, which regulate the level of immune activation and are vital in preventing the onset of autoimmunity. Group 1 exhibited increased expression of PDCD1 and LAG3.

Figure 2 Clinical significance and immune landscape of ICDRGs groups in the TCGA cohort. (a) Immune cell infiltration analysis between ICDRG groups. (b) Scores of stromal, immune, estimate, and tumor purity were compared across ICDRG groups. (c) Gene expression levels of HLA and immune checkpoint markers were analyzed across ICDRG groups. *Statistical significance indicated by *p <0.05, **p <0.01, ***p <0.001 and ****p <0.0001.

Enrichment Analysis

A total of 730 differentially expressed genes (DEGs) were identified between group 1 and group 2 using the limma method under the filtering conditions of FDR q value 0.01 and absolute value of logFC > 0.5 (Figure 3a, Supplementary Table 1). Pathway enrichment analysis ultimately demonstrated that the Estrogen Signaling Pathway and Phagosome were the key enrichments for the differentially expressed genes in group 1 and group 2 (Figure 3b).Gene enrichment analysis was conducted to investigate the biological pathways and processes associated with these groups. Figure 3c illustrates that the majority of biological processes related to immune response activation and molecule synthesis were significantly enriched. Furthermore, the cellular components such as the exterior side of the plasma membrane and the immunoglobulin complex were found to be significantly enriched (Figure 3d). These findings relate to the body’s active immune metabolic processes, reflecting the body’s immune response as cancer progresses. Analysis of Figure 4e revealed that the genes in molecular function were primarily enriched in antigen binding and immune receptor activation.

Figure 3 Differential gene expression and functional analysis of ICDRGs groups. (a) Volcano plot of differentially expressed genes (DEGs) between ICDRG groups. (b) KEGG pathway analysis of DEGs. (c) Biological process analysis of DEGs. (d) Cellular component analysis of DEGs. (e) Molecular function analysis of DEGs.

Figure 4 Construction of the risk model. (a) Partial likelihood deviance for the LASSO regression model during tenfold cross-validation, with vertical dotted lines indicating the optimal values determined by the minimum and 1-SE criteria. (b) LASSO coefficient profiles for 12 selected genes based on tenfold cross-validation, with vertical dotted lines representing the optimal values. (c–e) Kaplan-Meier survival curves for test and training sets in TCGA, GSE30219, and GSE37745 datasets, respectively, with a Log rank test (P < 0.001). (f–h) Survival status of patients in high-risk and low-risk subgroups based on median levels of ICDRGs in TCGA, GSE30219, and GSE37745 datasets, respectively. (i–l) Time-dependent ROC curves for training and test sets in TCGA, GSE30219, and GSE37745 datasets, respectively.

Construction and Validation of the Prognostic Model

The LASSO approach was utilized to construct risk models, which ultimately led to the identification of 13 prognostically significant genes (Figure 4a–c). Subsequently, risk score-based models using these genes were developed based on the training dataset (TCGA, n = 530), test set 1 (GSE30219, n = 307), and test set 2 (GSE37745, n = 196) obtained from the TCGA LUSC and GEO datasets, respectively. The survival study revealed that higher risk scores in both the training and test sets were correlated with poorer survival (P<0.0001) (Figure 4d, e, g, h, j and k). The sensitivity of the prognostic model was assessed using time-dependent ROC curves. The results of the AUCs indicated that the 1-, 3-, and 5-year AUCs for the training set were 0.707, 0.655, and 0.792, respectively (Figure 4f), while the corresponding AUCs for test set 1 were 0.638, 0.66, and 0.684 (Figure 4i) and for test set 2 were 0.609, 0.642, and 0.648, respectively (Figure 4l).

Independent Prognostic Ability of Risk Model

Further analysis explored the relationship between the risk model and the clinicopathological characteristics of TCGA-LUSC, utilizing univariate regression analysis and multivariate Cox regression analysis to evaluate the independent predictive power of the developed risk model. Clinically relevant variables, including Grade, Age, TNM Stage, Stage, and Risk Model, were selected for analysis. The results indicate that the prognosis of LUSC patients may be independently predicted by Grade, Age, and Risk Model (Figure 5a and b). A nomogram was constructed to serve as a therapeutically useful quantitative tool for estimating the mortality of specific BC patients, integrating the independent prognostic indicators (Figure 6a). Patients are assigned a total score based on the sum of the results from each prognostic parameter, with higher scores corresponding to worse prognosis. The performance of the nomogram is comparable to that of the ideal model, as evidenced by the c-exponential curve of various variables across time in the TCGA cohort (Figure 6b–d). The risk model demonstrates greater accuracy in predicting patient survival than age, as determined by receiver operating curve analysis (Figure 6e). The nomogram outperforms other single factors, as indicated by the c-exponential curve of various variables across time in the TCGA cohort (Figure 6f). Additionally, the DCA curve reveals that the nomogram’s net benefit curve in terms of age is stable and reliable compared to other clinical parameters (Figure 6g).

Figure 5 Correlation between the risk model and clinicopathological features of TCGA–LUSC samples. (a) Univariate analysis including the risk model and clinical factors. (b) Multivariate analysis including the risk model and clinical factors. ***p < 0.001.

Figure 6 Correlation between groups and clinicopathological features of TCGA–LUSC samples. (a) Nomogram for predicting the 1-, 3-, and 5-year overall survival (OS) of LUSC patients in the TCGA cohort. (b–d) Calibration plots for predicting the 1-, 3-, and 5-year OS for TCGA–LUSC samples. (e) Time-dependent ROC curve for the risk model and clinical factors. (f) Time-dependent C-index plot for the nomogram and clinical factors. (g) Decision curve analysis to assess the clinical utility of the nomogram and clinical factors for 1-, 3-, and 5-year risk assessment.

Risk Model Can Predict How Chemotherapy Work

Eight chemotherapeutic medications’ IC50 discrepancies were investigated using the “pRRophetic” package to predict their sensitivity to drug therapy. The drug sensitivity data for Sorafenib, Gefitinib, Bleomycin, Bosutinib, Etoposide, Lenalidomide, Camptothecin, and Methotrexate in the LUSC risk model were presented in Figures 7a–h, respectively. Statistical analyses revealed that the IC50 values for these chemotherapeutic agents were higher in high-risk individuals.

Figure 7 Drug sensitivity prediction and risk groups. Violin plots showing drug sensitivity prediction scores for (a) Sorafenib, (b) Gefitinib, (c) Bleomycin, (d) Bosutinib, (e) Etoposide, (f) Lenalidomide, (g) Camptothecin, and (h) Methotrexate across risk groups. *Statistical significance: *p < 0.05; **p < 0.01; ***p < 0.001.

Gene Expression Level Verification via Quantitative Reverse Transcription PCR (qRT‑PCR)

Through protein-protein interaction (PPI) analysis, we identified a core network of genes associated with the risk model, which included AKR1B1, LOX, SERPINA1, SERPINA5, and GPC3. These genes were selected based on their strong prognostic relevance to LUSC patient survival, as determined by LASSO regression analysis. To validate the expression levels of these genes, qRT-PCR was performed on the LUSC tumor cell line A549 (human lung squamous cell carcinoma) and the adjacent normal cell line BEAS-2B. As shown in Figure 8, all five genes were expressed at significantly higher levels in the A549 cell line compared to the BEAS-2B cell line, indicating their potential role in LUSC progression and supporting their inclusion in the risk model.

Figure 8 Gene expression in LUSC cell lines and adjacent normal cell lines. (A) mRNA expression of AKR1B1 in LUSC cell lines and adjacent normal cell lines. (B) mRNA expression of GPC3. (C) mRNA expression of LOX. (D) mRNA expression of SERPINA5. (E) mRNA expression of SERPINA5 (duplicate panel). (F) Core network diagram for risk model genes. *Statistical significance: *p < 0.001, ***p < 0.001.

Targeted Gene Knockdown Suppresses A549 Proliferation

The CCK-8 assay demonstrated that knockdown of AKR1B1, LOX, SERPINA1, SERPINA5, and GPC3 significantly impaired the proliferative capacity of A549 cells in a time- and gene-dependent manner (Figure 9). Compared to the negative control (NC) group, which exhibited a sigmoidal growth curve peaking at 120 hours (OD450 = 1.55 ± 0.02), all knockdown groups showed reduced proliferation rates. Notably, AKR1B1 knockdown (AKR1B1-KD) exerted the strongest suppression, with a 19.4% reduction in OD450 at 120 hours (1.25 ± 0.02, **p < 0.01), followed by LOX knockdown (LOX-KD), which displayed progressive inhibition over time (8% at 48 hours vs 22.6% at 120 hours; 1.20 ± 0.03, **p < 0.01). In contrast, SERPINA1 knockdown (SERPINA1-KD) induced rapid suppression, achieving statistical significance as early as 24 hours (OD450 = 0.30 ± 0.01 vs NC 0.35 ± 0.02, *p < 0.05), while SERPINA5 knockdown (SERPINA5-KD) exhibited a stepwise inhibition pattern, reaching 9.7% suppression at 96 hours (1.40 ± 0.02, *p < 0.05). Interestingly, GPC3 knockdown (GPC3-KD) showed minimal early effects but partial recovery at late stages (120-hour OD450 = 1.50 ± 0.02, ns, no significance), suggesting compensatory mechanisms may attenuate its long-term inhibitory role. Cross-comparative analysis further revealed divergent dynamics among genes: LOX-KD and SERPINA5-KD exhibited the largest discrepancy in suppression efficacy at 48 hours (8% vs 15%, **p < 0.01), whereas GPC3-KD uniquely approached NC levels by 120 hours, implicating its context-dependent regulatory function in proliferation (Figure 9). These findings collectively highlight AKR1B1 and LOX as potent suppressors of A549 growth, with potential therapeutic relevance in LUSC.

Figure 9 Inhibition rates and OD values at 120 hours for risk model genes. (a) AKR1B1 inhibition rates and OD values at 120 hours. (b) GPC3 inhibition rates and OD values. (c) LOX inhibition rates and OD values. (d) SERPINA5 inhibition rates and OD values.(e) SERPINA5 (duplicate panel) inhibition rates and OD values. *Statistical significance: *p < 0.05; **p < 0.01; ns, no significance.

Discussion

Immunogenic cell death (ICD) is a form of regulated cell death that triggers an adaptive immune response in an immunocompetent environment. It is characterized as a unique, albeit poorly defined, phenomenon that can be induced by various factors, including chemotherapeutic agents, oncolytic viruses, physicochemical treatments, photodynamic therapy, and radiation. Upon cellular death or stress, dying cells release chemicals that the immune system can utilize as adjuvants or warning signals. Collectively, these signals are referred to as damage-associated molecular patterns (DAMPs) and are recognized by innate pattern recognition receptors such as Toll-like receptors and NOD-like receptors, activating tumor-specific immune responses. This mechanism not only contributes to the direct killing of cancer cells but also enhances the long-term efficacy of anticancer drugs by establishing immunological memory.

Previous literature has identified 34 ICDRGs. In this study, we focused on examining the expression differences of these genes between LUSC tumor tissues and adjacent normal tissues, revealing 14 ICDRGs with significant expression variations.24 These ICDRGs were further analyzed through unsupervised consistency clustering, which identified two groups significantly associated with the prognosis of LUSC patients. By performing biological function enrichment analysis on the differentially expressed genes of these two groups, we developed a prognostic survival model and identified key ICDRGs that could serve as potential biomarkers for patient prognosis and treatment response.25

Importantly, our findings demonstrate that these ICD-related subtypes not only correlate with significant survival differences but also with variations in chemotherapy sensitivity based on patient risk profiles. This highlights the potential clinical value of ICDRGs in stratifying LUSC patients for more personalized treatment approaches.26 Given the observed differences in immune responses and prognosis, these ICD-related markers could guide therapeutic decisions, particularly in the context of immunotherapy.27

Recent advancements in cancer immunotherapy have focused on understanding how tumor cells evade immune detection and resist standard treatments.28 ICD represents a promising strategy to overcome this resistance by triggering T-cell adaptive immune responses and establishing long-lasting immunological memory.29,30 Many cancer therapies have the potential to induce ICD, offering a multi-dimensional framework for categorizing, verifying targets, and assessing drug sensitivity.31 As immunotherapy continues to evolve, our findings suggest that targeting ICD-related pathways could enhance the effectiveness of current therapies, particularly for patients with LUSC.

The clinical relevance of ICD in cancer treatment necessitates further validation through clinical data on LUSC patient responses. Future studies should examine the specific effects of ICD on treatment outcomes, focusing on how these pathways influence patient responses to targeted therapies. As we continue to refine treatment regimens, understanding the role of ICD in LUSC could lead to improved patient stratification, more effective treatment plans, and ultimately better clinical outcomes.

Conclusion

In this study, we summarized 34 immunogenic cell death (ICD)-related genes from the existing literature and explored their biological roles in lung squamous cell carcinoma (LUSC). We first analyzed the expression differences of these ICD genes between tumor and adjacent tissues in LUSC patients, identifying 14 ICD-related genes (ICDRGs) with significant expression changes. Using these ICDRGs, we applied unsupervised consistency clustering to the LUSC expression matrix, revealing two distinct clusters that were significantly associated with patient prognosis. Further biological function enrichment analysis of the differentially expressed genes between these clusters led to the construction of a prognostic survival model, and key ICDRGs were identified as potential prognostic biomarkers.

Our findings suggest that ICD can activate the adaptive immune response in the context of normal immune function, offering promising insights for cancer treatment. The two immunogenic cell death-related subtypes we identified showed significant differences in survival outcomes, and our risk-based model demonstrated varying chemotherapy sensitivities between different risk groups. These results underline the potential for applying ICDRG-based models to guide clinical immunotherapy strategies in LUSC patients.

However, to further substantiate these findings, additional experimental validation and clinical data are needed to confirm the clinical utility of the identified ICDRGs and to refine the prognostic model for broader application in clinical settings.

Abbreviations

ICD, immunogenic cell death; TME, tumor microenvironment; LUSC, lung squamous cell carcinoma; ICDRGs, immune cell dysfunction-related genes; GEO, Gene Expression Omnibus; TCGA, The cancer genome atlas; GO, Gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; GDSC, The Cancer Drug Sensitivity Genomics; qRT-PCR, Quantitative Reverse transcriptase PCR.

Data Sharing Statement

In accordance with the national legislative guidelines, Article 32 (1) and (2) of the Measures for the Ethical Review of Life Science and Medical Research Involving Human Subjects are exempted from approval. Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/. The names of the repository/repositories and accession number(s) can be found in the article Supplementary Material.

Ethics Approval and Consent to Participate

Our experiments do not involve animal or human testing and do not require ethics committee approval.

Consent for Publication

All authors approved the publication of the article.

Acknowledgments

Yuhan Wang, Shuang Wang are co-first authors for this study. Erli Zhang is the senior author for this study. All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Disclosure

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Kratzer TB, Bandi P, Freedman ND, et al. Lung cancer statistics, 2023. Cancer. 2024;130(8):1330–1348. doi:10.1002/cncr.35128

2. Oliver AL. Lung cancer: epidemiology and screening. Surg Clin North Am. 2022;102(3):335–344. doi:10.1016/j.suc.2021.12.001

3. Abu Rous F, Singhi EK, Sridhar A, Faisal MS, Desai A. Lung cancer treatment advances in 2022. Cancer Invest. 2023;41(1):12–24. doi:10.1080/07357907.2022.2119479

4. Pei Q, Luo Y, Chen Y, Li J, Xie D, Ye T. Artificial intelligence in clinical applications for lung cancer: diagnosis, treatment and prognosis. Clin Chem Lab Med. 2022;60(12):1974–1983. doi:10.1515/cclm-2022-0291

5. Li Y, Wu X, Yang P, Jiang G, Luo Y. Machine learning for lung cancer diagnosis, treatment, and prognosis. Genomics Proteomics Bioinf. 2022;20(5):850–866. doi:10.1016/j.gpb.2022.11.003

6. Chen JW, Dhahbi J. Lung adenocarcinoma and lung squamous cell carcinoma cancer classification, biomarker identification, and gene expression analysis using overlapping feature selection methods. Sci Rep. 2021;11(1):13323. doi:10.1038/s41598-021-92725-8

7. Kroemer G, Galassi C, Zitvogel L, Galluzzi L. Immunogenic cell stress and death. Nat Immunol. 2022;23(4):487–500. doi:10.1038/s41590-022-01132-2

8. Ando M, Ito M, Srirat T, Kondo T, Yoshimura A. Memory T cell, exhaustion, and tumor immunity. Immunol Med. 2020;43(1):1–9. doi:10.1080/25785826.2019.1698261

9. Zindel J, Kubes P. DAMPs, PAMPs, and LAMPs in immunity and sterile inflammation. Annu Rev Pathol. 2020;15:493–518. doi:10.1146/annurev-pathmechdis-012419-032847

10. Tomczak K, Czerwińska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol. 2015;19(1A):A68–77.

11. Rousseaux S, Debernardi A, Jacquiau B, et al. Ectopic activation of germline and placental genes identifies aggressive metastasis-prone lung cancers. Sci Trans Med. 2013;5(186):186ra66. doi:10.1126/scitranslmed.3005723

12. Botling J, Edlund K, Lohr M, et al. Biomarker discovery in non-small cell lung cancer: integrating gene expression profiling, meta-analysis, and tissue microarray validation. Clin Cancer Res. 2013;19(1):194–204. doi:10.1158/1078-0432.CCR-12-1139

13. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. doi:10.1093/nar/gkv007

14. Ding R, Wang Y, Fan J, et al. Identification of immunosuppressive signature subtypes and prognostic risk signatures in triple-negative breast cancer. Front Oncol. 2023;13:1108472. doi:10.3389/fonc.2023.1108472

15. Ding R, Liu Q, Yu J, et al. Identification of breast cancer subtypes by integrating genomic analysis with the immune microenvironment. ACS Omega. 2023;8(13):12217–12231. doi:10.1021/acsomega.2c08227

16. Tian Z, Yang Z, Jin M, et al. Identification of cytokine-predominant immunosuppressive class and prognostic risk signatures in glioma. J Cancer Res Clin Oncol. 2023;149:13185–13200. doi:10.1007/s00432-023-05173-4

17. Hu I, Zeng Y, Ge N, et al. Identification of the shared gene signatures and biological mechanism in type 2 diabetes and pancreatic cancer. Front Endocrinol. 2022;13:847760. doi:10.3389/fendo.2022.847760

18. Ping S, Wang S, Zhao Y, et al. Identification and validation of a ferroptosis-related gene signature for predicting survival in skin cutaneous melanoma. Cancer Med. 2022;11(18):3529–3541. doi:10.1002/cam4.4706

19. Li M, Ding R, Yang X, Ran D. Study on biomarkers related to the treatment of post-stroke depression and alternative medical treatment methods. Neuropsychiatr Dis Treat. 2022;18:1861–1873. doi:10.2147/NDT.S370848

20. Liu TT, Li R, Huo C, et al. Identification of CDK2-related immune forecast model and ceRNA in lung adenocarcinoma, a pan-cancer analysis. Front Cell Dev Biol. 2021;9:682002. doi:10.3389/fcell.2021.682002

21. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 2013;41(Database issue):D955–61. doi:10.1093/nar/gks1111

22. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One. 2014;9(9):e107468. doi:10.1371/journal.pone.0107468

23. Chan BKC. Data analysis using R programming. Adv Exp Med Biol. 2018;1082:47–122.

24. Fucikova J, Kepp O, Kasikova L, et al. Detection of immunogenic cell death and its relevance for cancer therapy. Cell Death Dis. 2020;11(11):1013. doi:10.1038/s41419-020-03221-2

25. Fucikova J, Moserova I, Urbanova L, et al. Prognostic and predictive value of DAMPs and DAMP-associated processes in cancer. Front Immunol. 2015;6:402. doi:10.3389/fimmu.2015.00402

26. Turhon M, Maimaiti A, Gheyret D, et al. An immunogenic cell death-related regulators classification patterns and immune microenvironment infiltration characterization in intracranial aneurysm based on machine learning. Front Immunol. 2022;13:1001320. doi:10.3389/fimmu.2022.1001320

27. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell mol Immunol. 2020;17(8):807–821. doi:10.1038/s41423-020-0488-6

28. Krysko DV, Garg AD, Kaczmarek A, Krysko O, Agostinis P, Vandenabeele P. Immunogenic cell death and DAMPs in cancer therapy. Nat Rev Cancer. 2012;12(12):860–875. doi:10.1038/nrc3380

29. Duan X, Chan C, Lin W. Nanoparticle-mediated immunogenic cell death enables and potentiates cancer immunotherapy. Angew Chem Int Ed Engl. 2019;58:670–680. doi:10.1002/anie.201804882

30. Ruan H, Leibowitz BJ, Zhang L, Yu J. Immunogenic cell death in colon cancer prevention and therapy. mol Carcinog. 2020;59(7):783–793. doi:10.1002/mc.23183

31. Li Z, Lai X, Fu S, et al. Immunogenic cell death activates the tumor immune microenvironment to boost the immunotherapy efficiency. Adv Sci. 2022;9(22):e2201734. doi:10.1002/advs.202201734