An immunoscore nomogram for predicting survival in patients with oesophageal cancer

Clinical characteristics of study cohorts

An overview of this study is delineated in Figure 1. The clinical information of the two oesophageal cancer cohorts was downloaded from the public Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases. Based on the inclusion criterion (CIBERSORT P < 0.05), a total of 374 samples were included in this study, including 104 normal samples from the GEO database and 270 cancer samples from GEO (147) and TCGA (123). The cancer samples from the GEO database were the training cohort to identify prognostic markers and construct the model, which was then validated with the TCGA data. In the training cohort, the median age at diagnosis was 60.0 (54.0-66.0) years, and 79.6% of the patients were men. In the validation cohort, the median age at diagnosis was 60.0 (54.0-73.0) years, and 82.9% of the patients were men. The detailed clinical information of these two groups is shown in Table 1.

Figure 1. Flowchart of the study design. GEO, Gene Expression Omnibus; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic; GSEA, gene set enrichment analysis.

Differential immune cell fraction

We first compared the relative abundance of 22 immune cell types between normal and oesophageal cancer tissues from the GEO database. In general, the top five immune cells in cancer and normal tissues were plasma cells, CD8⁺ T cells, dendritic cells (resting), regulatory T cells (Tregs), and mast cells (activated); the total fraction averaged at 64.8% in cancer tissues and 64.0% in normal tissues (Figure 2A). Moreover, dendritic cells (resting), M0 and M1 macrophages, mast cells (activated), memory B cells, CD4⁺ naive T cells, and CD4⁺ memory T cells (activated) were significantly enriched in cancer tissues relative to normal tissues, while mast cells (resting), M2 macrophages, naive B cells, CD8⁺ T cells, gamma and delta T cells, plasma cells, Tregs, and monocytes were significantly reduced in cancer tissues (Figure 2A, 2B, 2D). Immune cell fractions calculated with xCell showed similar results. The analysis validated the reliability of the CIBERSORT results and showed that the immune cell profile could be used as a potential prognostic marker, as it effectively distinguished normal tissue and cancer tissue (Figure 2D). Correlation analysis showed weak to moderate correlations between different immune cell types in the cancer group. CD8⁺ T cells were positively correlated with CD4⁺ memory T cells (activated) (r = 0.65, P < 0.001) and were negatively correlated with CD4⁺ memory T cells (resting) (r = -0.46, P < 0.001) (Figure 2C).

Figure 2. Different immune cell profiles between normal tissue and cancer tissue. (A) Relative immune cell fractions in the normal and cancer groups estimated with CIBERSORT based on gene expression profile data (GSE53625). (B) Heat map of differential immune cell fractions between the normal and cancer groups (FDR < 0.05). P values were calculated with the Mann–Whitney U test and adjusted for multiple testing (FDR). The left side bar shows the type of immune cells in A. (C) Correlation matrix of 22 immune cell types in the cancer group. Correlation coefficients were used to sort the cells by hierarchical clustering. (D) Violin plots of the abundance of CD8⁺ T cells and M1 macrophages calculated with CIBERSORT and xCell. The box plots in the violin indicate the median and interquartile range of the data distribution.

Feature selection and immunoscore calculation

For the training cohort, each immune cell fraction was dichotomized based on the optimal cut-off value generated by the survminer package (Supplementary Table 1). The LASSO Cox regression analysis identified 12 non-zero coefficient features for the calculation of prognostic immunoscores (Figure 3A, 3B). The immunoscore calculation formula is as follows: Immunoscore = (-0.294 × Memory B cells) + (-0.567 × Plasma cells) + (-0.257 × CD8⁺ T cells) + (0.068 × CD4⁺ Memory T cells (resting)) + (0.161 × Regulatory T cells (Tregs)) + (-0.043 × Gamma delta T cells) + (-0.410 × Monocytes) + (-0.251 × M0 Macrophages) + (-0.354 × M1 Macrophages) + (0.301 × Dendritic cells (resting)) + (-0.440 × Dendritic cells (activated)) + (-0.236 × Neutrophils), where 0 or 1 is assigned based on whether the immune cell fraction is below or above the corresponding cut-off value, respectively. For the training cohort, correlation analysis showed that a lower immunoscore was correlated with a higher survival rate of oesophageal cancer patients (Figure 4A). A time-dependent receiver operating characteristic (ROC) curve was used to evaluate the prognostic accuracy of 12 immune cell types over time. The AUC of the prognostic model was 0.733 at year 2, 0.736 at year 3, and 0.747 at year 5 (Figure 4B). Next, the training cohort was divided into a high-immunoscore group and a low-immunoscore group based on the median score (hazard ratio [HR], 2.70; 95% confidence interval [CI], 1.74-4.20; P < 0.001), which had 5-year survival rates of 23.7% (95% CI, 14.1-33.2) and 55.0% (95% CI, 43.0-67.0), respectively (Figure 4C). Stratification analysis of different clinical characteristics showed that immunoscore was significantly correlated with survival status (Figure 5 and Supplementary Table 2).

