Screening and identification of angiogenesis-related genes as potential novel prognostic biomarkers of hepatocellular carcinoma through bioinformatics analysis

Clinical characteristics of HCC patients

Figure 1 is a flow chart of this research. After excluding patients with a follow-up time of 0 days or lack of clinical data, our study finally included 365 HCC patients from the TCGA dataset as a training set, and 231 HCC patients from the ICGC dataset as a validation set. The demographic and specific clinical features of the HCC samples in both cohorts were listed in Table 1.

Identification of prognostic angiogenesis-related DEGs in the TCGA dataset

We screened out 52 genes from 79 angiogenesis-related genes with significant differential expression between adjacent normal tissues and tumor tissues (Supplementary Table 1). With the help of univariate Cox regression analysis, 25 genes were closely related to the prognosis of HCC (Figure 2A). By taking the intersection of both sets of genes, it was clear that 16 DEGs were tightly bound to the overall survival (OS) of HCC patients (Figure 2B, 2C). MMP9, PGF, TGFB1, VEGFA, ANGPT2, CTNNB1, PDCD10, AGGF1, ANGPT1, ITGAV had a higher hazard ratio, and were highly expressed in HCC tissues, indicating that the expression of these genes could promote the secretion of angiogenic factors, and angiogenesis was also accelerated with abnormalities of tumor blood vessels, leading to rapid tumor development and poor prognosis. On the contrary, the hazard ratio of TEK, ENG, COL18A1, IL1RN, PLG, and PON1 was less than 1. Moreover, these 6 genes were highly expressed in normal tissues, and had the function of normalizing blood vessels. High expression of these genes is conducive to prognosis. VEGFA and TEK were clearly identified as the core genes according to the interaction network (Figure 2D). Meanwhile, the correlation diagram based on these angiogenesis-related genes was exhibited in Figure 2E.

Establishment of a prognostic signature using the TCGA dataset

We conducted the least absolute shrinkage and selection operator (LASSO) regression to filter suitable genes from the above expression of the 16 genes to establish a prognostic signature. A 7-gene signature was identified on the basis of the optimal value of λ (Supplementary Figure 1A, 1B). We calculated the risk score as shown below: risk score = (0.123128 × expression of ANGPT1) + (–0.190501 × expression of ENG) + (0.282949 × expression of PDCD10) + (0.122334 × expression of PGF) + (–0.033803 × expression of COL18A1) + (0.016683 × expression of ITGAV) + (–0.066838 × expression of PON1). In terms of the median risk score, 365 HCC patients were dichotomized into a high-risk group (n = 182) and a low-risk group (n = 183) (Figure 3A). A Kaplan-Meier curve was created to indicate that patients in the low-risk group had a better OS probability than those in the high-risk group (P < 0.001) (Figure 3B). The OS rates at 1-, 2-, 3-year for the high- risk group were 69.7%, 53.0%, and 42.2%, whereas the corresponding rates for another group were 95.4%, 85.7%, and 80.9%, respectively. Other similar studies focused on the analysis of disease-specific survival (DSS), progression-free survival (PFS), and disease-free survival (DFS) of these HCC patients, whose outcomes of survival status and Kaplan-Meier curves were highly consistent with the results of OS (Supplementary Figure 2A–2C). To estimate its predictive performance, we displayed the time-dependent receiver operating characteristic (ROC) curve, with the area under the curve (AUC) reaching 0.785, 0.722 and 0.715 in the 1 year, 2 years, and 3 years, respectively, indicating that the prognostic signature exhibited an outstanding specificity and sensitivity (Figure 3C). To study this signature in-depth, we divided the patients into subgroups based on their clinical characteristics (age <=65 and >65, female and male, stage I and II and III and IV, grade G1 & 2 and 3 & 4), and compared the divergences of OS between the two groups in each subgroup. The Kaplan-Meier curves of OS in each subgroup had the same trends as the OS of all samples, suggesting that the signature could be applied to different clinically specific populations (Supplementary Figure 3A–3H).

