Identification of a Five Immune Term Signature for Prognosis and Therapy Options (Immunotherapy versus Targeted Therapy) for Patients with Hepatocellular Carcinoma

Bin, Xiaoyun; Luo, Zongjiang; Wang, Jianchu; Zhou, Sufang

doi:https://doi.org/10.1155/2023/8958962

Computational and Mathematical Methods in Medicine

On this page

Abstract Introduction Methods Results Discussion Data Availability Conflicts of Interest Supplementary Materials References Copyright Related Articles

Special Issue

Machine Learning Approaches Based on Multiscale Data in Diagnosis, Treatment and Prognosis of Diseases

View this Special Issue

Research Article | Open Access

Volume 2023 | Article ID 8958962 | https://doi.org/10.1155/2023/8958962

Identification of a Five Immune Term Signature for Prognosis and Therapy Options (Immunotherapy versus Targeted Therapy) for Patients with Hepatocellular Carcinoma

Xiaoyun Bin,¹Zongjiang Luo,²Jianchu Wang,²and Sufang Zhou^1,3

Academic Editor: Xiaohan Ren

Received25 Aug 2022

Revised15 Oct 2022

Accepted17 Oct 2022

Published02 Feb 2023

Abstract

Background. Immune microenvironment implicated in liver cancer development. Nevertheless, previous studies have not fully investigated the immune microenvironment in liver cancer. Methods. The open-access data used for analysis were obtained from The Cancer Genome Atlas (TCGA-LIHC) and the International Cancer Genome Consortium databases (ICGC-JP and ICGC-FR). R program was employed to analyze all the data statistically. Results. First, the TCGA-LIHC, ICGC-FR, and ICGC-JP cohorts were selected for our analysis, which were merged into a combined cohort. Then, we quantified 53 immune terms in this combined cohort with large populations using the ssGSEA algorithm. Next, a prognostic approach was established based on five immune principles (CORE.SERUM.RESPONSE.UP, angiogenesis, CD8.T.cells, Th2.cells, and B.cells) was established, which showed great prognostic prediction efficiency. Clinical correlation analysis demonstrated that high-risk patients could reveal higher progressive clinical features. Next, to examine the inherent biological variations in high- and low-risk patients, pathway enrichment tests were conducted. DNA repair, E2F targets, G2M checkpoints, HEDGEHOG signaling, mTORC1 signaling, and MYC target were positively correlated with the risk score. Examination of genomic instability revealed that high-risk patients may exhibit a higher tumor mutation burden score. Meanwhile, the risk score showed a strong positive correlation with the tumor stemness index. In addition, the Tumor Immune Dysfunction and Exclusion outcome indicated that high-risk patients could be higher responsive to immunotherapy, whereas low-risk patients may be higher responsive to Erlotinib. Finally, six characteristic genes DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11, and TRIP13 were identified for risk group prediction. Conclusions. In summary, our study identified a signature as a useful tool to indicate prognosis and therapy options for liver cancer patients.

1. Introduction

Every year, around 840, 000 new cases of liver cancer are diagnosed and 780, 000 people lose their lives as a result of this disease. Hepatocellular carcinoma is the most common pathological subtype among all cases, accounting for 90% of primary liver cancer [1]. As a multifactorial and multicause disease, the genesis and progression of liver cancer are linked to numerous risk variables, such as genetics, lifestyle, and environmental aspects [2]. At present, surgery is still the first-line therapy option for liver cancer patients in the early stage [3]. However, considering its insidious symptoms, a substantial portion of patients have been in an advanced stage when first diagnosed, therefore, leading to loss of the appropriate timing of surgery [4]. Meanwhile, due to the characteristic of rapid progression and early metastasis, liver cancer patients tend to have a poor prognosis [2]. Therefore, it is important to find novel molecules and directions with implications for liver cancer diagnosis and therapy.

