Differential expressions of IRGs and TFs in CCA
The main design and process of this study are shown in the following figure (Fig. 1). We used RNA-seq data of 45 samples, including 36 CCA and 9 non-tumor tissues, collected from TCGA-CHOL. Totally, 1755 DEGs were identified, which consist of 743 significantly down-regulated genes and 1012 significantly up-regulated genes, for subsequent bioinformatics analysis (Fig. 2A and 2B).
The mRNA levels of 2498 IRGs (ImmPort database) and 318 TFs (Cistrome Cancer database) were examined. This analysis eventually revealed 156 IRGs (73 upregulated and 83 downregulated) and 18 TFs (11 upregulated and 7 downregulated) according to the screening criteria. Then, the heatmap and volcano plot were visualized to show the expression pattern of differentially expressed IRGs (Fig. 2C and 2D) and TFs (Additional file 1: Figure S1A and B) between CCA and non-tumor tissues. The names of DEGs were listed on Additional file 8: Table S3.
Construction of TF regulatory network
To investigate the biological roles and pathways of the DEGs, we found that top enriched GO terms were associated with regulation of signaling in immune and inflammatory response, such as “immune response”, “inflammatory response”, and “cytokine activity” (Fig. 2E). The significant pathways identified by KEGG pathway analysis included “cytokine- cytokine receptor interaction”, “natural killer cell mediated cytotoxicity” and “chemokine signaling pathway”.
Furthermore, GO enrichment analyses of differentially expressed IRGs demonstrated that the differentially expressed IRGs were mainly primarily participated in the pathway about tumor immunity, such as “regulation of immune effector process” and “regulation of cell killing” (Fig. 2F).
Identification of prognostic immune-related genes signature
For observing the association between IRGs and OS, we further conducted a Cox regression analysis to screen the IRGs. We found 11 prognostic IRGs (RAET1E, CST4, CCL24, CCK, CGB, GUCA2A, TDGF1, TDGF3, THPO, AVPR1B and IL9R) that were significantly related to survival of CCA patients in the training group (Fig. 3A). To further explore potential regulatory mechanisms behind deregulation of 11 prognostic IRGs in CCA, we analyzed the correlation between TFs and IRGs expression, as showed on Additional file 1: Figure S1C.
To investigate significance of risk genes in estimating prognosis, the 11 survival-related IRGs were further submitted to Lasso method and multivariate Cox regression analysis (Fig. 3B-C; Table 1). Finally, only 5 candidates prognostic IRGs (RAET1E, CST4, TDGF1, AVPR1B and IL9R) which may serve as significant predictors were obtained for risk model. Among these genes, RAET1E, CST4 and TDGF1 were identified as high-risk genes (serving as risk factors), while AVPR1B and IL9R were identified as low-risk genes (serving as protective factors). In the TCGA training cohort, the expression distribution of five IRGs and K-M curves for OS were shown in Additional file 2: Figure S2.
A prognostic index was established based on the expression levels and estimated regression coefficients of the risk genes. The formula of prognostic model containing five IRGs is described as follows: Risk score= (0.5631 × expression of RAET1E) + (0.4147 × expression of CST4) + (0.3941 × expression of TDGF1) + (-1.2478 × expression of AVPR1B) + (-1.4883 × expression of IL9R).
Analysis of the immune-related signature in the training set
The survival status of each patient in the training cohort was marked on dot plot, which showed that mortality rate in the high-risk group was significantly higher than that of the low-risk group (Fig. 4A and 4B). Among high- risk patients, three risk genes (RAET1E, CST4 and TDGF1) were upregulated, while two protective genes (AVPR1B and IL9R) were downregulated (Fig. 4C). In low- risk patients, these genes displayed an opposite expression pattern. The survival risk curve showed that patients with high risk had a significantly shorter OS time than those with low risk (P <0.001, Fig. 4D). Interestingly, we also explored the difference of DFS between two groups, and found the similar results (Fig. 4E). Then ROC analysis revealed reliable performance in the survival prediction of the model, with the AUC value of ROC at 1-, 3-, 5- year were 0.920, 0.966 and 0.940, respectively (Fig. 4F).
Analysis of the immune-related signature in the validation set
The GSE107943 cohort including 30 CCA samples was used for validation of prognostic model. The risk score distribution, survival status and risk gene expression in the validation set were displayed in Fig. 5A-C. The K-M curves demonstrated that the OS (P <0.05; Fig. 5D) and DFS (P <0.01; Fig. 5E) was significantly poorer in high-risk group than low-risk group. The AUC at 1-, 3-, 5- year in the validation set were 0.689, 0.713 and 0.710, respectively (Fig. 5F). Taking together, these results indicated the risk score could well assess patients’ prognosis. K-M analysis of OS among patients with different expression of each IRGs were displayed in Additional file 3: Figure S3.
Independent prognostic value of the risk model
Next, combined with other clinicopathological factors, we performed Cox regression analysis to assess performance of our model, such as age, gender and AJCC stage. The univariate analysis and multivariate analysis revealed that risk score was an independent risk factor in the TCGA cohort (HR=16.459, 95%CI, 3.923-69.046, P <0.001) (Table 2).
Nomogram for the prediction of prognosis in the training cohort
To better estimate the prognosis of CCA patients, we established a nomogram that integrated prognostic risk score and clinical variables (age, gender and AJCC stage) (Fig. 6A). The nomogram could accurately predict OS, and demonstrated that risk score contributed much more risk points than other variables. C-index for OS was 0.814, and calibration plot of internal verification for probability of survival at 1-, 2- and 3-year showed better conformity between the predicted and realistic observation results (Fig. 6B-D).
The correlation between the risk score and tumor-infiltrating immune cells (TIICs)
To determine whether our IRGs signature could reflect the proportion of TIICs on tumor immune microenvironment, we investigated the correlation between risk score and 22 immune cells in CCA patients. The difference on proportion of 22 TIICs types between high- and low- risk group were investigated (Fig. 7A and Additional file 4: Figure S4). As shown on the violin plot, we found the low- risk group had more infiltrated macrophages M1 compared to high- risk group. In addition, we also analyzed the association between immune cells subpopulations, which indicated that memory resting CD4+ T cell is negatively related to CD8+ T cell, while active NK cell is negatively related to resting NK cell in CCA patients (Fig. 7B). To analyze whether prognostic IRGs might contribute to TIICs, we analyzed correlation between IRGs and tumor-immune cell proportions (Additional file 5: Figure S5). Interestingly, we found the expression of RAET1E is positively correlated to the proportion of B cells naïve; TDGF1 is positively correlated to memory active CD4+ T cells. IL9R is positively correlated to follicular helper T cells, and negatively correlated to M2 Macrophages (Fig. 7C-F).
Drug Sensitivity prediction
We used GDSC database to predict the likelihood of response to several common chemotherapy and targeted drugs. We estimated the IC50 of each sample and observed a significant difference of IC50 between high-risk and low-risk groups among anticancer drugs listed on the GDSC database. Patients in the high-risk group were more sensitive to commonly administered chemodrugs (P =0.037 for AMG 706, P =0.019 for Thapsigargin) (Fig. 8A-D). In contrast, the chemotherapeutic response of cisplatin and gemcitabine was not significantly different between both groups.
External cohort validation for 5- IRGs signature
Finally, the 5-gene immune related signature was also validated in our cohort. Compared with adjacent normal tissues, the expressions of RAET1E, CST4, AVPR1B and IL9R were significantly higher (Fig. 8E), while TDGF1 were significantly lower in CCA samples.