To further investigate the genes involved in the signature, we examined a human liver cell line (L02) and a hepatoma cell line (HepG2) for the expression of these genes through qRT-PCR. A series of specific primers of each gene were shown in Table 2. The relative expression of each gene between the L02 cell line and HepG2 cell line confirmed by qRT-PCR was not completely consistent with the corresponding RNA-sequence from the TCGA cohort, in light of the high heterogeneity of HCC. The expression of ANGTP1, COL18A, ITGAV, PGF in HepG2 cells were higher than those of L02 (Figure 4A–4D), while ENG and PON1 were highly expressed in hepatocytes (Figure 4E–4F). There was no difference in PDCD10 expression in both cell lines (Figure 4G). We further observed the immunohistochemistry of these 7 proteins in normal and pathological sections through the Human Protein Atlas (HPA) database, which lacked PGF data (Supplementary Figure 4). The comparison results were highly consistent with paired transcriptome sequencing differences (Figure 5A–5L). Survival analyses, on the basis of the optimal cut-off value of each gene expression, revealed that the differential expression of these genes was relevant to the prognosis of HCC patients. Particularly, the relatively high expression of ANGPT1, ITGAV, PDCD10, and PGF were in association with poor prognosis, whereas the higher expression of COL18A1, ENG, and PON1 had relation to a better prognosis, which was completely consistent with the sign of the corresponding gene coefficient in the signature (Supplementary Figure 5A–5G). There was an extremely obvious difference in survival analysis between various risk groups divided by the optimal cut-off configuration as well (Supplementary Figure 5H).

Meanwhile, Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis demonstrated that HCC patients in various risk groups were distributed in different directions (Figure 3D, 3E). Patients belonging to the latter had lower mortality as well as better prognosis; with the risk score increasing, the probability of death elevated obviously and the OS was significantly shortened (Figure 3F).

Validation of the 7-gene signature in the ICGC dataset

The HCC patients in the ICGC dataset were organized into a high-risk group (n = 115) and a low-risk group (n = 116) in the light of the median value, calculated with the above-mentioned model formula, to estimate the robustness of the signature what has been established (Figure 6A). Kaplan-Meier analysis suggested that there was a remarkable difference in the survival probability of the two groups (P < 0.001) (Figure 6B). The AUC of the 7-gene signature for 1-, 2-, 3-year OS were 0.764, 0.705, and 0.724 respectively (Figure 6C). Similarly, both t-SNE and PCA analysis illustrated that the distribution of these samples in the two risk groups presented discrete directions, which corroborated the results of the TCGA cohort (Figure 6D, 6E). Patients in the high-risk group were confronted with a shorter survival time and more deaths compared with their low-risk counterparts (Figure 6F). Apart from ANGPT1, ITGAV, and PGF, other genes and risk scores were fully related to OS and consistent with the trend of the TCGA outcomes (Supplementary Figure 6A–6H). We divided these HCC patients into different subgroups according to their clinical characteristics and performed survival analysis. And the signature was also applicable to the population with different clinical characteristics in the ICGC dataset (Supplementary Figure 7A–7F).

Independent prognostic value of the 7-gene signature

After removing samples with incomplete clinical information, the further analysis involved several clinical characteristics of 240 HCC patients, containing age, gender, grade, and stage. We applied univariate and multivariate Cox regression analysis to assess the effectiveness of independent prognostic predictions of risk scores and other clinical characteristics. Univariate analysis indicated that risk score was considerably correlated with the OS probability in the TCGA dataset (HR = 6.160, 95% CI = 3.568–10.633, P < 0.001) and the ICGC dataset (HR = 5.953, 95% CI = 2.774–12.775, P < 0.001). The risk score was still proved to be an independent predictor for OS in the TCGA dataset (HR = 5.430, 95% CI = 2.998–9.833, P < 0.001) and the ICGC dataset (HR = 4.376, 95% CI = 1.987–9.637, P < 0.001), after correction for other confounding factors through the multivariate Cox regression analysis. Besides, stage was also an independent prognostic factor for predicting OS (TCGA dataset: HR = 1.996, 95% CI = 1.361–2.927, P < 0.001; ICGC dataset: HR = 2.559, 95% CI = 1.330–4.925, P = 0.005). The detailed information of the independent prognostic analyses was presented in Table 3, 4.

