Establishment and validation of a 3-gene prognostic model
Using the TCGA and the GTEx resources, we first performed univariate and multivariate Cox analyses on the differentially expressed immune-related genes (IRGs) across the esophageal cancer patient samples to obtain prognostic-related genes (Figure S1 and Table S2), which were further verified by performing qRT-PCR. Using LASSO Cox regression in the Shihezi cohort, we identified a 3-gene risk score model based on the genes CCL25, IGLV10-54 and OSM (risk score = 0.0778×normalized expression level of CCL25-19.3202×normalized expression level of IGLV10-54+16.8585 × normalized expression level of OSM). We used the median level of the risk score to classify EC patients into low and high-risk groups. In the TCGA dataset, we found that the survival rates between low and high-risk patients were significantly different, with patients in the high-risk group having a significantly worse survival rate than their low-risk counterparts (p<0.01, Figure 1a). In addition to predicting survival, our gene model was markedly robust in the diagnostic accuracy as shown by the excellent receiver operating characteristic (AUC) values (risk score AUC = 0.735; Figure 1b). Separate analyses of esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC) yielded similar results (Figure S2).
In order to further test the performance of the 3-gene prognostic model, these findings were subsequently validated in the Shihezi datasets. Similar to the analysis of the TCGA cohort, high-risk patients had worse OS rates than those with low risk (p<0.001; Figure 1c). Likewise, in the Shihezi cohort, the diagnostic accuracy of the model was remarkable (AUC=0.851) and much higher than other clinicopathological characteristics (Figure 1d). Moreover, Kaplan-Meier curves and ROC analysis showed that the performance of the 3-gene model was better than the performance of any individual gene (Figure S3).
Next, the predictive value of the 3-gene model was compared to that of EC-related prognostic biomarkers published in the literature. The 3- gene model reached an AUC of 0.714 at 1 years, 0.635 at 3 year, and 0.766 at 5 years in the TCGA cohort, while AUC values for 1-, 3-, 5-year survival were 0.644, 0.696, 0.588 in a recent study by Wang et al. 21, suggesting that the 3-gene model is better at predicting 1 and 5-years OS than other biomarkers (Figure 2).Taken together, these results demonstrate the robustness of our genetic model in accurate diagnosis and prognostic prediction of EC.
The gene model serves as an independent predictor in EC patients
Next, we carried out univariate and multivariate Cox regression analyses by combining a variety of clinicopathological features in the TCGA and Shihezi datasets. In the univariate analysis, advanced stage (stage 1/2 vs. stage 3/4;TCGA p<0.001, Shihezi p<0.05), positive lymph node metastasis (N;TCGA p<0.001, Shihezi p<0.05 ), and the 3-gene signature (TCGA p=0.025,Shihezi p=0.002) were found to be reliable predictors of poor survival in both the TCGA and the Shihezi cohorts (Figure 3a and 3b). After correction for other confounding factors, subsequent multivariate analysis revealed that both the 3-gene model (TCGA p=0.009, Shihezi p=0.020) and positive lymph node metastasis(N)(TCGA p=0.026, Shihezi p=0.004) were independent predictors of poor survival in EC (Figure 3c and 3d). However, there was no obvious relationship between our model and clinical parameters (Figure S3). These results confirm that our genetic model can accurately predict overall EC survival and that its predictive value is comparable to or superior to the clinicopathological features currently used in clinic. These results also show that the model can be used as an independent predictor of survival in patients with EC.
Establishment of a nomogram for predicting survival in EC
To better predict the survival probability in EC patients, we developed a nomogram based on the Shihezi cohort to evaluate the probability of the 12- , 24‐, and 36-month OS by combining multiple significant univariate features. The nomogram contained 3 prognostic factors, including clinical stage, N and risk score (Figure 4a). In order to verify the model, the prediction performance of the nomogram was evaluated through the calibration chart of the model and calculation of the discrimination index. The prediction value of the nomogram was higher than that of either the gene model, N or clinical stage alone as the C-index of the nomogram reached 0.788, while it reached 0.741 for the gene-model, 0.609 for N and 0.564 for clinical stage. Calibration plots for the prediction of 12 and 24- month survival rates showed a good performance of the nomogram (Figure 4b and 4c). These findings suggest that the clinical significance of our genetic model can be improved by combining it with other factors such as N and clinical stage.
Enriched molecular pathways in high and low-risk EC subgroups
To elucidate the molecular mechanisms that might mediate the good performance of the 3-gene model, the TCGA discovery cohort was analyzed by GSEA. As shown in Figures 5a-d, BLADDER CANCER and INTESTINAL IMMUNE NETWORK FOR IGA PRODUCTION pathways were significantly enriched in the tumors of the high-risk group as defined by the 3-gene signature. In the low-risk group, the significantly enriched pathways were mainly related to PRIMARY IMMUNODEFICIENCY, T CELL RECEPTOR SIGNALING PATHWAY, THYROID CANCER, and VIRAL MYOCARDITIS. Analysis of the individual three genes in our model showed that the samples were mainly enriched in pathways such as various metabolic processes, immune-related pathways and tumor classic pathways(Figures 5a-d). In addition, we investigated the tumor mutational burden (TMB) of high and low-risk patients and found that high-risk patients had higher TMB than low-risk patients (p<0.05; Figure 5e).
Correlation between the genes model and infiltrating immune cells
To explore the tumor microenvironment (TME) of patients with EC, we used the ESTIMATE algorithm to estimate the overall immune status of the high-risk and low-risk groups. The results showed that the tumor purity of the high-risk group was higher than that of the low-risk. In contrast, the values of the ESTIMATE score and immune score were lower in the high-risk cluster than in the low-risk group (Figure 6a). To better understand the cell type composition of esophageal tumors, the deconvolution tool CIBERSORT was used to estimate the proportions of 22 tumor-infiltrating immune cells (TIIC). Similar to other tumors, naive B cells, CD8 T cells, CD4 memory resting T cells and Macrophages (M0, M1, M2) were of the most abundant infiltrating immune cells in esophageal carcinoma (Figure 6b). Moreover, naive B cells, CD4 memory resting T cells, CD4 memory activated T cells, regulatory T cells (Tregs), M0 Macrophages, M1 Macrophages, resting Dendritic cells, and activated Dendritic cells accounted for larger proportions in EC tumors than the normal tissue (Figure 6c). In order to determine whether risk scores based on our 3-gene model are related to the patient's TME status, we assessed the abundance of distinct immune cells in the low and high-risk groups. We found that naive B cells, CD8 T cells and plasma cells were predominant in the low-risk group, while M0 Macrophages, CD4 memory resting T cells and activated Mast cells were predominant in the high-risk group, the latter showing a distribution trend of partially infiltrating immune cells that was similar to that of tumor tissue (p<0.05; Figure 6d). These findings indicate that these subsets of immune cells may play an important role in the predictive ability of our model.
Choice of treatment methods for different risk groups
Since the high and low-risk groups have different immune scores and consist of different immune cell populations, and since immune checkpoints are key targets of immunotherapy, we investigated important immune checkpoint genes in the low and high-risk immune groups. We found that the expression levels of these checkpoint genes (including those encoding CTLA-4 and PD-1/PDL-1) were significantly increased in the low-risk group (p< 0.05; Figure 7a-c). Using the TIDE program to investigate the effectiveness of immune checkpoint (PD-1 and CTLA4) inhibitors in the two groups, we found that the high-risk patients had markedly lower TIDE scores compared to low-risk patients, indicating that low-risk patients may have a better outcome when receiving immune checkpoint (PD-1 and CTLA-4) inhibitors (p<0.05; Figure 7d).
To identify other potentially effective immune drugs for patients in the high-risk group, we used the pRRophetic algorithm to predict the IC50 of 70 agents in the high and low-risk patients. We found that several drugs (Bryostatin.1, FH535, PHA.665752, RDEA119 and X17. AAG ) had lower IC50 in low-risk patients relative to high-risk patients (p<0.05; Figure 7e-i and Table S3), and that the IC50 of Bryostatin.1 was the lowest among these drugs, indicating that the high-risk patients were more sensitive to these drugs, which may be used to treat patients at high risk.