Five-gene prognostic model based on autophagy-dependent cell death for predicting prognosis in lung adenocarcinoma

Lung cancer is a major global health issue, with around 2.2 million new cases and 1.8 million deaths each year. Adenocarcinoma is the most common subtype, accounting for approximately 40% of non-small cell lung cancer (NSCLC) cases, which make up over 85% of all lung cancer diagnoses¹. Research has revealed several key somatic mutations that drive lung cancer, such as EGFR, ALK, ROS1, and ERBB2. EGFR mutations are notably more common in East Asia (38–50%) compared to the Americas (24%) and Europe (14%), affecting approximately 38% of lung adenocarcinoma patients^2,3,4.

Recent clinical trials have demonstrated notable improvements in treatment outcomes with targeted therapies and immunotherapy. For example, the combination of osimertinib with EGFR-targeted chemotherapy achieved an objective response rate of 83%, though the 2-year progression-free survival rate was limited to 57%⁵. Immune checkpoint inhibitors targeting PD-1/PD-L1 and CTLA-4 have improved treatment outcomes for NSCLC patients. A recent Phase III trial (CheckMate 722, NCT02864251) evaluated the combination of Nivolumab with chemotherapy for patients who progressed after initial EGFR TKI treatment. The study found that this combination did not significantly improve progression-free survival (PFS) compared to chemotherapy alone in EGFR-mutated metastatic NSCLC, with median PFS of 5.6 versus 5.4 months⁶. Given these challenges, there is an urgent need to identify biomarkers that predict outcomes and guide treatment decisions, especially with immune checkpoint inhibitors (ICI). Autophagy, a key cellular process that involves the formation of autophagosomes to degrade and recycle damaged components, is essential for maintaining cellular homeostasis and responding to stress⁷. Autophagy-dependent cell death (ADCD) is a form of cell death associated with autophagy, affecting the balance between cell survival and death. This process helps cells survive stress and influences immune responses by releasing danger signals and Damage-Associated Molecular Patterns (DAMPs), which may enhance anti-cancer immune reactions⁸. Autophagy plays a protective role by helping cells manage stress and emergency situations. In ADCD, cells initiate autophagy through molecular processes involving Atg5, Beclin-1, and LC3, which can lead to apoptosis or other forms of cell death⁹. Apoptotic cell death releases internal cellular components, producing danger signals and Damage-Associated Molecular Patterns (DAMPs), such as organelles, proteins, and nucleic acids. These molecules can be recognized by surrounding cells and the immune system, triggering immune responses. This recognition recruits and activates antigen-presenting cells, initiating tumor-specific immune reactions¹⁰. Different stressors, including chemotherapy drugs, targeted therapies, radiotherapy, and intracellular pathogens, can induce autophagy-dependent cell death (ADCD). Research shows that combining chemotherapy with immunotherapy enhances tumor cell immunogenicity, stimulating anti-cancer immune responses and improving treatment outcomes¹¹. Therefore, inducing ADCD is considered a promising strategy to enhance immune responses against cancer and improve treatment efficacy. To gain a deeper understanding of ADCD’s role in the tumor immune microenvironment, single-cell RNA sequencing (scRNA-seq) is an essential tool for investigating cellular heterogeneity. This technology enables the analysis of gene expression profiles at the single-cell level, providing insights into the complex interactions and dynamic changes of various immune cells within the tumor microenvironment¹². An important area of research involves utilizing ADCD as a biomarker to stratify patient prognosis and predict responses to immunotherapy and chemotherapy.

In recent years, ADCD has garnered widespread attention in cancer research as a significant mode of cell death. ADCD plays a crucial role in regulating the survival and death of tumor cells, particularly in LUAD, where it holds substantial research value.

Autophagy plays a multifaceted role in regulating the tumor microenvironment (TME), involving immune modulation, antigen presentation, tumor resistance, and metabolic stress. For instance, studies have shown that the activation of autophagy in tumor-associated macrophages (TAMs) and dendritic cells (DCs) significantly influences tumor immune evasion and antigen presentation efficiency. Additionally, TME stressors such as hypoxia and nutrient deprivation activate autophagy, enhancing the survival of tumor cells in harsh conditions. Through interactions with stromal cells, autophagy further promotes tumor growth and metastasis¹³.

Despite the introduction of various LUAD prognostic models, most have overlooked the role of ADCD-related genes. Existing LUAD prognostic models typically rely on conventional biomarkers and clinical features, which have certain limitations in prediction accuracy and biological interpretation. Meanwhile, the tumor microenvironment (TME) is recognized for its critical role in tumor progression and prognosis. Recent studies have highlighted that immune cells and cytokines within the TME significantly impact tumor immune evasion and drug resistance¹⁴. Thus, integrating the mechanisms of ADCD with TME research may offer new perspectives for prognostic assessment in LUAD.

Further research into autophagy-dependent cell death (ADCD) is essential for understanding its mechanisms in cancer therapy. This knowledge could enable the development of personalized treatment strategies that might improve treatment outcomes and boost survival rates for cancer patients. However, research on the role of ADCD in immunotherapy and prognosis for lung adenocarcinoma (LUAD) remains limited. Using this information, we developed a distinct ADCD signature based on five genes (BIRC3, TAP1, SLAMF1, XCR1, and HLA-DMB), analyzed through both single-cell and bulk transcriptomic data from LUAD datasets. We established and validated the ADCD-related risk scoring system (ADCDRS), which proved effective in predicting disease prognosis. The system also correlated well with clinical characteristics, immune cell infiltration, the tumor microenvironment, and other pertinent factors. Additionally, we created a user-friendly nomogram for the ADCDRS to facilitate practical application. The objective of our research is to establish a foundational framework utilizing ADCD to assess immunotherapy effectiveness and predict outcomes for LUAD patients. This model aims to offer a more accurate basis for guiding immunotherapy and improving treatment strategies in LUAD.

Data collection

