Integrated transcriptomic analysis and machine learning identify immunogenic cell death genes as prognostic markers and therapeutic targets in non-small cell lung cancer

Introduction

Non-small cell lung cancer (NSCLC), accounting for approximately 85% of all lung cancer cases, represents the most prevalent histological subtype. Based on tissue characteristics, NSCLC is classified into three major variants: adenocarcinoma, squamous cell carcinoma, and large cell carcinoma, each exhibiting distinct biological behaviors, prognostic profiles, and therapeutic responses. The staging system endorsed by the Union for International Cancer Control (UICC) and the American Joint Committee on Cancer (AJCC) categorizes NSCLC progression through tumor-node-metastasis (TNM) classification, which evaluates tumor size (T), lymph node metastasis (N), and distant metastasis (M). While early-stage NSCLC (e.g., stage I–II) may achieve curative outcomes through surgical resection, advanced stages (e.g., stage III–IV) often require multimodal therapies combining radiotherapy, chemotherapy, and targeted agents. Despite advancements in targeted therapies [e.g., epidermal growth factor receptor (EGFR), tyrosine kinase inhibitors (TKIs)] and immune checkpoint inhibitors (ICIs), clinical benefits remain suboptimal. The 5-year survival rate for advanced NSCLC patients persists below 20%, with particularly dismal prognoses observed in specific subtypes such as KRAS-mutant tumors. Recent insights into NSCLC molecular heterogeneity have driven the evolution of personalized treatment strategies to improve survival outcomes and quality of life (1,2).

Immunogenic cell death (ICD) refers to a programmed cell death modality capable of eliciting antitumor immune responses. Distinct from conventional apoptosis, ICD involves the release of specific immune-activating molecules termed damage-associated molecular patterns (DAMPs), including high-mobility group box 1 (HMGB1), calreticulin, and adenosine triphosphate (ATP), which collectively enhance antitumor immunity (3). Beyond directly eliminating cancer cells, ICD remodels the tumor microenvironment (TME) to promote tumor-specific immune activation and augment systemic immune surveillance (4). The molecular mechanisms of ICD involve endoplasmic reticulum stress (ERS) and apoptosis-related signaling cascades. Chemotherapeutic agents such as doxorubicin and platinum-based drugs induce ICD by activating ERS pathways (5). Small molecules like (-)-guaiol have been shown to potentiate ICD through dual modulation of autophagy and apoptosis (6). Notably, antidepressants (e.g., sertraline, indatraline) enhance ICD by inhibiting cholesterol binding to lysosomal transport proteins Niemann-Pick type C1/C2 intracellular cholesterol transporter 1/2 (NPC1/NPC2) while upregulating their expression, thereby disrupting cholesterol trafficking and serving as novel immunostimulants (7). In NSCLC, cisplatin-induced ICD promotes surface exposure of calreticulin and HMGB1 release, facilitating dendritic cell (DC) activation and T cell infiltration (8). Targeted therapies such as high-dose crizotinib synergize with ICIs by triggering ICD (9). Inhibition of glucose-6-phosphate dehydrogenase (G6PD) reduces nicotinamide adenine dinucleotide phosphate (NADPH) synthesis, impairing oxidative stress tolerance and inducing ICD-associated HMGB1 release and calreticulin membrane translocation, suggesting G6PD targeting as a promising strategy to improve ICI responsiveness (10). These discoveries collectively advance ICD-based therapeutic paradigms for NSCLC.

Machine learning (ML), as an emerging computing technology with the ability to process complex datasets and discover potential patterns, has been widely applied in the medical field, especially in oncology. By analyzing large-scale clinical and genomic data, ML enables the identification of critical biomarkers associated with disease prognosis, thereby advancing personalized therapeutic strategies (11). For instance, ML-based models have successfully integrated genome-wide association studies (GWAS) and transcriptomic data to pinpoint genetic risk factors in NSCLC (12). Furthermore, ML applications in radiomics, proteomics, and metabolomics have enhanced the precision of early diagnosis and prognostic evaluation in NSCLC. A representative example includes ML-driven analysis of serum proteomic profiles, which achieves high-accuracy classification of early-stage NSCLC, significantly improving diagnostic reliability (13). These advancements underscore the transformative potential of integrating ML with conventional biomarker research to redefine prognostic paradigms in NSCLC. However, there are still limitations in existing prognostic models for predicting individual patients’ response to immunotherapy and long-term survival. A robust prognostic model that specifically integrates the characteristics of ICD—a key determinant of anti-tumor immunity—has not been fully established and validated. The development of such models is expected to meet the unmet needs of biomarkers reflecting tumor immunogenicity potential and guiding personalized immunotherapy.

Therefore, the primary objective of this study was to develop and validate a novel ICD-based prognostic model for NSCLC. We first systematically identified ICD-related genes (ICDRGs) at single-cell and bulk transcriptomic levels. Using an integrative ML framework, we then developed a robust ICD-related risk score (ICDRS). This model was rigorously validated in independent cohorts to assess its predictive performance for patient survival. Furthermore, we characterized the association between the ICDRS and features of the tumor immune microenvironment, as well as therapeutic sensitivity. This work aims to provide a clinically relevant tool for prognosis stratification and to offer insights for personalizing ICD-targeting immunotherapy strategies in NSCLC. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2682/rc).

Methods

Data collection

This retrospective cohort study utilized publicly available genomic data from two primary sources.

The data for model development were obtained from The Cancer Genome Atlas (TCGA) database (download date: April 2025). This cohort included gene expression profiles and matched survival data for 598 NSCLC samples.

For external validation, data were retrieved from the Gene Expression Omnibus (GEO) database (download date: April 2025). Two independent NSCLC cohorts, GSE11969 and GSE68465, were used, collectively comprising tissue-based RNA sequencing (RNA-seq) data and survival information from 441 NSCLC samples. For further validation on an independent RNA-seq platform, the GSE81089 dataset was obtained from GEO (download date: February 2026). This cohort includes 199 NSCLC samples profiled by Illumina HiSeq 2000, providing a suitable platform to assess the cross-platform generalizability of the ICDRS. Additionally, the single-cell RNA sequencing (scRNA-seq) dataset GSE189357 (from 3 NSCLC samples) was accessed for preliminary exploration of ICD heterogeneity and was not used for prognostic model building.

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

To identify ICDRGs, we curated 34 ICD-associated genes from prior publications (14-16).

scRNA-seq

Data processing

scRNA-seq data were processed using the Seurat package (v4.3.3) in R. Quality control involved filtering low-quality cells, retaining those with 200–7,000 detected genes and mitochondrial gene content <10%. Erythrocyte contamination was assessed via hemoglobin genes (e.g., HBA1, HBB), excluding cells with >1% hemoglobin gene expression. Data normalization employed the LogNormalize method (scale factor =10,000), followed by identification of highly variable genes (HVGs) using the FindVariableFeatures function (top 2,000 HVGs retained). Batch effects across three samples were corrected using the Harmony package. Cell clustering was performed via FindClusters (resolution =0.8) and FindNeighbors (dims =1:30) functions, with dimensionality reduction visualized via t-distributed stochastic neighbor embedding (t-SNE). Finally, cell populations were annotated using canonical marker genes (e.g., CD3E for T cells, CD19 for B cells, CD68 for macrophages).

ICD scoring

To investigate the role of ICD in the TME, we quantified cellular ICD activity using single-sample gene set enrichment analysis (ssGSEA) based on a literature-curated core gene set (e.g., HMGB1, CALR, STING), with median expression thresholds stratifying high/low-ICD subgroups. Differential expression analysis between groups was performed via the FindAllMarkers function (log₂ fold-change threshold =0.35, minimum detection rate =35%). Pathway enrichment analysis was conducted using the clusterProfiler package, complemented by gene set enrichment analysis (GSEA) against Hallmark gene sets to identify significantly enriched pathways. Results were visualized through dot plots and enrichment plots to highlight ICD-associated pathways. Cellular ICD score distributions across cell types were illustrated via Uniform Manifold Approximation and Projection (UMAP) projections, with violin plots comparing scores among annotated cell populations. Statistical significance was defined as P<0.05.

Cell-cell communication analysis

Cell-cell communication networks were reconstructed using the CellChat platform (v1.6.1), focusing on tumor-specific subpopulations and their crosstalk with stromal/immune cells. The workflow comprised three phases: (I) quantitative comparison of interaction strengths via the netVisual_diffInteraction module; (II) extraction of key communication modules using the identifyCommunicationPatterns algorithm; (III) ligand-receptor database-integrated pathway enrichment. Significant interactions were filtered via 1,000 permutation tests [false discovery rate (FDR) <0.05], with communication entropy matrices quantifying network heterogeneity.

Weighted gene co-expression network analysis (WGCNA)

WGCNA was implemented via the R package “WGCNA” using TCGA-NSCLC bulk RNA-seq data to identify co-variant gene modules. We first determined an optimal soft threshold power (β=4) to achieve scale-free network topology. Weighted adjacency matrices were transformed into topological overlap matrices (TOMs), and dissimilarity measures were calculated for hierarchical clustering. Gene modules were identified using the dynamic tree cut algorithm, with the most ICD-correlated module (r=0.82, P<0.001) selected for downstream mechanistic exploration.

