Comprehensive Analysis of Angiogenesis and Ferroptosis Genes for Predicting the Survival Outcome and Immunotherapy Response of Hepatocellular Carcinoma

Introduction

Ferroptosis is a novel form of iron-dependent cell death characterized by impaired intracellular lipid oxide metabolism due to the accumulation of unsaturated fatty acids and lipid reactive oxygen species (ROS) on the cell membrane.¹ In recent years, an increasing number of studies have found that ferroptosis plays an essential role in tumorigenesis, especially in hepatocellular carcinoma (HCC).^2,3 Cancer cells undergoing ferroptosis may suppress the body’s anti-tumor immunity and promote tumor growth by releasing some of the signaling substances such as oxidized lipid mediators.^4,5 Many molecules could modulate the susceptibility of cells to ferroptosis in terms of lipid metabolism, such as Glutathione peroxidase 4 (GPX4).⁶ Inactivation of GPX4 was associated with increased sensitivity to ferroptosis, which significantly inhibited the development of HCC.⁷ Furthermore, anti-PD-1 immunotherapy-activated CD8⁺ T cells could enhance ferroptosis in tumor cells by releasing interferon-gamma (IFN-γ) and further improve immune potency without harming normal cells.⁸ The above indicated that further exploration of the ferroptosis mechanism will contribute to the development of new therapeutic strategies for HCC.

Angiogenesis is a common process in solid tumors, contributing significantly to tumor development, infiltration, and vascular anomalies by giving the tumor tissue oxygen and nutrients for metabolism.⁹ Inhibition of tumor angiogenesis has become the key to modern tumor therapy, and therapies that normalize blood vessels have been approved for first-line treatment. For anti-angiogenesis treatment in advanced HCC, sorafenib and several additional angiogenesis inhibitors (also known as tyrosine kinase inhibitors, or TKIs) have been employed.^10,11 Finding biomarkers of therapy response, however, will aid in the development of individualized treatment plans and enhance patient outcomes as treatment response varies from person to person.

The fact that Evodiamine can inhibit angiogenesis by inducing iron metabolism in prostate cancer¹² shows us that combined ferroptosis and angiogenesis therapy have a promising future in the treatment of tumors. Unfortunately, no studies have reported the combined mechanism of ferroptosis and angiogenesis in HCC. Our goal in this study was to investigate the possible significance of prognostic genes linked to ferroptosis and angiogenesis in determining the prognosis and efficacy of immunotherapy in HCC patients.

Materials and Methods

Public Datasets Acquisition and Genes Associated with Angiogenesis and Ferroptosis

The transcriptome expression data and clinical features of HCC patients were obtained from the TCGA-LIHC, ICGC (LIRI-JP), and GSE14520 datasets, and the clinical characteristics of the patients are listed in Table S1. Samples with complete clinical survival data or survival times of more than one month were included. Genes involved in angiogenesis and ferroptosis pathways were derived from previous studies.^13–15

Molecular Subtypes Identification by Non-Negative Matrix Factorization (NMF) Algorithm

The 342 hCC samples in TCGA were grouped using the NMF algorithm with the criteria “brunet” and 50 iterations based on the angiogenesis and ferroptosis-related genes. The number of clusters (K) ranged from 2 to 6, with cophenetic, dispersion, and profile being used to determine the ideal number of clusters. Kaplan-Meier survival analysis was also undertaken to see if there were any changes in survival across the NMF subtypes.

DEG Identification and Prognostic Risk Score Model Construction

