A Novel Hypoxic-Angiogenesis-Immune-Related Gene Model for Prognostic and Therapeutic Effect Prediction in Hepatocellular Carcinoma Patients

Background Hepatocellular carcinoma (HCC) is one of the most heterogeneous malignant tumors that have been discovered so far, which makes the prognostic prediction difficult. The hypoxia, angiogenesis, and immunity-related genes (HAIRGs) are closely related to the development of liver cancer. However, the prognostic and treatment effect of hypoxia, angiogenesis, and immunity-related genes in HCC continues to be further clarified. Methods The gene expression quantification data and clinical information in patients with liver cancer were downloaded from the TCGA database, and HAIRG signature was built by using the least absolute shrinkage and selection operator (LASSO) technique. Patient from the ICGC database validated the model. Then, tumor immune dysfunction and exclusion (TIDE) algorithm was applied to estimate the clinical response to immunotherapy and the sensitivity of drugs was evaluated by the half-maximal inhibitory concentration (IC50). Result The HAIRGs were identified between the HCC patients and normal patients in the TCGA database. In univariate Cox regression analysis, seventeen differentially expressed genes (DEGs) were associated with overall survival (OS). An eight HAIRG signature model was constructed and was used to divide the patients into two groups according to the median value of the risk score base on the TCGA dataset. Patients in the high-risk group had a significant reduction in OS compared to those in the low-risk group (P < 0.001 in the TCGA, P < 0.001 in the ICGC). For TCGA and ICGC databases of univariate Cox regression analyses, the risk score was used as an independent predictor of OS (HR > 1, P < 0.001). Functional analysis showed that the relevant immune pathways and immune responses were enriched, cellular component analysis showed that the immunoglobulin complex and other related substances were enriched, and immune status existed a difference in the high- and low-risk groups. Then, the tumor immune dysfunction and exclusion (TIDE) algorithm presented differences in immune response in the high- and low-risk groups (P < 0.05), and based on drug sensitivity prediction, patients in the high-risk group were more sensitive to cisplatin compared to those in the low-risk group in both the TCGA and ICGC cohorts (P < 0.05). Conclusions HAIRG signature can be utilized for prognostic prediction in HCC, while it can be considered a prediction model for clinical evaluation of immunotherapy response and chemotherapy sensitivity in HCC.


Introduction
Liver cancer is the sixth most commonly diagnosed cancer in terms of morbidity and the fourth leading cause of cancer related to death [1]. According to the WHO prediction, more than one million people will die of liver cancer in 2030, based on an annual projection [1]. Factors such as chronic hepatitis B or C, alcohol addiction, metabolic liver disease (particularly NAFLD), and dietary toxicosis, for instance aflatoxin and aristolochic acid increase the risk of liver cancer development [2]. Being a highly heterogenous form of liver cancer, hepatocellular carcinoma (HCC) generates a low survival rate and increases the difficulty of overall survival prediction for HCC patients. Therefore, developing a precise prognostic and therapeutic model will help to clinical evaluation and prolong the survival time for HCC patients.
Tumor microenvironment (TME) can be affected by some important components, such as hypoxia, angiogenesis, and immune cells, which can influence tumor growth [3]. For the development of HCC, hypoxic stimulation can cause bone marrow cells enter to the TME; then, they differentiate into tumor-associated macrophages or neutrocytes and secrete cytokines and proangiogenic growth factors to promote tumor development [4]. Meanwhile, hypoxia can stimulate the release of hypoxia-inducible factors (HIF), which signal both natural immune cells and HCC cells. It is beneficial to the recruitment and maintenance of prototumor immune cells, inhibition of antitumor immune cells, and promotion of immune escape [5]. Immune escape plays an important role in the occurrence and metastasis of hepatic carcinoma [6]. In a hypoxic environment, some suppressive immune cells, like regulatory T cells and M2 macrophages, are frequently recruited to cancer tissues to form the immunosuppressive microenvironment in HCC, which can secret some procancer inflammatory cytokines and activate the STAT3 and NF-κB signaling pathways [6]. Therefore, hypoxia, angiogenesis, and immune response act as critical roles in the progression of HCC. Nevertheless, it is still unclear whether hypoxia, angiogenesis, and immunity-related genes (HAIRGs) were correlated with the prognosis and immune checkpoint therapy response.
In our study, we first established a multigene prognostic model base on HAIRGs in the TCGA database and validated it in the ICGC database. Finally, we conducted functional enrichment analysis and immune response prediction to explore the underlining mechanism of the prognostic model.

