Exploring immune-related signatures for predicting immunotherapeutic responsiveness, prognosis, and diagnosis of patients with colon cancer

Acquiring and preprocessing the relevant data

An anti-PD-L1 cohort (IMvigor210) was obtained using the IMvigor210CoreBiologies R package [8], which was treated with atezolizumab (anti-PD-L1 agent); only the samples with immunotherapeutic response information were included. The samples that displayed the complete response (CR) or partial response (PR) were categorized as responders and the samples that displayed stable (SD) or progressive disease (PD) were categorized as non-responders. In addition, the scores (estimate, stromal, and immune) of each sample were calculated using the estimate R package [9]. The differences in the clinical features and the scores between responders and non-responders were statistically compared using the Wilcoxon test, and the detailed information is shown in Table 1 and Supplementary Figure 1A–1C.

A total of 471 COAD samples were downloaded from UCSC Xena platform (https://xenabrowser.net/datapages/), including the mRNA expression data, mutation profiling data, and clinical and survival information. Subsequently, the samples (n = 471) were divided into the normal group (n = 39) and tumor group (n = 432); the tumor group was used as the training dataset, and the detailed information is listed in Supplementary Table 1.

Two independent GEO datasets (GSE39582 and GSE17536), considered testing datasets, were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE39582 included 556 colon cancer patients, and GSE17536 included 177 colon cancer patients, with available clinical and survival information. The detailed information was shown in Supplementary Tables 2 and 3.

A total of 1509 IRGs were obtained from the ImmPort database (https://immport.niaid.nih.gov/).

Identification of immune-related genes from the anti-PD-L1 cohort

The edgeR R package [10] was adopted to compare the gene expression differences between the responders and non-responders. Genes that met the cut-off criteria of false discovery rate (FDR) <0.05 and |log2 fold change (FC)|>1 were considered immunotherapeutic response related DEGs, and then intersected with IRGs to obtain the IRIGs. Subsequently, the Gene ontology (GO) functional, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed on IRIGs using the clusterProfiler R package [11].

Construction of estimate-stromal-immune scores model

The ESTIMATE algorithm calculated the immune, stromal and estimate scores for each colon cancer sample of the TCGA-COAD cohort [9]. Then, the colon cancer patients were divided into the high-score and low-score groups based on the median levels of the estimate, stromal, and immune scores, respectively. The association between the estimate-stromal-immune scores and OS was analyzed using the Kaplan-Meier curve. The pROC R package was used to plot the receiver operating characteristic (ROC) curves to evaluate the predictive capacities of the estimate-stromal-immune scores model.

Identification and validation of immune-related signatures for predicting immunotherapeutic responsiveness and prognosis of patients with colon cancer

The univariable Cox regression analysis was performed on the IRIGs using the training dataset with P < 0.05 as the criteria to identify the prognostic signature. Subsequently, we established the IRIRScore model for predicting the prognosis of colon cancer patients using the multivariate Cox regression analysis and the coefficient of the prognostic indicator was obtained. The IRIRScore for each patient was calculated according to the formula:

$$\text{IRIRScore}=\sum Coxcoefficientofgene{\chi }^i\times scaleexpressionvalueofgene{\chi }^i$$

Moreover, two independent GEO datasets with accession numbers GSE39582 and GSE17536 and anti-PD-L1 cohort were used to further validate the prognostic signatures of the constructed IRIRScore model. Besides, the Kruskal-Wallis test was used to compare the difference in IRIRScore between responders and non-responders in the anti-PD-L1 cohort.

In addition, a nomogram was designed to visualize the prognostic value of different patients’ features, including IRIRScore, StromalScore, ImmuneScore, age, gender, and tumor stage. The calibration curves were plotted to evaluate the predicted probabilities compared to the ideal predictive line using the rms R package. Moreover, the hazard ratios of the patients’ features were illustrated by a forest plot.

Statistical analysis

To verify and access the predictive capabilities of the constructed IRIRScore model, the ROC curves and Kaplan-Meier curves were plotted using the training dataset (tumor samples in TCGA-COAD cohort) and testing dataset (GSE39582, GSE17536, and anti-PD-L1 cohort). The optimal cutoff value for grouping patients were chosen at the points of the ROC curves where the difference between true positive and false positive was the most significant. Additionally, the associations between the risk score level and clinical features, including tumor stages and microsatellite status, were performed using the Kruskal-Wallis test.

Exploring potential immunotherapy related signatures in the tumor microenvironment of colon cancer

Firstly, the training dataset was divided into a high-risk group and a low-risk group based on the constructed IRIRScore model. Moreover, the proportions of 22 types of tumor-infiltrating immune cells of each sample in the TCGA-COAD cohort was estimated using the CIBERSORT algorithm [12]. Then, the statistical difference in immune landscape between the high-risk and low-risk score groups was compared using an unpaired t-test. The association between the significant differential immune cell types (P-value < 0.01) and OS were further accessed using Kaplan-Meier curve. Moreover, the correlations between tumor-infiltrating immune cells and IRIRScore were analyzed using the cor.test in R software (version 4.0.2, http://www.R-project.org).

In addition, the maftools R package was utilized to calculate the TMB value and visualize the mutation profiles of the high- and low-risk groups [13]. Subsequently, the TMB value between high- and low-risk groups was compared by the Kruskal-Wallis test, and the differences in the mRNA levels of immune checkpoints and their ligands between the high- and low-risk groups were statistically compared using the Wilcoxon test.

Preliminary exploration of IRIGs as predictive diagnostic indicators base on TCGA-COAD cohort

To explore the potential application of the IRIGs in the diagnosis prediction of colon cancer, we compared the expressions of the IRIGs between tumor and normal samples using the Wilcoxon test. Moreover, the IRIGs were employed as independent biomarkers to establish diagnosis prediction models using Random Forest (RF) algorithm. Besides, we tried to use the IRIGs as keywords to search immunohistochemical images in the Human Protein Atlas (HPA) (https://www.proteinatlas.org/), which may experimentally verify our findings.

Data availability statement

The immunotherapeutic cohort (IMvigor210) was available according to the guideline on http://research-pub.gene.com/IMvigor210CoreBiologies using the IMvigor210CoreBiologies R package. The TCGA-COAD cohort was available in the UCSC Xena (https://xenabrowser.net/datapages/). GSE39582 and GSE17536 were downloaded from NCBI-GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Identification of IRIGs

Figure 1 shows the whole analysis flow of this study. 1323 immunotherapeutic response-related DEGs (634 upregulated and 689 downregulated) were identified by comparing responders and non-responders in the anti-PD-L1 cohort (Figure 2A, Supplementary Table 4). After the intersection with 1509 IRGs, 108 IRIGs was obtained (Figure 2B, Supplementary Table 5). KEGG analysis results revealed that IRIGs were remarkably enriched in terms associated with Cytokine-cytokine receptor interaction, Neuroactive ligand-receptor interaction, JAK-STAT signaling pathway, cAMP signaling pathway, PI3K-Akt signaling pathway, and so on (Figure 2C). Furthermore, the GO enrichment analysis indicated that IRIGs mapped to hormone-related terms, including humoral immune response, antimicrobial humoral response, regulation of hormone secretion, and hormone transport (Figure 2D).

Construction of immune-related prognostic models

The AUC of estimate-stromal-immune scores model was less than 0.6 (ESTIMATEScore: 0.5; StromalScore: 0.513; ImmuneScore: 0.519, Supplementary Figure 1D), indicating its poor performance for predicting the prognosis of colon cancer. Additionally, the Kaplan-Meier curves revealed that the association between the scores (estimate, stromal, or immune) and OS were not significant (ESTIMATEScore: P = 0.61, StromalScore: P = 0.47, ImmuneScore: P = 0.94, Supplementary Figure 1E–1G).

Therefore, effective prognostic signatures for patients with colon cancer needed to be explored. Based on the training dataset, the univariate and multivariate Cox regression analyses were performed on the IRIGs, and 6 genes were found to be significantly related to OS status (Table 2). The coefficient and hazard ratio (HR) of each prognostic indicator is shown in Figure 2E and 2F, respectively. TDGF1, IL13RA2, and NRG1 are protective factors (HR <1), while CRLF1, GRP, and IL20RB are risk factors (HR >1). Then, the IRIRScore of each patient was calculated according to the formula described in the Materials and methods section. The scatter diagram in Figure 2G displayed that the survival time of patients in the low-risk group was longer than that in a high-risk group, which was also significantly supported by Figure 2H (P < 0.0001). The heatmap in Figure 2G reveals that the expression of risk factors was low in the low-risk group and high in the high-risk group, while the trend of protective factors was the opposite. The AUC of IRIRScore model at 3, 4, and 5 years of OS are 0.659, 0.677, and 0.653, respectively (Figure 2I). Patients with larger tumor sizes yielded higher IRIRScore than those with small tumor size (Kruskal-Wallis test, P = 3.9e⁻⁶, Figure 3A). Patients with metastases yielded higher IRIRScore than those without metastases (Kruskal-Wallis test, P = 0.012, Figure 3B). Patients with more distant or more lymph nodes are involved in patients with higher IRIRScore (Kruskal-Wallis test, P = 3.2e⁻⁵, Figure 3C). Patients with advanced stage yielded higher IRIRScore than those with early stage (Kruskal-Wallis test, P = 0.00016, Figure 3D). Figure 3E reveals that the patients with higher IRIRScore yielded more instability of the microsatellite status (Kruskal-Wallis test, P = 0.015, Figure 3E).

Validation and assessment of the identified immune-related signatures in predicting prognosis and immunotherapeutic responsiveness

To verify whether the identified immune-related signatures were suitable for different cohorts, The IRIRScore of each patient in two independent GEO datasets (GSE39582 and GSE17536), and the anti-PD-L1 cohort was calculated based on the same method used in the training dataset. Based on the testing datasets, the OS rate of patients in the high-risk group was significantly lower than that in low-risk group (P < 0.0001, Supplementary Figure 2A, Supplementary Figure 3A, and 3C), and the areas under the curves (AUCs) of IRIRScore model were calculated (GSE39582: 0.623 at 3-years, 0.607 at 4-years, and 0.563 at 5-years, Supplementary Figure 2B; GSE17536: 0.656 at 3-years, 0.733 at 4-years, and 0.657 at 5-years, Supplementary Figure 3B; Anti-PD-L1 cohort: 0.645, Supplementary Figure 3D).

Based on the testing data from GSE39582, the Kruskal-Wallis test revealed that higher IRIRScore was associated with higher T stage (P = 0.00038), metastasis (P = 0.011), N stages (P = 0.0015), and advanced pathological stage (P = 0.021) (Supplementary Figure 2C–2F).

In addition, based on the anti-PD-L1 cohort, we found that patients sensitive to immunotherapy had significantly lower risk scores than non-responders (P = 7.1e⁻⁵, Supplementary Figure 3E), indicating that responders will have better OS than non-responders.

Construction and validation of the nomogram

A prognostic nomogram was constructed to improve the accuracy of the performance of the IRIRScore model and to provide a quantitative and visual method for predicting 3-, 4-, and 5-years OS probability of patients with colon cancer. In the nomogram, it is easy to estimate the score for each variable on the point scale and calculate the probability of survival at 3-, 4-, 5-years (Figure 3F). Compared with the other variables, IRIRScore had the maximum score points. The forest plot shows that the variables, including the age (>60), advanced pathological stage (III and IV) and IRIRScore are significantly associated with OS (P < 0.05, Figure 3G), but the microsatellite status, StromalScore, and ImmuneScore are not (P > 0.05).

To verify the predictive performance of the nomogram, we plotted the calibration curves and observed that the predictive curves were close to the ideal curve (Figure 3H–3J). Moreover, the predictive accuracy (C-index) of the nomogram has been dramatically improved from 0.63 to 0.75, indicating good functioning.

Exploring potential immunotherapy related signatures

We estimated the proportions of 22 types of immune cells in patients with colon cancer in the training dataset using the CIBERSORT method. The profiles of high-risk group and low-risk group were shown in Supplementary Figure 3F and 3G. Subsequently, the composition of immune cell types between the high-risk and low-risk groups was compared using the Kruskal-Wallis test. Moreover, remarkable differences were found in T cells CD8, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), and Eosinophils (P < 0.01, Figure 4A). The Kaplan-Meier curve revealed that high level of T cells regulatory (Tregs) was significantly related to poor OS (P = 0.029, Figure 4C), while a high level of T cells CD4 memory resting was significantly related to better OS (P = 0.0071, Figure 4B). 10 of the 22 types of tumor-infiltrating immune cells were significantly related to IRIRScore (P < 0.05, Supplementary Figure 4). Among them, T cells regulatory (Tregs), Eosinophils, Macrophages M1, T cells follicular helper, T cells CD8, and Macrophages M2 were positively related to IRIRScore, while Dendritic cells activated, Plasma cells, T cells CD4 memory activated, and T cells CD4 memory resting were negatively related to IRIRScore.

The mutation profiles of the high-risk group and low-risk group were plotted and compared using maftools R package. The top 20 significantly mutated genes were APC, TP53, TTN, KRAS, SYNE1, MUC16, PIK3CA, FAT4, RYR2, OBSCN, ZFHX4, DNAH5, PCLO, CSMD3, ABCA13, DNAH11, LRP1B, FAT3, USH2A and CSMD1 (Figure 4D). Among them, the mutation rate of APC was significantly higher in the low-risk group (P = 0.000716, Figure 4E), while the mutation rate of MUC16 was significantly higher in the high-risk group (P = 0.0154, Figure 4E). However, there was no significant association between the mutation status of APC or MUC16 and patients’ OS (Figure 4F). Besides, the TMB value of each sample was calculated and visualized (median value of high-risk group: 2.66/MB; median value of low-risk group: 2/MB, Figure 4G). As shown in Figure 4H, the TMB value of patients in high-risk group was significantly higher than those in low-risk group (P = 0.00031). Moreover, the Wilcoxon test was applied to statistically compare the expression of immune checkpoint inhibitor targets between the two groups. The results revealed that the high-risk group had a significantly higher expression level (P < 0.001, Figure 5).

The identified IRIGs as potential diagnosis biomarkers

5 of 6 IRIGs were significantly differentially expressed between tumor and normal samples (P < 0.05, Figure 6A). Among them, TDGF1, and GRP were up-regulated in tumor samples, while NRG1, CRLF1, and IL20RB were down-regulated in tumor samples. The IRIGs were used to establish the prediction models with RF algorithm, and the AUC of each IRIG was calculated and plotted using the ggroc R package (IL13RA2: 0.579; TDGF1: 0.813; NRG1:0.933; CRLF1: 0.795; IL20RB: 0.625; GRP: 0.724, Figure 6B). Among them, TDGF1 and NRG1 revealed excellent diagnostic predictive performance (AUC >0.8). The distribution and expression of CRLF1, and NRG1 at the protein level are shown in Figure 6C–6F, whereas IL13RA2, TDGF1, IL20RB, and GRP remained inaccessible in HPA.