Multi−cohort validation based on coagulation-related genes for predicting prognosis of esophageal squamous cell carcinoma

Figure 1 exemplified the investigation process. In this study, four transcriptome cohorts encompassing complete overall survival rate (OS) information and clinical information (TCGA-ESCC, GSE53625↗, GSE53624↗, and GSE53622↗) were obtained from TCGA and GEO databases. Besides, two other transcriptome cohorts from the GEO database were included (GSE20347and GSE38129↗). All the six transcriptome cohorts undergo a log-2 transformed to ensure normalization. Subsequently, to eliminate the batch effect, the ComBat algorithm was employed, and the extent of correction was examined via principal component analysis (PCA, Supplementary Figures 1A, B).

The pre-preprocessing of ESCC single-cell data (with a sample size of n = 7 for ESCC)) from the GSE145370↗ dataset was conducted using the Seurat v4 and Harmony (version 1.0). Low-quality cells with a mitochondrial gene proportion >20% and a gene count of <200 and >6000 were removed. To visualize the cell clusters, the t-distributed random neighbor embedding (t-SNE) was employed.

Furthermore, a curated list of 203 CRGs (), such as SERPINE1, F3, and THBD, were obtained from the Gene Set Enrichment Analysis (GSEA) database. These genes originated from the hsa04610 pathway (complement and coagulation cascade) and the hsa04611 pathway (platelet activation). Both of these pathways are closely linked to coagulation and relevant processes. 1

Initially, a univariate cox analysis of 203 CRGs in the GSE53625↗ cohort was conducted to identify prognostic CRGs. Next, the ‘ConsensusCluster Plus’ R package was utilized for consensus clustering with these prognostic CRGs. The optimal cluster count was determined by making use of the CDF curve, consensus score matrix, and PAC score. To detect prognostic DEGs between clusters, the ‘limma’ package (with a |logFC| > 0.585 and a p-value < 0.05) and univariate cox regression were utilized (28). In this research, the samples from the GSE53625↗ cohort were randomly assigned to a training cohort with 90 samples and a testing cohort encompassing 89 samples through the utilization of the R package ‘caret’. The clinical information of the two cohorts was shown in Table 1. Here, the CRGs signature was formulated with GSE53625↗ training cohort and verified in diverse cohorts, including the GSE53625↗ testing, the entire GSE53625↗, the GSE53624↗, the GSE53622↗, and the TCGA-ESCC cohorts. In the GSE53625↗ training cohort, hub prognostic DEGs were screened out through univariate cox regression, LASSO regression, and multiple cox regression analysis with a stepwise approach. Each sample’s risk score was evaluated utilizing the formula: Risk score =∑i=EXP (i) ×Coef (i). Next, the median score of each cohort was tallied, and each ESCC sample was categorized as high- or low-risk based on the median score. We utilized the ‘Survival’ R package and employed the ‘survival ROC’ package to evaluate the predictive capability (29). Additionally, cox regression and nomogram analysis were performed in conjunction with the clinical features.

Researchers executed relevant enrichment analysis on DEGs utilizing ‘clusterProfiler’ R package (30), encompassing GO and KEGG analyses. Moreover, to analyze potential modifications in signaling pathways, Hallmark gene sets were applied in GSEA.

To analyze the differences between risk groups, the immune, estimate, stromal, and tumor purity scores were calculated by utilizing the ESTIMATE method. Moreover, ssGSEA algorithms were utilized to assess both the overall immune function and the immune cell infiltration.

Furthermore, the software package “maftools” was utilized to obtain tumor mutation burden (TMB) data (31). Meanwhile, the microsatellite instability (MSI) score was retrieved from the public data of TCGA (Supplementary Table 2). Besides, each patient’s TIDE scores were from the online website (32). Moreover, “oncoPredict” package was utilized to perform drug sensitivity analysis (p < 0.05) (33).