Genes with statistically significant differences between different subtypes identified by NMF were defined as differential expression genes (DEGs), when the P-value was less than 0.05 and |logFC| > 1. The link between overlapping DEGs and HCC patient survival outcomes was calculated using univariate Cox regression. Then, to explore the hub genes, LASSO-COX regression, random forest technique, and extreme gradient boosting (XGBoost) were used. The ten-fold cross-validated LASSO-Cox regression reduces the number of variables in the model and avoids overfitting and covariance compared to the traditional Cox regression model. The selection of genes corresponding to the min(λ) is clinically valuable for subsequent model construction. Moreover, through the XGBoost algorithm, we chose the fifteen most influential genes by ranking the genes based on their relevance. Similarly, the random forest algorithm ranked the genes based on their relevance and subsequently chose the genes corresponding to the number of trees with the smallest error rate for model construction. The intersecting genes among the three machine learning methods were finally put into multivariate Cox regression examination to establish a prognostic model. The multivariate Cox relapse coefficient (β) was used to create a risk score based on the concept of directly mixing the equation below with the mRNA expression level. Score = (Coef_mRNA1 × mRNA1) + (Coef_mRNA2 × mRNA2) +…+ (Coef_mRNAn × mRNAn). Patients in the TCGA and ICGC datasets were categorized into high-risk and low-risk groups based on the median value of the score. Kaplan-Meier survival analysis and receiver operating characteristic curves (ROCs) were subsequently performed to validate the predictive power of the model. Additionally, gene set enrichment analysis (GSEA) was carried out between the groups with high and low risk groups to distinguish the entirely cautious gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) items with FDR 0.05.

qRT-PCR and Immunohistochemistry (IHC) Assay in Clinical Samples

Fifty fresh frozen tumor tissues and 50 paracancerous normal tissues from HCC patients admitted to the Department of Pathology of the Henan Provincial People’s Hospital from 2023 to 2024 were chosen as clinical samples for qRT-PCR investigation. None of the patients received radiotherapy or chemotherapy before surgery, and those with incomplete data were excluded. The ethics committee of Zhengzhou University approved the enrollment of patients. Total RNA was isolated from tissues with TRIzol reagent and reverse-transcribed with a PrimeScript RT Master Mix Synthesis Kit for mRNAs. In a Real-Time PCR system, the relative expression levels of mRNAs were measured using Maxima SYBR Green qRT-PCR Master Mix. The relative mRNA levels of Gene expression were calculated using the delta Cycle threshold (Ct) technique, with β-ACTIN expression levels serving as endogenous controls. Table S2 contains the primer sequences. All participants signed a written informed consent form. After quenching endogenous peroxidase with 2% hydrogen peroxide and blocking with 5% bovine serum albumin, the sections were heated in citrate buffer for antigen retrieval (BSA, A602449, Sangon Biotech, China). The slices were then incubated with an anti-SLC10A1 (ab131084, Abcam) antibody and an isotype control for the remainder of the evening at 4 °C. After PBS washed the sections, DAB chromogen (A690009, Sangon Biotech, China) was used to stain them gradually (Phosphate Buffer Solution with Tween-20). The IHC results assessment was carried out by two skilled pathologists employing single-blind and unified criteria techniques as previously reported.¹⁶ A final score was calculated by combining the extent of expression score (no positive cells = 0, <10% = 1, 10%–50% = 2, >50% positive staining = 3) and intensity score (negative = 0, weak = 1, moderate = 2, strong = 3) to distinguish between low/loss (≤4) and high (>4) expression of SLC10A1 protein.

Analysis of Genetic Changes and Tumor-Infiltrating Immune Cells

We calculated the abundance ratios of tumor-infiltrating immune cells in the HCC immune microenvironment using the CIBERSORT¹⁷ databases. The R package “maftools” was used to assess changes in genetic variants across various subgroups using the mutation data from TCGA gathered from HCC patients.

Drug Susceptibility Analysis

We analyzed the correlation of 574 drugs participating in late-stage clinical trials and 216 FDA-approved drugs in the CellMiner database,¹⁸ with the expression of each gene in the risk model. An adjusted P value of less than 0.001 and a Pearson correlation coefficient greater than 0.5 for tumor-sensitive drugs were used as cut-off criteria. Subsequently, we also analyzed the differences in the activity Z-score for tumor-sensitive drugs in high- and low-risk samples.

Statistical Analysis

The independent-sample t-test was used to evaluate quantitative variables. R software (Version 4.0.3) was used to analyze the prediction performance of survival outcomes using ROC curve analysis and Kaplan-Meier survival analysis. The association between a prognostic classifier and survival outcomes, as well as other clinical variables, was investigated using a Cox proportional model. When the P-value was less than 0.05, the results were considered statistically significant.

Results

Molecular Subtype Identification Based on Ferroptosis-Associated Genes