Figure 3. Immune feature selection based on the LASSO Cox regression model. (A) Optimal feature subsets based on the minimum criteria with 10-fold cross-validation of the LASSO regression. The vertical dashed lines represent the optimal values using the minimum criteria and the 1-standard error (SE) of the minimum criteria. The λ selected with 10-fold cross-validation (minimum criteria) was 0.039, and log (λ) was -3.245. (B) Regression coefficient profiles of 22 immune signatures in the LASSO model. The vertical dotted line is the optimal λ selected with 10-fold cross-validation, which resulted in 12 non-zero coefficients.

Figure 4. The immunoscore for each patient, the time-dependent ROC curve, and the Kaplan–Meier curves in (A–C) the training cohort and (D–F) the validation cohort. The 2-, 3-, and 5-year AUCs were used in the training cohort, and the 1-, 2-, and 3-year AUCs were used in the validation cohort to evaluate the accuracy of the prognostic model. The log-rank test was performed to analyse survival status.

Figure 5. Kaplan–Meier survival analysis based on the immunoscore signature of oesophageal cancer patients in the training cohort stratified based on clinical characteristics. (A, B) Age. (C, D) Sex. (E, F) TNM stage. The log-rank test was performed to calculate P values.

Immunoscore validation in the TCGA cohort

To evaluate the prognostic value of the immunoscore-based prognostic model in different populations, we used the same formula to calculate the immunoscore in the validation cohort. As in the model above, a lower immunoscore was associated with a more favourable prognosis (Figure 4D). For the validation cohort, the AUCs were 0.649 at year 1, 0.664 at year 2, and 0.724 at year 3. As in the training cohort, the patients in the validation cohort were divided into a high-immunoscore group and a low-immunoscore group based on the median score (Figure 4E). The 3-year survival rates were 66.9% (95% CI, 49.7-84.1) in the low-immunoscore group and only 22.0% (95% CI, 6.1-37.9) in the high-immunoscore group (HR, 2.93; 95% CI, 1.56-5.50; P < 0.001) (Figure 4F).

Univariate and multivariate Cox regression analyses were performed with immunoscore as a continuous variable. After controlling for other clinical covariates, a significant between-group difference was still observed in the training cohort (HR, 4.78; 95% CI, 2.64-8.67; P < 0.001) and the validation cohort (HR, 2.15; 95% CI, 1.04-4.45; P = 0.040), suggesting that the immunoscore may be used as an independent predictor for the prognosis of oesophageal cancer (Table 2).

Nomogram plotting and validation

Cox regression analysis was performed to determine the immunoscore, which showed that age and tumour stage were independent predictors of the prognosis of oesophageal cancer. For quantitative prognostic predictions, we used the training cohort to build a model containing the above independent predictors, which is presented as a nomogram (Figure 6A). The calibration curve demonstrated good agreement between the prediction model and the ideal model for the prognosis of oesophageal cancer (Figure 6B). The C-index of the nomogram was 0.718 (95%, 0.668-0.768) in the training cohort. The C-index of the prediction model was 0.716 (95%, 0.612-0.820) in the validation cohort. Decision curve analysis showed that when the threshold probability was greater than 0.3, the use of the plotted nomogram for prognostic prediction had more net benefits than indiscriminate treatment or no treatment, suggesting that the nomogram was clinically useful (Figure 6C).

Figure 6. Nomogram plotting and evaluation. (A) Clinical characteristics and immunoscores were used to develop a nomogram for prediction of the 2-, 3-, and 5-year survival rates of oesophageal cancer patients. (B) The calibration curve demonstrated agreement between the predictive and observed outcomes for 2-, 3-, and 5-year survival. The 45-degree dashed line indicates a perfect prediction of the ideal model, while the dotted lines indicate the actual performance of the nomogram. The closer the dotted line matches the dashed line, the better the prediction accuracy. (C) Decision curve analysis of the nomogram and tumour stage for the 5-year risk among patients with oesophageal cancer.

Biological functions underlying the immunoscore

