Development and Validation of a Novel Immune-Gene Pairs Prognostic Model Associated with CTNNB1 Alteration in Hepatocellular Carcinoma.

BACKGROUND Immunotherapy is one of the research hotspots in the field of hepatocellular carcinoma (HCC). Successive clinical trials have shown that patients with CTNNB1 mutations are resistant to immunotherapy, but the mechanism is still unclear. MATERIAL AND METHODS We identified differentially expressed immune genes (DEIGs) in patients with and without CTNNB1 mutations in the Cancer Genome Atlas (TCGA) database and then paired them to explore any correlation with prognosis. Univariate Cox regression analysis and Lasso regression analysis were used to develop the prognostic model. We first divided the TCGA cohort into 29 subgroups for internal validation and then used the International Cancer Genome Consortium (ICGC) cohort to conduct external validation. We also used a CIBERSORT algorithm to quantify immune infiltration of the different risk groups. RESULTS The novel prognostic model consisted of 45 immune-gene pairs with general applicability. It was more accurate than the traditional prognostic signature, which is based on gene expression by comparison of area under the receiver operating characteristic curve (AUC) values. The infiltration proportion of B cells, CD8 T lymphocytes, activated natural killer cells, and M1 macrophages in the low-risk group was greater in the high-risk group, while the infiltration proportion of M0 and M2 macrophages was greater in the high-risk group. CONCLUSIONS In this study, a novel approach was proposed for evaluating HCC prognosis, which may be useful in evaluatingthe intensity of the immune response in the HCC microenvironment.

Background

Despite the great progress made in treatment of hepatocellular carcinoma (HCC), the prognosis for patients with the disease remains poor [1-3]. In recent years, immunotherapies and particularly immune checkpoint inhibitors, has revolutionized treatment of many malignancies. HCC is a typical inflammation-related cancer that occurs in the presence of immunosuppression, making immunotherapy a potentially attractive treatment option [1,4-6]. However, some recent clinical trials have found that patients with CTNNB1 mutations in HCC tissues are not sensitive to immunotherapy [7-9]. The underlying mechanism is not entirely clear.

CTNNB1 is an important component of the typical Wnt signaling pathway and plays an indispensable role in physiological embryology, zone formation and metabolic regulation of the liver [10]. Activating and inactivating mutations in linker b1 (CTNNB1), which encodes the protein β-catenin, are associated with tumorigenesis of the liver in humans and [11] Approximately 11% to 41% of all liver malignancies are HCC with mutated CTNNB1,which has specific metabolic, morphological, and clinicopathological characteristics [12]. For example, tumors with CTNNB1 mutations often occur in the setting of a liver without cirrhosis, the infection rate of hepatitis b virus is low, and the differentiation of the tumor is good, but it is often accompanied by microvascular invasion, tumor capsule invasion, and intratumor cholestasis [13]. Mutations in CTNNB1 also are associated with malignant transformation of hepatocellular adenoma [10].

A recent study found that activation of β-catenin promoted immune escape as well as resistance to anti-programmed cell death protein-1 therapy for HCC. The mechanism may be related to defective dendritic cells (DC) and antigen-specific T-cell recruitment [9]. In a mouse model, activation of β-catenin led to poor T-cell infiltration, promoting immune escape and leading to resistance to immunotherapy.7 Although the relationship between CTNNB1 mutations and prognosis remains controversial, there is increasing evidence suggesting that changes in the tumor immune response can affect prognosis in patients with HCC [14,15].

In this study, we attempted to explore the immunophenotype of HCC from the perspective of mutations in CTNNB1. We identified differentially expressed immune genes (DEIGs) in samples from patients with and without CTNNB1 mutations in The Cancer Genome Atlas (TCGA) database and paired them to explore any correlation with prognosis. Creation of our prognostic model was based on the above DEIGs. In addition, the International Cancer Genome Consortium (ICGC) dataset was used to validate the accuracy and universal applicability exhibited by the prognostic model.

Material and Methods

Data collection