By using the NMF algorithm, the optimal number of clusters was identified as two based on cophenetic, dispersion, and profile (Figures 1A and S1). Then, 342 hCC samples were separated into two groups with significant survival differences (Figure 1B). We identified 431 ferroptosis-related DEGs between the two clusters (Figure 1C). Interestingly, the top 10 gene mutation rates were completely different between our two clusters (Figure 1D), and there were also significant differences in the infiltration levels of various immune cells (Figure 1E). Overall, we identified two clusters with different survival outcomes, gene mutation, and infiltration levels of immune cells based on these ferroptosis-associated genes.

Figure 1 Molecular subtype identification based on ferroptosis-associated genes. (A) By using the NMF algorithm, the optimal number of clusters was identified as two. (B) HCC samples were separated into two groups with significant survival differences. (C) 431 ferroptosis-related DEGs between the two clusters. (D) The top 10 gene mutation rates in the two clusters. (E) Significant differences in the infiltration levels of various immune cells. ns, not significant; *p < 0.05; **p <0.01; ***p < 0.001.

Molecular Subtype Identification Based on Angiogenesis-Associated Genes

By using the NMF algorithm, the optimal number of clusters was identified as three based on cophenetic, dispersion, and profile (Figure S2). However, there was no difference in patient survival between the three clusters at this time. Therefore, we re-defined the number of clusters as two (Figure 2A). Then, 342 hCC samples were separated into two groups with significant survival differences (Figure 2B). We identified 1284 angiogenesis-related DEGs between the two clusters (Figure 2C). Interestingly, the top 10 gene mutation rates were completely different between our two clusters (Figure 2D), and there were also significant differences in the infiltration levels of various immune cells (Figure 2E). Overall, we identified two clusters with different survival outcomes, gene mutation, and infiltration levels of immune cells based on these angiogenesis-associated genes.

Figure 2 Molecular subtype identification based on angiogenesis-associated genes. (A) By using the NMF algorithm, the optimal number of clusters was identified as two. (B) HCC samples were separated into two groups with significant survival differences. (C) 1284 angiogenesis-related DEGs between the two clusters. (D) The top 10 gene mutation rates in the two clusters. (E) Significant differences in the infiltration levels of various immune cells. ns, not significant; *p < 0.05; **p <0.01; ***p < 0.001.

Development of a Prognostic Model Associated with Angiogenesis and Ferroptosis

Among these angiogenesis and ferroptosis-related DEGs, a total of 120 genes were identified as overlapping DEGs (Figure 3A) and were calculated by univariate Cox regression (Figure 3B). The 88 prognostically relevant DEGs were further screened by LASSO-COX, RF, and XGBoost analysis. We screened 15 significant genes from each machine-learning method (Figure 3C–E). Six overlapping genes among these genes were used to construct a prognostic model (Figure 3F). Risk score = (0.0466 × PAEP) + (0.0501 × DPYSL4) + (0.0893 × MSC) + (0.0575 × NQO1) + (0.0945 × CD24) - (0.0668 × SLC10A1).

Figure 3 Hub genes associated with angiogenesis and ferroptosis screening. (A) Overlap of DEGs associated with angiogenesis and ferroptosis. (B) Univariate Cox regression. (C) LASSO-Cox analysis. (D) XGBoost analysis. (E) RF analysis. (F) Venn Plot.

Validation of the Prognostic Model in the TCGA, ICGC, and GSE14520 Cohorts