Differential gene analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the R package clusterProfiler to elucidate functional annotations and biological pathways. Genome-wide variation analysis was conducted via the R package gene set variation analysis (GSVA), leveraging the c2. cp. kegg. v7.4. symbols. gmt gene set from the Molecular Signatures Database (MSigDB) database to characterize pathway dynamics. Visualization of enrichment results was achieved using ggplot2, with statistical significance defined as P<0.05 to ensure scientific rigor.

Integrated ML prognostic model

To construct a robust prognostic model, we first identified candidate ICDRGs. The expression levels of these candidates were log₂(x+1) transformed and standardized prior to analysis. Using the TCGA-NSCLC cohort as the development set, we employed a Cox proportional hazards model as the base framework. Feature selection was performed using least absolute shrinkage and selection operator (LASSO) regression with 10-fold cross-validation to identify the most predictive gene signature and prevent overfitting. For comprehensive evaluation, an ensemble ML framework integrating 10 algorithms (e.g., Ridge regression, CoxBoost) generated 101 combinatorial models; their performance was compared based on the average concordance index (C-index) across cross-validation folds. For external validation, the final model coefficients derived from the TCGA cohort were directly applied to calculate the ICDRS for each patient in the independent GEO cohorts (GSE11969 and GSE68465). Model performance was assessed using discrimination metrics, including the C-index and time-dependent receiver operating characteristic (ROC) curves, and calibration metrics via calibration plots.

The final multivariable prediction model (ICDRS) is based on the expression levels of eight genes identified through the ML framework: NPC2, PDE4B, FCGRt, CTSH, CD69, HLA-DQA1, PTPRC, RGS1. These predictors were measured using RNA-seq (for the TCGA cohort) or microarray platforms (for the GEO validation cohorts) on primary tumor tissues obtained at diagnosis, prior to any systemic therapy. The ICDRS for each patient was calculated as a linear combination of the normalized expression values of these eight genes, weighted by their respective regression coefficients derived from the Cox model in the training set.

Prognostic nomogram construction

The primary outcome predicted by the ICDRS prognostic model was overall survival (OS). OS was defined as the time interval (in months) from the date of initial pathological diagnosis to the date of death from any cause. Patients who were alive at the last follow-up were considered censored observations. This outcome information was extracted directly from the clinical metadata provided by the TCGA and GEO databases.

Patients in the TCGA cohort were stratified into high-/low-risk subgroups using median ICDRS cutoff. Kaplan-Meier survival analysis (survminer package) confirmed significant survival disparity [log-rank (P<0.05)]. Time-dependent ROC curves (timeROC algorithm) demonstrated robust predictive accuracy [1-/3-/5-year area under the curve (AUC) >0.7]. Multivariate Cox regression validated ICDRS as an independent prognostic factor [hazard ratio (HR) =2.43, 95% confidence interval (CI): 1.58–3.74]. A dynamic nomogram integrating ICDRS and TNM staging enabled personalized survival probability prediction.

Immune microenvironment profiling

Immune infiltration was quantified using the GSVA extension package with ssGSEA-based enrichment of 22 immune signatures. Key steps included: (I) gene expression profile normalization and feature ranking; (II) cumulative distribution function-based immune score calculation; (III) standardized enrichment score matrix generation. Differential immune cell infiltration between subgroups was visualized via limma and ggpubr (FDR <0.05, Benjamini-Hochberg correction).

ICDRS-drug sensitivity correlation

Drug sensitivity profiling was conducted using the pRRophetic package (v0.5) to calculate half-maximal inhibitory concentrations (IC₅₀). Comparative analysis revealed significant differences in therapeutic responsiveness between high- and low-risk subgroups.

Correlation analysis between glycolytic genes and ICD markers

To further explore the association between glycolysis and immunity/ICD, we performed correlation analysis in the TCGA-lung adenocarcinoma (LUAD), GSE68465, and GSE81089 cohorts. Expression data acquisition and preprocessing are described in “Data collection” and “Integrated ML prognostic model”; probes in GSE68465 were annotated using GPL96 and converted to gene symbols, with multiple probes averaged. All expression data were log₂-transformed (if the maximum value >50) and then Z-score normalized. We analyzed five glycolytic genes (HK2, PKM2, SLC2A1, LDHA, PGK1) in relation to two sets of target genes: (I) ICD molecules (CALR, HMGB1, HSP90AA1, ATG5, PDIA3); (II) T cell effector/exhaustion markers (CD8A, GZMB, PRF1, IFNG, TOX, EOMES). Spearman correlation coefficients and their significance were calculated. Heatmaps displayed the correlation coefficients and significance annotations, with clustering based on Euclidean distance and complete linkage method, plotted using pheatmap.

Statistical analysis

All analyses were performed in R (v4.3.1). Chi-squared tests compared clinical features between training and validation cohorts. Wilcoxon tests evaluated non-normally distributed variables. Statistical significance was set at (P<0.05).

The analysis was conducted on a complete-case basis. For the construction of the final prognostic model, only samples with complete data on OS (time and status) and expression values for all genes comprising the final ICDRS signature were included. No imputation was performed for missing data in the predictors or outcome, as the rate of missingness for these core variables in the curated datasets was negligible. This approach ensures the integrity and reliability of the model coefficients and predictions.

Results

Molecular characteristics of ICD at the single-cell level in LUAD

We obtained scRNA-seq data encompassing 50,115 cells from three NSCLC cases (GSE189357 dataset). Batch effects were harmonized using the Harmony package, achieving effective integration of all samples (Figure S1A,S1B). Dimensionality reduction was performed via principal component analysis (PCA) and t-SNE on the top 2,000 variable genes. Cluster resolution was optimized at 0.3 through dendrogram analysis (Figure S1C), yielding 14 distinct clusters (Figure S1D). Cellular distribution patterns across samples revealed elevated tumor cell proportions in GSM5699785 with reduced immune infiltration (Figure 1A-1C). Using lineage-specific marker genes, we annotated 10 major cell populations: T cells (n=15,272), natural killer cells (NK cells) (n=6,954), B cells (n=5,641), tumor cells (n=7,429), DCs (n=5,095), macrophages (n=2,344), monocytes (n=1,022), and endothelial cells (n=284) (Figure 1B). A heatmap displaying the top four marker genes per cluster confirmed annotation specificity (Figure 1D). The “AddModuleScore” function (Seurat) quantified ICD activity using a curated 32-gene set across cell types (Figure 1E). DCs and lymphocytes exhibited significantly elevated ICD activity compared to other populations (DCs: median score =0.82 vs. other cells: 0.31, P<0.001; lymphocytes: 0.78 vs. other cells: 0.31, P<0.001; Figure 1F,1G). Stratification into high/low ICD activity groups revealed most pronounced differential expression in DCs and lymphocytes (Figure S1E).

Figure 1 Single-cell transcriptome profiling of lung adenocarcinoma and quantification of ICD activity. (A) Distribution of 11 cell populations in three NSCLC samples. (B) Clustering map of 10 major cell types annotated based on marker genes (T cells, NK cells, tumor cells, etc.). (C) Bar chart of cell population distribution ratios in three samples. (D) Heatmap of the top 4 marker gene expressions in each cell population. (E) Activity scores of ICD in each cell. (F) Violin plot of expression scores of 32-gene modules related to ICD in eight cell types. (G) Scatter plot of interaction intensity among eight cell types between high/low ICD expression groups. DC, dendritic cell; ICD, immunogenic cell death; NK, natural killer; NSCLC, non-small cell lung cancer.

Differences in intercellular communication between low- and high-ICD groups

Leveraging single-cell transcriptomic data and the CellChat algorithm, we systematically identified ligand-receptor/co-receptor gene pairs between low-ICD and high-ICD groups. By quantifying expression levels of these genes, we inferred interaction strengths of specific ligand-receptor signaling pathways across cell types, constructing comprehensive intercellular communication networks. While both groups exhibited interaction networks among eight major cell clusters, the high-ICD group demonstrated significantly denser network connectivity, particularly between tumor cells and DCs/lymphocytes, indicating enhanced communication frequency and interaction complexity (Figure 2A,2B). Cell-cell communication frequency and signal intensity were markedly elevated in the high-ICD group (Figure 2C). Heatmap visualization revealed distinct interaction patterns through a red-to-blue gradient (red: higher relative values; blue: lower values), with tumor cells in the high-ICD group showing strengthened interactions with all other cell types, amplifying immune response efficacy. Notably, DC-lymphocyte interactions exhibited intensified signaling (Figure 2D). Pathway analysis identified complete loss of glycoprotein Ib platelet alpha subunit (GP1BA) and TNF-related apoptosis-inducing ligand (TRAIL) signaling in the high-ICD group, while CD99 signaling showed group-specific dominance (Figure S2A). Among hub gene interactions, CD99 and chemokine (C-C motif) ligand (CCL) pathways demonstrated pronounced activity: CD99 signaling enhanced DC-NK/T cell and lymphocyte-NK/T cell crosstalk (Figure 2E,2F, Figure S2B), whereas CCL signaling strengthened DC-lymphocyte/NK cell and lymphocyte-NK cell interactions (Figure 2G,2H, Figure S2C). These pathways potentially regulate immune synapse formation, T cell activation, and cancer cell survival.