Three ESCC cell lines, KYSE30, KYSE150, and KYSE410, were sourced from Pricella (Wuhan, China), and cultured in RPMI-1640 (Pricella) supplemented with 10% FBS (Pricella) and 1% penicillin-streptomycin (Pricella). These cells were cultured at 37 °C in 5% CO2. The siRNAs for SULT1B1 (si- SULT1B1-1/-2/-3/-4, details in), negative control (si-NC), and SULT1B1 overexpression plasmid were provided by GenePharma (Shanghai, China). And to transfect the cells, Lipofectamine 3000 regent (L3000015, Invitrogen, USA) was employed. Forty-eight hours post-transfection, cells were harvested to analyze the proteins levels, apoptosis, and cell cycle. Annexin V-FITC Apoptosis Detection Kit (Beyotime, C1062, China) and Cell cycle Analysis Kit (Beyotime, C1052, China) were used for apoptosis and cell cycle assay. 1

To prepare protein extracts, cells were lysed with a mixture of RIPA buffer (Beyotime, P0013B, China) and 1% PMSF (Beyotime, ST505, China). The protein samples were then denatured, and the proteins were resolved by SDS-PAGE gel. Subsequently, the proteins were transferred onto a PVDF membrane, which was blocked with 5% skim milk for 2 hours. Next, the membrane was incubated with primary antibodies, namely SULT1B1 (1:500; Proteintech, 16050-1-AP, China), E-cadherin (1:5000; Proteintech, 20874-1-AP, China), Vimentin (1:1000; Beyotime, AF1975, China), and GAPDH (1:1000; Servicebio, GB15004-100, China). After incubation with the secondary antibody, the bands were developed by means of the ECL developer (Beyotime, E422-01, China).

A density of 2000 cells per well was utilized for seeding in 96-well plates. After that, to each well, 10 µL of CCK-8 solution in a serum-free medium was added and incubated at 37 °C for one hour. Next, measurements of absorbance were conducted at a wavelength of 450 nm, with the optical density (OD) being recorded at the same daily interval.

A volume of 100 µL of serum-free medium, which contained 5 × 10⁴ cells, was seeded on the top chamber of the transwell. Concurrently, 600 µL of complete medium was inoculated into the bottom chamber. Following 24 hours of incubation, the cells were stained with a 0.5% crystal violet solution and then counted under a microscope.

Utilize the sterile plastic pipette to scrape the transfected cells. Subsequently, wash the cells twice with PBS, and then replenish the culture with fresh medium. Images of the scratch were captured under a microscope at 0 and 24 hours post-treatment.

For statistical analysis and graphing, R 4.3.1 and GraphPad Prism 8 were utilized. Differences between two groups were analyzed utilizing Student’s t-test or Wilcoxon’s rank sum test. For statistical analysis above two groups, a one-way ANOVA test was performed. A p < 0.05 was considered statistically significant.

Univariate cox regression was performed on 203 CRGs, and 19 prognostic CRGs were filtered out (Figure 2A). We then proceeded with a consensus cluster analysis, and identified the optimal number k was 2 (Figures 2B, C). The differences between the two consensus clusters (C1 and C2) in terms of 19 prognostic CRGs and clinical characteristics were illustrated in Figure 2D. According to the KM analysis, there were marked prognostic differences between the two clusters (Figure 2E). Next, 1025 DEGs (Figure 2F; Supplementary Table 4) and 215 prognostic DEGs were filtered out. No notable variations were observed between the two cohorts. Figure 2G illustrated the 36 prognostic DEGs in the training cohort. Subsequently, six model DEGs were screened out through the utilization of LASSO regression (Figures 2H, I) and multivariate Cox regression. Supplementary Table 5 presented the univariate and multivariate results of these six model DEGs. Thereafter, a CRGs signature was established (Figure 2J) according to the formula: Risk score = PTX3* 0.18815 + CILP * 0.14112 + CFHR4* (-0.17575) + SULT1B1* (-0.19359) + IL5RA * (-0.29789) + FAM151A * (-0.49759).