HCC samples were divided into high-risk and low-risk groups according to the risk score that was determined for each patient (Figure 4A). The Kaplan-Meier survival curves demonstrated that the high-risk group’s survival rate was lower than that of the low-risk group (Figure 4B). According to ROC analysis, the model’s area under the curve (AUC) at 1-, 2-, and 3 years was 0.725, 0.679, and 0.671 (Figure 4C), respectively. Moreover, this model might be used as an independent prognostic factor for HCC patients, according to univariable (HR=2.955, 95% CI 1.829-4.773, P <0.001) and multivariate Cox regression models (HR=2.566, 95% CI 1.544-4.265, P <0.001). Similarly, in the ICGC cohort, HCC samples were also divided into high- and low-risk subgroups (Figure 4D). Patients with a higher risk score had considerably worse survival outcomes, according to a Kaplan-Meier study (Figure 4E). ROC analysis demonstrated that this model had high prognostic performance, with AUCs of 0.696, 0.649, and 0.641 at 1-, 2-, and 3 years, respectively (Figure 4F). Moreover, this model might be used as an independent prognostic factor for HCC patients, according to univariable (HR=4.032, 95% CI 1.657-9.814, P = 0.002) and multivariate Cox regression models (HR=2.829, 95% CI 1.181-6.777, P = 0.019). Finally, HCC samples in the GSE14520 cohort were also divided into high- and low-risk subgroups (Figure 4G). Patients with a higher risk score had considerably worse survival outcomes, according to a Kaplan-Meier study (Figure 4H). ROC analysis demonstrated that this model had high prognostic performance, with AUCs of 0.617, 0.633, and 0.623 at 1-, 2-, and 3 years, respectively (Figure 4I). Unfortunately, this model was not an independent prognostic factor for HCC patients in the GSE14520 cohort according to the univariate model (HR=3.944, 95% CI 1.765-8.816, P=0.000) and the multivariate Cox regression model (HR=2.279, 95% CI 0.988-5.255, P=0.051), which may be related to individual variability of HCC patients. In summary, the prognostic model we constructed performed well in predicting the survival outcome of HCC patients.

Figure 4 Validation of the prognostic model in the TCGA, ICGC, and GSE14520 cohorts. (A) HCC samples were divided into high-risk and low-risk groups in the TCGA cohort. (B) The Kaplan-Meier survival curves in the TCGA cohort. (C) ROC analysis in the TCGA cohort. (D) HCC samples were divided into high-risk and low-risk groups in the ICGC cohort. (E) The Kaplan-Meier survival curves in the ICGC cohort. (F) ROC analysis in the ICGC cohort. (G) HCC samples were divided into high-risk and low-risk groups in the GSE14520 cohort. (H) The Kaplan-Meier survival curves in the GSE14520 cohort. (I) ROC analysis in the GSE14520 cohort.

The Correlation Analysis Between Risk Scores and Clinical Features

We then analyzed the relationship between risk scores and clinical case characteristics of HCC patients. We found that risk score was associated with later grade stage, later TNM stage, later T stage, tumor recurrence, cirrhosis, higher AFP values, microvascular infiltration, and HBV infection, independently of age and gender (Figure S3). These suggest that risk scores are strongly associated with disease progression in HCC patients.

Analysis of Immune Infiltrating Cells

We computed the four ESTIMATE indices to provide an extensive overview of tumor immunity in HCC. ESTIMATE is an algorithm implemented with the “ESTIMATE” package that uses the unique properties of a sample’s transcriptional profile to infer the content of tumor cells and the different infiltrating normal cells. A higher ESTIMATE score suggests that the sample has a higher proportion of infiltrating stromal and immune cells and a lower purity of tumor cells. We discovered that while tumor purity was lower in the high-score group, stroma, immunity, and ESTIMATE scores were all higher in the high-score group (Figure 5A). Furthermore, variations in the percentage of partly immunized cells between the two groups were revealed using CIBERSORT analysis (Figure 5B). Additionally, CIBERSORT is an inverse convolution analysis technique based on the linear support vector regression theory that employs roughly 22 distinct immune cell classes. Specifically, compared to samples from the low-score group, samples from the high-score group had a higher proportion of memory B cell, resting memory CD4 T cell, helper follicular T cell, Tregs, activated NK cell, M0 Macrophage, M1 Macrophage, M2 Macrophage, resting myeloid dendritic cell, and neutrophil. Finally, we found that there was a tendency for the high-score group to have higher levels of major histocompatibility complexes (MHCs) expression compared to the low-score group (Figure 5C). When combined, the immunological infiltration of the high-score group was often stronger than that of the low-score group.

Figure 5 Analysis of immune infiltrating cells. (A) While tumor purity was lower in the high-score group, stroma, immunity, and ESTIMATE scores were all higher in the high-score group. (B) The CIBERSORT analysis. (C) There was a tendency for the high-score group to have higher levels of MHCs expression compared to the low-score group. ns, not significant; *p < 0.05; **p <0.01; ***p < 0.001.

Immunotherapy Sensitivity and Genetic Alterations Analysis

