Multi-cohort validation based on a novel prognostic signature of anoikis for predicting prognosis and immunotherapy response of esophageal squamous cell carcinoma

The general study process is detailed in Figure 1. From the Gene Expression Omnibus (GEO) database, we obtained the transcriptome data and clinical characteristics. The GSE53624↗ cohort encompassed 119 tumor samples and 119 normal adjacent samples. The tumor samples from the GSE53624↗ cohort was split into training (60 samples) and testing cohorts (59 samples) using the R package "caret" in a random manner. For external validation of prognostic signature, transcriptome data and associated clinical characteristics were extracted from GEO (GSE53625↗ cohort, n = 179) and TCGA (TCGA-ESCC cohort, n = 93) cohorts. In addition, we obtained IMvigor210 database through R package IMvigor210CoreBiologies (28) and selected samples with complete treatment response information (n = 298). A list of ARGs was obtained from previously published literature (48), and 779 ARGs were included for analysis after excluding genes absent in the GSE53624↗ cohort.

ESCC single-cell data were obtained from the GSE188900↗ dataset (sample size: n = 4). The Seurat v4.1 package was used to perform a standard pre-processing workflow for single-cell sequencing. Low-quality cells (mitochondrial genes > 20 %, gene numbers < 200, and gene numbers > 6000) were excluded from the subsequent analysis. Harmony (version 1.0) was employed to correct batch effects in the dataset comprising four samples and to integrate the merged objects. The t-distributed stochastic neighbor embedding (t-SNE) was utilized to visualize cell clusters.

In the GSE53624↗ cohort, we initially conducted a comparison of the expression levels of 779 ARGs between tumor and normal adjacent esophageal tissues. Subsequently, we conducted univariate cox analysis of differentially expressed ARGs to identify those associated with prognosis. Based on these prognostic ARGs, consensus clustering was conducted by employing the "ConsensusCluster Plus" R package with 1000 repetitions, 80% re-sampling, clusterAlg = "km", and distance = "euclidean". Then, a combination of the consensus score matrix, the CDF curve, and the PAC score was employed to determine the optimal number of clusters. The differentially expressed genes (DEGs) between subgroups were detected using the "limma" package, with criteria of |logFC| > 1 and FDR < 0.05. An ARGs signature was established using the GSE53624↗ training cohort and validated across multiple cohorts including the GSE53624↗ testing cohort, the GSE53624↗ entire cohort, the GSE53625↗ cohort, and TCGA-ESCC cohort. Prognostic DEGs were identified through univariate Cox regression in the GSE53624↗ training cohort followed by LASSO regression. The most effective prognostic signature was constructed using stepwise regression methods via multiple cox regression analyses. Risk score = ∑i=EXP (i) × Coef (i). Based on the medium risk score, each ESCC patient was classified as high- or low-risk accordingly. We utilized Kaplan-Meier analysis to assess the overall survival (OS) differences between subgroups in order to validate and evaluate performance. We used the ROC curve to assess their predictive ability (49). In addition, principal component analysis (PCA) was used to evaluate the grouping effect of ARGs signature (50). Taking into account various ESCC patients' clinical feature, Cox regression analyses were performed to investigate the potential of the ARGs signature as an independent risk factor.

Utilizing the "limma" package, DEGs between risk groups were identified (|logFC| > 0.585 and FDR < 0.05) (51). Relevant analyses were implemented through "clusterProfiler" R package (52) to examine the potential function of DEGs, including GO and KEGG analyses. Additionally, to evaluate potential pathways for signaling and biological functional alterations between groups, Gene Set Enrichment Analysis (GSEA) were applied by utilizing KEGG and Hallmark gene sets (p < 0.05).

In this study, we utilized the ESTIMATE method to calculate the stromal, estimate, and immune scores for each ESCC sample, as well as the tumor purity score. Additionally, we evaluated the abundance of 22 immune cells infiltration between different risk groups was assessed through CIBERSORT algorithm (53). Furthermore, we utilized single sample gene set enrichment analysis (ssGSEA) algorithms to measure the infiltration of 22 immune cells and overall immune function. Correlations between risk score and 22 immune cells infiltration were assessed using Pearson correlation coefficient.

