Introduction

Malignant lung tumors are the most common malignant tumors of the respiratory system and the leading cause of cancer-related deaths in both men and women.1 The WHO classifies malignant lung tumors based on histological type into non-small cell lung cancer (NSCLC) and small cell lung cancer (SCLC),2 with the former accounting for approximately 85% and the latter 15%. NSCLC primarily includes adenocarcinoma and squamous cell carcinoma. Currently, lung adenocarcinoma (LUAD) is the most common type of NSCLC, comprising about 40% of lung cancer cases.3 Although surgery is the most effective treatment for stage I to II and some stage IIIA LUAD, the invasive nature of these tumors often makes complete surgical resection challenging. The advent of immunotherapy has ushered in a new era for LUAD patients, with immune checkpoint inhibitors significantly improving clinical outcomes and offering a neoadjuvant immunotherapy option for early resectable disease. Unfortunately, research indicates that while immunotherapy significantly improves the prognosis for some patients, the 5-year survival rate remains only 4%–10%.4 Therefore, there is an urgent need to elucidate the underlying molecular mechanisms and develop reliable molecular classification models to effectively assess prognosis and guide personalized treatment strategies for LUAD patients.

Programmed cell death (PCD), also known as regulated cell death, is a crucial physiological process that plays a pivotal role in maintaining tissue homeostasis and eliminating damaged or unnecessary cells. PCD can occur through various mechanisms, including apoptosis, anoikis, autophagy, alkaliptosis, cuproptosis, entosis, entotic cell death, immunogenic cell death, ferroptosis, lysosome-dependent cell death, methuosis, necroptosis, netotic cell death, NETosis, oxeiptosis, pyroptosis, parthanatos, and paraptosis. PCD involves a regulated program of intrinsic clearance mechanisms to remove damaged or unnecessary cells, thereby maintaining tissue health.5 Among these PCD mechanisms, several are associated with mitochondrial dysfunction. Mitochondria play a key role in providing energy for cellular functions, regulating cell signaling pathways, and controlling PCD. Recent studies suggest that cancer is primarily a mitochondrial metabolic disease, with mitochondrial dysfunction characterized by changes in structure, function, and dynamics being implicated in tumorigenesis.6 Apoptosis is a well-recognized PCD mechanism that plays a significant role in maintaining tissue homeostasis and eliminating damaged or unnecessary cells, characterized by a series of biochemical and morphological changes.4 Inhibition or resistance to apoptosis often contributes to cancer development. Pyroptosis is a form of PCD that occurs following inflammasome activation and caspase-1 cleavage, characterized by cell swelling, membrane rupture, and the release of pro-inflammatory cytokines. Ferroptosis is a recently discovered PCD characterized by iron-dependent cell death and lipid peroxidation. Autophagy is a cellular mechanism that maintains cellular homeostasis by degrading damaged proteins and organelles. Autophagy can promote cell survival or induce cell death, depending on the specific context in which it occurs. Previous studies have shown that USP15 may induce autophagy to negatively regulate lung cancer progression. Necroptosis is a form of PCD characterized by necrosis-like cell death, triggered by the activation of RIPK1 and RIPK3. Cuproptosis is a PCD induced by copper overload, characterized by lipid peroxidation and mitochondrial dysfunction. Entotic cell death occurs in viable cells and their neighboring regions, distinct from traditional apoptotic pathways, as it does not require the activation of apoptotic execution pathways. Netotic cell death is another form of PCD resulting from the release of neutrophil extracellular traps, commonly observed during infections or injuries. Parthanatos is a tightly regulated form of cell death induced by the overactivation of the nuclear enzyme PARP-1. Lysosome-mediated cell death involves the action of hydrolytic enzymes that permeate the cell membrane and enter the cytoplasm. Currently, ferroptosis, cuproptosis, and pyroptosis have garnered widespread discussion in the context of LUAD.7–9 Additionally, alkaliptosis is an emerging form of PDC regulated by intracellular alkalinization processes. Oxeiptosis, leveraging the reactive oxygen species sensing ability of KEAP1, is a recently identified cellular pathway that likely interacts with other cell death pathways.

In the context of LUAD, tumor cells can evade the activation of PCD mechanisms through various strategies. They may activate specific signaling pathways or alter their morphology to avoid being eradicated by the host. As our understanding of PCD mechanisms deepens, numerous drugs targeting these pathways have been developed. For instance, a novel HDAC1/2 inhibitor can effectively inhibit colorectal cancer by regulating the cell cycle and inducing apoptosis.10 GSDME-induced pyroptosis is a unique form of PCD that has shown potential as an antitumor immunotherapy. Additionally, research has indicated that obstructing ferroptosis can lead to resistance against programmed cell death protein-1 (PD-1)/programmed death-ligand 1 (PD-L1) therapy. These findings underscore the significance of PCD research in enhancing our understanding of LUAD and developing new anti-LUAD therapies.

