Immunotherapy: Open Access
Open Access

A Novel CCR7- Related Immune Prognostic Signature for Liver Hepatocellular Carcinoma (LIHC) Based on Immunogenomic Profiling

Author info »

Introduction

Liver Hepatocellular Carcinoma (LIHC) accounts for 85%-90% of primary liver cancer cases, and leads to poor survival with high morbidity and mortality [1]. Because of the delayed diagnosis and limitation of current therapeutic strategy, there are only 18% LIHC patients can survival 5 years after diagnosis [2]. LIHC is characterized by high degree of malignancy, insidious onset and rapid progression, and most patients are often diagnosed in advanced stage [2]. LIHC can be managed by various treatments, including surgical resection, radiofrequency ablation, chemotherapy and radiotherapy [3]. Targeted therapy and immunotherapy have made breakthroughs in recent years, but still with many unsolved problems and need further study [4].

Liver is the largest immune organ of human body, and there are various immune cells. LIHC is typically inflammatory associated tumor with autogenic immunogenicity, and immune tolerance and flight-avoidance mechanisms play an important role in this cancer progression [5]. The dynamical system was reported constituted by cancer cells, immune cell, non-tumor cells and extracellular matrix in tumor microenvironment, which promotes the progression of LIHC base on that composition of the system interact by producing growth factors, cytokines and chemokines [6]. Chemokines, one member of the cell super-factor family, work as attractants in cells migration, and involve in angiogenesis, metastasis, tumor cell growth and invasion [7]. Chemokines are reported promote tumor directly or perform its’ function of immune cell mediation, as well as prevent apoptosis of tumor cells [8]. In addition, the chemokine ligand-receptor interaction perform functions involve in surveillance and inflammation [9].

Chemokine receptor 7 (CCR7), one of G protein-coupled receptors, contributes to tumor neovascularization and accelerates tumor cells proliferation and migration in LIHC [10]. Study report that CCR7 play important roles in lymphatic metastasis, which is closely related to poor prognosis of LIHC patients [11]. CCR7 also was reported participate in epithelial-mesenchymal transition of LIHC cells via impeding histone deacetylase [10]. Besides, regulatory T-cells are key cells that play an immunosuppressive role in the immune microenvironment of LIHC, and it’s function is closely related to CCR7 and SLC (secondary lymphoid-tissue chemokine) [12]. CCR7 can also up-regulate Bcl-2 expression to protect CD8⁺ T-cell form apoptosis and promote CD3⁺ T-cells to inhibit proliferation of tumor cell [13]. Therefore, CCR7 was supported to influence immune microenvironment and correlate to prognosis of LIHC patients.

In this study, the LIHC expression data were downloaded from TCGA (The Cancer Genome Atlas) database, and the samples were divided into CCR7 high-expression group (CCR7_high) and CCR7 low-expression group (CCR7_L). The differences between those two groups were evaluated by calculating the level of ssGSEA (single-sample gene-set enrichment analysis) to 29 immune markers. The Differentially Expressed Immune Genes (DEIGs) were also identified between CCR7_H and CCR7_L groups. Then, the prognosis model was constructed based on LASSO regression for LIHC patients, as well as the nomogram model.

Materials and Methods

Gene expression datasets and data processing