We examined the possible sensitivity to immunotherapy by comparing the expression of numerous genes involved in immunomodulation, such as those related to PD1, CTLA4, and other agonists or antagonists of T-cell activation, between the two groups. The majority of the altered genes were found to be elevated in the high-score group, as shown in Figure 6A–C. We also assessed the gene mutation profile and illustrated the topographies of the two groupings (Figure 6D) as gene mutations could impact immunotherapy sensitivity. Tumor mutation burden (TMB), however, did not show a discernible variation between the two groups (Figure 6E). According to the aforementioned findings, the high-score group could be more sensitive to immunotherapy than the low-score group.

Figure 6 Immunotherapy sensitivity and genetic alterations analysis. (A) PD1 and CTLA4-related genes. (B) Agonists or antagonists of T-cell activation-related genes. (C) Antagonists of T-cell activation-related genes. (D) The gene mutation profile of the two groups. (E) TMB did not show a discernible variation between the two groups. ns, not significant; *p < 0.05; **p <0.01; ***p < 0.001.

Function Enrichment Analysis

To confirm the mechanism by which this prognostic model is related to tumor immunity, we performed the differential analysis of gene expression between high- and low-score groups. At the P-value was less than 0.05 and |logFC| > 1, we screened a total of 1190 upregulated and 183 downregulated genes (Figure 7A). Subsequently, we performed GO and KEGG enrichment analyses on these DEGs. GO enrichment analysis showed that these genes were enriched in immunity related to cellular lipid fatty acid metabolism (Figure 7B). In addition, KEGG analysis showed that these genes were enriched in ECM-receptor interaction, PPAR signaling pathway, and drug metabolism-related pathways (Figure 7C).

Figure 7 Function enrichment analysis. (A) 1190 upregulated and 183 downregulated genes between high- and low-score groups. (B) GO enrichment analysis. (C) KEGG analysis.

Establishment of a Nomogram Model

A nomogram model was created in the TCGA dataset to investigate the coefficient prediction effectiveness of this classifier, and the results suggested that the nomogram might help us give a quantitative way for properly predicting the 1-, 2-, and 3-year survival rates (Figure 8A). The calibration curves showed good agreement between the anticipated and actual probability of 1-, 2-, and 3-year survival rates no matter in the TCGA cohort (Figure 8B), ICGC cohort (Figure 8C), or GSE114520 cohort (Figure 8D).

Figure 8 Establishment of a nomogram model. (A) A nomogram model was created in the TCGA dataset. (B) The calibration curves in the TCGA cohort. (C) The calibration curves in the ICGC cohort. (D) The calibration curves in the GSE14520 cohort.

Tumor-Sensitive Drug Screening

Nine of the 790 drugs were found to be strongly associated with NQO1 and CD24 in the CellMiner database based on drug sensitivity studies (Figure 9A). We also compared the difference in the activity Z-score between these 9 drugs in the high- and low-risk groups and found that DACARBAZINE was more sensitive in the high-scoring group, whereas Sapitinib, Poziotinib, AZD-3759, and EGF-816 were more sensitive in the low-scoring group (Figure 9B).

Figure 9 Tumor-sensitive drug screening. (A) Nine of the 707 drugs were found to be strongly associated with NQO1 and CD24 in the CellMiner database. (B) The difference in the activity Z-score of the 9 drugs between the high- and low-risk groups. ns, not significant; **p <0.01.

Verification of the Expression of the Genes in the Prognostic Model in Clinical Samples

