Identification of esophageal cancer tumor antigens and immune subtypes for guiding vaccine development

Introduction

The primary clinical treatments for esophageal cancer (ESCA) are varied and encompass surgical resection, radiotherapy, chemotherapy, and targeted therapy, as well as agents that target surface antigens, signaling pathways, and immune checkpoints. However, the efficacy of these treatments is often limited due to late diagnosis. Consequently, it is crucial to identify early diagnostic indicators of ESCA to enhance patient prognosis (1). To enhance the management of ESCA and improve patient survival, it is essential to develop individualized treatment strategies that are informed by an understanding of both host and tumor genomics, as well as approaches that leverage the host immune system (2). Cancer vaccines predominantly utilize tumor-associated antigens and tumor-specific antigens to activate the patient’s immune system. In principle, these vaccines are designed to elicit targeted cellular and humoral immune responses that can inhibit tumor growth and ultimately lead to the eradication of tumor cells (3).

In recent times, messenger RNA (mRNA) vaccines have emerged as a potential substitute for traditional immunotherapy in the domain of anticancer treatment. Furthermore, through the progress made in mRNA synthesis, the chemical alteration of mRNA, and delivery system technology, substantial enhancements have been achieved in the stability and translational efficacy of mRNA (4).

Tumor immunotherapy and other biological treatments have the potential to be applied in controllable immunogenicity (5). In the clinical context, mRNA vaccines can be economically synthesized through the advancement of rapid in vitro transcriptional reactions, enabling the design and modification of mRNA sequences for encoding various pathological antigens. mRNA vaccines are highly suitable for targeting antigens specific to tumors. Additionally, there is increasing evidence supporting the use of immunophenotyping for assessing the combined immune status of a tumor and its immune microenvironment, which significantly impacts therapeutic response and vaccination potential. In clinical settings, mRNA vaccines can be efficiently manufactured due to the development of prompt in vitro transcriptional responses. Finally, the design and modification of mRNA sequences allow for the encoding of any pathological antigen (6).

The objective of this investigation was to identify potent antigens in ESCA for the advancement of mRNA vaccines and to subsequently distinguish immune subcategories of ESCA for the creation of an immune landscape for the identification of suitable candidates for vaccination. The immune landscape of ESCA was characterized by examining the dispersion of pertinent genetic characteristics in each patient. The findings indicated a multifaceted tumor immune microenvironment in each ESCA patient, offering a conceptual foundation for the formulation of mRNA vaccines and the identification of suitable candidates for vaccination. We present this article in accordance with the STREGA reporting checklist (available at https://jtd.amegroups.com/article/view/10.21037/jtd-2025-233/rc).

Methods

Sample source

ESCA-related expression data and clinical data were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) to remove overall survival (OS) deletion samples, resulting in 196 data samples, comprising 185 cancer samples and 11 normal control samples. GSE19417 was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), including 76 samples, with 52 of these being ESCA tumor samples. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Acquisition of overexpressed genes

Differentially expressed genes between cancer samples and normal samples were obtained using the “limma” R package (The R Foundation for Statistical Computing). The thresholds were fold change >1.5 and a P value <0.05, and visualization was conducted via the R packages “ggplot” (v.3.3.2) and “Ideogram” (v.02.2). Significantly different genes were then extracted as overexpressed genes and used for subsequent analysis.

Mutation analysis for identifying mutated genes

The sample mutation gene information was extracted from the mutation MAF file of TCGA-ESCA, and data processing and visual analysis were implemented using the R package “maftools” and “Upset”. High-frequency mutated genes were also summarized for subsequent analysis (7).

Identification of overexpressed and mutated genes as antigen candidates

The intersection of overexpressed and mutated genes (containing copy number variants) was used to identify potential antigenic candidates, and then enrichment analysis of the gene set was performed using the R package “clusterProfiler” (v.3.14.3); all of the enrichment analysis thresholds were a P value <0.05 and a q value <0.1.

Survival analysis for identifying potential prognosis-related antigens

Batch single-factor Cox regression analysis of TCGA-ESCA data was conducted based on OS data and progression-free interval (PFI) data via the R packages “survival” and “survminer”. Genes significantly associated with OS and disease-free interval (DFI) were screened according to a threshold P value <0.05. Intersecting genes from the prognostic analysis were taken as significant prognosis-associated antagonistic genes, and Kaplan-Meier (KM) survival curves were plotted, with the mean being used as the threshold for delineating the high- and low-expression sample groups (8).

Immunological subtype identification

Based on the previously obtained expression profiles of potential antigenic genes and genes associated with immune antigen-presenting cell (APC) infiltration, the ESCA cancer samples were subjected to unsupervised clustering via the R package “Consensusclusterplus”, with the KM clustering algorithm and Euclidean distance. Survival analysis was performed and KM curves were plotted. The same clustering algorithm and distances were used to perform unsupervised clustering of the GSE19417 samples, and after subtypes consistent with TCGA data were obtained, survival analysis was performed on the different subtypes and KM curves were plotted (9).

Analysis of intersubtype mutations and copy number variation (CNV)

Changes in CNV in different subgroups were analyzed by Genomic Identification of Significant Targets In Cancer (GISTIC) using the masked copy number segment of TCGA-ESCA samples. Finally, visualization was achieved using the R package “maftools”. Based on the ESCA mutation MAF file and subsequent combination with the subtype grouping, the mutation landscapes of the samples between subtypes were also plotted using the R package “maftools”.

Differential analysis of immune activity cycle scores and other factors among the subtypes