Figure 2 Intercellular communication networks in ICD subgroups. (A) Network diagram of interaction counts among 8 cell populations in the low ICD expression group. (B) Network diagram of interaction counts among 8 cell populations in the high ICD expression group. (C) Bar chart comparing communication frequency and signal intensity between high/low expression groups. (D) Heatmap of interaction differences, with left showing count differences and right showing intensity differences. (E,F) Network diagrams of CD99 signaling pathways in high/low expression groups. (G,H) Network diagrams of CCL signaling pathways in high/low expression groups. CCL, chemokine (c-c motif) ligand; DC, dendritic cell; ICD, immunogenic cell death; NK, natural killer.

Identification of key modules and genes related to ICD in batch RNA seq

We quantified ICD activity in the TCGA-NSCLC cohort using ssGSEA, calculating sample-specific enrichment scores based on a predefined ICD gene set to generate ICD activity indices for each NSCLC case. These indices were employed as core phenotypic features in WGCNA, constructing a scale-free topological network (Figure 3A) with an optimal soft thresholding power of 4 (scale-free R²=0.875; Figure S2D). Setting the minimum module gene count to 50 and a merge threshold of 0.15, we identified five co-expression modules (Figure 3B), among which the yellow module (MEyellow, n=150 genes) showed the strongest correlation with ICD scores in bulk RNA-seq (Pearson’s cor =0.87; Figure 3C). Module membership (MM) and gene significance (GS) analysis revealed a highly significant correlation (cor =0.9, P<0.001; Figure 3D), suggesting functional relevance of yellow module genes to ICD regulation. Using TCGA-NSCLC transcriptomic data, we systematically analyzed the regulatory roles of ICDRGs in NSCLC tumorigenesis. DESeq2-based differential expression analysis [FDR <0.05, |log₂fold change (FC)| ≥1]identified 4,173 significant differentially expressed genes (DEGs), including 1,891 tumor-upregulated and 2,282 downregulated genes (Figure 3E). Venn intersection with WGCNA co-expression modules prioritized 78 hub genes (Figure 3F), designated as ICD regulatory genes (ICDRGs) and validated via multi-omics integration (scRNA-seq and whole-transcriptome sequencing) as core regulators of ICD pathways. GO analysis revealed ICDRG enrichment in leukocyte adhesion activation, cytokine regulation, T cell activation, and antigen processing/presentation (biological processes); endocytic vesicles and major histocompatibility complex (MHC) class II complexes (cellular components); amide/peptide binding, MHC complex binding, and immune receptor activity (molecular functions) (Figure 3G). Univariate Cox regression identified 29 prognosis-significant ICDRGs (P<0.05; Figure S2E), with copy number variation (CNV) analysis showing ~10% amplification frequency for GMFG (Figure 3H). Protein-protein interaction (PPI) network mapping further characterized functional connectivity among these genes (Figure 3I).

Figure 3 Identification of key modules and genes related to ICD in bulk RNA-seq. (A) Construction of WGCNA co-expression network (based on ICD activity phenotype, soft threshold =4, R²=0.875, ensuring scale-free topology). (B) Module-ICD activity correlation analysis (yellow module significantly positively correlated with ICD score, cor =0.87). (C) GS analysis in the yellow module (MM highly correlated with GS, cor =0.9, P<0.001). (D) Scatter plot showing highly significant correlation between GS and MM in the MEblue module. (E) DEGs between tumor and normal tissues [DESeq2 analysis: 4,173 DEGs (1,891 upregulated, 2,282 downregulated), FDR <0.05, |log₂FC| ≥1]. (F) Screening of ICD-related hub genes (78 ICDRGs identified by intersecting yellow module genes of WGCNA with DEGs). (G) Functional enrichment analysis of ICDRGs. (H) CNV of key ICDRGs (CNV frequencies of 29 prognosis-related genes, GMFG amplification frequency ≈10%). (I) PPI network of ICDRGs (PPI core regulatory network of 29 prognosis genes). BP, biological process; CC, cellular component; CNV, copy number variation; DEG, differentially expressed gene; FDR, false discovery rate; FC, fold change; GS, gene significance; ICD, immunogenic cell death; ICDRG, immunogenic cell death-related gene; MF, molecular function; MHC, major histocompatibility complex; MM, module membership; PPI, protein-protein interaction; RNA-seq, RNA sequencing; WGCNA, weighted gene co-expression network analysis.

Construction of predictive signature based on integrated ML

We included 598 NSCLC patients from TCGA for model development and 441 patients from the combined GEO datasets (GSE11969 and GSE68465) for model validation. The key available clinical characteristics of the two groups of patients are summarized in Table 1.

To establish a precise and consistent ICD-related signature, we analyzed 29 prognostic genes derived from univariate Cox regression using 10 base algorithms [ridge regression, CoxBoost, random survival forest (RSF), etc.] ×10 feature selection methods ×1 resampling strategy, plus 1 LASSO-Cox stepwise model. The TCGA cohort served as the training set, with GSE11969 and GSE68465 as external validation cohorts. Error curve analysis via random forest modeling optimized feature selection (Figure 4A). Integration of LASSO and stepwise Cox regression achieved the highest mean C-index (0.69) across all validations. The CoxBoost + RSF model, incorporating only eight genes (NPC2, PDE4B, FCGRT, CTSH, CD69, HLA-DQA1, PTPRC, RGS1), demonstrated comparable predictive efficiency (Figure 4B). A refined 6-gene signature (Figure 4C) was used to calculate risk scores weighted by Cox regression coefficients, stratifying patients into high-/low-risk groups via median thresholds. Survival analysis confirmed significantly shorter OS in high-risk groups across training and validation cohorts (Figure 4D-4G). Confusion matrix evaluation revealed 62 true positives (TP) and 242 true negatives (TN) against 28 false negatives (FN) and 13 false positives (FP), yielding 88.1% accuracy [(62+242)/345] and 82.7% precision [62/(62+13)] (Figure 4H), demonstrating robust classification performance with high reliability for high-risk stratification.

Figure 4 Construction and validation of ICD prognostic signature based on integrated machine learning. (A) Error convergence curve of random forest (trend of OOB error with number of trees in 101-algorithm integrated modeling). (B) Comparison of prediction efficiency among multi-algorithm combinations (CoxBoost + RSF model achieved the highest average C-index =0.69 with 8 genes). (C) Screening of core prognostic genes (genes identified by LASSO + stepwise Cox joint screening). (D) Survival analysis in training set (TCGA). (E) Survival analysis in GSE68465 training set. (F) DSS analysis in training set (TCGA). (G) Survival analysis in GSE11969 training set. (H) Confusion matrix of risk prediction. C-index, concordance index; DSS, disease-specific survival; GBM, Gradient Boosting Machine; ICD, immunogenic cell death; OOB, out-of-bag; RSF, random survival forest; TCGA, The Cancer Genome Atlas.

To address the concern of platform heterogeneity, we further validated the ICDRS in an independent RNA-seq cohort GSE81089 (n=199). Applying the fixed Cox coefficients derived from the TCGA training set, we calculated the ICDRS for each patient and stratified them into high- and low-risk groups using the median cutoff. Kaplan-Meier analysis revealed that patients in the high-risk group had significantly shorter OS compared to those in the low-risk group (log-rank P=0.03; Figure S3A). Time-dependent ROC analysis yielded AUCs of 0.69, 0.71, and 0.68 for 1-, 2-, and 3-year survival, respectively (Figure S3B). These results demonstrate that the ICDRS maintains robust prognostic performance on an RNA-seq platform, effectively mitigating concerns about cross-platform generalizability.

Evaluation of ICDRS model

We evaluated the ICDRS model using ROC curves, demonstrating AUC values of 0.97 (1-year), 0.97 (3-year), and 0.98 (5-year) in the TCGA training cohort. External validation cohorts GSE11969 and GSE68465 showed AUCs of 0.74/0.72/0.69 and 0.80/0.84/0.83, respectively (Figure 5A-5C), confirming robust discriminative capacity. Heatmap visualization revealed elevated expression of all eight signature genes in the low-risk group across risk scores, clinical features, and molecular profiles (Figure 5D). Significant disparities in survival status, TNM staging, and nodal (N) distribution were observed between risk groups (Figure 5E). Stacked bar plots illustrated comparable T-stage (T1–T4) distributions between groups, with T2 predominance in both (Figure 5F). Patients with advanced T3–4 tumors, III–IV staging, or N2–3 nodal status exhibited significantly higher risk scores than those with T1–2, I–II staging, or N0–1 status (P<0.001), while age-stratified groups (≤60 vs. >60 years) showed no significant differences (Figure 5G). Kaplan-Meier analysis confirmed inferior survival probabilities for high-risk subgroups across all TNM categories: T1–2 vs. T3–4 (HR =2.8, P<0.001), M0 vs. M1 (HR =3.1, P<0.001), and early (I–II) vs. advanced (III–IV) staging (HR =2.5, P<0.001) (Figure 5H-5J). These results validate ICDRS as a potent prognostic stratifier across diverse clinical subtypes, providing robust decision support for personalized therapeutic strategies in NSCLC.