The Cancer Genome Atlas (TCGA) provided RNA-seq data and clinical information for lung adenocarcinoma (LUAD). Transcripts per kilobase million (TPM) values were used, excluding genes with average expression below 0.1 for data integrity. Samples lacking full clinicopathological data were excluded, resulting in a cohort of 598 TCGA-LUAD patients, including 539 tumor samples and 59 normal samples. Additionally, we accessed the LUAD single-cell RNA sequencing dataset (HRA) from the Genome Sequence Archive (GSA) database at [https://ngdc.cncb.ac.cn/gsa/]. Somatic mutation data were collected from TCGA using Mutation Annotation Format (MAF) and copy number variation (CNV) data were obtained from the UCSC Xena database. To assess the predictive accuracy of our signature for LUAD development, we utilized four datasets from the GEO database: GSE13213 (119 samples), GSE26939 (115 samples), GSE29016 (39 samples), and GSE30219 (86 samples). This study undertook a comprehensive systematic review and analysis of the existing literature to identify and curate genes implicated in autophagy-dependent cell death. We focused on 368 genes previously linked to ADCD for our analysis (Supplementary Table S1)¹⁵.

Single-cell RNA-Seq data collection and processing involved obtaining data from four LUAD patients through the Genome Sequence Archive (GSA) at the Beijing Institute of Genomics (BIG) Data Center, with accession number HRA 001130¹⁶. The analysis utilized the “Seurat” package¹⁷. Initially, quality assurance included filtering out cells with mitochondrial gene content less than 10%, erythrocyte gene content less than 3%, and ensuring a minimum of three genes expressed between 200 and 7500. Highly variable genes, totaling 2000, were identified for subsequent analysis. Gene expression data were normalized using the log-normalization method, where raw counts were scaled by the total expression per cell, multiplied by 10,000, and subsequently log-transformed. Following normalization, the data were further scaled and centered to ensure comparability across cells and genes. The “Harmony” package was employed to mitigate batch effects across the samples. Cell clustering was performed using the “FindClusters” and “FindNeighbors” functions, and results were visualized using the “UMAP” method. Cell types were categorized based on marker genes. The Seurat package’s ‘AddModuleScore’ function assessed the expression levels of specific gene sets in individual cells. Differential gene expression between groups was analyzed using the ‘FindMarkers’ function, employing the Wilcoxon test with an adjusted p-value threshold of < 0.05. Genes exhibiting differential expression associated with ADCD scores were identified, highlighting their role in ADCD at the single-cell transcriptome level. These genes were further analyzed using weighted gene co-expression network analysis (WGCNA). Additionally, cell interaction analysis was conducted using the “CellChat” R package¹⁸.

We first employed single-sample gene set enrichment analysis (ssGSEA) to measure the enrichment scores of specific gene sets in each sample. ssGSEA is a widely used technique for assessing the upregulation or downregulation of gene sets in individual samples. We used the ssGSEA method from the “GSVA” R package to calculate ADCD scores for each sample in the TCGA-LUAD dataset¹⁹. This analysis enabled the preliminary screening of ADCD-related gene sets, laying the foundation for further in-depth analysis. Next, we applied weighted gene co-expression network analysis (WGCNA) to explore co-expression patterns among genes and identify highly co-expressed gene modules²⁰. Using the “WGCNA” R package, we conducted WGCNA on the comprehensive RNA-seq data from TCGA-LUAD. We first determined an appropriate soft-threshold power (β) to ensure the network adhered to a scale-free topology. A weighted adjacency matrix was then constructed and transformed into a topological overlap matrix (TOM), followed by the calculation of gene dissimilarity (dissTOM). Gene clustering was performed using the dynamic tree cut algorithm, leading to the identification of several gene modules. We then conducted statistical analyses to identify the modules most strongly associated with ADCD scores. The genes within these modules were rigorously filtered to ensure their functional relevance to ADCD and their close association with the biological mechanisms of LUAD. These genes ultimately formed the core basis for our in-depth exploration of ADCD mechanisms and their role in LUAD.

Individuals were classified into distinct molecular clusters based on ADCDRs expression using the ConsensusClusterPlus R package²¹ and unsupervised clustering. In this study, univariate Cox regression analysis was conducted on 249 ADCD-related genes to identify those significantly associated with patient prognosis. Eighteen genes with p-values below 0.05 were selected, and a prognostic network was constructed to evaluate their interactions. Consensus clustering was then applied to LUAD patients from the TCGA cohort based on the expression levels of these 18 genes, achieving optimal clustering stability at k = 2. The clinical relevance of ADCDRs in LUAD was evaluated using the Kaplan-Meier method, and survival curves were generated using the survival and survminer packages in R. Principal component analysis was performed using the ggplot2 package. The proportions of immune and stromal cells in LUAD patients were quantified using CIBERSORT^22,23. Single-sample gene set enrichment analysis was employed to assess the enrichment scores of different immune cell types in each sample²⁴. We further explored the clinical implications of two ADCD-related subtypes by examining their association with key clinical and pathological factors such as age, T stage, and N stage, as well as prognosis in LUAD patients. Kaplan-Meier analysis was used to analyze differences in overall survival among these subtypes. Differentially expressed genes (DEGs) within ADCD-related subgroups were identified using the limma R package, with adjusted p-values less than 0.05 and fold changes exceeding 2. Enrichment analyses were conducted using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathways to identify significantly enriched pathways, demonstrated using the ClusterProfiler package²⁵.

To quantify ADCD patterns in individual tumors, we developed an ADCDRS (Autophagy-Dependent Cell Death-Related Score). Initially, differentially expressed genes (DEGs) underwent univariate Cox regression analysis to identify correlations with overall survival (OS) in LUAD patients. Subsequently, unsupervised clustering was used to explore prognosis-related ADCDR expression, categorizing patients into distinct subtype groups (A, B). LUAD patients were randomly divided into a training group (n = 258) and a testing group (n = 258) to compute a prognostic ADCDRS score. To avoid overfitting, Lasso Cox regression was performed using the ‘glmnet’ R package²⁶, with each independent variable’s path analyzed and model creation validated through 10-fold cross-validation. Before performing Lasso-Cox regression, continuous variables were standardized to a mean of 0 and a standard deviation of 1. The optimal penalty parameter (lambda) was determined through 10-fold cross-validation, with lambda.min (the value corresponding to the minimum cross-validated error) used for model fitting, and lambda.1se (the largest value within 1 standard error of the minimum) employed for comparison. To ensure the robustness of our findings, both internal and external validations were conducted. For internal validation, the dataset was randomly divided into a training set (50%) and a test set (50%). External validation employed independent datasets from the GEO database, including GSE13213 (119 LUAD samples), GSE26939 (115 LUAD samples), GSE29016 (39 LUAD samples), and GSE30219 (86 LUAD samples). These datasets underwent preprocessing and standardization similar to the training set, providing an unbiased evaluation of the model’s generalization capability. A prognostic ADCDRS score was derived from multivariate Cox analysis, where Coefi and Expi represent risk coefficients and gene expression levels. Patients in both training and testing sets were stratified into low-risk (ADCDRS score ≤ median) and high-risk (ADCDRS score > median) groups, followed by Kaplan-Meier survival analysis. External datasets from the GEO database were utilized to validate the model’s predictive performance through Kaplan-Meier survival analysis and receiver operating characteristic (ROC) curves. Additionally, to evaluate the performance of ADCDRS in comparison to other established signatures, we conducted a comprehensive analysis of previously published models. A total of 15 distinct features, including both mRNA and lncRNA characteristics, were selected for evaluation. These features were linked to various biological processes such as autophagy-driven cell death, ferroptosis, cuproptosis, pyroptosis, immune-related functions, anoikis, and thermoptosis. The concordance index (c-index) was employed to assess the predictive performance of ADCDRS relative to other signatures across multiple datasets. Statistical comparisons were then conducted to determine the efficacy of ADCDRS in outperforming these models.

Development and confirmation of a nomogram rating system

Using the rms package, we constructed a nomogram to predict overall survival by integrating clinically significant features with the ADCDRS score. The nomogram assigned scores to each relevant clinical feature, combined these scores, and generated a total predictive score. ROC curves were generated to assess the nomogram’s accuracy in predicting survival rates at 1, 3, and 5 years. Calibration of the model was conducted to evaluate the accuracy of predicted survival outcomes against observed survival rates during these specific time periods.

We employed the CIBERSORT algorithm to quantify levels of 22 different types of immune cells, examining the relationship between ADCDRS and immune cell infiltration in the LUAD tumor microenvironment. Additionally, the ESTIMATE algorithm was used to analyze stromal, immune, and ESTIMATE scores in groups stratified by high and low ADCDRS risk²⁷.

To assess personalized treatment options, we utilized the oncoPredict R package²⁸ to predict chemotherapy sensitivity in LUAD patients based on ADCDRS risk scores. oncoPredict aligns patient tissue gene expression profiles with cancer cell profiles to calculate the half-maximal inhibitory concentration (IC50). We compared drug IC50 variances between high-risk and low-risk groups using the Wilcoxon test, with a significance level set at p < 0.05.

Cell culture

BEAS-2B, a human bronchial epithelial cell line, and two human lung adenocarcinoma cell lines (A549 and PC9), were sourced from the Cell Bank of the Chinese Academy of Sciences. These cells were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin from Gibco, Thermo Fisher Scientific, Inc. Cultures were maintained under standard conditions at 37 °C with 5% CO2 and 95% humidity.

Isolation of RNA and real-time reverse transcription polymerase chain reaction for XCR1 gene expression

Total RNA was reverse-transcribed into complementary DNA (cDNA) using the PrimeScript™ RT Reagent Kit from Takara, following the manufacturer’s instructions. Real-time RT-PCR was conducted on a Bio-Rad CFX96 system to quantify XCR1 expression. Fold changes in gene expression were determined using the 2-ΔΔCT method, with GAPDH serving as the internal control gene.

Statistical analysis

R software (version 4.3.2) and GraphPad Prism 8 were employed for all statistical analyses in this study. Clinical characteristics differences between the training and internal validation sets were analyzed using the chi-square test. Non-normally distributed variables were assessed with the Wilcoxon test, a non-parametric approach. For differential expression gene (DEG) analysis, significance of DEGs was determined using FDR-corrected p-values. Survival analysis for assessing differences in overall survival (OS) among patient subgroups was performed using the Kaplan-Meier method and log-rank tests. Independent prognostic factors were identified through univariate and multivariate Cox regression analyses. Model performance was evaluated using ROC curve analysis and area under the curve (AUC) calculations, employing the “timeROC” package in R. Spearman correlation analysis was utilized to explore relationships between risk scores and immune cell infiltration. Statistical significance was set at a p-value less than 0.05, unless otherwise specified. For missing data, features with over 20% missing values were excluded from the analysis. For features with less than 20% missing data, mean imputation was applied for continuous variables, and mode imputation was used for categorical variables. All imputation was conducted using the “mice” package in R.

Patient information and tissue sample source

Tumor and healthy tissues were obtained from patients at Tianjin Medical University Cancer Institute and Hospital. Ethical approval for this study was granted by the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (Ethics No. E20220942). All procedures performed in this study involving human participants were conducted in accordance with the ethical standards of the institutional and national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Informed consent was obtained from all individual participants included in the study.

Autophagy-dependent cell death characteristics in single-cell transcriptome

We conducted single-cell RNA sequencing on 90,260 cells obtained from four patients with LUAD. Batch effects were corrected using the Harmony package (Supplementary Fig. 1a,1b), effectively integrating the data from all four samples. Principal component analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) were employed to reduce the dimensionality based on the top 2000 variable genes. Cell clustering into 50 groups at a resolution of 1.5 resulted in the identification of eleven major cell types using marker genes: T cells, AT2-like cells, NK cells, Macrophages, Ciliated cells, DCs, Granulocytes, AT2/AT1 cells, Club/Basal cells, Fibroblasts/Endothelial cells, and B cells. The bubble plot visualized the marker genes defining each cell cluster (Fig. 1a). The UMAP plot shows the clustering of cells (Fig. 1b).