The immune cycle activity data of TCGA-ESCA samples were downloaded, and then the difference in immune activity, homologous recombination deficiency, neoantigen load and chromosomal instability, and stemness index among different subtypes was visualized through boxplots. The Wilcoxon test was applied to determine the significance of differences (10).

Comparison of immune checkpoints and immunogenic cell death regulators between subtypes

Immune checkpoint-related genes and immunogenic cell death regulators were collected, and based on TCGA and GEO data, the expression differences of the related genes between different subtypes were statistically calculated (11).

Comparison of immune cell infiltration and tumor purity between subtypes

Stromal score, immune score, estimate score, and tumor purity were calculated for TCGA and GEO tumor samples with the R package “ESTIMATE”. The differences in correlation scores between subtypes was determined, from which box plots were drawn, and significant differences were tested via the Wilcoxon test. The CIBERSORT deconvolution method allows for the inference of the proportion of 22 immune cells within a sample from the characteristic matrix of different cells. Immune cell proportions were calculated for all samples in the TCGA and GEO databases using the R package “CIBERSORT” (v.103), and then heatmap and box plots were plotted to depict differences in immune cell proportions between high-and low-risk groups, with the Wilcoxon test being used to test for significant differences (12).

Comparison of subtypes with ESCA’s existing tumor biomarkers and existing pancancer immune subtypes

Existing ESCA tumor biomarkers include m6A Reader HNRNPA2B1 (13), CDK4/6 (14), and FABP3 (15). The expression of these markers is closely related to the prognosis of ESCA, and so the differences in the expression of related marker genes between different subtypes were calculated. The classification data of pancancer immune subtypes were downloaded, and the differences in the distribution of existing pancancer immune subtypes in the subtypes classified in this study were counted (16).

Characterization of the immune landscape

Based on the counts matrix of TCGA-ESCA immune-related genes, trajectory analysis was performed based on cancer samples using the R package “monocle” and the dimensionality reduction algorithm DDRTREE, with the maximum group score set to 3; the distribution of immune subtypes of the samples under different trajectory classifications was then visualized with a tree diagram. As there could arise a situation in which an immune subtype is distributed in two different immune trajectories, principal component 1 (PC1) and PC2, obtained based on the trajectory analysis, were used to subdivide the samples for the same immune subtype appearing in different branches; subsequently, the subdivided samples were displayed in a tree diagram. Survival analysis was performed on the grouped data and KM curves plotted to verify the extreme survival differences of the samples (17).

Weighted gene coexpression network construction and screening of prognostic modules

The purpose of weighted gene coexpression network analysis (WGCNA) is to identify co-expressed gene modules, clarify gene networks, determine the relationships between phenotypes, and mine the core genes in the network modules. WCGNA uses a soft threshold–based determination method, the value of which directly affects the construction of the modules and the delineation of genes within the modules. A WGCNA of immune-related genes was performed using the R package “WGCNA” to obtain two different modules based on a soft threshold of 0.8. After this, the signature genes of the modules were calculated, survival analysis (OS) of the modules was performed using the signature genes with P<0.05, and correlation scoring was conducted with the previously obtained immune-related first PCs. Pearson correlation coefficient analysis with P<0.05 as the threshold was implemented to obtain prognosis-related or immune-related modules and to determine the differences in the distribution of eigengenes of the related module between subtypes, with the Wilcoxon test being used to identify the significant differences (18,19).

Enrichment analysis within key modules

Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis was performed on signature genes within survival-significant or immune-related modules, and bubble maps were plotted after enriched pathways with a P value <0.05 and a q value <0.2 were obtained (20,21).

Identification of hub genes in modules and construction of risk models

The module most relevant to the OS or immune PC was selected, and the hub genes were screened based on the correlation between the module genes and the module eigengenes determined by a correlation test with a P value <0.01 being used as the threshold. Further multifactor Cox regression analysis, with P<0.05 as the cutoff, was performed on the hub genes to screen for prognosis-related genes, and in accordance with the Cox regression model factors, risk models were constructed. The risk score was calculated as follows: RiskScore_i = _ji ´ θ_j,where i denotes the sample, j denotes the gene, exp is the expression of the corresponding gene, and θ is the regression coefficient of the corresponding gene in the multifactor regression results. The risk score is obtained by multiplying the expression of the corresponding key gene in each sample by the coefficient of the corresponding gene. Risk model genes can be used as biomarkers to predict patient prognosis and to identify patients suitable for immunological drugs. The risk score of the sample was then used as the prediction result of the model, the survival data were used to calculate the area under the curve (AUC) value of the model score, and the receiver operating characteristic (ROC) curve was plotted.

Statistical analysis

The Wilcoxon test was employed to compare data between two groups, while the Kruskal-Wallis test was utilized for comparisons involving three or more groups. The univariable Cox regression analysis was conducted to define the hazard ratio (HR) for genes based on mRNA expression data. Survival analysis was performed using the Kaplan-Meier method and validated with the log-rank test. A P value of less than 0.05 was deemed statistically significant. All statistical analyses were carried out using R software.

Results

Screening of overexpressed genes

Based on the gene expression profiles of ESCA cancer samples and normal samples, significant differentially expressed genes [fold change (FC) >1.5; adjusted P value <0.05] were calculated (Figure 1A,1B), the expression matrices of 7,436 upregulated genes were extracted for subsequent analysis, and the location of differentially expressed genes in chromosomes was determined (Figure 1C).

