Lysosome-dependent cell death reveals a prognostic signature in colorectal cancer via integrated analysis of scRNA-seq and bulk RNA-seq data

Introduction

Colorectal cancer (CRC) ranks fourth in mortality rate among cancers, and nearly one million people are killed every year (1). Pathogenesis of CRC is highly complicated and involves several pathophysiologic mechanisms, like cell differentiation, abnormal cell proliferation, and cell apoptosis (2). With respect to molecular pathological mechanisms, abnormal activity of Wnt/β-catenin signaling increases cancer cell proliferation and differentiation speed and interferes with CRC development and treatment response (3). Epidermal growth factor receptor (EGFR)-mediated carcinogenic signals promote CRC malignant progression via a dual mechanism: on the one hand, they disrupt the cell cycle to stimulate cell proliferation and angiogenesis; on the other hand, they inhibit the apoptotic pathway through the phosphorylation of BAD protein (4). Recently, programmed cell death (PCD) dysregulation has been viewed as a pivotal node in CRC occurrence and development. Ferroptosis, a common cell death form, has drawn great attention due to successful elimination of cancer cells resistant to cell death in CRC (5). At present, known forms of PCD, like apoptosis, ferroptosis, necrotizing apoptosis, autophagy dependent, pyrocytosis and immunogenic cell death, are promising target systems for anticancer treatment and feasible new directions for novel combined CRC treatment strategies (6).

Lysosome-dependent cell death (LDCD), an important PCD form, is a distinct cell death form highly important for cellular activities (7). Studies have shown that lysosomal homeostasis imbalance is widespread in malignant tumor cells. In tumor microenvironment (TME), lysosomes affect tumor cell occurrence, growth, invasion and drug resistance by regulating activities of tumor-associated macrophages, fibroblasts, and dendritic and T cells (8). Immune escape and tumorigenesis are triggered by maintaining lysosomal activity (9). Furthermore, at the level of cell death regulation, lysosomes are involved in regulating cancer cell death type I (apoptosis), II (autophagy), or III (necrosis) (10). It activates sphingomyelinase through the interaction of RIPK1 and RIPK3, induces lysosomal membrane permeation, and thereby regulates apoptosis (11). Multiple studies have confirmed that specifically inducing LDCD can effectively inhibit tumor progression, which is important for tumor treatment (12). Addressing lysosomal impairments constitutes a potent therapeutic avenue for CRC, especially in the context of mitigating insensitivity to conventional drugs. The distinct addiction of tumor cells to lysosomal turnover creates a unique therapeutic opportunity that can be leveraged for personalized medical approaches (13).

In this study, to explore the regulatory mechanism of LDCD in CRC, we incorporated single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing (bulk RNA-seq) data to construct an LDCD activity scoring system and analyzed the dynamic changes in LDCD activity during differentiation trajectory of tumor cells and LDCD-mediated cell interaction network. A prognostic signature was established to clarify the relationships between key LDCD genes related to CRC and survival prognosis and to assess correlations between prognostic markers and clinical characteristics, the immune microenvironment and drug responses. In summary, our research strategy, ranging from single-cell resolution to the population level, delves deeply into the biological functions of LDCD-related genes and offers a basis theoretically and potential targets for developing individualized treatment plans. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-aw-2219/rc).

Methods

Data acquisition

