A six-gene prognostic signature for aflatoxin B1-associated hepatocellular carcinoma identified through integrated bioinformatics and network toxicology

Introduction

Liver cancer, as one of the most prevalent and lethal malignant tumors in the world, poses a severe threat to human life and health. Its pathogenesis is complex, involving the interplay of genetic, environmental, and other factors. Despite global regulatory efforts to reduce harmful exposure including maximum residue limits and improved agricultural practices, hepatocellular carcinoma (HCC) remains a major public health challenge with increasing incidence rates in many regions (1). The persistence of HCC problem highlights the need for better understanding of molecular mechanisms and improved risk prediction tools. Investigating the etiology of liver cancer, identifying effective therapeutic targets, and constructing precise prognostic models have remained critical yet challenging areas in oncology research.

Aflatoxin B1 (AFB1), a potent carcinogen widely found in mold-contaminated food, is closely associated with the development of liver cancer (2). AFB1, aflatoxin produced by Aspergillus flavus and Aspergillus parasiticus, is recognized as one of the most potent naturally occurring hepatocarcinogens. Its carcinogenic potential was first identified in the 1960s through studies linking groundnut consumption to increased liver cancer incidence in African populations, where improper storage conditions led to widespread fungal contamination. The International Agency for Research on Cancer (IARC) classified AFB1 as a group 1 carcinogen in 1993, based on overwhelming evidence from epidemiological studies in Chinese mainland, Taiwan, and sub-Saharan Africa showing strong dose-response relationships between dietary AFB1 exposure and HCC risk, particularly in populations with concurrent hepatitis B virus infection where synergistic effects dramatically increase cancer risk (2,3).

Extensive epidemiological studies indicate that long-term exposure to AFB1 significantly increases the risk of liver cancer. However, the specific molecular mechanisms by which AFB1 induces liver cancer are not yet fully understood, and the key genes, signaling pathways, and their associations with clinical features require further investigation (4).

With the advancement of bioinformatics and omics technologies, the Cancer Genome Atlas (TCGA) database has provided a vast resource for cancer research (5). TCGA integrates multi-omics data, including genomics, transcriptomics, and proteomics, along with clinical information, enabling a comprehensive and systematic analysis of the molecular characteristics and pathogenesis of cancer. Concurrently, network toxicology, an emerging interdisciplinary field, combines theories and methods from toxicology, systems biology, and bioinformatics. By constructing compound-target-disease networks, it reveals the toxicological mechanisms of chemicals at a systems level (6).