To illustrate the clinical application value of this signature, we compared the differences in risk scores for various clinical characteristics from the TCGA cohort. There were no remarkable differences in risk scores in patients of various ages and genders, indicating that the two have no additional contribution to the prognostic risk in this model (Figure 7A, 7B). It was worth noting that there are obvious differences in the risk scores of this model in different grades, except for G3 and G4, and as the grade increased step by step, the risk scores also increased (Figure 7C). This signature could effectively distinguish the different grades of HCC, so that could be a potential biomarker for HCC grades. The low-risk group was represented by G2, and the high-risk group corresponded to G3 (Figure 7D). Accompanied by the increase of stage and T classifications, risk scores also had a corresponding increase trend; notably, compared with stage II, stage III and other T classifications, the risk scores of stage I and T1were pronouncedly lower (Figure 7E–7H). The angiogenesis status of various stages in the ICGC cohort also had similar outcomes (Figure 7K, 7L). For tumor status, the risk scores of patients with tumor were higher than that of patients with tumor-free, suggesting that the differences in high and low scores of patients in this signature revealed that the total resection could be performed, which is also closely related to the angiogenesis (Figure 7I, 7J).

A personalized prognostic prediction model

We applied the nomogram to quantitatively estimate the individual’s actual clinical survival risk by integrating multiple risk factors [19]. We constructed a nomogram by integrating the 7-gene signature, age, gender and TNM classification to forecast the prognosis at 1, 2, and 3 years (Figure 8A). The length of each line represented its degree of influence on the prognosis, and different values corresponded to different results. Especially, the risk score corresponded to the longest line, indicating that it had the strongest predictive ability for OS. The 1-, 2-, and 3-year survival probability of the patient could be judged based on the total points obtained by adding the point of the risk score related to the patient’s prognosis and the points corresponding to the clinical characteristics. The calibration curves implied that the predicted survival probability matched with the actual one well (Figure 8B–8D).

Functional analysis of the angiogenesis-related 7-gene signature

We performed GO and KEGG analyses to further clarify the biological pathways and functions relevant to the 7-gene signature through the functional enrichment analysis of DEGs. Conspicuously, DEGs were enriched in several functions concerning cell proliferation and mitosis, as well as nuclear division, which were closed in association to angiogenesis in both datasets (adjusted P < 0.05) (Figure 9A, 9B) [20–22]. For biological processes (BP), the major enriched GO terms were nuclear division and organelle fission, and other immune-related pathways, such as neutrophil activation involved in the immune response. The most enriched cellular component (CC) was the spindle and chromosomal region. And it revealed that primary functional categories in molecular function (MF) were microtubule binding and single-stranded DNA binding. Therefore, the application of anti-mitotic drugs such as colchicine may effectively inhibit tumor growth and metastasis. KEGG pathway primarily involved in PI3K-Akt signaling pathway, cell cycle, microRNAs in cancer, ECM-receptor interaction, as well as other cancer-promoting pathways (adjusted P < 0.05) (Figure 9C, 9D).