Characteristics of Autophagy-Dependent Cell Death (ADCD) in Single-Cell Transcriptomes. (a) Cell annotation using bubble plots based on marker genes for different cell types in single-cell data. (b) UMAP-1 plot displaying cell types identified by marker genes. (c) ADCD activity scores in various cell types. (d)Grouping of ADCD activity scores in different cell types.

Full size image

To assess autophagy-dependent cell death (ADCD) levels across different cell populations, we utilized the ‘AddModuleScore’ function in Seurat software to evaluate the expression of 368 ADCD-associated genes across all cells (Fig. 1c). Elevated ADCD activity was notably observed in Macrophages, DCs, and AT2-like cells among the identified cell types. Based on this activity, cells were stratified into high-ADCD and low-ADCD groups, revealing 703 differentially expressed genes (DEGs) between these groups for further analysis (Fig. 1d).

The study also involved the identification of central modules and ADCD-related genes in bulk RNA-seq

The ssGSEA algorithm is frequently employed to evaluate changes in the activity of biological processes and pathways within individual samples. In our study, we used ssGSEA to calculate ADCD activity scores for each sample in the TCGA dataset, which were then used as phenotype data for subsequent Weighted Gene Co-expression Network Analysis (WGCNA). The UMAP plot illustrates the scores of various cell subpopulations (Fig. 2a). Applying WGCNA to the TCGA dataset, we identified modules that were strongly associated with ADCD scores. We constructed a co-expression network using 703 ADCD-related differentially expressed genes (DEGs) identified at the single-cell level, excluding outliers. The top six genes within this network were visualized (Fig. 2b and c).

