Function annotation and pathway enrichment result of CCR genes.
It was revealed in the gene functional enrichment analysis that CCR genes was enriched in chemotaxis, positive regulation of cytosolic calcium ion concentration, chemokine-mediated signaling pathway, immune response, dendritic cell chemotaxis, cellular defense response, and so on (figure 1A). The correspondence between CCRs and GO terms was shown in figure 1B. The details of the enriched Gene Ontology (GO) terms in molecular function (MF), biological process (BP) and cellular component (CC) categories and KEGG pathway for CCR genes was displayed in table S1.
Expression of CCRs in HCC and para-carcinoma tissues
It was observed in GSE14520 cohort that expression of CCR1, CCR2, CCR3, CCR5, CCR7 and CCR8 in HCC tissues were significantly lower than para-carcinoma tissues, whereas expression of CCR6 and CCR9 was higher in HCC tissues (figure 2A). CCR4 and CCR10 are the only two members of the CCR family that show no difference in expression between HCC and para-carcinoma tissues. Expression correlation analysis between any two members of the CCR family showed that there were strong correlations among expression of CCR1, CCR2, CCR5, and CCR7 in HCC (figure 2B).
Then, we further evaluated the expression characteristics of CCR family genes in HCC in TCGA cohort. It was observed that expression of CCR1, CCR2, CCR4, CCR5, CCR7 and CCR9 were significantly lower in HCC tissues, whereas expression of CCR3, CCR8 and CCR10 were significantly higher in HCC tissues (figure 2C). Expression correlation analysis indicated that there were relatively high expression correlations among CCR1, CCR2, CCR4, CCR5, CCR6, CCR7 and CCR8 in HCC (figure 2D).
Diagnostic significance of CCRs in HCC
After a preliminary exploration of the expression characteristics of CCR gene family members in HCC, we assessed the possibility of these genes as diagnostic markers of HCC using the area under the ROC curve. In GSE14520 cohort, CCR1 (AUC=0.731, figure 3A) and CCR5 (AUC=0.714, figure 3E) was observed to be with good diagnostic performance in HCC, while diagnostic significance of the other CCR family members (figure 3B-D, F-J) were not satisfactory. In TCGA cohort, CCR1 (AUC=0.833, figure 3K) and CCR9 (AUC=0.835, figure 3S) was found to be with good diagnostic performance in HCC, while diagnostic significance of the other CCR family members (figure 3L-R, T) were not satisfactory.
Survival analysis result in GSE14520 and TCGA
In addition to whole-transcriptome microarray data and prognostic data, clinical information on 212 HCC patients was obtained from the GSE14520 dataset. In order to adjust for the effect of clinical factors in subsequent survival analyses, we first investigated the relationship between clinical factors and prognosis. The baseline information about the 212 HCC patients was displayed in table S2. It revealed that tumor size, cirrhosis, BCLC stage, TNM stage and AFP were associated with the OS of HCC, and tumor size, gender, TNM stage and BCLC stage were associated with the RFS of HCC.
We analyzed the relationship between CCR family members and RFS in GSE14520 and TCGA, respectively. In GSE14520 cohort, none of CCR gene was observed to be associated with RFS of patients in HCC, neither in univariate survival analysis nor after adjustment for clinical factors in Cox proportional hazards model (table 1, figure 4A-J); However, In TCGA cohort, CCR1, CCR2, CCR4, CCR5, CCR6, CCR7, CCR8 and CCR9 were observed to be associated with RFS of patients in HCC (figure 4K, L, N-S), while prognostic significance was not found for CCR3, CCR10 (figure 4M, T).
Then we evaluated the relationship between CCR family members and OS in GSE14520 and TCGA, respectively. Prognostic significance of CCR1 in OS (P=0.189, table 1, figure 5A) was not observed in univariate survival analysis; however, it (adjusted P=0.044, table 1) was observed to be associated with OS in Cox proportional hazards model after adjusted for clinical factors. CCR5 (P=0.022, adjusted P=0.021, table 1, figure 5E) and CCR7 (P=0.021, adjusted P=0.039, table 1, figure 5G) were both significantly correlated with OS in either Cox proportional hazards model or Kaplan Meier method in GSE14520 cohort. Other members of the CCR gene family were found to be associated with OS in HCC (figure 5B-D, F, H-J).
In TCGA cohort, CCR1, CCR2, CCR3, CCR4, CCR5 and CCR7 were observed to be associated with OS (figure 5K-O, Q), while prognostic significance for CCR6, CCR8, CCR9 and CCR10 were not observed (figure 5P, R-T).
Nomogram and prognostic signature
Based on the prognostic significance of CCR1, CCR5 and CCR7 found in our above study, in order to optimize our discovery and produce a better predictive prognostic model for patients with HCC, we respectively performed combined effect survival analysis, nomogram and prognostic signature in terms of the data of GSE14520. Combined analysis of CCR1 and CCR5 in HCC showed that patients in the group with low expression of CCR1 and CCR5 had the best outcome (figure 6A). Similarly, in other combined analyses, the patient in the group Ⅲ, the patients in the group c and the patients in group 3 all had the longest survival in their respective comparisons (figure 6B-D). The grouping protocols and outcomes were listed in table 2. We observed that the differences between the best and worst groups were more significant in the combined analysis than in the single gene survival analysis.
We established a nomogram and a prognosis signature based on the expression levels of CCR1, CCR5 and CCR7 in GSE14520. In nomogram, the length of corresponding line segment of each variable represents its contribution degree for prognosis. The parameter with the highest prognostic contribution was BCLC stage, followed by the degree of cirrhosis. The contribution of CCR1, CCR5 and CCR7 in predicting prognosis were similar (figure 6E). We evaluated the predictive power of the histogram by the match degree between the training group and the validation group. In the nomogram of GSE14520, there was a high degree of superposition between the self-validation cohort (red line) and training group (gray line) for predicting a 1-, 3 -, or 5-year prognosis (figure 6F-H).
The risk score formula for prognosis signature in GSE14520 was: risk score = expression value of CCR1 x -0.278 + expression value of CCR5 x -0.348 + expression value of CCR7 x -0.306. A total of 212 patients with HCC in GSE14520 were classified as high-risk group or low-risk group. Ranking patients by risk score from left to right (figure 6I, K), we observed that patients in the high-risk group had a higher concentration of individuals who reach terminal event in short term (figure 6J). The difference between the high and low risk groups in OS was statistically significant (P=0.025, figure 6L). Besides, the ROC curve also revealed that the prognostic signature worked well in predicting 1-, 2-, 3-, 4-, and 5-year outcome (figure 6M).
Validation for clinic significance of CCR1,5,7 in Guangxi cohort
Twenty-five patients from the Department of Hepatobiliary Surgery of Guangxi Medical University were enrolled as a validation cohort. The baseline data of the HCC patients in Guangxi cohort are listed in Table S3. In Guangxi cohort, IHC assay and qPCR assay showed that CCR1, CCR5 and CCR7 expression were significantly decreased in HCC tissues (figure 7A, B). Meanwhile, we observed that the expression levels of CCR1, CCR5 and CCR7 were strongly correlated (figure 7C). Besides, it was observed in Guangxi cohort that CCR1, CCR5 and CCR7 performed well in HCC diagnosis (figure 7D-F). In full agreement with the results in GSE14520 and TCGA database, CCR1 (P=0.045, table 3, figure 7G), CCR5 (P=0.013, table 3, figure 7H) and CCR7 (P=0.029, table 3, figure 7I) were significantly associated with prognosis of HCC, and their up-regulated expression predicting a good prognosis.
Nomogram and Prognostic signature construction in Guangxi cohort.
Based on the expression of CCR1, CCR5 and CCR7, we constructed the prognostic signature and the nomogram for HCC patients of Guangxi cohort. The risk score formula in Guangxi cohort was: risk score = expression value of CCR1 x -0.845 + expression value of CCR5 x -0.117 + expression value of CCR7 x -0.129. The risk score and the time of outcome event in HCC patients of Guangxi cohort of were displayed in the scatter plots (figure 8A, B), and the CCR1, CCR5 and CCR7 expression profile of these patients was shown using heat map. We observed that patients in the high-risk group had a shorter survival compared to those in the low-risk group. The results of survival analysis in the high and low risk groups indicated that the difference in prognosis was statistically significant (figure 8D, P=0.006). the survival ROC curve indicated that the prognostic signature worked well in predicting 1 year OS (figure 8E).
In the nomogram constructed in Guangxi cohort, the parameter with the highest prognostic contribution was tumor size, followed by the AFP (figure 8F). Predictive power of the nomogram was assessed using the match degree between the training group and the validation group. In the nomogram of Guangxi cohort, there was a high degree of superposition between the self-validation cohort (red line) and training group (gray line) for predicting a 1- or 2-year prognosis (figure 8G-H).
After intersecting the GSEA result of GSE14520 cohort with the GSEA result of TCGA cohort, it was observed that the results of these two datasets were very similar. Some of the more representative results were presented. It revealed that CCR1 (figure 9A, B) was associated with B cell receptor signal pathway, chemokine signaling pathway, nod-like receptor signal pathway, T cell receptor signal pathway, JAK-STAT signaling pathway, etc. CCR5 (figure 9C, D) was associated with B cell receptor signal pathway, chemokine signaling pathway, cytokine- cytokine receptor signal pathway, T cell receptor signal pathway, toll-like receptor signaling pathway, etc. CCR7 (figure 9E, F) was associated with B cell receptor signal pathway, chemokine signaling pathway, natural killer mediated cytotoxicity, nod-like receptor signal pathway, toll-like receptor signaling pathway, etc. We observed that these CCR genes were enriched in very similar pathways in HCC data sets, which suggested that there might be synergy between them.
Tumor-Infiltrating Immune Cells
We found a significant association between CCR1,5,7 and immune cell infiltration. CCR1 was positively correlated with the degree of B cells (Cor=0.498), CD8+ T cell (Cor=0.500), CD4+ T cell (Cor=0.389) and macrophage (Cor=0.629) infiltration in HCC tissues (figure 10A). CCR5 was also positively correlated with the degree of B cells (Cor=0.634), CD8+ T cell (Cor=0.680), CD4+ T cell (Cor=0.477) and macrophage (Cor=0.552) infiltration in HCC tissues (figure 10B). Similarly, the HCC tissues with high expression of CCR7 was companied with high degree infiltration of B cells (Cor=0.456), CD8+ T cell (Cor=0.405), CD4+ T cell (Cor=0.429) and macrophage (Cor=0.302) (figure 10C).