mRNA Sequence data from 374 HCC samples came from TCGA (). Corresponding clinical data – including survival time, survival status, sex, age, race, body mass index (BMI), histological grading, alpha fetoprotein (AFP) level, vascular invasion, American Joint Committee on Cancer tumor, node, metastasis (AJCC-TNM) stage, family history of cancer, prior malignancy history, and new tumors presenting following initial treatment, as well as tumor status – were taken from the University of California Santa Cruz’s Xena Funcational Genomics Exporer (). The CTNNB1 mutation sample list was obtained from the cBioPortal (). We downloaded another independent mRNA sequence dataset and corresponding clinical data from the ICGC () to conduct external validation. The two sequence datasets came from the Illumina HiSeq RNA-Seq platform. We downloaded the immune-related gene list from the ImmPort database (). It should be stated that acquisition of the above data followed the data access policies and release guidelines for these databases. Therefore, there was no requirement to obtain the approval from local ethics committee.

Identification of DEIRGs

The “edgeR” R package was used to identifying differentially expressed immune-related genes (DEIRGs) in 276 HCC samples free of CTNNB1 mutations and 96 HCC samples with CTNNB1 mutations in the TCGA cohort. A false discovery rate (FDR) of <.05 was considered significant.

Annotations of DEIRG biological processes

The “clusterProfiler” R package was used to annotate the function of the DEIRGs and conduct enrichment analysis of their involvement in biological processes (gene ontology). The threshold value was an FDR <.05.

Construction and validation of immune gene-pairs prognostic model

We then paired the DEIRGs. In each gene pair, if the expression level of the former gene was higher than that of the latter gene, the value assigned was 1. On the contrary, if the expression level of the former gene was lower than that of the latter gene, the value assigned was 0. Immune-gene pairs with values of 0 or 1 that comprised less than 20% of the TCGA cohort were excluded. A total of 14 485 immune gene pairs were obtained after screening. Univariate Cox regression analysis was applied to identify the immune gene pairs associated with prognosis (P<.001). The “|glmnet” R package was used for Least absolute shrinkage and selection operator (Lasso) regression analysis of the prognostic immune gene pairs. The penalty coefficient (log Lambda value) with the smallest cross-validation error helped to confirm the gene pair coefficient. The formula for the risk score was the sum of each coefficient of immune-gene-pair multiplied by each value in the immune-gene pair. The risk score corresponding to the maximum value for sensitivity and specificity under the receiver operating characteristic (ROC) curve was used to classify the groups with high and low risks. The R packages “survminer” and “survivalROC” were used to generate Kaplan-Meier survival curves and ROC curves regarding the risk score as a means of assessing the predictive power of the model. The log-rank test was used to compare the two groups with regard to survival curves, and P<.05 was considered statistically significant. Patients in TCGA were divided into several subgroups based on their clinical characteristics for internal verification, and external verification was conducted by using the ICGC cohort.

Independence validation of the prognostic model

Univariate and multivariate Cox regression analysis assisted in verifying whehter risk score could independent predict HCC prognosis in the TCGA and ICGC cohorts.

Analysis of association between risk score and clinicopathology

We used the Wilcoxon test or Kruskal’s algorithm (multiple groups) to analyzing the association between the risk score and clinicopathology (including AJCC-TNM stage, histologic grade, vascular tumor cell type, tumor status, and new tumor event after initiation of treatment). P<0.05 was deemed to be statistically significant. The “beeswarm” package in the R software was used to generate the boxplot.

Analysis of infiltration of immune cells in different risk groups

The CIBERSORT algorithm assesses cell composition in complex tissues based on standardized gene presentation data, quantifying the richness of specific kinds of cells [16]. We used this method (with a characteristics matrix of 547 genes standing for 22 kinds of infiltrating immune cells) to quantify immune cell infiltration in tumor tissues from groups with high and low risks, respectively. P<.05 was used as the standard by which to judge accuracy of the estimation.

Statistical analysis

Statistical analyses were conducted with R software v3.6.1 (R Foundation for Statistical Computing, Vienna, Austria) as well as with GraphPad Prism v7.00 (GraphPad Software Inc., United States). Analysis of qualitative variables was performed with Pearson’s χ2 test or Fisher’s exact test; analysis of quantitative variables was done with a non-parametric Wilcoxon test for unpaired samples in the appropriate mode. Analysis of various normalized data groups was done with Kruskal’s algorithm. P<.05 was considered statistically significant.

Results

Identification of differentially expressed immune genes