There are 374 RNA-seq data and 377 clinical data of LIHC samples were downloaded from TCGA website (https://portal.gdc.cancer.gov/repository). RNA-seq data was normalized for FPKM values by DEseq2 package of R software. We divided those samples into CCR7_H and CCR7_L groups based on the expression of CCR7.

Immunogenomic analysis and survival analysis

The enrichment levels of the 29-immune signature in each sample were quantized using the ssGSEA score [14,15]. Then, the immune cell infiltration level (immune score), stromal content (stromal score), ESTIMATE score and tumor purity were calculated for each LICH samples [16]. The comparison of expression of HLA genes and immune checkpoint genes between CCR7_H and CCR7_L were conducted by ANOVA test for further analysis. In addition, the overall survival of CCR7_H and CCR7_L in LIHC patients was compared. The log-rank test was used in order to calculate the significance of survival time differences with a threshold of P-value<0.05.

Differential expression analysis and functional enrichment analysis

The differential expression analysis between CCR7_H and CCR7_L genes were analyzed using limma package of R software with the threshold of |log fold change (FC)|>1 and adjusted P values<0.05 [17]. Then, DEIGs were pooled into online tools Metascape (https://metascape.org/gp/index.html#/main/step1) and DAVID (https://david.ncifcrf.gov/) for gene annotation, visualization and providing gene attributes ADDIN EN.CITE [18-21].

Construction and test of the Immune Prognostic Signature (IPS)

The LIHC samples were divided into train set and test set stochastically. The DEIGs in the train set were pooled into LASSO Cox regression to establish immune prognosis signature [19,20]. To estimate a risk score for each patient, an IPS was established by weighting the Cox regression coefficients. Using survminer package of R software, LIHC patients were divided into low-risk and high-risk group. Then, the survival curve and risk curve, especially ROC curve were calculated to evaluate the sensitivity and specificity of this model [21]. The Area Under the Curve (AUC) values from the ROC curves was calculated and testified the prognostic prediction ability of this IPS in the test set synchronously. The concordance index (C-index) was used to evaluate model discrimination.

Development of the nomogram

The independent prognostic ability of the IPS was assessed by univariate and multivariate Cox analyses. Based on the application of the result of Cox analyses, we used the rms package of R software to develop a novel nomogram. On the purpose of determining correction range, calibration plots of observed vs. predicted probabilities of 1-3 and 5-year Survival were also figured.

Results

Immunogenomic analyses between CCR7_H and CCR7_L

The ssGSEA score was conducted to analyze the 29 sets of immune-associated genes of LIHC samples, which represent levels of immune cells, functions or pathways. The ssGSEA results show that in CCR7_H group immune cells, functions and pathways enriched. The Immune score, Stromal score and ESTAMATE score are higher in CCR_H group, while tumor purity higher in CCR_L group (P<0.001). In general, the level of immune status contraries to tumor purity.

The expression of HLA genes and immune checkpoint genes were also analyzed for understand immunogenomic of LIHC. The HLA genes expressed higher in CCR7_H group, as well as immune chick point genes (P<0.001). The survival curve between CCR7_H and CCR7_L was also printed, and the results showed favorable survival rate in CCR7_H group, which has higher immune score and lower tumor purity than CCR7_L group. This means immune compromising condition of LIHC patients in CCR7_H group and potent positive effect of immune checkpoint inhibitors in LIHC patients in CCR7_H group (Figure 1).

Figure 1: Immunogenomic analyses between CCR7_H and CCR7_L. (A): The enrichment levels of the 29-immune signature by ssGSEA score in each GBM sample. ESTIMATE was used to evaluate Tumor purity, Stromal score and Immune score; (B): Comparison of the expression levels of HLA genes between CCR7_H and CCR7_L; (C): Comparison of the Immune score, Stromal score, ESTIMATE score, Tumor purity between CCR7_H and CCR7_L; (D): Comparison of immune checkpoint genes expression levels between CCR7_H and CCR7_L; (E): Survival analysis between CCR7_H and CCR7_L. .

Identification of differential expressed immune genes

After analyzed the expression levels of immune genes, there are 903 DEIG’s were identified by limma package. Then, the functions and pathways enrichment of DEIG’s was achieved through pooling those genes into Metascape and DAVID website. Results showed that the DEIG’s mainly enriched in complement activation, adaptive immune response, T-cell activation, lymphocyte differentiation and cytokine-cytokine receptor interaction. We also made a network of those terms to clarify the relationship among them, the terms with a similarity of >0.3 was chosen and connected by edges. Besides, when the criteria of no more than 15 items per cluster and no more than 250 items in total were met, the item with the best P value was selected from each of the 20 clusters (Figure 2).

Figure 2: Volcano plot of immune genes differentially expressed between CCR7_H and CCR7_L; (B) GO terms and KEGG pathways enrichment of DEIG’s; (C) coloured by cluster ID, where nodes that share the same cluster ID are typically close to each other in DEIG’s; (D) coloured by p-value, where terms containing more genes tend to have a more significant P-value in DEIG’s. .

Construction of the CCR7-related immune prognostic signature

The samples in train set were performed by LASSO Cox regression and the risk score for patients can be calculated through coefficients in Table 1. The prognostic signature has good model recognition ability with the C-index (0.755). According to the risk score, patients were divided into high risk group and low risk group both in train set and test set. The survival rate of LIHC patients with low risk score have favorable survival rate in train set, as well as in test set. The ROC curve analysis of IPS showed hopeful prognostic ability for LIHC patients with AUC of 0.763 in train set and 0.706 in test set. In addition, the risk curve and expression pattern of DEIGs in IPS were conducted, the result showed patients with high risk score suffer from shorter survival time and higher mortality rate (Figure 3).

Table 1: The coefficients of genes in IPS.

Figure 3: Construction of the CCR7-related immune prognostic signature. (A, B): LASSO Cox analysis identified four genes most correlated to overall survival in train set; (C): Kaplan–Meier curves of overall survival based on the IPS in train set and test set. (D) ROC curve analysis of the IPS; (E,F): Risk scores distribution, survival status of each patient, and heat maps of prognostic four-gene signature in train set and test set: (1): In the first picture, a red dot means a patient with high risks, and a green dot means a patient with low risks; (2): In the second picture, a red dot means a dead patient and a green dot means a patient who survived; (3): In the third picture, red means high risk, and green means low risks.

Establishment an IPS-based nomogram model

Before construct the nomogram model, Univariate Cox analysis and Multivariate Cox analysis were performed. The result showed that the IPS correlated with OS of LIHC patients (Hazard ratio: 1.059 (1.040-1.079), P<0.001) and is an independent prognostic factor (Hazard ratio: 1.056 (1.034-1.078), P<0.001) (Table 2). The nomogram model was constructed based on Multivariate Cox analysis. A good consistency was showed by the 1, 3 and 5 year calibration diagram (Figure 4).

Table 2: Univariate and multivariate Cox analyses indicated that IPS was an independent prognostic factor.

Figure 4: Construction of the nomogram model. (A): Univariate and multivariate Cox analyses indicated that IPS was significantly associated with OS; (B): Nomogram model for predicting the probability of 1, 3 and 5 years OS in LIHC; (C): Calibration plots of the nomogram for predicting the probability of OS at 1, 3 and 5 years; (D): ROC curve analysis of the nomogram.

Discussion

LIHC is characterized by high degree of malignancy, insidious onset and rapid progression, and contributes to poor survival with high morbidity and mortality [2]. The immune tolerance and flight-avoidance mechanisms were reported play important roles in the progression of LIHC [5], especially the functions of chemokine ligand-receptor interaction [22]. CCR7 was reported contributes to tumor neovascularization and accelerates tumor cells proliferation and migration in LIHC [10]. Therefore, CCR7 was supported to influence immune microenvironment and correlate to prognosis of LIHC patients.

In this study, the ssGSEA results pointed that the CCR7_H group enriched with more immune cells, functions and pathways. The higher immune score, stromal score and ESTIMATE score also were identified in CCR7_H group, as well as lower tumor purity. The expression levels of HLA genes and immune checkpoints genes also perform higher in CCR7_H group. Besides, the survival analysis showed that the CCR7_H group patients have favorable survival rate than CCR7_L group. This means higher immune activity related to CCR7_H, and it can result to a better prognosis. The differential expression analysis was conducted to identify 903 DEIG’s, and the functions and pathways of those DEIGs significantly enriched in complement activation, adaptive immune response, T-cell activation, lymphocyte differentiation and cytokine-cytokine receptor interaction. That explored the molecular mechanism of CCR7 in progression of LIHC.

The IPS of LIHC based on CCR7 was constructed by applying LASSO Cox analysis, and the IPS was evaluated as a good prognostic model with C-index, survival analysis, ROC curve and risk curve analysis. In addition, the Univariate and Multivariate Cox analysis were conducted and fund that the IPS was an independent prognostic factor for LIHC patients. There are 10 genes were identified as hub genes of the IPS, including IGHV4-59, SCML4, AKR1B10, LINC00426, TRGC1, F2RL2, TRBV10-3, SAMD9L, SLC4A10 and ROR2. IGHV4-59 was reported was the most mutation frequent IGHV gene and over expression in chronic lymphocytic leukemia and contributes to progression of this cancer [23] and SCML4 also participate in the mechanism of acute and chronic lymphocytic leukemia [24]. AKR1B10 can protect hepatocytes from damage induced by oxidoreduction and detoxification, the deficiency of AKR1B10 play important roles in development of LIHC especially inflammation related hepatocarcinogenesis [25]. Zhu et al. reported that LINC00426 contributes to immune related signature construction to predict LIHC patients prognosis [26]. TRGC1, T-cell receptor gamma 1 region, was pointed paly roles in functions of T-cells [27]. F2RL2 was reported associated with development of preneoplastic foci in liver and promote initiation of liver cancer [28]. TRBV10-3 was one of subtypes of TRBV10, which related to CD4⁺ T-cells of primary biliary cholangitis patients [29]. The inactivation of SAMD9L was pointed contributes to progression of B virus-associated hepatocellular carcinoma via facilitation G1-S [30]. SLC4A10 is one of the SLC4 members, which had been proved related to cancer, hypertension and drug addiction [31]. ROR2, the receptor tyrosine kinase-like orphan receptor 2, can promote secretion of Wnt5a to active Wnt signaling pathway in colon cancer cells [32], which associated with proliferation, migration and invasion of LIHC closely [33].

This study elaborates molecular mechanism and constructs IPS of LICH based on CCR7 expression, provide new ways of thinking for LIHC immune microenvironment studies. However, this study is a retrospective research and still has limitations such as lack of foresight. Besides, the specific functions and action mechanisms of those hub genes of the IPS still need further study.

Conclusion

Based on the differences of immunogenomic profiling between CCR7_H group and CCR7_L group, this study constructed the CCR7-related IPS for LIHC patients. The expression level of CCR7 related to prognosis of LIHC patients. The IPS is an independent prognostic factor for LIHC. Besides, the nomogram was constructed to predict survival rate of LIHC patients.

Funding

This study was supported by the foundation of the Henan Medical Science And Technology Public Relations Project (LHGJ202105534).

References

1. Zhou J, Sun HC, Wang Z, Cong WM, Wang JH, Zeng MS, et al. Guidelines for diagnosis and treatment of primary liver cancer in China (2017 Edition). Liver Cancer. 2018;7:235-260.

   [Crossref][Google Scholar][PubMed].

2. Forner A, Reig M, Bruix J. Hepatocellular carcinoma. Lancet. 2018;391:1301-1314.

   [Crossref][Google Scholar][PubMed].

3. Llovet JM, Montal R, Villanueva A. Randomized trials and endpoints in advanced HCC: Role of PFS as a surrogate of survival. J Hepatol. 2019;70:1262-1277.

   [Crossref][Google Scholar][PubMed]. 

4. Li W, Ni CF. Current status of the combination therapy of transarterial chemoembolization and local ablation for hepatocellular carcinoma. Abdom Radiol. 2019;44:2268-2275.

   [Crossref][Google Scholar][PubMed].  

5. Shlomai A, de Jong YP, Rice CM. Virus associated malignancies: The role of viral hepatitis in hepatocellular carcinoma. Semin Cancer Biol. 2014;26:78-88.

   [Crossref][Google Scholar][PubMed].  

6. Hato T, Goyal L, Greten TF, Duda DG, Zhu AX. Immune checkpoint blockade in hepatocellular carcinoma: Current progress and future directions. Hepatology. 2014;60:1776–82.

   [Crossref][Google Scholar][PubMed].  https://doi.org/10.1002/hep.27246.

7. Valeta-Magara A, Gadi A, Volta V, Walters B, Arju R, Giashuddin S, et al. Inflammatory breast cancer promotes development of M2 tumor-associated macrophages and cancer mesenchymal cells through a complex chemokine network. Cancer Res. 2019;79:3360-3371.

   [Crossref][Google Scholar][PubMed].  

8. Yanru W, Zhenyu B, Zhengchuan N, Qi Q, Chunmin L, Weiqiang Y. Transcriptomic analyses of chemokines reveal that down-regulation of XCR1 is associated with advanced hepatocellular carcinoma. Biochem Biophys Res Commun. 2018;496:1314-1321.

   [Crossref][Google Scholar][PubMed].  

9. Anders HJ, Vielhauer V, Schlöndorff D. Chemokines and chemokine receptors are involved in the resolution or progression of renal disease. Kidney Int. 2003;63:401-415.

   [Crossref][Google Scholar][PubMed].  

10. Yang L, Chang Y, Cao P. CCR7 preservation via histone deacetylase inhibition promotes epithelial-mesenchymal transition of hepatocellular carcinoma cells. Exp Cell Res.

    [Crossref][Google Scholar][PubMed]. 2018;371:231-237.

11. Yan C, Zhu ZG, Yu YY, Ji J, Zhang Y, Ji YB, et al. Expression of vascular endothelial growth factor C and chemokine receptor CCR7 in gastric carcinoma and their values in predicting lymph node metastasis. World J Gastroenterol. 2004;10:783-790.

    [Crossref][Google Scholar][PubMed].

12. Smigiel KS, Richards E, Srivastava S, Thomas KR, Dudda JC, Klonowski KD, et al. CCR7 provides localized access to IL-2 and defines homeostatically distinct regulatory T cell subsets. J Exp Med. 2014;211:121-136.

    [Crossref][Google Scholar].

13. Vignali DAA, Collison LW, Workman CJ. How regulatory T cells work. Nat Rev Immunol. 2008;8:523-532.

    [Crossref][Google Scholar][PubMed].

14. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009;462:108-112.

    [Crossref][Google Scholar][PubMed].

15. Hänzelmann S, Castelo R, Guinney J. GSVA: Gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics. 2013;14:14-17.

    [Crossref][Google Scholar][PubMed].

16. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612

    [Crossref][Google Scholar][PubMed].

17. Servant N, Gravier E, Gestraud P, Laurent C, Paccard C, Biton A, et al. EMA - A R package for Easy Microarray data analysis. BMC Res Notes. 2010;3-277.

    [Crossref][Google Scholar][PubMed].

18. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: Tool for the unification of biology. Nat Genet. 2000;25:25-29.

    [Crossref][Google Scholar][PubMed].

19. Tian MX, Liu WR, Wang H, Zhou YF, Jin L, Jiang XF, et al. Tissue-infiltrating lymphocytes signature predicts survival in patients with early/intermediate stage hepatocellular carcinoma. BMC Med. 2019;17:1341-1346.

    [Crossref][Google Scholar][PubMed].

20. Tibshirani R. The lasso method for variable selection in the cox model. Stat Med. 1997;16:385-395.

    [Crossref][Google Scholar][PubMed].

21. Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56:337-344.

    [Crossref][Google Scholar][PubMed].

22. Li W, Gomez E, Zhang Z. Immunohistochemical expression of stromal cell-derived factor-1 (SDF-1) and CXCR4 ligand receptor system in hepatocellular carcinoma. J Exp Clin Cancer Res. 2007;26:527-533.

    [Google Scholar][PubMed].

23. Stanganelli C, Travella A, Bezares R, Slavutsky I. Immunoglobulin gene rearrangements and mutational status in argentinian patients with chronic lymphocytic leukemia. Clin Lymphoma, Myeloma Leuk. 2013;13(4):447-457.

    [Crossref][Google Scholar][PubMed].

24. Zhang S, Zhang Q, Yin J, Wu X. Overlapped differentially expressed genes between acute lymphoblastic leukemia and chronic lymphocytic leukemia revealed potential key genes and pathways involved in leukemia. J Cell Biochem 2019;120:15980-15988.

    [Crossref][Google Scholar][PubMed].  

25. Liu Y, Zhang J, Liu H, Guan G, Zhang T, Wang L, et al. Compensatory upregulation of aldo-keto reductase 1B10 to protect hepatocytes against oxidative stress during hepatocarcinogenesis. Am J Cancer Res. 2019;9:2730-2748.

    [Google Scholar][PubMed].

26. Zhu R, Guo W, Xu X-J, Zhu L. An Integrating Immune-Related Signature to Improve Prognosis of Hepatocellular Carcinoma. Comput Math Methods Med. 2020;2020:1-13.

    [Crossref][Google Scholar][PubMed].  

27. Lefranc MP, Chuchana P, Dariavach P, Nguyen C, Huck S, Brockly F, et al. Molecular mapping of the human T cell receptor gamma (TRG) genes and linkage of the variable and constant regions. Eur J Immunol. 1989;19:989-994.

    [Crossref][Google Scholar][PubMed].

28. Lupu DS, Orozco LD, Wang Y, Cullen JM, Pellegrini M, Zeisel SH. Altered methylation of specific DNA loci in the liver of Bhmt-null mice results in repression of Iqgap2 and F2rl2 and is associated with development of preneoplastic foci. FASEB J. 2017;31:2090-2103.

    [Crossref][Google Scholar][PubMed].

29. Nakagawa R, Muroyama R, Saeki C, Oikawa T, Kaise Y, Koike K, et al. CD4+ T cells from patients with primary biliary cholangitis show T cell activation and differentially expressed T-cell receptor repertoires. Hepatol Res. 2019;49:653-662.

    [Crossref][Google Scholar][PubMed].  

30. Wang Q, Zhai YY, Dai JH, Li KY, Deng Q, Han ZG. SAMD9L inactivation promotes cell proliferation via facilitating G1-S transition in hepatitis B virus-associated hepatocellular carcinoma. Int J Biol Sci. 2014;10:807-816.

    [Crossref][Google Scholar][PubMed].

31. Aalkjaer C, Boedtkjer E, Choi I, Lee S. Cation-coupled bicarbonate transporters. Compr Physiol. 2014;4:1605-1637.

    [Crossref][Google Scholar][PubMed].

32. MacLeod RJ, Hayes M, Pacheco I. Wnt5a secretion stimulated by the extracellular calcium-sensing receptor inhibits defective Wnt signaling in colon cancer cells. Am J Physiol - Gastrointest Liver Physiol. 2007;293.

    [Crossref][Google Scholar][PubMed].  

33. Xu M, Hu J, Zhou B, Zhong Y, Lin N, Xu R. TRIM29 prevents hepatocellular carcinoma progression by inhibiting Wnt/β-catenin signaling pathway. Acta Biochim Biophys Sin (Shanghai). 2018;51:68-77.

    [Crossref][Google Scholar].