Figure 5 Clinical validation and subgroup analysis of ICDRS prognostic model. (A) ROC curve of training set (TCGA) showing specificity and sensitivity for predicting 1-, 3-, and 5-year OS. (B,C) ROC analysis in external validation sets (B: GSE11969 AUC =0.74/0.72/0.69; C: GSE68465 AUC =0.80/0.84/0.83). (D) Heatmap of correlations among risk scores, gene expressions, and clinical features [8 core genes (NPC2, PDE4B, etc.) showed significantly higher expressions in low-risk group]. (E) Distribution of clinicopathological features between high/low risk groups [higher proportions of advanced stages (III–IV), T3–4, and N2–3 in high-risk group, P<0.001]. (F) Comparison of T stage distribution (stacked bar chart: both groups dominated by T2, but higher proportion of T3–4 in high-risk group). (G) Differences in risk scores across clinical subgroups. (H-J) Survival validation in clinical subgroups (H: T stage subgroups; I: TNM stage subgroups; J: M stage subgroups. OS significantly shorter in high-risk patients across all subgroups, P<0.01). *, P<0.05; ***, P<0.001. AUC, area under the curve; CI, confidence interval; ICDRS, immunogenic cell death-related risk score; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; TNM, tumor-node-metastasis.

Establishment and validation of nomograms combined with clinical features

To evaluate ICDRS as an independent prognostic factor in NSCLC, we performed univariate and multivariate Cox regression analyses on TCGA-LUAD patients for OS, disease-specific survival (DSS), and progression-free survival (PFS). ICDRS demonstrated robust prognostic independence across all endpoints in both univariate (HR =2.94, P<0.001) and multivariate analyses (HR =2.67, P<0.001) after adjusting for age, gender, and TNM stage (Figure 6A-6C). External validation in the GSE68465 cohort replicated these findings, confirming ICDRS as an independent predictor of OS (multivariate HR =2.31, P=0.003; Figure 6D). A nomogram integrating risk scores with key clinical variables (age, gender, stage) achieved high concordance between predicted and observed 1-/3-/5-year OS probabilities (C-index =0.82; Figure 6E,6F). Calibration curves showed minimal deviation from ideal predictions (slope =0.96, intercept =0.02), with narrow 95% CIs across all timepoints. Comparative analysis revealed superior predictive performance of the risk score (AUC =0.89) over individual clinical features (age: AUC =0.62; gender: AUC =0.55; stage: AUC =0.75; Figure 6G). The nomogram’s prognostic accuracy significantly outperformed standalone clinical parameters (C-index difference =0.21 vs. stage-only model; Figure 6H), validating its clinical utility for NSCLC risk stratification.

Figure 6 Construction of clinical integration nomogram for ICDRS and independent prognostic validation. (A-C) Multivariate Cox regression analysis of TCGA-LUAD across multiple endpoints. (D) Independent validation in GSE68465 external cohort. (E) Nomogram for predicting 1/3/5-year survival probabilities (integrating ICDRS score, TNM stage, age, and gender). (F) Calibration curve of the nomogram. (G) Comparison of prediction efficiency among clinical features (risk score vs. stage/age/gender: 5-year AUC =0.84 vs. 0.72/0.58/0.53). (H) Comprehensive comparison of C-index [nomogram C-index =0.78 (0.74–0.82), significantly higher than single clinical features]. AUC, area under the curve; CI, confidence interval; C-index, concordance index; ICDRS, immunogenic cell death-related risk score; OS, overall survival; TCGA-LUAD, The Cancer Genome Atlas-lung adenocarcinoma; TNM, tumor-node-metastasis.

Molecular mechanisms of ICDRS in batch transcriptome

To elucidate the molecular mechanisms linking ICDRS with NSCLC prognosis, we performed GSEA and GSVA comparing high- vs. low-ICDRS groups. Low-risk patients exhibited enrichment in immune response pathways, including HALLMARK_ALLOGRAFT_REJECTION (immune recognition), HALLMARK_INTERFERON_GAMMA_RESPONSE (immune cell activation), and HALLMARK_IL6_JAK_STAT3_SIGNALING (inflammatory regulation) (Figure 7A-7C). High-risk patients showed dominant activation of cell cycle (HALLMARK_G2M_CHECKPOINT), transcriptional regulation (HALLMARK_E2F_TARGETS), and metabolic pathways (HALLMARK_GLYCOLYSIS). Correlation analysis confirmed strong associations between ICDRS scores and hallmark pathway activities (Figure 7D). Kaplan-Meier analysis revealed favorable prognosis for low-risk patients with elevated HALLMARK_IL6_JAK_STAT3_SIGNALING (HR =0.62, P=0.008) and HALLMARK_HEME_METABOLISM (HR =0.58, P=0.003) pathway scores (Figure 7E,7F), whereas high-risk patients with upregulated HALLMARK_E2F_TARGETS (HR =2.1, P<0.001) and HALLMARK_MYC_TARGETS_V1 (HR =1.9, P=0.002) exhibited significantly worse survival (Figure 7G,7H), demonstrating pathway-specific prognostic relevance of ICDRS stratification.

Figure 7 Enrichment of ICDRS-related molecular pathways and their prognostic values. (A-C) Hallmark pathway enrichment maps between high-/low-risk groups. (D) Correlation network between ICDRS and characteristic pathways. (E) Kaplan-Meier curve for the HALLMARK_IL6_JAK_STAT3_SIGNALING pathway. (F) Kaplan-Meier curve for the HALLMARK_HEME_METABOLISM pathway. (G) Kaplan-Meier curve for the HALLMARK_E2F_TA RGETS pathway. (H) Kaplan-Meier curve for the HALLMARK_MYC_TARGETS_V1 pathway. ICDRS, immunogenic cell death-related risk score.

We observed mixed correlations between glycolytic genes and ICD-related DAMPs (e.g., CALR, HSP90AA1), which may reflect the complex adaptive mechanisms of tumor cells under metabolic stress. Critically, however, a highly consistent and robust finding across all three cohorts was that glycolytic genes were significantly negatively correlated with T-cell effector molecules (CD8A, GZMB, PRF1, IFNG) and positively correlated with the immune exhaustion marker TOX (Figure S3C). This pattern strongly suggests that enhanced glycolytic activity is closely linked to the suppression of anti-tumor immune function, thereby providing a molecular-level explanation for the poor prognosis observed in the high ICDRS group.

Correlation between immune microenvironment and immune characteristics and ICDRS

We calculated immune scores, stromal scores, and ESTIMATE scores across ICDRS risk subgroups, with low-risk patients exhibiting significantly elevated immune (mean =1,852 vs. 982, P<0.001), stromal (mean =1,243 vs. 688, P<0.001), and ESTIMATE scores (mean =3,095 vs. 1,670, P<0.001) compared to high-risk counterparts (Figure 8A-8C). GSEA revealed enhanced activity of immune effector pathways in low-risk groups, including complement cascade (NES =2.1, FDR =0.02), FcγR-mediated phagocytosis (NES =1.9, FDR =0.03), and leukocyte transendothelial migration (NES =2.3, FDR =0.01) (Figure 8D). CIBERSORT analysis identified distinct immune infiltration patterns: low-risk tumors were enriched in resting lymphocytes (CD8⁺ T cells: 18.7% vs. 9.2%), quiescent DCs (6.5% vs. 2.1%), and NK cells (7.8% vs. 3.4%), while high-risk groups showed M0 macrophage predominance (23.6% vs. 8.9%) (Figure 8E). ssGSEA and xCell algorithms validated these infiltration patterns (concordance correlation coefficient =0.86, P<0.001; Figure 8F). Correlation analysis demonstrated negative associations between eight ICDRS signature genes and M0 macrophage infiltration (mean ρ=−0.62, P<0.001), alongside positive correlations with resting DC populations (mean ρ=0.71, P<0.001) (Figure 8G).

Figure 8 Tumor microenvironment characteristics and immune cell landscape in ICDRS risk subgroups. (A-C) Comparison of microenvironment scoring systems (A: immune score; B: stromal score; C: ESTIMATE score). (D) GSEA enrichment of immune-active pathways. (E) Differences in immune cell infiltration by CIBERSORT. (F) Validation of infiltration signatures by multiple algorithms. (G) Correlation network between hub genes and immune cells. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant. GSEA, gene set enrichment analysis; ICDRS, immunogenic cell death-related risk score.

Correlation analysis between ICDRS and drug sensitivity and validation of gene expression

In advanced NSCLC, we assessed ICDRS-based risk stratification for predicting responses to TKIs and chemotherapeutic agents. Low-risk patients exhibited significantly lower IC₅₀ for gefitinib (mean =1.2 vs. 3.8 µM, P<0.001) and sunitinib (mean =0.9 vs. 2.5 µM, P<0.001), with risk scores showing positive correlations to drug IC₅₀ values (ρ=0.72 and ρ=0.68, respectively; Figure 9A,9B), indicating enhanced TKI sensitivity in low-risk subgroups. No significant intergroup differences were observed for vinorelbine, a cell cycle-specific agent (Figure 9C). Intriguingly, high-risk patients demonstrated greater therapeutic benefit from docetaxel, evidenced by lower IC₅₀ values (mean =0.6 vs. 1.4 µM, P=0.002; Figure 9D). These findings validate ICDRS as a predictive biomarker for guiding targeted therapy in NSCLC.

