Prognostic significance of DNA damage response-related markers in esophageal squamous cell carcinoma using machine learning approaches

Introduction

Esophageal squamous cell carcinoma (ESCC) represents one of the most prevalent histological subtypes of esophageal cancer worldwide and is associated with high incidence and mortality rates, particularly in regions including China, Southeast Asia, and East Africa (1). Despite advances in diagnostic and therapeutic strategies, the 5-year survival rate for patients with ESCC remains low, as most cases are diagnosed at an advanced stage with limited availability of effective targeted treatment options (2). Consequently, further investigation into the molecular mechanisms underlying ESCC, as well as the identification of potential prognostic biomarkers and therapeutic targets, remains a key priority in this field.

Advances in genomics have established a key foundation for precise molecular classification of tumors and the development of personalized therapeutic strategies. The Cancer Genome Atlas (TCGA), a large-scale, multidimensional cancer genomics project, comprises whole-genome sequencing data alongside clinical data across multiple cancer types and serves as a valuable resource for investigating molecular mechanisms and identifying potential biomarkers in oncology research (3). Utilization of gene expression profiles and clinical data from the TCGA-ESCC dataset enables bioinformatics analyses that characterize molecular features associated with ESCC progression and prognosis, thereby providing a basis for the development of individualized treatment approaches.

Homologous recombination deficiency (HRD) has emerged as an important indicator of genomic instability in cancer research and has been closely associated with malignant progression, treatment resistance, and adverse prognosis across multiple cancer types (4). The HRD score reflects the degree of genomic instability within tumors and serves as a key biomarker for predicting patient responses to therapies targeting DNA damage response (DDR) and repair (5). However, the use of HRD scoring in ESCC remains insufficiently investigated, and evaluation of its association with prognosis in patients with ESCC may provide further insight into precision therapeutic strategies.

Immune cell infiltration within the tumor microenvironment is recognized as a key factor influencing ESCC progression and prognosis (6). Complex interactions between immune cells and tumor cells contribute to immune evasion and may affect patient responses to immunotherapy (7). Analysis of immune cell infiltration in ESCC samples may therefore clarify its association with prognosis and support the identification of potential immunotherapeutic targets.

This study used ESCC samples obtained from the TCGA database and integrated HRD scores with DDR-related genes to construct a molecular model associated with prognosis in patients with ESCC. Differentially expressed genes (DEGs) analysis, gene set enrichment analysis (GSEA), mutation profiling, and immune infiltration assessment were performed to assess the impact of HRD scores on ESCC prognosis. Multiple machine learning approaches were used to identify key genes and to develop prognostic models. Given that the incidence and mortality of ESCC vary across geographic regions and socioeconomic contexts, the importance of assessing potential health disparities was acknowledged. As the TCGA database does not include detailed sociodemographic variables, future investigations incorporating more diverse populations will be necessary to assess the generalizability of these findings across different sociodemographic groups and to support equitable clinical application of the HRD-based prognostic model. These findings are intended to provide further insight and a theoretical basis for molecular classification, prognostic evaluation, and personalized treatment strategies in ESCC. We present this article in accordance with the TRIPOD+AI reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0588/rc).

Methods

Study design

The technology roadmap is presented in Figure 1.

Figure 1 Technology roadmap. DEG, differentially expressed gene; ESCC, esophageal squamous cell carcinoma; GBM, gradient boosting machine; GSEA, gene set enrichment analysis; HRD, homologous recombination deficiency; KM, Kaplan-Meier; LASSO, least absolute shrinkage and selection operator; ORA, over‑representation analysis; PCA, principal component analysis; QC, quality control; RSF, random survival forest; Survival-SVM, survival support vector machine; TCGA, The Cancer Genome Atlas.

Data download