In the TCGA cohort, lower levels of NQO1, CD24, and PAEP and higher levels of SLC10A1 were found in normal tissues compared with HCC samples (Figure 10A); whereas in the ICGC cohort, lower levels of NQO1, CD24, and higher levels of SLC10A1 were found in normal tissues compared with HCC samples (Figure 10B). We also explored the association of NQO1, CD24, and SLC10A1 with survival outcomes in HCC patients. We also explored the association of NQO1, CD24, and SLC10A1 with survival in HCC patients. We found that in the TCGA cohort, low expression of SLC10A1 and high expression of NQO1 were associated with worse prognosis in HCC patients, while CD24 did not affect the prognosis of HCC patients (Figure 10C); whereas in the ICGC cohort, we only found that low expression of SLC10A1 was associated with worse prognosis in HCC patients, while CD24 and NQO1 did not affect the prognosis of HCC patients (Figure 10D). These suggest that SLC10A1 is the most critical gene in the prognostic model we constructed. According to the Human Protein Atlas database (HPA),¹⁹ the protein level of SLC10A1 was differentially expressed in tumor and normal tissues (Figure 10E). To validate the above bioinformatics analysis results, we validated the SLC10A1 expression in 50 clinical samples by qRT-PCR and IHC. SLC10A1 mRNA levels were lower in HCC tissues compared to normal samples (Figure 10F). We also found that the SLC10A1 protein showed differential expression in HCC tissues (Figure 10G), and the positivity rate in HCC tissues was significantly lower than that in normal tissues (Figure 10H).

Figure 10 Verification of the expression of the genes in the prognostic model in clinical samples. (A) In the TCGA cohort, lower levels of NQO1, CD24, and PAEP and higher levels of SLC10A1 were found in normal tissues compared with HCC samples. (B) In the ICGC cohort, lower levels of NQO1, CD24, and higher levels of SLC10A1 were found in normal tissues compared with HCC samples. (C) In the TCGA cohort, low expression of SLC10A1 and high expression of NQO1 were associated with worse prognosis in HCC patients, while CD24 did not affect the prognosis of HCC patients. (D) In the ICGC cohort, we only found that low expression of SLC10A1 was associated with a worse prognosis in HCC patients, while CD24 and NQO1 did not affect the prognosis of HCC patients. (E) SLC10A1 expression in the Human Protein Atlas database (HPA). (F) SLC10A1 mRNA levels were lower in HCC tissues compared to normal samples. (G and H) The positivity rate of SLC10A1 protein in HCC tissues was significantly lower than that in normal tissues. *p < 0.05; ***p < 0.001.

Discussion

Ferroptosis plays an essential role in the pathophysiology of cancer.² Considering that the extensively mutated p53 gene in tumors is strongly associated with susceptibility to ferroptosis,²⁰ ferroptosis, which could effectively inhibit tumor progression, might be an underlying mechanism that resists tumorigenesis. Furthermore, direct crosstalk between ferroptosis and antitumor immunity has been demonstrated in a previous study. Enhanced induction of ferroptosis by CD8+ T cells could facilitate the anti-tumor effect of immunotherapy,⁸ such as PD-L1 inhibitor treatment. Meanwhile, tumor angiogenesis, which enhances invasion, progression, and drug resistance, is one of the features of HCC progression, and many treatments for HCC revolve around the inhibition of angiogenesis and vascular normalization.²¹ Therefore, antiangiogenic therapy combined with ferroptosis induction has the potential to improve the prognosis of HCC patients.

To date, with the widespread use of RNA-seq technology, machine learning methods have become increasingly popular among researchers in the medical field, shining in feature extraction, diagnostic classification, or patient survival analysis of tumors or non-tumors.^22–24 Clinicians may employ big data analysis and machine learning to evaluate vast amounts of clinical data, make precise prognostic predictions for their patients, and help identify workable treatment strategies and symptomatic treatments.²⁵ Machine learning can be used to comprehensively analyze 24 clinical variables commonly found in critical care settings to accurately predict the need for vasopressor administration to patients and successfully apply it to clinical care.²⁶ Additionally, machine learning may be utilized as a decision-support tool for diagnosing illnesses, identifying candidates for surgery, and forecasting possible consequences after surgery.²⁷ In cancers, machine learning is not only applied to the diagnosis of tumors but can also be used to predict the efficacy of immunotherapy for patients, which greatly improves the survival outcomes of tumor patients.^28–30 Our study explored HCC transcriptomic data in detail through a series of multiple machine learning methods including NMF clustering, LASSO-COX, RF, and XGBOOST in an attempt to identify a novel ferroptosis- and angiogenesis-related gene signature to screen high-risk HCC patients. Finally, we successfully constructed a six-gene prognostic model, which not only can classify HCC patients into two subgroups with significant survival differences, mutation profiles, and immune cell infiltration status but also can predict the efficiency of HCC receiving immunotherapy. Overall, through the use of powerful algorithms, machine learning can optimize the subjectivity and human error that traditional screening procedures are prone to, identify key biomarkers associated with cancer, assist in assigning personalized treatment choices to improve patient prognosis, and provide important insights into the underlying course of the disease.