Immune cell infiltration in tumor tissues in the wild-type (WT) CTNNB1 group was significantly different from that in the mutated CTNNB1 group (Figure 1A). Among the 434 differentially expressed immune genes (the heat map shown in Figure 1B), 231 were upregulated in the non-CTNNB1-mutated group (n=276), and their biological processes mainly involved activation and regulation of leukocytes and lymphocytes (Figure 1C). We found that 203 genes were upregulated in the CTNNB1-mutated group (n=98), and their biological processes were mainly related to innate immune response and to regulation of the Fc receptor signaling pathway. (Figure 1D). This may help to explain the insensitivity of CTNNB1-mutated HCC to immunotherapy.

Figure 1

Identification of differential expressed immune-related genes (DEIRGs) between with and without CTNNB1 mutation HCC in TCGA. (A) Boxplot of immune cell infiltration in with and without CTNNB1 mutation HCC patients (P-value significant codes: 0< *** <0.001< ** <0.01< * <0.05). (B) Heatmap of DEIRGs in the significantly enriched gene sets in CTNNB1-mutant HCC. (C) The biological process Enrichment Map of DEIRGs upregulated in the WT-CTNNB1. (D) The biological process Enrichment Map of DEIRGs upregulated in the mutated-CTNNB1.

Construction of a robust prognostic model consisting of 45 immune gene pairs

Univariate Cox regression analysis showed that among the 14 485 gene pairs, 480 i were associated with prognosis (P<.001). After the overfitting was removed using Lasso regression, a robust prognostic model consisting of 45 immune gene pairs was established (Table 1). According to the ROC curve, the optimal critical value for dividing the group in to highand low risk was.886. Kaplan-Meier survival curve analysis showed that the group with high risk had lower rates of overall (OS) and disease-specific survival (DFS) relative to the group with low risk (Figure 2A, 2C). AUC values for 1-, 2-, 3-, 4-, and 5-year OS reached.899,.905,.892,.915, and.894, respectively (Figure 2B) and those for DFS were.956,.945,.936,.935, and.903, respectively (Figure 2D). As displayed in the survival status map, a risk score increase meant an obvious increase in the number of patient deaths (Figure 3A). As illustrated in univariate and multivariate Cox regression analyses, only risk score was an independent predictor of prognosis (Figure 3B, 3C).

Table 1

The immune-gene pair list and coefficient.

Immune-gene pair	Coef
RXRB|ROBO1	−0.37609
CDH1|GRB2	−0.15355
NEDD4|ITGAV	−0.20133
FABP3|CCL5	0.060678
TGFB1|CCL19	0.153498
S100A1|BIRC5	−0.02978
S100A9|JAK1	0.220469
S100A9|PSMD6	0.051975
MX1|PAK1	−0.19625
NFKBIE|GHR	0.01206
PTK2|AHNAK	0.055809
PLAUR|PIK3CD	0.181683
PLAUR|INPP5D	0.358956
IL1R1|BIRC5	−0.012
CLEC11A|CD48	0.132153
OAS1|ENG	0.275178
TNFAIP3|RORC	0.053256
TLR2|CD3E	0.031311
MMP9|RAC2	0.118758
MMP9|HLA-DQA1	0.094582
HLA-E|RHOA	−0.08382
IL15RA|CCL19	0.036259
IL15RA|RELB	0.258749
IL15RA|ITGAL	0.146002
IL15RA|IRF1	0.157378
IL15RA|AKT2	0.02658
SEMA3F|FYN	0.134619
PGF|IL7R	0.174782
SDC3|BLNK	0.195853
SDC3|CCL5	0.046152
SDC3|CD3E	0.141212
PSMD1|IL6ST	0.034187
GHR|LTBP1	−0.13369
GHR|BMP2	−0.07265
GHR|PROCR	−0.00744
MMP12|PDCD1	0.059706
EPOR|BIRC5	−0.04532
CCL21|EPO	−0.03773
BMP2|PSMD5	0.132738
RELB|CSF1	−0.23732
TAP2|CSF1	−0.01425
PSMD6|ENG	0.167224
MSR1|CCL5	0.00684
ENG|PSMD5	−0.20019
PLSCR1|FYN	0.280502

Figure 2

Survival analysis of the prognostic model. (A, B) The Kaplan-Meier survival curve of overall survival (OS) and disease-specifici survival (DSS) in TCGA cohort (C, D) The time-dependent ROC curve of overall survival (OS) and diseasespeciFIC survival (DSS) in the TCGA cohort

