Prognostic implications of a programmed cell death-related long non-coding RNA signature and its relevance to immune features in esophageal carcinoma

Introduction

Esophageal carcinoma (ESCA) is a prevalent global malignancy, occupying the seventh and sixth positions in incidence and mortality, respectively. Disparities in demographic and histologic characteristics, such as sex ratio, geography, and tumor subtype, contribute to its epidemiologic diversity (1). Although significant progress has been made with surgery, chemoradiotherapy, and targeted approaches, prognosis remains unfavorable, with most patients not surviving beyond 5 years (2,3). This may be attributed to late-stage diagnoses, treatment resistance, and complex tumor biology. Furthermore, exposure to key risk factors such as alcohol and tobacco for squamous carcinoma, and reflux-related inflammation for adenocarcinoma, may promote oncogenic pathways that limit therapeutic efficacy and worsen survival (4,5). Furthermore, although programmed cell death protein 1 (PD-1)/programmed death-ligand 1 (PD-L1) blockade offers benefit to some, its overall efficacy remains limited, with many patients showing insufficient responses (6). These challenges emphasize the necessity of developing reliable biomarkers and innovative therapeutic targets to enhance patient-specific management (7,8).

The role of programmed cell death (PCD) in cancer is under increasing investigation due to its potential effects on how tumors develop, progress, and resist to therapy (9). It comprises a diverse array of regulated mechanisms, such as apoptosis, necroptosis, pyroptosis, ferroptosis, entosis, NETosis, parthanatos, lysosome-dependent death, autophagy-related cell death, alkaliptosis, and oxeiptosis (10). Among these, apoptosis is one of the most extensively studied forms, characterized by caspase-mediated DNA fragmentation and generally associated with tumor suppression (11,12). Ferroptosis, driven by oxidative damage in the presence of iron, has been associated with cellular stress responses and may influence therapy outcomes (13). Under certain conditions, autophagy can shift from a survival mechanism to a self-destructive process, potentially resulting in cell death (14). Pyroptosis, a highly inflammatory form of PCD mediated by gasdermin proteins, has been implicated in both antitumor immunity and immunopathology (15,16).

Long non-coding RNA (lncRNA), which exceeds 200 nucleotides in length, can influence gene expression patterns under both normal and disease conditions through a range of mechanisms—including interaction with chromatin, RNAs, and proteins—and are increasingly recognized for their regulatory functions in gene expression (17). Emerging evidence suggests that lncRNA is involved in modulating distinct cell death mechanisms beyond apoptosis, such as ferroptosis, pyroptosis, and necroptosis (18). For instance, LINC00336 has been shown to inhibit ferroptosis in lung cancer via modulation of iron metabolism and microRNA (miRNA) activity (19), while MEG3 is downregulated in esophageal squamous cell carcinoma (ESCC) via endoplasmic reticulum stress activation (20). Despite these findings in several cancer types, studies focusing on PCD-related lncRNAs in ESCA remain limited.

In this investigation, using data from The Cancer Genome Atlas (TCGA) database, we constructed a prognostic signature for patients with ESCA using eight PCD-related lncRNAs. Univariate Cox regression analysis was first applied to screen for prognostic lncRNAs. Subsequently, the least absolute shrinkage and selection operator (LASSO) regression followed by multivariate Cox regression was conducted to establish the prognostic value of the signature by comparing overall survival (OS) between distinct risk groups. To further elucidate the biological significance of the signature, we performed weighted gene co-expression network analysis (WGCNA), gene set enrichment analysis (GSEA), and Gene Ontology (GO) enrichment analysis. These analyses indicated that high-risk patients were associated with suppressed immune-related pathways and decreased infiltration of immune cell populations. Additionally, drug sensitivity prediction suggested differential response potential between the risk groups.

This investigation introduces a novel lncRNA-based prognostic signature associated with PCD in ESCA, which may aid in patient risk stratification and offer theoretical support for individualized treatment planning and future mechanistic research. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1598/rc).

Methods

Data acquisition

