Innovative prognostic modeling in ESCC: leveraging scRNA-seq and bulk-RNA for dendritic cell heterogeneity analysis

This study obtained original scRNA-seq data of 9 cases of ESCC and 9 normal esophageal tissues from the National Center for Biotechnology Information Sequence Read Archive (SRA) (Project Number: PRJNA777911, URL: https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA777911↗). Furthermore, we acquired dataset GSE53624↗ from the Gene Expression Omnibus (GEO), which includes sequencing data and clinical information of 119 ESCC samples and 119 normal samples, for the purpose of training our model. Subsequently, for model validation and further analysis, we obtained sequencing data (FPKM format), clinical information, and genetic mutation details of 198 ESCA patients from The Cancer Genome Atlas (TCGA) database (URL: https://portal.gdc.cancer.gov/repository↗). In order to better match the TCGA data with GEO data, we converted the gene expression information from FPKM format to TPM format. Moreover, aligning with previous studies, we retrieved information about 10 cancer-related biological pathways (17).

In the scRNA-seq analysis, we first utilized the hg38 reference genome and employed the default parameters of Cell Ranger v.7.1.0 software for gene alignment (18). Following this, we conducted standard single-cell RNA sequencing data analysis using the ‘Seurat’ package in R. After filtering out cells with mitochondrial gene percentage (pMT) over 20%, hemoglobin gene percentage (pHB) over 1%, and those expressing fewer than 500 genes, while retaining genes expressed in at least five cells, we successfully obtained 91,810 high-quality cell samples. Data normalization was achieved using the NormalizeData function, followed by the use of the FindVariableFeatures function to identify 2000 highly variable genes, and data scaling was accomplished through the ScaleData function. Next, Principal Component Analysis (PCA) was utilized for data dimensionality reduction. To further ensure comparability between different sequencing data, we employed the ‘Harmony’ R package to eliminate batch effects among samples (19). Utilizing t-distributed Stochastic Neighbor Embedding (t-SNE) technology, we visualized the above results, and ultimately identified 30 cell clusters. The annotation of each cell cluster was conducted by combining the ‘singleR’ package (20) and CellTypist (https://github.com/Teichlab/celltypist↗), based on the expression patterns of known marker genes. Moreover, we extracted DC cells using the subset function and subjected them to similar analytical processing as mentioned above. Subsequently, the analysis of intercellular communication was conducted using the ‘CellChat’ R package (21). Following this, we performed single-cell regulatory network inference and clustering (SCENIC) analysis on DC cells using the ‘Scenic’ R package (22) Additionally, the malignancy level of DC cells was inferred using the copykat R package (23), and pseudotime analysis was conducted with the ‘monocle2’ R package (24). Finally, we calculated marker genes between three types of dendritic cells using the FindallMarkers function, with selection parameters including logfc.threshold = 1, min.pct = 0.25, only.pos = T. Furthermore, the ‘SCP’ R package was used during the data visualization process.

In the GSE53624↗ cohort, enrichment scores for each type of DCs macrophage in every sample were calculated using the Single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm (25), based on marker genes of dendritic cell. Initially, the differences in enrichment scores between normal and tumor samples were assessed, followed by dividing the ESCC tumor samples into two groups based on the median enrichment score, to conduct Kaplan-Meier survival analysis and assess survival differences between the groups.

In this phase of our study, the ‘limma’ R package (26) was initially utilized to analyze the GSE53624↗ dataset, specifically aiming to identify differentially expressed genes (DEGs) between normal esophageal tissues and ESCC samples. The selection criteria were established as a false discovery rate (FDR) less than 0.05 and an absolute log2 fold change (|log2(FoldChange)|) greater than 1. After this evaluation, we analyzed the correlation between the DEGs and markers of three types of dendritic cells, selecting genes with a correlation coefficient greater than 0.4 and a p-value less than 0.05 for subsequent univariate Cox regression analysis to identify genes significantly impacting ESCC prognosis. Thereafter, a prognostic model was constructed using the identified genes in conjunction with least absolute shrinkage and selection operator (LASSO) Cox proportional hazards regression and multivariate regression. The model’s formula is: RiskScore=∑k=1ncoef(k)×Expr(K) where cof (k) is an abbreviation for the regression coefficient, and Expr (k) represents the expression level of the genes used in model construction. Based on the median risk score, all patients were categorized into high and low-risk groups. To enhance the confirmation of the model’s prognostic effectiveness, the ‘survivalROC’ R package was applied on two datasets: the training set GSE53624↗ and the validation set TCGA-ESCC. This was done to construct Receiver Operating Characteristic (ROC) curves and compute the Area Under the Curve (AUC). The aim of this procedure was to ascertain the risk model’s precision and evaluate its viability for clinical application.

Initially, both univariate and multivariate Cox regression analyses, integrating clinical features (such as age, alcohol consumption, lesion location, gender, and pathological staging) along with risk scores, were employed to identify factors significantly affecting ESCC prognosis. Utilizing these analysis outcomes, a nomogram was developed using the ‘rms’ package (27, 28) for predicting the survival probability of patients at 1, 3, and 5 years. To confirm the diagnostic and prognostic capabilities of the nomogram, decision curves and calibration curves were generated. These methods were utilized to assess the clinical benefits of the model at different risk thresholds and to evaluate the concordance between predicted survival probabilities and actual observed survival probabilities, thereby ensuring the accuracy and practicality of our model in clinical applications.

Mutation data of ESCC patients were retrieved from the TCGA database using the ‘TCGAbiolinks’ R package and then uniformly decompressed. Using the read.maf function of the ‘maftools’ R package, mutation data and clinical information were read into MAF files. The oncoplot function was employed to create a heatmap that combines clinical and mutation information, showcasing the mutation profiles of high and low-risk groups, and the somaticInteractions function was used to analyze the co-mutation patterns of hub genes and the top 10 most frequently mutated genes in TCGA-ESCC.

Initially, the correlation between each Hub gene and common signaling pathways was assessed, followed by the application of the ‘GSEA’ algorithm (29) to analyze significantly enriched Gene Ontology (GO) pathways (30) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (31) between high and low-risk groups.

Initially, we assessed the expression differences of immune checkpoint-related genes and major histocompatibility complex genes between high and low-risk groups. Subsequently, we validated the model’s prognostic and immunotherapy efficacy prediction capabilities in the IMvigor210 and GSE78220↗ cohorts, respectively.

Using TRIzol reagent from Thermo Fisher Scientific (Waltham, MA, USA), total RNA was extracted from ESCC cells and tissues. Following the protocol provided by the manufacturer, this RNA was then reverse-transcribed into complementary DNA (cDNA) using the RevertAid™ First Strand cDNA Synthesis Kit, also from Thermo Fisher Scientific. Quantitative real-time PCR (qRT-PCR) analyses were conducted using a Takara Bio’s SYBR Green PCR kit (Otsu, Japan) on Thermo Fisher Scientific’s StepOne Real-Time PCR system. For quantifying the levels of gene expression, the 2-△△CT method was employed.

Following approval by the Ethics Committee, 16 paraffin-embedded ESCC tissue sections, including tumors and corresponding peritumoral tissues, were obtained from the Pathology Department of the Second Hospital of Hebei Medical University. The immunohistochemical staining commenced with incubating the slides at 60°C for four hours to fix the tissue, followed by deparaffinization in xylene and rehydration through a graded series of alcohol. Subsequently, antigen retrieval was conducted in citrate buffer under high pressure for 20 minutes before the slides were cooled to room temperature and washed with PBS. Endogenous peroxidase activity was quenched with a 30-minute treatment of 3% hydrogen peroxide, followed by further PBS washes. The sections were then blocked with PBS containing 10% bovine serum. Primary antibody incubation proceeded overnight, after which the slides were washed with PBS, treated with secondary antibody for one hour, and washed again. Color development was achieved by applying DAB chromogen and timing the reaction carefully. Following this, the sections were rinsed with water, counterstained with hematoxylin for three minutes, washed again, and finally subjected to a dehydration, clearing, and mounting process. This sequence of steps completed the staining of the tissue sections for subsequent microscopic observation and analysis.

Firstly, strict quality control measures were implemented on the data, resulting in the elimination of low-quality cells that did not meet the requirements (Supplementary Figures 1A, B). Subsequently, during the analysis phase, batch effects in the original dataset were identified, a critical factor when explaining differences between cells (Supplementary Figure 1C). Following this, the non-biological effects were effectively mitigated using the Harmony algorithm, thereby optimizing sample distribution and minimizing the impact of experimental condition variability on the analysis results (Supplementary Figure 1D). Thereafter, following data normalization and dimensionality reduction clustering, the single-cell data were visually depicted in a t-SNE plot. Figure 1A, depicting the cell distribution by sample origin, reveals that cells from different samples are evenly mixed in space, unaffected by significant batch effects. Continuing with the analysis, 30 cell subpopulations were identified from 91,810 high-quality cells, representing various cell states in the TME (Figure 1B). Through the integration of manual and algorithmic annotations, these subpopulations were further refined into 10 major cell clusters (Figure 1C). The functional annotation of these cell clusters highlighted that dendritic cells (DCs) are primarily associated with pathways such as positive regulation of leukocytes, allograft rejection, and T-cell regulation in pancreatic cancer (Figure 1D). Figure 1E presents the top 10 most highly expressed genes in each cell cluster, further elucidating the biological characteristics of each cluster. Furthermore, the expression patterns of marker genes in Figure 1F, Supplementary Figure 1E illustrate the precision of our clustering. Ultimately, Figure 1G, depicting the distribution of cell tissue origins, further confirms the uniform mixture of various cell types in ESCC, apart from epithelial cells, emphasizing the cellular heterogeneity in the TME.

Our cellular communication analysis revealed a marked escalation in both the number and intensity of intercellular communications within tumor tissues, as opposed to normal tissues, with Figure 2A visually supporting this finding. This trend likely indicates a heightened level of interaction between tumor cells and adaptive responses within the TME. Elaborating on these insights, the heatmap in Figure 2B uncovers noticeable disparities in communication probabilities among different tumor cell subpopulations, particularly highlighting the intensified interactions between dendritic and epithelial cells in tumor contexts. Delving into signal pathway analysis, illustrated in Figure 2C, we discerned notable variances in the frequency and intensity of cellular communications between tumor and normal tissues. This analysis brought to light the prominent activation of pathways such as GALECTIN, MHC−I, and MHC−II within tumors, hinting at their potential pivotal roles in tumor evolution. Following this trajectory, our subsequent examination, as demonstrated in Figure 2D, brought to the fore significant alterations in ligand-receptor pairs within tumor tissues. These specific molecular interactions suggest critical regulatory roles in intercellular communication, thereby unveiling potential new molecular targets for therapeutic intervention. A deeper dive into the GALECTIN and MHC-I pathways’ roles in cellular communication revealed distinct dynamics. The GALECTIN pathway, predominantly initiated by dendritic cells within the TME and targeting T lymphocytes, is illustrated in Figure 2E. Conversely, the MHC−I pathway, prominently active in T lymphocytes and functioning as a signal emitter in dendritic cells, is depicted in Figure 2F. Intriguingly, the comparative analysis of key molecules in the GALECTIN pathway between tumor and normal samples, as seen in Figure 2G, did not indicate significant expression differences. This suggests that their role in the TME might be mediated through mechanisms beyond mere expression level alterations. Concluding our analysis, Figure 2H elucidates the pivotal receptor-ligand pairs in the MHC-I pathway, particularly highlighting the dominance of HLA-B - CD8A, HLA-A - CD8A, and HLA-C - CD8A in tumor communication. This finding accentuates their significance in the tumor immune landscape, potentially positioning them as central targets for future therapeutic exploration. These comprehensive analytical results not only demonstrate the intricacies and dynamic nature of cellular communication within the TME but also pinpoint potential key regulatory nodes in tumor progression. This provides novel insights for developing ESCC treatment strategies and lays a solid foundation for future advancements in precision medicine research.

In order to gain a deeper understanding of DC subgroups, we initially utilized the subset function within the Seurat package, specifically focusing on isolating DC cells for enhanced scrutiny. This approach led to the identification of 3,877 high-quality DC cells following the crucial steps of data normalization and dimensionality reduction clustering. These cells were distinctly visualized in a t-SNE plot, where they were further subdivided into 25 unique subgroups, as illustrated in Figure 3A. Subsequently, by employing specific marker genes for cell annotation (32), we meticulously differentiated these cells into three DC subtypes: cDCs, tDCs, and pDCs. Each subtype, with its unique features, is concisely represented in Figure 3B. Furthermore, the analysis of the sample origin distribution of these DC cells, showcased in Figure 3C, revealed a uniform distribution across various samples. This uniformity is indicative of minimal batch effect influence, thereby reinforcing the reliability of our clustering approach. Advancing further into our investigation, the expression analysis of marker genes was conducted, as depicted in Figure 3D. This analysis not only validated the precision of our cell clustering but also enriched our understanding of the cellular characteristics. Moreover, the functional enrichment analysis of the three DC subgroups, detailed in Figure 3E, unveiled distinct functional pathways associated with each subgroup. In Figure 3F, we demonstrate the significantly overexpressed and underexpressed differential genes across three dendritic cell types. We discovered that cDC cells are predominantly involved in pathways like ‘response to molecule of bacterial origin’ and ‘positive regulation of cytokine production’, while tDC cells are linked to the ‘regulation of leukocyte proliferation’. In contrast, pDC cells are primarily associated with ‘ribonucleoprotein complex biogenesis’. These insights elucidate the diverse and significant roles that different DC subgroups play within the immune system, further contributing to our comprehensive understanding of their functionalities in various physiological contexts.

To elucidate the evolutionary relationships among the three types of DCs, we embarked on a pseudotime analysis using the Monocle2 package. This analysis revealed that DC cells progress through nine distinct differentiation states. Intriguingly, cDCs were predominantly found in the early stages of differentiation, while tDCs appeared to advance towards later stages. pDCs, however, were present in both early and late stages, as illustrated in Figures 4A-C. This finding points towards the dynamic and complex nature of DC cell differentiation. Further examining the number and percentage of DC cells among different patients, as shown in Figures 4D, E, we uncovered significant heterogeneity. This variability suggests that DC cells may assume diverse roles in different individuals, highlighting the complexity of their functions in the immune response. Subsequently, we compared the cell cycle distribution of the dendritic cell subpopulations, presented in Figures 4F, G. This comparison revealed that pDCs had a higher proportion of cells in the G2/M and S phases, indicating a more active cell cycle status, whereas cDCs and tDCs predominantly occupied the G1 phase. Delving deeper, our analysis, as demonstrated in Figure 4H, focused on the expression patterns of various transcription factors across dendritic cell types. A detailed examination unveiled distinct expression trends for transcription factors such as IRF1, NFKB1, and RELB within the DC subgroups. Notably, IRF1 expression was markedly higher in cDC cells compared to pDC cells. NFKB1 exhibited relatively high expression in both cell types, albeit more pronounced in cDCs. Conversely, an increased expression of RELB was observed predominantly in pDC cells. Additionally, POLR2A was found to be highly expressed across all dendritic cell types, particularly in cDC cells, which might indicate its broad and pivotal role in dendritic cell functions. These results shed light on potential transcriptional regulatory differences between dendritic cell subgroups, which are integral to comprehending their distinct functions in immune responses. The unique expression patterns of these transcription factors likely mirror the specialized roles of dendritic cells in immune surveillance, antigen presentation, and inflammatory responses, thereby contributing to our understanding of the intricate dynamics within the immune system.

In the GSE53624↗ dataset, we conducted a thorough analysis of the differences in the behavior of cDC, tDC, and pDC dendritic cells between normal esophageal tissue and ESCC tissue. Figures 5A-C revealed significant differences in the enrichment scores of the three types of DCs between normal and cancerous tissues. In particular, tDC and pDC showed significantly higher enrichment scores in ESCC tissues compared to normal tissues, reflecting their possible activated state in the TME. Further survival analysis using the Kaplan-Meier curve method explored the correlation between these enrichment scores and patient prognosis. Figures 5D, E showed that high levels of cDC and tDC enrichment fractions were positively correlated with worse prognosis in ESCC patients. While the correlation for cDC did not reach statistical significance (p=0.31), tDC exhibited a statistically significant correlation (p=0.03).In stark contrast, high levels of pDC enrichment scores were significantly correlated with better prognosis in patients (Figure 5F). These results revealed the unique roles of different dendritic cell subgroups in the pathogenesis of esophageal cancer and their potential mechanisms affecting patient prognosis. The significant association of tDC and pDC enrichment scores with patient survival probabilities highlights their importance as potential biomarkers in the TME and in future clinical decision-making. These findings provide valuable molecular targets for future therapeutic strategies targeting dendritic cells and offer new perspectives for clinical prognosis assessment.

In our in-depth molecular characterization of DC subgroups, we employed the Copykat algorithm to estimate the benign or malignant state of each cell. t-SNE visualization results indicated a clear distinction between malignant and benign cells among the 3,877 cells analyzed, with 1,141 malignant and 2,536 benign cells identified (Figure 6A). Comparing the proportion of malignant cells across three different types of dendritic cells revealed variability in malignancy rates within each cell type, potentially reflecting their distinct functions or pathological states in the TME (Figure 6B). Specifically, cDCs had the lowest proportion of malignancy (17%), followed by pDCs (29%), and tDCs had the highest (48%). The activity of ten tumor-related signaling pathways was assessed using GSVA, with the results presented in a heatmap format (Figure 6C). This analysis demonstrated a more significant association of malignant cells with these pathways, revealing specific correlations between each dendritic cell type and certain signaling pathways, providing clues to their potential roles in tumor development. Further analysis of functional state differences between benign and malignant cells within each subgroup showed significant scoring differences in multiple tumor-related signaling pathways for cDCs (Figure 6D), tDCs (Figure 6E), and pDCs (Figure 6F). This disparity may reveal different mechanisms by which malignant dendritic cells promote tumor growth and modulate tumor immune responses. Integrating these findings, we conclude that different dendritic cell subgroups exhibit distinct signaling pathway activation patterns and functional states in malignant tumors, which is crucial for understanding their roles in the tumor immune microenvironment.

In the process of conducting an in-depth analysis of the GSE53624↗ dataset, we initially calculated the differentially expressed genes (DEGs), selecting genes with an absolute log fold change (|lgFC|) greater than 1 and a p-value less than 0.05 (Figure 7A). Next, we analyzed the correlation between tDC and pDC marker genes and these DEGs, selecting genes with a correlation coefficient greater than 0.4 and a p-value less than 0.05. Based on these genes, univariate Cox regression analysis identified 46 genes with significant prognostic value for esophageal squamous cell carcinoma (ESCC), of which 28 were risk factors and 18 were protective factors (Figure 7B). To further narrow down the gene pool, we employed Lasso regression analysis and selected 13 genes at the optimal cutoff value (lambda=0.0735) (Figure 7C). Ultimately, we constructed a prognostic model comprising six genes using stepwise multivariate Cox regression, including Coiled-Coil Domain Containing 50 (CCDC50), ETS Variant 5 (ETV5), Neuralized E3 Ubiquitin Protein Ligase 3 (NEURL3), Lysosomal-Associated Membrane Protein Family Member 5 (LAMP5), Complement Receptor 2 (CR2), and Serine Dehydratase (SDS). The model formula is: Risk = 0.684 * CCDC50 + 0.221 * ETV5 + 0.158 * LAMP5 + 0.203 * NEURL3 - 0.62 * SDS - 0.15 * CR2. Kaplan-Meier survival analysis was used to evaluate the prognosis of patients in high-risk and low-risk groups, with results indicating a significantly worse prognosis for the high-risk group. In the training set, the prognostic model demonstrated good predictive performance with area under the curve (AUC) values of 0.81, 0.76, and 0.74 for 1-year, 3-year, and 5-year survival predictions, respectively (Figure 6D). Additionally, we validated the prognostic assessment and diagnostic capabilities of the model in the TCGA-ESCA dataset, obtaining similar positive results (Figure 7E). The risk score distribution chart (Figure 7F) clearly differentiated between high-risk and low-risk group patients, while the survival status distribution chart (Figure 7G) showed a noticeably higher number of deaths in the high-risk group compared to the low-risk group. Lastly, a heatmap (Figure 7H) detailed the expression patterns of the six core genes in the training set, further supporting their critical roles in esophageal cancer. In summary, the gene signature model we constructed is not only statistically significant but also holds potential clinical application value. It provides important molecular markers for the personalized treatment of esophageal cancer patients.

To comprehensively assess the prognosis of ESCC patients and its applicability in clinical decision-making, we conducted a series of statistical analyses and model validations. Initially, we identified factors significantly affecting the prognosis of ESCC patients through univariate and multivariate Cox regression analyses, combined with clinical features and risk scoring (Figures 8A, B). The analysis indicated that pathological staging (Stage) and risk scoring are important factors affecting prognosis. Based on these findings, we constructed a nomogram (Figure 8C) combining risk scores and pathological staging to predict 1-year, 3-year, and 5-year survival probabilities. This tool aims to provide a quantitative method to assist physicians in treatment decision-making, estimating patients’ survival probabilities by calculating a total score for each patient. To validate the predictive accuracy of this nomogram, we plotted calibration curves (Figure 8D) to assess the concordance between the predicted and actual survival probabilities. The calibration curves demonstrated the model’s accuracy in predicting 1-year, 3-year, and 5-year survival, providing an intuitive validation of the model’s predictive capability. From the calibration curves in Figure 8D, our nomogram prediction model demonstrated good accuracy in forecasting 1-year, 3-year, and 5-year survival probabilities of ESCC patients in the training set. The calibration curves closely followed the ideal line, indicating a match between predicted survival probabilities and actual observed survival rates, confirming good calibration performance of the model at different time points. In the decision curve analysis, when applying our nomogram prediction model, particularly at moderate threshold ranges, the model showed higher net benefits, indicating strong clinical applicability in differentiating medium to high-risk groups of ESCC patients (Figure 8E). This result emphasizes the important value of the model in accurately stratifying patient risk and assisting in the formulation of corresponding treatment strategies. Overall, our analysis revealed a powerful prognostic assessment tool that combines clinical features and biomarker scoring of patients, effectively predicting the survival probabilities of ESCC patients. The model’s accuracy and clinical utility in predicting the prognosis of esophageal cancer patients were confirmed through the validation by calibration curves and decision curve analysis.

An in-depth bioinformatics analysis of gene expression data from ESCC patients revealed significant differences in biological processes and metabolic pathways between different risk groups, offering new perspectives on the molecular characteristics and pathological mechanisms of these patient groups. As shown in Figures 9A, B, we first analyzed potential signaling pathways related to the genes used in model construction. A total of 30 pathways were significantly related to these genes, including the B cell receptor signaling pathway, primary immunodeficiency, colorectal cancer, etc., which play key roles in the biology and clinical characteristics of tumors. Subsequent GO and KEGG enrichment analyses further emphasized the molecular-level differences between different risk groups. As illustrated in Figure 9C, GO enrichment analysis showed that gene expression features of patients in the high-risk group were primarily focused on processes related to tumor invasiveness and metastasis, such as external encapsulating structure organization and collagen fibril organization. This finding suggests that tumors in these patients may have a stronger tendency for invasiveness and deterioration. In contrast, the low-risk group showed enrichment in processes related to immune response, such as keratinization and T cell receptor complex (Figure 9D), possibly reflecting a stronger immune response and lower tumor invasiveness in these patients. KEGG enrichment analysis results (Figures 9E, F) further highlighted these differences. The high-risk group was significantly enriched in ecm receptor interaction and pathways in cancer, suggesting that the tumor microenvironment might be more conducive to tumor growth and spread. Conversely, the low-risk group was mainly enriched in drug metabolism-cytochrome P450 and linoleic acid metabolism pathways, revealing unique characteristics of these patients in drug and lipid metabolism. In summary, our analysis disclosed marked differences in molecular traits and biological processes between high-risk and low-risk groups among ESCC patients.

Initially, we compared the expression differences between immune checkpoint-related genes and major histocompatibility complex genes in high and low-risk groups (). The results indicated that these genes, including classic molecules like PDCD1 and HAVCR2, were predominantly expressed at higher levels in the low-risk group. Next, we compared the survival differences between high and low-risk groups in the IMvigor210 cohort. Consistent with previous analyses, there were significant prognostic differences between the different risk groups (). Additionally, in the group with better immunotherapy responses (CR+PR group), there were lower risk scores (). Concurrently, the proportion of patients with poorer immunotherapy responses (PD+SD group) was significantly lower in the low-risk group (). Lastly, we validated the prognostic ability of the model in another immunotherapy cohort, obtaining similar results (). 1 1 1 1 1

Given that our analysis was primarily based on bioinformatics, it is possible that certain biases existed. To confirm the accuracy of our analysis, we initially downloaded the expression data of six model genes across various cancers from the Timer2.0 database. We then conducted qPCR and IHC validations of these model genes’ expression levels. As shown in Figure 10, our qPCR analysis (involving 8 pairs of ESCC patients and their corresponding peritumoral tissues) revealed significant overexpression of the six model genes in tumor tissues. While the pan-cancer analysis did not demonstrate statistical significance for CCDC50 and CR2, the trend of gene overexpression was consistent. This discrepancy might be attributed to the genetic differences between Eastern and Western populations. We further validated the expression of CCDC50 and CR2 at the protein level using IHC, which confirmed their significant overexpression in tumor tissues (Figure 11). These analyses substantiate the accuracy and credibility of our bioinformatics finding.