Increasing evidence suggests that the tumor microenvironment (TME) plays a major role in the occurrence and progression of HCC.³¹ TME can affect vascular invasion, early recurrence, overall survival outcomes, and even the efficacy of immunotherapy in HCC. In this study, we used the CIBERSORT algorithms to evaluate immune cell infiltration in HCC patients and found that patients in the high-score group had a higher proportion of various immune cell infiltration, including memory B cell, resting memory CD4 T cell, helper follicular T cell, Tregs, activated NK cell, M0 Macrophage, M1 Macrophage, M2 Macrophage, resting myeloid dendritic cell, and neutrophil, suggesting that these immune cells may be involved in immune evasion. Of course, bioinformatics analysis alone is not enough, and single-cell sequencing techniques can be used to confirm this result in the future.

The high-risk and low-risk groups showed a variety of medications with varying drug sensitivity levels; those with a significant difference were Dacarbazine, Sapitinib, Poziotinib, AZD-3759, and EGF-816. This suggested that choosing alternative medications for individuals with various risk ratings would enhance the result. Although all five drugs can treat tumors,^32–35 only Dacarbazine has been reported to be associated with HCC,³⁶ but its mechanism is unknown. To verify this finding, future enrollment of HCC patients on these medications is necessary.

Most of the genes that make up the risk model are strongly associated with tumors. PAEP, a glycoprotein with a molecular weight of 28 kDa belonging to the lipocalin superfamily, is over-expressed in melanoma, and PAEP small interfering RNA significantly inhibits cell migration and cell invasion ability of melanoma cells.³⁷ PAEP is also an independent significant risk factor for death in patients with clear cell renal cell carcinoma.³⁸ PAEP may also serve as a novel therapeutic target to weaken the defenses of the NSCLC immune system by targeting glycoproteins.³⁹ DPYSL4 is a collagen response-mediated protein that can regulate oxidative phosphorylation and cellular energy supply through association with the mitochondrial supercomplex thereby participating in the regulation of cancer invasion and development.⁴⁰ Unfortunately, no studies have yet demonstrated a link between PAEP and HCC. In addition to being a component of the body’s antioxidant defense system, which shields cells from carcinogenesis and mutation, NQO1 may be a crucial biomarker for the early detection, diagnosis, and treatment of gastrointestinal malignancies.⁴¹ NQO1 can regulate ros-induced apoptosis to mediate lenvatinib resistance in HCC.⁴² Inhibiting macrophage phagocytosis and transmitting immunosuppressive signals, CD24 is a new therapeutic target for tumor treatment. It interacts with the inhibitory receptor Siglec-10 on tumor-associated macrophages. In HCC, oxaliplatin and CD24 expression suppression dramatically decreased tumor invasion, migration, and proliferation while increasing apoptosis.⁴³ It has been suggested that the gene muscledin, which controls lymphocyte growth, might be an immunotherapeutic target for HCC patients.⁴⁴ Lower expression of SLC10A1 is linked to an earlier recurrence and a worse prognosis for HCC patients.⁴⁵ By inhibiting aerobic glycolysis through the regulation of SLC10A1, LINC00659 prevents the malignant development of HCC.⁴⁶ In the present study, we found that SLC10A1 was involved in the development and progression of HCC as an angiogenesis and ferroptosis-associated gene, and its specific molecular mechanism is worthy of further investigation in future work.

Our research has certain limitations, to be sure. The biological processes underlying the crosstalk between angiogenesis and ferroptosis in HCC must be clarified by biological experimental research. Furthermore, the diversity and individual variability of HCC patients may alter the expression of this signature, necessitating future prospective multicenter randomized controlled trials to assess this signature. Finally, additional in vivo and in vitro experiments are needed to better investigate the particular mechanisms of the angiogenesis-ferroptosis-related signature in HCC.

Conclusions

In conclusion, we established and validated an individualized angiogenesis-ferroptosis-related prognostic signature in HCC patients, which could potentially guide the selection of optimal individualized immunotherapy for HCC patients in future clinical practice.