The tumor microenvironment (TME) significantly affects the biological process in cancer progression [5]. Immune cells and status are the essential components in TME. Different factors, such as cytokines, chemokines, and others, can help cancer cells shape their microenvironment to support their growth [6]. Meanwhile, the reprogramming of other cells surrounding cancer cells plays a decisive role in tumor survival and progression [7]. An example is that TME can enhance immunosuppressive M2 monocyte-derived macrophages by secreting cytokines such as IL-4, which allows the tumor to grow and progress because monocyte-derived macrophages can account for 50% of the tumor mass [8]. Also, through a CXCL13/CXCR5/NFB/p65/miR-934 positive feedback mechanism, Zhao et al. demonstrated that tumor-derived exosomal miR-934 might stimulate macrophage M2 polarizing action to increase liver metastasis of colorectal malignancy [9]. The crosstalk between immune factors and other systems can make TME more complex [10]. In TME, cancer cells showed different metabolic modes from normal cells, like the “Warburg” effect, which has a broader significance in regulating tumor immunity [10]. This energetic interaction involving tumor and immune cells results in metabolic conflict in the tumor ecosystem, restricts the number of nutrients, and causes microenvironmental acidity, which inhibits the activity of immune cells ¹⁰. In addition, metabolic reprogramming is essential for the maintenance of immune cell stability and body balance [11]. Presently, increasing studies have pointed out that metabolic reprogramming occurs in the mechanism of immune cell growth, development, and functional activity, which is crucial for immune reaction [12]. Nowadays, immunotherapy presents promising therapeutic effects in specific cancer populations. Cancer immune status can also influence the intensity and time of the anticancer reaction of immunotherapy [13]. Therefore, it is meaningful to explore the immune microenvironment in liver cancer.

Here, we firstly quantified 53 immune terms in a combined cohort with large populations using the ssGSEA algorithm. Then, a prognostic modeling was created based on five immunological elements (CORE.SERUM.RESPONSE.UP, angiogenesis, CD8.T.cells, Th2.T.cells, and B.cells), which showed great prognosis prediction efficiency. Analysis of clinical correlations revealed that high-risk patients may exhibit higher clinical characteristics progression. Subsequently, a pathway enrichment analysis was conducted to investigate intrinsic biological variations between high- and low-risk patients. Assessment of genomic instability revealed that high-risk individuals could possess a greater TMB score. The risk score showed a high positive correlation with tumor stemness index. In addition, the Tumor Immune Dysfunction and Exclusion (TIDE) outcome indicated that high-risk patients may show higher responsiveness to immunotherapy, whereas low-risk patients may show higher responsiveness to Erlotinib. Finally, six characteristic genes DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11, and TRIP13 were identified for risk group prediction.

2. Methods

2.1. Data Acquisition