This study integrates TCGA data with network toxicology approaches to explore the intrinsic link between AFB1 and liver cancer. We extracted sequencing and clinical data of 377 patients from TCGA, systematically identified key targets of AFB1 in liver cancer, and constructed a robust prognostic model to reveal the potential molecular mechanisms by which AFB1 promotes liver cancer progression, aiming to provide a valuable prognostic prediction tool for AFB1-induced liver cancer. The findings help to reveal theoretical foundations and novel insights for the precise diagnosis, personalized treatment of liver cancer, and prevention of AFB1-associated liver cancer. We present this article in accordance with the TRIPOD reporting checklist (available at https://tgh.amegroups.com/article/view/10.21037/tgh-2026-0053/rc).

Methods

Core genes and clinical data of liver cancer

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Gene expression profile data (RNA-seq data) and clinical data (including patient age, gender, tumor stage, survival status, etc.) of liver cancer “LIHC” were obtained from the TCGA database (7). The data were downloaded on 1 September 2025. Differential analysis was performed using the limma package of R software to screen differentially expressed genes, lncRNAs, and mRNAs closely related to the occurrence and development of liver cancer, and the results were visualized graphically. Figure 1 provides a schematic representation of the research conceptualization and technical implementation strategy employed in this study.

Figure 1 Flow chart of integrating bioinformatics and network toxicology to identify a prognostic signature for AFB1-associated HCC. AFB1, aflatoxin B1; AIC, Akaike Information Criterion; GO, Gene Ontology; HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; LIHC, liver hepatocellular carcinoma; PPI, protein-protein interaction; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Compound target screening

The SMILES name and 3D protein structure of AFB1 were obtained from the PubChem database (8). The SMILES name was input into the SwissTarget database (http://www.swisstargetprediction.ch/) (9), the ChEMBL database (10), and the SEA database (Similarity Ensemble Approach) (11) for toxicant target prediction. A confidence threshold (≥95%) was applied for screening, and after merging and deduplication, the relevant target proteins of AFB1 were identified.

Disease target screening

The keywords “liver cancer” and “hepatic carcinoma” were used to search in the GeneCard (12), Online Mendelian Inheritance in Man (OMIM) (13), and Comparative Toxicogenomics Database (CTD) (14) databases to identify targets associated with liver cancer. After merging and deduplication, the relevant target proteins of liver cancer were obtained.

Venn intersection genes

The differentially expressed genes obtained from TCGA differential analysis, the predicted target genes of AFB1, and the liver cancer-related genes were visualized using the ggvenn package in R software to generate a Venn diagram, revealing the genes through which AFB1 influences liver cancer.

Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes pathway analysis

To investigate the biological roles of the potential targets of AFB1 in liver cancer pathogenesis, GO analysis and KEGG pathway enrichment analysis were performed. The GO analysis covered biological process (BP), cellular component (CC), and molecular function (MF) assessments to elucidate the core biological functions of the relevant targets. Through KEGG pathway enrichment analysis, significant pathways associated with the potential targets of AFB1 in liver cancer were screened. An adjusted P value <0.05 was used as the criterion to identify the major target pathways. The clusterProfiler package in R software was employed for GO functional enrichment and KEGG pathway analysis of the intersection genes (15), and the results were visualized.

Protein-protein interaction (PPI) network construction and core target screening

The intersection targets of AFB1 and liver cancer were uploaded to the STRING platform (16), with the species limited to “Homo sapiens” and the confidence set to high (0.4). Disconnected nodes in the network were hidden, resulting in a PPI network diagram of the intersection genes. The file was imported into Cytoscape 3.9.1 software (17) in “TSV” format, and the Hubba algorithm in the CytoHubba plugin (18) was used to screen the top 10/20 genes with the highest degree values as core genes.

Prognosis-related gene screening

The core genes from the PPI network were correlated with survival data and clinical data obtained from the TCGA database. The Cox proportional hazards regression model was used to screen genes associated with the prognosis of liver cancer patients. The survival package in R software was utilized for this analysis (19).

Construction of disease prognostic model

Based on the screened prognosis-related genes, a multivariate Cox regression model was employed to construct a prognostic model for liver cancer. The rms package in R software was used to construct and optimize the model (19).

SHapley Additive exPlanations (SHAP) analysis

SHAP analysis was performed on the constructed prognostic model to interpret its predictions. The SHAP values for each feature (prognosis-related gene) were calculated, and visualizations (such as SHAP waterfall plots, summary plots, and dependence plots) were generated to demonstrate the impact of each gene on the model’s predictions.

Model accuracy validation

Survival analysis: the survminer and survival packages in R software were used to plot survival curves (Kaplan-Meier curves) to compare survival differences among patient groups stratified by the prognostic model. Receiver operating characteristic (ROC) curve plotting: the survminer package was used to plot ROC curves and calculate the AUC to evaluate the predictive performance of the model. An AUC value closer to 1 indicated better predictive accuracy. Risk curve/heatmap plotting: based on the predicted risk scores, risk curves were plotted to reveal the changes in patient risk at different time points.

Independent prognostic model analysis

Possible clinical covariates (such as age, gender, tumor stage, etc.) were incorporated into the Cox regression model to construct a multivariate Cox regression model. The independence of the prognosis-related genes was analyzed after adjusting for these clinical covariates (19).

Molecular docking

Molecular docking analysis was conducted to evaluate the binding affinity and interaction patterns between liver cancer and its targets. The protein structures encoded by the prognosis-related genes CCNB1, CCNB2, CHEK1, MMP1, ESR, and TTK were obtained from the Protein Data Bank (PDB) (20). The online platform https://cadd.labshare.cn/cb-dock2/index.php was used to perform molecular docking with the receptor compound AFB1 and visualize the results (21). The binding affinity and binding mode of AFB1 with the proteins were assessed. Lower binding free energy indicated a more stable docking conformation.

Statistical analysis

This article employed the following statistical methods: analysis of variance, Cox regression, SHAP analysis, survival analysis, and ROC curve analysis, using the R software version 4.3.2.

Results

Differential gene expression of liver cancer

Gene expression profile data (RNA-seq data) and corresponding clinical data for liver cancer [liver hepatocellular carcinoma (LIHC)] were successfully obtained from the TCGA database. The clinical data included information from 377 patients, such as age, gender (255 males, 122 females), tumor stage (175 cases in stage I, 87 in stage II, 86 in stage III, 5 in stage IV, and 24 unclassified), and survival status (245 alive, 132 deceased). Differential analysis of gene expression data between liver cancer tissues and normal liver tissues was performed using the limma package of R software (22). The screening criteria were set at |log₂FC| >1 and adjusted P value <0.05, resulting in the identification of 4,506 differentially expressed mRNAs and 2,806 differentially expressed lncRNAs. The pheatmap package of R software was used to generate heatmaps and volcano plots of the differential genes (as shown in Figure 2A), illustrating the expression pattern differences between liver cancer and normal liver tissues. The results clearly distinguished the two sample types, confirming that the screened differential genes effectively reflect the gene expression characteristics of liver cancer tissues.

Figure 2 Identification of key targets for AFB1-related HCC. (A) Volcano plot and heatmap showing the DEGs in HCC from TCGA. (B) Venn diagrams displaying the number of targets obtained from AFB1 prediction and the HCC disease database. (C) Venn diagram identifying 35 overlapping genes among DEGs, AFB1 targets, and HCC-related genes. (D) GO functional enrichment of the 35 overlapping genes. (E) KEGG pathway analysis of the 35 overlapping genes. (F) Protein-protein interaction (PPI) network of the overlapping genes, highlighting the top 19 hub genes. (G) GO functional enrichment and KEGG pathway analysis of the top 19 hub genes are highlighted. AFB1, aflatoxin B1; DEGs, differentially expressed genes; GO, Gene Ontology; HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; TCGA, The Cancer Genome Atlas.

AFB1 toxicity target screening

The SMILES structure (COC1=C2C3=C(C(=O)CC3)C(=O)OC2=C4[C@@H]5C=CO[C@@H]5OC4=C1) and 3D protein structure of AFB1 were retrieved from the PubChem database. The SMILES structure was input into the SwissTargetPrediction, ChEMBL, and SEA databases for toxicity target prediction. SwissTargetPrediction predicted 104 potential toxicity targets with a probability >0. The ChEMBL database screened 89 potential toxicity targets with a probability >0. The SEA database predicted 1 potential toxicity target with a probability >0. After merging and deduplicating the potential toxicity targets predicted by the three databases, a total of 186 AFB1-related potential toxicity targets were ultimately obtained (as shown in Figure 2B, 1).

Liver cancer disease target retrieval results

Based on the previous discussion, the screening process involved searching for liver cancer-related targets in the GeneCard, OMIM, and CTD 14 databases using the keywords liver cancer and hepatic carcinoma. The results were as follows: GeneCard: 4,662 targets (Relevance score >10), OMIM: 205 disease-associated genes, CTD: 731 targets linked to liver cancer development. After merging and deduplicating the data, a total of 5,037 liver cancer-related targets were identified (as shown in Figure 2B, 2).

Venn intersection gene analysis results

Using the R package ggvenn, a Venn diagram analysis was performed on three gene sets: 4,506 differentially expressed mRNA genes (screened from the TCGA database), 186 AFB1-related toxicity targets, 5,037 liver cancer-related targets (as previously identified). The results (Figure 2C) revealed 35 overlapping genes at the intersection of all three datasets. These 35 genes represented key candidate genes through which AFB1 might exert its carcinogenic effects on liver cancer, providing critical targets for further mechanistic studies on AFB1-induced hepatocarcinogenesis.

GO functional enrichment analysis of intersection genes

Using the clusterProfiler package of R software, GO functional enrichment analysis was performed on the 35 intersection genes identified above, with a significance threshold of corrected P value <0.05. BP: the genes were significantly enriched in processes such as protein autophosphorylation (GO:0046777, P=9.8×10⁻⁷), nuclear division (GO:0000280, P=2.4×10⁻⁶), and mitotic cell cycle phase transition (GO:0044772, P=2.4×10⁻⁶). CC: enrichment was observed in cyclin-dependent protein kinase holoenzyme complex (GO:0000307, P=6.9×10⁻⁸), serine/threonine protein kinase complex (GO:1902554, P=6.2×10⁻⁶), and protein kinase complex (GO:1902911, P=7.3×10⁻⁶). MF: the genes were primarily enriched in protein tyrosine kinase activity (GO:0004713, P=5.5×10⁻⁷), protein serine/threonine/tyrosine kinase activity (GO:0004712, P=5.5×10⁻⁷), and protein serine kinase activity (GO:0016301, P=5.5×10⁻⁷) (Figure 2D). These results suggest that the 35 intersection genes are critically involved in cell cycle regulation and kinase-mediated signaling pathways, which may be key mechanisms underlying AFB1-induced hepatocarcinogenesis.

KEGG pathway analysis of intersection genes

KEGG pathway enrichment analysis was performed on the former 35 intersection genes, revealing 58 significantly enriched pathways (corrected P value <0.05). The top pathways included: Cell cycle (hsa04110, P=2.2×10⁻⁹), p53 signaling pathway (hsa04115, P=1.9×10⁻⁷), Cellular senescence (hsa04218, P=9.1×10⁻⁷) (Figure 2E). These pathways were closely associated with cell growth, proliferation, and apoptosis regulation, suggesting that AFB1 may induce liver cancer by disrupting the normal function of these pathways.

Results of PPI network construction and core target screening

The former 35 intersection genes were uploaded to the STRING platform, with the species restricted to Homo sapiens and a confidence score threshold set to 0.4. After hiding disconnected nodes, the resulting PPI network (Figure 2F, 1) comprised 35 nodes and 131 edges, with an average node degree of 7.49 (a measure of node importance).

The PPI network file (in TSV format) was imported into Cytoscape 3.9.1, and the Hubba algorithm (via the CytoHubba plugin) was used to calculate node degree values. This identified 19 core genes with the highest degree values, ranked as follows: SRC (degree =21), MMP9 (degree =19), ESR1 (degree =16), AURKA (degree =14), CDK1 (degree =14), CCNA2 (degree =13), CCNB1 (degree =13), CCNE1 (degree =13), AURKB (degree =12), MCL1 (degree =12), CHEK1 (degree =11),CCNE2 (degree =10), TTK (degree =9), PDGFRB (degree =9), CCNB2 (degree =9), PDGFRA (degree =8), MMP3 (degree =8), CXCR2 (degree =7), AKT1 (degree =6), MMP1 (degree =6) (Figure 2F, 2,3). These core genes (e.g., SRC, MMP9, ESR1) are central to the PPI network and are likely critical in mediating the oncogenic effects of AFB1 on liver cancer, as they are enriched in pathways related to cell cycle regulation, proliferation, and apoptosis (as previously discussed in the KEGG analysis).

GO functional enrichment and KEGG pathway analysis results of core genes

GO functional enrichment analysis

The 19 core genes identified from the PPI network underwent GO functional enrichment analysis (23). The results revealed significant enrichments in the following categories: BP: nuclear division (GO:0000280, P=4.0×10⁻¹²), mitotic cell cycle phase transition (GO:0044772, P=5.6×10⁻¹²), cell cycle G2/M phase transition (GO:0044839, P=7.5×10⁻¹⁰). CC: cyclin-dependent protein kinase holoenzyme complex (GO:0000307, P=1.4×10⁻¹¹), serine/threonine protein kinase complex (GO:1902554, P=2.6×10⁻⁹), protein kinase complex (GO:1902911, P=4.6×10⁻⁹). MF: cyclin-dependent protein serine/threonine kinase regulator activity (GO:0016538, P=9.0×10⁻¹⁰), phospholipase activator activity (GO:0016004, P=9.9×10⁻⁷), lipase activator activity (GO:0060229, P=1.8×10⁻⁶).

KEGG pathway enrichment analysis

KEGG pathway analysis of the core genes identified 46 significantly enriched pathways (corrected P value <0.05) (24). In addition to the pathways enriched in the intersection genes (e.g., cell cycle, p53 signaling pathway), the core genes were also significantly enriched in: cellular senescence (hsa04218, P=1.3×10⁻⁸), prostate cancer (hsa05215, P=4.2×10⁻⁸) (Figure 2G). These results further suggest that the core genes play a critical role in the cellular abnormalities induced by AFB1 during hepatocarcinogenesis.

Univariate cox regression

The 19 core genes were jointly analyzed with the survival data and clinical data of liver cancer patients from the TCGA database. The R software survival package was used to perform univariate cox regression. The results showed that a total of 13 core genes significantly affected the prognosis of liver cancer patients (P<0.05). Among them, the genes whose high expression was associated with poor prognosis included MMP9, CDK1, AURKA, CCNB1, CCNA2, AURKB, CHEK1, CCNE2, CCNB2, TTK, MMP3, and MMP1, while the gene whose high expression was associated with good prognosis was ESR1 (Figure 3A).

Figure 3 Construction and validation of the prognostic prediction model. (A) Forest plot of univariate Cox regression analysis for 13 prognosis-related hub genes. (B) Coefficients of the 6 genes included in the final multivariate Cox proportional hazards model. (C) Kaplan-Meier survival curves for high-risk and low-risk groups stratified by the median risk score. (D) Time-dependent receiver operating characteristic (ROC) curves at 1, 3, and 5 years. (E) Distribution of risk scores and patient survival status. (F) Heatmap of the expression levels of the 6 prognostic genes in high- and low-risk groups. AUC, area under the curve; CI, confidence interval.

Multivariate cox regression

Based on the 13 prognosis-related genes identified from univariate cox regression, a multivariate cox regression model was constructed using the rms package in R. The final predictive model incorporated 6 genes, with the following risk score formula:

Risk Score = 0.37 × CCNB1 expression + 0.28 × CHEK1 expression + 0.14 × MMP1 expression + 0.31 × TTK expression − 0.68 × CCNB2 expression − 0.09 × ESR1 expression (Figure 3B).

Model accuracy validation results

Survival analysis results

Based on the risk scores calculated by the above prognostic model, a cohort of 377 patients with HCC, comprising 343 cases with available survival data from the TCGA database, was stratified into a high-risk group (risk score > median, n=163) and a low-risk group (risk score ≤ median, n=180).Using the R packages survminer and survival, Kaplan-Meier survival curves were plotted (Figure 3C). The results showed that the median survival time of patients in the high-risk group was significantly shorter than that of the low-risk group, with a statistically significant difference in survival rates between the two groups (P<0.001). This indicates that the prognostic model can effectively distinguish liver cancer patients with different prognosis risks.

ROC curve results

Using the survminer R package, ROC curves were plotted to evaluate the prognostic model’s predictive performance at 1-year, 3-year, and 5-year time points (Figure 3D). The results showed that the AUC values were 0.791 (1-year), 0.745 (3-year), and 0.656 (5-year). Notably, the 1-year and 3-year AUC values were both above 0.74, indicating strong predictive accuracy for these time frames. The 5-year AUC value (0.656) suggests a moderate long-term predictive capability, which aligns with the model’s ability to stratify patients by risk (as shown in the Kaplan-Meier analysis).

Risk curve and heatmap results

The risk curve (Figure 3E) was plotted based on patient risk scores, showing that over time, the high-risk group had a significantly higher mortality risk than the low-risk group, with earlier and more frequent death events. Additionally, a heatmap of the 6 prognosis-related genes in the high- vs. low-risk groups (Figure 3F) revealed: Higher expression of poor-prognosis genes (e.g., CCNB1, CCNB2, CHEK1) in the high-risk group. Lower expression of favorable-prognosis genes (e.g., ESR1) in the high-risk group. These findings align with previous analyses, further validating the model’s ability to stratify patients by risk.

SHAP analysis results

SHAP analysis was performed on the constructed liver cancer prognostic model using R software (25). The SHAP summary plot (Figure 4A) revealed that CCNB1, CCNB2, and TTK had the greatest contributions to the model’s predictions. Specifically, the SHAP values for CCNB1 and TTK were positive, indicating that their high expression increases the risk of liver cancer in patients. In contrast, the SHAP value for CCNB2 was negative, suggesting that its high expression may reduce the risk of liver cancer, consistent with the Cox regression results. A randomly selected high-risk patient was used to generate a SHAP waterfall plot (Figure 3B), which visually illustrates the contribution of each gene to the individual risk prediction. Additionally, the SHAP dependence plot (Figure 4B) further demonstrated a positive correlation between CCNB1 expression and SHAP values. As CCNB1 expression increased, the SHAP value gradually rose, leading to an elevated risk score, further validating the critical role of CCNB1 in the prognostic model.

Figure 4 Interpretation of the prognostic prediction model and independent predictive analysis. (A) SHAP summary plot showing the contribution of each gene to the model output. (B) SHAP dependence plots illustrating the relationship between expression levels and SHAP values (impact on risk score). (C) Forest plot of the multivariate Cox regression model incorporating clinical covariates and the 6-gene signature, demonstrating its independent prognostic value. AIC, Akaike Information Criterion; CI, confidence interval; SHAP, SHapley Additive exPlanations.

Independent prognostic prediction model analysis results

A multivariate Cox regression model incorporating clinical covariates (age, sex, tumor stage) and 6 prognosis-related genes was constructed. The likelihood ratio test (Log-Rank) showed high global significance (χ²=7.8636e−07, P=7.8636e−07), indicating at least one clinical variable significantly affected survival outcomes. Akaike Information Criterion (AIC) comparison revealed the model with prognostic genes (AIC =1050.71) had significantly better fit than the clinical-only model (AIC =1910.9). The model’s Concordance Index (C-index =0.71) demonstrated moderate-to-high predictive discrimination capability (Figure 4C).

Molecular docking results

Protein structures of prognosis-related genes (CCNB1, CCNB2, CHEK1, MMP1, ESR1, TTK) were retrieved from RCSB PDB. AFB1 was docked with these proteins using CB-Dock2. Results showed strong binding affinity (binding free energy <−5.0 kcal/mol) between AFB1 and all six proteins. The lowest binding free energy was observed with CHEK1 binding to AFB1 (−9.4 kcal/mol). Visualization of the docking complex revealed that AFB1 primarily binds to CHEK1’s active site via hydrogen bonds (Asn132, Lys133) and hydrophobic interactions (Leu129, Val130), suggesting potential functional interference in HCC development (Figure 5).

Figure 5 Molecular docking validation of AFB1 with proteins encoded by genes in the prognostic prediction model. Three-dimensional structures of molecular docking between AFB1 and CCNB1, CCNB2, CHEK1, MMP1, ESR1, and TTK proteins. AFB1, aflatoxin B1.

Discussion

AFB1, a highly potent hepatocarcinogen produced by Aspergillus species, remains a critical global health concern due to its strong association with HCC (26). Epidemiological studies have consistently demonstrated a dose-dependent relationship between dietary AFB1 exposure and HCC risk, particularly in regions with high prevalence of hepatitis B virus (HBV) infection, where synergistic effects amplify carcinogenesis (27). The International Agency for Research on Cancer (IARC) classified AFB1 as a Group 1 carcinogen, underscoring its established role in liver cancer development (28).

Recent advances in omics technologies and bioinformatics have significantly enhanced our understanding of AFB1-induced hepatocarcinogenesis (29). Integrated analysis of multi-omics data (genomics, transcriptomics, proteomics) from platforms like TCGA has identified key molecular pathways disrupted by AFB1, including cell cycle regulation (e.g., through CCNB1/CCNB2), p53 signaling, and DNA damage response (30). Network toxicology approaches, combining toxicology with systems biology, have revealed complex interactions between AFB1 and cellular targets, highlighting critical genes such as CHEK1, MMP1, and TTK as potential therapeutic targets (31).

Despite the well-established carcinogenic role of AFB1 in HCC, the development of prognostic prediction models specifically tailored to AFB1-related HCC using network toxicology approaches remains relatively underexplored (32). While traditional epidemiological studies and clinical risk factors have been extensively utilized for HCC prognosis, they often fail to capture the intricate molecular interactions between AFB1 exposure and hepatic carcinogenesis. Network toxicology offers a promising framework to dissect these complex relationships (33). However, its application in predicting the prognosis of AFB1-induced HCC is still in its infancy, with several limitations:

• Limited integration of multi-omics data: most existing studies focus on single-omics analyses rather than integrating multi-omics data to comprehensively identify AFB1-related prognostic biomarkers (34).
• Sparse network-toxicology-based models: few studies have employed network toxicology to systematically identify AFB1-specific target genes and pathways and integrate them into prognostic models (31,35).
• Lack of clinical validation: many proposed models are based on computational predictions or small-scale datasets, with limited validation in independent clinical cohorts (36).
• Overlooked synergistic effects: AFB1’s interaction with other risk factors (e.g., hepatitis B virus infection) is often neglected in model construction, despite its significant impact on HCC progression (27).

Prognostic models based on AFB1-related gene signatures have shown promising predictive accuracy (37). Molecular docking studies further support direct binding of AFB1 to these proteins, particularly CHEK1, with strong binding affinity, suggesting a mechanistic link between AFB1 exposure and HCC progression (38).

Despite these insights, challenges remain in translating findings into clinical practice. Current research focuses on developing early diagnostic biomarkers, targeted therapies (e.g., inhibitors of CHEK1 or p53 pathways), and preventive strategies (39). Future studies aim to integrate single-cell sequencing and spatial transcriptomics to elucidate cell-specific responses to AFB1 and enhance personalized treatment approaches for AFB1-associated HCC (40).

This study aims to systematically identify key targets of AFB1 in HCC through integrated analysis of TCGA data and network toxicology approaches, and subsequently construct a robust prognostic model. The model demonstrated good predictive performance in the TCGA cohort (1-year and 3-year AUC >0.74). Survival analysis showed a statistically significant difference between the two groups (P<0.001). SHAP analysis of the constructed prognostic model demonstrated that CCNB1, CCNB2, and TTK contributed significantly to the model’s predictions. Further inclusion of clinical covariates in the model was performed to evaluate model performance. The model’s Concordance Index (C-index =0.71) demonstrated moderate-to-high predictive discrimination ability (41).

The SHAP analysis revealed that CCNB1 and TTK contribute most significantly to risk prediction through their roles in mitotic regulation and spindle assembly checkpoint control, while CCNB2 exhibits protective effects-findings that are supported by existing literature on these genes’ roles in cell cycle regulation (42). Molecular docking studies provided a mechanistic underpinning for the model by demonstrating strong binding affinities between AFB1 and the proteins encoded by all six prognostic genes. The exceptionally high affinity for CHEK1 (binding free energy =−9.4 kcal/mol) is particularly noteworthy. CHEK1 is a critical kinase in the DNA damage response (43). The potential disruption of its function by direct AFB1 binding could compromise the G2/M checkpoint, permitting the proliferation of cells with AFB1-induced DNA damage and fostering genomic instability—a hallmark of carcinogenesis (44).

The clinical validation of our model through independent Cox analysis (C-index =0.71, P<0.001) underscores its potential for clinical translation into risk stratification protocols (45). The model’s ability to incorporate both genetic markers and clinical parameters addresses a critical need in HCC management, where current staging systems often fail to predict patient outcomes with sufficient accuracy, particularly in early-stage disease (45). These findings collectively support a model where AFB1 induces HCC through multi-hit mechanisms. The identification of these actionable targets may inform development of precision prevention strategies for high-risk populations in endemic regions, including targeted chemoprevention approaches and enhanced surveillance protocols, particularly where vaccination against hepatitis B virus has been implemented with varying success due to vaccine coverage limitations and the persistent threat of AFB1 exposure (46).

Conclusions

Based on integrated bioinformatics and network toxicology analyses, this study established a six-gene prognostic signature (CCNB1, CCNB2, CHEK1, MMP1, TTK, and ESR1) specifically for AFB1-associated HCC. The model demonstrated satisfactory predictive performance (1-year and 3-year AUC >0.74) and effectively stratified patients into high- and low-risk groups with significantly distinct survival outcomes (P<0.001). Molecular docking confirmed strong direct binding between AFB1 and all six encoded proteins, with the highest affinity for CHEK1 (−9.4 kcal/mol). This prognostic signature provides a clinically applicable tool for risk stratification in AFB1-exposed HCC patients and offers mechanistic insights into AFB1-driven hepatocarcinogenesis, warranting further external validation and prospective evaluation.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tgh.amegroups.com/article/view/10.21037/tgh-2026-0053/rc

Peer Review File: Available at https://tgh.amegroups.com/article/view/10.21037/tgh-2026-0053/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tgh.amegroups.com/article/view/10.21037/tgh-2026-0053/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Loomba R, Sanyal AJ, Kowdley KV, et al. Factors Associated With Histologic Response in Adult Patients With Nonalcoholic Steatohepatitis. Gastroenterology 2019;156:88-95.e5. [Crossref] [PubMed]
3. Puvača N, Avantaggiato G, Merkuri J, et al. Occurrence and Determination of Alternaria Mycotoxins Alternariol, Alternariol Monomethyl Ether, and Tentoxin in Wheat Grains by QuEChERS Method. Toxins (Basel) 2022;14:791. [Crossref] [PubMed]
4. Fraquelli M, Paggi S. Risk of hepatocellular carcinoma development in long-term nucleos(t)ide analog-suppressed patients with chronic hepatitis B. Hepatoma Res 2023;9:3.
5. Shah P, Sethuraman A. A novel machine learning approach for tumor detection based on telomeric signatures. Biol Methods Protoc 2025;10:bpaf069. [Crossref] [PubMed]
6. Mei S, Ma H, Chen X. Anticancer and anti-inflammatory properties of mangiferin: A review of its molecular mechanisms. Food Chem Toxicol 2021;149:111997. [Crossref] [PubMed]
7. Tomczak K, Czerwińska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (Pozn) 2015;19:A68-77. [Crossref] [PubMed]
8. Kim S, Chen J, Cheng T, et al. PubChem in 2021: new data content and improved web interfaces. Nucleic Acids Res 2021;49:D1388-95. [Crossref] [PubMed]
9. Daina A, Michielin O, Zoete V. SwissTargetPrediction: updated data and new features for efficient prediction of protein targets of small molecules. Nucleic Acids Res 2019;47:W357-64. [Crossref] [PubMed]
10. Mendez D, Gaulton A, Bento AP, et al. ChEMBL: towards direct deposition of bioassay data. Nucleic Acids Res 2019;47:D930-40. [Crossref] [PubMed]
11. Keiser MJ, Roth BL, Armbruster BN, et al. Relating protein pharmacology by ligand chemistry. Nat Biotechnol 2007;25:197-206. [Crossref] [PubMed]
12. Stelzer G, Rosen N, Plaschkes I, et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinformatics 2016;54:1.30.1-1.30.33.
13. Amberger JS, Bocchini CA, Schiettecatte F, et al. OMIM.org: Online Mendelian Inheritance in Man (OMIM®), an online catalog of human genes and genetic disorders. Nucleic Acids Res 2015;43:D789-98. [Crossref] [PubMed]
14. Davis AP, Wiegers TC, Johnson RJ, et al. Comparative Toxicogenomics Database (CTD): update 2023. Nucleic Acids Res 2023;51:D1257-62. [Crossref] [PubMed]
15. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
16. Szklarczyk D, Kirsch R, Koutrouli M, et al. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res 2023;51:D638-46. [Crossref] [PubMed]
17. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498-504. [Crossref] [PubMed]
18. Chin CH, Chen SH, Wu HH, et al. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol 2014;8:S11. [Crossref] [PubMed]
19. Shen PS. Equivalence tests under the Cox-Aalen model and the partly Aalen model. J Biopharm Stat 2022;32:789-801. [Crossref] [PubMed]
20. Burley SK, Bhikadiya C, Bi C, et al. RCSB Protein Data Bank (RCSB.org): delivery of experimentally-determined PDB structures alongside one million computed structure models of proteins from artificial intelligence/machine learning. Nucleic Acids Res 2023;51:D488-508. [Crossref] [PubMed]
21. Liu Y, Yang X, Gan J, et al. CB-Dock2: improved protein-ligand blind docking by integrating cavity detection, docking and homologous template fitting. Nucleic Acids Res 2022;50:W159-64. [Crossref] [PubMed]
22. Xie Z, Bailey A, Kuleshov MV, et al. Gene Set Knowledge Discovery with Enrichr. Curr Protoc 2021;1:e90. [Crossref] [PubMed]
23. Gene Ontology Consortium. The Gene Ontology knowledgebase in 2023. Genetics 2023;224:iyad031.
24. Kanehisa M, Furumichi M, Sato Y, et al. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res 2023;51:D587-92. [Crossref] [PubMed]
25. Mangalathu S, Hwang SH, Jeon JS. Failure mode and effects analysis of RC members based on machine-learning-based SHapley Additive exPlanations (SHAP) approach. Eng Struct 2020;219:110927.
26. Rushing BR, Selim MI. Aflatoxin B1: A review on metabolism, toxicity, occurrence in food, occupational exposure, and detoxification methods. Food Chem Toxicol 2019;124:81-100. [Crossref] [PubMed]
27. Niu Y, Fan S, Luo Q, et al. Interaction of Hepatitis B Virus X Protein with the Pregnane X Receptor Enhances the Synergistic Effects of Aflatoxin B1 and Hepatitis B Virus on Promoting Hepatocarcinogenesis. J Clin Transl Hepatol 2021;9:466-76. [Crossref] [PubMed]
28. Ostry V, Malir F, Toman J, et al. Mycotoxins as human carcinogens-the IARC Monographs classification. Mycotoxin Res 2017;33:65-73. [Crossref] [PubMed]
29. Wang S, Yang X, Liu F, et al. Comprehensive Metabolomic Analysis Reveals Dynamic Metabolic Reprogramming in Hep3B Cells with Aflatoxin B1 Exposure. Toxins (Basel) 2021;13:384. [Crossref] [PubMed]
30. Gao YN, Yang X, Wang JQ, et al. Multi-Omics Reveal Additive Cytotoxicity Effects of Aflatoxin B1 and Aflatoxin M1 toward Intestinal NCM460 Cells. Toxins (Basel) 2022;14:368. [Crossref] [PubMed]
31. Gao J, Zhang M, Chen Q, et al. Integrating machine learning and molecular docking to decipher the molecular network of aflatoxin B1-induced hepatocellular carcinoma. Int J Surg 2025;111:4539-49. [Crossref] [PubMed]
32. Pan Y, Zhang D, Chen Y, et al. Development and validation of robust metabolism-related gene signature in the prognostic prediction of hepatocellular carcinoma. J Cell Mol Med 2023;27:1006-20. [Crossref] [PubMed]
33. Peng Y, Zhu G, Ma Y, et al. Network Pharmacology-Based Prediction and Pharmacological Validation of Effects of Astragali Radix on Acetaminophen-Induced Liver Injury. Front Med (Lausanne) 2022;9:697644. [Crossref] [PubMed]
34. Chen L, Zhang YH, Wang S, et al. Prediction and analysis of essential genes using the enrichments of gene ontology and KEGG pathways. PLoS One 2017;12:e0184129. [Crossref] [PubMed]
35. Ji H, Song N, Ren J, et al. Systems Toxicology Approaches Reveal the Mechanisms of Hepatotoxicity Induced by Diosbulbin B in Male Mice. Chem Res Toxicol 2020;33:1389-402. [Crossref] [PubMed]
36. Yang L, Li C, Qin Y, et al. A Novel Prognostic Model Based on Ferroptosis-Related Gene Signature for Bladder Cancer. Front Oncol 2021;11:686044. [Crossref] [PubMed]
37. Long J, Wang A, Bai Y, et al. Development and validation of a TP53-associated immune prognostic model for hepatocellular carcinoma. EBioMedicine 2019;42:363-74. [Crossref] [PubMed]
38. Zhu X, Liu S, Pei H, et al. Study on Dihydromyricetin Improving Aflatoxin Induced Liver Injury Based on Network Pharmacology and Molecular Docking. Toxics 2023;11:760. [Crossref] [PubMed]
39. Tawfik RTM, Abd El-Azeem EM, Elsonbaty SM, et al. Green-synthesized selenium-hydroxytyrosol nanocomposites attenuate hepatocellular carcinoma in rats by modulating oxidative stress, inflammation, and apoptosis. Naunyn Schmiedebergs Arch Pharmacol 2025;398:12381-403. [Crossref] [PubMed]
40. Ma L, Hernandez MO, Zhao Y, et al. Tumor Cell Biodiversity Drives Microenvironmental Reprogramming in Liver Cancer. Cancer Cell 2019;36:418-430.e6. [Crossref] [PubMed]
41. Wang W, Wang E, Hong K, et al. m6A regulators-based gene expression pattern is associated with immune microenvironment characteristics in hepatocellular carcinoma. Sci Rep 2025;15:42530. [Crossref] [PubMed]
42. Otto T, Sicinski P. Cell cycle proteins as promising targets in cancer therapy. Nat Rev Cancer 2017;17:93-115. [Crossref] [PubMed]
43. Saldivar JC, Cortez D, Cimprich KA. The essential kinase ATR: ensuring faithful duplication of a challenging genome. Nat Rev Mol Cell Biol 2017;18:622-36. [Crossref] [PubMed]
44. Engin AB, Engin A. DNA damage checkpoint response to aflatoxin B1. Environ Toxicol Pharmacol 2019;65:90-6. [Crossref] [PubMed]
45. Reig M, Forner A, Rimola J, et al. BCLC strategy for prognosis prediction and treatment recommendation: The 2022 update. J Hepatol 2022;76:681-93. [Crossref] [PubMed]
46. Liu Y, Wu F. Global burden of aflatoxin-induced hepatocellular carcinoma: a risk assessment. Environ Health Perspect 2010;118:818-24. [Crossref] [PubMed]