Figure 3

Independent validation of the prognostic model. (A) Survival status of patients in the TCGA cohort. (B) Univariate cox regression analyses of the association between clinicopathological factors and OS. (C) Multivariate Cox regression analyses of the association between clinicopathological factors and OS.

Internal validation based on clinical features of TCGA

To observe whether the model constructed was applicable to different populations, we divided the patients into 29 subgroups according to their clinical characteristics (including sex, age, race, BMI, alpha-fetoprotein levels, histological grade, vascular invasion, AJCC-TNM stage, family history of cancer, malignancy history, new tumors presenting following initial treatment, and individual tumor status) for internal verification. The survival curve showed that in each subgroup, the group with a high risk had lower OS relative to the group with low risk (Figure 4A–4L). This indicated that our prognostic model had universal applicability.

Figure 4

Internal validation of the prognostic model in the TCGA cohort according to clinical features. (A) AFP. (B) Age. (C) BMI. (D) Family cancer history. (E) Gender. (F) Tumor status. (G) New tumor event after treatment initiation. (H) Prior malignancy. (I) Histologic grade. (J) AJCC stage. (K) Vascular tumor cell type. (L) Race.

Correlation analysis between risk score and clinicopathology

With increasing histological grade, AJCC stage, and vascular invasion, patients’ median risk scores gradually increased. Median risks score in patients with new tumors after initial treatment were higher than in those without new tumors, and median risk scores in patients with a tumor was higher than those in individuals without tumors (Figure 5). We also carried out Chi-squared tests on the two groups, which revealed an obvious difference in clinical stage, histological grade, vascular invasion degree, and tumor histological grade between them, which may help to explain the differences in their prognoses (Tables 2, 3).

Figure 5

Correlation of risk score with clinicopathological characteristics. (A) Histopathological grade. (B) Vascular tumor cell type. (C) AJCC stage. (D) New tumor event after treatment initiation. (E) Cancer status.

Table 2

The chi-square test of the relation between risk score and clinical features in TCGA.

Clinical feature	Risk score		χ2	P
	High risk n(%)	Low risk n(%)
AFP			1.482	0.223
>300 ng/ml	22 (35.48%)	40 (64.52%)
≤300 ng/ml	54 (27.41%)	143 (72.59%)
Age			0.079	0.779
>65	46 (36.22%)	81 (63.78%)
≤65	75 (34.72%)	141 (65.28%)
BMI			4.461	0.035
>25	42 (27.45%)	111 (72.55%)
≤25	63 (38.65%)	100 (61.35%)
Family cancer history			0.149	0.699
No	29 (23.58%)	94 (76.42%)
Yes	19 (26.03%)	54 (73.97%)
Gender			0.282	0.595
Female	54 (42.90%)	69 (56.10%)
Male	80 (34.33%)	153 (65.67%)
Tumor status			14.215	0.001
Tumor free	45 (25.14%)	134 (74.86%)
With tumor	67 (44.97%)	82 (55.03%)
New tumor event after initiate treatment			9.012	0.003
No	43 (26.54%)	119 (73.46%)
Yes	71 (42.26%)	97 (57.74%)
Prior malignancy			2.406	0.121
No	114 (36.53%)	198 (63.46%)
Yes	7 (22.58%)	24 (77.42%)
Histologic grade			9.308	0.025
G1	10 (18.87%)	43 (81.13%)
G2	56 (34.78%)	105 (65.22%)
G3	48 (42.86%)	64 (57.14%)
G4	5 (41.67%)	7 (58.33%)
AJCC stage			26.373	0.001
I	37 (22.98%)	124 (77.02%)
II	29 (37.66%)	48 (62.34%)
III–IV	45 (56.25%)	35 (43.75%)
Vascular tumor cell type			10.5	0.005
None	50 (26.60%)	138 (73.40%)
Micro	29 (34.12%)	56 (65.88%)
Macro	10 (62.50%)	6 (37.50%)
Race			3.835	0.147
White	55 (32.54%)	114 (67.46%)
Asian	59 (39.86%)	89 (60.14%)
Blank & Indian	3 (18.75%)	13 (81.25%)

Table 3

The chi-square test of the relation between risk score and clinical features in ICGC.