Figure 9 Predictive efficacy of ICDRS risk score for sensitivity to targeted and chemotherapeutic drugs in NSCLC. (A) Gefitinib IC₅₀ significantly lower in low-risk group of molecular subtypes (risk score positively correlated with IC₅₀). (B) Sunitinib sensitivity analysis showed better drug response in low-risk group (IC₅₀ positively correlated with risk score). (C) No significant sensitivity difference of Vinorelbine between risk subgroups (no statistical significance in IC₅₀). (D) Docetaxel IC₅₀ significantly higher in low-risk group (high-risk patients showed more significant chemotherapy benefit). **, P<0.01; ****, P<0.0001; ns, not significant. IC₅₀, half-maximal inhibitory concentration; ICDRS, immunogenic cell death-related risk score; NSCLC, non-small cell lung cancer.

Discussion

NSCLC remains a leading cause of cancer-related mortality worldwide, with its high incidence and mortality rates posing a major public health challenge. Despite significant advancements in early screening, surgical techniques, and targeted therapies over the past decade, most patients are diagnosed at advanced stages, resulting in suboptimal treatment outcomes and shortened survival (17). Current therapeutic strategies—including surgical resection, radiotherapy, chemotherapy, and targeted therapies—face limitations such as drug resistance and high recurrence rates, which restrict their ability to sustainably prolong patient survival (18). These challenges underscore the urgent need to develop novel therapeutic approaches, particularly through leveraging insights into tumor immunogenicity and microenvironmental dynamics to advance personalized therapeutic paradigms.

Cell death encompasses critical physiological processes including apoptosis, necrosis, autophagic cell death, and ferroptosis. Apoptosis, a programmed cell death mechanism, operates through caspase activation via intrinsic/extrinsic pathways to achieve orderly cellular dismantling without inflammation (19). Necrosis involves unregulated cell demise with membrane rupture and inflammatory cascades. Autophagic cell death degrades cellular components via lysosomal pathways under nutrient deprivation or stress (20). Ferroptosis, an iron-dependent lipid peroxidation-driven death modality, disrupts membrane integrity (21). These mechanisms collectively shape TME dynamics, influencing cancer progression. ICD, conceptualized by Guido Kroemer’s team in 2005, redefined apoptotic signaling by demonstrating chemotherapy-induced release of DAMPs that activate DCs and T cells, converting tumors into “in situ vaccines” to enhance ICI efficacy, particularly in cold tumors (22). Recent studies highlight an iridium(III) complex inducing ICD in NSCLC, triggering durable anti-tumor CD8⁺ T cell responses and Foxp3⁺ T cell depletion in murine models (23). Metabolic reprogramming studies reveal G6PD, a pentose phosphate pathway enzyme maintaining NADPH/redox homeostasis, as an ICD modulator—its inhibition amplifies immunotherapy efficacy by inducing oxidative stress-mediated ICD in lung cancer (10). Prognostic biomarker research has identified ICD-related signatures (e.g., IL-1β, IFN-γ) correlating with improved survival (24). ML-based ICD gene models enable personalized treatment stratification (25). However, current studies predominantly rely on bulk transcriptomic single-omics analyses with simplistic methodologies, failing to comprehensively decode ICD’s molecular landscape in gastric cancer (GC). Our work delineates ICD heterogeneity across TME cell subtypes and links ICD activity to survival outcomes in NSCLC, providing a framework for precision immunotherapy to improve therapeutic efficacy and survival rates (26,27).

Our study employed a novel multi-omics integrative approach to advance ICD research. Leveraging single-cell transcriptomic profiling, we dissected the spatial heterogeneity of ICD activity and identified core regulatory networks through differential gene enrichment algorithms. WGCNA on TCGA cohorts constructed scale-free topological structures, pinpointing ICD-associated hub gene clusters with correlations exceeding r=0.9. To enhance analytical rigor, we developed an ensemble algorithmic framework integrating 101 combinatorial variants across 10 ML paradigms, systematically optimizing prognostic signatures validated for cross-platform robustness in independent cohorts (GSE11969/GSE68465, C-index =0.82). Importantly, the prognostic value of ICDRS was not only validated in two microarray-based GEO cohorts (GSE11969 and GSE68465) but also in an independent RNA-seq cohort (GSE81089). The consistent HRs and time-dependent AUC values across these technologically distinct platforms underscore the robustness of the eight-gene signature to platform-specific technical variations. This is particularly relevant for clinical translation, where diagnostic assays may be implemented on different platforms. The successful replication in GSE81089, which uses Illumina HiSeq technology identical to the training set, directly addresses concerns regarding platform heterogeneity and reinforces the generalizability of our findings. Multidimensional clinical validation demonstrated superior discriminative power in TNM stage stratification (P<0.001), underscoring translational potential. Mechanistic interrogation via functional enrichment revealed signature gene clusters governing immune response modulation and transmembrane receptor signaling. Ultimately, we established a predictive system integrating immunological infiltration scoring and drug sensitivity profiling, pioneering a precision oncology paradigm to accelerate personalized therapeutic development.

We identified nine distinct cell types from NSCLC scRNA-seq data, with DCs and monocytes exhibiting markedly elevated ICD scores, highlighting their pivotal roles in NSCLC pathogenesis. This underscores a robust association between ICD activity and these immune subsets, reflecting their functional specialization within the TME. Prior single-cell analyses revealed that ICD-high tumor cells upregulate co-stimulatory molecules (CD86, HLA-DR) on DC surfaces and promote DC migration to lymph nodes via the CCL5-CCR7 axis (28). Stratifying cells into high- and low-ICD groups using ssGSEA, we observed significant transcriptional divergence in the high-ICD cohort. These findings illuminate potential links between ICD activity and tumor progression while nominating novel biomarkers for immunotherapy optimization. Pathway enrichment analysis implicated apoptosis (e.g., CASP3/BAX-mediated mitochondrial pathways) and immune response activation (e.g., TLR4/STING/CXCL9 signaling) in high-ICD phenotypes. Mechanistically, apoptotic DAMPs released via CASP3/BAX activation synergize with HMGB1-TLR4-MYD88-driven DC maturation and CD8⁺ T cell recruitment, amplifying ICD effects to enhance anti-tumor immunity (29,30).

We identified an 8-gene ICD-related signature (NPC2, PDE4B, FCGRT, CTSH, CD69, HLA-DQA1, PTPRC, RGS1) through integrated LASSO and stepwise Cox regression, with the CoxBoost + RSF model achieving a C-index of 0.69 across training and validation cohorts. This signature provides novel biomarkers for ICD stratification and a foundation for personalized therapeutic strategies. Transcriptomic profiling of ICDRS risk subgroups revealed significant associations between ICDRS scores and Hallmark pathway GSVA scores, with low-risk patients enriched in immune recognition/activation pathways (e.g., interferon-γ response) and high-risk patients in cell cycle regulation (e.g., E2F targets) and glycolysis. Tumor microenvironmental inflammation, mediated by tumor-associated macrophages (TAMs) and cytokines like IL-1β/TNF-α, amplifies ICD efficacy by enhancing immune cell infiltration (31-33). Cell cycle phase-specific ICD modulation was observed: S/G2-phase DNA damage triggers non-immunogenic mitotic catastrophe (34) while G1-phase damage enables error-prone non-homologous end joining (NHEJ), potentiating immunogenic interferon responses (32). Dysregulated cyclin-dependent kinases (CDKs) (CDKs 1/2/5) exhibit dual roles—CDK inhibition (e.g., Dinaciclib) induces ICD to overcome interferon-γ-mediated resistance (35), while CDK5 inhibitors reduce PD-L1 expression to attenuate immune evasion and ICD (36). These findings elucidate cell cycle-ICD crosstalk and therapeutic vulnerabilities in NSCLC.

Our ICDRS model identified eight functionally diverse genes (NPC2, PDE4B, FCGRT, CTSH, CD69, HLA-DQA1, PTPRC, RGS1) with distinct roles in ICD and antitumor immunity. NPC2, a cholesterol transporter, correlates positively with tumor-infiltrating macrophages and DCs, modulating immune responses via DAMPs release (37,38). PDE4B (phosphodiesterase 4B) degrades cAMP to regulate immune cell activation, where its inhibition enhances macrophage M1 polarization and antitumor activity (39,40). FCGRT (neonatal Fc receptor) influences IgG metabolism and tumor immune evasion, with overexpression linked to immunosuppressive microenvironments and reduced T-cell infiltration (41,42). CTSH (Cathepsin H) is an important lysosomal protease that mainly participates in the degradation and recycling of intracellular proteins. The activity of CTSH is related to the sensitivity of tumor cells to immunotherapy. CTSH can enhance the immunogenicity of tumor cells by promoting the release and processing of intracellular antigens. HLA-DQA1 is a part of the human MHC class II molecule, primarily responsible for antigen presentation and activation of T cells. Research has shown that high expression of HLA-DQA1 is associated with the immune escape mechanism of tumor cells, which may affect tumor progression and patient prognosis by regulating immune cell infiltration and function in the TME. For example, in high-grade soft tissue sarcoma, the expression of HLA-DQA1 is significantly positively correlated with patient survival, suggesting that it may serve as a prognostic marker (43). In melanoma patients, high expression of HLA-DQA1 is associated with a good response to ICIs (44). Through single-cell analysis, it can be seen that these genes are mainly highly expressed in DCs and lymphocytes, which is consistent with the expression characteristics observed in previous ICD evaluations. These findings further demonstrate the close relationship between DCs, ICD, and non-small cell carcinoma tumor cells. In the TME, DCs activate and enhance the anti-tumor activity of CD8⁺ T cells by recognizing DAMPs released by tumor cells (45). Research has shown that the function of DCs is influenced by the TME, and the immune escape mechanism of tumors is often achieved by inhibiting the function of DCs (46). Therefore, enhancing the activity and function of DCs, or restoring their anti-tumor ability through immunotherapy, has become an important strategy to improve the efficacy of tumor immunotherapy (47). For example, the combination therapy strategy of ICIs and ICD inducers has shown the potential to significantly enhance anti-tumor immune responses in various tumor models (38).