Disruption of mitochondrial morphology, such as changes in shape, size, or cristae organization, can impair normal mitochondrial function and trigger PCD. Structural abnormalities may affect the release of pro-apoptotic factors from mitochondria, leading to caspase activation and subsequent apoptosis. Mitochondrial function is also closely associated with PCD mechanisms. Dysfunctional mitochondria with impaired oxidative phosphorylation and ATP production can induce cellular stress and initiate PCD pathways. To better link mitochondrial function with the aforementioned PCD mechanisms, this study screened a series of relevant biomarkers. Currently, several molecular markers have been identified as clinically significant for the diagnosis and prognosis of LUAD.11 Markers such as PD-L1 expression levels, tumor mutational burden (TMB), and hematological biomarkers play a crucial role in predicting treatment efficacy and are closely related to the prognosis of LUAD. Besides genetic factors, various environmental and lifestyle factors contribute to the development of LUAD. Carcinogens in tobacco smoke can induce differentiation of bronchial epithelial cells, thereby promoting tumor formation.12

In the process of screening and validating myriad biomarkers, various confounding factors—including sample selection bias, intratumoral and intertumoral heterogeneity, analytical variability, and publication bias—can significantly impact the accuracy and generalizability of results. Consequently, the identification of survival-associated genes through transcriptomic profiling is paramount for prognostic stratification and the selection of targeted therapeutic modalities. Over recent decades, a substantial body of evidence has elucidated the pivotal role of mitochondrial dysfunction and PCD mechanisms in the pathogenesis and metastatic progression of malignant neoplasms. Neoplastic cells must circumvent diverse modalities of cellular demise and mitigate mitochondrial dysfunction to sustain proliferation and invasion. Notwithstanding these advances, a comprehensive elucidation of the intricate interplay between mitochondrial perturbations and PCD pathways in LUAD remains elusive. Moreover, there is a paucity of detailed functional studies delineating these processes specifically within the context of LUAD tumorigenesis and progression. To bridge these knowledge gaps, we introduced a mitochondrial-related PCD signature (MPCDS) to predict the efficacy of therapeutic interventions and prognosis in LUAD. Through our research, we have uncovered the high heterogeneity13 among LUAD patients and evaluated their clinical outlooks.

MethodDataset source