Previous studies had shown that tumor anti-angiogenesis therapy went hand in hand with immunotherapy [23–25]. With the normalization of tumor vasculature (or anti-angiogenesis), tissue perfusion was enhanced, and the infiltration of immune effector cells was improved, leading to immunotherapy potentiation; meanwhile, the stimulation of immune cell functions also contributed to normalizing tumor blood vessels. Thus, ssGSEA was utilized to perform enrichment analysis on different immunocyte subsets, related pathways and functions to investigate the relationship between risk scores and immune cells infiltration with and related components. The scores of adaptive immune cells, including CD8+ T cells, B cells, Neutrophils, and TIL, were different between the two groups in the TCGA cohort; notably, the expression values of these cells in the high-risk group were lower than those in the low-risk group (all P < 0.05) (Figure 10A). Some molecules and signals related to the antigen presentation process, such as MHC class I, aDCs, DCs, pDCs, and APC co-inhibition, also had gaps between various groups in the TCGA dataset, but it was worth noting that their regulation directions were not consistent (all P < 0.05) (Figure 10A, 10B). What’s more, the high-risk group had lower scores of NK cells, cytolytic activity, type I and II IFN response, and mast cells, while the score of macrophages was just the opposite (all P < 0.05). Similarly, the differences of B cells, Neutrophils, macrophages, NK cells, cytolytic activity, type I and II IFN response between the two groups were verified in another dataset (all P < 0.05) (Figure 10C, 10D).

To further investigate the relation between the model and immunotherapy, the immunophenoscore was utilized to evaluate the distinction in the effect of immunotherapy between the two groups under various immune checkpoint conditions. For HCC patients with both CTLA4 and PD-1 double-positive or double-negative, or PD-1 positive but CTLA4 negative, the IPS of patients in the high-risk group was noticeably lower than that of the low-risk group, that is, under the above-mentioned conditions, patients in the former had a worse effect on immunotherapy (Figure 11A–11C). Yet, a higher IPS was exposed to the patients (with CTLA4 positive but PD-1 negative) in the high-risk group (Figure 11D), resulting from the higher expression of CTLA4 in this group (Figure 11E). Besides, we evaluated the risk scores of the HCC patients corresponding to different immune subtypes (C1: Wound Healing, C2: IFN-gamma Dominant, C3: Inflammatory, C4: Lymphocyte Depleted, C5: Immunologically Quiet, C6: TGF-beta Dominant) (Figure 11F). The immune subtype C1, which is characterized by wound healing, had a remarkably higher risk score than other subtypes. The immune subtype C3 had a relatively low risk score, resulting in a better prognosis, which was completely consistent with existing research conclusions [26].

In general, the immune cell infiltration and corresponding immune response of samples in the low-risk group were dramatically more striking, indicating that relatively normal tumor blood vessels are more likely to be immune cell infiltrated and even beneficial for immunotherapy and combination therapy of anti-angiogenesis therapy and immunotherapy.

Transcriptome data and clinical data collection