To further explore the potential associations between glycolytic pathway activation, immune suppression, and ICD regulation in the high-risk group, we systematically evaluated expression correlations between core glycolytic genes (HK2, PKM2, SLC2A1, LDHA, PGK1), ICD key molecules (CALR, HMGB1, HSP90AA1, ATG5, PDIA3), and T cell effector markers (CD8A, GZMB, PRF1, IFNG) across three independent cohorts (TCGA-LUAD, GSE68465, GSE81089). The results showed that glycolytic genes were significantly positively correlated with ICD-related DAMPs (e.g., CALR, ATG5, HSP90AA1; ρ=0.2–0.5, P<0.001) but significantly negatively correlated with T cell effector molecules (e.g., CD8A, GZMB, IFNG; ρ=−0.2 to −0.5, P<0.001), a pattern consistent across all three cohorts, revealing a transcriptomic-level quantitative coupling between glycolytic reprogramming, ICD molecule expression, and T cell dysfunction. As a key glycolytic enzyme, high PGK1 expression promotes lactate accumulation, TME acidification, CD8⁺ T cell exhaustion, and reduced anti-PD-1 efficacy (48,49), while PGK1 inhibition reverses these effects; additionally, Ito et al. demonstrated that co-positive GLUT1 and PKM2 expression correlates with worse disease-free survival (HR =2.1, P<0.001) and invasive features in resected NSCLC (50), highlighting glycolytic enzymes as drivers of malignant progression. Notably, glycolytic genes were positively correlated with ICD positive regulators (e.g., CALR, ATG5), suggesting glycolytic activation may induce ICD-related molecule expression, but the immunogenic functions of these DAMPs depend on subcellular localization (e.g., CALR membrane translocation) and post-translational modifications (e.g., phosphorylation, acetylation) (51-53), meaning transcriptional upregulation alone may not generate effective immunogenic signals—particularly in a glycolysis-driven immunosuppressive microenvironment—explaining why the high-risk group exhibits immunosuppression despite ICD molecule expression. Collectively, our multi-cohort correlation analysis, combined with existing functional evidence, constructs a mechanistic model linking metabolic reprogramming to immune escape in the high-risk ICDRS subgroup: enhanced glycolytic activity (characterized by upregulated HK2, PGK1, SLC2A1, etc.) leads to lactate accumulation and microenvironmental acidification, inducing CD8⁺ T cell exhaustion and impaired effector function, such that tumor cells expressing ICD-related molecules fail to activate effective anti-tumor immunity, ultimately resulting in immune surveillance escape and poor prognosis; subsequent studies should validate this model in vitro and in vivo via CRISPR-Cas9-mediated knockout of core ICDRS genes (e.g., PDE4B, RGS1), combined with glycolytic stress testing, T cell-tumor cell co-culture, and ICD function assays.

In addition to glycolysis, the high-risk group exhibited significant enrichment of the G2M checkpoint pathway, consistent with elevated proliferative activity and unfavorable prognosis. Mounting evidence supports crosstalk between cell cycle progression and ICD. CDK4/6 kinases directly phosphorylate the p53 family member p73 at Thr86, sequestering it in the cytoplasm; CDK4/6 inhibition restores p73 nuclear translocation and upregulates death receptor 5 (DR5), thereby promoting ICD via the extrinsic apoptotic pathway. DR5 depletion abolishes the ICD-sensitizing effects of CDK4/6 inhibitors in combination with chemotherapy or immunotherapy, confirming a direct mechanistic link between CDK4/6 activity and ICD competence (54). Moreover, pharmacological inhibition of CDK12/13 induces hallmarks of ICD, including HMGB1 release, ATP secretion, and CALR membrane translocation, while enhancing anti-tumor immunity (55,56). These observations indicate that CDKs modulate ICD through distinct substrate-specific mechanisms. We therefore propose that heightened G2M checkpoint activity in the high-risk subgroup reflects hyperactive CDK4/6 and CDK12 signaling. On one hand, sustained CDK activity drives unchecked cell cycle progression and tumor proliferation; on the other, it suppresses DR5-mediated ICD via p73 phosphorylation, thereby diminishing tumor immunogenicity and enabling immune escape. This model aligns with the enriched interferon-γ response and increased effector T-cell infiltration observed in the low-risk group, where restrained CDK activity preserves ICD functionality and anti-tumor immune surveillance. Future studies will validate this model in NSCLC models via CDK4/6 or CDK12 ablation, combined with ICD functional assays and T-cell co-culture systems. Targeting CDKs in combination with ICD inducers may represent a promising therapeutic strategy for high-risk patients.

Our study validated the clinical relevance of ICDRS stratification, with high-ICDRS subgroups demonstrating superior survival outcomes compared to low-ICDRS counterparts. High-ICDRS tumors exhibited elevated infiltration of tumor-infiltrating lymphocytes (TILs) and DCs, correlating positively with prolonged survival, while low-ICDRS groups showed higher tumor burden and immunosuppressive microenvironments characterized by reduced immune infiltration. Single-cell resolution confirmed ICDRS-associated genes predominantly localized to lymphocytes and DCs. Tumor-mediated DC dysfunction—via immunosuppressive cytokine secretion—compromises antigen presentation and T-cell activation, facilitating immune evasion (57), though DC abundance remains prognostically favorable (58). Low-ICDRS tumors harbored enriched M0 macrophages, traditionally considered “resting” precursors to polarized M1/M2 subtypes. Emerging evidence implicates M0 macrophages as functional components of TAM populations, potentially driving protumoral processes. In hepatocellular carcinoma (HCC), M0 macrophage infiltration significantly exceeds normal liver tissue, with 35 M0 macrophage-related genes (M0RGs) linked to poor prognosis (59). Prognostic modeling revealed M0RG signatures as independent predictors of adverse outcomes, suggesting M0 macrophages may promote tumor progression through gene regulatory networks. These findings illuminate NSCLC progression mechanisms and provide a scientific foundation for developing precision immunotherapies targeting ICD-related pathways, DC functional restoration, and macrophage polarization dynamics.

ICDRS demonstrates unique clinical utility in personalizing NSCLC treatment through drug sensitivity correlations, offering novel strategies to overcome therapeutic resistance. In advanced NSCLC, TME heterogeneity and genetic evolution frequently reduce conventional therapy efficacy, while ICDRS enables molecular subtype-driven therapeutic optimization by quantifying ICDRG expression patterns. TKIs gefitinib and sunitinib exhibit significantly lower IC₅₀ in low-risk subgroups (P<0.001), with risk scores positively correlating to IC₅₀ values (ρ=0.72–0.68). Mechanistically, elevated FCGRT (immunoglobulin transport receptor) expression in low-risk tumors enhances antibody-dependent cellular cytotoxicity (ADCC), while NPC2 potentiates TKI efficacy by promoting lysosomal cholesterol metabolism-mediated drug endocytosis. Conversely, taxane chemotherapeutic docetaxel demonstrates superior efficacy in high-risk subgroups (P=0.003), attributable to enrichment of cell cycle pathways (e.g., HALLMARK_G2M_CHECKPOINT) and CDK1/2 hyperactivity. High-risk tumors exhibit heightened vulnerability to cycle-specific agents, whereas active immunometabolic networks (e.g., IL-6/JAK/STAT3) in low-risk subgroups may attenuate docetaxel cytotoxicity via mitochondrial apoptosis suppression. Notably, vinorelbine shows comparable efficacy across subgroups, suggesting tubulin inhibitor mechanisms operate independently of ICD status, supporting ICDRS-guided chemotherapy selection.