To further investigate the biological functions underlying the immunoscore, we divided the oesophageal cancer patients into the high-immunoscore group and the low-immunoscore group and then performed gene set enrichment analysis (GSEA) to identify potential biological pathways based on the gene expression data of oesophageal cancer. Many immune-related pathways, such as the inflammatory response, IL-6-JAK-STAT3 signalling pathway, interferon-gamma response, chemokine-signalling pathway, T-cell receptor signalling pathway, and cytokine–cytokine receptor interaction pathways, were significantly enriched in the low-immunoscore group (false discovery rate (FDR) < 0.05 for all; Figure 7A, 7B). Correlation analysis of the immunoscore signature and immune-related gene expression showed that a high immunoscore was significantly correlated with the expression of antigen processing and presentation and B-cell receptor signalling pathway genes, while a low immunoscore was correlated with the expression of T-cell receptor signalling pathway genes and cytokine receptors (Figure 7C, 7D and Supplementary Table 3). These high correlations with immune-related pathways provide clues for further molecular mechanism studies related to the immunoscore.

Figure 7. Biological functions underlying the immunoscore. GSEA of (A) h.all.v7.1.symbols and (B) c2.cp.kegg.v7.1.symbols identified the most significant pathways and processes related to a low immunoscore. Correlation analysis of (C) a high immunoscore and (D) a low immunoscore with immune-related gene expression. The red line indicates a positive correlation, and the green line indicates a negative correlation. A darker colour corresponds to a stronger correlation.

Data source

The expression profile and clinical data of oesophageal cancer were retrieved and downloaded from the public (https://www.ncbi.nlm.nih.gov/geo/) database with a query of GSE53625 [24]. The expression data were generated with the Agilent human lncRNA + mRNA array v2.0 platform. Before the analysis, we re-annotated Agilent array probes. The probe sequences were aligned with SeqMap [25] from the GENCODE database (GRCh38, release 29) [26]. Only probes that were uniquely mapped to the genome with no mismatch were retained. Validation was performed on an independent cohort of oesophageal cancer patients from TCGA. The high-throughput sequencing (HTSeq)-fragments per kilobase of transcript per million mapped reads (FPKM) transcriptome data of the TCGA-ESCA project and survival information were obtained from Genomic Data Commons through the R package TCGAbiolinks [27].

Estimation of immune cell abundance

We used CIBERSORT to estimate the fractions of different immune cells in tissues based on the reference signature (LM22) matrix of 547 genes [10]. We uploaded the gene expression data prepared based on the example mixture file to the CIBERSORT web portal (https://cibersort.stanford.edu/) and ran the program with the default LM22 feature matrix at 1000 permutations. Only samples with significant CIBERSORT global deconvolution (P < 0.05) were included in subsequent analyses. For a given mixture sample, the relative fraction of immune cells was calculated with CIBERSORT and may be used directly for comparisons between immune cells or between studies.

The main results from CIBERSORT were validated with xCell (https://xcell.ucsf.edu/), a new proven method for counting cell subpopulations from a tissue expression profile [28]. xCell is a gene signature-based method and performs single-sample gene set enrichment analysis (ssGSEA) from the gene expression data for 64 immune and stromal cells.

Gene set enrichment analysis

To investigate biological pathways related to the immunoscore signature, we performed GSEA using the R clusterProfiler package on the gene expression profile [29]. The h.all.v7.1.symbols.gmt and c2.cp.kegg.v7.1.symbols reference gene set files used in this study were downloaded from the Molecular Signature Database (https://www.gsea-msigdb.org/gsea/index.jsp). Gene sets with an adjusted P < 0.05 after 1000 permutations were considered significantly enriched.

Statistical analysis

SPSS 22.0 and R 3.6.2 were used for statistical analysis. All statistical tests were two-tailed, and P < 0.05 was considered statistically significant. Count data were analysed with the chi-squared test or Fisher’s exact test. Normally distributed measurement data were analysed with the independent t-test, and non-normally distributed measurement data were analysed with the Mann–Whitney U test. To control for the FDR, the Benjamini-Hochberg method was used to adjust the P value when comparing immune cell fractions between the cancer group and the control group. Pearson correlation analysis was performed to identify any correlations between subsets of immune cells, and the corrplot package was used to visualize the resulting correlation matrix. The least absolute shrinkage and selection operator (LASSO) Cox regression method was used to screen the most significant immune cells for the prognosis of oesophageal cancer [30]. In the LASSO model, all immune cell fractions were dichotomized based on the optimal cut-off value calculated by the survminer package. The Kaplan–Meier method was used to analyse survival rates, and the log-rank test was used to analyse survival time. The time-dependent ROC curve was used to analyse the sensitivity and specificity of immunoscore-based survival predictions, and the area under the ROC curve (AUC) was quantified with the timeROC package [31]. The Cox proportional hazards model was used for univariate and multivariate prognostic analyses, and the regplot package was used to plot a nomogram for statistically independent predictors. The performance of the nomogram was evaluated with the calibration curve and Harrell's concordance index (C-index). In addition, decision curve analysis was performed by quantifying the net benefits of different threshold probabilities in evaluating clinical usefulness [32, 33].