The transcriptomic, copy number variation (CNV), mutation, H&E slides, and clinical data of LUAD were successfully retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) and used as the training set for model construction. Six transcriptome datasets were obtained from the GEO database for model validation, including GSE13213 dataset14 (n=119), GSE26939 dataset15 (n=115), GSE29016 dataset16 (n=39), GSE30219 dataset17 (n=86), GSE31210 dataset18 (n=227), and GSE42127 dataset19 (n=134). Three LUAD immunotherapy cohorts (OAK dataset,20 POPLAR dataset,21 NG dataset22) and three pan-cancer immunotherapy cohorts (GSE91061 dataset,23 phs000452 dataset,24 and Braun_2020 dataset25) were used to validate the model’s ability to predict the prognosis of immunotherapy patients. A pan-cancer cohort covering 21 tumors was used to evaluate the model’s broad applicability (data sourced from UCSC, https://xena.ucsc.edu/). Two single-cell datasets were sourced from GEO (GSE189357 dataset26) and the Genome Sequence Archive at the BIG Data Center (HRA001130 dataset27). All datasets used in this study are sourced from those listed in online supplemental table 1. Mitochondrial and PCD genes were derived from previous studies28; see online supplemental table 2 for details.

To ensure data uniformity and comparability, gene expression data were converted to transcripts per million format. The “combat” function in the “sva” package29 was applied to mitigate potential batch effects. Additionally, all datasets obtained from TCGA and GEO underwent log transformation to establish a standardized data format from the onset of the analysis. Principal component analysis (PCA) was used to evaluate batch effects between datasets.

Identification of prognostic MPCDS

The limma package was used to analyze differentially expressed mitochondrial PCD-related genes (MPRGs) between normal and tumor samples, using selection criteria of false discovery rate (FDR) <0.05 and log2 fold change (FC)>1. Univariate Cox regression analysis was conducted to evaluate genes with prognostic value for LUAD patients. Subsequently, using 10-fold cross-validation, we applied 101 combinations of 10 machine learning algorithms, including stepwise Cox, Lasso, Ridge, partial least squares regression for Cox (plsRcox), CoxBoost, random survival forest (RSF), generalized boosted regression modeling (GBM), elastic net (Enet), supervised principal components (SuperPC), and survival support vector machine (survival-SVM). The goal was to identify the most valuable MPCDS characterized by the highest concordance index (C-index).

Our machine learning approach is similar to the combination of 10 machine learning algorithms used in Liu et al’s article.30 We used the same algorithms: stepwise Cox, Lasso, Ridge, plsRcox, CoxBoost, RSF, GBM, Enet, SuperPC, and survival-SVM. Both studies employed 101 algorithm combinations and used the C-index as the metric for evaluating model performance. The main difference lies in the cross-validation method: we used 10-fold cross-validation, while Liu et al used leave-one-out cross-validation.

The accuracy of MPCDS was evaluated using receiver operating characteristic (ROC) curves and PCA. Subsequently, we comprehensively retrieved 114 prognostic signatures from published papers, including lncRNA and mRNA, and assessed the performance of MPCDS using the C-index as the evaluation metric.

Investigating potential mechanisms and pathways

To investigate the biological functions and pathway processes associated with MPRGs, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the “clusterProfiler” package. Specifically, we used MPRGs as inputs, converted them to Entrez IDs, and conducted GO and KEGG enrichment analyses with an adjusted p value <0.05 as the criterion. Additionally, the potential mechanisms were investigated using Gene Set Variation Analysis (GSVA). The cancer immunity cycle and pathways predicted to respond to immunotherapy were examined according to previously established methods.31 32 A comprehensive collection of pathway gene sets was sourced from the Molecular Signatures Database.

Comprehensive analysis of genomic alterations and immune cell infiltration

Genomic alterations, including recurrently amplified and deleted regions, were identified using GISTIC 2.0 analysis. The TMB was calculated using the R package “maftools.” Previous studies have evaluated immune cell infiltration abundance in TCGA data, with IFNG scores sourced from the TIDE website (http://tide.dfci.harvard.edu/). The Cancer Immunome Atlas (https://tcia.at/home) assesses the immunophenoscore (IPS) of LUAD patients and identifies populations suitable for immunotherapy. Additionally, the ssGSEA algorithm evaluates the abundance of immune cell infiltration and the activity of immune-related pathways in tumor samples. The TIMER2.0 database provides a summary of immune cell infiltration abundance in TCGA, incorporating results from multiple algorithms (see online supplemental table 3 for details). The IOBR package33 includes a variety of immune infiltration algorithms that are used to infer the abundance of immune cells within tumor samples.

Preprocessing and integration of single-cell gene expression data

The initial gene expression matrix underwent preprocessing using the Seurat R package (V.4.2.0).34 Genes expressed in at least 10 cells within each sample were included in the analysis. Subsequently, low-quality cells were excluded based on specific criteria: cells with more than 5000 or fewer than 200 expressed genes, or cells with more than 10% of unique molecular identifiers originating from the mitochondrial genome. Sample integration was performed using the harmony R package. Next, a set of highly variable genes was selected for PCA, and the top 30 significant principal components were chosen for Uniform Manifold Approximation and Projection dimension reduction to visualize gene expression patterns. Differentially expressed genes within each cell subpopulation were identified using the “FindAllMarker” function and annotated based on previously published reliable cell type markers.27

Cell–cell communication

We used CellChat35 to integrate gene expression data and analyze differences in cell–cell communication modules. The default CellChatDB was used as the ligand-receptor database, following the standard CellChat protocol. By detecting overexpressed ligands or receptors within cell groups, we inferred cell type-specific interactions and identified enhanced ligand-receptor interactions when ligands or receptors were highly expressed.

Analysis of drug sensitivity in cancer cell lines and LUAD patients

Expression profile data of human cancer cell lines were sourced from the Broad Institute Cancer Cell Line Encyclopedia database (https://sites.broadinstitute.org/ccle). CERES scores for 18333 genes in 739 cell lines were obtained from the DepMap portal (https://depmap.org/portal/). CERES scores measure the dependency of genes in specific cell lines, with lower scores indicating a higher likelihood that the gene is essential for cell growth and survival. Drug sensitivity data for cancer cell lines were obtained from the Cancer Therapeutics Response Portal (CTRP,36 https://portals.broadinstitute.org/ctrp.v2.1/) and the PRISM database37 (https://www.theprismlab.org/). These datasets provide the area under the dose-response curve (AUC) as a measure of drug sensitivity, with lower AUC values indicating greater sensitivity. Missing AUC values were imputed using the K-nearest neighbor method, and compounds with more than 20% missing data, as well as hematopoietic and lymphoid tissue cell lines, were excluded.38 The sensitivity of each LUAD patient to chemotherapeutic agents was assessed based on the IC50 values quantified using the pRRophetic package from the Genomics of Drug Sensitivity in Cancer database.39 A ridge regression model was applied to generate drug sensitivity estimates for each patient.

Cell culture and small interfering RNA transfection

The A549 and H1299 LUAD cell lines were obtained from the Institute of Biochemistry and Cell Biology at the Chinese Academy of Sciences in Shanghai, China. The culture medium, RPMI 1640, was supplemented with 10% FBS and 1% antibiotics (100 U/mL penicillin and 100 mg/mL streptomycin) and was used to maintain the A549 and H1299 cells, respectively. Small interfering RNA (siRNA) transfection was performed using the Lipo2000 reagent (Invitrogen, Shanghai, China) strictly following the manufacturer’s protocol. Typically, A549 and H1299 cells were seeded onto coverslips in six-well plates, and siRNA transfection was carried out the following day. qRT-PCR was used to assess knockdown efficiency. The primer information for nucleoside diphosphate kinase 4 (NME4) and GAPDH is provided in online supplemental table 4.

Colony formation

5000 cells were seeded into each well of a six-well plate for the colony formation assay, with conventional growth medium added and replaced after 1 week. After 2 weeks, once the colonies had matured, they were fixed with methanol for 15 min and then stained with 0.1% crystal violet (Sigma) for 30 min. Following staining, the colonies were counted to assess their formation capability.

Invasion and migration assays

Invasion and migration assays were conducted using Corning’s Transwell system (24 wells, 8 µm pore size, New York, USA). For migration assays, 5×104 transfected cells were placed in the upper chambers with 350 µL of serum-free medium, while 700 µL of medium containing 10% FBS was added to the lower chambers. The invasion assays used Transwell membranes precoated with Matrigel (Sigma-Aldrich). After a 16-hour incubation, cells on the upper surface were removed, and those that migrated to the lower surface were stained with methanol and 0.1% crystal violet. Images were captured using an Olympus inverted microscope (Tokyo, Japan).

Clinical and pathological characteristics of LUAD patients

All paraffin-embedded tissue sections were obtained from the Department of Pathology at Tianjin Medical University Cancer Institute and Hospital. All case slides were histopathologically confirmed as LUAD, and none of the patients received radiotherapy, chemotherapy, or targeted therapy prior to surgery. The clinical and pathological characteristics of the patients are shown in online supplemental table 5.

Multiplex immunohistochemistry protocol for NME4 and immune cell marker detection

According to the standard protocol,40 multiplex immunohistochemistry (mIHC) was used to examine the expression relationship between NME4 and immune cell markers, with simultaneous detection of DAPI. The primary antibodies used for mIHC included: NME4 (1:200 dilution, Cat# ab228005, Abcam, USA), CD3 (1:150 dilution, Cat# ab21703, Abcam, USA), and CD20 (1:100 dilution, Cat# ab64088, Abcam, USA). The immunohistochemistry protocol involved deparaffinizing the LUAD tissue sections using xylene, followed by rehydration in a graded series of ethanol solutions. To inhibit peroxidase activity, the sections were treated with 3% hydrogen peroxide and then blocked with 5% normal goat serum. Subsequently, the sections were incubated overnight at 4°C with primary antibodies, including NME4 (1:200 dilution, Cat# ab228005, Abcam, USA) and CD8A (1:2000 dilution, Cat# ab217344, Abcam, USA). Following this, HRP-conjugated secondary antibodies were applied and incubated for 30 min at the same temperature. The sections were then stained with DAB (3,3′-diaminobenzidine) and counterstained with hematoxylin for visualization. Two independent pathologists independently evaluated all stained sections. For the semi-quantitative evaluation of NME4 staining, the H-score criterion was used. Additionally, tumors were classified into three phenotypes based on the spatial distribution of CD8+ T cells: inflamed, excluded, and desert subtypes.

Subcutaneous xenograft experiments in C57BL/6 mice

C57BL/6 mice were purchased from Jiangsu Gempharmatech (Jiangsu, China). For subcutaneous xenograft experiments, 2×106 LLC cells (sh-NME4) were resuspended in 100 µL PBS and injected subcutaneously into the mice. Once the tumors reached approximately 100 mm3, the mice were treated with either mouse PD-1 monoclonal antibody (mAb) (BioXcell, BE0146; New Haven, CT, USA) or IgG isotype control (BioXcell, BE0089) via intraperitoneal injection (100 mg per mouse in 100 µL D-PBS buffer). Tumor size was measured every 2–3 days. Tumor weight was recorded, and volumes were estimated using the formula 1/2×(length×width2).

Flow cytometry analysis of tumor-infiltrating immune cells in C57BL/6 mice

After euthanizing the C57BL/6 mice, fresh tumor tissue was excised and digested with collagenase IV (Gibco 17104019) and DNase I (Roche 10104159001) for 2 hours. The digested tissue was then filtered through 70 µm cell strainers to obtain a single-cell suspension. The cells were washed, resuspended in 40% Percoll, and centrifuged against 70% Percoll to isolate lymphocytes. For intracellular cytokine staining, the cells were stimulated for 4 hours with a leukocyte activation cocktail (BD Pharmingen, 550583; San Diego, CA, USA). Non-specific antibody binding was blocked with CD16/CD32 antibody (553141), followed by surface staining with anti-CD45 (553080), anti-CD3e (553064), and anti-CD8α (551162) antibodies (BD Pharmingen). After fixation and permeabilization with a Fixation/Permeabilization kit (BD Biosciences, 554714), intracellular GZMB was stained with an anti-GZMB antibody (Biolegend, 372204; San Diego, CA, USA). Stained cells were analyzed using a BD FACS Canto II flow cytometer, and the data were processed with FlowJo software.

Statistical analysis

Data manipulation, statistical analyses, and visualization were performed using R software V.4.2.0. Kaplan-Meier analysis and log-rank tests were used to estimate and compare subtype-specific overall survival (OS). Differences in continuous variables between two groups were assessed using the Wilcoxon test or t-test, while categorical variables were evaluated with χ2 tests or Fisher’s exact tests. P values were adjusted using the FDR method. Pearson correlation analysis was employed to explore relationships among variables. All p values were calculated using a two-tailed test, with p<0.05 considered statistically significant.

ResultElucidating the role of mitochondrial and PCD-related biomarkers

Tumor cells can evade PCD mechanisms through various strategies, affecting the tumor microenvironment and thereby limiting immune response, which has become a significant obstacle in current cancer treatments. Numerous studies have demonstrated the close connection between mitochondrial function and PCD. To better elucidate this relationship, our study screened a series of mitochondrial and PCD-related biomarkers specifically for LUAD. Figure 1 presents the overall framework of this study. To ensure the consistency of data across different databases, we removed batch effects from normal samples in both the GTEx and TCGA databases (online supplemental figure S1A). We integrated normal lung tissue data from GTEx and LUAD samples from TCGA to perform differential expression analysis and identify differentially expressed MPRGs (figure 2A). Next, we conducted univariate Cox analysis to identify MPRGs associated with prognosis (figure 2B), ultimately identifying 92 prognosis-related MPRGs (p<0.05, online supplemental table 6). Further GO and KEGG enrichment analyses revealed that MPRGs are mainly enriched in mitochondrial components, regulation of autophagy, apoptosis, PI3K-Akt, other signaling pathways, and so on (figure 2C,D). The chromosomal locations and expressions of MPRGs are shown in figure 2E. To ensure comparability between datasets, we removed batch effects from seven LUAD transcriptome cohorts (figure 2F). CNV status analysis indicated frequent changes in MPRGs. Notably, STYXL1 and STAT1 exhibited the most widespread CNV amplification (figure 2G).

Figure 1

Workflow summary. This figure illustrates the workflow of the study, exploring the roles of PCD and mitochondrial function in LUAD. The study used machine learning to construct the MPCDS for predicting prognosis and immunotherapy response in LUAD patients. NME4 was identified as a substitute marker for MPCDS, and its role was demonstrated through in-house cohorts and in vivo/vitro experiments. PCD, programmed cell death; LUAD, lung adenocarcinoma; MPCDS, mitochondrial-related PCD signature; NME4, nucleoside diphosphate kinase 4; PD-1, programmed cell death protein-1; ROC, receiver operating characteristic; TME, tumor microenvironment.

Figure 2

Identification and analysis of prognostic MPRGs. (A) Identification of differentially expressed MPRGs in tumor compared with normal lung tissue. (B) Univariate Cox analysis to determine potential prognostic MPRGs. (C and D) GO and KEGG enrichment analyses revealing the potential pathways enriched by prognostic MPRGs. (E) Chromosomal locations of prognostic MPRGs, with red dots representing higher expression in tumors relative to normal lung tissue (depicted by log2FC values; the larger the value, the more peripheral the position), and black dots representing the opposite. The central pie chart shows all LUAD datasets used for modeling. (F) Principal component analysis plot after batch effect removal in seven LUAD datasets. (G) Chromosomal locations and copy number variations of prognostic MPRGs. MPRGs, mitochondrial programmed cell death related genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LUAD, lung adenocarcinoma.

Development and validation of an MPCDS for prognosis and immunotherapy suitability in patients with cancer

Using the expression profiles of 92 prognosis-related MPRGs, we developed an MPCDS through a machine learning combinatorial algorithm. The TCGA dataset served as the training cohort, while six GEO datasets were used for validation. The average C-index across the six validation cohorts was the criterion for model selection. Ultimately, the Lasso+Stepcox (forward) algorithm emerged as the optimal model (figure 3A). The MPCDS score distinguished patient prognosis across all seven cohorts (figure 3B–H). Patients in the high MPCDS group exhibited worse outcomes compared with those in the low MPCDS group. Interestingly, we extrapolated the MPCDS scores for the immunotherapy cohorts using the model’s formula and found consistent results in three LUAD immunotherapy cohorts (figure 3I–L) and three pan-cancer cohorts (including two melanoma cohorts and one clear cell renal cell carcinoma cohort, figure 3M–P). The low MPCDS group appeared to be more suitable for immunotherapy, as they demonstrated better prognoses following treatment (figure 3I–P). Furthermore, IPSs calculated using the TCIA website also suggested that the low MPCDS group is more appropriate for treatment with PD-1, CTLA-4, or their combination therapy (figure 3Q–T).

Figure 3

Construction and validation of MPCDS using machine learning. (A) Construction of MPCDS using various machine learning combinations, with values in the heatmap representing the C-index of corresponding models for predicting prognosis; the bar graph on the right shows the average C-index across multiple datasets. (B–H) MPCDS effectively stratifies prognosis in LUAD patients across seven datasets, with low MPCDS group showing better outcomes. (I–P) Evaluation of MPCDS in immunotherapy cohorts, including three LUAD immunotherapy cohorts and three pan-cancer cohorts (two melanoma and one clear cell renal cell carcinoma cohorts), showing consistent stratification of prognosis. (Q-T) Prediction of IPS for TCGA-LUAD patients using The Cancer Immunome Atlas, consistently indicating higher IPS and greater sensitivity to immunotherapy in the low MPCDS group. MPCDS, mitochondrial-related programmed cell death signature; C-index, concordance index; LUAD, lung adenocarcinoma; IPS, immunophenoscore; TCGA, The Cancer Genome Atlas.

Superior predictive performance of MPCDS in prognosis

To evaluate the predictive efficacy of MPCDS, we integrated clinical features from seven datasets. MPCDS demonstrated higher C-index values than any other clinical feature (such as age, gender, stage, EGFR status, etc) (figure 4A). Subsequently, we conducted a PCA based on the expression levels of the model genes of MPCDS across all datasets. The results showed that the expression levels of the model genes could effectively distinguish LUAD patients, with high and low MPCDS groups forming distinct clusters (figure 4B). ROC curves further confirmed that MPCDS scores could reliably predict the prognosis of LUAD patients, which was validated across all datasets (with 1-, 3-, and 5-year AUC values generally above 0.65, figure 4C). Next, we compared MPCDS against 114 previously published LUAD signatures and found that MPCDS consistently demonstrated the best predictive performance across all datasets, achieving the highest C-index values (figure 4D). Further exploration of the correlation between MPCDS and clinical indicators revealed that patients with high MPCDS scores had a higher number of deaths, and tended to have more advanced T stage, N stage, and overall stage (online supplemental figure S1B). Additionally, the model genes LDHD, CTSG, and ADRB2 were highly expressed in the low MPCDS group, while the remaining model genes were highly expressed in the high MPCDS group. These observations indicate that MPCDS is an effective prognostic predictor and is significantly associated with clinical characteristics.

Figure 4

Comparison of MPCDS with other clinical indicators for prognostic prediction. (A) Comparison of MPCDS with other clinical indicators for prognostic prediction. (B) Principal component analysis based on model gene expression in seven datasets, showing clear separation into two clusters (red dots represent high MPCDS group, green dots represent low MPCDS group). (C) ROC curves evaluating the prognostic prediction capability in seven datasets. (D) Comparison of MPCDS with publicly published LUAD-related signatures, demonstrating superior prognostic capability of MPCDS. MPCDS, mitochondrial-related programmed cell death signature; ROC, receiver operating characteristic; LUAD, lung adenocarcinoma.

Complex relationship between MPCDS and key cancer immunity and biological pathways

To explore the potential pathway mechanisms underlying the differences between different MPCDS groups, GSVA enrichment analysis was conducted. The primary pathways enriched in the high MPCDS group include DNA replication checkpoint signaling, regulation of mitotic cytokinesis, positive regulation of DNA-templated DNA replication, DNA unwinding involved in DNA replication, phosphorylation of EMI1, and MYC gene upregulation-related target genes (online supplemental figure S2A–G). The abnormal activation or suppression of these pathways leads to uncontrolled cell cycle progression, increased DNA replication stress, accumulation of genetic variations, and significant chromosomal structural changes. These collectively promote the unlimited proliferation of tumor cells, genomic instability, and enhanced invasive and metastatic potential. Next, we explored the associations between MPCDS scores and immune therapy-related pathways as well as cancer immunity cycle steps. First, MPCDS showed significant positive correlations with various cancer immunity cycle steps (online supplemental figure S2H), such as cancer antigen release (step 1) and the recruitment of immunosuppressive cells (eg, MDSCs and neutrophils). This suggests that high MPCDS score groups may exhibit heightened activity in these immune steps, leading to more effective antigen release and recruitment of immunosuppressive cells. Second, MPCDS demonstrated strong correlations with several key biological pathways, such as DNA replication, cell cycle regulation, and mismatch repair. These pathways may be abnormally activated or suppressed in high MPCDS score groups, thereby affecting tumor progression and immune response. For example, abnormally active DNA replication and cell cycle regulation can result in rapid tumor cell proliferation and genomic instability, while alterations in mismatch repair mechanisms may increase mutation rates. Surprisingly, when we extended the analysis of MPCDS to a pan-cancer level, GSVA confirmed our previous results. MPCDS scores were significantly positively correlated with cell cycle-related pathways such as E2F targets, MYC targets V1/V2, and the G2M checkpoint across 21 types of tumors (online supplemental figure S5A,B). Overall, patients with high MPCDS scores exhibit more pronounced activity in these key pathways and immunity cycle steps, which are closely related to tumor invasiveness and prognosis. This indicates a complex relationship between MPCDS and multiple cancer-related biological pathways and immunity cycle steps, where high MPCDS scores may predict poorer prognosis and higher tumor aggressiveness.

Exploring the relationship between tumor mutations and patient prognosis

Given the significant correlation between tumor mutations, immune response, and patient prognosis, this analysis explores these relationships in detail. Online supplemental figure S3A,B vividly shows distinct chromosomal alterations between high- and low-MPCDS groups. The heatmap highlights a markedly elevated TMB in the high-MPCDS group, with TP53, TNN, CSMD3, RYR2, and LRP1B being the most frequently mutated genes (online supplemental figure S3C). Analyses reveal more frequent chromosomal deletions or amplifications within the high-MPCDS group (online supplemental figure S3D). Furthermore, online supplemental figure S3E indicates the worst prognosis in the L-TMB+ high-MPCDS subgroup.

Immune infiltration and regulatory gene expression in high and low MPCDS groups

Genomic-level comparisons revealed that the proportions of leukocytes, lymphocytes, and NK cells were higher in the low MPCDS group than in the high MPCDS group (p<0.05) (figure 5A–D). Additionally, we used the tumor-infiltrating lymphocyte (TIL) fraction data provided by Saltz et al41, who employed deep learning methods to estimate TILs on H&E-stained slides. Strikingly, the H&E-estimated TIL fraction was higher in the low MPCDS group (p<0.05) (figure 5E). The IFNG levels estimated by the TIDE website were also higher in the low MPCDS group (figure 5F). Immune regulatory genes play a significant role in tumor immune response, and we further compared the relative expression levels of immune regulatory genes between the two groups. The low MPCDS group appeared to favor the activation of MHC class II molecules and co-stimulatory molecules, whereas the high MPCDS group showed activation of MHC class I molecules and certain co-inhibitory molecules(figure 5G). Furthermore, we used the ssGSEA algorithm to calculate the abundance of immune infiltrating cells and immune-related pathway activity between the high and low MPCDS groups. The results were consistent with previous findings: the low MPCDS group had higher immune cell abundance, such as B cells and helper T cells, and the activation of type II interferon response and antigen presentation-related immune pathways (figure 5H). The resultsN calculated using the TIMER2.0 database further corroborated these findings (figure 5I). To confirm the above results, we examined H&E slides from the TCGA database for the high and low MPCDS groups and found more lymphocyte infiltration in the low MPCDS group (figure 5J).

Figure 5

Immune cell infiltration and genomic differences between high and low MPCDS groups. (A–E) Genomic and H&E staining evaluations of immune cell infiltration differences between different MPCDS groups. (F) TIDE website prediction of IFNG differences between MPCDS groups. (G) Expression analysis of immune regulatory genes between the two groups. (H) ssGSEA evaluation of immune cell infiltration and immune-related pathway activity between the groups. (I) TIMER2.0 database assessment of immune cell abundance differences between high and low MPCDS groups. (J) Representative slices of TCGA-LUAD high and low MPCDS groups. MPCDS, mitochondrial-related programmed cell death signature; TCGA, The Cancer Genome Atlas; LUAD, lung adenocarcinoma.

Detailed exploration of MPCDS distribution and cellular interactions in different cell types

We explored the detailed distribution of MPCDS across different cell types using single-cell RNA transcriptome data (GSE189357 and HRA001130). We annotated the major cell types in these datasets (figure 6A,D) and examined the distribution and intensity of MPCDS scores across different cell populations (figure 6B,E). Based on the median MPCDS score, we divided the cell populations into high and low MPCDS groups to observe cell type proportions (figure 6C,F). We found that the proportion of cycling cells was higher in the high MPCDS group, while immune cells such as B cells, plasma cells, and NK/T cells were lower. This finding is consistent with the bulk transcriptome analysis results. Violin plots also demonstrated that cycling cell types had the highest MPCDS scores (figure 6G,H), indicating that higher MPCDS scores may be associated with greater proliferative and stem cell potential. Furthermore, we examined cell interactions between high and low MPCDS groups (figure 6I–K) and found stronger cell–cell communication and more robust outgoing and incoming signals in the high MPCDS group (figure 6L,M). The high MPCDS group is likely to have stronger pro-tumor signaling pathways, including VEGF, EGF, and PDGF pathways (online supplemental figure S4).

Figure 6

Single-cell analysis of MPCDS in different cell types. (A–F) Annotation of cell types after dimensionality reduction clustering using two single-cell datasets, observing MPCDS distribution across different cell types. Median values used to distinguish high and low MPCDS groups, with differences in cell type proportions observed between the groups. (G and H) MPCDS differences across different cell types, with cycling cells consistently having the highest MPCDS scores. (I–K) Analysis of cell communication differences between different MPCDS groups. (L and M) Comparison of outgoing and incoming signal strength of different cell types between high and low MPCDS groups. MPCDS, mitochondrial-related programmed cell death signature; UMAP, Uniform Manifold Approximation and Projection.

Identification of therapeutic targets and candidate drugs for high MPCDS score patients

Proteins that show a high positive correlation with MPCDS may have potential therapeutic implications for patients with high MPCDS scores. However, most human proteins remain undruggable because they lack clear active sites for small molecule compounds to bind, or they are located within cells, making them inaccessible to biological agents. Therefore, to identify potentially druggable therapeutic targets for patients with poor prognosis in the high MPCDS group, we calculated the correlation coefficients between the expression levels of druggable proteins and MPCDS, identifying protein targets with correlation coefficients greater than 0.30 (p<0.05). Subsequently, using the CERES scores and MPCDS scores from LUAD cell lines, we identified genes with Spearman’s r<−0.3 and p<0.05. Through these two analyses, we identified seven genes: PARP1, HK2, P4HTM, PLOD2, LDHA, GPX8, and RFK (online supplemental figure S6A–O). Notably, the CERES scores of PLOD2 were greater than zero in most LUAD cell lines, suggesting that PLOD2 may not be essential in LUAD. The remaining six genes are considered potential therapeutic targets, indicating that inhibiting these genes in patients with high MPCDS scores may achieve favorable treatment outcomes.

To identify candidate drugs with higher sensitivity in patients with high MPCDS scores, we employed two different methods, using drug response data derived from CTRP and PRISM. First, we conducted differential drug response analysis between the high MPCDS group (top decile) and the low MPCDS group (bottom decile) to identify compounds with lower estimated AUC values in