While our multi-omics and ML approach systematically characterized ICD features in NSCLC, several limitations warrant consideration: (I) although we have validated the ICDRS across multiple independent cohorts encompassing both microarray (GSE11969, GSE68465) and RNA-seq (GSE81089) platforms, all validations were performed on retrospective publicly available datasets. Prospective studies with standardized protocols and homogeneous patient populations are warranted to confirm the clinical utility and generalizability of the signature before its translation into practice; (II) mechanistic interpretations of ICDRG functional enrichments (e.g., glycolysis pathway associations) lack experimental validation, with causal relationships to ICD remaining unconfirmed via metabolomic profiling or gene knockout studies; (III) immune microenvironment analyses omitted spatial transcriptomics to resolve spatial cellular dynamics; (IV) although the ICDRS shows robust prognostic value and has been validated across multiple independent cohorts, our assessment of ICD activity is based exclusively on transcriptomic data. It is important to note that functional ICD is largely governed by non-transcriptional events, including post-translational modifications and epigenetic regulation, which control the subcellular localization, translocation, and release of DAMPs (60-62). For instance, the extracellular secretion of HMGB1 is tightly regulated by acetylation at specific lysine residues, which determines its nuclear export and immunostimulatory activity (63-65). Similarly, surface exposure of calreticulin—the canonical “eat me” signal for DC activation—depends on phosphorylation cascades mediated by protein kinase C and other kinases, rather than simply mRNA expression levels (66). Epigenetic mechanisms, including DNA methylation and histone modifications, further modulate the activity of ICD-related molecules, adding another layer of regulation beyond transcript abundance (67,68). Therefore, while a transcriptome-based ICD signature enables convenient and scalable quantification in clinical cohorts, it may not fully reflect the actual immunogenicity of dying tumor cells. This constraint should be considered when interpreting ICDRS performance across datasets. Future studies integrating proteomic, post-translational modification, and functional DAMP measurements will help establish more comprehensive and mechanistically accurate ICD evaluation systems for NSCLC. These gaps highlight the need for orthogonal multi-omics validation and functional interrogation to refine ICD-driven therapeutic strategies.

Conclusions

This study identified a cluster of immune-related differentially expressed genes (ICDRGs) in the TCGA-NSCLC cohort through WGCNA and differential expression analysis. By integrating ML approaches, we developed an 8-gene immune-related differentially expressed gene signature (ICDRS) that demonstrated robust prognostic predictive performance across training and validation cohorts. The ICDRS demonstrated robust prognostic performance across both microarray and RNA-seq platforms in multiple independent cohorts, supporting its potential for broad clinical application. Comprehensive multi-omics investigations—including genomic, single-cell transcriptomic, bulk transcriptomic, and drug sensitivity profiling—further validated the potential of ICDRS as a novel molecular framework to guide precision immunotherapy and personalized therapeutic strategies for NSCLC. Future studies incorporating expanded cohorts and functional validation of ICDRG mechanisms will optimize this model for clinical translation, ultimately advancing patient-tailored therapeutic paradigms in oncology.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2682/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2682/prf