RNA-sequencing (RNA-seq) data and corresponding clinical information for 174 esophageal tissue samples were obtained from TCGA database (https://www.cancer.gov/ccg/research/genome-sequencing/tcga). Among these, 162 tumor samples were derived from patients diagnosed with ESCC or esophageal adenocarcinoma (EAC), along with 12 adjacent non-tumor samples. Clinical data related to ESCA patients were also retrieved from the TCGA database. All datasets underwent standardized preprocessing, including normalization and log₂ transformation, performed using the DESeq2 R package (https://bioconductor.org/packages/DESeq2/) to ensure data consistency and quality. The overall study workflow is illustrated in Figure 1. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Figure 1 The flowchart of the whole study. ESCA, esophageal carcinoma; GSEA, gene set enrichment analysis; LASSO, least absolute shrinkage and selection operator; TCGA, The Cancer Genome Atlas; WGCNA, weighted gene co-expression network analysis.

Identification of PCD-related lncRNAs

To identify lncRNAs associated with PCD, we compiled a set of genes associated with PCD based on data from Zou et al.’s previous study (21). We then performed Pearson correlation analysis on the annotated lncRNAs expression profiles and these genes associated with PCD. Criteria for screening relevant lncRNAs include absolute correlation coefficient (|Cor|) exceeds 0.6 and the P value is below 0.05.

Construction and validation of a prognostic signature based on lncRNAs associated with PCD

Using univariate Cox regression analysis, combined with the LASSO method implemented through R package glmnet (https://CRAN.R-project.org/package=glmnet), we identified the lncRNAs associated with PCD-related that are associated with OS. The TCGA esophageal cancer data were randomly divided into a training set and a validation set, each containing an equal number of samples. A risk score was then calculated for each patient by multiplying the expression level of each selected lncRNA with its corresponding multivariate Cox regression coefficient and summing the results. The formula used was: risk score = ∑ (expression level of lncRNA × its multivariate Cox regression coefficient). Patients were subsequently categorized into high-risk and low-risk groups using the median risk score from the training set as the threshold. Differences in survival between groups were analyzed using the Kaplan-Meier method, and significance was assessed by a log-rank test. The signature’s ability to predict survival was evaluated through receiver operating characteristic (ROC) curve analysis using the pROC R package (22), with a focus on 1-, 2-, and 3-year survival outcomes. Signature performance was quantified using indicators such as the area under the curve (AUC), sensitivity, and specificity.

Association between prognostic signature constructed from PCD-related lncRNAs and clinicopathological features of patients

Using the ggplot2 (https://CRAN.R-project.org/package=ggplot2), ggpubr (https://CRAN.R-project.org/package=ggpubr), and pheatmap (https://CRAN.R-project.org/package=pheatmap) libraries available from CRAN in R software, the relationship between PCD-related lncRNAs characteristics and multiple clinical pathology characteristics was visualized. Survival in different patient subgroups was assessed using stratified OS analysis, conducted with the “survival” (https://CRAN.R-project.org/package=survival) and “survminer” (https://CRAN.R-project.org/package=survminer) R packages, which support univariate and multivariate Cox regression modeling to identify independent prognostic indicators.

Tumor mutation burden (TMB) analysis

To analyze mutation patterns and TMB scores in the high- and low-risk groups, we employed the R software packages maftools (23) and ggplot2, along with the forestPlot package (https://CRAN.R-project.org/package=forestplot).

Nomogram construction

To develop a prediction tool, a nomogram was built in the R statistical environment using the “rms” package (available at https://CRAN.R-project.org/package=rms). This nomogram integrates key prognostic variables, including age, clinical tumor (T), node (N), and metastasis (M) stages, as well as the derived risk scores. The reliability of the signature was then evaluated using a calibration plot comparing predicted survival probabilities with actual outcomes. The discriminative power and predictive accuracy of the feature were assessed by calculating the concordance index (C-index) and AUC, using the riskRegression package (available at https://CRAN.R-project.org/package=riskRegression) and the timeROC package (24).

WGCNA

A weighted gene co-expression network was constructed using the WGCNA R package (25) to identify gene modules associated with the PCD-related lncRNAs signature. In order to maintain a scale-free network in the TCGA-ESCA dataset, a soft threshold power of 12 was selected. Genes with similar expression patterns were grouped into modules, and the module showing the strongest correlation with PCD-related lncRNAs characteristics was selected for subsequent functional enrichment analysis.

Functional enrichment analyses

In order to understand the biological significance of the gene module associated with this trait, we conducted a series of enrichment analyses, including GO, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and hallmark pathway enrichment analyses. These analyses were performed using the clusterProfiler, org.Hs.eg.db (https://bioconnector.org/packages/org.Hs.eg.db/), and enrichplot (https://bioconnector.org/packages/enrichplot/) R packages.

GSEA

GSEA was performed to detect hallmark pathways that were preferentially active in high-risk and low-risk groups. Genes were ranked according to their differential expression between groups, and the “clusterProfiler” R package (26) was used to analyze pathway enrichment. Pathways with adjusted P values below 0.05 were considered significantly enriched.

Immune infiltration analysis

To assess immune infiltration, the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) R package (27) was used to calculate immune scores, stromal scores, and ESTIMATE composite scores, providing an in-depth understanding of the tumor immune microenvironment. To quantify the infiltration levels of 28 immune cell types in tumor tissues, we applied single-sample GSEA (ssGSEA) using the GSVA package in R (28). Associations between immune cell infiltration and risk scores were analyzed using both Spearman’s rank and Pearson’s correlation methods.

Drug sensitivity

Using the oncoPredict package (29) to identify potential chemotherapeutic drugs, using the GDSC2 dataset containing 198 drug data from the oncoRespond project hosted on Open Science Framework (OSF) (available at https://www.cancerrxgene.org/), using the calcPhenotype function to estimate drug responses [half maximal inhibitory concentration (IC₅₀) values], and setting the significance cutoff to P value <0.005 to screen for relevant drugs.

Statistical analysis

All statistical assessments were performed using R software version 4.2.2. Results were considered statistically significant when P values <0.05. The analysis was based on publicly available data sets, and non-parametric tests were used to compare groups when appropriate. The significance levels are expressed as follows: *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001.

Results

Screening of PCD-related lncRNAs in ESCA

Transcriptome datasets from the TCGA were analyzed to accurately lncRNA associated with PCD in ESCA. Initially, PCD-related genes were derived from a comprehensive review by Zou et al. (21). To identify lncRNAs significantly correlated with these genes, Pearson correlation analysis was performed. As shown in Figure 2A, a Sankey diagram was created to visually illustrate the associations between PCD-related genes and their corresponding lncRNAs in ESCA. Subsequently, in order to identify lncRNAs with prognostic value from lncRNAs associated with PCD, we performed a univariate regression analysis. The Venn diagram (Figure 2B) shows the overlap between PCD-associated lncRNAs and lncRNAs with univariate Cox P values below 0.05. This process identified 146 lncRNAs that were both associated with PCD and significantly associated with OS in ESCA patients (Table S1).

Figure 2 Identification and prognostic analysis of PCD-related lncRNAs in ESCA. (A) The Sankey diagram shows the relationship between PCD-related genes and lncRNAs. (B) The Venn diagram shows the overlap between PCD-related lncRNAs and prognostic lncRNAs (univariate Cox regression P<0.05). ESCA, esophageal carcinoma; lncRNAs, long non-coding RNAs; PCD, programmed cell death.

Construction and validation of the PCD-related lncRNAs signature

From TCGA, we initially selected 160 cases of ESCA. The cases were randomly divided into two groups: a training set (n=80) and a validation set (n=80). Subsequently, both univariate Cox regression and LASSO Cox regression were applied to investigate 146 PCD-related lncRNAs. Through cross-validation and examination of regression coefficient trajectories, 14 candidate lncRNAs were selected for subsequent modeling (Figure 3A,3B). After optimization by multivariate Cox regression analysis, 8 lncRNAs were finally selected AC109347.1, BLACE, AP001527.2, AP001001.1, LINC00402, AC087289.5, FAM83C-AS1, and AL132655.2 (Figure 3C,3D). A risk assessment signature was constructed using these eight lncRNAs, and the risk score for each patient was calculated using a formula, risk score = ∑ (Expression × corresponding regression coefficient), where the expression level reflects the number of each lncRNA, and the coefficient is from a multivariate Cox analysis. To distinguish patient groups, the median risk score was used to classify patients into high- and low-risk groups. The distribution pattern of risk scores, patient survival status, and classification between the two groups is shown in (Figure 3E), which combines risk score visualization, survival results, and a heatmap of gene expression. Kaplan-Meier analysis (Figure 3F) showed that patients classified as high-risk had significantly worse OS than patients in the low-risk group (P<0.001), emphasizing the prognostic utility of this signature. The stability of this signature was further verified by time-dependent ROC curve analysis, with AUC values reaching 0.902 at 1 year, 0.871 at 2 years, and 0.860 at 3 years for OS prediction (Figure 3G), indicating strong predictive accuracy. The chromosomal locations of eight prognostic lncRNAs are annotated in Figure 3H. To study their functional significance, correlation heatmaps were created to examine the relationship between these lncRNAs and genes known to be involved in PCD. As shown in Figure 3I, several lncRNAs show strong correlations with key PCD regulators, implying their possible role in regulating the PCD mechanism in ESCA.

Figure 3 Development of a prognostic signature based on PCD-related lncRNAs. (A,B) LASSO-Cox regression analysis was performed to develop the prognostic signature. (C,D) Forest plot displays HR with 95% CI for prognostic variables, while horizontal bar chart visualizes LASSO-selected gene coefficients, distinguishing protective (green) and risk-associated (red) markers. (E) Ranked dot, heatmap and scatter plots of the signature gene expressions in Train coherent. (F) Kaplan-Meier survival curve of OS of patients in high- and low-risk groups. (G) ROC curves and AUC for 1-, 2-, and 3-year survival rates. (H) A circular plot illustrating gene location distribution across chromosomes. (I) Heatmap of the correlation between PCD-related genes and 8 prognostic PCD-related lncRNAs. Statistical significance: *, P<0.05; **, P<0.01; ***, P<0.001. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; lncRNAs, long non-coding RNAs; OS, overall survival; PCD, programmed cell death; ROC, receiver operating characteristic.

To test the reliability and predictability of this signature, survival analysis was performed using two separate internal datasets, the TCGA training dataset and the TCGA entire dataset. The risk score distribution and survival results in the TCGA-test set are presented in (Figure 4A), and the Kaplan-Meier survival curve is presented in (Figure 4B), confirming that the survival of low-risk patients is significantly better than that of high-risk individuals (P<0.0001). Similarly, Figure 4C,4D presents the analysis results of the entire TCGA dataset, showing highly significant differences in survival between the low- and high-risk groups (P<0.0001), consistent with previous results. The persistent differences in lncRNA expression patterns observed among these groups highlight the reliability of this signature in classifying ESCA patients based on risk levels, as well as the ability to derive stable prognostic information from risk assessments and associated lncRNA expression profiles.

Figure 4 The evaluation and validation of prognostic signatures in the test group and the entire group. (A,C) The risk score distribution, survival status, and heatmap of prognostic PCD-related lncRNAs expression in the corresponding set. (B,D) Kaplan-Meier survival analysis curves and ROC curves, with curves based on patients’ risk scores. Green, blue, and red indicate 1-, 2-, and 3-year AUC, respectively. AUC, area under the curve; lncRNAs, long non-coding RNAs; PCD, programmed cell death; ROC, receiver operating characteristic.

Association between PCD-related lncRNAs signature and clinicopathological characteristics

To assess the clinical utility of PCD-related lncRNAs signatures, we conducted extensive analysis of their correlation with various clinical characteristics. Specifically, as shown in Figure 5A, we investigated the relationship between the risk score distribution and different clinical pathology factors in patients with ESCA. The box plots in Figure 5B-5H show that patients with higher T stages (Figure 5E), patients with lymph node metastases (Figure 5F), patients with distant metastases (Figure 5G), and patients who died (Figure 5D) tended to have higher risk scores. In contrast, there were no significant differences in risk scores when grouped by age (Figure 5B), gender (Figure 5C), or overall clinical stage (Figure 5H). These observations show that this prognostic characteristic is closely related to tumor severity and progression indicators, highlighting its clinical significance for ESCA. To further confirm its predictive accuracy, we performed stratified survival analyses in various clinical subgroups. Figure 6 show the results of the Kaplan-Meier survival analysis, consistently demonstrating a continuing trend across all subgroups: patients classified as high-risk had significantly shorter OS times than those in the low-risk group. The prognostic ability of this signature has proven to be stable and reliable under different clinical parameters, supporting its potential application value in personalized risk stratification of ESCA.

Figure 5 Association between the PCD-related lncRNAs signature and clinicopathological characteristics. (A) Heatmap revealing the correlation between clinical features and risk scores in ESCA patients, including age, gender and clinical stage. (B-H) Evaluating variations in risk scores across different clinical features. Statistical significance: *, P<0.05; **, P<0.01; ***, P<0.001. ESCA, esophageal carcinoma; lncRNAs, long non-coding RNAs; M, metastasis; N, node; ns, not significant; PCD, programmed cell death; T, tumor.

Figure 6 Stratified survival analysis of the PCD-related lncRNAs signature across different clinical subgroups. Patients were categorized into subgroups based on age (A,B), gender (C,D), clinical T stage (E,F), clinical N stage (G,H), clinical M stage (I,J), and clinical stage (K,L). Survival probabilities are displayed for each group (high-risk vs. low-risk), with P value highlighting significant differences between the groups. LncRNAs, long non-coding RNAs; M, metastasis; N, node; PCD, programmed cell death; T, tumor.

Genomic alterations and TMB analysis

Somatic mutation analysis revealed that the mutation frequencies were nearly identical in the high-risk (98.73%) and low-risk (98.75%) groups (Figure 7). However, despite this similarity in overall frequency, distinct differences in mutation patterns were observed between the two groups. Mutations were most common in individuals with high TMB in TP53 (82%), TTN (42%), SYNE1 (20%), MUC16 (20%), LRP1B (19%), and ARID1A (16%) (Figure 7A). Notably, the high-risk group showed significantly higher mutation frequencies in genes such as ARID1A, MRC1, and SCN9A (P<0.05), whereas FLG mutations were more prevalent in the low-risk group, as illustrated in Figure 7B. In contrast, the low TMB group mainly showed mutations in TP53 (89%), TTN (35%), FLG (21%), CSD3 (20%), MUC16 (19%) and SYNE1 (19%) (Figure 7C). In addition to somatic mutation, we also studied the association of TMB with clinical prognosis in ESCA. Kaplan-Meier analysis (Figure 7D) revealed that higher TMB levels have significantly shortened OS compared with patients with lower TMB levels. In addition, as shown in Figure 7E, individuals with high TMB and high-risk scores have the worst survival outcomes, while patients with low TMB and low risk scores have the most favorable prognosis.

Figure 7 Genomic alterations analysis in high- and low-risk ESCA patients. (A,C) Waterfall plots showing the somatic mutation landscape in the high- and low-risk groups. (B) Forest plot displaying differentially mutated genes between high- and low-risk groups. (D) Kaplan-Meier survival curves comparing ESCA patients’ overall survival stratified by risk score. (E) Kaplan-Meier survival curves showing the combined effect of TMB and risk score on overall survival. CI, confidence interval; ESCA, esophageal carcinoma; OR, odds ratio; TMB, tumor mutation burden.

Nomogram construction

We constructed a predictive nomogram (Figure 8A) to determine whether features of our risk score could function as independent predictors of patient prognosis. This tool integrates multiple clinical factors, such as age, clinical T stage, clinical N stage, clinical M stage, and risk score, into a single numerical model, providing a direct and quantitative risk assessment method. A nomogram is an easy-to-use statistical tool that can be used to estimate survival prospects and simplify the clinical application of risk scores. Our results show that this nomogram has reliable predictive accuracy for patients’ OS at 1, 2, and 3 years (Figure 8B). Its consistency index is significantly higher than that of a single predictor, highlighting the enhanced predictive power of combined features (Figure 8C). Based on Kaplan-Meier analysis, patients with higher nomogram scores exhibited significantly worse OS (Figure 8D). The predictive performance measured by AUC was 0.691, 0.713, and 0.794 for 1-, 2-, and 3-year survival, respectively (Figure 8E), reinforcing the predictive stability and clinical importance of this prognostic nomogram.

Figure 8 Construction and validation of the nomogram. (A) Nomogram incorporating significant prognostic factors for predicting 1-, 2-, and 3-year OS. (B) Calibration curves assessing the predictive accuracy of the nomogram for 1-, 2-, and 3-year OS. (C) Time-dependent C-index comparison between the nomogram and individual prognostic factors. (D) Kaplan-Meier survival analysis based on nomogram-derived risk stratification. (E) ROC curve analysis for 1-, 2-, and 3-year OS predictions using the nomogram. AUC, area under the curve; C-index, concordance index; M, metastasis; N, node; OS, overall survival; ROC, receiver operating characteristic; T, tumor.

WGCNA and GSEA

To explore the functional significance of lncRNA associated with PCD, we used their expression data to conduct a WGCNA network. To construct a co-expression network, a soft threshold power of 12 was applied, and 14 independent gene modules were identified through hierarchical clustering, each module having a unique color distinction (Figure 9A). Among the modules, the red one exhibited the strongest positive correlation with the risk score, with a correlation coefficient of 0.49 and a P value of 5.5×10⁻⁷², while the turquoise module displayed the greatest inverse correlation, with a correlation coefficient of −0.38 and an exceptionally low P value of 1.2×10⁻¹¹³ (Figure 9B). Further inspection showed that these modules were most closely associated with PCD-related phenotypes (Figure 9C). GO enrichment analysis revealed that the genes in these modules are predominantly engaged in immune-related biological processes, such as immune response regulation, leukocyte immunity, and cytokine receptor activity (Figure 9D). KEGG pathway analysis shows that these genes are significantly concentrated on immune-related pathways, including cytokine-cytokine receptor interactions, antigen processing and presentation, and autoimmune disease-related pathways (Figure 9E). Additionally, hallmark pathway analysis indicated that these genes were primarily involved in immunological and inflammatory processes, including interferon gamma response, allograft rejection, and the IL-2/STAT5 signaling pathway (Figure 9F).

Figure 9 WGCNA analysis of PCD-related lncRNAs signature. (A) Cluster dendrogram of gene co-expression modules detected by WGCNA, with colors indicating module assignment. (B) Module-trait relationship heatmap showing correlations between module eigengenes and risk scores (correlation coefficients and P value are shown). (C) Correlation between module membership and gene significance in the red and turquoise modules. (D) GO enrichment analysis (BP, CC, MF) of genes associated with PCD-lncRNAs. (E) KEGG pathway enrichment analysis of these mRNAs. (F) Hallmark pathway enrichment analysis revealing key functional signatures in esophageal carcinoma. BP, biological process; CC, cellular component; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNAs, long non-coding RNAs; lncRNAs, long non-coding RNAs; MF, molecular function; MHC, major histocompatibility complex; PCD, programmed cell death; PCD, programmed cell death; WGCNA, weighted gene co-expression network analysis.

We performed a GSEA to gain insight into the biological mechanisms that may explain the observed differences in survival rates between high- and low-risk groups, as defined by the prognostic signature based on PCD-related lncRNAs. Enrichment of hallmark pathways showed that PCD-related lncRNA signatures were associated with pathways in the dataset such as epithelial-mesenchymal transition (EMT), the interleukin 6 (IL-6)/Janus kinase (JAK)/signal transducer and activator of transcription 3 (STAT3) signaling pathway, the tumor necrosis factor α (TNF-α) signaling pathway through nuclear factor kappa B (NF-κB), and allograft rejection (Figure 10). The high-risk group exhibited greater activation of immune-associated pathways, including EMT, IL-6/JAK/STAT3 signaling pathway, TNF-α/NF-κB signaling, and allograft rejection (Figure 10). Notably, these pathways showed substantial overlap with those enriched in the red and turquoise gene modules, both of which were previously linked to immune-related functions. In contrast, several of these immune pathways appeared downregulated in the low-risk group (Figure 10C).

Figure 10 GSEA of hallmark pathways associated with PCD-lncRNAs signature. (A) Bubble plot showing hallmark pathways significantly enriched in the high-risk and low-risk groups based on GSEA. (B) GSEA enrichment plot of activated pathways. (C) GSEA enrichment plot of suppressed pathways. GSEA, gene set enrichment analysis; lncRNAs, long non-coding RNAs; PCD, programmed cell death.

Immune infiltration analysis

At the outset, to explore the immune landscape of ESCA and the characteristics of its tumor microenvironment, we evaluated immune, stromal, and ESTIMATE scores in patients stratified by risk level. The results showed that compared with the low-risk group, individuals in the high-risk group had significantly lower scores on these three indicators and were statistically significant (Figure 11A). We performed ssGSEA to investigate the association between risk scores and tumor-infiltrating immune cells, and found that immune infiltration patterns differed notably between the two risk groups. Specifically, high-risk patients exhibited increased levels of activated dendritic cells, type 17 helper T cells, as well as various CD8⁺ and memory T cell populations. By contrast, patients categorized as low risk exhibited increased levels of neutrophils, immature dendritic cells, and plasmacytoid dendritic cells (Figure 11B). Moreover, significant differences were observed in the infiltration levels of 19 distinct immune cell types between the high- and low-risk cohorts (Figure 11C).

Figure 11 Immune infiltration analysis in high- and low-risk ESCA patients. (A) Box plots comparing ESTIMATE score, immune score, and stromal score between high- and low-risk groups. (B) Heatmap displaying the infiltration levels of tumor-infiltrating immune cells in high- and low-risk groups based on ssGSEA. (C) Violin plots showing the distribution of immune cell infiltration levels between high- and low-risk groups. Statistical significance: *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. ESCA, esophageal carcinoma; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; MDSC, myeloid-derived suppressor cells; ns, not significant; ssGSEA, single-sample gene set enrichment analysis.

PCD-related lncRNAs correlate with therapeutic drugs

To further explore the possible therapeutic significance of our prognosis signature, we used the we used the OncoPredict R package to conduct drug sensitivity assessments combined with data from the GDSC2 database. Our analysis (Figure 12) revealed that the predicted IC₅₀ values of SB505124, TAF1_5496, and Ulixertinib were significantly lower in the high-risk ESCA group compared to the low-risk group (Wilcoxon’s test, P<0.05). These findings imply that these drugs could potentially exhibit enhanced therapeutic effects in patients with high-risk classification, although additional experimental studies are needed to confirm this. Overall, these results highlight the potential role of PCD-related lncRNAs signatures in guiding personalized treatment options by identifying promising drug candidates for ESCA treatment.

Figure 12 Potential drugs for patients in the high- and low-risk groups. (A-C) Spearman correlations between risk score and predicted sensitivity scores for SB505124, TAF1_5496, and Ulixertinib. (D-F) Boxplots comparing drug sensitivity scores between high- and low-risk groups. Statistical significance: *, P<0.05; **, P<0.01; ***, P<0.001. ns, not significant.

Discussion

ESCA is a highly aggressive malignancy with poor prognosis and limited treatment options. Given its heterogeneous clinical presentation, identifying reliable prognostic markers is essential for effective patient stratification and personalized therapy. PCD, such as apoptosis, ferroptosis, pyroptosis, and disulfidptosis, plays a pivotal role in tumor progression, immune regulation, and treatment resistance (10). LncRNAs have emerged as key regulators of gene expression through chromatin remodeling, RNA stabilization, and protein interactions (30). Increasing evidence suggests that lncRNA are involved in non-apoptotic PCD pathways such as ferroptosis, pyroptosis, and necroptosis, influencing tumor behavior and therapeutic responses (18). Prognostic signatures based on PCD-related lncRNAs have been proposed in cancers such as hepatocellular carcinoma, colorectal cancer, and lung adenocarcinoma (31-33). However, their role in ESCA remains largely uncharacterized. To fill this gap, we developed and validated a novel signature consisting of PCD-associated lncRNAs aimed at predicting the prognosis of ESCA and analyzing its immune microenvironment.

Our identification of a PCD-related lncRNA signature enabled the effective stratification of ESCA patients into high- and low-risk groups, with significantly poorer survival observed in the high-risk cohort. Bridging the gap between bioinformatics findings and clinical utility is a critical goal. In this context, the PCD-related lncRNA signature demonstrates potential value as a prognostic candidate biomarker for patient risk stratification and personalized management. The distinct high- and low-risk classifications provide a testable hypothesis for future clinical investigation, particularly regarding whether these groups respond differentially to various treatment modalities, such as immunotherapy or chemotherapy. Although rigorous experimental validation is essential to establish the clinical relevance of these lncRNAs, our findings lay a valuable bioinformatics foundation for identifying high-risk subsets and guiding future mechanism research in ESCA personalized management. Compared with traditional prognostic factors such as tumor-node-metastasis (TNM) staging or pathological grading, the PCD-related lncRNA signature captures tumor heterogeneity at the molecular level. This transcriptome-based approach may reveal high-risk subgroups not identifiable by standard clinicopathological features, thereby providing additional prognostic value. Correlation analyses further showed that patients with lymph node metastases, distant metastases, and poor survival outcomes had significantly higher risk scores. These associations indicate that the PCD-related lncRNA signature reflects clinically relevant tumor progression features, and may complement established prognostic parameters. This aligns with previous findings that underscore the role of lncRNAs in immune regulation and tumor progression (34). Functional enrichment analysis of genes associated with the PCD-related lncRNAs, identified via WGCNA, revealed significant involvement in immune-related biological processes and signaling pathways. GO and KEGG analyses highlight their role in cytokine receptor function, antigen processing and presentation, and autoimmune responses. GSEA confirmed the enrichment of IL-6/JAK/STAT3 and TNF-α/NF-κB pathways, both of which are critical in mediating inflammatory responses and promoting immune escape (35).

The IL-6/JAK/STAT3 signaling cascade is recognized because it promotes tumor growth, immune escape, and metastasis spread, mainly through STAT3 that increases the transcription of genes related to cell survival and proliferation (36). A variety of lncRNAs affect this pathway, usually by acting as competing endogenous RNAs (ceRNAs) to isolate specific miRNAs or regulate the expression of upstream cytokines. The lncRNA maternal expressed gene 3 (MEG3) regulates STAT3 activity by sponging miR-708, thereby inhibiting the growth and stemness of colorectal cancer cells (37). Similarly, TNF-α/NF-κB signaling is the fundamental axis of chronic inflammation associated with cancer. LncRNA such as LINC02820 have been shown to augment this pathway, remodel the cytoskeleton, and potentiate metastasis in ESCC (38). These findings underscore the need for further exploration into whether the PCD-related lncRNAs identified in our signature may similarly modulate IL-6/JAK/STAT3 or TNF-α/NF-κB signaling, potentially influencing inflammatory responses and tumor progression in ESCA.

Tumor microenvironment analysis revealed that low-risk patients generally exhibited higher immune, stromal, and ESTIMATE scores, suggesting that their tumor microenvironment favors immune involvement. In comparison, high-risk patients demonstrated greater infiltration of activated dendritic cells and various T cell subsets, which may reflect altered immune activation states or dysfunction. Previous studies have shown that certain long-chain non-coding RNAs can affect T cell depletion and dendritic cell maturation, which can hinder effective antigen presentation and immune checkpoint function (39).

From a therapeutic perspective, drug sensitivity analyses highlighted three promising compounds, SB505124 (a TGF-β pathway inhibitor), TAF1_5496 (a TAF1 inhibitor) and Ulixertinib (an ERK pathway inhibitor), that may be more effective in high-risk patients who have lower IC₅₀ values for these drugs. These findings suggest that lncRNAs associated with PCD may regulate drug responses by affecting specific pathways or directly changing drug targets, which is consistent with previous evidence from other gastrointestinal malignancies (40). Beyond conventional chemotherapy, immunotherapy has emerged as an important treatment strategy across various malignancies. Recent evidence indicated that PD-1, PD-L1, NY-ESO-1, and MAGE-A4 expression are interrelated and may contribute to tumor aggressiveness (41). In the present study, the PCD-related lncRNA signature was closely associated with immune-related pathways and distinct immune cell infiltration patterns, suggesting that it may help stratify ESCA patients for immunotherapy. Future research should further explore whether patients in different risk groups respond differently to immune checkpoint inhibitors or antigen-targeted strategies, thereby linking molecular risk stratification with potential clinical decision-making.

Crucially, our prognosis signature simultaneously integrates multiple PCD processes. Although early research focused on single cell death pathways, emerging research emphasized the complex interactions between apoptosis, autophagy, iron death and pyroapoptosis. LINC00402 has previously been recognized as a ferroptosis-related lncRNA and linked to inflammatory responses within the tumor microenvironment (42), while the roles of other PCD-related lncRNAs such as AC109347.1, BLACE, and AP001527.2 remain to be experimentally validated in the context of apoptosis and autophagy. Our multi-PCD lncRNAs-based framework offers an expanded perspective on tumor susceptibility and immune participation, reflecting the potential interplay between different cell death mechanisms-an improvement over single-pathway model.

Nonetheless, certain limitations should be acknowledged. Its retrospective design and reliance solely on the TCGA database limit the breadth of the study results. In order to determine the causal significance of these lncRNAs, an independent prospective set is needed to conduct additional verification and conduct laboratory-based in vitro and in vivo studies. In addition, spatial transcriptomics and single cell sequencing techniques can provide more detailed insights into the role of these lncRNAs in the tumor microenvironment. In short, the lncRNAs characteristics associated with PCD provide a valuable basis for risk assessment, immunotherapy prediction and mechanism research of ESCA, highlighting the need for further verification and improvement.

Conclusions

In conclusion, this study systematically identified a PCD-related lncRNA signature (AC109347.1, BLACE, AP001527.2, AP001001.1, LINC00402, AC087289.5, FAM83C-AS1, and AL132655.2) in ESCA based on transcriptomic and clinical data from the TCGA database. The signature showed promising prognostic relevance by distinguishing patients with different risk levels, suggesting its possible value in supporting personalized prognosis assessment. In addition, functional enrichment analyses suggested its association with immune suppression and inflammation-related pathways, while drug sensitivity analysis pointed to several compounds that may merit further investigation as targeted therapeutic candidates. These findings hold promise in guiding clinical management and facilitating individualized therapy approaches for ESCA patients in the future.

Acknowledgments

We acknowledge the TCGA database for providing the platform and its contributors for uploading meaningful datasets.

Footnote

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

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

Funding: This work was supported by grants from the China Health and Medical Development Foundation (No. chmdf2024-xrzx01-10), the National Natural Science Foundation of China (No. 82203096), and the Chinese Postdoctoral Science Foundation (No. 2022M711717).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1598/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. Bray F, Ferlay J, Soerjomataram I, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68:394-424. [Crossref] [PubMed]
2. Huang J, Koulaouzidis A, Marlicz W, et al. Global Burden, Risk Factors, and Trends of Esophageal Cancer: An Analysis of Cancer Registries from 48 Countries. Cancers (Basel) 2021;13:141. [Crossref] [PubMed]
3. He S, Xu J, Liu X, et al. Advances and challenges in the treatment of esophageal cancer. Acta Pharm Sin B 2021;11:3379-92. [Crossref] [PubMed]
4. Rustgi AK, El-Serag HB. Esophageal carcinoma. N Engl J Med 2014;371:2499-509. [Crossref] [PubMed]
5. Thrift AP. Esophageal Adenocarcinoma: The Influence of Medications Used to Treat Comorbidities on Cancer Prognosis. Clin Gastroenterol Hepatol 2015;13:2225-32. [Crossref] [PubMed]
6. Zhou X, Ren T, Zan H, et al. Novel Immune Checkpoints in Esophageal Cancer: From Biomarkers to Therapeutic Targets. Front Immunol 2022;13:864202. [Crossref] [PubMed]
7. Sundar R, Smyth EC, Peng S, et al. Predictive Biomarkers of Immune Checkpoint Inhibition in Gastroesophageal Cancers. Front Oncol 2020;10:763. [Crossref] [PubMed]
8. Noori M, Yousefi AM, Zali MR, et al. Predictive value of PD-L1 expression in response to immune checkpoint inhibitors for esophageal cancer treatment: A systematic review and meta-analysis. Front Oncol 2022;12:1021859. [Crossref] [PubMed]
9. Sazonova EV, Yapryntseva MA, Pervushin NV, et al. Cancer Drug Resistance: Targeting Proliferation or Programmed Cell Death. Cells 2024;13:388. [Crossref] [PubMed]
10. Tang D, Kang R, Berghe TV, et al. The molecular machinery of regulated cell death. Cell Res 2019;29:347-64. [Crossref] [PubMed]
11. Carneiro BA, El-Deiry WS. Targeting apoptosis in cancer therapy. Nat Rev Clin Oncol 2020;17:395-417. [Crossref] [PubMed]
12. Plati J, Bucur O, Khosravi-Far R. Dysregulation of apoptotic signaling in cancer: molecular mechanisms and therapeutic opportunities. J Cell Biochem 2008;104:1124-49. [Crossref] [PubMed]
13. Li J, Cao F, Yin HL, et al. Ferroptosis: past, present and future. Cell Death Dis 2020;11:88. [Crossref] [PubMed]
14. Dikic I, Elazar Z. Mechanism and medical implications of mammalian autophagy. Nat Rev Mol Cell Biol 2018;19:349-64. [Crossref] [PubMed]
15. Xia X, Wang X, Cheng Z, et al. The role of pyroptosis in cancer: pro-cancer or pro-”host”? Cell Death Dis 2019;10:650. [Crossref] [PubMed]
16. Yu P, Zhang X, Liu N, et al. Pyroptosis: mechanisms and diseases. Signal Transduct Target Ther 2021;6:128. [Crossref] [PubMed]
17. Statello L, Guo CJ, Chen LL, et al. Gene regulation by long non-coding RNAs and its biological functions. Nat Rev Mol Cell Biol 2021;22:96-118. [Crossref] [PubMed]
18. Huang K, Yu L, Lu D, et al. Long non-coding RNAs in ferroptosis, pyroptosis and necroptosis: from functions to clinical implications in cancer therapy. Front Oncol 2024;14:1437698. [Crossref] [PubMed]
19. Wang M, Mao C, Ouyang L, et al. Long noncoding RNA LINC00336 inhibits ferroptosis in lung cancer by functioning as a competing endogenous RNA. Cell Death Differ 2019;26:2329-43. [Crossref] [PubMed]
20. Huang ZL, Chen RP, Zhou XT, et al. Long non-coding RNA MEG3 induces cell apoptosis in esophageal cancer through endoplasmic reticulum stress. Oncol Rep 2017;37:3093-9. [Crossref] [PubMed]
21. 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]
22. Robin X, Turck N, Hainard A, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics 2011;12:77. [Crossref] [PubMed]
23. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
24. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med 2013;32:5381-97. [Crossref] [PubMed]
25. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
26. Xu S, Hu E, Cai Y, et al. Using clusterProfiler to characterize multiomics data. Nat Protoc 2024;19:3292-320. [Crossref] [PubMed]
27. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
28. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
29. 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]
30. Marchese FP, Raimondi I, Huarte M. The multidimensional mechanisms of long noncoding RNA function. Genome Biol 2017;18:206. [Crossref] [PubMed]
31. Song Z, Cao X, Wang X, et al. A disulfidptosis-related lncRNA signature for predicting prognosis and evaluating the tumor immune microenvironment of lung adenocarcinoma. Sci Rep 2024;14:4621. [Crossref] [PubMed]
32. Liu Y, Meng J, Ruan X, et al. A disulfidptosis-related lncRNAs signature in hepatocellular carcinoma: prognostic prediction, tumor immune microenvironment and drug susceptibility. Sci Rep 2024;14:746. [Crossref] [PubMed]
33. Xiao L, Yin W, Chen X, et al. A disulfidptosis-related lncRNA index predicting prognosis and the tumor microenvironment in colorectal cancer. Sci Rep 2023;13:20135. [Crossref] [PubMed]
34. Park EG, Pyo SJ, Cui Y, et al. Tumor immune microenvironment lncRNAs. Brief Bioinform 2022;23:bbab504. [Crossref] [PubMed]
35. Johnson DE, O’Keefe RA, Grandis JR. Targeting the IL-6/JAK/STAT3 signalling axis in cancer. Nat Rev Clin Oncol 2018;15:234-48. [Crossref] [PubMed]
36. Wang Y, Wang Z, Li S, et al. Deciphering JAK/STAT signaling pathway: A multifaceted approach to tumorigenesis, progression and therapeutic interventions. Int Immunopharmacol 2024;131:111846. [Crossref] [PubMed]
37. Rahbar Farzam O, Najafi S, Amini M, et al. Interplay of miRNAs and lncRNAs in STAT3 signaling pathway in colorectal cancer progression. Cancer Cell Int 2024;24:16. [Crossref] [PubMed]
38. Wang J, Huang TJ, Mei Y, et al. Novel long noncoding RNA LINC02820 augments TNF signaling pathway to remodel cytoskeleton and potentiate metastasis in esophageal squamous cell carcinoma. Cancer Gene Ther 2023;30:375-87. [Crossref] [PubMed]
39. Li K, Wang Z. Non-coding RNAs: Key players in T cell exhaustion. Front Immunol 2022;13:959729. [Crossref] [PubMed]
40. Song X, Hou L, Zhao Y, et al. Metal-dependent programmed cell death-related lncRNA prognostic signatures and natural drug sensitivity prediction for gastric cancer. Front Pharmacol 2022;13:1039499. [Crossref] [PubMed]
41. Hashimoto K, Nishimura S, Ito T, et al. Immunohistochemical expression and clinicopathological assessment of PD-1, PD-L1, NY-ESO-1, and MAGE-A4 expression in highly aggressive soft tissue sarcomas. Eur J Histochem 2022;66:3393. [Crossref] [PubMed]
42. Guo S, Chen J, Yi X, et al. Identification and validation of ferroptosis-related lncRNA signature as a prognostic model for skin cutaneous melanoma. Front Immunol 2022;13:985051. [Crossref] [PubMed]