A total of 1466 genes (Supplementary Table 1) that are associated with twelve PCD procedures were collected from available gene sets, recent papers, and manual search [9]. Firstly, the PCD genes with prognostic value were identified in the Cancer Genome Atlas (TCGA)-LUAD cohort (prognosis cohort). The prognostic CDI was constructed and then validated in the Clinical Proteomic Tumor Analysis Consortium (CPTAC)-LUAD cohort (validation cohort). The correlation of CDI with clinical outcomes, functional enrichment, tumor microenvironment and drug sensitivity prediction was further investigated. The workflow diagram of this research is delineated in Fig. 1.

Fig. 1.

Overview of the main workflow of this study. LUAD, lung adenocarcinoma; TCGA, the Cancer Genome Atlas; CPTAC, the Clinical Proteomic Tumor Analysis Consortium; PCD, programmed cell death; DEGs, differentially expressed genes; GSVA, gene set variation analysis.

Raw transcriptome counts data from LUAD and normal lung tissues were subjected to comparison within the TCGA database. The DESeq2 package was used to identify differentially expressed genes (DEGs), with the criterion criteria set at adjusted p 0.05 and $|$log${}_2$fold change (FC)$|$ $>$1. The biological function enrichment of these DEGs was elucidated through the Database for Annotation, Visualization and Integrated Discovery (DAVID).

To pinpoint the PCD genes with significant implications for the prognosis of LUAD, a series of analyses including univariate Cox regression, Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression and multivariate Cox regression were executed. The ultimate prognostic CDI of each patient was calculated by the sum up the product of risk coefficient and the expression of each prognostic PCD gene. Subsequent exploration delved into the association of CDI with clinical status, encompassing survival, tumor stage, etc.

LUAD patients were stratified into high- and low-CDI groups based on the median CDI score of the TCGA-LUAD cohort. Gene set variation analysis (GSVA) was utilized to explore the altered biological signaling pathways among patients with high- and low-CDI through the R packages GSVA.

The independent prognostic efficacy of CDI score was appraised in TCGA-LUAD cohort and subsequently validated in CPTAC-LUAD cohort. The Kaplan-Meier curves depicting the high- and low-CDI groups of LUAD patients were constructed and compared. Employing the CDI score as the predictive variable, Receiver Operating Characteristic (ROC) curves were delineated at 1-, 3- and 5-year, with the area under the receiver operating characteristic curve (AUC) calculated using the timeROC R package. Calibration curves and dot plots were generated to generated the congruence of CDI-predicted survival and actual survival.

Based on the expression matrix of CDI genes, the ConsensusClusterPlus package was used to investigate the unidentified subtypes of LUAD. The parameter of consensus clustering was chosen as follows: maxK = 8, clusterAlg = hc and distance = pearson. The parameter k means that each sample is partitioned into up to k cluster by a specified clustering algorithm of agglomerative hierarchical clustering.

To verify the independent prognostic significance of CDI, univariate Cox regression and stepwise multivariate Cox regression were conducted, incorporating CDI score along with clinical parameters. To further developed a prognostic nomogram, the stepwise regression was conducted to select the final model factors. The predictive performance of the nomogram was evaluated using the Kaplan-Meier curve, AUC and calibration analysis as described above.

The mRNA levels of 70 immunomodulators which were involved in the functions including antigen presentation, cell adhesion, coinhibitor, costimulator, ligand, receptor, etc., were analysed to identify their correlation with CDI level. Then, we used CIBERSORT (Cell-type Identification by Estimating Relative Subsets of RNA Transcripts) to measure the infiltration of 22 types of immune cells in the LUAD tumors of the high- and low-CDI groups. CIBERSORT is a computational method for estimating immune cell composition in complex tissues base on their gene expression profiles.

Inhibitory concentration (IC50) of nearly two hundred common anti-tumor agents for each patient were predicted by using the oncoPredict package in R. To assess the potential response to immunotherapy, Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was employed between high- and low-CDI groups.