Identification of ADCD-Related Module Genes. (a) Cells are divided into high ADCD group and low ADCD group. (b) Violin plots of the top six ADCD-related DEGs identified at the single-cell sequence level. (c) UMAP plots of the top six ADCD-related DEGs identified at the single-cell sequence level. (d) Module-trait heatmap showing that the MEblue module is closely related to ADCD traits. (e) Scatter plot showing the relationship between gene significance (GS) and module membership (MM) in the yellow module.

Full size image

To achieve a scale-free network topology, we selected a soft threshold power of 12 (R² = 0.865) during the thresholding process. With a minimum gene count of 60 and a MEDissThres of 0.25, we identified nine distinct modules (Fig. 2d). Our results revealed a strong correlation between the MEyellow module and ADCD scores in bulk RNA-seq data (correlation coefficient = 0.76; Fig. 2d). Additionally, the scatter plot for the yellow module showed a strong correlation between gene significance (GS) and module membership (MM) (correlation coefficient = 0.87, p-value = 8.1e-78; Fig. 2e), suggesting the functional relevance of genes within the yellow module for autophagy-dependent cell death. The 249 genes within the yellow module were identified as ADCD-related genes (ADCDRgenes), which are involved in ADCD as evidenced by both bulk and single-cell transcriptomic analyses.

Patient typing based on differential expression of ADCD-related genes

An analysis was conducted using univariate Cox regression on 249 ADCDRgenes, identifying 18 genes with significant p-values below 0.05 (see Fig. 3a). A prognostic network graph was constructed for these genes (Fig. 3b), revealing PIK3AP1, SH2B3, SLCO2B1, FPR3, and PIK3R5 as protective factors, while the remaining ADCD-related genes acted as risk factors. The network diagram demonstrated strong correlations among these 18 genes.

Biological and Clinicopathological Characteristics of ADCD Subtypes. (a) Forest plot of univariate Cox regression analysis for ADCDR genes, identifying 18 significant genes with p-values less than 0.05. (b) Interaction between ADCD subtypes (red and blue indicate positive and negative correlations, respectively; the intensity of the correlation is indicated by the color shade). (c) Consensus matrix heatmap for two clusters (k = 2). (d) Kaplan-Meier OS curves for specific subtypes. (e) PCA analysis showing significant transcriptomic differences between the two subtypes. (f) Scatter plot showing the relationship between gene significance (GS) and module membership (MM) in the yellow module. (g) Scatter plot showing the relationship between gene significance (GS) and module membership (MM) in the yellow module.

Full size image

Consensus clustering analysis was performed on LUAD patients from the TCGA cohort with survival data to explore the relationship between expression levels of the 18 ADCD-related genes and LUAD subtypes. Optimal clustering stability was achieved at k = 2 (Fig. 3c). Patients were stratified into two clusters based on differential gene expression. Kaplan-Meier analysis showed significantly improved overall survival in cluster 2 compared to cluster 1 (Fig. 3d). A heatmap (Supplementary Fig. 2) depicted the clinical characteristics of patients in each cluster, while PCA highlighted distinct features separating cluster 1 from cluster 2 (Fig. 3e).

Using the CIBERSORT algorithm, we quantified immune cell infiltration in each sample to assess differences between cluster 1 and cluster 2 (Fig. 3f). Cluster 1 exhibited higher levels of plasma cells, CD8 + T cells, memory-activated CD4 + T cells, follicular helper T cells, and regulatory T cells, whereas CD56dim NK cells predominated in cluster 2, known for antibody-dependent cell-mediated cytotoxicity (ADCC). GSVA analysis using ‘c2.cp.kegg.v7.4’ revealed statistically enriched immune-related biological processes and pathways in the high-expression group of ADCDRgenes (Fig. 3g).

To explore pathway differences between the two subtypes, we analyzed differentially expressed genes (DEGs) based on high and low ADCDRgene expression scores (|logFC| > 0.5 and adjusted p < 0.05). GO and KEGG enrichment analyses highlighted significant pathways and molecular functions. Biological process analysis from the GO database indicated DEGs were primarily involved in immune response regulation, leukocyte interaction, cell adhesion control, signaling receptor complexes on the plasma membrane, T cell receptor complexes, and immune receptor activity signaling pathways (Fig. 4a and b). KEGG analysis identified DEGs associated with cytokine-cytokine receptor interactions, chemokine signaling pathways, hematopoietic cell lineage, and cell adhesion molecules signaling pathways (Fig. 4c and d).