Utilizing the "maftools" package, we conducted an analysis of TMB data obtained from TCGA database (54). Furthermore, TIDE scores for each patient were acquired through the online website (55). Additionally, we used the R package "oncoPredict" to predict the IC50 value of potential therapeutic drugs for ESCC between risk groups (p < 0.05) (56).

We obtained a tissue microarray containing 80 pairs of ESCC tissues and their corresponding adjacent normal esophageal tissues from Changsha Xiangya Biotechnology Co., Ltd. After dewaxing, antigen repair, and blocking endogenous peroxidase, the microarray was left to incubate overnight at 4°C with anti-F2RL2 primary antibody (bs-9510P, Bioss, China). Subsequently, secondary antibody incubation and visualization of immunoreactivity with DAB (G1211, Servicebio, China) were performed. In order to ensure an impartial evaluation, two pathologists who were unaware of the clinical information independently assessed the IHC results and resolved any discrepancies through discussion. The H-score method was used to evaluate staining intensity and extent by taking into account both intensity (ranging from 0 to 3+) and the percentage of tumor cells that were positively stained (ranging from 0% to 100%). The level of F2RL2 expression was equal to the product of these two estimates.

To conduct the statistical analyses, GraphPad and R 4.3.1 were used. According to the specific circumstances, we employed either the student t-test or Wilcoxon rank sum test for intergroup comparisons. Pearson correlation coefficient was utilized to assess correlation between variables. Chi-square test was implemented to evaluate whether there are differences in clinical characteristics. All statistical tests were considered significant at p < 0.05.

After conducting differential analysis between tumor and normal adjacent samples, 631 differentially expressed ARGs were identified (Supplementary Table S1). Through univariate Cox regression analyses, 66 prognostic ARGs were identified (Figure 2A). Subsequently, we conducted a consensus cluster analysis (k = 2-6) and found that the optimal number was obtained when k = 2 (Figure 2B). KM analysis showed significant differences in prognosis between two clusters (Figure 2C). Next, we conducted differential analysis between two clusters and identified 319 DEGs (Figure 2D; Supplementary Table S2). Subsequently, 119 ESCC patients were divided into two cohorts: training (n = 60) and testing (n = 59), maintaining an approximate 1:1 ratio. Table 1 displays the baseline characteristics of both the training and testing cohorts, showing no significant disparities. Through univariate Cox regression analysis, 36 prognostic DEGs were selected (p < 0.05) in the training cohort (Figure 2E). Subsequent LASSO regression analysis identified a minimum lambda value of 11 (Figures 2F, G). After multivariate Cox regression screening, an ARGs signature was developed based on four DEGs (Figure 2H) (RAMP1, F2RL2, FOXL1, and SLCO1B3). Following formula: Risk score = RAMP1 * 0.29303 + F2RL2 * 0.25188 + FOXL1 * 0.20157 + SLCO1B3 * (-0.14595), each ESCC patient's risk score was computed.

In accordance with the median score of training cohort, each ESCC patient was classified as either the high-risk or low-risk within their respective cohorts for GSE53624↗ training, GSE53624↗ testing, and GSE53624↗ entire. The heatmap of four modeling genes expression and the distribution of OS status and risk score were presented in Figures 3B–D. After conducting KM analysis, it was noted that the low-risk ESCC group exhibited better OS in comparison to the high-risk ESCC group across all cohorts (p < 0.05) (Figure 3A). Following internal validation, GSE53625↗ and TCAG-ESCC were selected as external cohorts to reconfirm the predictive effects of the signature. According to the results, both external validation and internal validation cohorts exhibit good cross-validation effects (Figures 3E, F). ROC analysis revealed that risk scores exhibited high levels of specificity and sensitivity in both the internal and external cohorts (Figures 3G–I).

The PCA results indicate that the ARGs signature has good grouping effect (Figures 4A–C). In addition, to explore the influence of the ARGs signature on patient prognosis in different clinical subgroups, we assessed the prognostic signature within various clinical features of ESCC. As shown in Figures 4D–I, in Age, Gender, and N subgroups, the low-risk ESCC group exhibited a higher survival rate than the high-risk ESCC group. Furthermore, to further validate the prognostic performance of the ARGs signature, we incorporated 32 published prognostic signatures and compared the C-index in the GSE53624↗, GSE53625↗, and TCGA-ESCC cohorts (Figures 4J–L). Our ARGs signature outperformed the majority of other published signatures in the GSE53624↗ and GSE53625↗ cohorts, and showed an intermediate performance in the TCGA-ESCC cohort.