The OS status and risk score distribution for GSE53625↗ training, testing, and entire cohorts were depicted in Figures 3D–I. Besides, within the aforementioned three cohorts, notably poorer prognosis was detected in high-risk group (Figures 3A–C). Additionally, the validity of signature prediction was corroborated by the TCGA-ESCC, GSE53624↗, and GSE53622↗ cohorts (Figures 3M–O). The time-dependent ROC of the signature was plotted (Figures 3J–L, P–R). The results indicated that across all cohorts, the AUC values at 1–5 years all exceeded 0.6. This observation evidenced high specificity and sensitivity. Furthermore, a random-effects meta-analysis was conducted on the hazard ratios (HR) across the above four cohorts (GSE53625↗, TCGA-ESCC, GSE53624↗, and GSE53622↗), and the results showed that the CRGs signature was a risk factor for OS in ESCC (HR = 1.73, 95% CI = 1.52 - 1.97, I² = 0, illustrated in Supplementary Figure 1C).

Figures 4A–C illustrated that the CRGs signature shows a clear grouping effect as revealed by PCA. Moreover, we investigated the influence of clinical features on the CRGs signature, as depicted in Figures 4D–K, within subgroups including age, T, N and stage, the prognosis of the low-risk group was better. Furthermore, we conducted a comparison of the CRGs signature with 32 previously published prognostic signatures. The findings indicated that in the GSE53625↗, TCGA-ESCC, GSE53624↗, and GSE53622↗ cohorts, the CRGs signature was more effective in comparison to other published signatures (Figures 4L–O). In conclusion, these above results emphasized the remarkable predictive accuracy of CRGs signature in forecasting the prognosis of ESCC patients.

We utilized univariate and multivariate cox analyses in the GSE53625↗ (Figures 5A, B), GSE53624↗ (Figures 5C, D), and GSE53622↗ (Figures 5E, F) cohorts to explore the prognostic implications of CRGs signature alongside various clinical features. In the cohorts mentioned above, the CRGs signature was demonstrated to be an independent prognostic risk factor (p < 0.05), emphasizing its significant prognostic potential. Next, a nomogram was developed that integrates risk scores with several clinical features, illustrated in Figure 5G. Time-dependent ROC analysis indicated that this nomogram possessed high sensitivity and specificity (Figure 5H). The calibration curve verified the feasibility of this nomogram in practical settings (Figure 5I). The outcomes of DCA indicated and the C index indicated that, in comparison to other clinical features, this nomogram possessed stronger predictive power and a higher net clinical benefit (Figures 5J, K). Overall, this nomogram holds the potential to emerge as an effective instrument for the accurate prognosis of patients with ESCC.

The DEGs (illustrated in Figure 6A and Supplementary Table 6) were mainly concentrated in “Inflammatory mediator regulation of TRP channels” and “Complement and coagulation cascades” (as shown in Figure 6B). GO analysis reveals a notable enrichment within the domain of biological process (BP), specifically in relation to “cellular developmental process” and “cell differentiation” (Figure 6C). Additionally, the GSEA analysis revealed that the high-risk groups predominantly display the activation of multiple cancer-associated signaling pathways, such as “UV RESPONSE DN” and “EPITHELIAL MESENCHYMAL TRANSITION”. Alternatively, the low-risk group was predominantly marked by “KRAS SIGNALING DN” and “P53 PATHWAY” (Figure 6D). The correlation analysis between CRGs scores and hallmarks pathway scores further indicated that CRGs scores was strongly associated with cancer-related biological processes and metabolic pathways (Figure 6E). To sum up, these results suggested that the activation or inhibition of these pathways might give rise to distinct prognostic outcomes observed in different CRGs signature subgroups.