ADCD Subtypes Associated with TME Infiltration. (a-b) DEG enrichment analysis between the two ADCD subtypes using GO. *p < 0.05, **p < 0.01, ***p < 0.001. (c-d) DEG enrichment analysis between the two ADCD subtypes using KEGG. *p < 0.05, **p < 0.01, ***p < 0.001. (e)C onsensus matrix heatmap defining the two clusters (k = 2). (f) Kaplan-Meier OS curves for the two gene subtypes. (g) Differences in clinicopathological characteristics between the two gene subtypes. (h) Expression changes of eighteen ADCDR genes in the two gene subtypes. ***p < 0.001.

Full size image

Gene subtypes identified using DEGs

Consensus clustering was employed to classify 516 LUAD patients into two molecular genetic groups based on prognostic genes, aiming to explore the relationship between distinct gene expression profiles and molecular subtypes of LUAD. When set at k = 2, we identified subtypes A (n = 307) and B (n = 209), indicating two distinct subclasses of LUAD cases (Fig. 4e and Supplementary Table S2). We further investigated the association between clinical characteristics and gene subtypes among LUAD patients (Fig. 4g). Kaplan-Meier survival analysis revealed significantly reduced overall survival in patients belonging to genetic subtype A compared to those in genetic subtype B (p < 0.01; Fig. 4f).The expression of ADCDRgenes varied significantly between the two ADCD gene subtypes, consistent with the ADCD pattern observed (Fig. 4h).

Construction of a predictive model

The 516 participants were split into two cohorts, a training cohort and a validation cohort, each consisting of 258 individuals, using the R caret package with a 1:1 ratio. In the training cohort dataset, we performed LASSO-Cox analysis using 67 differentially expressed genes (DEGs) associated with subtype prognosis, identifying additional prognostic genes and developing a model.

Through ten-fold cross-validation and by minimizing the deviance of the partial likelihood, we constructed a prognostic model comprising five genes selected based on the optimal λ value from LASSO analysis (see Fig. 5a and b). The model is represented by the following equation for calculating risk scores, considering the coefficients and gene expression levels:

Construction and Validation of the Predictive Model. (a-b) LASSO regression visualization to obtain the optimal λ when partial likelihood deviance reaches its minimum. (c-d) Distribution of patient survival status and ADCDRS scores in the training and testing sets. (e-p) Kaplan-Meier OS curves and ROC curves for the training set, testing set, validation set GSE13213, validation set GSE26939, validation set GSE29016, and validation set GSE30219, respectively, to verify the differences between high-risk and low-risk groups and predict 1-year, 3-year, and 5-year OS in the cohorts. *p < 0.05, ***p < 0.001.

Full size image

Risk score = (BIRC3 × 0.217464258073396) - (SLAMF1 × 0.422658700457643) + (TAP1 × 0.380601938911013) - (XCR1 × 0.294615486530056) - (HLA-DMB × 0.304131811097357).

Subsequently, individuals in the training group were classified into high-risk (n = 129) and low-risk (n = 129) categories using the median risk score. Notable variations in gene expression levels between the high-risk and low-risk subgroups were observed in eleven marker genes (Supplementary Fig. 3a).

Supplementary figures illustrate the distribution of patients across the two risk score groups, ADCD subtypes, and gene subtypes (Supplementary Fig. 3b). This figure provides a comprehensive overview of patient categorization based on these variables, highlighting distinct distribution patterns and relationships between them. Additionally, significant differences in risk scores were observed between ADCD subtypes, with cluster 2 exhibiting significantly higher risk scores than cluster 1 (Supplementary Fig. 3c). This figure further underscores the disparity in risk levels associated with different ADCD subtypes, demonstrating that cluster 2 is characterized by a notably higher risk profile compared to cluster 1. This differential risk scoring highlights the importance of ADCD subtype classification in assessing patient prognosis and stratification.

Evaluation of the ADCDRS model

We evaluated our model by comparing the outcomes between high-risk and low-risk categories. The mortality rate among LUAD patients increased proportionally with their risk scores (Fig. 5c and d), a trend consistently observed in both the test and entire datasets (Fig. 5d, Supplementary Fig. 4a). Specifically, BIRC3 and TAP1 expression levels were significantly elevated in the high-risk subgroup across the training, test, and entire datasets, suggesting these ADCD-related genes as potential indicators of adverse outcomes. Conversely, SLAMF1, XCR1, and HLA-DMB expression levels were significantly higher in the low-risk subgroup, indicating their potential protective roles (Supplementary Fig. 4b).

We evaluated our model by comparing outcomes between high-risk and low-risk categories. The mortality rate among LUAD patients increased proportionally with their risk scores (Fig. 5c and d), a trend consistently observed in both the test and full datasets (Fig. 5d, Supplementary Fig. 4a). Notably, BIRC3 and TAP1 expression levels were significantly elevated in the high-risk subgroup across the training, test, and full datasets, suggesting these ADCD-related genes as potential indicators of adverse outcomes. Conversely, SLAMF1, XCR1, and HLA-DMB expression levels were significantly higher in the low-risk subgroup, indicating their potential protective roles. Supplementary Fig. 4b provides a detailed view of these expression patterns, highlighting the consistent elevation of protective genes in the low-risk group and risk-associated genes in the high-risk group, further reinforcing the biological relevance of these ADCD-related genes in stratifying patient outcomes.

We evaluated our model by comparing the outcomes between high-risk and low-risk categories. The mortality rate among LUAD patients increased proportionally with their risk scores (Fig. 5c and d), a trend consistently observed in both the test and entire datasets (Fig. 5d, Supplementary Fig. 4a). Specifically, BIRC3 and TAP1 expression levels were significantly elevated in the high-risk subgroup across the training, test, and entire datasets, suggesting these ADCD-related genes as potential indicators of adverse outcomes. Conversely, SLAMF1, XCR1, and HLA-DMB expression levels were significantly higher in the low-risk subgroup, indicating their potential protective roles. We provides a detailed visualization of these differential expression patterns, highlighting the pronounced differences between high-risk and low-risk subgroups (Supplementary Fig. 4b). This figure underscores the distinct expression profiles associated with each risk category and reinforces the potential of these genes as reliable biomarkers for prognostic stratification in LUAD.