Open access data retrieval was carried out based on The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/), Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/gds), and International Cancer Genome Consortium (ICGC; https://dcc.icgc.org/) databases (Access time: 2022/06/19). For the TCGA database, the expression profile information was obtained in a manner of “STAR-Counts” and then collated using the author’s code. Clinical information was obtained in the “xml” form. For the GEO database, GSE14520 and GSE76427 were first included in the study for complete prognosis and transcriptional profiling data. After data quality evaluation, the GSE14520 was eliminated for the reason that the number of probes is less than 20,000; the GSE76427 was eliminated for a large number of missing values (NaN). For the ICGC, the ICGC-JP and ICGC-FR projects were included in the study. Sva package was applied for the data combination with batch effect reduction. First, by taking the logarithm, the order of magnitude of the data (TCGA and ICGC) reached the same range. Then, the combat function in sva package was used to detect and reduce the batch effect between different cohorts. Baseline information of the included patients was shown in Tables 1–3.

2.2. Immune Term Quantification

Immune elements were quantified using the single sample gene set enrichment analysis (ssGSEA) technique, which is an in-built algorithm of Gene Set Variation Analysis (GSVA) [14]. The reference immune terms set was obtained from the previous study, which was used to quantify the enrichment score of 53 immune terms [15]. The advantage of ssGSEA algorithm is the high freedom, in which you can quantify the enrichment score according to the given gene set. However, considering that gene sets can be freely defined, potential quality bias is inevitable.

2.3. Prognosis Analysis and Model Construction

Based on the assessment of 53 immunological elements from the ssGSEA algorithm, a univariate Cox regression analysis was conducted to discover the prognosis-related terms (). Afterward, the random survival forest variable hunting (RSFVH) method was processed to reduce the number of dimensions and filter genes. Finally, a multivariate Cox regression analysis was conducted for the development of a prognostic model. Kaplan-Meier (KM) analysis was employed to detect the best gene combination or final signature by analyzing log-rank values. Each enrolled patient with complete prognosis and expression profile data was assigned a risk score with the formula of “” [15]. If the risk score was greater than the median value, the patients were classified into the high-risk or low-risk group, accordingly. ROC and KM survival curves were utilized to assess the accuracy of our model’s prognostic projections. Univariate and multivariate Cox regression were also employed to validate the independence of the prediction model. The 1-, 3-, and 5-year survival can reflect the short-, medium-, and long-term prognosis of patients, and therefore, were selected as the time node in prognosis analysis.

2.4. Pathway Enrichment and Genomic Analysis

The GSVA and GSEA algorithms were utilized for pathway enrichment analysis [16]. Hallmark was used as the standard gene set for the GSVA algorithms, whereas metabolism-related gene sets (41 metabolism terms) were acquired from the website https://www.gsea-msigdb.org/. Standard gene sets for the GSEA algorithm were c2.cp.kegg.v7.5.1.symbols.gmt and c5.go.v7.5.1.symbols.gmt. The TCGA database was accessed using genomic mutation information, including the tumor mutation burden (TMB) and microsatellite instability (MSI) score. Based on the expression profile and utilizing the one-class logistic regression machine learning (OCLR) machine-learning technique, the tumor stemness index was computed [17].

2.5. Immunotherapy and Drug Sensitivity Analysis

Using the TIDE methodology [18], patients were evaluated for immunotherapy sensitivity. The parameter of “Cancer type” was set as “Other.” The parameter of “Previous immunotherapy” was set as “No.” The analysis of drug responsiveness was conducted using the database of Genomics of Drug Responsiveness in Cancer [19].

2.6. Feature Gene Identification

To identify the feature genes for the risk group, the LASSO regression and SVM-RFE (support vector machine recursive feature elimination) algorithm were applied to find the best variable [20].

2.7. Western Blot

Total proteins were extracted using a total protein extraction kit (Beyotime, China). Western blot was conducted based on the standardized process (10% SDS-PAGE gel). The primary antibody of CALCRL (1 : 2000) and GAPDH (1 : 50000) was purchased from Proteintech.

2.8. Statistical Analysis

This research was analyzed using R software version 4.2.1. Two-sided values < 0.05 were considered statistically significant. For continuous variables with normal distributions, an independent -test was applied, and for continuous variables with skewed distributions, a Wilcoxon rank-sum test was conducted. The study of differentially expressed genes (DEGs) was conducted using the limma program with the criteria and .

3. Results

3.1. Immune Term Quantification

Figure 1 displays the flowchart of the entire investigation. Three distinct liver cancer cohorts, TCGA-LIHC, ICGC-JP, and ICGC-FR, were chosen for our research (Figure 2(a)). The sva package was employed to combine data and decrease the batch effect. Then, a significant batch impact decline was noted (Figure 2(b)). The ssGSEA method was employed in the pooled cohorts to quantify 53 immunological elements (Figure 2(c)).

(a)

(b)

(c)

3.2. Prognosis Model Construction

First, the patients were randomly divided into training and validation cohorts according to the 1 : 1 ratio. Based on the 53 immunological elements, a univariate Cox regression analysis was conducted to find the prognosis-related variables with (Figure 3(a) and Table 4). The random forest approach was then used to reduce the dimensionality, and the top ten important terms were CORE.SERUM.RESPONSE.UP, angiogenesis, CD8.T.cells, Cytotoxic.cells, CSR.activated, Th2.cells, IL13.score, TcClassII.score, B.cells, and T.cells.receptors.score (Figures 3(b) and 3(c)). Through multivariate Cox regression analysis and permutations, the five immune terms were used for prognosis model construction, including CORE.SERUM.RESPONSE.UP, angiogenesis, CD8.T.cells, Th2.cells, and B.cells (Figure 3(d)). The risk score was calculated with the formula of “Risk score = CORE.SERUM.RESPONSE.UP 0.281 + angiogenesis -0.214 + CD8.T.cells -0.071 + Th2.cells 0.031 + B.cells -0.027”. The overview of our model (training group) was shown in Figure 3(e), whereas a greater number of fatalities were seen in the high-risk group. The KM survival curve revealed that high-risk individuals tend to have a poorer outcome (Figure 3(f), , , and ). ROC curves for 1-, 3-, and 5-year patients demonstrated a high predictive accuracy (Figures 3(g)–3(i), 1-year: , 3-year: , and 5-year: ). Also, in the validation group, the same trend was observed (Figure 3(j)). KM survival curves indicated that the high-risk patients might have a worse prognosis performance (Figure 3(k), , , and ). Meanwhile, the performances of ROC curves for 1-, 3-, and 5-year-old patients are still satisfactory (Figures 3(l)–3(n), 1-year: , 3-year: , and 5-year: ).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

Figure 3

Prognosis model construction. Notes: (a) univariate Cox regression analysis was utilized to find prognosis-related elements with a significance level of ; (b, c) the random forest technique was used for dimensionality reduction; (d) five elements, including CORE.SERUM.RESPONSE.UP, angiogenesis, CD8.T.cells, Th2.cells, and B.cells, were selected for model creation; (e) the summary of our model in the training cohort; (f) the KM curves indicated that high-risk patients may have a poorer outcome (training cohort); (g–i) ROC curves for one, three, and five years (training cohort); (j) the summary of our model in the validation cohorts; (k) the KM curves indicated that high-risk patients may have a poorer outcome (validation cohort); and (l–n) ROC curves for one, three, and five years (validation cohort).

3.3. Clinical Correlation Analysis

Furthermore, we explore the clinical correlation of our model. Univariate Cox regression and multivariate Cox regression analysis revealed because our model is not related to other clinical characteristics (Figures 4(a) and 4(b), univariate: , ; multivariate: , ). Clinical correlation analysis showed that high-risk patients could possess more aggressive clinical characteristics, such as clinical stage, grade, and T classification (Figure 4(c)). Interestingly, we found that the patients with mild adjacent hepatic tissue inflammation might have a higher risk score compared to the severe group (Figure 4(d)). Moreover, the patients with had a higher risk score than those with (Figure 4(e)). No significant difference was found in different gender patients. Asian populations might have a higher risk score than White populations (Figure 4(f)). A negative correlation was found between height, weight, and BMI (Figures 4(h)–4(j), weight, , , height, , , BMI, , ).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

3.4. Pathway Enrichment Analysis

We investigated the biological differences between high-and low-risk patients. For Hallmark pathways, we found that risk score was positively correlated with DNA repair, E2F targets, G2M checkpoints, HEDGEHOG signaling, mTORC1 signaling, and MYC targets (Figure 5(a)). For metabolism-related pathways, we observed that riskscore was positively correlated with purine and pyrimidine metabolism (Figure 5(a)). The GSEA method illustrated that in the high-risk group, the terms of DNA-dependent DNA replication, chromosome, cell-cell junction assembly, nuclear chromosome, leukocyte transendothelial migration, ubiquitin-mediated proteolysis, tight junction, actin cytoskeleton regulation, and MAPK signaling were upregulated (Figures 5(b) and 5(c)). Also, based on the DEGs discovered between high- and low-risk patients, a network of protein-protein linkages was developed (Figure 5(d)). ClueGO analysis revealed that these nodes were mainly enriched in fat-soluble vitamin catabolic process, proximal/distal pattern formation, regulation of neuronal synaptic plasticity, and membrane depolarization (Figure 5(e)).

(a)

(b)

(c)

(d)

(e)

3.5. Genomic Instability Analysis

The genomic feature can also affect the tumor biological process. Genetic variations are between high- and low-risk patients. Mutation information was obtained from the TCGA database (Figure 6(a)). The outcome indicated that high-risk patients may possess a greater TMB score (Figure 6(b)). Meanwhile, we found that all mutant counts, synonymous mutation counts, and nonsynonymous mutation counts were elevated in patients at high risk (Figure 6(c)). No significant variation in the MSI score was observed between high- and low-risk patients (Figure 6(d)). Moreover, we discovered that TP53 was the gene greatest substantially mutated across patients at high and low risk (Figure 6(e)). The KM curve indicated that patients with TP53 mutations may suffer a poorer outcome (Figure 6(f)). In addition, a positive significant correlation was identified between risk score and tumor stemness index (Figures 6(g) and 6(h), mDNAsi, , ; mRNAsi, , ).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 6

Genomic instability analysis. (a) A summary of the distribution of TMB in TCGA pan-cancer; (b) comparison of TMB between high-risk and low-risk patients; (c) the disparity between high-risk and low-risk patients’ mutant count, nonsynonymous mutation count, and synonymous mutation count; (d) the MSI distinction between patients at high and low risk; (e) TP53 and MUC16 were the most frequently mutated genes comparing patients at high and low risk; (f) KM survival curve of the TP53 wild and mut patients; and (g, h) the correlation of riskscore with tumor stemness index.

3.6. Immunotherapy and Drug Sensitivity Analysis

Immune checkpoint modules contribute substantially to the course of cancer; hence, we analyzed the correlation between the risk score and numerous checkpoint modules. Correlation analysis indicated significant variations between the high- and low-risk groups for various immune checkpoint modules, such as PD-L1 and CTLA-4 (Figures 7(a)–7(e)). TIDE research revealed that high-risk patients may show a lower TIDE score and a greater probability of immunotherapy responses (Figures 7(f) and 7(g)). The submap algorithm suggested that high-risk patients could be highly responsive to PD-1 and CTLA-4 therapy (Figure 7(h)). Analysis of drug responsiveness revealed that low-risk patients could be higher responsive to Erlotinib (Figure 7(i)).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 7

The risk score was associated with the sensibility to immunotherapy and chemotherapy. (a) Several disparities in immunological checkpoints between high- and low-risk groups; (b–e) variations in PD-L1, CTLA4, PD-1, and PD-L2 between high- and low-risk patients; (f, g) High-risk patients may have a lower TIDE score and a higher percentage of immunotherapy responders; (h) high-risk individuals may respond better to CTLA4 and PD-1 treatment; and (i) low-risk patients may respond better to Erlotinib.

3.7. Identification of the Characteristic Genes of the Risk Group

Considering the prognosis and therapy sensitivity difference between high- and low-risk patients, we try to identify the characteristic gene that could robustly indicate the risk group. Utilizing LASSO regression and the SVM-RFE technique, distinctive genes were identified (Figures 8(a) and 8(b)). The intersection of these two algorithms identified six genes, including DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11, and TRIP13 (Figure 8(c) LASSO logistic regression: DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11, and TRIP13; SVM-RFE: CALCRL, MCM10, SCN4A, NGFR, DEPDC1, PRR11, DEPDC1B, RAD51, ANLN, PYGM, CD5L, DCN, IYD, GLI2, TRIP13, TNMD, and SLC12A1). Among these genes, DEPDC1, DEPDC1B, PRR11, and TRIP13 were increased in the high-risk group, while NGFR and CALCRL were decreased (Figure 8(d)). All these characteristic genes showed great prediction efficiency in the patients’ risk group (Figure 8(e), training cohort, DEPDC1, ; DEPDC1B, ; TRIP13, ; CALCRL, ; PRR11, ; NGFR, ). Logistic regression was performed to combine these genes with the formula of “”, which showed extremely great prediction efficiency in patients risk group (Figure 8(f), ). In the validation cohort, these characteristic genes also showed great prediction efficiency in the patients’ risk group (Figure S1A–F, DEPDC1, ; DEPDC1B, ; TRIP13, ; CALCRL, ; PRR11, ; NGFR, ), as well as the logistic regression model (Figure S1G, ). Among these characteristic genes, DEPDC1, DEPDC1B, PRR11, and TRIP13 were risk factors for liver cancer patients (Figure 8(g)). Meanwhile, DEPDC1, DEPDC1B, CALCRL, PRR11, and TRIP13 were significantly upregulated in liver cancer tissue, yet NGFR was downregulated (Figure 8(h)). In liver cancer, the genes DEPDC1, DEPDC1B, PRR11, and TRIP13 have been explored in previous studies [21–24]. Therefore, we selected CALCRL for further validation. Western blot indicated that the protein level of CALCRL was upregulated in liver cancer tumor tissue (Figure S2).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8

Identification of the characteristic genes of the risk group. Notes: (a) the LASSO regression algorithm was employed to determine the distinctive genes; (b) the SVM-RFE algorithm was utilized to determine the genes features; (c, d) DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11, and TRIP13 were identified as the characteristic genes; (e) ROC curves were utilized to assess the prediction accuracy of characteristic genes on risk group (training cohort); (f) the logistic regression model showed an extremely high prediction efficiency in the risk group (training cohort); (g) the prognosis effect of DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11, and TRIP13; and (h) expression level of these six genes in tumor and normal tissue.

4. Discussion

Liver cancer is still a serious public health concern worldwide [25]. Recently, researchers have focused on the immune microenvironment in liver cancer, which can significantly affect the progression of the disease. Thus, a comprehensive investigation of the liver cancer immune microenvironment can contribute to its diagnosis and therapy options.

In our study, we firstly quantified 53 immune terms in a combined cohort with large populations using the ssGSEA algorithm. Afterward, a prognostic model based on five immune elements (CORE.SERUM.RESPONSE.UP, angiogenesis, CD8.T.cells, Th2.cells, and B.cells) was established, which showed great prognosis prediction efficiency in both training and validation cohorts. Analysis of clinical correlations revealed that high-risk patients could possess higher clinical progression characteristics. Angiogenesis plays an important role in tumor metastasis and progression. From the comprehensive review conducted by Morse et al., abnormal angiogenesis in liver cancer often leads to poor prognosis and facilitated progression [26]. Wolf et al. found that the cross-talk between intrahepatic CD8+ T cells and NKT cells contributes to nonalcoholic steatohepatitis and liver cancer [27]. Xu et al. found that the Th2 response could affect liver fibrosis, which is a risk factor for liver cancer [28]. Garnelo et al. indicated that the interaction between tumor-infiltrating B cells and T cells significantly affects the progression of liver cancer [29].

Next, a pathway enrichment analysis was conducted to investigate the intrinsic biological variations between high-and low-risk patients. Analysis of genomic instability revealed that high-risk patients could possess a higher TMB score. The risk score showed a high positive correlation with the tumor stemness index. In addition, the TIDE outcome indicated that high-risk patients could be better responsive to immunotherapy, whereas low-risk patients could be better responsive to Erlotinib. Finally, six characteristic genes DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11, and TRIP13 were identified for risk group prediction.

Pathway enrichment analysis showed that in high-risk patients, the Hallmark pathway of DNA repair, E2F targets, G2M checkpoints, Hedgehog signaling, mTORC1 signaling, and MYC targets were elevated. Malignant tumor cells have higher genomic damage compared to normal cells. Meanwhile, abnormal DNA repair procedures are difficult to meet the needs of DNA damage, therefore, leading to a more aggressive phenotype [30]. G2/M checkpoint is a critical phase in the cell cycle, which might directly affect cell proliferation [31]. Hedgehog signaling has been reported and extensively included in liver cancer progression. For instance, Gu et al. revealed that circular RNA circIPO11 could drive self-renewal and cancer progression of liver cancer through Hedgehog signaling, which might be an underlying therapeutic target [32]. In addition, Wu et al. discovered that CK2 might enhance stemness and chemotherapy resistance via the Hedgehog signaling pathway in liver cancer [33]. These findings suggested that patients in the high-risk category may have a higher level of activity in the aforementioned pathways, thus resulting in a poor prognosis.

Genomic analysis revealed that high-risk patients might have a greater level of genomic instability. In cancer, genomic instability is a prominent feature, which is also a key marker for separating cancerous cells from normal ones [34]. Genomic instability is defined as the increase in mutation frequency in the genome and its potential source is the sum of defects in the DNA damage and repair pathway, which could lead to uncontrolled proliferation of cancer cells [34]. In the pathway enrichment analysis, we have found the abnormal activity of DNA repair in high-risk patients. Therefore, higher genomic instability in high-risk patients seems to be reasonable. Moreover, a strong correlation was observed between the risk score and tumor stemness. Liver cancer has a high recurrence rate, making it one of the most highly deadly malignant tumors, partly due to the existence of cancer stem cells (CSC) [35]. Tumor stemness could also affect the immune microenvironment of liver cancer. For example, M2 TAMs could promote the expression of stemness proteins [36]. Another aspect, the CXCL11, a cytokine involved in the recruitment of activated T cells to inflammatory sites, was found significantly affect the stemness genes, sphere formation, and tumorigenicity in liver cancer [37].

Immunotherapy and chemotherapy are the best choices for the advanced liver cancer stage [4]. Low-risk patients could be higher responsive to Erlotinib, whereas high-risk patients could be more responsive to PD-1 and CTLA-4 treatment. According to the LASSO regression and SVM-RFE algorithm, we identified six characteristic genes, including DEPDC1, DEPDC1B, NGFR, CALCRL, PRR11 and TRIP13. The logistic regression model showed an extremely great prediction efficiency in the patients risk group. In the clinical application, measuring the expression level of these six genes might accurately identify the risk group of a patient, which might have the potential to guide the prognosis and therapy options of liver cancer patients.

Although our study was performed based on high-quality bioinformatics analysis, some limitations should be noticed. First, the populations included in our study were mainly White populations. Therefore, the underlying race bias might reduce the credibility of our conclusion. Second, some specific clinical features, for example, laboratory examination, lifestyle, and other basic diseases were not provided in most patients. However, our approach proved an excellent pattern for identifying liver cancer patients’ survival time as well as their immunotherapy and chemotherapy responsiveness.

Data Availability

Based on reasonable request, data are available from the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Supplementary Materials

Supplementary 1. Figure S1: the prediction accuracy of characteristic genes on risk group in validation cohort. Notes: (a–g): ROC curves of characteristic genes.

Supplementary 2. Figure S2: Western blot of CALCRL.

References

F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, and A. Jemal, “Global cancer statistics 2018: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries,” CA: a Cancer Journal for Clinicians, vol. 68, no. 6, pp. 394–424, 2018.
View at: Publisher Site | Google Scholar
A. Marengo, C. Rosso, and E. Bugianesi, “Liver cancer: connections with obesity, fatty liver, and cirrhosis,” Annual Review of Medicine, vol. 67, no. 1, pp. 103–117, 2016.
View at: Publisher Site | Google Scholar
A. Kanazawa, T. Tsukamoto, S. Shimizu, S. Yamamoto, A. Murata, and S. Kubo, “Laparoscopic hepatectomy for liver cancer,” Digestive Diseases, vol. 33, no. 5, pp. 691–698, 2015.
View at: Publisher Site | Google Scholar
D. Anwanwan, S. K. Singh, S. Singh, V. Saikam, and R. Singh, “Challenges in liver cancer and possible treatment approaches,” Biochimica et Biophysica Acta Reviews on Cancer, vol. 1873, no. 1, article 188314, 2020.
View at: Publisher Site | Google Scholar
G. Wang, Q. Wang, N. Liang et al., “Oncogenic driver genes and tumor microenvironment determine the type of liver cancer,” Cell Death & Disease, vol. 11, no. 5, p. 313, 2020.
View at: Publisher Site | Google Scholar
X. Mao, J. Xu, W. Wang et al., “Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives,” Molecular Cancer, vol. 20, no. 1, p. 131, 2021.
View at: Publisher Site | Google Scholar
P. McGuirk and K. H. Mills, “Pathogen-specific regulatory T cells provoke a shift in the Th1/Th2 paradigm in immunity to infectious diseases,” Trends in Immunology, vol. 23, no. 9, pp. 450–455, 2002.
View at: Publisher Site | Google Scholar
J. Kim and J. S. Bae, “Tumor-associated macrophages and neutrophils in tumor microenvironment,” Mediators of Inflammation, vol. 2016, Article ID 6058147, 11 pages, 2016.
View at: Publisher Site | Google Scholar
S. Zhao, Y. Mi, B. Guan et al., “Tumor-derived exosomal mir-934 induces macrophage M2 polarization to promote liver metastasis of colorectal cancer,” Journal of Hematology & Oncology, vol. 13, no. 1, p. 156, 2020.
View at: Publisher Site | Google Scholar
L. Xia, L. Oyang, J. Lin et al., “The cancer metabolic reprogramming and immune response,” Molecular Cancer, vol. 20, no. 1, p. 28, 2021.
View at: Publisher Site | Google Scholar
H. J. Hurley, H. Dewald, Z. S. Rothkopf et al., “Frontline science: Ampk regulates metabolic reprogramming necessary for interferon production in human plasmacytoid dendritic cells,” Journal of Leukocyte Biology, vol. 109, no. 2, pp. 299–308, 2021.
View at: Publisher Site | Google Scholar
S. Ricciardi, N. Manfrini, R. Alfieri et al., “The translational machinery of human CD4⁺ T cells is poised for activation and controls the switch from quiescence to metabolic remodeling,” Cell Metabolism, vol. 28, no. 6, pp. 895–906.e5, 2018.
View at: Publisher Site | Google Scholar
D. S. Chen and I. Mellman, “Elements of cancer immunity and the cancer-immune set point,” Nature, vol. 541, no. 7637, pp. 321–330, 2017.
View at: Publisher Site | Google Scholar
S. Hänzelmann, R. Castelo, and J. Guinney, “GSVA: gene set variation analysis for microarray and Rna-Seq data,” BMC Bioinformatics, vol. 14, no. 1, p. 7, 2013.
View at: Publisher Site | Google Scholar
X. Ren, X. Chen, X. Zhang et al., “Immune microenvironment and response in prostate cancer using large population cohorts,” Frontiers in Immunology, vol. 12, article 686809, 2021.
View at: Publisher Site | Google Scholar
A. Subramanian, P. Tamayo, V. K. Mootha et al., “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 43, pp. 15545–15550, 2005.
View at: Publisher Site | Google Scholar
T. M. Malta, A. Sokolov, A. J. Gentles et al., “Machine learning identifies stemness features associated with oncogenic dedifferentiation,” Cell, vol. 173, no. 2, pp. 338–54.e15, 2018.
View at: Publisher Site | Google Scholar
J. Fu, K. Li, W. Zhang et al., “Large-scale public data reuse to model immunotherapy response and resistance,” Genome Medicine, vol. 12, no. 1, p. 21, 2020.
View at: Publisher Site | Google Scholar
W. Yang, J. Soares, P. Greninger et al., “Genomics of drug sensitivity in cancer (Gdsc): a resource for therapeutic biomarker discovery in cancer cells,” Nucleic Acids Research, vol. 41, pp. D955–D961, 2013.
View at: Publisher Site | Google Scholar
R. Y. Choi, A. S. Coyner, J. Kalpathy-Cramer, M. F. Chiang, and J. P. Campbell, “Introduction to machine learning, neural networks, and deep learning,” Translational Vision Science & Technology, vol. 9, no. 2, p. 14, 2020.
View at: Publisher Site | Google Scholar
C. Tian, M. Abudoureyimu, X. Lin, X. Chu, and R. Wang, “Linc-Ror facilitates progression and angiogenesis of hepatocellular carcinoma by modulating Depdc1 expression,” Cell Death & Disease, vol. 12, no. 11, p. 1047, 2021.
View at: Publisher Site | Google Scholar
X. W. Dang, Q. Pan, Z. H. Lin et al., “Overexpressed Depdc1b contributes to the progression of hepatocellular carcinoma by Cdk1,” Aging, vol. 13, no. 16, pp. 20094–20115, 2021.
View at: Publisher Site | Google Scholar
W. Qiao, H. Wang, X. Zhang, and K. Luo, “Proline-rich protein 11 silencing inhibits hepatocellular carcinoma growth and epithelial-mesenchymal transition through Β-catenin signaling,” Gene, vol. 681, pp. 7–14, 2019.
View at: Publisher Site | Google Scholar
M. X. Zhu, C. Y. Wei, P. F. Zhang et al., “Elevated Trip13 drives the Akt/Mtor pathway to induce the progression of hepatocellular carcinoma via interacting with Actn4,” Journal of Experimental & Clinical Cancer Research: CR, vol. 38, no. 1, p. 409, 2019.
View at: Publisher Site | Google Scholar
C. Ayuso, J. Rimola, R. Vilana et al., “Diagnosis and staging of hepatocellular carcinoma (Hcc): current guidelines,” European Journal of Radiology, vol. 101, pp. 72–81, 2018.
View at: Publisher Site | Google Scholar
M. A. Morse, W. Sun, R. Kim et al., “The role of angiogenesis in hepatocellular carcinoma,” Clinical Cancer Research: An Official Journal of the American Association for Cancer Research, vol. 25, no. 3, pp. 912–920, 2019.
View at: Publisher Site | Google Scholar
M. J. Wolf, A. Adili, K. Piotrowitz et al., “Metabolic activation of intrahepatic Cd8+ T cells and Nkt cells causes nonalcoholic steatohepatitis and liver cancer via cross-talk with hepatocytes,” Cancer Cell, vol. 26, no. 4, pp. 549–564, 2014.
View at: Publisher Site | Google Scholar
F. Xu, R. Cheng, S. Miao et al., “Prior Toxoplasma gondii infection ameliorates liver fibrosis induced by Schistosoma japonicum through inhibiting Th2 response and improving balance of intestinal flora in mice,” International Journal of Molecular Sciences, vol. 21, no. 8, 2020.
View at: Publisher Site | Google Scholar
M. Garnelo, A. Tan, Z. Her et al., “Interaction between tumour-infiltrating B cells and T cells controls the progression of hepatocellular carcinoma,” Gut, vol. 66, no. 2, pp. 342–351, 2017.
View at: Publisher Site | Google Scholar
A. K. Basu, “DNA damage, mutagenesis and cancer,” International Journal of Molecular Sciences, vol. 19, no. 4, 2018.
View at: Publisher Site | Google Scholar
L. L. Liu, J. M. Zhu, X. N. Yu et al., “Ube2t promotes proliferation via G2/M checkpoint in hepatocellular carcinoma,” Cancer Management and Research, vol. 11, pp. 8359–8370, 2019.
View at: Publisher Site | Google Scholar
Y. Gu, Y. Wang, L. He et al., “Circular Rna Circipo11 drives self-renewal of liver cancer initiating cells via Hedgehog signaling,” Molecular Cancer, vol. 20, no. 1, p. 132, 2021.
View at: Publisher Site | Google Scholar
D. Wu, Y. Q. Yin, Y. Li, L. Zhang, Y. H. Jiang, and Z. Wang, “Ck2α causes stemness and chemotherapy resistance in liver cancer through the Hedgehog signaling pathway,” Hepatobiliary & Pancreatic Diseases International, vol. 20, 2021.
View at: Publisher Site | Google Scholar
H. Li, S. E. Zimmerman, and U. Weyemi, “Genomic instability and metabolism in cancer,” International Review of Cell and Molecular Biology, vol. 364, pp. 241–265, 2021.
View at: Publisher Site | Google Scholar
M. Huch and L. Dollé, “The plastic cellular states of liver cells: are Epcam and Lgr5 fit for purpose?” Hepatology, vol. 64, no. 2, pp. 652–662, 2016.
View at: Publisher Site | Google Scholar
F. Zhu, X. Li, S. Chen, Q. Zeng, Y. Zhao, and F. Luo, “Tumor-associated macrophage or chemokine ligand Ccl17 positively regulates the tumorigenesis of hepatocellular carcinoma,” Medical Oncology, vol. 33, no. 2, p. 17, 2016.
View at: Publisher Site | Google Scholar
Y. Zhang, W. Zhao, S. Li et al., “CXCL11 promotes self-renewal and tumorigenicity of α2δ1⁺ liver tumor-initiating cells through CXCR3/ERK1/2 signaling,” Cancer Letters, vol. 449, pp. 163–171, 2019.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2023 Xiaoyun Bin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

333

Downloads

436

Citations