Patients with high-risk ESCC exhibited notably elevated immune, stromal, and estimate scores, in conjunction with decreased tumor purity score (Figures 7A, C, E, G). Besides, the risk score was positively correlated with immune, estimate, and stromal scores (Figures 7B, D, F). In contrast, a negative correlation was observed in tumor purity (Figure 7H). By utilizing the ssGSEA algorithm and the Wilcoxon test, notable differences were found (Figure 7I). Specifically, activated dendritic cells, fibroblasts, and macrophages M0 were notably more common in patients with high-risk ESCC. In contrast, in low-risk ESCC patients, there was a higher abundance of activated mast cells, naive CD4 T cells, and plasma cells. Besides, Pearson correlation analysis pinpointed 8 immune cells (p < 0.05, Figure 7J). Moreover, five intersecting immune cell types were identified (Figure 7K). Then, researchers evaluated the associations among the 6 CRGs and immune cells (Figure 7L). Thereafter, the two risk groups exhibited disparities in HLA enrichment (Figure 7M). Finally, an analysis of immune checkpoints (ICs) was performed, and seventeen ICs exhibited notable (p < 0.05, Figure 7N). These findings underscored the distinctions in immune cell infiltration between CRGs signature subgroups.

Initially, in the TCGA-ESCC cohort, the waterfall plot was employed to depict the somatic mutation landscape (Figure 8A). Subsequent analysis unveiled patients within low-risk subgroup possessed higher TMB levels (Figure 8B), and the risk score was inversely related to TMB (R = -0.25, p = 0.025, Figure 8C). Moreover, the results suggested that patients in the “high risk + high TMB” category exhibited the poorest prognosis (Figure 8D, p = 0.031). Furthermore, we assessed the MSI score between different risk groups. No significant disparity in MSI score between different risk groups was detected (Supplementary Figure 1D), and no significant correlation existed between the risk score and the MSI score (Supplementary Figure 1E).

We employed TIDE to assess the potential of the immunotherapy response between different risk groups. Analyses of the GSE53625↗, GSE53624↗, GSE53622↗, and TCGA-ESCC cohorts revealed that, when compared with low-risk ESCC patients, the TIDE score in high-risk ESCC patients was relatively higher. (Figures 8E, I, M, Q). Moreover, as opposed to the responder group, the average risk score of the non-responder group exhibited an upward tendency (Figures 8G, J, N, R). Additionally, a positive correlation was identified (Figures 8F, K, O, S). Furthermore, in the aforementioned four cohorts, a higher proportion of patients in the low-risk group were predicted to respond effectively to immunotherapy (Figures 8H, L, P, T). These findings indicated that the CRGs signature hold substantial potential in predicting the response to immunotherapy.

Through Pearson correlation analysis, nine drugs with the strongest correlation between IC50 values and risk scores were screened out. Patients with high-risk ESCC showed higher sensitivity to Erlotinib, Acetalax, Gefitinib, Afatinib, Ibrutinib, Sapitinib, AZD3759, Lapatinib, and SCH772984 (Figures 9A–I). The results of the Wilcoxon test revealed that notable statistical differences were identified (p < 0.05). The findings may provide insights into the selection of potential treatment options to inhibit the malignant progression of cancer.

Figure 10A depicted the integration results of 7 ESCC patients after eliminating the batch effect. Subsequently, the cells were classified into nine clusters (Figure 10B), and Figure 10D illustrated the three marker genes in each cell clusters. Moreover, the histogram presented the proportions of cell clusters in the 7 samples (Figure 10C). Finally, the expression patterns of five model CRGs were analyzed (Figures 10E–J), and FAM151A was not detected in this dataset. The research findings showed that the expressions of five model CRGs exhibited significant differences. Specifically, SULT1B1 was predominantly expressed in fibroblasts.