Kaplan-Meier curves (Fig. 5e and f) and time-dependent ROC curves (Fig. 5g and h) were plotted to assess the sensitivity of our prognostic model. Significant differences in survival curves were observed between high-risk and low-risk groups in both the training and test datasets. In the training set, the AUC values were 0.68, 0.73, and 0.71 for 1, 3, and 5 years, respectively, while in the test set, they were 0.73, 0.7, and 0.66. Additionally, our model was validated using four external datasets obtained from the GEO database (GSE13213, GSE26939, GSE29016, GSE30219). Kaplan-Meier curves (Fig. 5i-l) and ROC curves (Fig. 5m-p) demonstrated the superior predictive capabilities of our model.

To assess the performance of ADCDRS relative to other established signatures, we conducted a comprehensive evaluation of published models. Ultimately, we incorporated 15 distinct features, including both mRNA and lncRNA characteristics (Supplementary Table S3). These features are associated with various biological processes such as autophagy-driven cell death, ferroptosis, cuproptosis, pyroptosis, immune-related functions, anoikis, and thermoptosis. The c-index comparison of ADCDRS against other signatures demonstrated that ADCDRS consistently outperformed nearly all models across multiple datasets (Fig. 6a-g). This suggests that our signature holds significant promise for broad application.

A comparison of lung adenocarcinoma prognostic signatures based on gene expression. The C-index analysis of ADCDRS and 15 published signatures was conducted across the TCGA, GSE13213, GSE26939, GSE29016, GSE30219, GSE31210, and GSE42127 cohorts. Statistical test: two-sided z-score test. * p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001.

Full size image

Creation of a nomogram to predict survival

To enhance the predictive model and facilitate efficient communication, we integrated various factors visually and utilized collected data to develop a nomogram predicting the expected survival of LUAD patients at 1, 3, and 5 years using rms software (Supplementary Fig. 4c). Each patient’s total score was determined based on gender, age, risk level (high or low), and TNM stage, with higher scores indicating poorer prognosis. Calibration curves showed that the accuracy of predictions was much higher than that of the theoretical model (Supplementary Fig. 4d).

Analysis of the ADCDRS model and tumor microenvironment

We utilized the CIBERSORT algorithm to assess immune infiltration in LUAD samples from the TCGA cohort. Notable differences were observed in immune cell infiltration patterns between high-risk and low-risk cohorts. Specifically, the high-risk group exhibited elevated levels of CD8 + T cells, CD4 + T cells, and macrophages, suggesting their potential involvement in tumor progression. Conversely, the low-risk group showed higher levels of B cells, NK cells, and dendritic cells, indicating a possible role in anti-tumor immunity (Fig. 7b-l). Additionally, building on our previous analysis, macrophages and dendritic cells exhibited significantly elevated expression levels of ADCD-related genes. Therefore, ADCD-related genes may play a critical role in either promoting antitumor immunity or facilitating tumor progression, depending on the context. Further supporting this notion is the pivotal role of resting CD4 + memory T cells within the immune system. As a key component of immunological memory, these cells possess the ability to recognize prior antigen exposure and initiate a robust immune response upon reactivation. These cells play a key role in immune surveillance, identifying and destroying tumor cells, and assisting in immune clearance when situated in tumors²⁹. However, some tumors can evade immune detection due to reduced antigen presentation, altered antigen processing, and enhanced immunosuppressive pathways, which dampen or evade the response of resting CD4 + memory T cells. By enhancing their activity and boosting their interaction with other immune cells, these cells could potentially improve immunotherapy outcomes, allowing the immune system to fight cancer more effectively³⁰. Therefore, an inverse relationship with resting CD4 + memory T cells suggests a more favorable prognosis. Our study also investigated the correlation between immune cell levels and five specific genes. These genes significantly influenced the majority of immune cells (Fig. 7a). A significant correlation was found between Tumor Mutation Burden (TMB) levels and ADCDRS scores (P < 0.001, Fig. 8b). ESTIMATE analysis revealed that the low ADCDRS score group exhibited significantly higher stromal, immune, and overall ESTIMATE scores compared to the high ADCDRS score group (Fig. 8a), indicating lower tumor purity and a greater presence of immune cells in the low ADCDRS score group. Additionally, TMB analysis showed a significant correlation between expected TMB levels and ADCD gene subtypes (R = 0.25, P < 0.001; Fig. 8c).

Relationship Between ADCDRS Scores and Immune Cell Quantities. (a) The correlation between the quantity of immune cells and the five genes within the model. (b-i) The correlation between immune cell types and ADCDRS scores.

Full size image

Relationship Between ADCDRS Scores, Tumor Microenvironment, Tumor Mutation Burden, and Drug Sensitivity Analysis. (a) The relationship between ADCDRS scores and stromal and immune cells in the tumor microenvironment. (b) Violin plot depicting differences in Tumor Mutation Burden (TMB) scores between high and low ADCDRS risk groups. (c) Spearman correlation analysis between ADCDRS scores and TMB. (d) The correlation between CSC index and CRG scores. (e-g) The relationship between ADCDRS scores and drug sensitivity (BI-2536, docetaxel, paclitaxel, BMS-754801, lapatinib, tamoxifen).

Full size image

Drug sensitivity analysis and CSC index