We initiated univariate and multivariate Cox analyses on the GSE53624↗ and GSE53625↗ cohorts to analyze the prognostic importance of ARGs signature in relation to various clinical features. In GSE53624↗ cohort, as depicted in Figures 5A, N, stage, and risk score were identified as prognostic risk factors for ESCC patients through univariate regression Cox analysis. Subsequently, multivariate Cox regression analysis revealed risk score (p < 0.001) as an independent prognostic risk factor for ESCC patients. Furthermore, in the GSE53625↗ cohort (Figure 5B), we observed that risk score continued to be an independent prognostic risk factor, indicating its robust prognostic ability in ESCC patients.

We developed a nomogram utilizing risk score and various clinical features (Figure 5C). ROC analysis revealed that nomogram exhibited a high level of specificity and sensitivity (Figure 5D). The calibration curve showed high consistency between the findings of the nomogram and the observed probability of OS in practical application (Figure 5E). The results from the C index and DCA indicated that the nomogram has a more robust and strong predictive capability as well as net clinical benefit than other clinical features (Figures 5F, G), indicating that this nomogram has the potential to be utilized as a precise prognostic tool for ESCC patients.

Figure 6A and Supplementary Table S3 presents the DEGs between risk groups. KEGG pathway analysis indicated that DEGs were predominantly associated with “Cell adhesion molecules”, “Wnt signaling pathway”, and “Cytokine-cytokine receptor interaction” (Figure 6B). GO analysis, especially in the field of biological process (BP), indicated that there was a notable enrichment in terms of “immune system process”, “cell death”, “programmed cell death”, and “immune response” (Figure 6C). Furthermore, the GSEA analysis indicated a notable difference between subgroups (Figures 6D, E). According to “c2.cp.kegg_legacy.v2024.1.Hs.symbols.gmt” and “H.all.v2024.1.Hs.symbols.gmt”, we found that the high-risk group primarily showed activation of various cancer-related and immune-related signaling pathways, such as

“HALLMARK_WNT_BETA_CATENIN_SIGNALING“, "HALLMARK_KRAS_SIGNALING_UP",

"KEGG_ANTIGEN_PROCESSING_AND_PRESENTATION",“KEGG_TGF_BETA_SIGNALING_PATHWAY“, and “KEGG_INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION“, and the low-risk group mainly exhibits activation of the following signaling pathways, such as “KEGG LINOLEIC ACID METABOLISM“, “KEGG RETINOL METABOLISM“, and “HALLMARK_KRAS_SIGNALING_DN“.

According to the Wilcoxon test, high-risk ESCC patients showed significantly higher immune, estimate, and stromal scores, as well as lower tumor purity scores (Figures 7A–D). According to the findings in the cibersort (Figures 7E, F), macrophages M0, macrophages M2, mast cell resting, and T cell gamma delta were more abundant in high-risk ESCC patients, while in the low-risk group, plasma cells, monocytes, and mast cell activated were more abundant. By applying the ssGSEA algorithm, we found significant differences between the two risk groups (Figure 7G). Among them, macrophages M0, macrophages M2, and T cell gamma delta were more abundant in high-risk ESCC patients, which is consistent with previous results. Besides, we employed Pearson correlation analysis to identify 6 immune cell types that are significantly correlated with risk scores (p < 0.05, Figure 7H). Ultimately, we identified two intersecting tumor microenvironment cell types (macrophages M0 and T cells gamma delta, Figure 7I). Next, we obtained immune function scores using the ssGSEA algorithm. Compared with low-risk ESCC patients, high-risk ESCC samples exhibited greater enrichment in co stimulation of antigen-presenting cells (APCs), check point, cytolytic activity, HLA, MHC class I, and T cell co-inhibition (Figure 7J). Finally, we compared the immune checkpoints (ICs) between the high and low groups, and found that 24 ICs had significant differences between the two groups (Wilcox test, p < 0.05, Figure 7K).