R software (v.4.3.0, https://www.r-project.org/) was used to perform all statistical analysis. Cumulative survival curve was presented by Kaplan-Meier plot and the differences of survival were evaluated using the log-rank test. p value less than 0.05 was considered as statistical significance.

In the prognosis cohort, a total of 236 DEGs encompassing 136 up-regulated genes and 100 down-regulated genes were delineated (Fig. 2A and Supplementary Table 2). The expression heatmap of these DEGs showed pronounced disparities between LUAD and normal lung tissues (Fig. 2B). To elucidate the primary functions associated with these DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed (Fig. 2C–F). The results indicated that these PCD DEGs are intricately linked to diverse functional pathways, including apoptotic process, inflammatory response, pathways in cancer and others.

Fig. 2.

Exploration of Differentially Expressed Genes (DEGs) in PCD in LUAD. (A) Volcano plot depicting Log fold change (FC) changes and p-values of DEGs in PCD relative to normal controls in LUAD. (B) Heatmap showing the expression levels of DEGs in PCD in both LUAD and normal controls. (C–F) Representation of functional enrichment pathways for these DEGs in terms of gene ontology-biological process (GO BP), cellular component (GO CC), molecular function (GO MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG), respectively.

Univariate COX regression analysis identified 49 PCD genes that were significantly linked to the survival of LUAD patients in the prognosis cohort with p value 0.05 (Supplementary Table 3). Using the LASSO COX regression, the collinearity between PCD genes as well as those with less pronounced survival significance were eliminated (Fig. 3A,B), and 15 genes were further selected. These 15 genes were then integrated into a multivariable COX regression model, and a stepwise selection analysis were performed. Consequently, the ultimate prognostic model comprised 10 PCD genes, namely ERO1A, CDK5R1, TRIM6, DNASE2B, ITPRIP, MRGPRX2, FGA, NDUFA13, NLRP2, and CD68 (Fig. 3C). The CDI for each patient can be computed based on the expression level of these 10 genes using the following formula: CDI = 0.2052 $\times$ERO1A + 0.2821 $\times$CDK5R1 + 0.3238 $\times$TRIM6 + (–0.1625) $\times$DNASE2B + (0.3568) $\times$ITPRIP + 0.3277 $\times$MRGPRX2 + 0.0672 $\times$FGA + (–0.2602) $\times$NDUFA13 + (–0.1133) $\times$NLRP2 + (–0.4214) $\times$CD68.

Fig. 3.

Establishment of cell death index (CDI). (A,B) Illustration of the Least Absolute Shrinkage and Selection Operator (LASSO) regression process for penalty coefficients, resulting in the inclusion of 10 genes for modeling. (C) Display of hazard ratio values, 95% Confidence Interval (CI), and p-values for the 10 genes in the CDI model. (D) Association between CDI scores in LUAD patients and clinical indicators. (E,F) Differentiation of high and low CDI groups in LUAD patients based on CDI levels, presented as heatmap and ridge plot, highlighting the top 10 significantly enriched KEGG pathways for each group. *, p 0.05; **, p 0.01; ***, p 0.001; ****, p 0.0001; ns, non-significant.

The relationship of CDI with clinical features were depicted in Fig. 3D. LUAD patients in the prognosis cohort who experienced mortality during follow-up demonstrated a notably higher CDI score. Furthermore, an observable inclination towards elevated CDI was noted in the context of more advanced tumor stages, encompassing American Joint Committee on Cancer (AJCC) stage, T stage, and N stage.

In furthering our comprehension of the distinctions in biological function status between high- and low-CDI groups, KEGG pathways analysis demonstrated that the cell cycle, P53 signaling pathway and several cancer pathways were notably prominent (Fig. 3E,F).

Employing the CDI as a prognostic factor, the predictive performance was evaluated in prognosis cohort and validation cohort, respectively. LUAD patients with high-CDI showed worse survival probability in both cohorts (Fig. 4A,B). The AUC values of CDI in the prognosis cohort at 1-, 3- and 5-year of follow-up were 0.784, 0.796 and 0.805, respectively, with the C-index value of 0.801 (Fig. 4C). The CDI also demonstrated commendable performance in the validation cohort with 1-, 3- and 5-year AUC values of 0.803, 0.787 and 0.778, respectively, as well as the C-index value of 0.794 (Fig. 4C). The calibration analysis of CDI in predicting the survival also exhibited highly favorable fitting with the actual survival outcomes in both LUAD cohorts (Fig. 4D,E).

Fig. 4.

Evaluation of the Prognostic Value of CDI. (A) Survival Kaplan-Meier (KM) curves for high and low CDI groups in LUAD patients from both TCGA and CPTAC datasets. (B) Correlation between individual patient survival risk and outcomes in high and low CDI groups from TCGA and CPTAC datasets. (C) Receiver Operating Characteristic (ROC) curves and area under the curve (AUC) values for CDI as a prognostic indicator in 1/3/5 years for LUAD patients in TCGA and CPTAC datasets. (D) Calibration curves for CDI as a prognostic indicator in 1/3/5 years for LUAD patients in TCGA and CPTAC datasets. (E) Calibration scatter plots for CDI as a prognostic indicator in 1/3/5 years for LUAD patients in TCGA and CPTAC datasets.

To investigate the unidentified subsets of LUAD, 10 PCD signature genes were utilized to conduct the consensus clustering procedure. According to the consensus index and Cumulative Distribution Function (CDF) curve, it was observed that the most significant differences among subsets were evident when k = 3 (Fig. 5A–C). Subsequently, there was a significant difference in survival among the subsets of LUAD patients (Fig. 5D). The distribution of survival and CDI risk groups revealed that cluster 3 was associated with a lower CDI and better survival outcome (Fig. 5E).

Fig. 5.

Unsupervised Clustering Analysis Based on CDI Genes. (A–C) Clustering analysis process demonstrating three distinct patient clusters. (D) Survival KM curves for the three patient clusters. (E) Scatter plot showing the distribution of CDI levels and survival outcomes for patients in the three clusters. CDF, Cumulative Distribution Function.

COX regression analysis showed that CDI is a significant risk factor for LUAD survival independent of the clinical factors including age, sex, stages of tumor (T), lymph node (N), metastases (M) and the AJCC stage (Fig. 6A,B). When adjusted for potential confounding variables which were non-significant in the multivariable COX model such as age, sex, etc., the hazard ratio for CDI was 2.655 (95% Confidence Interval (CI), 2.009–3.508, p 0.001). A prognostic nomogram model was established after multivariate COX regression and stepwise selection. AJCC stage, age and CDI were finally used for construction of the nomogram (Fig. 6C). The red lines, arrows and dots were plotted on the nomogram to show an example of how to use the nomogram. For a LUAD patient of 60 years old, stage II and with CDI score of 7.5, he gets a total risk score of 121 points which indicates that the 1-, 3- and 5-year death probability were 0.0619, 0.272 and 0.515, respectively. The Kaplan-Meier curves of the nomogram-predicted high- and low-points groups showed very obvious separation in both cohorts (Fig. 6D). The AUC values of 1-, 3- and 5-year ROC curves were all beyond 0.75, suggesting that the nomogram had a relatively high accuracy (Fig. 6E). Besides, the calibration curves also indicate favorable fitting of the predicted survival with actual survival outcomes (Fig. 6F).

Fig. 6.

Construction of Nomogram Based on CDI. (A,B) Forest plots displaying HR values, 95% CI, and p-values for CDI in single-factor Cox analysis and multi-factor Cox analysis, respectively. (C) Nomogram based on AJCC, age, and CDI, explaining how to predict a patient’s 1/3/5-year survival probability. (D) KM curves for high-risk and low-risk groups predicted by nomogram for TCGA and CPTAC datasets. (E) ROC curves for 1/3/5-year survival prediction based on nomogram for TCGA and CPTAC datasets. (F) Calibration curves for 1/3/5-year survival prediction based on nomogram for TCGA and CPTAC datasets. HR, hazard ratio.

The correlation analysis of CDI with the expression of immune-related genes suggested that higher CDI level might had lower antigen presentation, cell adhesion and stronger coinhibitor, indicating a relatively immuno-tolerant microenvironment in these LUAD tissues (Fig. 7A). The estimated proportion of immune cells in the tumor microenvironment was also calculated and shown in Fig. 7B. Elevated CDI score was correlated with marked increase in the infiltration of neutrophils, activated mast cells, activated dendritic cells, M0 macrophages, resting natural killer (NK) cells, $\gamma$$\delta$T cells, and activated memory CD4${}^{+}$T cells. Conversely, there was a reduction in the infiltration of resting mast cells, resting dendritic cells, monocytes, and activated NK cells (Fig. 7C,D).

Fig. 7.

Assessment of CDI-Related Immune Microenvironment. (A) Bar chart illustrating the correlation between CDI and various immune regulatory genes in LUAD patients from TCGA, with panel descriptions of each immune regulatory gene’s properties and type. (B) Heatmap displaying the proportion levels of 22 immune cells in each LUAD patient from TCGA, stratified into high and low CDI groups. (C) Heatmap further showing the correlation between 22 immune cells, CDI, and the 10 CDI modeling genes, with asterisks indicating statistically significant correlations (p 0.05). (D) Box plots summarizing the average levels and differences in 22 immune cells between high CDI and low CDI groups. *, p 0.05; **, p 0.01; ***, p 0.001; ****, p 0.0001; ns, non-significant.

To investigate the correlation between the CDI score and treatment sensitivity, we obtained IC50 values of each agent in LUAD tissues to identify those showing significant differences. The top 20 drugs with the most significant correlation of IC50 value with CDI score were recognized. As depicted in Fig. 8A, several commonly used chemotherapeutic drugs, such as cisplatin and targeting drugs, including osimertinib in NSCLC, were found to be highly associated with CDI score. Notably, patients with lower CDI scores appeared to be resistant to traditional chemotherapeutic drugs such as cisplatin, docetaxel and gemcitabine (with higher IC50), but sensitive to the recent developed targeting drug osimertinib (Fig. 8B).

Fig. 8.

Prediction of CDI-Related Drug Sensitivity. (A) Bubble plot demonstrating the correlation prediction between CDI score, CDI genes, and sensitivity to various anticancer drugs in LUAD patients from TCGA, with bubble colors indicating the level of correlation and bubble size representing statistical significance. (B) Box plots presenting IC50 values for common chemotherapy and targeted drugs in high CDI and low CDI groups. (C) Box plots showing TIDE score, Exclusion score, and Dysfunction score in high CDI and low CDI groups. (D) Scatter plots describing the correlation between CDI score and Tumor Immune Dysfunction and Exclusion (TIDE) score, Exclusion score, and Dysfunction score. **, p 0.01; ****, p 0.0001; ns, non-significant.

Regarding immunotherapy prediction (Fig. 8C,D), CDI showed no significant difference in relation to the TIDE score. However, patients with higher CDI possessed elevated exclusion score and dysfunction score, suggesting a positive association between CDI and T cell immune exclusion, and a negative association with T cell function.