Research indicates a strong connection between the ADCDRS score and the CSC index (R = 0.39, P < 0.001) (Fig. 8d). This suggests that lung adenocarcinoma cells with higher gene retention scores demonstrate increased stem cell characteristics and reduced differentiation. We analyzed drug sensitivity using the “oncoPredict” R package to explore potential therapeutic targets. High-risk patients exhibited increased sensitivity to multiple chemotherapy drugs, including paclitaxel, cisplatin, and docetaxel, compared to low-risk patients. Additionally, high-risk patients showed heightened responsiveness to targeted treatments such as EGFR inhibitors and ALK inhibitors (Fig. 8e-j). These findings suggest that our ADCDRS model may aid in identifying patients who could benefit from specific therapeutic strategies.

Analysis of XCR1 gene expression discrepancies in tumor and normal tissues

XCR1 is crucial in modulating anti-tumor immune responses and shows a significant correlation with prognosis in our model. The role of XCR1 gene in cancer is not well studied. In order to delve into this subject, we examined the mRNA expression levels of the XCR1 gene in a normal lung epithelial cell line and two lung cancer cell lines. Our findings revealed that the expression of XCR1 was significantly elevated in the normal lung epithelial cell line in comparison to the lung cancer cell lines (Fig. 9a). To further confirm the differential expression of the XCR1 gene in normal versus cancer tissues, we collected RNA samples from 10 sets of corresponding normal and cancer tissues belonging to individuals diagnosed with lung cancer. Through an examination of the mRNA expression levels of the XCR1 gene, we made a notable finding. In the majority of cases (7 out of 10), the expression of XCR1 was noticeably reduced in cancerous tissues in comparison to their healthy tissue counterparts (Fig. 9b). These findings align perfectly with the predictions made by our model.

Differential Expression of XCR1 Gene in Normal and Tumor Tissues. (a) Expression of XCR1 gene in normal cells (BEAS-2B) and tumor cells (A549, PC9). (b) Expression of XCR1 gene in 10 pairs of normal and tumor tissues.

Full size image

Lung adenocarcinoma (LUAD) is the most common histological subtype of non-small cell lung cancer (NSCLC)³¹. Surgical removal is often the preferred treatment for many early-stage LUAD patients, whereas immunotherapy offers potential benefits for those with advanced-stage disease, especially stage IV³¹. Despite notable progress in lung cancer research, including advancements in targeted therapies, immunotherapies, multidisciplinary approaches, and novel drugs, challenges in treatment and prognosis remain. Consequently, there is an urgent need to explore new prognostic biomarkers and therapeutic strategies for LUAD³².

Regulation of cell apoptosis remains a central focus in cancer therapy, with significant advancements achieved in targeted anti-cancer drug development. Recent research highlights the complex interactions among non-apoptotic cell death pathways, which impact cancer progression and treatment outcomes across diverse types of human cancers. Mechanisms of controlled cell death, including apoptosis, autophagy-dependent cell death (ADCD), ferroptosis, necroptosis, pyroptosis, and immunogenic cell death, offer promising strategies for advancing cancer treatment.

The discovery of autophagy-mediated cell death has significantly advanced cancer immunotherapy research, emerging as a pivotal focus in oncology. However, previous studies have largely overlooked the relationship between LUAD and ADCD. Previous research has predominantly focused on autophagy-dependent cell death (ADCD) in cancers like colorectal and liver cancer^33,34. Our study, however, addresses a critical gap by investigating the interplay between ADCD and immune function in lung adenocarcinoma (LUAD). Instead of experimental approaches to examine autophagy pathways, we leverage advanced bioinformatics techniques, integrating multi-omics data analyses—such as single-cell RNA sequencing and bulk transcriptomic profiling—to construct an ADCD-related prognostic model. Our findings suggest that prognostic models focusing on immune system function can effectively predict outcomes for LUAD patients, including immune cell infiltration levels in lung cancer. Furthermore, these models can aid healthcare professionals in tailoring personalized treatment strategies and provide a foundational basis for comprehensive research on lung adenocarcinoma and its tumor microenvironment.

Liu Ming’s research team at West China Hospital, Sichuan University, has elucidated that autophagy-dependent cell death (ADCD) is a specific type of cell death regulated exclusively by autophagy, which can be modulated through both chemical and genetic interventions³². The mechanisms of ADCD involve several pathways, including AMPK, PI3K/AKT/mTOR, and the Vps34 complex. AMPK regulates autophagy by inhibiting mTOR, activating ULK1, and phosphorylating beclin-1. The PI3K/AKT/mTORC1 pathway, frequently activated in various cancer cells, consistently promotes ADCD upon inhibition of AKT/mTOR. The VPS34-beclin-1 complex is essential for controlling autophagy, with elevated beclin-1 levels facilitating autophagic cell death. Furthermore, protein kinases such as JNK, MEK, and ERK can also induce ADCD, underscoring its role as “autophagy-dependent cell death”^35,36. Research demonstrates that ADCD occurs through distinct mechanisms in various organisms. In Dictyostelium discoideum (slime mold), ADCD occurs independently of traditional apoptotic pathways. During the transition from a nutrient-rich state to fruiting body formation, programmed cell death (PCD) is essential and can be inhibited by genetic modifications that target autophagy. Here, starvation-induced autophagy triggers cell death through an additional signal, differentiation-inducing factor (DIF-1), highlighting autophagy’s role in cell death beyond apoptosis. Similarly, in Drosophila (fruit flies), during the transition from larva to pupa, the hormone ecdysone facilitates the degradation of old larval tissues, such as the midgut and salivary glands. Although apoptotic genes are upregulated and caspase activity is high, midgut degradation continues normally, though it is delayed if autophagic components are depleted. These observations, when integrated with the role of autophagy in the tumor microenvironment, suggest a more intricate role for autophagy in tumor development than previously understood³⁷.

Despite this, ADCD remains understudied in lung cancer, with little exploration into its correlation with LUAD.

To bridge this research gap, our study developed a predictive model for ADCD to forecast LUAD outcomes and investigate its association with the tumor immune microenvironment (TIME) and response to immunotherapy. Initially, we screened 18 ADCD-associated genes and derived a predictive signature termed ADCDRS to assess post-treatment outcomes in LUAD patients. This signature comprises five genes: BIRC3, TAP1, SLAMF1, XCR1, and HLA-DMB. Using LASSO, univariate, and multivariate Cox regression analyses across training, testing, and comprehensive datasets, we identified BIRC3 and TAP1 as prognostic risk factors, while SLAMF1, XCR1, and HLA-DMB emerged as protective factors.