We acquired The Cancer Genome Atlas (TCGA)-COAD mRNA-seq data, clinical data and mutation data, and TCGA-READ mRNA-seq data from UCSC Xena website (https://xenabrowser.net/). TCGA-COAD cohort comprised 512 specimens, namely, 471 cancer specimens and 41 normal specimens. TCGA-READ cohort comprised 177 specimens, i.e., 167 cancer specimens and 10 normal specimens. Furthermore, we gained RNA-seq data for CRC patients from the GSE18105 (14), GSE21501 (15), GSE17536 (16) and GSE39582 (17) datasets, together with related clinical data from GSE17536 and GSE39582 datasets, from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). GSE18105 cohort includes 111 specimens, including 94 cancer specimens and 17 normal specimens. GSE21501 cohort comprised 148 specimens, including 123 cancer specimens and 25 normal specimens. TCGA-COAD cohort included 448 cancer patients with valid survival data, while GSE17536 and GSE39582 cohorts, adopted as validation cohorts, included 177 and 585 CRC patients with survival data, respectively. Regarding to CRC scRNA-seq data, we evaluated GSE144735 and GSE132465 cohorts from GEO, comprising a total of 91,103 cells (18). A total of 220 LDCD genes were identified by Zou et al. (19). We acquired transcriptomic and matched clinical data of patients undergoing anti-programmed cell death 1 ligand 1 (PD-L1) therapy from IMvigor210 cohort from R package IMvigor210CoreBiologies to survey infiltration of three cell clusters during immunotherapy, including 348 immunotherapy patients, 25 complete response (CR) patients, 43 partial response (PR) patients, 63 stable disease (SD) patients and 167 progressive disease (PD) patients (20). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Processing of scRNAseq data and calculation of LDCD activity

We preprocessed scRNA-seq data via “Seurat” R package (21) and utilized findVariableFeatures function in R package Seurat for scaling data with top 2,000 most variable genes. Variable genes were employed for principal component analysis (PCA), FindNeighbors in Seurat was utilized to obtain nearest neighbors for graph clustering founded upon principal components (PCs), whereas FindCluster in Seurat was adopted to gain cell subtypes, with uniform manifold approximation and projection (UMAP) algorithm utilized for visualizing cells. To eliminate batch impact, batch correction was removed via the Harmony algorithm in Harmony R package, followed by clustering analysis (22), with FindNeighbors and FindCluster in Seurat adopted for obtaining cell types.

For identifying marker genes that were expressed specifically in every cluster, we performed Wilcoxon tests via “FindAllMarkers” and “FindMarkers” algorithms and compared varied cell types. Furthermore, we annotated varied cell types founded upon marker genes from adoptive cell therapy (ACT) database (23). A total of 220 genes were analyzed according to the expression profiles pertaining to LDCD. To evaluate the enrichment scores of the CRC single-cell sequencing data, AUCell algorithms were utilized. Single-cell pseudotime trajectories for tumor cells were built via the “Monocle” R package (24). Furthermore, we surveyed cell-cell communications with union sets of ligand-receptor pairs and single-cell transcriptomes as inputs via CellChat (25).

Generation of differentially expressed genes (DEGs)

We adopted the limma package to conduct differential expression analysis. |Log2 (fold change) FC| >1 and adjusted P<0.05 was regarded as cutoff to screen for DEGs.

Construction of the prognostic signature and calculation of risk scores

We utilized univariate and multivariate Cox regression analyses to screen for prognosis-related candidate genes among the key LDCD genes associated with CRC. Feature selection was carried out via least absolute shrinkage and selection operator (LASSO) regression analysis, with optimal λ value measured through 10-fold cross-validation. We computed each CRC patient’s risk score as per selected genes’ expression level and regression coefficient. Formula is risk score = ∑ⁿ (coef_i × expr_i). coef_i denotes the LASSO Cox i =1 coefficient of gene_i, and expr_i represents gene_i expression.

Construction and validation of predictive nomogram

First, we screened independent prognostic factors via univariate and multivariate Cox regression analyses. After that, we integrated significant factors to establish a nomogram for forecasting patients’ 1-, 5- and 10-year survival rates. Finally, we verified the nomogram’s predictive impact via a calibration curve.

Functional enrichment and gene set level analysis

We employed R package clusterProfiler for conducting DEGs’ Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. We downloaded gene set files of “h all.v7.4.symbols” containing 50 critical gene sets from MSigDB of Broad Institute (http://www.gsea-msigdb.org/) and compared disparities in pathway activity between high- and low-risk assessment groups via gene set enrichment analysis (GSEA) method.

Abundance of infiltrating immune cells in TME

CIBERSORT provides an evaluation of member cell type abundances in a blended cell population, utilizing gene expression data. Expression abundances of 22 sorts of immune cells were estimated in all specimens.

Sensitivity analysis of anticancer drugs

The transcriptome-drug response pairing data of the Genomics of Drug Sensitivity in Cancer (GDSC) database acted as training set, with half-maximal inhibitory concentration (IC50) area under the curve (AUC) values for each compound computed via calcPhenotype function in oncoPredict R package. We subsequently utilized Spearman correlation analysis for evaluating relationships among prognostic characteristics and drug sensitivity (26).

Statistical analysis

We executed all statistical analyses in R version 4.3.2 and performed Kruskal-Wallis and Wilcoxon rank sum tests for analyzing relationships among risk score and clinicopathological parameters and conducted correlation analysis between two groups of variables via Spearman correlation coefficient. Significance of differences between survival curves was measured via log-rank test.

Results

Dissection of the TME founded upon LDCD genes

To elucidate the heterogeneity of TME in CRC and explore LDCD genes’ in CRC patients, we obtained single-cell transcriptome data from 29 CRC patients, of which 88,295 cells were analyzed after quality control and batch correction (Figure 1A). Seven major cell types were revealed by graph-based clustering and canonical cell marker annotation: T cells (n=28,973), which expressed CD3D, CD3E and CD3G; epithelial cells (n=23,563), which were identified by KRT18, CDH1 and EPCAM expression; B cells (n=12,746), marked by CD79A, CD79B and MS4A; fibroblasts (n=10,424), identified by COL1A1, DCN and COL3A1; neutrophils (n=9,070), which were specifically upregulated S100A8, CCL3 and CCL4; endothelial cells (n=3,115), defined by the classical markers VWF, CLDN5 and PECAM1; and mast cells (n=404), positive for KIT, MS4A2 and CPA3 (Figure 1B,1C). To further clarify the tumor cell population, we screened 21,453 malignant epithelial cells via the Cancer-Finder algorithm (27). Next, we computed proportions of these 8 cell types in different patients. Notably, proportions of these 8 cell types differed across patients. These cell populations were uniformly distributed among different patients, indicating that our cell clustering results have good reliability (Figure 1D,1E).

Figure 1 Classification of cell types and gene expression scores pertaining to lysosome-dependent cell death genes in colorectal cancer. (A) UMAP showing cell clusters’ dimensionality reduction visualizations. (B) All cell types’ UMAP plot visualization. (C) Dot plots exhibiting mean expression of marker genes in major cell types. (D) Patient distribution UMAP plot visualization. (E) Bar plots of cell proportions per CRC patient. (F) Heatmap showing lysosome-dependent cell death genes’ expression levels in major cell types. (G) LDCD activity distribution UMAP plot visualization. (H) The histogram shows distribution of high-LDCD-activity cells and low-LDCD-activity cells. (I) Bar plots of the proportions of high-LDCD-activity cells and low-LDCD-activity cells. CRC, colorectal cancer; LDCD, lysosome-dependent cell death; UMAP, uniform manifold approximation and projection.

Next, we assessed the level of LDCD in scRNA-seq data on 220 LDCD genes’ expression (Figure 1F). For example, LAPTM4B was highly expressed in tumor cells [log2FC =1.87, false discovery rate (FDR) <0.001] and endothelial cells (log2FC =1.76, FDR <3.69e−234). Researches have exhibited that LAPTM4B is increased in most cancer types and is associated with drug resistance, survival, cancer cell proliferation, along with unfavorable patient outcomes (28). Compared with that in other cells, RAC2 was dramatically increased in T cells (log2FC =2.63, FDR <0.001) and was associated with T-cell infiltration (29). Using AUCell algorithms, we scored gene sets to assess LDCD activity (Figure 1G). We obtained 15,731 cells with high LDCD activity and 72,564 cells with low LDCD activity (Figure 1H). Additionally, we contrasted distributions of various cell types between high- and low-LDCD-producing cells, revealing that the proportions of neutrophils and fibroblasts were greater in low-LDCD-producing cells and that T, tumor, and B cell distributions were greater in high-LDCD-producing cells (Figure 1I). Next, we performed differential analyses and GO and KEGG enrichment analyses of low-LDCD-activity cells and high-LDCD-activity cells in varied cell types. For example, in tumor cells, genes whose expression was specifically upregulated in high-LDCD-producing cells were enriched in lysosomes, antigen processing and presentation, and ECM-receptor interactions. In neutrophils, the genes whose expression was specifically upregulated in high-LDCD-treated cells were enriched in the terms “lysosome organization” and “glycolipid catabolic process”. These results revealed the underlying impact of LDCD in the CRC TME (Figure S1). Our findings highlight LDCD’s significance for regulating TME at single-cell level, particularly its potential association with antigen presentation system and metabolic reprogramming.

LDCD is involved in tumor cell differentiation and interactions

For exploring LDCD crucial mechanism in tumor cells, tumor cells were further reclustered into 4 subpopulations (tumor cell0-tumor cell3) (Figure 2A). We found that compared with other tumor cell subpopulations, tumor cells had the highest LDCD activity and highest proportion of cell subpopulations with high LDCD activity (Figure 2B-2D). Single-cell trajectory reconstruction through pseudotemporal ordering has become an essential methodology for decoding differentiation cascades, ontogenetic pathways, and spatiotemporal remodeling of immune landscapes in neoplastic microenvironments. In our study, trajectory analysis was used to decipher the cell trajectory of CRC from tumor cell subpopulations. The pseudotime trajectory axis revealed that tumor cells represented the lowest differentiated tumor cell subpopulation along this trajectory (Figure 2E). We also found that cells with high LDCD activity were mainly concentrated in the initial stage of the differentiation trajectory. Furthermore, through RNA-seq analysis, we systematically profiled the transcriptional dynamics of genes associated with LDCD pathogenesis across distinct cellular states in the developmental trajectory of tumor cells (Figure 2F). For example, LAPTM4B was markedly upregulated during terminal tumor cell development, and IL4R was markedly upregulated during the initial stage of tumor cell development. Past researches have demonstrated that IL4R is crucial to early CRC metastasis by activating the ERK pathway (30). These findings indicate that LDCD activity is of close relation with tumor cells’ differentiation status.

Figure 2 Cell developmental trajectory analysis and cell communication. (A) Tumor cell subpopulations’ UMAP plot visualization. (B) Distribution UMAP plot visualization of LDCD activity in tumor cells. (C) Tumor cell subpopulation proportions’ bar plots in high-LDCD-activity cells and low-LDCD-activity cells. (D) Distribution of LDCD activity in tumor cell subpopulations. (E) Tumor cell subpopulations’ cell trajectory and pseudotime analysis. (F) Heatmap depicting genes expression patterns related to LDCD that are differentially expressed during cell development. (G,H) The landscape of cell communication in LDCD high-activity and LDCD low-activity cells. (I) Bar plot displaying specific enrichment of signaling pathways between LDCD high-activity cells and LDCD low-activity cells. (J) Bubble plots displaying the communication probability of specific signaling pathways in LDCD high-activity and LDCD low-activity groups. LDCD, lysosome-dependent cell death; UMAP, uniform manifold approximation and projection.

Cell communication and signaling between tumor cells and other TME cells are important for cancer progression and metastasis. We split the cells into LDCD high- and low-activity groups in accordance with LDCD activity score. We herein conducted a systematic analysis of intercellular communication networks within the CRC TME via the CellChat algorithm to delineate the functional roles of distinct cellular subsets during dynamic remodeling of the TME. By constructing comprehensive cell-cell interaction maps, we elucidated the signaling mechanisms between microenvironmental components and their regulatory effects on tumor ecosystems. These results indicated that the interaction number between cell types was greater in low-LDCD activity group than in high-LDCD activity group, but interaction strength among cell types was superior in high-LDCD activity group than in high-LDCD activity group (Figure 2G,2H). Afterward, we identified specific ligand-receptor pairs that belong to the LDCD high-activity group and the LDCD low-activity group (Figure 2I). Among them, specific signaling pathways (ApoA, NPR1, PVR, CD6, SPP1, RA, ApoE, VCAM, ADGRA, ANNEXIN) were present in the LDCD high-activity group, whereas signaling pathways (MHC-II, PLAU, PSAP, CD40, IL6, BAFF, ACTIVIN, RELN, THY1, CD86, CSF, EPO, NPR2) were specific to the LDCD low-activity group. Furthermore, CellChat analysis revealed specific ligand-receptor interactions among different cell types in LDCD high-activity and LDCD low-activity groups (Figure 2J). In LDCD high-activity group, the ANNEXIN pathway was activated between tumor cells and neutrophils via ligand-receptor pairs (ANXA1-FPR1 and ANXA1-FPR3). Studies have shown that ANXA1 regulates cancer cell invasion by binding to the receptors FPR1 and FPR3 (31). The ApoA signaling pathway (APOA1-ABCA1, APOA1-TREM2) between tumor cells and neutrophils was specifically enriched in the LDCD high-activity group. APOA1 promotes tumor growth and worsens CRC prognosis by regulating ABCA1 expression (32). Furthermore, immune-associated ligand receptor pairs, including the MHC-II signaling pathway (HLA-DMA-CD4, HLA-DRA-CD4, and HLA-DRB1-CD4), the CD86 signaling pathway (CD86-CD28 and CD86-CTLA4) and the CD40LG (ITGA5 + ITGB1), CD40LG (ITGAM + ITGB2), and CD40LG-CD40, were enriched in LDCD low-activity group. Such outcomes illustrate heterogeneity of the ligand-receptor pairs we identified across different LDCD activity groups and reveal the latent control networks orchestrating cellular cross-talk in the tumor niche.

Identification of key LDCD genes related to CRC and establishment of prognostic signatures

For identifying key genes pertaining to CRC, we integrated and analyzed four bulk RNA-seq datasets from TCGA database (TCGA-COAD and TCGA-READ) and GEO database (GSE18105 and GSE21501) between cancer and normal specimens (Figure S2A). Differential expression analysis revealed 886 increased genes and 1,191 decreased genes in TCGA-COAD cohort. In TCGA-READ cohort, 944 genes were increased, and 1,405 genes were decreased. GSE18105 dataset contained 1,679 upregulated genes and 1,078 downregulated genes. The GSE21510 dataset revealed 1,577 genes whose expression was upregulated and 1,352 genes whose expression was downregulated (Figure S2B). Through multicohort cross-analysis, we identified 617 codifferentially expressed genes, like 220 increased and 397 decreased genes, which were defined as key CRC-related genes. GO and KEGG analyses demonstrated that key CRC-related genes were pertaining to primary alcohol metabolic processes, the response to chemokines, nitrogen metabolism, the stress response to metal ions and other functions (Figure S2C). Further analysis implied that 1,765 DEGs were defined as LDCD-related tumor genes between tumor cells with high LDCD activity and those with low LDCD activity. By intersecting LDCD-related tumor genes with key CRC-related genes, we identified 95 common genes that act as key LDCD genes associated with CRC (Figure 3A). Notably, such genes are involved mainly in biological processes (BPs), like chemokine signal transduction and metabolic regulation (Figure 3B).

Figure 3 Calculation of risk scores pertaining to LDCD and prognostic signature development. (A) Venn diagram displaying 95 intersection genes between key CRC-related genes and LDCD-related tumor genes. (B) 95 intersecting genes’ GO and KEGG enrichment analyses. (C) High- and low-risk score groups Kaplan-Meier curves in cohorts. (D) ROC curve in training set. (E) OS of patients in high- and low-risk score groups in GSE17536 and GSE39582 cohorts. (F) Cluster heatmap of key LDCD genes related to CRC. (G) HIGD1A and PARM1 expression levels according to clinicopathological features (pathological T, N, M and tumor stages). (H) Risk score according to clinicopathological characteristics. AUC, area under the curve; CRC, colorectal cancer; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LDCD, lysosome-dependent cell death; M, metastasis; N, node; OS, overall survival; ROC, receiver operating characteristic; T, tumor.

To screen prognostic genes pertaining to CRC progression, we firstly conducted univariate Cox regression analysis about 95 key LDCD genes associated with CRC and identified 8 key LDCD genes related to CRC as significantly related to patient survival (P<0.05). We conducted further multivariate Cox regression analysis in combination with clinicopathological features [like pathological tumor (T), node (N), metastasis (M) and tumor stages], and four key LDCD genes associated with CRC had independent prognostic value (P<0.05). By employing LASSO-regularized Cox proportional hazards modeling of four gene candidates, we derived a multiparametric prognostic algorithm. Through rigorous algorithmic filtering, HIGD1A and PARM1 emerged as cardinal biomarkers with predictive supremacy in survival stratification. We adopted formula as below to calculate risk score: HIGD1A expression * (−0.1276) + PARM1 expression * (−0.0903). We stratified patients into high- and low-risk cohorts as per median risk thresholds. Survival probability curves generated through the Kaplan-Meier method demonstrated markedly poorer survival outcomes in high-risk stratum than in high-risk stratum (log-rank test, P=0.02; Figure 3C). Furthermore, time-dependent receiver operating characteristic (ROC) analysis implied that such a risk score has excellent predictive efficacy in TCGA-COAD cohort, and AUC values of total survival at 1, 5, and 10 years all exceeded 0.6 (Figure 3D). Two independent validation cohorts (GSE17536, n=177; GSE39582, n=579) were used to evaluate the robustness of the signature by conducting validation analyses. Mirroring TCGA training set findings, log-rank testing confirmed marked divergence in survival trajectories across risk-stratified demarcations. The validation cohort consistently demonstrated an attenuated survival probability in elevated-risk populations compared with their baseline counterparts, confirming the cross-dataset predictive efficacy of this prognostic signature (GSE17536, log-rank test, P=0.025; GSE39582, log-rank test, P=0.01; Figure 3E). Further analysis implied that core genes PARM1 and HIGD1A were dramatically downregulated in high-risk score group (Figure 3F). PARM1 expression was negatively related to tumor progression. PARM1 expression in patients with stage III/IV disease was dramatically lower than that in patients with stage I/II disease (P=0.02; Figure 3G), and it was also significantly higher than that in patients with stage N0 disease than in patients with stage N1/N2 disease (P=0.02; Figure 3G). Furthermore, a high-risk score strongly depends on advanced clinical features. Risk score of patients with stage III/IV disease was greater than that of patients with stage I/II disease (P=0.02), and that of patients with stage N1/N2 disease was greater than that of patients with stage N0 disease (P=0.03, Figure 3H). These findings imply that this prognostic feature can not only effectively distinguish the survival risk of patients but also that its score changes strongly pertain to tumor malignant progression and has potential clinical value.

Independent predictive value of clinical models’ prognostic signature and construction for determining whether prognostic signature is a CRC independent factor, through univariate and multivariate Cox regression analyses, including prognostic signatures and clinicopathological factors like TNM stage, we confirmed that prognostic signature was a CRC independent prognostic factor (Figure 4A). This finding was further verified in two independent validation cohorts (GSE17536 and GSE39582), and univariate analysis implied that prognostic signature had significant independent predictive value (P<0.05, Figure 4B,4C). Time-dependent ROC analysis implied that model maintained stable efficacy in forecasting 10-year survival (AUC >0.6, Figure S2D). Based on the above findings, we integrated the prognostic signature, age and pathological stage to establish a clinically practical prognostic nomogram (Figure 4D). This model visually forecasts patients’ 1-, 5- and 10-year survival probabilities via quantifying degree of contribution of each risk factor. Calibration curves demonstrated powerful concordance between projected survival probabilities and empirical clinical results, validating algorithm’s predictive ability (Figure 4E). This comprehensive assessment system can predict CRC patient prognosis more accurately.

Figure 4 Independent prognostic analysis of risk scores and clinicopathologic factors. (A) Univariate Cox analysis and multivariate Cox analysis forest plot in TCGA-COAD cohort. (B) Univariate Cox analysis and multivariate Cox analysis results forest plot for GSE17536 cohort. (C) Univariate Cox analysis and multivariate Cox analysis results forest plot for GSE39582 cohort. (D) Nomogram founded upon independent prognostic factors. (E) Calibration curves for nomogram. **, P<0.01; ***, P<0.001. CI, confidence interval; HR, hazard ratio; LDCD, lysosome-dependent cell death; M, metastasis; N, node; OS, overall survival; T, tumor; TCGA, The Cancer Genome Atlas.

Function and variant landscape of prognostic signature in CRC patients

In TCGA-COAD cohort investigation, we split CRC specimens into high- and low-risk cohorts, and median risk score was used as threshold. Subsequent differential expression analysis identified 188 dramatically altered genes distinguishing these two prognostic groups (Figure 5A,5B). We subsequently conducted GSEA in two risk score groups. We found that epithelial mesenchymal transition, myogenesis, apical junction, hedgehog signaling and cancer progression-related functions were dramatically enriched in high-risk score group (Figure 5C). For instance, dysregulation of hedgehog signaling is associated with developmental abnormalities and colon cancer (33). Low-risk score group was dramatically enriched in metabolic and cell cycle functions, like fatty acid metabolism, glycolysis, oxidative phosphorylation, and the E2F and MYC targets (Figure 5D). Given the pervasive mutational landscape of CRC, we conducted comparative profiling of the 20 most frequently altered genes across defined risk strata, delineating distinct mutation patterns between high- and low-risk cohorts (Figure 5E). TP53, RAR1 and DST mutation frequencies were greater in high-risk score groups than in low-risk score groups. In accordance with our findings, TP53 mutations in CRC generally promote cancer cell dryness, cell proliferation, invasion, and metastasis, thus contributing to cancer progression (34). In comparison with those in high-risk group, KRAS, CSMD1 and PCLO mutation frequencies were increased in low-risk group. KRAS mutations lead to KRAS protein activation, stimulating cell proliferation and survival and leading to tumorigenesis (35). These results collectively indicate that different risk score groups have unique molecular characteristics and biological behaviors.

Figure 5 Biological characteristics of various risk score groups in TCGA cohort. (A) Heatmap of DEGs between high- and low-risk score groups in TCGA-COAD cohort. (B) Volcano plot of DEGs between high- and low-risk score groups in TCGA-COAD cohort. (C) GSEA outcomes in high-risk score group. (D) GSEA outcomes in low-risk score group. (E) Tumor somatic mutation distribution waterfall plot of the top 20 highly mutated genes in high- and low-risk score groups. DEGs, differentially expressed genes; GSEA, gene set enrichment analysis; TCGA, The Cancer Genome Atlas.

Prognostic signature as a predictive biomarker for immunotherapy in CRC

To elucidate mechanistic connections among risk stratification frameworks, lymphocyte infiltration dynamics, and tumor-immune ecosystem remodeling, we applied the CIBERSORT computational deconvolution platform for systematic immune profiling. These findings suggest considerable disparities between high- and low-risk score groups regarding to M0 macrophages, M2 macrophages, naive B, memory B, plasma, CD8⁺ T, resting CD4⁺ memory T, activated CD4⁺ memory T, regulatory T and gamma delta T cells, resting dendrites, and activated dendritic and resting mast cells (Figure 6A). Moreover, Spearman correlation analysis revealed a strong relationship between the risk score and such three differential immune cell types (Figure 6B). Furthermore, immune checkpoint genes (CD40LG, LAIR1, and TIMGD2) were strongly different between the two groups, with CD40LG, LAIR1, and TIMGD2 strongly related to the risk score (Figure 6C,6D). Next, we conducted correlation analysis of survival time and the prognostic signature in bladder cancer patients dealt utilizing anti-PD-L1 therapy via IMvigor210 dataset (20). We found that high-risk patients exhibited impaired OS responding to PD-L1 antibody treatment (log-rank test, P<0.001; Figure 6E). Anti-PD-L1-treated patients with high-risk scores were of shorter survival times. Importantly, patients in high-risk score group had more PD and SD than patients in low-risk score group did (Figure 6F), meaning that patients with high-risk scores in the prognostic signature had an especially inferior response to anti-PD-L1 therapy than the other patients.

Figure 6 Prediction of immunotherapy efficacy and drug sensitivity effects. (A) Box plot of immune cell percentages in high- and low-risk score groups. (B) Heatmap of the relationships among immune cells, the prognostic signature, and model genes. (C) Box plot of immune checkpoints in high- and low-risk score groups. (D) Heatmap of the relationships among immune checkpoints, the prognostic signature, and model genes. (E) Survival analyses of high- and low-risk score groups in anti-PD-L1 immunotherapy cohort utilizing Kaplan-Meier curves (IMvigor210 cohort). (F) Bar plots depicting that patients in high- and low-risk score groups are more probably to respond to anti-PD-L1 treatment. (G) Heatmap of relationships among drugs, prognostic signatures, and model genes. (H) IC50 value comparison of drugs between high- and low-risk score groups and correlation between IC50 values and prognostic signature in TCGA-COAD cohort. ns, P>0.05; *, P<0.05; **, P<0.01; ***, P<0.001. CR, complete response; IC50, half-maximal inhibitory concentration; LDCD, lysosome-dependent cell death; PD, progressive disease; PD-L1, programmed cell death 1 ligand 1; PR, partial response; SD, stable disease; TCGA, The Cancer Genome Atlas.

To systematically investigate therapeutic correlations between molecular risk stratification and pharmacogenomic vulnerabilities, we conducted quantitative pharmacodynamic profiling via IC50 metrics across colorectal carcinoma specimens, identifying differential chemotherapeutic response patterns stratified by the risk cohort. The landscape of the relationship between drug sensitivity and the prognostic signature is shown in Figure 6G. Among them, six drugs had the strongest association with the prognostic signature (|R|≥0.35, Figure 6H). We found that the IC50 values of panobinostat, dacinostat, erlotinib, navitoclax, the JNK inhibitor VIII and WZ-1-84 were dramatically greater in high-risk score group. Researches have exhibited that panobinostat promotes autophagy and leads to CRC cell death (36). Combining bevacizumab and erlotinib in metastatic CRC maintenance treatment apparently enhances patient survival (37). Therefore, we identified multiple potential drugs, which can target high-risk score group of the prognostic signature.

Discussion

Regulatory cell death (RCD), a precisely regulated cell death form, has dual functions in tumor biology. This programmed death process is involved in maintaining tissue homeostasis, and its dysregulation has emerged as a key feature of tumorigenesis and development, providing a vital target for the development of novel anticancer strategies (38). Lysosomes, as important degradation centers within cells, regulate various physiological and pathological processes through the hydrolase system in which they are rich (39). In the TME, lysosomes affect cell fate through dual mechanisms. On the one hand, lysosomes take a cellular protective effect via autophagy pathway to help tumor cells cope with metabolic stress (40). On the other hand, when lysosomal function is impaired, cell death can be induced through mechanisms like disrupting autophagosome-lysosomal fusion and disrupting autophagic flux (41); this makes targeted LDCD a highly promising new direction for anticancer treatment. Single-cell sequencing technology offers a revolutionary research perspective for CRC studies. Analysis of TME heterogeneity characteristics at single-cell resolution significantly advances understanding of disease occurrence mechanisms. Research has revealed that specific connective tissue formation characteristics in the TME may be closely related to immunotherapy resistance (42). More importantly, in-depth analysis based on single-cell transcriptomics enables us to precisely characterize patients according to the molecular characteristics of the TME and discover that different subtypes of tumor cells adopt unique immune escape strategies (43). These breakthrough discoveries have improved our insight into CRC pathogenesis and offer a major theoretical foundation for proposing targeted treatment regimens and realizing individualized medicine.

Using single-cell transcriptome sequencing technology, our study systematically analyzed the expression profiles of LDCD-related genes in CRC. Research has shown that LDCD-related genes exhibit significantly heterogeneous expression characteristics in different cell subsets. By establishing the LDCD activity assessment system, cells in the TME were split into high/low LDCD activity subgroups, and highly active cells were enriched mainly at the initial stage of the tumor differentiation trajectory. Further analysis of cellular interactions revealed significant remodeling of communication networks between high-LDCD-active cells and TME components, particularly cancer-associated fibroblasts (CAFs) and neutrophils. Interaction among cancer cells and CAFs affects tumor growth and invasion (44). These findings indicate that tumor cells with different LDCD activities have unique communication patterns, offering a novel perspective for understanding regulatory part of LDCD in CRC occurrence and development.

A novel risk-scoring system utilizing 2 key LDCD genes pertaining to CRC (HIGD1A and PARM1) was subsequently created to predict therapy responses and stratify CRC patients from TCGA and GEO datasets into different risk groups. We split patients into high- and low-risk score groups according to gene expression patterns, with those in high-risk score group exhibiting poorer prognoses. Researches have implied that HIGD1A overexpression inhibits CRC cell proliferation, invasion and migration (45). PARM1 low expression in CRC tumor tissues might be a promising novel prognostic biomarker for CRC (46). These key LDCD genes associated with CRC are closely related to cancer development and regulation of TME. Through multicohort validation, we confirmed that prognostic signature has stable prognostic predictive efficacy and could be adopted as a CRC-independent prognostic marker. Significant biological differences were shown between high-risk and low-risk score groups. The former was enriched in tumor progression-related pathways, whereas the latter was characterized by metabolic reprogramming and cell cycle regulation. Importantly, among patients treated utilizing anti-PD-L1 therapy, those in high-risk score group had a poorer prognosis. Such outcomes offer a major foundation for optimizing clinical immunotherapy regimens.

We incorporated scRNA-seq and bulk RNA-seq data for analyzing LDCD-related gene expression characteristics systematically in CRC. Research has revealed not only cell type-specific expression patterns of such genes but also their critical roles in the regulation of tumor differentiation and intercellular communication. Therefore, we established a prognostic risk score signature with clinical predictive value, which was significantly correlated with treatment response, genomic variation and immune microenvironment characteristics. Although this research provides a novel perspective for assessing prognosis and regulating CRC microenvironment, there are still limitations. The spatial distribution characteristics of gene expression in the TME have not yet been clarified. The scRNA-seq trajectory inference and CellChat analysis are computational predictions based on transcriptomic snapshots. They reflect correlational patterns rather than definitive lineage tracing or functional protein interactions. Tissue dissociation may introduce selection bias, while trajectory analysis infers a developmental trend, future studies utilizing spatial transcriptomics (ST) are warranted to validate the precise in situ localization of high-LDCD cells and rule out dissociation artifacts. Experimental validation (e.g., lineage tracing in mice) is required to confirm these causal relationships and knockout experiment is performed to validate that the identified key genes modulate CRC progression through the LDCD mechanism (47). Larger-scale clinical cohorts are needed to verify our model’s diagnostic efficacy. We should focus on clarifying the molecular mechanism and regulatory network of LDCD in CRC in the future, which will provide important clues for better understanding its occurrence and development.

Conclusions

Through integrating scRNA-seq and bulk RNA-seq, a prognostic signature for CRC founded upon LDCD was successfully established. This prognostic signature not only demonstrates stable clinical predictive efficacy but also reveals the regulatory role of LDCD in tumor progression. Such outcomes offer important tools for clinical prognosis judgment and lay theoretical basis for novel therapeutic strategies targeting LDCD development.

Acknowledgments

All authors thank all colleagues and collaborators contributing to our work. We also acknowledge the GEO and UCSC databases for their datasets, which made our work possible.

Footnote

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

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

Funding: None.

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

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

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

References

1. Dekker E, Tanis PJ, Vleugels JLA, et al. Colorectal cancer. Lancet 2019;394:1467-80. [Crossref] [PubMed]
2. Ionescu VA, Gheorghe G, Bacalbasa N, et al. Colorectal Cancer: From Risk Factors to Oncogenesis. Medicina (Kaunas) 2023;59:1646. [Crossref] [PubMed]
3. Yamamoto D, Oshima H, Wang D, et al. Characterization of RNF43 frameshift mutations that drive Wnt ligand- and R-spondin-dependent colon cancer. J Pathol 2022;257:39-52. [Crossref] [PubMed]
4. Pabla B, Bissonnette M, Konda VJ. Colon cancer and the epidermal growth factor receptor: Current treatment paradigms, the importance of diet, and the role of chemoprevention. World J Clin Oncol 2015;6:133-41. [Crossref] [PubMed]
5. Wang Y, Zhang Z, Sun W, et al. Ferroptosis in colorectal cancer: Potential mechanisms and effective therapeutic targets. Biomed Pharmacother 2022;153:113524. [Crossref] [PubMed]
6. Li R, Wu Y, Li Y, et al. Targeted regulated cell death with small molecule compounds in colorectal cancer: Current perspectives of targeted therapy and molecular mechanisms. Eur J Med Chem 2024;265:116040. [Crossref] [PubMed]
7. Gemingnuer A, Yin R, Liu Y, et al. Lysosome-dependent cell death: disease implications and potential therapeutic targets. Mol Biol Rep 2025;52:881. [Crossref] [PubMed]
8. Zhang Z, Yue P, Lu T, et al. Role of lysosomes in physiological activities, diseases, and therapy. J Hematol Oncol 2021;14:79. [Crossref] [PubMed]
9. Yan R, Zhang P, Shen S, et al. Carnosine regulation of intracellular pH homeostasis promotes lysosome-dependent tumor immunoevasion. Nat Immunol 2024;25:483-95. [Crossref] [PubMed]
10. Piao S, Amaravadi RK. Targeting the lysosome in cancer. Ann N Y Acad Sci 2016;1371:45-54. [Crossref] [PubMed]
11. Vandenabeele P, Galluzzi L, Vanden Berghe T, et al. Molecular mechanisms of necroptosis: an ordered cellular explosion. Nat Rev Mol Cell Biol 2010;11:700-14. [Crossref] [PubMed]
12. Guo Z, He L, Wang W, et al. FUT2-dependent fucosylation of LAMP1 promotes the apoptosis of colorectal cancer cells by regulating the autophagy-lysosomal pathway. Cancer Lett 2025;619:217643. [Crossref] [PubMed]
13. Wang X, An F, Wang B, et al. From "lysosomal addiction" to targeted therapies: exploiting novel windows in colorectal cancer. Eur J Med Res 2025;30:957. [Crossref] [PubMed]
14. Matsuyama T, Ishikawa T, Mogushi K, et al. MUC12 mRNA expression is an independent marker of prognosis in stage II and stage III colorectal cancer. Int J Cancer 2010;127:2292-9. [Crossref] [PubMed]
15. Stratford JK, Yan F, Hill RA, et al. Genetic and pharmacological inhibition of TTK impairs pancreatic cancer cell line growth by inducing lethal chromosomal instability. PLoS One 2017;12:e0174863. [Crossref] [PubMed]
16. Chen MS, Lo YH, Chen X, et al. Growth Factor-Independent 1 Is a Tumor Suppressor Gene in Colorectal Cancer. Mol Cancer Res 2019;17:697-708. [Crossref] [PubMed]
17. Marisa L, de Reyniès A, Duval A, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med 2013;10:e1001453. [Crossref] [PubMed]
18. Lee HO, Hong Y, Etlioglu HE, et al. Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. Nat Genet 2020;52:594-603. [Crossref] [PubMed]
19. Zou Y, Xie J, Zheng S, et al. Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery. Int J Surg 2022;107:106936. [Crossref] [PubMed]
20. Mariathasan S, Turley SJ, Nickles D, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 2018;554:544-8. [Crossref] [PubMed]
21. Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell 2021;184:3573-3587.e29. [Crossref] [PubMed]
22. Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods 2019;16:1289-96. [Crossref] [PubMed]
23. Quan F, Liang X, Cheng M, et al. Annotation of cell types (ACT): a convenient web server for cell type annotation. Genome Med 2023;15:91. [Crossref] [PubMed]
24. Borcherding N, Voigt AP, Liu V, et al. Single-Cell Profiling of Cutaneous T-Cell Lymphoma Reveals Underlying Heterogeneity Associated with Disease Progression. Clin Cancer Res 2019;25:2996-3005. [Crossref] [PubMed]
25. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
26. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
27. Zhong Z, Hou J, Yao Z, et al. Domain generalization enables general cancer cell annotation in single-cell and spatial transcriptomics. Nat Commun 2024;15:1929. [Crossref] [PubMed]
28. Ji X, Ma H, Du Y. Role and mechanism of action of LAPTM4B in EGFR-mediated autophagy. Oncol Lett 2022;23:109. [Crossref] [PubMed]
29. Xu Y, Cai Y, Deng Y, et al. RAC2 as a Tumor-Suppressive Biomarker Associated with T Cell Infiltration in Breast Cancer. Cancer Biother Radiopharm 2025;40:62-77. [Crossref] [PubMed]
30. Jin Q, Yang H, Jing Z, et al. IL4/IL4R signaling promotes the osteolysis in metastatic bone of CRC through regulating the proliferation of osteoclast precursors. Mol Med 2021;27:152. [Crossref] [PubMed]
31. Cheng TY, Wu MS, Lin JT, et al. Annexin A1 is associated with gastric cancer survival and promotes gastric cancer cell invasiveness through the formyl peptide receptor/extracellular signal-regulated kinase/integrin beta-1-binding protein 1 pathway. Cancer 2012;118:5757-67. [Crossref] [PubMed]
32. Aguirre-Portolés C, Feliu J, Reglero G, et al. ABCA1 overexpression worsens colorectal cancer prognosis by facilitating tumour growth and caveolin-1-dependent invasiveness, and these effects can be ameliorated using the BET inhibitor apabetalone. Mol Oncol 2018;12:1735-52. [Crossref] [PubMed]
33. Skoda AM, Simovic D, Karin V, et al. The role of the Hedgehog signaling pathway in cancer: A comprehensive review. Bosn J Basic Med Sci 2018;18:8-20. [Crossref] [PubMed]
34. Liebl MC, Hofmann TG. The Role of p53 Signaling in Colorectal Cancer. Cancers (Basel) 2021;13:2125. [Crossref] [PubMed]
35. Zhu G, Pei L, Xia H, et al. Role of oncogenic KRAS in the prognosis, diagnosis and treatment of colorectal cancer. Mol Cancer 2021;20:143. [Crossref] [PubMed]
36. Zamani M, Safari F, Siri M, et al. Epigenetic modulation of autophagy pathway by small molecules in colorectal cancer: a systematic review. J Cancer Res Clin Oncol 2024;150:474. [Crossref] [PubMed]
37. Kaveh S, Ebrahimi P, Rezapour A, et al. Bevacizumab and erlotinib versus bevacizumab for colorectal cancer treatment: systematic review and meta-analysis. Int J Clin Pharm 2019;41:30-41. [Crossref] [PubMed]
38. Peng F, Liao M, Qin R, et al. Regulated cell death (RCD) in cancer: key pathways and targeted therapies. Signal Transduct Target Ther 2022;7:286. [Crossref] [PubMed]
39. Berg AL, Rowson-Hodel A, Wheeler MR, et al. Engaging the Lysosome and Lysosome-Dependent Cell Death in Cancer. In: Mayrovitz HN. editor. Breast Cancer. Brisbane: Exon Publications; 2022.
40. Galluzzi L, Vitale I, Aaronson SA, et al. Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death Differ 2018;25:486-541. [Crossref] [PubMed]
41. Sabatini DM. mTOR and cancer: insights into a complex relationship. Nat Rev Cancer 2006;6:729-34. [Crossref] [PubMed]
42. Qi J, Sun H, Zhang Y, et al. Single-cell and spatial analysis reveal interaction of FAP(+) fibroblasts and SPP1(+) macrophages in colorectal cancer. Nat Commun 2022;13:1742. [Crossref] [PubMed]
43. Chu X, Li X, Zhang Y, et al. Integrative single-cell analysis of human colorectal cancer reveals patient stratification with distinct immune evasion mechanisms. Nat Cancer 2024;5:1409-26. [Crossref] [PubMed]
44. Beatty GL, Chiorean EG, Fishman MP, et al. CD40 agonists alter tumor stroma and show efficacy against pancreatic carcinoma in mice and humans. Science 2011;331:1612-6. [Crossref] [PubMed]
45. Lv JL, Tan YJ, Ren YS, et al. Procyanidin C1 inhibits tumor growth and metastasis in colon cancer via modulating miR-501-3p/HIGD1A axis. J Adv Res 2024;60:215-31. [Crossref] [PubMed]
46. Liu YR, Hu Y, Zeng Y, et al. Neurexophilin and PC-esterase domain family member 4 (NXPE4) and prostate androgen-regulated mucin-like protein 1 (PARM1) as prognostic biomarkers for colorectal cancer. J Cell Biochem 2019;120:18041-52. [Crossref] [PubMed]
47. Liu Z, Cai L, Liao B, et al. Comprehensive characterization of lysosome-dependent cell death reveals prognostic significance and immune landscape in colon adenocarcinoma. Sci Rep 2025;15:33132. [Crossref] [PubMed]