Figure 1 Screening of overexpressed and mutant gene sets. (A) Volcano map of differentially expressed genes. Red dots represent significantly differentially upregulated genes, and blue dots indicate significantly differentially downregulated genes. (B) Heatmap of differentially expressed genes. The color of the heatmap changes from blue to red as the number of expression increases. (C) Location of differential genes on chromosomes. Green indicates low expression, and red indicates high expression. (D) Waterfall plot of mutated genes. (E) Grouping statistics for change in gene CNV values. (F) The distribution of the top 10 genes in terms of CNV change values in the different subgroups. (G) Statistics on the number of mutated genes grouped together. (H) The distribution of the top 10 genes with the highest number of mutations in the different subgroups. CNV, copy number variation.

Identification of mutant genes in ESCA via mutation analysis

Information on the mutated genes in the samples was extracted from the MAF compiled file of TCGA-ESCA. Of the 183 cancer samples, 172 were mutated, with TP53 having the highest mutation rate (69%). A total of 9,064 genes were mutated (Figure 1D), with the highest proportion of genes having a mutation number of 1 and the lowest number of genes having a mutation number of >20 (Figure 1E). The 10 genes with the largest CNV changes were ANO1, CCND1, CTTN, FADD, FGF19, FGF3, FGF4, MIR548K, ORAOV1, and PPFIA1 (Figure 1F). Subsequently, the average CNV change value of each gene was calculated, indicating that most genes had a change value between 0 and 0.1 and the least number of genes had a change value greater than 0.4 (Figure 1G). The 10 genes with the highest number of mutations were TP53, TTN, MUC16, SYNE1, CSMD3, FLG, DNAH5, LRP1B, HMCN1, and PCLO (Figure 1H).

Enrichment analysis of antigen candidate genes