BIRC3, or Baculoviral IAP Repeat Containing 3, belongs to the inhibitor of apoptosis protein family, pivotal in apoptosis regulation, cell survival promotion, and intracellular signaling maintenance³⁸. Dysregulated BIRC3 expression significantly correlates with tumorigenesis, progression, and treatment resistance, notably through NF-κB pathway regulation, influencing inflammatory and immune responses crucial in tumor progression³⁹. BIRC3 has been integrated into prognostic models for breast and pancreatic cancers.

TAP1, also known as Transporter Associated with Antigen Processing 1, encodes a protein essential for antigen processing and presentation critical for immune system function⁴⁰. TAP1 expression levels closely correlate with tumor initiation, growth, treatment response, and prognosis, highlighting its role in regulating tumor immune evasion⁴¹. Targeting TAP1 expression or function could potentially enhance immune recognition and tumor cell clearance, thereby improving immunotherapy efficacy.

XCR1, or X-C motif chemokine receptor 1, plays a pivotal role in modulating the immune microenvironment of tumors, influencing immune cell infiltration and aiding in tumor immune evasion⁴². Studies suggest XCR1 as a promising candidate for cancer immunotherapy by enhancing dendritic cell (DC) performance, leading to heightened activation and proliferation of tumor-specific T cells, thereby boosting immunotherapy effectiveness^42,43.

SLAMF1, also referred to as CD150, is an immunoglobulin superfamily member potentially linked to tumor prognosis based on its expression levels. Low SLAMF1 expression often correlates with poor tumor prognosis, suggesting its utility as a prognostic biomarker⁴⁴.

HLA-DMB, a component of the human leukocyte antigen (HLA) complex, is crucial for antigen processing and presentation⁴⁵. While its direct role is limited, HLA-DMB expression levels may impact tumor prognosis. High HLA-DMB expression correlates with favorable outcomes, whereas low or abnormal expression may associate with tumor malignancy and poor prognosis.

Our study constructed the ADCDRS score based on these genes, demonstrating its robust predictive power for LUAD patient prognosis through survival and ROC analyses across various datasets. Additionally, we explored the differential response to common anticancer medications between high and low ADCDRS score groups, revealing potential implications for personalized treatment strategies.

This study developed a prognostic scoring tool, ADCDRS, and confirmed its excellent performance in predicting patient prognosis trends. Comparing low-risk and high-risk ADCDRS values revealed significant disparities in survival rates, clinical features, genetic mutation patterns, tumor microenvironment characteristics, cancer stem cell levels, and resistance to medication. A nomogram was devised using patient characteristics and ADCDRS scores to enhance the tool’s effectiveness and user-friendliness, providing insights into the molecular-level pathological mechanisms of lung adenocarcinoma and offering new perspectives for future cancer therapies. Additionally, given the prevalence of chemotherapy and targeted therapy in treating advanced LUAD, this research investigated how low-risk and high-risk groups respond to commonly prescribed anticancer medications. High ADCDRS scores correlated with elevated IC50 values for BI-2536, docetaxel, paclitaxel, BMS-754801, lapatinib, and tamoxifen, suggesting potentially better treatment responses in these individuals⁴⁶.

Our observations indicate that XCR1 expression is significantly elevated in tumor cells⁴⁷. We established a correlation between XCR1 expression and both ferroptosis and autophagy-dependent cell death, leading us to concentrate on XCR1 for further studies. qRT-PCR analyses of cell lines and human tissues demonstrated significant differences in XCR1 expression between tumor and normal tissues, highlighting XCR1 as a promising new therapeutic target for LUAD patients. This study introduces an innovative approach, differing from previous research that primarily emphasizes predictive features. We conducted a comprehensive analysis that integrates genomics, single-cell transcriptomics, and high-throughput transcriptomics to examine the characteristics of ADCDRS. Our investigation unveiled the molecular basis and potential mechanisms of ADCDRS across different levels and confirmed its significant association with the development and prognosis of LUAD. These results provide a strong biological foundation for incorporating ADCDRS into personalized medicine strategies. Furthermore, we employed an innovative bioinformatics method combining AddModuleScore, ssGSEA, and WGCNA algorithms to identify ADCD-related genes at both single-cell and bulk transcriptomic levels. This approach facilitates the identification of potential therapeutic targets within the ADCD environment, providing new insights and research directions for investigating ADCD in LUAD.

However, the study has several inherent limitations. Given our reliance on bioinformatics analysis, further experimental validation is necessary. Our in vitro cell experiments are preliminary, and additional in vivo studies are required to validate our findings. Due to sample size limitations, data may exhibit some bias. While we conducted external validation with independent datasets from the GEO database, several potential limitations must be considered. Differences in clinical characteristics and data quality between external datasets and our primary cohort could introduce batch effects or sampling bias. Variations in data collection protocols and processing methods might impact the reproducibility and generalizability of our findings. Moreover, external datasets may exhibit selection bias and may not fully capture the diversity of patient populations or tumor subtypes present in our primary cohort. To address these limitations, future research should incorporate multiple independent datasets and perform sensitivity analyses to evaluate the robustness of results across various settings.

In summary, the use of ADCDRS features can independently predict the prognosis of LUAD patients, providing valuable insights for guiding immunotherapy in LUAD patients.

Through extensive bioinformatics analyses, this study developed a prognostic model utilizing five key ADCD-related genes (BIRC3, TAP1, SLAMF1, XCR1, HLA-DMB) to predict disease outcomes and treatment efficacy in LUAD. We also explored the potential of the autophagy-dependent cell death-related signature (ADCDRS) as a predictive biomarker for therapeutic response. Our findings underscore the clinical importance of ADCDRS and offer a strong basis for future research focused on developing personalized diagnostic and treatment strategies for LUAD patients.