Clinical feature	Risk score		χ2	P
	High risk n(%)	Low risk n(%)
Gender			0.538	0.463
Male	47 (27.81%)	122 (72.19%)
Female	20 (32.79%)	41 (67.21%)
Age			0.382	0.537
>65	39 (27.66%)	102 (72.34%)
≤65	28 (31.46%)	61 (68.54%)
Stage			30.071	0.001
I–II	23 (16.20%)	119 (83.80%)
III–IV	44 (50.00%)	44 (50.00%)
Prior malignancy			0.562	0.453
Yes	7 (23.33%)	23 (76.67%)
No	60 (30.00%)	140 (70.00%)

External validation regarding the prognostic model in the ICGC cohort

Patient risk scores in the ICGC cohort were calculated based on the risk score calculation formula from the TCGA cohort. The cutoff values of the groups with high and low risks were the same as those for the TCGA cohort. As with the TCGA cohort, the group with a high risk had a poor prognosis (Figure 6A). AUC values for 1-, 2-, 3-, 4-, and 5-year OS reached.831,.727,.730,.786 and.786, respectively (Figure 6B). As revealed by univariate and multivariate Cox regression analyses, risk score was an independent predictor of prognosis (Figure 6C, 6D).

Figure 6

External validation of the prognostic model in the ICGC cohort. (A) Kaplan-Meier curve of overall survival (OS) in the ICGC cohort. (B) Time-dependent ROC analysis in the ICGC cohort. (C) Univariate cox regression analyses of the association between clinicopathological factors and OS in the ICGC cohort. (D) Multivariate Cox regression analyses of the association between clinicopathological factors and OS in the ICGC cohort.

Infiltration of immune cells in tumor tissue between different risk groups

After screening (P<.05), 42 tumor tissues from the low-risk group and 38 tumor tissues from the high-risk group were included. The infiltration ratio for B cells, CD8 T lymphocytes, activated natural killer cells, and M1 macrophages in the group with low risk appeared higher relative to the group with high risk, while the latter group had a higher infiltration ratio for M0 and M2 macrophages relative to the former group (Figure 7).

Figure 7

The immune infiltration landscape in HCC patients with high- and low-risk. (A) The barplot of the proportion of immune cell infiltration. (B) Heatmap of the proportion of immune cell infiltration. (C) Violin plot of 22 types of immune cell infiltration abundances in groups with high and low risks (red and blue denote high and low risks).

Discussion

The prognosis for patients with HCC is still not satisfactory, mainly because most cases are only diagnosed when the disease is advanced [1]. Adjuvant therapy, such as immunotherapy, is one of the primary treatments for advanced HCC. However, several recent clinical trials have found that patients who have HCC with CTNNB1 mutations are not sensitive to immunotherapy [5,8,9]. The mechanism of CTNNB1 mutation affecting immunophenotypic regulation and prognosis of HCC remains unclear. Some immune signatures related to prognosis of HCC have been identified. For example, Long [17] identified a dual-immune-gene prognostic model, Wang [18] identified a 9-immune-gene prognostic signature, and Li [19] identified a 6-immune-gene prognostic signature. These reports provide theoretical evidence on which to base assessment of HCC prognosis, taking into account genes related to the immune system.

In this study, we found that several immune cells in tumor tissues with CTNNB1 mutation had a significantly decreased infiltration level, which is similar to previous reports [9,15]. In HCC without CTNNB1 mutations, the upregulated immune genes were significantly involved in regulation of lymphocyte activation. Therefore, immunotherapy resistance associated with a CTNNB1 mutation may be related to low lymphocyte activity.

After it was clear that the differentially expressed genes could significantly affect the immune response of the tumor, we combined them in pairs. This is a different approach than that in previous studies. Most previous studies have assessed the association between a single gene and prognosis based on its level of expression. With gene pairing, we were able to analyze the effect of the relative expression level of two genes on prognosis, transforming continuous variables into a dichotomous variable. With this methods, other centers using the model would not have to carry out batch normalization.

We screened 45 gene pairs from 480 immune gene pairs related to prognosis by means of lasso regression. The results showed that risk stratification of patients with based on these 45 immune gene pairs was very accurate, especially for evaluating the specific survival rate for HCC. It was important that the prognostic model had universal applicability and could be applied to people with different clinical characteristics. In addition, risk scores were associated with factors correlated with poor prognosis, such as low tumor differentiation, vascular invasion, and advanced clinical stage. This may help to explain the difference in prognosis between different groups.