Funding: This work was supported by the Clinical Research Fund of Shandong Medical Association (grant number: yxh2022zx02031).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2682/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Hendriks LEL, Remon J, Faivre-Finn C, et al. Non-small-cell lung cancer. Nat Rev Dis Primers 2024;10:71. [Crossref] [PubMed]
2. Yan Y, Sun Y, Guo X, et al. Immune Evasion Mechanism of Neurotropic Viruses. Rev Med Virol 2024;34:e2589. [Crossref] [PubMed]
3. Frelaut M, Besse B. Immunotherapy for non-small-cell lung cancer. La Revue du praticien 2021;71:1057-9.
4. Su S, Chen F, Xu M, et al. Recent advances in neoantigen vaccines for treating non-small cell lung cancer. Thorac Cancer 2023;14:3361-8. [Crossref] [PubMed]
5. Xie D, Wang Q, Wu G. Research progress in inducing immunogenic cell death of tumor cells. Front Immunol 2022;13:1017400. [Crossref] [PubMed]
6. Yang X, Yang J, Gu X, et al. (-)-Guaiol triggers immunogenic cell death and inhibits tumor growth in non-small cell lung cancer. Mol Cell Biochem 2023;478:1611-20. [Crossref] [PubMed]
7. Alvarez-Valadez K, Sauvat A, Diharce J, et al. Lysosomal damage due to cholesterol accumulation triggers immunogenic cell death. Autophagy 2025;21:934-56. [Crossref] [PubMed]
8. Flieswasser T, Van Loenhout J, Freire Boullosa L, et al. Clinically Relevant Chemotherapeutics Have the Ability to Induce Immunogenic Cell Death in Non-Small Cell Lung Cancer. Cells 2020;9:1474. [Crossref] [PubMed]
9. Liu P, Zhao L, Senovilla L, et al. In Vivo Imaging of Orthotopic Lung Cancer Models in Mice. Methods Mol Biol 2021;2279:199-212. [Crossref] [PubMed]
10. Nakamura M, Magara T, Yoshimitsu M, et al. Blockade of glucose-6-phosphate dehydrogenase induces immunogenic cell death and accelerates immunotherapy. J Immunother Cancer 2024;12:e008441. [Crossref] [PubMed]
11. Joo MS, Pyo KH, Chung JM, et al. Artificial intelligence-based non-small cell lung cancer transcriptome RNA-sequence analysis technology selection guide. Front Bioeng Biotechnol 2023;11:1081950. [Crossref] [PubMed]
12. Feng X. Integrative analysis of GWAS and transcriptomics data reveal key genes for non-small lung cancer. Med Oncol 2023;40:270. [Crossref] [PubMed]
13. Wang M, Dai X, Yang X, et al. Serum Protein Fishing for Machine Learning-Boosted Diagnostic Classification of Small Nodules of Lung. ACS Nano 2024;18:4038-55. [Crossref] [PubMed]
14. Pan H, Liu P, Zhao L, et al. Immunogenic cell stress and death in the treatment of cancer. Semin Cell Dev Biol 2024;156:11-21. [Crossref] [PubMed]
15. Fang K, Yuan S, Zhang X, et al. Regulation of immunogenic cell death and potential applications in cancer therapy. Front Immunol 2025;16:1571212. [Crossref] [PubMed]
16. Wang X, Wu S, Liu F, et al. An Immunogenic Cell Death-Related Classification Predicts Prognosis and Response to Immunotherapy in Head and Neck Squamous Cell Carcinoma. Front Immunol 2021;12:781466. [Crossref] [PubMed]
17. Skřičková J, Kadlec B, Venclíček O. Non-small cell lung cancer. Vnitrni lekarstvi 2018;63:861-74.
18. Chen P, Liu Y, Wen Y, et al. Non-small cell lung cancer in China. Cancer Commun (Lond) 2022;42:937-70. [Crossref] [PubMed]
19. Park W, Wei S, Kim BS, et al. Diversity and complexity of cell death: a historical review. Exp Mol Med 2023;55:1573-94. [Crossref] [PubMed]
20. Yang J, Hu S, Bian Y, et al. Targeting Cell Death: Pyroptosis, Ferroptosis, Apoptosis and Necroptosis in Osteoarthritis. Front Cell Dev Biol 2021;9:789948. [Crossref] [PubMed]
21. Dou J, Liu X, Yang L, et al. Ferroptosis interaction with inflammatory microenvironments: Mechanism, biology, and treatment. Biomed Pharmacother 2022;155:113711. [Crossref] [PubMed]
22. Casares N, Pequignot MO, Tesniere A, et al. Caspase-dependent immunogenicity of doxorubicin-induced tumor cell death. J Exp Med 2005;202:1691-701. [Crossref] [PubMed]
23. Wang L, Guan R, Xie L, et al. An ER-Targeting Iridium(III) Complex That Induces Immunogenic Cell Death in Non-Small-Cell Lung Cancer. Angew Chem Int Ed Engl 2021;60:4657-65. [Crossref] [PubMed]
24. Wen Y, You D, Bai Y, et al. Development and validation of immunogenic cell death-applied prediction model for esophageal carcinoma. Am J Cancer Res 2023;13:2104-15.
25. Liu Z, Luo Y, Su L, et al. Identification of immunogenic cell death-related subtypes used for predicting survival and immunotherapy of endometrial carcinoma through a bioinformatics analysis. Medicine (Baltimore) 2023;102:e34571. [Crossref] [PubMed]
26. Garg AD, Dudek-Peric AM, Romano E, et al. Immunogenic cell death. Int J Dev Biol 2015;59:131-40. [Crossref] [PubMed]
27. Troitskaya OS, Novak DD, Richter VA, et al. Immunogenic Cell Death in Cancer Therapy. Acta Naturae 2022;14:40-53. [Crossref] [PubMed]
28. Dai Z, Wang Y, Sun N, et al. Characterizing ligand-receptor interactions and unveiling the pro-tumorigenic role of CCL16-CCR1 axis in the microenvironment of hepatocellular carcinoma. Front Immunol 2023;14:1299953. [Crossref] [PubMed]
29. Ma X, Shan H, Chen Z, et al. Programmed cell death-related prognostic genes mediate dysregulation of the immune microenvironment in triple-negative breast cancer. Front Immunol 2025;16:1563630. [Crossref] [PubMed]
30. Zhao F, Jiang X, Li Y, et al. Characterizing tumor biology and immune microenvironment in high-grade serous ovarian cancer via single-cell RNA sequencing: insights for targeted and personalized immunotherapy strategies. Front Immunol 2024;15:1500153. [Crossref] [PubMed]
31. Romero Fernandez J, Cordoba Largo S, Benlloch Rodriguez R, et al. The Effects of Gynecological Tumor Irradiation on the Immune System. Cancers (Basel) 2024;16:2804. [Crossref] [PubMed]
32. Tian Z, Li X, Jiang D. Analysis of immunogenic cell death in atherosclerosis based on scRNA-seq and bulk RNA-seq data. Int Immunopharmacol 2023;119:110130. [Crossref] [PubMed]
33. Lai JM, Chen PL, Shi QY, et al. A Self-Delivery Nanodrug Simultaneously Inhibits COX-2/PGE(2) Mediated Inflammation and Downregulates PD-L1 to Boost Photoimmunotherapy. Adv Healthc Mater 2024;13:e2400367. [Crossref] [PubMed]
34. Szmyd R, Casolin S, French L, et al. Homologous recombination promotes non-immunogenic mitotic cell death upon DNA damage. Nat Cell Biol 2025;27:59-72. [Crossref] [PubMed]
35. Janssens S, Rennen S, Agostinis P. Decoding immunogenic cell death from a dendritic cell perspective. Immunol Rev 2024;321:350-70. [Crossref] [PubMed]
36. Liu J, Zhao Z, Qiu N, et al. Co-delivery of IOX1 and doxorubicin for antibody-independent cancer chemo-immunotherapy. Nat Commun 2021;12:2425. [Crossref] [PubMed]
37. Liu W, Liu L, Kuang T, et al. Cholesterol metabolism-related genes predict immune infiltration and prognosis in gastric cancer patients. J Cancer 2025;16:2087-102. [Crossref] [PubMed]
38. Shi F, Huang X, Hong Z, et al. Improvement strategy for immune checkpoint blockade: A focus on the combination with immunogenic cell death inducers. Cancer Lett 2023;562:216167. [Crossref] [PubMed]
39. Hernández ÁP, Juanes-Velasco P, Landeira-Viñuela A, et al. Restoring the Immunity in the Tumor Microenvironment: Insights into Immunogenic Cell Death in Onco-Therapies. Cancers (Basel) 2021;13:2821. [Crossref] [PubMed]
40. Luo Y, Linghu M, Luo X, et al. Remodeling tumor immunosuppressive microenvironment through dual activation of immunogenic panoptosis and ferroptosis by H(2)S-amplified nanoformulation to enhance cancer immunotherapy. Acta Pharm Sin B 2025;15:1242-54. [Crossref] [PubMed]
41. Song DH, Yang J, Kim CH, et al. FcRn Expression in Endometrial Cancer and Its Association with Clinicopathologic Features. Diagnostics (Basel) 2023;13:3660. [Crossref] [PubMed]
42. Xue F, Yang L, Dai B, et al. Bioinformatics profiling identifies seven immune-related risk signatures for hepatocellular carcinoma. PeerJ 2020;8:e8301. [Crossref] [PubMed]
43. Bae JY, Choi KU, Kim A, et al. Evaluation of immune-biomarker expression in high-grade soft-tissue sarcoma: HLA-DQA1 expression as a prognostic marker. Exp Ther Med 2020;20:107. [Crossref] [PubMed]
44. Pahkuri S, Katayama S, Valta M, et al. The effect of type 1 diabetes protection and susceptibility associated HLA class II genotypes on DNA methylation in immune cells. HLA 2024;103:e15548. [Crossref] [PubMed]
45. Chen Y, Xiong T, Gu M, et al. Pyrazolone-Modified Photosensitizers for Precise Cell Membrane Rupture to Enhance Cancer Immunotherapy. J Am Chem Soc 2025;147:17915-25. [Crossref] [PubMed]
46. Ribatti D. Different subpopulations of macrophages, neutrophils, mast cells, and fibroblasts are involved in the control of tumor angiogenesis. Front Med (Lausanne) 2024;11:1481609. [Crossref] [PubMed]
47. Ji J, Liu S, Liao X, et al. Research progress on ferroptosis in tumor immunity. Xi Bao Yu Fen Zi Mian Yi Xue Za Zhi 2024;40:362-6.
48. Zhang K, Sun L, Kang Y. Regulation of phosphoglycerate kinase 1 and its critical role in cancer. Cell Commun Signal 2023;21:240. [Crossref] [PubMed]
49. Wang P, Wang Q, Yang X, et al. Targeting the Glycolytic Enzyme PGK1 to Inhibit the Warburg Effect: A New Strategy for Keloid Therapy. Plast Reconstr Surg 2023;151:970e-80e.
50. Ito R, Yashiro M, Tsukioka T, et al. GLUT1 and PKM2 may be useful prognostic predictors in patients with non-small cell lung cancer following curative R0 resection. Oncol Lett 2023;25:129. [Crossref] [PubMed]
51. Betancourt Ponce M, Boonpattrawong N, Fagan AJ, et al. Atovaquone-induced oxidative stress activates the pentose phosphate pathway and Immunogenic cell death in ovarian cancer. Sci Rep 2025;15:41702. [Crossref] [PubMed]
52. Luo J, Huang K, Yi X, et al. Bimetallic Nanozyme Amplifier for Synergistic Ferroptosis-Cuproptosis and Metabolic Reprogramming to Reshape Immunosuppressive Tumor Microenvironment. Adv Sci (Weinh) 2026;13:e12764. [Crossref] [PubMed]
53. Zhang Z, Tang S, Qi M, et al. Mitochondria-targeting natural product rhein conjugated with dichloroacetate as the dual inhibitor of glycolysis and oxidative phosphorylation to off energize cancer cells and induce ROS storm. Theranostics 2025;15:4909-29. [Crossref] [PubMed]
54. Tong J, Tan X, Hao S, et al. Inhibition of multiple CDKs potentiates colon cancer chemotherapy via p73-mediated DR5 induction. Oncogene 2023;42:869-80. [Crossref] [PubMed]
55. Orhan E, Velazquez C, Tabet I, et al. CDK inhibition results in pharmacologic BRCAness increasing sensitivity to olaparib in BRCA1-WT and olaparib resistant in Triple Negative Breast Cancer. Cancer Lett 2024;589:216820. [Crossref] [PubMed]
56. Houles T, Boucher J, Lavoie G, et al. The CDK12 inhibitor SR-4835 functions as a molecular glue that promotes cyclin K degradation in melanoma. Cell Death Discov 2023;9:459. [Crossref] [PubMed]
57. Tang Y, Zhang X, An L, et al. Diverse role of NMDA receptors for dendritic integration of neural dynamics. PLoS Comput Biol 2023;19:e1011019. [Crossref] [PubMed]
58. Bhattacharjee S, Iyer EPR, Iyer SC, et al. The Zinc-BED Transcription Factor Bedwarfed Promotes Proportional Dendritic Growth and Branching through Transcriptional and Translational Regulation in Drosophila. Int J Mol Sci 2023;24:6344. [Crossref] [PubMed]
59. Zhang Y, Zou J, Chen R. An M0 macrophage-related prognostic model for hepatocellular carcinoma. BMC Cancer 2022;22:791. [Crossref] [PubMed]
60. Feng Y, Geng Y, Liu Z, et al. QRICH1, as a key effector of endoplasmic reticulum stress, enhances HBV in promoting HMGB1 translocation and secretion in hepatocytes. Immunobiology 2025;230:152913. [Crossref] [PubMed]
61. Singharajkomron N, Seephan S, Iksen I, et al. CAMSAP3-mediated regulation of HMGB1 acetylation and subcellular localization in lung cancer cells: Implications for cell death modulation. Biochim Biophys Acta Gen Subj 2024;1868:130614. [Crossref] [PubMed]
62. Zhang M, Ma Z, Cui H, et al. Involvement of circadian clock protein PER2 in controlling sleep deprivation induced HMGB1 up-regulation by targeting p300 in the cortex. Sci Rep 2025;15:12253. [Crossref] [PubMed]
63. Cui S, Hu H, Chen A, et al. SIRT1 activation synergizes with FXR agonism in hepatoprotection via governing nucleocytoplasmic shuttling and degradation of FXR. Acta Pharm Sin B 2023;13:559-76. [Crossref] [PubMed]
64. Wang P, Wang L, Liu C, et al. YAP K236 acetylation facilitates its nucleic export and deprived the protection against cardiac hypertrophy in mice. Pharmacol Res 2025;211:107573. [Crossref] [PubMed]
65. Seo Y, Kim DK, Park J, et al. A Comprehensive Understanding of Post-Translational Modification of Sox2 via Acetylation and O-GlcNAcylation in Colorectal Cancer. Cancers (Basel) 2024;16:1035. [Crossref] [PubMed]
66. Montico B, Nigro A, Lamberti MJ, et al. Phospholipid scramblase 1 is involved in immunogenic cell death and contributes to dendritic cell-based vaccine efficiency to elicit antitumor immune response in vitro. Cytotherapy 2024;26:145-56. [Crossref] [PubMed]
67. Shoaib S, Ansari MA, Ghazwani M, et al. Prospective Epigenetic Actions of Organo-Sulfur Compounds against Cancer: Perspectives and Molecular Mechanisms. Cancers (Basel) 2023;15:697. [Crossref] [PubMed]
68. Manna S, Mishra J, Baral T, et al. Epigenetic signaling and crosstalk in regulation of gene expression and disease progression. Epigenomics 2023;15:723-40. [Crossref] [PubMed]