In the GSE53625↗ cohort, consisting of 179 tumor samples and 179 normal adjacent samples, the capacity of six model CRGs to predict tumor states was assessed utilizing ‘pROC’ R package. The ROC results indicated that among the six model CRGs, SULT1B1 exhibited the highest accuracy in predicting the tumor status (Figure 11A). Consequently, we selected it as the focal point and investigated its role in ESCC. Initially, by accessing various online databases, such as GEPIA (http://gepia.cancer-pku.cn/↗), TNMplot (https://tnmplot.com↗), and Kaplan-Meier plotter (http://kmplot.com↗), we analyzed the expression level and prognosis of SULT1B1 in ESCC patients. The analysis revealed that, in comparison with normal tissues, a notable down-regulation of SULT1B1 was detected in ESCC cancer tissues (Figures 11B, D). Additionally, when compared with the SULT1B1-lower group, higher SULT1B1 patients expression presented better OS (Figure 11C). Furthermore, we conducted a comparison of the expression of SULT1B1 across five transcriptome cohorts retrieved from the GEO databases. The findings indicated that, in comparison with normal tissues, SULT1B1 was significantly down-regulated in ESCC cancer tissues (Figures 11E–I). Moreover, in the four prognostic cohorts, it was noted that, in compared to the SULT1B1-lower group, ESCC patients with higher SULT1B1 expression demonstrated better OS (Figures 11J–M). Subsequently, we examined three pairs of clinical ESCC tissue samples. The results showed that compared to the paired normal tissues, ESCC cancer tissues demonstrated decreased SULT1B1 expression (Figure 12A). The collected data proposed that SULT1B1 might have a tumor-suppressing function in the progression of ESCC.

To investigate the role of SULT1B1, we initially assessed the basal protein expression levels of SULT1B1 in KYSE150, KYSE30, and KYSE410 (Figure 12B). The results indicated that SULT1B1 exhibited high expression in KYSE150 cell line, moderate expression in KYSE410 cell line, and low expression in the KYSE30 cell line. Subsequently, SULT1B1 was knocked down in the KYSE150 and KYSE410 cell lines (Figure 12C), while it was overexpressed in the KYSE30 and KYSE410 cell lines (Figure 12D). The efficiency of overexpression and knockdown was substantiated by western Blot analysis. Here, the subsequent experiments utilized si-SULT1B1#4, as it exhibited a higher level of knockdown efficacy. The CCK8 assay results demonstrated that the knockdown of SULT1B1 was capable of promoting the tumor cell proliferation of KYSE150 and KYSE410 cells (Figures 12E, H). Conversely, the overexpression of SULT1B1 was able to suppress the proliferation of KYSE30 and KYSE410 cells (Figures 12N, Q). Flow cytometry analysis revealed that upon knockdown of SULT1B1 in KYSE150 and KYSE410 cells, the proportion of cells in the G0/G1 phase showed a decline, whereas the proportion of cells in the S phase exhibited an increase (Figures 12F, I). Additionally, the apoptosis rate decreased (Figures 12G, J). Conversely, after the overexpression of SULT1B1 in KYSE30 and KYSE410 cells, the proportion of cells in the G0/G1 phase showed an upward trend, whereas the proportion of cells in the S phase declined (Figures 12O, R). Moreover, the apoptosis rate increased (Figures 12P, S). Results from the scratch assay and transwell invasion experiments confirmed that knockdown of SULT1B1 promoted the migration of KYSE150 and KYSE410 cells (Figures 12K, L), while overexpression of SULT1B1 inhibited the migration of KYSE30 and KYSE410 cells (Figures 12T, U). The findings of western Blot analysis indicated that upon the knockdown of SULT1B1, the expression of E-cadherin significantly decreased, while that of Vimentin significantly increased (Figure 12M). Conversely, when SULT1B1 was overexpressed, an opposite outcome was observed (Figure 12V). In summary, the findings of our research indicated that SULT1B1 is capable of effectively suppressing the proliferation and migration of ESCC cells. Further mechanistic investigations confirmed that its tumor-suppressing function is achieved through promoting cell cycle arrest at the G0/G1 phase, inducing apoptosis, and suppressing the epithelial-mesenchymal transition (EMT) process.