Finally, we used the ICGC dataset to conduct external verification of the prognostic model. Conforming to the findings of the TCGA cohort, the prognostic model was highly accurate for predicting survival in the short and long term, and it was an independent predictor of prognosis. We compared AUC values to determine the accuracy of the prognostic model relative to other prognostic signatures. AUC values for Long’s double immune gene prognostic signature were.738,.727, and.655 at 1, 3, and 5 years, respectively [17]. AUC values for Wang’s 9-immune-gene prognostic signature were.811,.711, and.734 at 1, 3, and 5 years, respectively [18]. AUC values for Li’s 6-immune-gene prognostic signature were.761,.681, and.692 at 1, 3, and 5 years, respectively [19]. AUC values for our prognostic model at 1, 3, and 5 years were.899,.892, and.894, respectively. Obviously, our prognostic model has more value than any of the traditional prognostic immune signatures.

Tumor microenvironment (TME) is the cellular environment of tumor cells, and includes tumor stromal cells, the extracellular matrix, and soluble molecules. Once the TME is shaped, a lot of immune cells, such as T cells, macrophages, and myelogenic inhibitory cells, are drawn by chemotaxis to constitute it [20]. In the TME, immune and stromal cells, which are major types of non-tumor components, may have a role in diagnosis and prognostication of tumors [21,22]. As found in this study, prognosis for patients in the group with low risk was dramatically better than in the group with high risk. Garnelo [23] found that patients with higher proportions of tumor-infiltrating T and B cells had stronger local immune activity and better prognosis. This was consistent with our research results. In addition, patients whose activated NK-cell infiltration level was high had better prognosis. That is in accord with conclusions by Zhu and other researchers [24].

Tumor-associated macrophages (TAMs) are macrophages that infiltrate tumor tissues and they make up the largest number of immune cells in the TME [25,26]. M1 macrophages have an antitumor effect, which can distinguish tumor cells from normal cells and eventually kill tumor cells by identifying them. M2 macrophages mainly promote tumor development, invasion, and metastasis [27,28]. In patients who had a high level of infiltration of M2 macrophages, the prognosis was poor, while for those with high levels of M1 macrophage infiltration, the prognosis was better. Therefore, regulating the phenotype of TAM may function as a positive factor in improving prognosis in patients with HCC.

Some genes in our model have been proven to be related to tumor prognosis and immune regulation. For example, AHNAK can promote tumor metastasis via transforming growth factor-β-mediated epithelial-mesenchymal transition [29], and CCL21 can promote immune activity in the TME through co-localization of dendritic cells (DC) and play a role in T-cell programming of ectopic lymph node architectural structures, which has been correlated to cancer prognosis [30]. CD48 is indispensable as an environmental sensor for regulating HSC and progenitor cell numbers, and it also inhibits tumor development [31]. Differential activation of the transcription factor IRF1 laid the foundation for the distinct immune responses elicited by Type I and Type III interferons [32]. Levels of Januse kinase 1 mRNA have been correlated with prognosis and degree of immune infiltration of breast cancer [33]. Matrix metalloproteinase-12 (MMP-12) over-presentation is associated with poor prognosis in HCC [34]. Stimulating MSR1 can enhance c-Jun N-terminal kinase -mediated inflammation in interleukin-4-activated macrophages [35]. NEDD4 exerts a key effect in enhancing HCC multiplication and metastasis by activating the PTEN/phosphatidylinositol-3 kinase/Akt signaling pathway [36]. SEMA3F can accelerate HCC metastasis by activating the focal adhesion pathway [37].CCL5 is a mediator for CD40-driven CD4+ T-cell tumor infiltration as well as immunity in pancreatic cancer.

Conclusions

Our study provides a new way to evaluate prognosis in patients with HCC, but there are still some limitations. Although our prognostic model has been verified internally and externally, this was a retrospective study, so further prospective clinical trials are necessary. The relative expression level exhibited by the two immune genes had a significant impact on prognosis and further study of its mechanism is warranted.

This research suggests a novel approach for evaluating the prognosis of HCC, which may be beneficial for guiding individualized treatment of patients.