Transcriptome data, containing gene expression, together with clinical characteristics of American HCC samples (n = 371) and Japanese HCC samples (n = 231) were respectively collected from the TCGA portal (https://portal.gdc.cancer.gov/repository) (up to July 1, 2020) and the ICGC portal (https://dcc.icgc.org/projects/LIRI-JP) (up to August 1, 2020). These Japanese samples were mainly derived from patients with HCC caused by HBV or HCV infection [50]. These extracted RNA-sequence data were normalized with formula log₂ (x + 1) by the “limma” R package. After removing samples with a follow-up time of 0 days and missing clinical characteristics, 365 HCC samples in the TCGA portal, with intact survival status, OS time, age, gender, histological grade, and TNM stage were included. The samples in the ICGC database had complete clinical data. Then, 79 angiogenesis-related genes with a correlation score > 10, provided in Supplementary Table 2, were extracted from the GeneCards database (https://www.genecards.org/).

Establishment and validation of a potential prognostic angiogenesis-related gene signature

The “limma” R package was used to screen out the DEGs between the tumor tissues and the adjacent normal tissues in the HCC samples in the TCGA data with the false discovery rate (FDR) < 0.05. Then, we screened angiogenesis-related genes that could be utilized to judge the prognosis of HCC patients through univariate cox analysis. The intersection of these two groups of genes was utilized to identify the differentially expressed angiogenesis-related genes with prognostic value in HCC. The STRING database generated an interactive network that overlaps the prognostic DEGs [51]. The interactions with high confidence (0.70) were considered statistically significant. LASSO regression was performed with the “glmnet” R package to establish a prognosis signature to reduce the dimension and enhance its generalization ability [52, 53]. The response variables in the LASSO regression were OS probability and survival status, and the independent variables were the normalized expression matrix of the above candidate prognostic DEGs. Based on tenfold cross-validation following the minimum criteria, we identified the penalty parameter (λ) for the signature. The following computational formula was then established for further analyses:

$$\text{Risk score = }{\displaystyle \sum _{i=1}^n[ExpressionofGen{e}_i\times {\beta }_i]}$$

Where n represents the number of finally enrolled genes, β indicates the coefficient of Gene_i obtained after the LASSO regression.

The patients in both portals were divided into two risk groups for survival analysis in line with the median value of the risk score. And Kaplan–Meier survival analysis was utilized to compare the OS, DFI, PFI, DSS between these two groups. In addition to the above survival analysis, we also utilized the “survminer” R package to identify the optimal cut-off expression to perform survival analysis for each gene and risk score. PCA was carried out to explore whether the difference between the two groups was significantly identified on the basis of unsupervised learning through the “stats” R package. Such dimensionality reduction can extract the portions of the principal components of various groups, so as to clarify whether there were remarkable differences between the two groups. In order to determine the complex structural relationship between features, we performed t-SNE to investigate the distribution of the two groups with the “Rtsne” R package. The time-dependent ROC curves of 1-, 2-, 3-year was implemented to assess the predictive power by the “survivalROC” R package. The paired expression value of angiogenesis-related genes in the signature, along with immunohistochemistry of the expressed proteins were obtained from the HPA database.

Detection of the expression of each gene in the cell lines by qRT-PCR

We extracted total RNA extracted from the human hepatocarcinoma cell line HepG2 and normal liver cell line L02 using the RNA-Quick Purification Kit (Shanghai Yishan Biotechnology Co., Ltd, China). The A260/A280 absorption (1.9–2.2) was used to evaluate the quality of the extracted total RNA. Then, the HiScript III All-in-one RT SuperMix Perfect for qPCR kit (Vazyme, China) was utilized to reverse transcribe the total RNA as following steps: 50°C for 15 min and 85°C for 5s. The Taq Pro Universal SYBR qPCR Master Mix kit (Vazyme, China) was utilized to perform PCR according to the following reaction procedure: 95°C for 30s, 95°C for 10s with 40 cycles, and 60°C for 30s. The specific primer sequences of GAPDH and these genes were displayed in Table 2. We calculated the relative expression of target genes using the 2^−ΔΔCT method.

Functional enrichment analysis

GO and KEGG analyses were conducted to clarify the enrichment of related functional pathways on the basis of the DEGs (FDR < 0.05, |log2FC| ≥ 1) using the “clusterProfiler” R package. The ssGSEA could be utilized to calculate the infiltrating score of 16 immune cells and the activity of 13 immune-related pathways [54]. Supplementary Table 3 displays the annotated gene set file.

Development of nomogram

The nomogram can predict the probability of a certain clinical outcome on the basis of the values of multiple variables [55]. To establish a nomogram for the TCGA dataset, we integrated age, gender, TNM classification and risk score with the “survival” and “rms” R package. Then, we plotted calibration curves to estimate the concordance between predicted and actual survival probability. Combining clinical characteristics and risk scores could predict 1-, 2-, 3-year survival rates of HCC patients.

Statistical analysis

All statistical analyses were performed with R software version 4.0.1 (https://www.R-project.org). Unless otherwise noted, p < 0.05 was considered to be statistically significant. We used Chi-square tests to evaluate the differences in proportions, and independent t-tests to assess the differences in gene expression levels between adjacent normal tissues and tumor tissues of HCC samples. Mann-Whitney test was utilized to evaluate the ssGSEA scores of pathways and immune cells between different groups for comparison. Kaplan-Meier curves were usually utilized to analyze the survival probability between different groups. Wilcoxon test and Kruskal-Wallis test were performed to compare the differences in risk scores and IPS between different groups. We implemented Cox regression to distinguish independent predictors of OS probability.