Materials and Methods
2.1. Collecting the Data. Quantitative gene expression data and clinical information of liver cancer patients were downloaded from the TCGA database (https://portal.gdc.cancer; containing 374 liver cancer samples and 50 normal tissue samples) and ICGC database (http://www.ncbi.nlm.nih .gov/geo/; containing 231 liver cancer samples). The genetic information about hypoxia, angiogenesis, and immunity was downloaded from the GeneCards database (https://www .genecards.org/). The TCGA and ICGC databases were used for training queue and validation queue, respectively.

Screening and
Identifying HAIRG Signature Associated with LC Prognosis. The HAIRGs were matched (the top 500 genes of the three gene sets). DEGs between hypoxia, angiogenesis, and immunity-related genes were recognized by "limma" R package with the error discovery rate of < 0.05. The overlapping prognostic DEGs were incorporated into the LASSO Cox regression using the "glmnet" R package. Univariate Cox analysis was accomplished for OS using the "survival" R package to screen HAIRGs with prognostic potential. According to the minimum criteria, the tenfold cross-validation was used in the penalty parameter (λ). Based on the expression of each gene and the corresponding regression coefficient, a risk score was calculated for each patient. The formula: Risk score = SUM ðexpression of each gene × corresponding coefficientÞ. Patients with liver cancer were divided into the high-and low-risk groups according to the median value of the risk score. Principal components analysis (PCA) was carried out by using the "prcomp" function of the "stats" R package based on signature gene expression in the the TCGA database. In addition, t-SNE was implemented through the "Rtsne" R package to investigate the distribution of the two groups.

The Predictive Nomogram Construction and Evaluation.
Univariate and multivariate Cox regression analyses were executed to determine whether the risk score was an independent prognostic predictor for OS compared to other clinical features in the TCGA database. The "rms" R package was utilized to construct a predictive nomogram and corresponding calibration maps based on independent predictive factors. To evaluate the predictive power of the nomogram, a time-dependent receiver operating characteristic (ROC) curve analysis was performed by using the "time ROC" R package. Patients from ICGC were analyzed by using the same formula as that for the TCGA database. The sensitivity and specificity of the nomogram were tested by the ROC curves.

Functional Enrichment Analysis and Immunotherapy
Response Predictions. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the DEGs were analyzed by using the STRING database ð|log 2FC | ≥1, FDR < 0:05) in the high-and low-risk groups. The P value was regulated by the BH method. The single-sample gene set enrichment analysis (ssGSEA) was utilized to calculate the infiltrating score of 16 immune cells and the activity of 13 immune-related pathways [7] in the "gsva" R package.
2.5. Immunotherapy Response Predictions for HCC Patients. TIDE (http://tide.dfci.harvard.edu/) is a combination model of T cell dysfunction and exclusive expression characteristic that can calculate and simulate the tumor immune escape [8]. The TIDE algorithm was utilized for predicting the clinical response to immunotherapy in the TCGA and ICGC cohorts.
2.6. Evaluation of the Sensitivity of Drugs. The Genomics of Drug Sensitivity in Cancer (GDSC; https://www .cancerrxgene.org/) database was used to assess the sensitivity of chemotherapy drugs [9]. The half-maximal inhibitory concentration (IC50) was assessed by using the pRRophetic package in R.

Statistical
Analysis of the Data. All statistical analyses were processed by the R software (version 4.0.3). Student's two-sided t-test was used to compare the expression of genes in HCC tissues and nearby nontumorous tissues. The Kaplan-Meier analysis and log-rank test were used to compare the difference of OS among groups. Univariate and multivariate Cox regression analyses were used to determine independent predictors of OS. SsGSEA score of immune

Disease Markers
pathways or cells were compared between the two groups by Mann-Whitney test. P value < 0.05 were considered statistically significant if no specific requirement is made, and all P values were double-tailed. Since the TCGA, ICGC, and Gen-eCards databases are publicly available, this study strictly followed access policies for databases and publication guidelines; thus, ethical approval from a local ethics committee was not required.

Results
3.1. Identification of Prognostic HAIRGs in the TCGA Cohort. All 44 DEGs between hypoxia, angiogenesis, and immune-related genes were identified (Figure 1   7 Disease Markers CASP8, and LGALS3) were significantly related to the OS for HCC patients, sixteen of which were high hazard ratio genes (HR > 1) and one was protective genes (HR < 1) (Figure 1(b)). Then, we picked out 8 HAIRGs (VEGFA, CTNNB1, PPARG, HSP90AA1, HMOX1, LGALS3, SPP1, and RAC1) from the 17 HAIRG model through the LASSO Cox regression, upregulated genes accounted for 7/8 in tumor tissue, which was shown by a heat map (Figure 1(c)). The correlation between above the 8 HAIRGs was shown in Figure 1(d).

Verifying the HAIRG Signature in the ICGC Database.
To examine the stability of the model established in the TCGA database, we used the same formula as the TCGA database to calculate the risk score for each patient in the ICGC database. Similar to the TCGA database, patients were also classified into high-risk group (n = 115) or low-risk group (n = 116) in the ICGC according to the median value of the risk score ( Figure 3(a)). The results were similar to the TCGA cohort, PC and T-SNE analyses verified the reliable aggregation ability of risk score in ICGC database (Figures 3(c) and 3(d)). It was known that patients in the high-risk group were likely to die earlier and had a shorter lifespan compare to those in the low-risk group (Figures 3(b)-3(e)). Besides, after completing the ROC analysis, AUC values of 1, 2, and 3 years were 0.654, 0.673, and 0.657, respectively (Figure 3(f)).

Establishment and Validation of Nomogram.
Univariable Cox regression analysis showed the higher cancer stage and risk scores were independent risk factors for prognosis (P < 0:05). Based on the multivariate analysis, tumor stage and risk score were also served as independent predictors for OS (P < 0:05). Later, these two variables were utilized to build a nomogram of OS, including the tumor stage and risk score ( Figure 5 3.6. Functional Analysis in the TCGA and ICGC Queues. To explore the potential molecular mechanism of the signature, we used the DEGs to perform GO enrichment and KEGG pathway analyses in the high-and low-risk groups. As we expected, the DEGs were significantly enriched in many immune-related biological processes, such as phagocytosis and leukocyte migration in the TCGA cohort ( Figure 6(a)). DEGs were also gathered in a few hypoxia and immunerelated molecular functions in the TCGA and ICGC databases, for instance, heme binding and antigen binding (P < 0:05, Figures 6(a) and 6(c)). For KEGG analysis, the differentially expressed HAIRGs were gathered in essential pathways associated with immune and cancer progression, for example, proteoglycans in cancer, phagosome, complement and coagulation, and PI2K-Akt signaling pathway (P < 0:05, Figures 6(b) and 6(d)).
We quantified the enrichment score of different immune cell subsets, associated functions or pathways using ssGSEA to further investigate the correlation between risk score and  11 Disease Markers immune status. In the TCGA cohort, there were significant differences in the antigen recognition process (including the score of T cell coinhibition, HLA, APC coinhibition, IDCs, and ADCs) between the high-risk group and the low-risk group (Figures 7(a) and 7(b)). The immunobiological processes in the GO analyses scored higher in the highrisk group in the TCGA cohort (P < 0:05, Figure 6(a)). However, scores of the type I IFN response, type II IFN response, NK cells, and mast cells were lower in the highrisk group, compared to that in the low-risk group. In contrast, scores of macrophages, Treg cell, TIL, and HLA were higher in the high-risk group (P < 0:05, Figures 7(a) and  7(b)). The differences of aDCs, DCs, iDCs, pDCs, Th2 cells, Treg, APC costimulation or inhibition, checkpoint, and type II IFN response between the two groups through comparisons in the ICGC database (P < 0:05, Figures 7(c) and 7(d)). In the TCGA and ICGC cohorts, the macrophage and ADC score existed the greatest statistical difference in the high-and low-risk groups, which were also consistent with the results from the GO analysis.

Prediction of Immunotherapeutic Response and Immune
Checkpoint Expression Pattern in HCC Patients. There are six main immune checkpoints PD1, PDL1, PDL2, CTLA4, CD80, and CD86; their expression level of HCC were detected and compared in the TCGA and ICGC databases. All six immune checkpoints were over expression in the high-risk group in the TCGA database (Figures 8(a)-8(f)). Besides, the six immune checkpoints, PD1, PDL1, PDL2, CTLA4, CD80, and CD86, were also over expression in the high-risk group in the ICGC database (Figures 8(g)-8(l)). So, the expression of immune checkpoint PD1, PDL1, PDL2, CTLA4, CD80, and CD86 in the high-risk group was significantly higher than the one in the low-risk group.
Later, the possibility of predicting immunotherapeutic responses is based on the tumor immune dysfunction and         Then, we found that high TIDE value was related to the high-risk group, while low TIDE value was related to the low-risk group in the TCGA database (P < 0:05, Figure 9(a)) and ICGC cohort (P < 0:05, Figure 9(c)). The proportion of immunotherapy response was significantly higher among patients in the low-risk group than the patients in the high-risk group in TCGA (P < 0:05, Figure 9(b)) and ICGC cohort (P < 0:05, Figure 9(d)), which indicated that immunotherapy has a better response rate on the low-risk group.
3.8. Sensitivity of Chemotherapy Drugs between TCGA and ICGC Cohorts. Adjuvant chemotherapy is an alternative treatment of HCC. Thus, it is important to predict the sensitivity of chemotherapy, which can help clinicians to choose the best chemotherapy regimen. The calculated IC 50 levels of Cisplatin (P = 4:3e − 05) were lowered in the high-risk group than in the low-risk group. AMG 706 (P = 3:8e − 16), gefitinib (P = 2:3e − 05) and docetaxel (P = 1:9e − 11) was lowered in the low-risk group in the TCGA database (Figures 10(a)-10(d)), indicating the high-risk group was more susceptible to cisplatin while the low-risk group was more sensitive to AMG 706, gefitinib, and docetaxel drugs. The estimated IC 50 levels of cisplatin (P = 0:026) were lowered in the high-risk group, but AMG 706 (P = 2:7e − 13) and docetaxel (P = 2:6e − 06) was lowered in the low-risk group in the ICGC database (Figures 10(e)-10(h)). In brief, the high-risk group was more susceptible to cisplatin, but the low-risk group was more susceptible to AMG 706 and docetaxel.

Discussion
In our study, 17 HAIRGs were firstly identified to be related to overall survive in HCC patients. Then, basing on the

20
Disease Markers cancer patients and normal controls. The GO results showed differentially expressed HAIRGs were significantly enriched in phagocytosis, leukocyte migration, immune responseactivating cell surface receptor signaling pathway, neutrophil activation, heme binding, and antigen binding (in terms of BP and MF). These results suggested that differentially expressions of HAIRGS were related to the development of liver cancer. Moreover, for KEGG analysis, it showed differentially expressed HAIRGs were involved in the development of liver cancer, which also indicated HAIRGs may regard as a potential biomarker for HCC. Then, immunotherapy response predication exhibited that the low-risk group was associated with low TIDE score and had a better immunotherapy effect in the TCGA and ICGC cohorts. Fur-thermore, for drug sensitivity analysis in HCC patients, we found that the high-risk group was more susceptive to cisplatin, compared to the low-risk group in the TCGA and ICGC cohorts. Therefore, our result firstly provides a HAIRG signature model which can be utilized for prognostic, immunotherapy response, and chemotherapy sensitivity prediction in HCC.
Currently, the most commonly use of the method for predicting liver cancer is Okuda System Staging. However, this model has some limitations. Over time, the early diagnosis of HCC has changed due to improved diagnostic methods. At the same time, Okuda staging is insufficient to stratify patients prior to radical or palliative treatment [10]. A new and reliable prognostic model is an urgent needed

21
Disease Markers to predict OS in HCC patients. In our study, we found the potential role of the HAIRGs in HCC patients and the possibility of constructing a prognostic model with these HAIRGs. Besides, the 3-year predictive ability of the nomogram for OS was 0.708 in the TCGA cohort and 0.722 in the ICGC database for the model.
The prognostic model consisted of HAIRGs in the present study, including VEGFA, CTNNB1, SPP1, PPARG, HMOX1, RAC1, HSP90AA1, and LGALS3. VEGFA, also called vascular endothelial growth factor-A, is not the only, major factor driving tumor vascular bed dilation [11]. Angiogenesis in HCC depends mainly on VEGFA-driven response, which leads to vascular system dysfunction to a large extent. The reason for this is not clear, although it seems that some aspects of the angiogenic environment stimulated by VEGFA-stimulated angiogenic milieu (high levels of microvascular permeability and density) are capable of promoting tumor expansion [11]. CTNNB1 may play a vital role in metabolic reprogramming and cell proliferation in HCC. Phosphorylation sites associated with CTNNB1 mutations were confirmed on key metabolic enzymes, including ALDOA, and the function of phosphate-ALDOA about promoting metabolic reprogramming and cell proliferation was demonstrated [12]. Secreted phosphoprotein-1 (SPP1) is a secreted arginine-glycine-aspartate (RGD) containing phosphoprotein [13]. SPP1 may play a role as a miR-181c-targeted growth promoter of HCC [14]. PPARG had three subtypes, namely, PPARG1, PPARG2, and PPARG3 [15]. The susceptibility to HCC may be affected by PPARG gene polymorphism [16]. The expression of PPARG was upregulated in lung, prostate, colorectal, bladder, and breast tumors [17]. Heme oxygenase-1 (HMOX-  , an important catalytic enzyme in heme degradation, is increased under stressful conditions [18]. The elevated expression of HMOX-1 in a variety of malignant tumors is related to the tumor microenvironment resistance to tumor cell growth, angiogenesis, metastasis, chemotherapy, and radiotherapy [19]. It suggested that HMOX-1 is a protective gene in liver cancer. HSP90AA1 had the vital functions in the process of the assembly, manipulation, folding, and degradation of its customer proteins [20]. It had been shown to be overexpressed in multifarious human cancer, including liver, breast, endometrial, ovarian, colon, lung, and prostate cancers [21][22][23][24]. LGALS3 (galectin-3) is a multifunctional protein, which has a variety of biological functions, including tumor cells proliferation and differentiation, angiogenesis, tumor progression, and metastasis [25]. Galectin-3 is associated with a lot of cancers, such as mesothelioma, breast, HCC, and colon cancers [26][27][28]. In summary, seven of the genes (VEGFA, CTNNB1, SPP1, PPARG, RAC1,

24
Disease Markers HSP90AA1, and LGALS3) were highly unregulated in patients with liver cancer. Only one gene (HMOX-1) was downregulated, hinting which may be a protective gene for HCC. Whether these genes affect the prognosis of HCC patients by influencing hypoxia, angiogenesis, and immune processes remained to be elucidated.
In functional analysis, we found that different HAIRG signatures of HCC shown a significantly distinguished immune microenvironments, such as immune infiltration levels which including aDCs, DCs, iDCs, pDCs, Th2 cells, and Treg. Moreover, the functional differences that T cell costimulation or inhibition, checkpoint, and type II IFN response were also found in the different groups of HAIRG signature in HCC. The IFN pathway was closely related to the progression of HCC patients. The interferons are divided into three types, type I (IFN-α and IFN-β), type II (IFN-γ), and type III. Currently, it is known that interferon-(IFN-) α is one of the vital treatment options for patients with liver cancer [29]. IFN-α activates interferon-stimulating gene (ISG) transcription by binding to the receptor and mediating its signal transduction. These genes determine the biological consequences of STAT1 signaling and mediate immune functions, inhibit cell proliferation, and induce apoptosis [30]. Moreover, a recent study showed that IFN-α can inhibit growth and induce apoptosis of HCC [31]. IFN-γ is mainly released by T cells that are recognized and activated by antigens [32], which can induce the expression of B7-H1 gene in lung cancer cells, bile duct cancer cells, head and neck cancers, and HCC through JAK/STAT1 pathway [32][33][34][35]. More mechanisms between HAIRGs and IFN pathway regulation in HCC need to be further explored.
Our study also found out different HAIRG signatures correlated with immune checkpoints expression in HCC, which include PD1, PDL1, PDL2, CTLA4, CD80, and CD86. The high expression of CTLA-4 on Tregs in HCC patients was negatively correlated with the cytolytic granzyme B produced by CD8+ T cells [36]. Meanwhile, both tumors and peritumoral cells (LSECs and HSCs) express the ligands PDL1 and PDL2, resulting in inactivation of CD8+ T cells that adhere to hepatic sinusoidal cells, which promotes immune tolerance. Besides, PDL1 expression in HCC also leads to follicular helper T cell failure and impairs the expression of cytokines and the help of B cells, therefore, promoting the development process to advanced tumor stage [36]. PD-1 may play an important function in promoting cancer development. PD-1 blocking combined with targeting mTOR pathway may enhance the antitumor curative effect in cancer [37,38]. Erin et al. found that PD-1 blocking and IL-2 combined therapy could synergistically increase CD8+ T cell response [39]. Thus, there may be a mechanism between calcineurin inhibitors and PD-1 blockades that can serve as therapeutic sites for anti-HCC immunity [36]. However, due to the overexpression of miR-221, the expression of CD86 and CD40 on the surface of DCs downregulated, and miR-221 could delay the maturation of DCs in the microenvironment of HCC cells [40]. As the HAIRG signatures can predict immune checkpoint expres-sion in HCC, our result also found that HAIRG signatures could predict the immunotherapy response in HCC patients. Therefore, our finding indicated the HAIRG signature could be a potential biomarker for the prediction of immunotherapy in HCC.
At last, we also explored the chemotherapeutic sensitivity of HCC patients based on HAIRG signature with traditional chemotherapeutic agents (cisplatin, AMG 706, gefitinib, and docetaxel). This study indicated that the high-risk patients might do better with cisplatin drug than the low-risk patients. However, when using the other drugs (AMG 706 and docetaxel), the low-risk patients did better as well as the high-risk patients. Cisplatin is a frequently used first-line chemotherapy drug for the treatment of liver cancer, ovarian cancer, small cell lung cancer, and other cancers [41]. Docetaxel has been widely utilized to relieve symptoms of breast, prostate, bladder, and ovarian cancers [42][43][44]. Furthermore, docetaxel has been recognized for its low toxicity and high efficacy in the treatment of liver cancer. Compared with previous studies [45][46][47], although had reported the hypoxic signatures, which included genes that PDSS1, SLC7A11, CDCA8, and angiogenesis-related immune signatures, such as BIRC5, KITLG, PGF, SPP1, and SHC1. These signatures of HCC had been reported can be used as potential biomarkers for diagnosis, prognosis, and recurrence of HCC. In our model, which had included more key tumor biological characteristics, included hypoxic, angiogenesis, and immune. More key tumor biological characteristics genes included in our model may construct a more accurate model for predicting the prognosis and treatment sensitivity in HCC. Our results also conducted immunotherapy response and chemotherapy sensitivity prediction in HCC. Therefore, we believe our model may be more novel and accurate. Therefore, our study firstly provides an effective and novel model to evaluate the chemosensitivity for HCC

Conclusion
In a word, our study defined a new prognostic model of HAIRGs in HCC patients. In the TCGA and ICGC cohorts, this model was shown to be independently associated with OS and can be considered a biomarker for prognosis prediction, clinical immunotherapy evaluation, and chemotherapy selection for HCC.

Data Availability
The data and materials used to support the findings of this study are available from the corresponding author upon request.