In order to investigate the TMB between different risk groups, we conducted a mutation landscape within TCGA-ESCC patients (Figures 8A, B). Our investigation identified different mutation lineages in the high- and low-risk groups. For example, tumor suppressor genes NFE2L2 and NOTCH1 exhibited a higher frequency of mutations in the high-risk group at 22% and 18%, respectively, compared to 13% and 11% in the low-risk group. Moreover, an exploration into the association between risk scores and TMB revealed that individuals classified in the low-risk category exhibited elevated levels of TMB (Figure 8C, p = 0.041). When combining TMB with risk scores, it was observed that patients in the "high risk + high TMB" category experienced poorer outcomes (Figure 8D, p = 0.019).

To assess the potential of immunotherapy responses based on ARGs signature, we utilized the TIDE. Analysis of GSE53624↗, GSE53625↗, and TCGA-ESCC cohorts revealed higher TIDE scores in high-risk ESCC patients (p = 0.0031, Figure 8E; p < 0.001, Figure 8I; p = 0.0024, Figure 8L), with the mean risk score being significantly elevated in the non-response group compared to the response group (p = 0.035, Figure 8F; p < 0.001, Figure 8J; p = 0.02, Figure 8M). In the GSE53624↗ cohort, Pearson correlation analysis demonstrated a positive association between risk score and TIDE (Figure 8G, R = 0.23, p = 0.011), with a higher proportion of low-risk ESCC patients responding to immunotherapy (36%) than high-risk ESCC patients (25%) as shown in Figure 8H. Furthermore, compared to individuals in the high-risk group identified by risk score, analysis of both GSE53625↗ and TCGA-ESCC cohorts indicated a greater percentage of immunotherapy recipients in the low-risk group (Figures 8K, N). We further validated the predictive capability of ARGs signature for ICI responses using IMvigor210 cohort by calculating a risk score for each patient based on the coefficients and expression of the four modeling genes. Patients were then categorized into high-risk or low-risk groups according to their median risk score. Reassuringly, we observed that the average risk score was higher in patients belonging to the SD/PD group compared to those in the CR/PR group (Figure 8O, p = 0.016), and a lower proportion of high-risk patients achieved CR/PR compared to low-risk patients (Figure 8P, p = 0.011).

This involved correlating IC50 values with risk score of various drugs and comparing drug sensitivity scores between risk groups. Figures 9A–I displayed the nine compounds with the most significant correlation between IC50 values and risk score, as determined by Pearson correlation and Wilcox tests (p < 0.05).

To mitigate batch effects, we used the Harmony package to effectively integrate the four patients with ESCC (Figure 10A). Subsequently, the top 2000 variant genes underwent dimensionality reduction using principal component analysis and t-SNE. The cells were grouped into 26 clusters using a resolution of 1 during the clustering process. We categorized the cells into nine major clusters using marker genes for different cell types: myeloid cells, fibroblasts cells, T cells, endothelial cells, epithelial cells, B cells, and mast cells (Figure 10B). The heatmap shows the three most significant marker genes for every cell population (Figure 10D). Additionally, we calculated the proportions of cell clusters in each sample and presented the results as histograms (Figure 10C). Furthermore, we analyzed the expression patterns of the four modeling genes in various cell types (Figures 10E–I). The results indicated the RAMP1 and F2RL2 were predominantly expressed in fibroblasts cells, while SLCO1B3 and FOXL1 had lower expression levels in various cells.

The findings indicated that F2RL2 exhibited superior accuracy in predicting tumor status (tumor versus normal) compared to other variables (Figure 11D). Consequently, we proceeded with further validation of F2RL2 in ESCC. We collected 80 pairs of ESCC tissues and adjacent normal esophageal tissues for IHC staining. The protein level of F2RL2 in ESCC tissues was significantly higher than that in adjacent non-tumor tissues (Figures 11A–C). Paired t-test results indicated that the average expression of F2RL2 in ESCC tissues was higher compared to adjacent normal tissues (Figure 11E, p = 0.035). Based on the median expression of F2RL2, the 80 samples were stratified into high- and low-F2RL2 expression groups. KM survival analysis revealed that the high-F2RL2 expression groups exhibited a significantly poorer prognosis (p = 0.036, Figure 11F). These findings suggested that F2RL2 could serve as a valuable prognostic biomarker related to anoikis in ESCC.

The original contributions presented in the study are included in the article/. Further inquiries can be directed to the corresponding author. 1