Genes with mutations equal to or greater than two are selected as candidate mutant genes. The intersection with overexpressed genes is then determined, resulting in a total of 3,970 mutant genes and 1,248 intersection genes (Figure S1A). GO and KEGG enrichment analyses were performed on the antigen candidate genes, and the main pathways of biological process enrichment were found to be regulation of DNA metabolic process, covalent chromatin modification, and histone modification. The main pathways for cellular component (enrichment were found to be spindle, nuclear chromatin, and nuclear envelope. The main pathways of molecular function enrichment were found to be cell adhesion molecule binding, ATPase activity, and protein serine. In further KEGG pathway analysis, the obtained pathways were mainly in human papillomavirus infection, PI3K-AKT signaling pathway, and focal adhesion (Figure S1B-S1E).

Survival analysis for identifying prognosis-related antigens

We further combined the OS data and DFI data, conducting a batch Cox univariate regression analysis on the antigen candidate genes. This analysis identified 24 genes significantly associated with OS and 55 genes significantly associated with DFI. By intersecting these gene sets, we identified six significant prognosis-related antigens (Figure 2A-2G), all of which exhibited a HR greater than 1, indicating that high expression levels of these genes are associated with reduced survival.

Figure 2 Survival analysis identified prognostic associated potential antigens. (A) Forest map of genes significantly associated with shared survival. Significantly related genes based on OS and PFS. (B-F) KM curves for potential prognosis-related antigens. (G) Venn diagram of genes significantly associated with OS and genes significantly associated with RFS. KM, Kaplan-Meier; Mut, mutation; OS, overall survival; PFS, progression-free survival; RFS, recurrence-free survival.

Correlation analysis of the potential prognosis-related antigens and immune cell infiltration

The correlation of prognosis-related potential antigens with APCs was analyzed using CIBERSORT pancancer data (22,23). Scatter plots showed significant correlations between PHF8, UBQLN2, MID1, and APCs such as B cells, suggesting that related genes may be involved in the antigen presentation process (Figure S2).

Identification of immune subtypes based on consistency clustering

Based on the expression matrix of six immune-related genes from TCGA data, unsupervised clustering of cancer samples was performed to obtain three immune subtypes (Figure 3A-3C). The expression heat map of immune-related genes showed a difference in the expression of related genes between the three subtypes (Figure 3D). PC analysis showed a clear distinction between the three subtypes (Figure 3E). Survival analysis indicated that the difference in OS between the three subtypes was not significant, but the difference in PFI was, with that in cluster A samples decreasing significantly faster than that of cluster B and C in the PFI survival curve. The same dimensionality reduction algorithm DDRTREE was used to cluster GSE19417, again yielding three subtypes, and survival analysis showed equally significant differences in survival between the three subtypes (Figure 3F-3H).

Figure 3 Immuno-subtyping of TCGA samples. (A) Heatmap of unsupervised clustering of samples. (B) Cumulative distribution graph. (C) Plot of consensus clustering cumulative distribution function. (D) Heatmap of expression of immune-associated antigen genes. (E) Survival analysis KM curves for each subtype based on PCA results for different subtypes. (F) Subtype analysis of TCGA-OS. (G) PFI. (H) OS survival analysis curve for GSE19417 data. CDF, cumulative distribution function; KM, Kaplan-Meier; OS, overall survival; PCA, principal component analysis; PFI, progression-free interval; TCGA, The Cancer Genome Atlas.

Comparison of subtypes in term of clinical features

We conducted separate counts of the differences among various subtypes of pathologic_stage, pathologic_N, and pathologic_T in TCGA samples. In the heatmap, darker colors indicate a higher −log(P) value of the correlation test. The results indicated that the immunological subtype composition of N3 samples differed significantly from that of N0 samples, with N3 samples having more cluster A subtype samples; there were significant differences between T4 and T1 and between T2 and T3 in terms of immunological subtype composition, but essentially all were composed of cluster A. Stage I samples were almost free of cluster A and differed significantly in immune-subtype composition from the other subtypes. Stage III and stage IV samples had a higher proportion of cluster A subtypes. The KM survival curves for cluster A curves decreased more sharply, and the analysis of clinical characteristics among immune subtypes indicated that the number of samples with cluster A increased as the stage and N grade increased, indicating that the survival of immune subtype subgroups was consistent with the severity of some clinical characteristics (Figure 4A-4C).

Figure 4 Clinical characteristics of the subtypes and presentation of CNVs. TCGA: (A) pathologic N, (B) pathologic stage, and (C) pathologic T. GEO: (D) tumor differentiation and (E) tumor volume (m3). (F-H) Immunoclusters A to C. The Y-axis is the G-Score, where red indicates amplification and blue indicates deletions. (I) Frequency statistics for the CNV locus of the clusters. CNV, copy number variation; GEO, Gene Expression Omnibus; TCGA, The Cancer Genome Atlas.

The GSE19417 sample was utilized to further investigate the differences in clinical characteristics among immune subtypes. In the heatmap, a darker color indicates a higher −log(P) value resulting from the correlation test (Figure 4D,4E). The results suggest that cluster B decreased proportionally with diminished differentiation and increased with increasing tumor volume, consistent with the significant results for OS.

Intersubtype mutations and the characterization of CNV

The CNV variation of the three immune subtypes was further characterized (Figure 4F-4I) through use of the mutated MAF files and combination with subtype groupings to map the differences in mutation landscapes of the samples between subtypes (Figure 5A-5C). In Figure 5D,5E, the violin plot indicates no significant difference in the degree of tumor mutation burden and number of mutant genes between the subtypes (24). The three subtypes differed in multiple immune cycle steps, and the final overall scores were also significantly different, with the cluster C sample having significantly higher overall immune activity than the cluster B sample (Figure 5F). The cluster C sample also had a higher number of homologous recombination defects (HRD). In addition, the cluster B sample had a higher mRNA expression-based stemness index (mRNAsi) (Figure 5G-5K).

Figure 5 Intersubtype mutation, immunological activity, HRD, neoantigen burden and chromosomal instability, and stemness index differences. (A-C) Mutation landscape of immunoclusters A to C. (D) TMB differences between subtypes. (E) Differences in the number of mutated genes between subtypes and mutations between subtypes. (F) Kruskal-Wallis test of intersubtype immunoreactivity. (G-K) Differences in aneuploidy score, HRD, SNV neoantigens, non-silent mutation rate, and mRNAsi between the subtypes. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant. HRD, homologous recombination deficiency; mRNAsi, mRNA-based stemness index; SNV, single nucleotide variant; TMB, tumor mutation burden.

Intersubtype differences in immune checkpoint and immunogenic cell death regulators

Immune checkpoint-related genes and regulators of immunogenic cell death were further collected, and then the differences in expression of the related genes between subtypes in TCGA and GEO data were noted. In TCGA database, consisting of 4 immune checkpoint genes and 26 immune prototype cell death regulator genes, the vast majority of immune checkpoints (ICPs) and immunogenic cell death modulators (ICDs) were significantly differentially expressed between the different immune subtypes. Similarly, in the GEO data, the majority of ICPs and ICDs were also significantly differentially expressed between the different immune subtypes (Figure S3).

Intersubtype immune cell analysis

Differences in immune cell proportions, stromal score, immune score, ESTIMATE score, and tumor purity were further calculated for all samples in TCGA dataset and the GEO dataset. The differences in correlation scores between subtypes were determined and visualized in box plots (Figure 6A-6L).

Figure 6 Intersubtype immune cell analysis of TCGA and GEO data. TCGA: (A) differences in immune cell ratios; (B) heatmap of immune cell ratios; and (C-F) stromal score, tumor purity, immune score, and cytolytic activity between the subtypes. GEO: (G) differences in immune cell ratios in GEO; (H) heatmap of immune cell ratios; and (I-L) stromal score, immune score, tumor purity, and cytolytic activity between the subtypes. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant. ESCA, esophageal cancer; GEO, Gene Expression Omnibus; NK, natural killer; TCGA, The Cancer Genome Atlas.

Analysis of subtypes and existing biomarkers for this tumor

The expression of HNRNPA2B1, CDK4/6, and FABP3 was strongly correlated with ESCA prognosis. The expression of HNRNPA2B1 and CDK4 was significantly lower in cluster C samples than in clusters A and B samples (Figure 7A-7D) in TCGA data and GEO data. The expression of FABP3 in cluster C samples was significantly higher than that in cluster A and cluster B (Figure 7E-7G).

Figure 7 Subtype and biomarker analysis for the tumor and immune landscape. Expression of other prognostic markers in (A-D) TCGA data and (E-G) GEO data. (H) Comparison of immune subtypes with existing pancancer immune subtypes. (I) Immune landscape of TCGA samples, with different colors representing different immunological subtypes. (J) Further breakdown of immune subtypes of the landscape. (K,L) Correlation of PC1 and PC2 with 22 immune cells and P values. The colors represent the magnitude of -log2(P), with darker colors representing smaller P values and greater significance. (M) Sample regroupings in the immune landscape are shown. (N) KM survival analysis of the resampled subgroups. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001; ns, not significant. GEO, Gene Expression Omnibus; KM, Kaplan-Meier; OS, overall survival; PC, principal component; TCGA, The Cancer Genome Atlas.

Comparison of immune subtypes with existing pancancer immune subtypes

In Figure 7H, the heatmap illustrates the P value [−log₁₀(P)] for the difference in detection, revealing a significant difference in the distribution of pancancer subtypes within the cluster C samples. Notably, the pancancer subtype C3 and C4 samples exhibit a greater tendency to be distributed within the immune subtype cluster C.

Characterization of the immune landscape

Trajectory analysis of cancer samples was conducted based on the count matrix of TCGA immune-related genes and was followed by construction of a tree diagram to visualize the distribution of immune subtypes in different trajectory samples. From the tree diagram in Figure 7I,7J, it can be seen that cluster C samples were further divided into different branches. The correlation between dendrogram PC1 and PC2 with the proportion of 22 immune cell types was further calculated. In Figure 7K,7L, the heat map shows that the abundance of gamma delta T cells, M0 macrophages, M1 macrophages, and resting dendritic cells had a significant positive correlation with PC1, while T regulatory cells (Tregs) and resting memory CD4 T cells had a significant positive correlation with PC1. Since there were cases in which an immune subtype was distributed to several different branches, PC1 and PC2 were used to subdivide the cluster, with the subdivided tree being displayed in Figure 7M. Subclassified samples of the three branches in the tree were obtained, survival analysis was performed on the subclassified samples, and KM curves were plotted, but the survival difference between subgroups of the subclassified samples was not significant (Figure 7N). Differences in the proportion of immune cells in subgroups within the cluster were then further examined via the Kruskal-Wallis test (Figure S4).

Construction of the weighted gene coexpression network and prognostic module

After soft threshold screening, the weighted gene coexpression network was constructed with a power of 10 and finally, two modules were obtained. Subsequently, the eigengenes of each module corresponding to each sample were calculated to detect whether there were differences in the distribution of eigengenes in different subtypes, and the results showed that the distribution of eigengenes varied in the different modules of the immune subtypes. Moreover, the distribution of eigengenes in the immunological subtypes was significantly different for each module. Further survival analysis of the modules using eigengenes showed that the gray module was significantly associated with prognosis. In addition, the eigengenes of the gray module were significantly correlated with PC1 of the immune landscape, suggesting that its internal genes are closely associated with both survival and immunity (Figure 8A-8G). The results of KEGG enrichment analysis for the gray module genes are significantly enriched in cytokine-cytokine receptor interaction, cell adhesion molecules and natural killer cell mediated cytotoxicity (Figure 8H).

Figure 8 Weighted gene coexpression network construction and prognostic module screening. (A) Sample classification tree diagram. (B) Soft threshold screening. (C) Network nodule clustering. (D) The number of genes in the module. (E) Differences in the distribution of module eigengenes between subtypes. (F) Single-gene regression analysis based on the module eigengenes. (G) Correlation analysis of the brown module with PC1. (H) KEGG functional analysis of genes within the gray module. (I) Risk score graphs for all samples. (J) Scatter plot of survival time for all samples. (K) Multifactor Cox regression for model construction and KM curve validation. (L) ROC curve validation. (M) Heatmap of gene expression in the model. ***, P<0.001; ****, P<0.001; ns, not significant. AUC, area under the curve; CI, confidence interval; FP, false positive; KM, Kaplan-Meier; KEGG, Kyoto Encyclopedia of Genes and Genomes; PC, principal component; ROC, receiver operating characteristic; TP, true positive.

Identification of module hub genes and construction of risk models

We divided the samples related to the module into high-risk and low-risk groups, and plotted the distribution of risk scores for the samples as well as the survival status of each sample (Figures 8I,8J). The correlations between genes within modules and module eigengenes were calculated, and 19 hub genes that were highly correlated with module eigengenes were obtained. Further multifactorial Cox regression analysis was performed on hub genes to obtain three genes for constructing the risk score model: BTN2A1 (0.070884346), HHLA2 (0.0059687), and MICA (–0.041615424); the risk scores of the samples were calculated using the above factors, and high- and low-risk groups were delineated using the mean as the node. The survival differences between high and low-risk groups were significant, and the 1-, 3-, and 5-year ROC curves were further plotted (Figure 8K-8M).

Discussion

Enhancements in surgical and chemotherapy techniques have significantly improved the treatment of ESCA over the years. However, challenges pertaining to efficacy and negative side effects still require attention. This is particularly crucial for patients with advanced ESCA, as their prognosis remains relatively unfavorable. Given the complexities associated with tumor heterogeneity and the intricate tumor immune microenvironment, the use of mRNA vaccines in clinical patients is not widespread. Instead, an alternative approach involves the implementation of unsupervised hierarchical cluster analysis, which is based on a comprehensive set of immune-related genes. This innovative strategy offers novel insights into selecting patients suitable for vaccination, substituting the development of a supervised learning model for patient risk stratification (25,26).

Tumors involve a complicated, diverse, and integrated ecosystem of relatively differentiated cancer cells, cancer stem cells, endothelial cells, tumor-associated fibroblasts, and infiltrating immune cells, among other cell types. The mRNAsi reflects gene expression, and there exists a robust correlation between the stemness index and new oncogenic pathways, somatic changes, microRNA, and regulatory networks for transcription. The progression of cancer involves complex immunological activities within the tumor microenvironment, where continuous and dynamic interactions occur between tumor cells and stromal cells. The tumor microenvironment includes ongoing and dynamic interactions between tumor cells and various types of stromal cells, such as fibroblasts, immune cells, and endothelial cells. Specifically, the tumor microenvironment is a highly complex and constantly changing ecosystem in which tumor cells interact closely and variably with multiple types of stromal cells. Therefore, a deep understanding of the complexity of immunological activities in the tumor microenvironment is crucial for developing more effective cancer therapies. Researchers are exploring how to leverage these complex interactions in order to identify new therapeutic targets or improve the effectiveness of existing treatments (27,28). As is well known, invasive immune cells in ESCA play a complex role in the tumor microenvironment, encompassing both anti-tumor immune cells (such as CD8⁺ T cells, NK cells) and immunosuppressive cells that facilitate tumor immune escape (such as Tregs, M2 macrophages, MDSCs). Through targeting these immune cells and their associated immunosuppressive mechanisms, more efficacious immunotherapy strategies can be developed to enhance the treatment outcome and survival rate of patients with ESCA (29).

The risk, progression, and treatment response of cancer are regulated by the pathway of DNA damage repair (DDR), with a loss of DDR function being a crucial factor in determining the risk, progression, and response to treatment of cancer. DDR genes can be classified based on functional pathways determined by genetic, biochemical, and mechanistic criteria. Proteins that belong to the same pathway often collaborate to repair specific types of DNA damage (30). Cancer treatment has been revolutionized by the emergence of cancer immunotherapy, which heavily relies on the development and activation of immune cells in the microenvironment of the host. The success of cancer immunotherapy depends on a series of stepwise events, collectively known as the cancer-immune cycle. These events include the release of cancer antigens, presentation of cancer antigens, initiation and activation, transfer of immune cells to the tumor, infiltration of immune cells into the tumor, recognition of cancer cells by T cells, and killing of cancer cells. The activity status of these seven-step, anti-cancer immune responses, along with the presence of immune cells infiltrating the tumor, contribute to the formation of a complex tumor immune phenotype within the tumor microenvironment. Thus, it is crucial to track tumor immunophenotypes to gain a comprehensive understanding of cancer immune mechanisms and to identify immunotherapeutic response biomarkers (31,32).

Unsupervised hierarchical clustering analysis relies on constructing a supervised learning model founded on an exhaustive collection of genes associated with the immune system. Its main purpose is to furnish novel perspectives regarding the identification of eligible patients for vaccination. The likelihood of a patient belonging to cluster A increases with the pathological grade of ESCA. Besides, the occurrence of cluster B diminishes as the degree of pathological differentiation diminishes. Notably, in our study, the immune activity exhibited by cluster C samples surpasses that of cluster B samples. Moreover, the immunotherapy strategies based on the molecular typing of ESCA can be integrated with the three stable immune subtypes (cluster A, B, and C) identified in this study. Firstly, the Cluster A subtype is likely to have relatively low immune cell infiltration and immune activity. The treatment strategy could consider the application of immune checkpoint inhibitors, such as programmed cell death protein 1 (PD-1) inhibitors (pembrolizumab, nivolumab, etc.), to enhance the recognition and killing ability of T cells against tumor cells. Meanwhile, chemotherapy or radiotherapy can be combined to augment the effect of immunotherapy. Whereas the Cluster C subtype is possibly characterized by higher immune cell infiltration and immune activity. Given its inherently high immune activity, immune checkpoint inhibitors may have better therapeutic effects. Additionally, the combination of adoptive cell immunotherapy (such as CAR-T cell therapy) can be contemplated to further strengthen the immune system’s attack on tumors (2,33).

In our study, our initial step involved the identification of potential antigens associated with ESCA. This was accomplished by identifying those mutated genes whose overexpression correlated with survival. Subsequently, we organized the immune gene sets into distinct clusters, creating novel immune subtypes. Additionally, we developed immune landscapes for TCGA samples based on these immune subtypes. Additionally, we performed prognostic analyses using WGCNA to identify correlates and build predictive models. The identified genes in the risk model can serve as valuable biomarkers for the prognosis of patients. To characterize the diversity within the ESCA samples, we implemented unsupervised clustering, resulting in the selection of three subtypes. Subsequent analyses exhibited significant distinctions between these three subtypes in various aspects. Importantly, our WGCNA-driven prognostic model accurately identified three potential antigens (BTN2A1, HHLA2, and MICA) that hold promise for mRNA vaccine development.

The regulation of T-cell function is attributed to HHLA2, a newly identified member of the B7 family (34). Several tumor types, such as clear-cell renal cell carcinoma, demonstrate the frequent expression of HHLA2. Interestingly, one study has shown that the expression patterns of HHLA2 and PDL1 in clear-cell renal cell carcinoma are not correlated, indicating that HHLA2 operates through a distinct mechanism unrelated to PDL1-mediated tumor immune evasion (35). Moreover, both the T-cell receptor (TCR) and CD28 signaling pathways are regulated by HHLA2, resulting in the inhibition of T-cell activation and proliferation (36).

The enhanced expression of major histocompatibility complex class I chain-associated proteins A and B (MICA/B) is attributable to cellular stress, and the release of MICA/B by neoplastic cells results in evasion of NKG2D recognition, thus facilitating oncogenesis (34). Additionally, exposure to low-dose chemotherapy induces senescence in neuroblastoma cells, and the acquisition of a senescent phenotype by tumor cells stimulates the secretion of the NKG2D ligand MICA/B via the MALAT1/miR-92a/ADAM10 pathway, thereby fostering the establishment of an inhibitory immune microenvironment and facilitating immune evasion (37).

Research suggests a correlation between BTN2A1 expression in cancer cells and the cytotoxicity induced by bisphosphonates in Vγ9Vδ2 T cells. The remarkable cytotoxic activity displayed by Vγ9Vδ2 T lymphocytes against tumor cells suggests their potential as candidates for anticancer treatment (38). The elimination of cancer cells by Vγ9Vδ2 T cells is governed by the action of anti-BTN2A1 monoclonal antibodies (mAbs), which relies on inhibiting the binding of BTN2A1 to Vγ9Vδ2 TCR. This supports the therapeutic prospects of BTN2A1 and enriches our understanding of the collaborative pathway of the butyric acid family in T-cell activation (39). Despite the relatively small sample size of the data, the survival results demonstrated robust statistical significance. In ESCA, specific molecules such as HER-2 and EGFR serve as potential tumor-associated antigens. These antigens are expressed on the surface of tumor cells and can be recognized by the immune system, thereby initiating an anti-tumor immune response (40). APCs, particularly dendritic cells (DCs), play a critical role in antigen processing and presentation to T cells, thus activating the immune response. However, in ESCA, tolerogenic dendritic cells (tDCs) are highly enriched, which suppress T cell activation and promote tumor immune evasion. MHC class II molecules, expressed on APCs, are essential for presenting antigens to CD4⁺ T cells. The expression levels of MHC class II in ESCA may be modulated by the tumor microenvironment. Therefore, the interplay between tumor-associated antigens, APCs, and MHC class II is crucial for understanding the immunological landscape of ESCA. Blocking immune checkpoints such as PD-1 can restore the function of CD8⁺ T cells, thereby enhancing the anti-tumor immune response (41,42). A comprehensive analysis of the interactions among these factors could provide a theoretical basis for developing more effective immunotherapy strategies, ultimately improving treatment outcomes and survival rates in patients with ESCA.

Conclusions

In this study, BTN2A1, MICA, and HHLA2 were found to be potent antigens that can be used in the development of mRNA vaccines targeting ESCA. Additionally, we identified three ESCA immune subtypes that demonstrate stability and reproducibility, which hold promise in predicting the outcome of mRNA vaccines. By shedding light on the role of mRNA vaccines, our study emphasizes the significance of BTN2A1, MICA, and HHLA2 as diagnostic factors and valuable therapeutic targets for ESCA. This study marks a pioneering effort to integrate the immune subtypes of ESCA with mRNA vaccine targets, providing a novel framework for personalized immunotherapy. Future research should explore the dynamic interaction networks within the tumor microenvironment across different subtypes to further refine and optimize combination treatment strategies.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the STREGA reporting checklist. Available at https://jtd.amegroups.com/article/view/10.21037/jtd-2025-233/rc

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

Funding: This study was supported by grants from the Nantong Science and Technology Project (No. JC22022068).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jtd.amegroups.com/article/view/10.21037/jtd-2025-233/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. Jiao R, Jiang W, Xu K, et al. Lipid metabolism analysis in esophageal cancer and associated drug discovery. J Pharm Anal 2024;14:1-15. [Crossref] [PubMed]
2. Yang H, Wang F, Hallemeier CL, et al. Oesophageal cancer. Lancet 2024;404:1991-2005. [Crossref] [PubMed]
3. Sayour EJ, Boczkowski D, Mitchell DA, et al. Cancer mRNA vaccines: clinical advances and future opportunities. Nat Rev Clin Oncol 2024;21:489-500. [Crossref] [PubMed]
4. Taina-González L, de la Fuente M. The Potential of Nanomedicine to Unlock the Limitless Applications of mRNA. Pharmaceutics 2022;14:460. [Crossref] [PubMed]
5. Jiang Y, Zhang H, Wang J, et al. Targeting extracellular matrix stiffness and mechanotransducers to improve cancer therapy. J Hematol Oncol 2022;15:34. [Crossref] [PubMed]
6. Zheng X, Xu H, Yi X, et al. Tumor-antigens and immune landscapes identification for prostate adenocarcinoma mRNA vaccine. Mol Cancer 2021;20:160. [Crossref] [PubMed]
7. Mai Z, Liu Q, Wang X, et al. Integration of Tumor Heterogeneity for Recurrence Prediction in Patients with Esophageal Squamous Cell Cancer. Cancers (Basel) 2021;13:6084. [Crossref] [PubMed]
8. Zhao F, Gao S, Qin X, et al. Comprehensive Analysis of TRP Channel-Related Genes for Estimating the Immune Microenvironment, Prognosis, and Therapeutic Effect in Patients With Esophageal Squamous Cell Carcinoma. Front Cell Dev Biol 2022;10:820870. [Crossref] [PubMed]
9. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-3. [Crossref] [PubMed]
10. Xu L, Deng C, Pang B, et al. TIP: A Web Server for Resolving Tumor Immunophenotype Profiling. Cancer Res 2018;78:6575-80. [Crossref] [PubMed]
11. Huang X, Tang T, Zhang G, et al. Identification of tumor antigens and immune subtypes of cholangiocarcinoma for mRNA vaccine development. Mol Cancer 2021;20:50. [Crossref] [PubMed]
12. Guo W, Tan F, Huai Q, et al. Comprehensive Analysis of PD-L1 Expression, Immune Infiltrates, and m6A RNA Methylation Regulators in Esophageal Squamous Cell Carcinoma. Front Immunol 2021;12:669750. [Crossref] [PubMed]
13. Guo H, Wang B, Xu K, et al. m(6)A Reader HNRNPA2B1 Promotes Esophageal Cancer Progression via Up-Regulation of ACLY and ACC1. Front Oncol 2020;10:553045. [Crossref] [PubMed]
14. Du Q, Guo X, Wang M, et al. The application and prospect of CDK4/6 inhibitors in malignant solid tumors. J Hematol Oncol 2020;13:41. [Crossref] [PubMed]
15. Liu T, Fang P, Han C, et al. Four transcription profile-based models identify novel prognostic signatures in oesophageal cancer. J Cell Mol Med 2020;24:711-21. [Crossref] [PubMed]
16. Huang TX, Fu L. The immune landscape of esophageal cancer. Cancer Commun (Lond) 2019;39:79. [Crossref] [PubMed]
17. Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 2014;32:381-6. [Crossref] [PubMed]
18. Jiang C, Liu Y, Wen S, et al. In silico development and clinical validation of novel 8 gene signature based on lipid metabolism related genes in colon adenocarcinoma. Pharmacol Res 2021;169:105644. [Crossref] [PubMed]
19. Peng XY, Wang Y, Hu H, et al. Identification of the molecular subgroups in coronary artery disease by gene expression profiles. J Cell Physiol 2019;234:16540-8. [Crossref] [PubMed]
20. Liu B, Gao Q, Liu B, et al. Application of Transcriptome Analysis to Understand the Adverse Effects of Hypotonic Stress on Different Development Stages in the Giant Freshwater Prawn Macrobrachium rosenbergii Post-Larvae. Antioxidants (Basel) 2022;11:440. [Crossref] [PubMed]
21. Huang Z, Wang S, Zhang HJ, et al. Characteristics of hypoxic tumor microenvironment in non-small cell lung cancer, involving molecular patterns and prognostic signature. Transl Lung Cancer Res 2021;10:2132-47. [Crossref] [PubMed]
22. Popi AF, Longo-Maugéri IM, Mariano M. An Overview of B-1 Cells as Antigen-Presenting Cells. Front Immunol 2016;7:138. [Crossref] [PubMed]
23. Jia Q, Wu W, Wang Y, et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun 2018;9:5361. [Crossref] [PubMed]
24. Rooney MS, Shukla SA, Wu CJ, et al. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 2015;160:48-61. [Crossref] [PubMed]
25. Liu L, Wu J, Shi M, et al. New Metabolic Alterations and A Predictive Marker Pipecolic Acid in Sera for Esophageal Squamous Cell Carcinoma. Genomics Proteomics Bioinformatics 2022;20:670-87. [Crossref] [PubMed]
26. Bregni G, Beck B. Toward Targeted Therapies in Oesophageal Cancers: An Overview. Cancers (Basel) 2022;14:1522. [Crossref] [PubMed]
27. Mai H, Xie H, Luo M, et al. Implications of Stemness Features in 1059 Hepatocellular Carcinoma Patients from Five Cohorts: Prognosis, Treatment Response, and Identification of Potential Compounds. Cancers (Basel) 2022;14:563. [Crossref] [PubMed]
28. Liu Z, He J, Han J, et al. m6A Regulators Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization In Nasopharyngeal Carcinoma. Front Immunol 2021;12:762243. [Crossref] [PubMed]
29. Schiffmann LM, Plum PS, Fuchs HF, et al. Tumor Microenvironment of Esophageal Cancer. Cancers (Basel) 2021;13:4678. [Crossref] [PubMed]
30. Oberdoerffer P, Miller KM. Histone H2A variants: Diversifying chromatin to ensure genome integrity. Semin Cell Dev Biol 2023;135:59-72. [Crossref] [PubMed]
31. Mills KHG. Innate lymphoid cells recruit T cells to turn up the heat on tumors. Cancer Cell 2022;40:362-4. [Crossref] [PubMed]
32. Kalaora S, Nagler A, Wargo JA, et al. Mechanisms of immune activation and regulation: lessons from melanoma. Nat Rev Cancer 2022;22:195-207. [Crossref] [PubMed]
33. Baba Y, Nomoto D, Okadome K, et al. Tumor immune microenvironment and immune checkpoint inhibitors in esophageal squamous cell carcinoma. Cancer Sci 2020;111:3132-41. [Crossref] [PubMed]
34. Oliviero B, Varchetta S, Mele D, et al. MICA/B-targeted antibody promotes NK cell-driven tumor immunity in patients with intrahepatic cholangiocarcinoma. Oncoimmunology 2022;11:2035919. [Crossref] [PubMed]
35. Bhatt RS, Berjis A, Konge JC, et al. KIR3DL3 Is an Inhibitory Receptor for HHLA2 that Mediates an Alternative Immunoinhibitory Pathway to PD1. Cancer Immunol Res 2021;9:156-69. [Crossref] [PubMed]
36. Rieder SA, Wang J, White N, et al. B7-H7 (HHLA2) inhibits T-cell activation and proliferation in the presence of TCR and CD28 signaling. Cell Mol Immunol 2021;18:1503-11. [Crossref] [PubMed]
37. Zhang Y, Hu R, Xi B, et al. Mechanisms of Senescence-Related NKG2D Ligands Release and Immune Escape Induced by Chemotherapy in Neuroblastoma Cells. Front Cell Dev Biol 2022;10:829404. [Crossref] [PubMed]
38. Laplagne C, Ligat L, Foote J, et al. Self-activation of Vγ9Vδ2 T cells by exogenous phosphoantigens involves TCR and butyrophilins. Cell Mol Immunol 2021;18:1861-70. [Crossref] [PubMed]
39. Cano CE, Pasero C, De Gassart A, et al. BTN2A1, an immune checkpoint targeting Vγ9Vδ2 T cell cytotoxicity against malignant cells. Cell Rep 2021;36:109359. [Crossref] [PubMed]
40. Gu YM, Zhuo Y, Chen LQ, et al. The Clinical Application of Neoantigens in Esophageal Cancer. Front Oncol 2021;11:703517. [Crossref] [PubMed]
41. Zhou H, Ma Y, Liu F, et al. Current advances in cancer vaccines targeting NY-ESO-1 for solid cancer treatment. Front Immunol 2023;14:1255799. [Crossref] [PubMed]
42. Jia Y, Zhang B, Zhang C, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Esophageal Squamous Cell Carcinoma. Adv Sci (Weinh) 2023;10:e2204565. [Crossref] [PubMed]

(English Language Editor: J. Gray)