The R package TCGAbiolinks was used to access esophageal cancer [esophageal carcinoma (ESCA)] data from the TCGA database (https://portal.gdc.cancer.gov/) (TCGA-ESCA) (8). Subsequently, ESCC-specific data were extracted and recorded as TCGA-ESCC. After exclusion of samples lacking clinical data or with an overall survival time of less than 1 day, sequencing data in fragments per kilobase of transcript per million mapped reads (FPKM) format and corresponding clinical data from 78 patients diagnosed with ESCC were retained. A random sampling approach was then used, whereby 30% of the samples were assigned to the test cohort and 70% to the training cohort for construction of a prognostic model based on DDR-related genes.

The single-cell dataset GSE199619, comprising 10 ESCC samples and 7 normal samples, was used (9). Sequencing was conducted using the GPL24676 and GPL20838 platforms, with all samples derived from Homo sapiens. DDR-related genes and HRD scores reported in published literature were obtained from the PubMed database (10).

Microsatellite instability, mutation count, and tumor mutation burden data were retrieved from the cBioPortal database (11). The fraction of genome altered (FGA), defined as the proportion of genomic regions exhibiting chromosomal copy number alterations (CNAs) outside the region, was obtained. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

HRD score acquisition and processing

HRD scores for TCGA samples were obtained from a previously published pan-cancer study conducted by Rempel et al., which provides comprehensive genomic scar scores across 33 cancer types (10). Affymetrix single nucleotide polymorphisms (SNPs) 6.0 array data from TCGA were processed using the PennCNV algorithm to identify CNA. HRD scores were calculated as the unweighted sum of three genomic instability metrics: (I) loss of heterozygosity (LOH), defined as regions of LOH greater than 15 Mb but smaller than whole chromosomes (12); (II) telomeric allelic imbalance, defined as regions with allelic imbalance extending to telomeric regions (13); and (III) large-scale state transitions, defined as chromosomal breaks between adjacent regions of at least 10 Mb (14). Detailed calculation protocols are described in the original publication.

For this study, HRD scores specific to ESCC samples were extracted by matching TCGA patient identifiers. Following the exclusion of samples with missing clinical data or survival time of less than 1 day, HRD scores for 78 patients with ESCC were retained for analysis. The distribution of HRD scores in this cohort ranged from 0 to 89 (median =34, mean =36.2).

HRD score grouping

To achieve biologically meaningful stratification, maximally selected rank statistics (MSRS) were applied using the maxstat R package to determine the optimal HRD score cutoff based on overall survival outcomes. This approach evaluates all possible cut points and identifies the value that maximizes the separation of survival curves, as determined by the standardized log-rank statistic. The optimal cutoff was identified as 38 (Figure 2C). Patients were subsequently categorized into high-HRD (score ≥38, n=32) and low-HRD (score <38, n=46) groups for all comparative analyses, including differential expression, pathway enrichment, and mutation profiling.

Figure 2 HRD score. (A) Histogram illustrating the distribution of HRD scores. The blue and red regions represent the optimal grouping of HRD scores. (B) Comparative plots of partial genomic changes between high- and low-HRD score groups, with asterisks indicating levels of statistical significance (***, P<0.001). (C) MSRS plot demonstrating the optimal cutoff value for HRD score stratification. (D) Volcano plot displaying DEGs, with statistical significance and logFC. Each dot represents an individual gene; red indicates significant upregulation, blue indicates significant downregulation, and gray indicates no statistically significant difference. (E) Heatmap illustrating gene expression patterns across samples. Color intensity represents expression levels, with red indicating higher expression and blue indicating lower expression. (F) Kaplan-Meier survival curves comparing overall survival between high- and low-HRD score groups, with the P value indicating statistical significance. (G) Time-dependent ROC curves demonstrating the relationship between sensitivity and 1 − specificity for model predictions, along with AUC values at different time points to assess predictive performance. AUC, area under the curve; DEG, differentially expressed gene; FC, fold change; HRD, homologous recombination deficiency; MSRS, maximally selected rank statistics; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Differential gene analysis of HRD score in bulk dataset

Differential expression analysis between HRD score groups in the dataset was performed using the limma package in R to identify DEGs (15). Genes meeting the criteria of |log fold change (logFC)| >1 and P<0.05 were selected as significant DEGs for further analysis. Genes with logFC >1 and P<0.05 were classified as upregulated DEGs, whereas those with logFC <1 and P<0.05 were classified as downregulated DEGs. A volcano plot and heat map were subsequently generated to visualize these results.

HRD score GSEA

GSEA was performed using the clusterProfiler R package to calculate normalized enrichment scores for each gene set and to identify signaling pathways enriched in high- and low-HRD score groups within the TCGA-ESCC dataset (16). GSEA, a computational approach, was used to determine whether predefined gene sets exhibited statistically significant differences between two biological states, thereby facilitating the assessment of pathway and biological process activity changes within gene expression datasets. The gene set “c2.cp.kegg.symbols.gmt” was obtained from the MSigDB database for GSEA, enabling evaluation of the effects of high- and low-HRD score groups on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with tumors. Results meeting the criteria of P<0.05 and false discovery rate (FDR) <0.25 were retained, and seven representative pathways were selected for visualization.

HRD score mutation analysis

SNPs in patients with ESCC stratified by high and low HRD scores were analyzed using the maftools package (17). The top 20% and bottom 20% of samples based on HRD scores were selected to assess high-frequency mutated genes among patients. Summary plots and waterfall plots were subsequently generated to visualize these findings.

DDR gene prognostic screening

Univariate Cox regression analysis was conducted for batch prognostic screening to identify genes significantly associated with prognosis (P<0.05). The identified prognostic-related genes were subsequently intersected with DDR-related genes to obtain DDR-associated prognostic genes for subsequent machine learning analyses.

Multi-machine learning to achieve one-stop feature gene screening and prognostic model construction

To construct a stable prognostic model for TCGA-ESCC, prognostic machine learning analyses were performed using DDR-associated prognostic genes as follows: (I) Ten classical algorithms were integrated, including random survival forest (RSF), least absolute shrinkage and selection operator (LASSO), gradient boosting machine, survival support vector machine (Survival-SVM), supervised principal component, ridge regression, Cox partial least squares regression, CoxBoost, stepwise Cox, and elastic net. Among these, RSF, LASSO, CoxBoost, and stepwise Cox provide dimensionality reduction and variable selection functions, enabling their integration with other algorithms to generate 112 machine learning algorithm combinations. (II) The randomly assigned training cohort was used as the training cohort. These 112 algorithm combinations were applied to identify key genes and to construct prognostic models based on previously identified feature genes (18). In the test cohort, features derived from the training cohort were used to calculate risk scores (RSs) for each cohort. Models that failed to converge or reduced the number of features to zero were excluded based on the average concordance index (C-index) across the two test cohorts. The optimal consensus prognostic model for ESCC and its corresponding RS was subsequently selected. Genes included in this model were defined as hub genes.

Survival analysis and multivariate Cox regression analysis were subsequently conducted to assess the independent prognostic significance of the RS. In combination with the RS, clinical variables, including stage, age, and sex, were incorporated into additional survival and multivariate Cox regression analyses. Finally, based on the median RS, patients were stratified into high-risk (high) and low-risk (low) groups for subsequent analyses (19).

Immune infiltration analysis

The CIBERSORT algorithm, based on linear support vector regression, was applied to deconvolute the transcriptomic expression matrix and to estimate the composition and abundance of immune cells within mixed cell populations. Data were filtered to retain samples with immune cell enrichment scores greater than zero, and an immune cell infiltration matrix was generated using the CIBERSORT algorithm in conjunction with the LM22 feature gene matrix. The ggplot2 R package was subsequently used to generate group comparison plots, illustrating differences in LM22 immune cell expression in TCGA-ESCC test cohort samples stratified according to the median RS of the clinical model. In addition, the pheatmap R package was used to construct correlation heatmaps, providing a visual representation of correlation analyses among LM22 immune cells and between hub genes and LM22 immune cells.

Single-cell quality control (QC) and normalization

Prior to analysis of single-cell gene expression data, it was ensured that all unique molecular identifier (UMI) data corresponded to viable cells. Cell QC was conducted based on three covariates: UMI count depth per cell, the number of identified genes per cell, and the proportion of mitochondrial gene counts per cell. Abnormal UMI profiles may correspond to dying cells, cells with compromised membrane integrity, or doublets. For example, UMIs characterized by low count depth, a limited number of identified genes, and a high mitochondrial gene fraction may indicate leakage of cytoplasmic mRNA due to membrane damage, whereas UMIs with high count depth and many identified genes may represent doublets. Therefore, a high-count depth threshold was used to exclude potential doublets. For outlier handling, a method based on the median absolute deviation (MAD) was adopted. Cells with a MAD greater than 3 for the three QC covariates were excluded, along with cells exhibiting a mitochondrial gene proportion exceeding 20%. For bimodal filtering, the Scrublet (20) function from the Python package scanpy (21) was used to identify and remove doublets, defined as events in which two or more cells are captured within the same droplet. Doublet prediction was performed for each sample, and the expression matrix was filtered based on the predicted_doublet attribute. Normalization of the filtered single-cell data was performed using the scran R package (22), which applies a deconvolution approach based on pooled size factor estimation to reduce technical variability between cells while preserving biological variation.

Single-cell deconvolution normalization

Normalization of single-cell sequencing data was conducted using the scran R (22) package. This package applies an inverse convolution normalization approach based on pooled size factor estimation, which reduces technical variability between cells while preserving biological variation. The method involves clustering cells into similar subsets, estimating a shared size factor within each subset, and subsequently deriving size factors for individual cells through inverse convolution. Initially, scanpy in Python was used to perform log₂ transformation and Leiden clustering. The computeSumFactors function in R was then used to calculate size factors for each cell, followed by normalization using these size factors. The resulting normalized expression matrix was used for subsequent analyses.

Single-cell batch correction

The scvi-tools Python package was utilized for single-cell gene expression data analysis. This package provides models and methods based on deep learning. An autoencoder-based approach was used for batch effect correction. The autoencoder, an unsupervised neural network, enables learning of low-dimensional representations of the data while reconstructing the input. Within scvi-tools, the autoencoder was used to mitigate batch effects in single-cell data, which arise from technical variability associated with different experimental conditions. Batch correction was performed on the original data, with the batch variable specified according to sample source and grouping.

Single-cell clustering

Clustering analysis of single-cell sequencing data was conducted using the scanpy Python package. Initially, the sc.pp.neighbors function was used to construct a neighborhood graph of cells, in which distances and similarities between cells were estimated using the Uniform Manifold Approximation and Projection (UMAP) algorithm. Subsequently, the leiden function was used to perform clustering based on the Leiden algorithm, which optimizes modularity to achieve high-quality clustering. The clustering results were visualized to illustrate the distribution of cells across different clusters within UMAP coordinates.

Single-cell annotation

Initial annotation of single-cell sequencing data was performed using the celltypist Python package. Celltypist is a deep learning-based method that uses pre-trained models to enable rapid, accurate, and interpretable annotation of single-cell data derived from diverse sources and tissues. Multiple models are available and can be selected according to specific research objectives. Two immune cell annotation models with different resolutions, Immune_All_Low.pkl and Immune_All_High.pkl, were downloaded using the models.download_models function, and annotation was conducted using the annotate function. The annotation results were visualized using the celltypist.dotplot function, illustrating the distribution of cell types and their corresponding confidence levels.

Over-representation analysis (ORA) based enrichment annotation was performed using the reference dataset OmniPath, which integrates extensive prior knowledge, including PanglaoDB, a database providing marker genes for various cell types. Data from PanglaoDB were accessed through the OmniPath interface using the decoupler package. The top 5% of the most highly expressed genes in each sample were selected as the gene set of interest (S). Set operations were then conducted between S and each reference gene set to construct contingency tables (23). Subsequently, ORA was applied using one-sided Fisher’s exact tests to assess the statistical significance of overlap between S and each gene set. The resulting P values were transformed logarithmically to derive functional enrichment scores, with higher scores indicating greater enrichment.

By comparing the two annotation approaches along with the previously obtained clustering results, the ORA-based annotation method was selected as the cell type assignment for subsequent analyses.

Analysis of single-cell communication

CellPhoneDB (v2), a ligand-receptor-based inference method, was used to examine signal transmission between different cell types in single-cell data. This approach utilizes a curated database of known ligand-receptor interactions in combination with gene expression profiles derived from single-cell transcriptomic data to predict potential intercellular communication networks. Analysis was restricted to the GSM7847137 sample, which contained the largest number of cells from colorectal cancer liver metastasis cases, as this method is designed to infer intercellular communication events under steady-state conditions. Therefore, it is not intended for comparative analysis across different samples or conditions but rather for analysis within a single sample or condition. The CellPhoneDB method was implemented using the liana Python package with the CellPhoneDB framework, and the cell type exhibiting the greatest perturbation was selected for visualization.

To increase confidence in inferred ligand-receptor interactions underlying intercellular communication, a consensus analysis strategy based on multiple methods was used. Results from different ligand-receptor inference approaches were compared to identify interactions consistently predicted across methods. The rank_aggregate function from the liana Python package was utilized, which integrates outputs from multiple methods to generate a consensus ranking. This approach facilitates the identification of ligand-receptor interactions demonstrating consistent significance across methods. The integrated methods included CellPhoneDB; Connectome, a network-based approach for inferring functional connections between cells using single-cell transcriptomic data; log₂FC, which identifies ligand-receptor pairs based on gene expression differences; NATMI, a machine learning based method for predicting intercellular ligand-receptor interactions; SingleCellSignalR, which reconstructs intercellular signaling networks based on pathway data; and CellChat, which applies probabilistic graph models to infer intercellular communication patterns from single-cell transcriptomic data.

Single-cell temporal sequence analysis

Pseudotemporal analysis was performed to arrange cells along a trajectory based on temporal patterns of gene expression. From gene expression states, samples were classified into multiple cell clusters representing distinct differentiation states, enabling the construction of a lineage development trajectory. This approach allows the inference of cellular differentiation and developmental trajectories. The results of pseudotemporal analysis were interpreted in conjunction with the distribution of cell types along the trajectory and the expression dynamics of characteristic genes to determine the starting and terminal points of differentiation.

Statistical analysis

All statistical analyses were conducted using R software (version 4.3.3) and Python software (version 3.12). Differences between two groups were assessed using the Wilcoxon rank-sum test, whereas the Kruskal-Wallis test was applied for comparisons among more than two groups. Spearman correlation analysis was used to assess correlations. A P value <0.05 was considered statistically significant.

Results

Association of HRD score with genomic instability and prognosis

The HRD score was defined as the unweighted sum of LOH, telomeric allelic imbalance, and large-scale state transition scores (5,12-14,24). To stratify patients, the optimal cutoff value for the HRD score was determined using MSRS, which identified the cut point that most effectively separated patients according to survival outcomes (Figure 2A,2C).

From this cutoff, patients were stratified into high- and low-HRD score groups. The biological characteristics associated with HRD were subsequently examined. Analysis of genomic alteration patterns indicated that the low-HRD score group exhibited a significantly higher frequency of CNA, as measured by the FGA score (Figure 2B), suggesting distinct genomic instability profiles between the groups.

To further characterize molecular differences, differential expression analysis was performed. Compared with the low-HRD score group, the high-HRD score group exhibited a distinct gene expression profile, with numerous genes significantly upregulated or downregulated (Figure 2D). The top 20 DEGs are presented in the heatmap, demonstrating consistent transcriptional changes associated with high HRD scores (Figure 2E).

Finally, the clinical impact of HRD-based stratification was assessed. Kaplan-Meier survival analysis indicated that patients in the high-HRD score group exhibited significantly poorer prognosis compared with those in the low-HRD score group (P<0.05; Figure 2F). The prognostic performance of the HRD score was further assessed using time-dependent receiver operating characteristic (ROC) analysis, which demonstrated area under the curve (AUC) values >0.7 for predicting 1-, 2-, and 3-year overall survival, indicating robust predictive accuracy (Figure 2G).

Pathway enrichment analysis in high- and low-HRD score groups

To assess the impact of gene expression differences between high- and low-HRD risk groups within the TCGA-ESCC dataset, GSEA was performed to assess associations with biological processes, cellular components, and molecular functions (Figure 3A). The results demonstrated significant enrichment of gene sets associated with the following pathways and functions: Blanco Melo coronavirus disease 2019 (COVID-19) bronchial epithelial cells severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection up (Figure 3B), Blanco Melo COVID-19 SARS-CoV-2 low multiplicity of infection (MOI) infection A594 angiotensin-converting enzyme 2 (ACE2)-expressing cells up (Figure 3C), KEGG glutathione metabolism (Figure 3D), KEGG metabolism of xenobiotics by cytochrome P450 (Figure 3E), Phong TNF targets up (Figure 3F), Turashvili breast lobular carcinoma vs. ductal normal up (Figure 3G), and Wieland upregulated by hepatitis B virus (HBV) infection (Figure 3H), among other biologically relevant pathways and signaling processes.

Figure 3 GSEA analysis of HRD scores. (A) GSEA enrichment plot (mountain plot) illustrating overall pathway enrichment patterns. (B-H) Enriched gene sets identified by GSEA, including: Blanco Melo COVID-19 bronchial epithelial cells SARS-CoV-2 infection up (B); Blanco Melo COVID-19 SARS-CoV-2 low MOI infection A594 ACE2-expressing cells up (C); KEGG glutathione metabolism (D); KEGG metabolism of xenobiotics by cytochrome P450 (E); Phong TNF targets up (F); Turashvili breast lobular carcinoma vs. ductal normal up (G); and Wieland upregulated by HBV infection (H). Selection criteria were adjusted P<0.05 and FDR <0.25. P value adjustment was performed using the Benjamini-Hochberg method. ACE2, angiotensin-converting enzyme 2; COVID-19, coronavirus disease 2019; FDR, false discovery rate; GSEA, gene set enrichment analysis; HBV, hepatitis B virus; HRD, homologous recombination deficiency; KEGG, Kyoto Encyclopedia of Genes and Genomes; MOI, multiplicity of infection; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2.

Although certain pathways identified by GSEA are associated with viral infections (e.g., SARS-CoV-2 and HBV), their enrichment in the high-HRD group reflects activation of shared immune programs, including interferon signaling and antigen presentation, which represent key components of the tumor immune microenvironment in patients with ESCC (25,26).

Mutation landscape in high- and low-HRD score groups

To assess frequently mutated genes in patients within the TCGA-ESCC cohort, the top 20% (HRD top) and bottom 20% (HRD tail) of samples based on HRD scores were selected. A summary map, waterfall plot, and mutation correlation heatmap were subsequently generated.

A mutation analysis summary map for the top 20% of samples was first constructed (Figure 4A), presenting mutation classifications, mutation types, single nucleotide variant (SNV) categories, the number of variants per sample, an overview of mutation classifications, and the top 10 most frequently mutated genes. The waterfall plot (Figure 4B) illustrated mutated genes, mutation types, and the proportion of mutations across samples. A mutation correlation heatmap was generated to depict relationships between mutations (Figure 4C).

Figure 4 HRD top mutation analysis. (A) Mutation analysis summary plot of samples in the top 20% of HRD scores, including variant classification, variant type, SNV categories, number of variants per sample, summary of variant classifications, and the top 10 most frequently mutated genes. (B) Waterfall plot illustrating mutation profiles in the top 20% of samples based on HRD scores, displaying mutated genes, mutation types, and the proportion of mutations across samples. (C) Heatmap depicting mutation associations in the top 20% of samples with HRD scores, presenting correlations between mutations. DEL, deletion; HRD, homologous recombination deficiency; INS, insertion; SNP, single nucleotide polymorphism; SNV, single nucleotide variant; TMB, tumor mutational burden.

Subsequently, a mutation analysis summary map was generated for the bottom 20% of samples (Figure 5A), presenting variant classifications, variant types, SNV categories, the number of variants per sample, a summary of variant classifications, and the top 10 most frequently mutated genes. The corresponding waterfall plot (Figure 5B) depicted mutated genes, mutation types, and mutation proportions across samples, accompanied by a mutation correlation heatmap (Figure 5C) depicting relationships between mutations.

Figure 5 Analysis of HRD tail mutations. (A) Mutation analysis summary plot of the bottom 20% of samples based on HRD scores, including variant classification, variant type, SNV categories, number of variants per sample, summary of variant classifications, and the top 10 most frequently mutated genes. (B) Waterfall plot illustrating mutation profiles in the bottom 20% of samples based on HRD scores, displaying mutated genes, mutation types, and the proportion of mutations across samples. (C) Heatmap depicting mutation associations in the bottom 20% of samples with HRD scores, presenting correlations between mutations. DEL, deletion; HRD, homologous recombination deficiency; INS, insertion; SNP, single nucleotide polymorphism; SNV, single nucleotide variant; TMB, tumor mutational burden.

Multi-machine learning model based on DDR score feature genes

Following evaluation of the prognostic significance of HRD scores in esophageal cancer, a prognostic model was constructed using multiple machine learning algorithms based on DDR-related genes. The TCGA-ESCC dataset was randomly divided into training (70%) and test (30%) cohorts. Within the training cohort, univariate Cox regression analysis was conducted for batch prognostic screening to identify genes significantly associated with prognosis (P<0.05). These genes were subsequently intersected with DDR-related genes to obtain DDR-associated prognostic genes for further machine learning analyses.

A total of 112 algorithm combinations were used to construct prognostic models based on DDR-associated prognostic genes, resulting in 27 successfully trained models (Figure 6A). Although the LASSO-RSF model demonstrated the highest average C-index (0.751) across cross-validation folds, the Survival-SVM model was selected as the final model due to its superior generalization performance, achieving both the highest training C-index (0.751) and test C-index (0.712). This selection prioritized robustness against overfitting, as the Survival-SVM model maintained stable prognostic accuracy during independent validation, whereas the LASSO-RSF model exhibited a substantial decline in test cohort performance (ΔC-index =−0.15, P<0.01). Genes with non-zero coefficients in the selected model, defined as hub genes, included PARP1, MBD4, TELO2, NSMCE3, SMUG1, and BABAM1. A heatmap was generated to depict correlations among these genes (Figure 6B). Cox regression analysis was applied to the linear predictor value (RS) of the model to assess associations with overall survival and survival status. Time-dependent ROC analysis for 3-year prognosis demonstrated AUC values exceeding 0.7 at 1, 2, and 3 years, indicating strong predictive performance (Figure 6C).

Figure 6 Multi-machine learning model based on DDR score feature genes. (A) Heatmap of C-index values for 112 machine learning model combinations in the training and validation cohorts. (B) Correlation heatmap of six hub genes. (C) Time-dependent ROC curves for the Survival-SVM model. *, P<0.05 indicates statistical significance; **, P<0.01 indicates high statistical significance; ***, P<0.001 indicates very high statistical significance. AUC, area under the curve; C-index, concordance index; DDR, DNA damage response; Lasso, least absolute shrinkage and selection operator; ROC, receiver operating characteristic; RSF, random survival forest; Survival-SVM, survival support vector machine.

Validation of the Survival-SVM-based prognostic model

From C-index comparisons among 112 machine learning model combinations, the Survival-SVM algorithm was identified as the optimal approach for constructing the prognostic prediction model for patients with ESCC. Samples in the training cohort were stratified into high- and low-RS groups using the median RS as the cutoff value.

The final prognostic model was established by incorporating additional clinical variables, including age and stage, and a multivariate Cox regression forest plot and nomogram were generated (Figure 7A). The final model-derived RS was used to generate a three-factor risk diagram (Figure 7B), as well as time-dependent ROC curves (Figure 7C) and Kaplan-Meier survival curves (Figure 7D).

Figure 7 Construction of risk model. (A) Prognostic model forest map and nomogram (*, P<0.05; **, P<0.01). Gender: 0, female; 1, male. Age: 0, ≤65 years; 1, >65 years. (B) Triple plot of risk factors. (C) Time ROC curve. (D) RS Kaplan-Meier curve. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; ROC, receiver operating characteristic; RS, risk score; TCGA, The Cancer Genome Atlas.

When compared to the time-dependent ROC curve based on the HRD score, the AUC values of the newly developed prognostic model were higher at 1, 2, and 3 years.

Immune microenvironment heterogeneity between risk groups

To examine differences between high- and low-risk groups defined by the Cox model within the TCGA-ESCC dataset, correlations between 22 immune cell types and risk groups [high-risk (high) and low-risk (low)] were calculated using the CIBERSORT algorithm based on gene expression data from the TCGA-ESCC training cohort (19). A bar chart depicting the proportions of immune cell types in TCGA-ESCC was generated from the immune infiltration analysis (Figure 8A). Boxplots were constructed to compare differences in immune cell infiltration abundance between groups (Figure 8B). The analysis demonstrated statistically significant differences (P<0.05) in the infiltration levels of plasma cells and neutrophils between high- and low-risk groups. A correlation heatmap was generated to depict relationships among immune cell infiltration levels in TCGA-ESCC tumor samples (Figure 8C). Finally, an additional correlation heatmap was constructed to depict associations between key genes and immune cell infiltration levels in TCGA-ESCC tumor samples (Figure 8D).

Figure 8 CIBERSORT immune infiltration analysis. (A) Proportional bar graph of immune infiltration analysis. (B) Group comparison plot of immunoassays. (C) Heatmap of immune cell correlations. (D) Heatmap of immune cell and gene correlation. *, P<0.05. NK, natural killer.

Single-cell QC and normalization

The scanpy Python package was used to import the count matrix for all samples from the single-cell RNA sequencing (scRNA-seq) dataset GSE199619 and to construct an AnnData object. QC procedures were subsequently applied to the AnnData object (Figure 9A,9B). A total of 88,954 cells were initially included. Using a MAD-based approach, cells with a MAD greater than 3 across the three QC covariates were excluded for each sample. Cells with mitochondrial gene content exceeding 10% were removed. After these filtering steps, 76,405 cells were retained. Subsequent doublet filtering resulted in a final dataset comprising 75,407 cells (Figure 9C,9D). Normalization of single-cell sequencing data was conducted using the scran R package (22). This package applies an inverse convolution normalization approach based on pooled size factor estimation, which reduces technical variability between cells while preserving biological variation. The method involves clustering cells into similar subsets, estimating a shared size factor within each subset, and subsequently deriving size factors for individual cells through inverse convolution. The count distributions before and after normalization are shown in Figure 9E-9H.

Figure 9 Single-cell data QC and normalization. (A,B) QC of single‑cell data based on gene count, total UMI counts, and mitochondrial gene proportion. Blue dots represent retained cells, and orange dots represent filtered cells. (C,D) Doublet detection based on doublet scores. Blue dots represent cells not predicted as doublets, and orange dots represent cells predicted as doublets. (E-G) Distribution of total counts in single‑cell data: (E) before normalization, (F) after shifted logarithmic transformation, and (G) after shifted logarithmic transformation with subsequent size factor estimation and log1p transformation using Scran. (H) Count distribution after size factor estimation and log1p transformation. QC, quality control; UMI, unique molecular identifier.

Single-cell data integration, clustering, and cell-type annotation

Batch correction was performed on the original data using the variational autoencoder (VAE) framework implemented in the scvi-tools Python package (20,27), with the batch variable specified according to sample origin. Nonlinear dimensionality reduction was subsequently applied to the corrected data. UMAP visualization demonstrated distinct clustering patterns according to Leiden clustering labels (Figure 10A), substantial overlap among samples when labeled by sample origin (Figure 10B), and substantial overlap between groups (ESC, tumor group; ESN, normal group) (Figure 10C). PCA and t-distributed stochastic neighbor embedding (t-SNE) visualizations of the batch-corrected data are provided in Figure S1.

Figure 10 Single-cell data integration, clustering, and cell-type annotation. (A-C) UMAP visualization of integrated single‑cell data, colored by Leiden clustering (A), group labels (B), and individual sample origin (C). (D) Distribution of coarse cell type annotations performed by CellTypist projected onto UMAP space. (E) Confidence scores for CellTypist cell type annotations. (F) Heatmap showing the correspondence between Leiden clusters and candidate cell type annotations based on marker enrichment scores. (G) UMAP visualization colored by Leiden clustering. (H) UMAP visualization colored by ORA‑based cell type annotations. ESC, tumor group; ESN, normal group. DC, dendritic cell; ILC, innate lymphoid cell; NK, natural killer; ORA, over‑representation analysis; UMAP, Uniform Manifold Approximation and Projection.

Clustering analysis of single-cell sequencing data was performed using the scanpy Python package (22). The neighbors function was used to construct a neighborhood graph of cells, in which distances and similarities between cells were estimated using the UMAP algorithm. Subsequently, the Leiden function was used to perform Leiden clustering (28) (Figure 10G).

Cell type annotation was performed using two complementary approaches. First, the celltypist Python package (29) was applied for automated cell type annotation, classifying cell clusters into 12 cell types (Figure 10D), with confidence levels shown in Figure 10E. Second, ORA-based enrichment annotation was performed using the decoupler Python package (30) with reference to the PanglaoDB dataset (23), which provides marker genes for different cell types. The top 5% of the most highly expressed genes in each sample were selected as the gene set of interest. ORA was applied using one-sided Fisher’s exact tests to assess the statistical significance of overlap between the gene set and each reference gene set. The normalized enrichment scores were visualized to present the three most probable cell types for each cluster (Figure 10F). By integrating the results from both annotation approaches and manual verification, the ORA-based annotation was considered reliable and was therefore selected for subsequent analyses, as shown in the UMAP plot with cell type annotations (Figure 10H).

Cell communication analysis

The CellPhoneDB ligand-receptor inference method was implemented using the liana Python package (Figure 11A). Subsequently, consensus cell-cell communication analysis was performed using LIANA, integrating multiple communication inference methods (Figure 11B).

Figure 11 Cell communication analysis. (A) Bubble plot of cell-cell communication inferred using the CellPhoneDB method. (B) Bubble plot of consensus cell-cell communication analysis performed using LIANA. Point size represents confidence level, and color intensity (dark to light) indicates the strength of communication interactions. Lr, ligand and receptor.

Single-cell composition, correlation, and temporal sequence analysis

Compositional data analysis (CoDA) is a statistical approach for analyzing compositional data, which represents proportional relationships between components and the whole. This method involves transformation of compositional data into an unconstrained log-ratio space, followed by application of appropriate statistical models to assess data characteristics and differences. CoDA reduces spurious correlations and biases arising from the compositional structure of the data, thereby improving analytical reliability. Proportional distributions were initially visualized using a conventional proportion comparison plot (Figure 12A) and a stacked bar chart (Figure 12B) to assess differences in cell composition between tumor and non-tumor groups. These analyses indicated that T cells and fibroblasts were the predominant cell populations. Subsequently, CoDA was applied to obtain standardized relative abundance measures (Figure 12C).

Figure 12 Single-cell composition, correlation, and temporal sequence analysis. (A) Group comparison chart of cell type proportions. (B) Stacked bar chart of cell type proportions. (C) Relative abundance of cells under CoDA algorithm, comparing tumor and primary groups. (D) Heatmap of cell interactions. (E) Temporal trajectory UMAP map of T cell subgroups. ESC, tumor group; ESN, normal group. CoDA, compositional data analysis; DC, dendritic cell; ILC, innate lymphoid cell; UMAP, Uniform Manifold Approximation and Projection.

The scanpy Python package was used to calculate correlations between cell types in the single-cell dataset using the dendrogram function. The correlation_matrix function was subsequently used to generate a correlation heatmap illustrating relationships among cell types (Figure 12D).

Pseudotemporal analysis enables the arrangement of cells along trajectories based on temporal patterns of gene expression, allowing classification of samples into multiple cell populations representing distinct differentiation states. This approach facilitates the construction of lineage trajectory diagrams and inference of cellular differentiation and developmental pathways. Hub genes were predominantly highly expressed in T cells, which constituted the largest proportion of the cell population. A trajectory plot of the T-cell subpopulation was generated (Figure 12E).

Discussion

This study demonstrates the potential role of the HRD score in patients with ESCC, highlighting its significant associations with genomic instability, the tumor immune microenvironment, and clinical outcomes through comprehensive bioinformatics analysis (31). Multiple machine learning models were applied to analyze DDR-related genes, thereby improving the accuracy of prognostic prediction.

Among the 112 algorithm combinations tested, the Survival-SVM model was selected as the final model due to its superior generalization performance (training C-index =0.741; test C-index =0.708; Table S1), achieving an optimal balance between model complexity and predictive accuracy (32). Machine learning approaches have been increasingly applied to biomarker discovery and prognostic modeling in ESCC, as summarized in recent comprehensive reviews (33). First, a significant association between HRD score and poor prognosis was observed in patients with ESCC. The findings indicate that patients with elevated HRD scores tend to have shorter survival times, suggesting that the HRD score may serve as an independent prognostic marker. This observation is consistent with findings in other cancer types, in which higher HRD scores have been associated with poorer prognosis, likely reflecting underlying genomic instability (34,35). Similar associations have been reported in breast, ovarian, and pancreatic cancers, where HRD scores correlate with genomic instability and adverse clinical outcomes (5,14,36,37). This cross-cancer consistency indicates that the HRD score may have broader biological relevance in solid tumors. In ovarian cancer, Telli et al. demonstrated a strong association between HRD score and sensitivity to platinum-based therapies (12). Recent studies have further demonstrated the utility of machine learning-based prognostic models in ESCC, with applications ranging from predicting survival in patients receiving immunochemotherapy (38) to developing artificial intelligence algorithms for postoperative prognosis (39). The HRD score cutoff of 38 was derived using MSRS, a method designed to optimize separation of survival outcomes. Although this threshold is lower than that reported in breast cancer (≥42), it is consistent with thresholds reported in ovarian (≥33) and pancreatic (≥35) cancers, supporting cross-cancer consistency in the prognostic role of HRD (5,12). The elevated baseline genomic instability observed in ESCC (e.g., frequent TP53 mutations) may necessitate a lower threshold to adequately capture HRD-driven clinical effects. Further validation of this cutoff in larger cohorts is required, along with evaluation of its utility in predicting responses to PARP inhibitors. However, given the distinct pathological characteristics of ESCC, additional validation and refinement of HRD score application in this cancer type remain necessary. In a recent study of 96 ESCC patients, Wang et al. similarly demonstrated that HRD score was an independent prognostic factor associated with tumor stage and recurrence, further supporting the clinical relevance of HRD assessment in ESCC (40).

Second, a strong association was observed between the HRD score and genomic instability. Patients with elevated HRD scores exhibited more pronounced indicators of genomic instability, including CNA and chromosomal aberrations, supporting the potential utility of the HRD score as a marker for genomic instability assessment. The HRD score reflects deficiencies in the homologous recombination repair (HRR) pathway, which is responsible for repairing double-strand DNA breaks; when HRR function is impaired, genomic changes accumulate, resulting in increased genomic instability and enhanced tumor cell proliferation, aggressiveness, and poorer prognosis in patients with ESCC (12,41,42). The HRD score may serve as an indicator of tumor genetic variation and provide a basis for identifying patients who may benefit from therapies targeting DNA damage repair mechanisms. Gene expression analysis indicated that elevated HRD scores were associated with dysregulation of key pathways, including cell cycle regulation, DNA repair, and chromosome segregation. These molecular changes may contribute to tumor progression and adverse prognosis. In tumors with high HRD scores, impaired DNA repair capacity may result in cell cycle dysregulation, thereby promoting tumor cell proliferation and metastasis (43-45). This observation supports the potential utility of targeted therapeutic strategies directed at HRD-associated pathways, such as PARP inhibitors, although their clinical efficacy in patients with ESCC requires further validation in clinical trials. Notably, the hub genes identified in the prognostic model (PARP1, MBD4, TELO2) are mechanistically associated with ESCC pathogenesis. PARP1 is a key mediator of DNA damage repair that has been implicated in multiple cancer types. Along with the established role of PARP1, the biological functions of other hub genes further support the mechanistic relevance of the prognostic signature. MBD4 is involved in DNA mismatch repair and plays an important role in suppressing CpG site mutations, with dysfunction potentially contributing to the accumulation of driver gene changes in ESCC (46). TELO2 functions as a regulator of phosphatidylinositol 3-kinase-related kinase family members, including ATM and ATR, thereby contributing to the maintenance of genomic stability under replication stress conditions commonly observed in aggressive tumors (47). NSMCE3, as a component of the SMC5/6 complex, is involved in HRR and replication fork stabilization, processes that are frequently dysregulated in genomically unstable cancers (48). SMUG1 functions in the base excision repair pathway to remove oxidative DNA damage, and its alteration may interact with HRD-associated phenotypes to influence tumor progression (49). Notably, BABAM1 is a core subunit of the BRCA1-A complex, directly linking this gene set to HRD, the central phenotype investigated in this study (50). Collectively, these genes represent distinct yet interconnected DDR pathways, and their coordinated dysregulation reflects widespread impairment of genomic maintenance mechanisms in aggressive ESCC, thereby providing biological plausibility for the machine learning-derived prognostic model. Beyond these six hub genes, other DDR-related molecules have also been implicated in ESCC progression and immune modulation. For instance, APE1, a key DNA repair enzyme, has been validated as an independent prognostic marker for postoperative chemotherapy in ESCC, with high expression linked to immunosuppressive microenvironment remodeling (51). Similarly, KIN17, a DNA/RNA-binding protein, has been shown to promote ESCC progression through R-loop resolution and noncanonical STING pathway regulation, further highlighting the interplay between DDR and tumor immunity (52).

The HRD score demonstrates a strong association with the tumor immune microenvironment. Immune infiltration analysis revealed significant differences in plasma cell and neutrophil infiltration between HRD-defined risk groups, suggesting that genomic instability may be accompanied by remodeling of the tumor immune microenvironment. Dai et al. recently characterized chromosome instability in ESCC and found that HRD-associated tumors exhibited enrichment of lipid-associated macrophages, an immunosuppressive cell population, further supporting the link between genomic instability and immune microenvironment remodeling in this malignancy (53). The observed alterations in plasma cells may reflect changes in B-cell-mediated humoral immunity, whereas abnormal neutrophil infiltration is potentially linked to myeloid-derived immunosuppression and tumor immune evasion (54). Meanwhile, HRD may enhance immune recognition through increased neoantigen exposure and altered inflammatory signaling (55,56). These immune features may therefore have implications for responsiveness to immune checkpoint blockade, as tumors with higher mutational burden and an inflamed microenvironment have been associated with improved immunotherapy outcomes in multiple cancer types (57). As therapeutic strategies such as chemoimmunotherapy combined with radiotherapy are increasingly explored in ESCC, the identification of immune biomarkers that can aid in patient stratification and treatment response assessment has become clinically relevant (58). Beyond immune cell abundance, the functional status of immune cells may also influence antitumor immunity. Studies in solid tumors have indicated that terminally exhausted CD8⁺ T cells hold potential as prognostic and immunotherapeutic biomarkers, although their significance is modulated by tumor type and the immune microenvironment (59). However, we acknowledge that the present study did not directly assess CD8⁺ T cell exhaustion status, nor did it include immunotherapy response data. Accordingly, the current immune infiltration findings should be interpreted as providing clues for understanding HRD-associated immune remodeling, rather than as direct predictors of immune checkpoint inhibitor efficacy. Further investigations incorporating immunotherapy-treated cohorts and functional studies are required to clarify the relationship between HRD-associated immune remodeling and therapeutic responsiveness in ESCC. The enrichment of viral infection-related pathways in high-HRD tumors indicates activation of an immune microenvironment characterized by interferon signaling, consistent with evidence indicating that HRD influences immune landscapes through increased neoantigen burden (57,60). This inflamed phenotype, supported by immune infiltration findings, indicates that patients with high HRD scores in ESCC may derive benefit from combined therapeutic strategies involving PARP inhibitors and immunotherapy (26,61).

Despite these findings, several limitations should be acknowledged. First, the modest sample size (n=78) reduces statistical power and increases the risk of overfitting, particularly given that 112 algorithm combinations were evaluated using the same dataset. Although the Survival-SVM model was selected based on superior generalization performance, the imbalance between predictors and sample size remains a concern, as limited observations may not adequately capture the underlying data distribution. Testing multiple algorithm combinations without correction for multiple comparisons may introduce selection bias, potentially inflating reported performance estimates. The determination of the optimal HRD score cutoff (38) using MSRS within the same dataset may introduce optimal cutpoint bias; although the consistency of this cutoff with genomic instability, transcriptomic profiles, and immune infiltration supports its biological relevance, external validation in larger cohorts is required. Furthermore, the model was developed and internally tested using only the TCGA-ESCC dataset, with no independent external validation cohort. Publicly available ESCC cohorts present substantial challenges for direct validation, including platform heterogeneity (microarray vs. RNA-sequencing), inconsistent gene annotation, incomplete coverage of the six hub genes, and limited survival data, making careful harmonization a prerequisite for future validation efforts. Second, the retrospective design relying on TCGA data introduces inherent limitations including potential selection bias, as TCGA samples may not reflect the general ESCC population. The cohort lacks detailed treatment data (e.g., specific chemotherapy regimens, radiation protocols) and lifestyle factors (e.g., smoking history, alcohol consumption), all of which are established prognostic determinants in ESCC. Unmeasured confounding factors, including comorbidities and performance status, may affect the observed associations between HRD scores and survival outcomes. The retrospective, non-randomized design precludes causal inference, and the absence of immunotherapy response data precludes direct evaluation of HRD scores as predictors of treatment efficacy. Third, although associations between HRD scores, prognosis, genomic instability, and the tumor immune microenvironment were identified, the underlying molecular mechanisms remain to be elucidated. Experimental studies are required to determine how HRD deficiency influences immune evasion and therapeutic response in ESCC. Fourth, as a static biomarker, the HRD score does not capture dynamic genomic changes occurring during tumor progression. Future studies may improve prognostic accuracy by integrating HRD scores with dynamic biomarkers, such as circulating tumor DNA. Fifth, the six hub genes were identified through computational analyses without independent validation. The current evidence is derived from TCGA-ESCC cohort analyses and previously documented roles in DNA damage repair pathways. However, no independent transcriptomic cohort has validated their expression and prognostic performance, and no clinical specimens or functional assays have confirmed their biological roles in ESCC. Future studies incorporating independent cohorts and functional experiments are needed to systematically evaluate these candidates.

Conclusions

The application of the HRD score in ESCC demonstrates substantial potential. Assessment of genomic instability through HRD scoring enables prediction of prognosis in patients with ESCC, characterization of tumor genomic features, and evaluation of interactions within the tumor immune microenvironment. These findings provide a theoretical basis for the development of personalized treatment strategies, including the combined use of PARP inhibitors and immune checkpoint inhibitors. With further validation in multicenter studies and additional functional investigations, these results may contribute to improved prognostic assessment and therapeutic decision-making in patients with ESCC.

Acknowledgments

We would like to acknowledge the hard and dedicated work of all the staff that implemented the intervention and evaluation components of the study.

Footnote

Reporting Checklist: The authors have completed the TRIPOD+AI reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0588/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0588/prf

Funding: This study was partially funded by the Medical Science Research Project of Hebei Provincial Health Commission (No. 20250663).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0588/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, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Pennathur A, Gibson MK, Jobe BA, et al. Oesophageal carcinoma. Lancet 2013;381:400-12. [Crossref] [PubMed]
3. Cancer Genome Atlas Research Network. Integrated genomic characterization of oesophageal carcinoma. Nature 2017;541:169-75.
4. Alexandrov LB, Kim J, Haradhvala NJ, et al. The repertoire of mutational signatures in human cancer. Nature 2020;578:94-101. [Crossref] [PubMed]
5. Abkevich V, Timms KM, Hennessy BT, et al. Patterns of genomic loss of heterozygosity predict homologous recombination repair defects in epithelial ovarian cancer. Br J Cancer 2012;107:1776-82. [Crossref] [PubMed]
6. Galon J, Mlecnik B, Bindea G, et al. Towards the introduction of the 'Immunoscore' in the classification of malignant tumours. J Pathol 2014;232:199-209. [Crossref] [PubMed]
7. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell 2011;144:646-74. [Crossref] [PubMed]
8. Colaprico A, Silva TC, Olsen C, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016;44:e71. [Crossref] [PubMed]
9. Song H, Lou C, Ma J, et al. Single-Cell Transcriptome Analysis Reveals Changes of Tumor Immune Microenvironment in Oral Squamous Cell Carcinoma After Chemotherapy. Front Cell Dev Biol 2022;10:914120. [Crossref] [PubMed]
10. Rempel E, Kluck K, Beck S, et al. Pan-cancer analysis of genomic scar patterns caused by homologous repair deficiency (HRD). NPJ Precis Oncol 2022;6:36. [Crossref] [PubMed]
11. Gao J, Aksoy BA, Dogrusoz U, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal 2013;6:pl1. [Crossref] [PubMed]
12. Telli ML, Timms KM, Reid J, et al. Homologous Recombination Deficiency (HRD) Score Predicts Response to Platinum-Containing Neoadjuvant Chemotherapy in Patients with Triple-Negative Breast Cancer. Clin Cancer Res 2016;22:3764-73. [Crossref] [PubMed]
13. Birkbak NJ, Wang ZC, Kim JY, et al. Telomeric allelic imbalance indicates defective DNA repair and sensitivity to DNA-damaging agents. Cancer Discov 2012;2:366-75. [Crossref] [PubMed]
14. Popova T, Manié E, Rieunier G, et al. Ploidy and large-scale genomic instability consistently identify basal-like breast carcinomas with BRCA1/2 inactivation. Cancer Res 2012;72:5454-62. [Crossref] [PubMed]
15. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
16. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
17. 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]
18. Benowitz NL, Schultz KE, Haller CA, et al. Prevalence of smoking assessed biochemically in an urban public hospital: a rationale for routine cotinine screening. Am J Epidemiol 2009;170:885-91. [Crossref] [PubMed]
19. Sudlow C, Gallacher J, Allen N, et al. UK biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age. PLoS Med 2015;12:e1001779. [Crossref] [PubMed]
20. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol 2018;19:15. [Crossref] [PubMed]
21. Wolock SL, Lopez R, Klein AM. Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data. Cell Syst 2019;8:281-291.e9. [Crossref] [PubMed]
22. Lun AT, McCarthy DJ, Marioni JC. A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor. F1000Res 2016;5:2122. [Crossref] [PubMed]
23. Franzén O, Gan LM, Björkegren JLM. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data. Database (Oxford) 2019;2019:baz046. [Crossref] [PubMed]
24. Teslovich TM, Musunuru K, Smith AV, et al. Biological, clinical and population relevance of 95 loci for blood lipids. Nature 2010;466:707-13. [Crossref] [PubMed]
25. Zhou Y, Mo S, Cui H, et al. Immune-tumor interaction dictates spatially directed evolution of esophageal squamous cell carcinoma. Natl Sci Rev 2024;11:nwae150. [Crossref] [PubMed]
26. Xu Q, Wen Y, Huang T, et al. The distinct landscape of tumor immune microenvironment in homologous recombination deficient cancers. Biomark Res 2025;13:108. [Crossref] [PubMed]
27. Gayoso A, Lopez R, Xing G, et al. A Python library for probabilistic analysis of single-cell omics data. Nat Biotechnol 2022;40:163-6. [Crossref] [PubMed]
28. Traag VA, Waltman L, van Eck NJ. From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep 2019;9:5233. [Crossref] [PubMed]
29. Domínguez Conde C, Xu C, Jarvis LB, et al. Cross-tissue immune cell analysis reveals tissue-specific features in humans. Science 2022;376:eabl5197. [Crossref] [PubMed]
30. Badia-I-Mompel P. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinform Adv 2022;2:vbac016. [Crossref] [PubMed]
31. Zhao D, He X, Guo Y, et al. Advances in multi-omics for esophageal squamous cell carcinoma: diagnostic, prognostic, and therapeutic perspectives. Protein Cell 2026;17:491-527. [Crossref] [PubMed]
32. Wang P, Feng Z, Chen W. Establishment of prognostic prediction model based on lipid metabolism related genes in esophageal squamous cell carcinoma by machine learning algorithms. BMC Gastroenterol 2026;26:409. [Crossref] [PubMed]
33. Yang W, Duan L, Zhao X, et al. Integration of machine learning in biomarker discovery for esophageal squamous cell carcinoma: Applications and future directions. Pathol Res Pract 2025;272:156083. [Crossref] [PubMed]
34. Alexandrov LB, Nik-Zainal S, Wedge DC, et al. Signatures of mutational processes in human cancer. Nature 2013;500:415-21. [Crossref] [PubMed]
35. Al Assaad M, Hadi K, Levine MF, et al. Whole genome sequencing approach to assess homologous recombination deficiency in a pan-cancer cohort. Commun Med (Lond) 2026;6:5. [Crossref] [PubMed]
36. Liao H, Pei W, Zhong J, et al. Impact of homologous recombination deficiency biomarkers on outcomes in patients with early breast cancer: a systematic review protocol. BMJ Open 2022;12:e059538. [Crossref] [PubMed]
37. Davies H, Glodzik D, Morganella S, et al. HRDetect is a predictor of BRCA1 and BRCA2 deficiency based on mutational signatures. Nat Med 2017;23:517-25. [Crossref] [PubMed]
38. Wang Y, Wu Q. Development and Validation of an Interpretable ML Model for Survival Prediction in Unresectable ESCC with Immunochemotherapy. Cancer Manag Res 2025;17:2609-19. [Crossref] [PubMed]
39. Wang Z, Xiao Z, Zhang T, et al. Development and validation of a novel artificial intelligence algorithm for precise prediction the postoperative prognosis of esophageal squamous cell carcinoma. BMC Cancer 2025;25:134. [Crossref] [PubMed]
40. Wang Y, Ding B, Tao Y, et al. Homologous recombination deficiency score is an independent prognostic factor in esophageal squamous cell carcinoma. J Pathol Clin Res 2024;10:e70007. [Crossref] [PubMed]
41. Chabanon RM, Muirhead G, Krastev DB, et al. PARP inhibition enhances tumor cell-intrinsic immunity in ERCC1-deficient non-small cell lung cancer. J Clin Invest 2019;129:1211-28. [Crossref] [PubMed]
42. Lord CJ, Ashworth A. BRCAness revisited. Nat Rev Cancer 2016;16:110-20. [Crossref] [PubMed]
43. Ding S, Liang H, Wang J, et al. Aly as a Key Regulator of DNA Repair and Immune Evasion in Esophageal Squamous Cell Carcinoma Radioresistance. Int J Radiat Oncol Biol Phys 2026; Epub ahead of print. [Crossref]
44. Takaya H, Nakai H, Takamatsu S, et al. Homologous recombination deficiency status-based classification of high-grade serous ovarian carcinoma. Sci Rep 2020;10:2757. [Crossref] [PubMed]
45. Tateno K, Okuda K, Haruna S, et al. Esophageal Cancer Cells Exhibit Heterogeneity in DNA Double-Strand Break Repair and G2/M Checkpoint Arrest Associated With Cell Viability After Ionizing Radiation. Adv Radiat Oncol 2026;11:101987. [Crossref] [PubMed]
46. Sjöblom T, Jones S, Wood LD, et al. The consensus coding sequences of human breast and colorectal cancers. Science 2006;314:268-74. [Crossref] [PubMed]
47. Hurov KE, Cotta-Ramusino C, Elledge SJ. A genetic screen identifies the Triple T complex required for DNA damage signaling and ATM and ATR stability. Genes Dev 2010;24:1939-50. [Crossref] [PubMed]
48. van der Crabben SN, Hennus MP, McGregor GA, et al. Destabilized SMC5/6 complex leads to chromosome breakage syndrome with severe lung disease. J Clin Invest 2016;126:2881-92. [Crossref] [PubMed]
49. Kavli B, Otterlei M, Slupphaug G, et al. Uracil in DNA--general mutagen, but normal intermediate in acquired immunity. DNA Repair (Amst) 2007;6:505-16. [Crossref] [PubMed]
50. Wang B, Matsuoka S, Ballif BA, et al. Abraxas and RAP80 form a BRCA1 protein complex required for the DNA damage response. Science 2007;316:1194-8. [Crossref] [PubMed]
51. Gao H, Rao J, Xiao H, et al. APE1 mediates chemoresistance in esophageal squamous cell carcinoma by remodeling the immunosuppressive microenvironment. Front Immunol 2025;16:1689468. [Crossref] [PubMed]
52. Wei Z, Zhao N, Kuang L, et al. DNA/RNA-binding protein KIN17 supports esophageal cancer progression via resolving noncanonical STING activation induced by R-loop. Signal Transduct Target Ther 2025;10:256. [Crossref] [PubMed]
53. Dai W, Ko JM, Yu VZ, et al. Characterizing chromosome instability reveals its association with lipid-associated macrophages and clonal evolution of lymph node metastasis in esophageal squamous cell carcinoma. Cancer Lett 2025;628:217874. [Crossref] [PubMed]
54. Kandoth C, McLellan MD, Vandin F, et al. Mutational landscape and significance across 12 major cancer types. Nature 2013;502:333-9. [Crossref] [PubMed]
55. Coffelt SB, Wellenstein MD, de Visser KE. Neutrophils in cancer: neutral no more. Nat Rev Cancer 2016;16:431-46. [Crossref] [PubMed]
56. Joyce JA, Fearon DT. T cell exclusion, immune privilege, and the tumor microenvironment. Science 2015;348:74-80. [Crossref] [PubMed]
57. Benci JL, Johnson LR, Choa R, et al. Opposing Functions of Interferon Coordinate Adaptive and Innate Immune Responses to Cancer Immune Checkpoint Blockade. Cell 2019;178:933-948.e14. [Crossref] [PubMed]
58. Ma W, Baran N. Expanding horizons in esophageal squamous cell carcinoma: The promise of induction chemoimmunotherapy with radiotherapy. World J Clin Oncol 2025;16:104959. [Crossref] [PubMed]
59. Guo X, Ma S, Wang J, et al. Terminally exhausted CD8(+) T cells in solid tumors: biology, biomarker potential and translational tools for precision oncology. Front Immunol 2025;16:1709852. [Crossref] [PubMed]
60. Lei M, Gai J, McPhaul TJ, et al. Homologous recombination-DNA damage response defects increase TMB and neoantigen load, but not effector T cell density and clonal diversity in pancreatic cancer. Exp Hematol Oncol 2025;14:86. [Crossref] [PubMed]
61. Swoboda A, Nanda R. Immune Checkpoint Blockade for Breast Cancer. Cancer Treat Res 2018;173:155-65. [Crossref] [PubMed]