The raw counts of RNA-Seq data and corresponding clinical information of TCGA-LIHC patients were collected as the training cohorts. Data were downloaded from UCSC Xena (http://xena.ucsc.edu/). We also obtained two independent validation cohorts of HCC patients from the ICGC portal (ICGC-LIRI-JP, https://dcc.icgc.org/projects/LIRI-JP) and the GEO (GSE144269, https://www.ncbi.nlm.nih.gov/geo/), with transcriptomic and clinical data available.

We first collected genes with annotations related to “angiogenesis” in NCBI. Then the signature gene sets of “angiogenesis” from the MSigDB database were also obtained (http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp). The final AFs gene set studied here was the combination of genes from the two resources (Supplementary Table S1).

In order to obtain the differential expressed genes between tumor and normal samples, Log-transformed fold change (FC) and FDR of each gene was analyzed by DESeq2 (version 1.34.0) package (Love et al., 2014). The two “recurrent solid tumor” samples were removed from analysis. Genes with FDR < 0.05 and | log2 FC|>1 were defined as differentially expressed genes.

Univariate Cox regression was performed for each AF genes to obtain the prognostic genes with p-value < 0.01 and HR > 1 or HR < 0.5. Kaplan-Meier analysis was also performed to screen the prognostic candidate AF genes using R package survival (version 3.2-13) (Borgan and Therneau, 2001). Multivariate Cox regression was used to assess the performance of risk score under the effects of other clinical factors.

The prognostic candidate AF genes were screened by univariate cox regression and log-rank test. Tumor samples of the TCGA-LIHC were used as the training cohort to establish the LASSO model. A lasso penalty was used to find the best gene model utilizing an R package glmnet. The risk score for each sample can be calculated with the final LASSO model coefficient as follows: R i s k   S c o r e =   ∑ i = 1 n exp ⁡ r e s s i o n   o f   g e n e   i ∗ l a s s o   c o e f f i c i e n t   o f   g e n e   i

The median risk score was used as cutoff for high-risk group (with higher risk score) and low-risk group (with lower risk score).

CIBERSORT algorithm was used to evaluate the infiltration of 22 immune cell types (Newman et al., 2015). The immune score and stromal score were tested by the R package ESTIMATE (Yoshihara et al., 2013). TIDE was used to predict the potential of patients to response for immunotherapy (Jiang et al., 2018). We also collected tumor stemness score, TMB and HRD score of TCGA tumors from previous studies (Chen et al., 2021). mDNAsi, EREG-mDNAsi, DMPsi and ENHsi data of TCGA-LIHC tumor samples were collected from exist studies (Malta et al., 2018).

Mutation comment file (MAF) of TCGA-LIHC cohort was downloaded from the GDC client. Differential analysis and visualization of somatic mutations were performed using maftools package. The Fisher’s exact test was used on all genes between two groups to detect differentially mutated genes. Segment file of the TCGA-LIHC cohort was downloaded from FIREHOSE and analyzed using the GISTIC 2.0 pipeline (Mermel et al., 2011).

AFs-derived risk scores, TNM stages, clinical stage, gender, age and grade were used as independent prognostic factors through univariate cox regression and AFs-derived prognosis risk model. Nomogram was finished based on the results of multivariate cox regression analysis. The calibration curves of the nomogram were constructed to test consistency between 1-, 3- and 5-years survival probability prediction and actual observation. The performance of the nomogram was evaluated using the concordance index (C-index) and time-dependent receiver operating characteristic (ROC) curves. Nomograms analysis and visualization were performed using R packages rms (version 6.2-0) (JrHarrell, 2021) and survival (version 3.2-13) (Borgan and Therneau, 2001) with default parameters.

First, the large population of liver hepatocellular carcinoma (LIHC) patients from the TCGA database was used as the training cohort. We downloaded the transcriptomic and clinical data for 371 tumor samples and 50 adjacent normal samples. A total of 8,250 differentially expressed genes were found between tumor and normal samples (FDR <0.01, |log2FC| > 1, DEseq2). We then systematically collected AF genes from NCBI and MSigDB (see Methods) and noticed 1,038 AFs were differentially expressed (hereafter termed DE-AFs, Figures 1A,B). Among them, 361 DE-AFs were downregulated and 677 DE-AFs were upregulated in tumors.

Then we investigated the prognostic ability and regression characteristics of DE-AFs using univariate Cox regression analysis. A total of 36 genes with significant differences were discovered in univariate regression model (p-value < 0.01, hazard ratio (HR) > 1 or HR < 0.5, Supplementary Table S2). Furthermore, the Kaplan-Meier curve with log-rank test was also used to filter important DE-AFs and finally screened 17 genes with confirmed prognostic efficacy (log-rank p-value < 0.05, Supplementary Figures S1 and S2), including EGF, GRM8, TRPM6, SLC38A6, BLM, BARD1, CLSPN, PRIM2, MSH2, FAM72A, SPC25, IGF2BP3, CENPP, GTF2IRD1, TMC7, FSD1L, and SLC39A10. The distribution of expression levels and Kaplan-Meier curves of the four genes with top significance were displayed (Figures 1C–F). It was obvious that higher expression of them indicated poor prognosis of patients.

Although these genes showed certain predictive efficacy by means of the intersection of univariate Cox regression analysis and log-rank test, we preferred combining the informative genes to obtain a more optimized prognostic model. The LASSO regression analysis was subsequently performed to remove redundant factors and also filtered the factors with less contribution. The 6 AFs with the most predictive value were selected (Supplementary Figures S3 and S4), including GRM8, SPC25, FSD1L, SLC38A6, FAM72A, and SLC39A10. In addition, several publications also supported these genes in cancer with experimental evidence (Supplementary Table S3), such as siRNA approach. Then a prognostic risk scoring model of AFs with the coefficients from LASSO regression analysis was constructed (Figure 2A). Their expression was accordantly upregulated to define the high-risk group of HCC (Figure 2B).

Based on the 6-AFs gene prognostic model, patients of TCGA-LIHC were stratified into high-risk group (n = 158) and low-risk group (n = 159) according to the median cut-off value of risk score. The overall survival (OS) time of patients in the high-risk group was remarkably decreased (Figure 2C). The Kaplan-Meier and C-index analysis showed the capacity of the AFs prognostic model (log-rank test p-value < 0.001, C-index = 0.7, Figure 2D). The predictive performance of the prognostic risk model was further evaluated by time-dependent ROC curves, and the area under the ROC curve (AUC) reached 0.802 at 1-year, 0.694 at 2-years, and 0.672 at 3-years (Figure 2E).

After performing the univariate Cox regression analysis (Figure 2F), we explored the relationship between clinical characteristic factors and AFs risk score. Age, gender, stage, grade, and the risk score of the prognostic model were included in the multivariate Cox regression model. The risk score was found to be an independent predictor for OS, with HR = 1.93, 95% CI: 1.26–2.94, p-value = 0.002 (Figure 2G, Table 1). Taken together, the AFs gene prognostic model was confirmed as a credible and independent predictor of OS in HCC.

To assess the robustness and generalizability of the AFs-derived prognostic risk model, the validation data sets from ICGC-LIRI-JP (N = 212) and GSE144269 (N = 68) were collected. Patients were all separated into high- or low-risk groups according to the risk score. The high-risk group of the ICGC-LIRI-JP validation cohort also showed significantly lower survival rate than the low-risk group (log-rank p-value = 0.03). The predictive capacity was proved as AUC reaching 0.602 at 1-year, 0.632 at 2-years, and 0.709 at 3-years in ICGC-LIRI-JP cohort (Figure 3A). Consistently, another validation set from GSE144269 (N = 68) also supported the poor prognosis of the high-risk group (log-rank p-value = 0.015, Figure 3B). Likewise, in the ICGC-LIRI-JP and GSE144269 cohorts, the risk score still proved to be an independent predictor for OS after correction for other confounding clinical factors (multivariate Cox regression analysis, ICGC-LIRI-JP cohort: HR = 2.87, 95%CI = 1.38-5.99, p = 0.005; GSE144269 cohort: HR = 13.14, 95%CI = 1.32-130.43, p = 0.028; Figures 3C,D).

In order to explore the underlying molecular mechanisms of this AFs-derived prognostic risk model, we assessed the associations between the risk score and the typical clinical characteristics including age, gender, TCGA molecular typing (iclust1, iclust2, iclust3), tumor stage, virus infection status, etc. (Cancer Genome Atlas Resea, 2017). Chronic infection of Hepatitis B virus (HBV) has been commonly considered as a major risk factor in the initiation and development of HCC (Chan et al., 2016). We observed the AFs risk score was positively associated with the risk of HBV infection (Wilcoxon rank-sum test p-value = 0.0018, Figure 4A), which further confirmed its predictive value of severe disease status. Besides, the risk score varied in different TCGA molecular typing groups (Kruska-Wallis test p-value = 0.00036, Figure 4B). Higher risk score indicated a higher tumor stage, and the exception of stage IV was probably due to the limited samples size (Kruska-Wallis test p-value = 3.2e-6, Figure 4C). Previous analysis proved that cancer stem cells promoted angiogenesis by secreting factors such as vascular endothelial growth factor (VEGF) and stromal cell-derived factor 1 (SDF1). Therefore, we obtained two measurements depicting the tumor stemness (Malta et al., 2018; Chen et al., 2021), and verified a positive correlation between the tumor stemness and AFs risk score in TCGA-LIHC cohorts (Figures 4D,E).

Next, we investigated the associations between the AFs risk score and the immune response in tumors. The mutation burden (TMB) was not differed between high- and low-risk groups (Figure 4F), while the homologous recombination deficiency (HRD) score was positively correlated with the AFs risk score (Spearman R = 0.53, p-value < 2.2e-16, Figure 4G). This was in line with the moderate potential to induce adaptive immunity in high-risk group (Figure 4H). As for individual immune cell types, we calculated the relative fraction of 22 immune cell types in five cohorts (TCGA-LIHC, ICGC-LIRI-JP, GSE10141, GSE144269 and GSE10186) by the CIBERSORT algorithm. And found that the fraction of immunosuppressive regulatory T cells (Tregs) cells, M0 macrophages, and resting dendritic cells showed positive correlations with the risk score, while several adaptive immune cells showed the trends of negative correlations (Figure 4I).

Furthermore, we compared the genomic aberrations between the different risk groups in TCGA-LIHC cohorts. The difference of copy number variation (CNV) was analyzed through maftools and GISTIC 2.0. As shown in Figures 5A,B, the high-risk group had significantly more deletion events and higher CNV frequencies than the low-risk group. Then, the differentially mutated AFs genes between the high-risk and low-risk groups were detected (Figure 5C, chi-squared test, p-value <0.05). Among them, the mutation frequency of HCC driver gene TP53 was enriched in the high-risk group (45 versus 13%, Fisher’s exact test p-value = 6.60e-10), this observation suggested the classic role of TP53 in cell-cycle regulation and guarding genome stability might also contribute to the malignant progression of HCC (Gao et al., 2019). While the CTNNB1 mutation was higher in the low-risk group than the high-risk group (29 versus 20%), the p-value slightly failed to reach statistical significance (Fisher’s exact test p-value = 0.063). Previously, multi-omics integration analysis revealed three HCC subtypes, one of which exhibited few CTNNB1 mutations companied by poor prognosis (Cancer Genome Atlas Resea, 2017). This was consistent with our observation that relatively lower mutation frequency of CTNNB1 in the high-risk group. In addition, the previous analysis also found that microvascular invasion was significantly reduced in the subtype with increased CTNNB1 mutation. We expected the mechanism of somatic mutations on angiogenesis in hepatocellular carcinoma to be further explored and verified.

Accumulating evidence supported that sorafenib was effective in extending the time of progression in HCC (Vitale et al., 2010). In order to investigate the effect of our risk scoring model in HCC with sorafenib as adjuvant treatment, GSE109211 data set was used for analysis. The proportion of samples responding to sorafenib in the high-risk group was significantly higher than that in the low-risk group (0.53 versus 0.09, Figure 6A). Accordingly, the risk score was found extremely higher in responders (Wilcoxon rank-sum test p-value < 0.0001, Figure 6B). In other words, the risk score of the AFs prognosis model could effectively predict the patient’s response to sorafenib. The results of ROC curve analysis confirmed the good sensitivity and specificity of risk score (AUC = 0.8416, Figure 6C).

TIDE is the tumor immune dysfunction and rejection score, depicting the primary mechanisms of tumor immune evasion. It was proved to predict the clinical response and outcomes of patients following immunotherapy (Jiang et al., 2018). We used TIDE to evaluate the potential of risk score as a predictor of immunotherapy. Interestingly, the high-risk group showed significantly lower TIDE score (Wilcoxon rank-sum test p-value < 0.05, Figure 6D).

Finally, a graphic prognostic nomogram based on the 6-AFs genes was developed for 1-, 3-, 5- and 10-years prediction of OS for HCC patients from TCGA. The tumor stage, grade, age, and gender were also included (Figure 6E). Meanwhile, the calibration plot showed that the prediction by the nomogram had good agreement with actual observation (Figure 6F).