Machine learning‑based identification of key genes underlying sex differences in hepatocellular carcinoma and targeted drug screening

Introduction

Hepatocellular carcinoma (HCC) is one of the leading causes of cancer-related mortalities worldwide. It occupies a dominant position among primary types of liver cancer (1). Due to the insidious onset of HCC in its early stages and the lack of typical symptoms and signs, the disease is delayed, and diagnosis is made at an advanced stage (2). Patients with advanced HCC are frequently unsuitable for surgical treatment due to tumor metastasis, compromised liver function and poor systemic status, leading to a generally poor prognosis (3,4).

Epidemiologic data show that there is a considerable sex difference in the incidence of HCC, with a global prevalence of 2.8-3 times higher in men when compared with women (3,4). In the present study, all analyses regarding these disparities are strictly based on biological sex (men vs. women) rather than sociocultural factors. Studies show that sex hormones carry out an important role in the development of HCC. Estrogens inhibit HCC by regulating the cytokine IL-6, whereas androgens promote HCC by inhibiting lipocalin secretion or directly activating pro-cancer signaling pathways (5,6). Testosterone in particular, leads to reduced levels of lipocalin by inhibiting its secretion, a mechanism that is considered to be one of the most important reasons for the increased incidence of HCC in men (7). However, the effect of these hormones on the sex differences in mechanisms of HCC has not been fully clarified at the molecular level.

Machine learning, as an emerging bioinformatics analysis method, has shown great potential for application in the medical field and has been widely used in the analysis of clinical datasets (8). By constructing robust risk models to aid clinical decision-making and simultaneously redefining patient categories based on data characteristics, precision medicine can be realized in the future (9). Machine learning has been used to analyze high-dimensional transcriptomic data with promising results in identifying key genes for diseases, providing key clues for the study of disease mechanisms (10).

The present study combined differential expression analysis, weighted gene co-expression network analysis (WGCNA) and four machine learning algorithms, including random forest (RF), support vector machine (SVM), generalized linear model (GLM) and extreme gradient boosting (XGBoost), to systematically screen the core genes associated with sex differences in the development of HCC. This was then combined with pan-cancer analysis, functional cell-based assays and virtual drug screening to explore the potential molecular associations of the key genes and identify putative candidate compounds. Specifically, the screening in the present study highlighted Saikosaponin A as a promising therapeutic candidate. This will provide a theoretical basis for understanding the sex disparity in HCC incidence and proposing candidate compounds for further evaluation.

Materials and methods

Cell lines and reagents

Human liver cancer cell lines HepG2 and Huh7, and 293T were purchased from Beijing Beina Chunglian Institute of Biotechnology and The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences. The authenticity of the cell lines used in the present study was verified via short tandem repeat profiling. The plasmids used in the experiments [including the CYP17A1 knockout plasmid (cat. no. L25930), pLenti-Control-sgRNA (cat. no. L00011) and lentiviral packaging vector set A (cat. no. L00002S)] were obtained from Beyotime Biotechnology. According to the manufacturer, the sgRNA incorporated in the commercial CYP17A1 knockout plasmid (cat. no. L25930) was generated using Beyotime's CRISPR/Cas9 sgRNA rapid screening and validation system, and its efficacy had been pre-validated by a T7 endonuclease I (T7EI) assay. The primary reagents included DMEM containing 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.), 0.25% trypsin-EDTA digest (Beijing Solarbio Science & Technology Co., Ltd.), penicillin-streptomycin solution (100X; Beyotime Biotechnology), BCA Protein Quantification Kit (cat. no. 500T; Beyotime Biotechnology) and puromycin (cat. no. ST551; Beyotime Biotechnology). Cell culture conditions were maintained at 37˚C in a humidified incubator with 5% CO2.

Data acquisition and pre-processing

The mRNA expression profiling data of HCC and corresponding paracancerous tissues were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/; accession nos. GSE19665, GSE54236, GSE84402 and GSE121248), which included samples of both sexes and were generated in the original studies by Deng et al (11), Villa et al (12), Wang et al (13) and Wang et al (14). Batch effects were corrected using Surrogate Variable Analysis [sva package, v3.54.0(15)] in R (v4.4.2; R Foundation for Statistical Computing), and the expression matrix of HCC vs. normal tissue for sex groupings was extracted. Differentially expressed genes (DEGs) were screened using limma (v3.62.2) based on |log2 fold change (FC)|≥0.5 and an adjusted P<0.05 after Benjamini-Hochberg correction.

WGCNA

A WGCN was constructed based on the expression matrix of DEGs, and a dynamic tree-cutting algorithm was used to divide the gene modules and screen for those significantly associated with the HCC phenotype (P<0.05 for module significance). The ‘WGCNA’ R package (v1.72) was utilized to construct the network. The network parameters were determined by setting a soft threshold of β=6 (scale-free topological fit index R²>0.85). Modules significantly associated with HCC phenotypes were screened (Pearson correlation; P<0.05) and the top 10% of genes with high intra-module connectivity were extracted as candidate hub genes.

Machine learning-based core gene screening

The WGCNA hub genes were intersected with sex-specific DEGs (yielding a total of 36 candidate genes). Given the relatively limited sample size and the small, highly pre-filtered feature set (only 36 candidate genes), further dividing the dataset into separate training and testing cohorts would severely compromise the training process and the statistical power. Instead, to maximize data utility while strictly preventing overfitting, the present study relied entirely on a rigorous 10-fold cross-validation strategy across the entire dataset. A total of four machine learning algorithms, RF, SVM, GLM and XGBoost, were used to assess gene importance. During the 10-fold cross-validation process, the data were iteratively partitioned into training and validation folds, ensuring that the predictive performance of each model was robustly assessed on unseen data within each fold. The predictive performance was evaluated by calculating the area under the receiver operating characteristic curve (AUC), with all models achieving an AUC >0.95. The overlapping genes (CYP17A1 and IRX3) ranked in the top ten for importance across all four algorithms were ultimately selected as the core genes. Their sex-specific differential expression profiles were further validated using the GEO database.

Bioinformatics validation

For pan-cancer expression analysis, the TIMER2.0 platform (http://timer.cistrome.org/) was used to analyze the expression differences of the core genes across 30 tumors in The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/; Wilcoxon rank-sum test; with P<0.05 as the criterion for determining statistically significant differences). For survival and clinical correlation analyses, TCGA clinical data were integrated and the association between core gene expression and patients' overall survival was evaluated using the Kaplan-Meier method (R package ‘survival’; v3.3.1). Correlations with tumor grade and stage were analyzed using the Spearman rank correlation test. Gene set enrichment analysis was performed to identify the biological pathways associated with the core genes. For immune infiltration analysis, the CIBERSORT algorithm (R package ‘CIBERSORT’; v1.06) was applied to quantitatively analyze the infiltration levels of 22 immune cells and explore the correlation between core gene expression and the immune microenvironment (P<0.05).

Collection of samples

Tumors and corresponding paracancerous tissues of liver cancer patients (n=18; 12 men and 6 women; age range, 40-75 years) were included in the present study. The samples were collected between January 2019 and June 2023. Specifically, samples dating from January 2019 to November 2022 were retrospectively accessed from historical archival formalin-fixed paraffin-embedded blocks, while samples from November 2022 to June 2023 were prospectively collected. Based on the study design, the precise inclusion criteria were: i) Patients who underwent radical resection for liver cancer; ii) pathologically confirmed primary liver cancer; iii) availability of matched paraffin-embedded tumor and paracancerous tissues (with paracancerous tissue strictly defined as being >10 mm away from the tumor margin). The exclusion criteria were: i) Patients who received preoperative radiotherapy, chemotherapy or targeted therapy; ii) presence of other primary malignancies. Following surgical resection, fresh tissues were snap-frozen and stored at -80˚C within 30 min, while matched tissue blocks were routinely formalin-fixed and paraffin-embedded for subsequent analysis. The present study was conducted with the approval of the Ethics Committee of the First Affiliated Hospital of Dali University (Approval no. OFY20221122001), which explicitly covered both the retrospective use of archived material and the prospective collection of new samples. Written informed consent was obtained from all patients prior to their participation; notably, for the newer prospectively collected samples, dedicated written informed consent for research use was strictly obtained prior to surgery. Strict confidentiality of patient identities was maintained. Detailed clinicopathological characteristics, including age, biological sex, tumor stage and viral hepatitis status, were retrospectively collected from medical records.

Immunohistochemistry (IHC)

Paraffin sections were deparaffinized with xylene and rehydrated with descending gradient ethanol. Antigen retrieval was performed in sodium citrate buffer (pH 6.0) at 100˚C for 20 min. After blocking with goat serum blocking solution (cat. no. C0265; Beyotime Biotechnology) at room temperature for 10 min, the sections were incubated in a humidified chamber with rabbit anti-CYP17A1 primary antibody (1:200; cat. no. bsm-54306R; Bioss) overnight at 4˚C, followed by an HRP-labeled goat anti-rabbit secondary antibody (1:1,000; cat. no. G1213-100UL; Wuhan Servicebio Technology) at 37˚C for 30 min. Visualization was achieved using DAB, followed by hematoxylin counterstaining at room temperature for 15 sec. Staining intensity (0-3) and the percentage score of positive cells (1-4) were independently evaluated by two pathologists in a double-blind manner; the final IHC score was calculated as intensity x percentage score, and discordant cases were reviewed jointly to reach consensus.

Knockout cell line construction and functional experiments

To construct stable CYP17A1 knockout cells, 293T cells were transfected at 37˚C using a third-generation lentiviral packaging system consisting of the CYP17A1 knockout plasmid (cat. no. L25930; Beyotime Biotechnology), pMDLg, Rev and VSV-g plasmids. For each 10-cm dish, 10 µg transfer plasmid, 6.5 µg pMDLg, 2.5 µg Rev and 3.5 µg VSV-g were mixed in 500 µl serum- and antibiotic-free DMEM and transfected using 40 µl Lipo6000 (cat. no. C0526; Beyotime Biotechnology) at 37˚C for 6 h. The medium was replaced after 4-6 h, and viral supernatants were collected at 24 and 48 h, filtered through a 0.45-µm membrane, and used to infect HepG2 cells. Because a fixed MOI was not determined, transduction conditions were optimized empirically using 0.5-2.0 ml viral supernatant/well according to transduction efficiency; 1 ml viral supernatant plus Polybrene (8 µg/ml; cat. no. ST1380; Beyotime Biotechnology) was used for routine infection in 6-well plates. Cells were transduced for 48 h, after which puromycin selection was initiated at 2 µg/ml and maintained at the same concentration for >1 week until stable knockout cells were obtained. The CYP17A1 sgRNA sequence was 5'-GCACCAGGGCACCTTCTCTT-3', and an empty vector was used as the negative control (NC).

Reverse transcription-quantitative PCR (RT-qPCR)

Total RNA was extracted from cells using TRIzol® reagent (cat. no. 15596026; Invitrogen; Thermo Fisher Scientific, Inc.). Reverse transcription was performed using HiScript III All-in-one RT SuperMix Perfect for qPCR (cat. no. R333-01; Vazyme Biotech Co, Ltd.) according to the manufacturer's instructions (37˚C for 15 min and 85˚C for 5 sec). qPCR was performed using ChamQ Universal SYBR qPCR Master Mix (cat. no. Q311-02; Vazyme Biotech Co, Ltd.) on cDNA derived from the indicated cells, and relative mRNA expression was calculated using the 2-ΔΔCq method (16) with GAPDH as the internal reference. The thermocycling conditions were 95˚C for 30 sec, followed by 40 cycles of 95˚C for 10 sec and 60˚C for 30 sec. The forward and reverse sequences of all primers used were as follows: CYP17A1 Forward: 5'-AGAAGTTATCATCAATCTGTGGGC-3', Reverse: 5'-CTGCTCCGAAGGGCAAATA-3'; GAPDH Forward: 5'-GGAAGCTTGTCATCAATGGAAATC-3', Reverse: 5'-TGATGACCCTTTTGGCTCCC-3'.

Western blotting

Total protein was extracted from cells using RIPA lysis buffer (cat. no. R0010; Beijing Solarbio Science & Technology) supplemented with PMSF. Protein concentrations were determined using a BCA Protein Assay kit. Equal amounts of protein (20 µg/lane) were separated by 10% SDS-PAGE and transferred onto PVDF membranes (cat. no. IPVH00010; MilliporeSigma) at 200 mA for 90 min. The membranes were blocked with 5% non-fat milk in TBST (0.1% Tween-20) at room temperature for 2 h and incubated with rabbit anti-CYP17A1 antibody (1:1,000; cat. no. bsm-54306R; Bioss) and mouse anti-GAPDH antibody (1:1,000; cat. no. GB12002-100; Wuhan Servicebio Technology) overnight at 4˚C, followed by HRP-conjugated goat anti-rabbit IgG (1:5,000; cat. no. bs-0295G-HRP; Bioss) or HRP-conjugated goat anti-mouse IgG (1:5,000; cat. no. bs-0296G-HRP; Bioss) for 1 h at room temperature. Signals were visualized using BeyoECL Plus chemiluminescence reagent (cat. no. P0018S; Beyotime Biotechnology).

Cell proliferation

Cell Counting Kit-8 (CCK-8) was used to measure cell proliferation (Dojindo Laboratories, Inc.). Cells were seeded in 96-well plates (5x10³ cells/well; six replicates per group) and cell viability was assessed at 0, 24, 48 and 72 h by measuring the OD450. In the functional cell-based assays, the KO-CYP17A1-HepG2 cell line was used as the experimental group, and the KO-Control-HepG2 cell line was used as the control group, to explore the effects of CYP17A1 knockout on the proliferation, migration and invasion ability of the cells.

Cell migration

Cell migration was evaluated using a scratch wound-healing assay. For CYP17A1 knockout experiments, HepG2 NC and KO cells were used; for Saikosaponin A experiments, HepG2 and Huh7 cells were treated with 20 µM Saikosaponin A (Aladdin Biochemical Technology Co., Ltd.). Cells were seeded in 6-well plates and allowed to reach ~100% confluence, scratched with a 200-µl pipette tip, washed three times with PBS, and then cultured in complete medium. Images were captured at 0 and 24 h under an inverted light microscope (magnification, x100), and wound-healing area was quantified using ImageJ (v1.53; National Institutes of Health). For invasion assays, Matrigel (cat. no. M8370; Beijing Solarbio Science & Technology Co., Ltd.) was diluted to ~200 µg/ml, 100 µl was added to each upper chamber and allowed to polymerize at 37˚C for 4 h. Transwell inserts (8-µm pore size; cat. no. WG3422; Wuhan Servicebio Technology) were used. Cells were resuspended at 5x105 cells/ml in serum-free DMEM, and 200 µl cell suspension was added to the upper chamber, while 600 µl complete medium containing 10% FBS was added to the lower chamber. After incubation for 24 h at 37˚C, cells on the lower membrane surface were fixed with 4% paraformaldehyde for 20 min at room temperature, stained with 0.1% crystal violet for 10 min at room temperature, and counted in five random fields under an inverted light microscope (magnification, x200).

Virtual drug screening and molecular docking

The crystal structure of the CYP17A1 protein was retrieved from the PDB database (https://www.rcsb.org/). Water molecules and co-crystallized ligands were removed using PyMOL (v2.5.2; https://pymol.org/2/). A Traditional Chinese Medicine monomer library containing 4,080 compounds was obtained from TargetMol (https://www.targetmol.com/compound-library/traditional-chinese-medicine-monomers-library.html) and virtually screened using AutoDock Vina (v1.2.3). For compounds demonstrating binding energies ≤-10 kcal/mol, molecular docking interactions with the CYP17A1 protein were further analyzed using PyMOL and PLIP software (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index). Saikosaponin A (Aladdin Biochemical Technology Co., Ltd.) exhibited the highest binding affinity for CYP17A1 and was therefore selected for subsequent in vitro validation. The concentration gradient of Saikosaponin A was 0-30 µM. In each experiment, the changes in cell proliferation, migration and invasion ability of HepG2 and Huh7 cells were observed under different concentrations of Saikosaponin A, and the dose-dependent inhibitory effects were analyzed.

Statistical analysis

Experimental data were analyzed using R (v4.4.2; R Foundation for Statistical Computing), GraphPad Prism (v6.0; Dotmatics) and ImageJ (v1.53; National Institutes of Health). Data are presented as mean ± SD. Two-group comparisons were performed using unpaired, two-tailed Student's t-tests. For >2 groups, one-way ANOVA followed by Tukey's post hoc test was used. Time-course proliferation data (CCK-8) were analyzed using two-way ANOVA with Bonferroni post hoc correction. Differential expression analyses were performed using limma (v3.62.2) with Benjamini-Hochberg false discovery rate correction. Pearson correlation analysis was used for WGCNA module-trait associations, Wilcoxon rank-sum tests were used for pan-cancer expression comparisons, and Spearman rank tests were used for correlations with tumor grade and stage. Survival analyses were performed using Kaplan-Meier curves and two-sided log-rank tests. Because late-stage crossover/convergence was observed in Fig. 3C, sensitivity analyses were additionally conducted using i) administrative truncation at 7 years, ii) weighted Gehan-Breslow testing (rho=1) and iii) a two-stage test [TSHRC package (v0.1-6)]. P<0.05 was considered statistically significant.

Figure 3 Pan-cancer analysis of CYP17A1 and single-gene analysis in HCC. (A) Expression of CYP17A1 in various tumors. (B) Expression of CYP17A1 in HCC and paracancerous tissues. (C) Kaplan-Meier survival curve of CYP17A1 high- and low-expression groups. (D) Correlation of CYP17A1 expression with M1 macrophages. (E) Correlation of CYP17A1 expression with activated dendritic cells. HCC, hepatocellular carcinoma; **P<0.01; ***P<0.001.

Results

Identification of core genes for sex differences in HCC

A total of 36 candidate genes were screened from samples from male and female patients with HCC by differential expression analysis combined with WGCNA (Fig. 1A and B). The evaluation results showed that the AUC of all models was >0.95 (Fig. 1C) and the residual analysis showed that the model residual values were <0.25 (Fig. 1D), which suggested that the constructed models achieved high discriminatory ability within the cross-validation framework. An intersectional analysis of the top ten genes ranked in importance by the four algorithms identified CYP17A1 and IRX3 as the core genes differentially expressed in HCC samples from men (Fig. 1E).

Figure 1 Screening for genes central to sex differences in hepatocellular carcinoma incidence by four machine learning algorithms. (A) Genes that are differentially expressed only in male or female patients intersect with the hub genes of WGCNA. (B) A total of four machine learning algorithms score the importance of 36 candidate core genes. (C) Receiver operating characteristic curve of four machine learning algorithms on model prediction accuracy. (D) Residual values of predictive accuracy of four machine learning algorithms. (E) Venn diagram of the top 10 most important genes for four machine learning algorithms. WGCNA, weighted gene co-expression network analysis; Diffgene, differentially expressed genes; RF, random forest; SVM, support vector machine; XGB, extreme gradient boosting; GLM, generalized linear model.

Validation of sex-specific expression of CYP17A1 and IRX3

To verify the differences in the expression of CYP17A1 and IRX3 in HCC tissues of different sexes, their mRNA expression levels were further evaluated using the comprehensive GEO datasets. The results showed that CYP17A1 and IRX3 mRNA expression levels in HCC tissues from men were significantly higher when compared with those in paired paracancerous tissues (P<0.001, Fig. 2A and C). However, there was no significant difference in CYP17A1 and IRX3 mRNA expression between HCC tissues from women and their paired paracancerous tissues (P>0.05; Fig. 2B and D). Subgroup analysis further strictly compared tumor tissues to their sex-matched normal controls. The results showed that the expression levels of CYP17A1 and IRX3 were elevated in the male patient HCC group compared with the normal male group (both P<0.001, Fig. 2E and F). These data verified that CYP17A1 and IRX3 were differentially expressed only in male patient HCC, suggesting a sex-specific regulatory pattern. Notably, while CYP17A1 was subjected to further experimental validation in the present study, IRX3 lacks functional validation and currently serves as a hypothesis-generating candidate gene warranting future investigation.

Figure 2 Validation of sex-specific expression of CYP17A1 and IRX3. (A) Expression of CYP17A1 in male patients with HCC. (B) Expression of CYP17A1 in female patients with HCC. (C) Expression of IRX3 in male patients with HCC. (D) Expression of IRX3 in female patients with HCC. (E) Differences in CYP17A1 expression among four groups. (F) Differences in IRX3 expression among four groups. Con, control; Treat, treated; HCC, hepatocellular carcinoma; ns, not significant; *P<0.05; ***P<0.001.

Pan-cancer exploratory analysis and HCC-specific clinical relevance of CYP17A1

To explore the broader context of CYP17A1, a pan-cancer expression analysis was performed. It showed that CYP17A1 was generally lowly expressed in the majority of tumors, such as cholangiocarcinoma and renal cancer, but exhibited specific high expression in HCC, rectal adenocarcinoma and gastric cancer (P<0.001; Fig. 3A). In a subsequent, strictly independent HCC-specific validation using the TCGA-LIHC dataset, the expression of CYP17A1 in HCC tissues was significantly increased compared with that in normal tissues (P<0.001; Fig. 3B). The survival analysis with the standard log-rank test indicated no statistically significant association with overall survival (log-rank P>0.05; Fig. 3C). Considering the late-stage curve crossover, additional sensitivity analyses were performed: A 7-year truncated log-rank (P>0.05), weighted Gehan-Breslow (P<0.05) and two-stage test (P=>0.05), indicating a trend-level but non-robust prognostic association. Furthermore, the correlation analysis of CYP17A1 expression and immune cell infiltration showed a positive correlation with M1 macrophage infiltration (R=0.33; P<0.05; Fig. 3D) and a negative correlation with activated dendritic cell infiltration (R=-0.3; P<0.05; Fig. 3E).

Validation of protein expression and functional mechanisms for CYP17A1

Validation of CYP17A1 protein expression by IHC suggested that CYP17A1 protein was significantly higher in male patient HCC tissues when compared with that in paracancerous tissues (P<0.001; Fig. 4A and B). By contrast, there was no significant difference in female patient samples (P>0.05; Fig. 4C and D). This protein-level result was consistent with the mRNA expression patterns. To investigate the functional mechanisms of CYP17A1 in vitro, stable knockout of CYP17A1 was established using lentivirus specifically in the HepG2 cell line. The mRNA (Fig. 5A) and protein (Fig. 5B and C) expression levels of HepG2 cells were significantly reduced following CYP17A1 knockout (P<0.001). Functional cell-based assays demonstrated that knockout of CYP17A1 inhibited HepG2 cell proliferation at 72 h (P<0.001; Fig. 5D; Table SI), reduced migration capacity by 78% in the scratch healing assay (P<0.001; Fig. 5E and F) and attenuated invasive capacity by 65% in the Transwell assay (P<0.001; Fig. 5G and H). These results suggest that CYP17A1 potentially contributes to the malignant phenotypes of HCC cells in vitro.

Figure 4 Immunohistochemical validation of CYP17A1 protein in HCC tissues. (A) Immunohistochemical images from HCC and paired paracancerous tissues from a male patient. (B) Difference in CYP17A1 expression between HCC and paired paracancerous tissues from a male patient. (C) Immunohistochemical images from HCC and paired paracancerous tissues from a female patient. (D) Difference in CYP17A1 expression between HCC and paired paracancerous tissues from a female patient. HCC, hepatocellular carcinoma; ns, not significant; ***P<0.001.

Figure 5 A preliminary investigation of the role of CYP17A1 in hepatocellular carcinoma based on CYP17A1 knockout and functional cell-based assays. (A) Relative expression of CYP17A1 mRNA in the control group and the KO group. (B) CYP17A1 protein expression in the control group and the KO group. (C) Difference analysis of CYP17A1 relative expression levels in the control group and the KO group. (D) Proliferation of the control group and the CYP17A1 gene KO group over 72 h. (E) Effect of CYP17A1 gene knockout on HepG2 migration ability under the microscope at 0 and 24 h (magnification, x100). (F) Difference analysis of scratch healing rates between the control group and the KO group. (G) Transwell results of the control group and the KO group under a microscope. (H) Analysis of the difference in the relative number of invasive cells in the control group and the KO group. NC, the control group; KO, the knockout group; OD, optical density; ***P<0.001.

Virtual screening and in vitro validation of putative CYP17A1-targeting compounds

In a virtual screening of 4,080 herbal monomers using the CYP17A1 protein as a receptor, molecular docking predicted that Saikosaponin A possessed the highest binding affinity for the CYP17A1 protein, with a binding energy of -12.5 kcal/mol and the formation of 6 pairs of hydrogen bonds (Table I; Fig. 6A). Based on these predictive docking data, Saikosaponin A was selected as a putative compound targeting CYP17A1 for further cellular validation. In vitro experiments evaluated the effects of Saikosaponin A on HCC cell phenotypes. Cell proliferation assays demonstrated that treatment with different concentrations of Saikosaponin A significantly inhibited the proliferation of HepG2 and Huh7 cells in a dose- and time-dependent manner (Fig. 6B and C; Tables SII and SIII). Furthermore, after treatment with 20 µM of Saikosaponin A for 24 h, the scratch healing rates of HepG2 and Huh7 cells were reduced by 63 and 58%, respectively (P<0.001; Fig. 6D-G). In response to the image-quality query, representative wound-healing control images were replaced with fields showing continuous near-confluent monolayers outside the wound area, and quantification was recalculated from qualified fields (Fig. 6F). Similarly, Transwell assays revealed that the number of invading HepG2 and Huh7 cells was reduced by 52 and 48%, respectively (P<0.01; Fig. 6H-J). These inhibitory effects on migration and invasion were also dose-dependent. This suggests that Saikosaponin A exhibits a significant inhibitory effect on the malignant phenotypes of HCC cells, highlighting it as a potential therapeutic candidate, although direct CYP17A1-dependency requires further rescue validations.

Figure 6 Screening and in vitro validation of drug candidates targeting CYP17A1. (A) Saikosaponin A interacts with the CYP17A1 protein. (B) Cell proliferation of HepG2 cells treated with different concentrations of Saikosaponin A for 24, 48 and 72 h. (C) Cell proliferation of Huh7 cells treated with different concentrations of Saikosaponin A for 24, 48 and 72 h. (D) Scratch healing of HepG2 control and Saikosaponin A-treated cells at 0 and 24 h under the microscope (magnification, x100). (E) Difference analysis of the scratch healing rate between the Control group and the drug-treated group. (F) Scratch healing of Huh7 Control and Saikosaponin A-treated cells at 0 and 24 h under the microscope (magnification, x100). (G) Difference analysis of scratch healing rates between the Control group and the drug-treated group. (H) Cell invasion of HepG2 and Huh7 cells in the Control group and the drug-treated group under the microscope. (I) Difference analysis of the relative cell number of HepG2 cells in the Control group and the drug-treated group. (J) Difference analysis of relative cell numbers of Huh7 cells in the Control group and the drug-treated group. *P<0.05; **P<0.01; ***P<0.001.

Table I Virtual screening of interactions between the top 10 drugs with the highest affinity and CYP17A1 protein.

Table I

Virtual screening of interactions between the top 10 drugs with the highest affinity and CYP17A1 protein.

Ligand	Chemical formula	Pubchem no.	Molecular weight	Affinity (kcal/mol)	Interactions
Wilforlide A	C30H46O3	158477	454.7	-13.7	5 hydrophobic interactions
Epifriedelanol	C30H52O	119242	428.7	-13.4	7 hydrophobic interactions
δ-Amyrenone	C30H48O	5318261	424.7	-13.2	6 hydrophobic interactions
Arnidiol	C30H50O2	10478550	442.7	-13.1	7 hydrophobic interactions
β-Amyrin	C30H50O	73145	426.7	-12.8	7 hydrophobic interactions
Taraxasterol	C30H50O	115250	426.7	-12.7	7 hydrophobic interactions
Saikosaponin A	C42H68O13	167928	781	-12.5	4 hydrophobic interactions, 6 H-bonds
Taraxastery acetate	C32H52O2	13889352	468.8	-12.4	7 hydrophobic interactions
Glycyrrhetinic acid	C30H46O4	10114	470.7	-12.3	6 hydrophobic interactions
Saikosaponin D	C42H68O13	107793	781	-12.3	4 hydrophobic interactions, 5 H-bonds

Discussion

HCC is considerably more prevalent in men than in women (3,4), and exploring the underlying molecular mechanisms of this sex difference is a prominent research topic in oncology. In the present study, by integrating multi-omics data and machine learning algorithms, CYP17A1 and IRX3 were computationally identified as key candidate genes associated with sex differences in HCC. The present study explored their potential mechanisms in influencing the progression of HCC, possibly through the modulation of the androgen signaling pathway. Virtual drug screening and in vitro functional validation revealed that the traditional Chinese medicine monomer Saikosaponin A could inhibit the malignant phenotypes of HCC cells, highlighting a putative new candidate for the treatment of HCC.

Cytochrome P450 enzyme (CYP17A1) is a key steroidogenic enzyme that carries out an important role in the development of several tumors. It catalyzes the conversion of pregnenolone to dehydroepiandrosterone and promotes testosterone synthesis (17,18). Niu et al (19) found that CYP17A1 knockout in vitro considerably inhibits the proliferation and invasive ability of glioma cells, while promoting apoptosis. Moreover, in women, it has also been hypothesized that the androgen receptor (AR) can substitute for estrogen-dependent signaling to stimulate transcription of steroid-responsive genes that drive breast cancer (20). Altered expression levels of genes involved in the androgen expression pathway have been suggested as one of the key factors that may be associated with sex differences in the development of these tumors (21-23). The present study found that CYP17A1 was specifically highly expressed in male patient HCC tissues. This result aligns with the mechanism that CYP17A1 drives tumor progression by promoting androgen synthesis in prostate cancer. At the same time, previous studies have found that cell culture medium containing certain levels of androgens promotes the proliferative capacity of glioma cells and that average serum testosterone levels are markedly higher in patients with glioma (24,25). Taken together, the present study speculates that the higher incidence of HCC in men may be associated with CYP17A1-mediated elevations of local androgen levels, which potentially promote tumor cell proliferation through the activation of AR signaling pathways. Notably, while IRX3 was identified alongside CYP17A1 as a core gene, to the best of our knowledge, its association with sex differences has not yet been previously reported. IRX3 regulates cell differentiation in embryonic development, and its function in HCC may involve aberrant activation of the Wnt/β-catenin pathway. However, because no functional validation was performed for IRX3 in the present study, it must be viewed as a hypothesis-generating candidate gene that requires rigorous future experimental validation.

In the field of tumor therapy, natural products have been an important source for the development of chemotherapeutic drugs (26). In order to explore potential interventions, the present study found that Saikosaponin A, an important bioactive triterpene glycoside compound extracted from Chaihu (Radix Bupleuri), has a predicted binding affinity with CYP17A1(27). Du et al (28) found that Saikosaponin A triggered apoptosis and inhibited the PI3K/Akt signaling pathway, exerting an anti-cervical cancer effect. Wang et al (29) showed that Saikosaponin A may induce apoptosis in gastric cancer cells. Zhang et al (30) found that Saikosaponin A has potent anti-angiogenic activity and inhibits tumor growth mainly by blocking the VEGFR2 signaling pathway. In the present study, molecular docking predicted a high affinity (binding energy=-12.5 kcal/mol) with the CYP17A1 protein, and in vitro experiments confirmed that it dose-dependently inhibited HCC cell migration and invasion. However, molecular docking alone does not establish direct enzymatic inhibition, presenting a limitation of the present study. Without target-dependency or rescue experiments, it remains inconclusive whether CYP17A1 is the primary functional target mediating the observed anti-migration and anti-invasion effects, as Saikosaponin A may act through pleiotropic mechanisms. In the future, enzymatic activity assays and rescue experiments in CYP17A1-knockout cells are necessary to comprehensively reveal its precise antitumor mechanism.

The sex-specific expression of CYP17A1 in HCC provides a new direction for the development of personalized treatment regimens for male patients. Abiraterone, an inhibitor of CYP17A1, has achieved significant efficacy in the treatment of prostate cancer, suggesting that targeting this gene axis might also be applicable to HCC (31). Despite these promising results, several limitations must be addressed. First, the clinical IHC sample size was small (n=18), and the subgroup analyses did not adequately control for potential confounding factors such as age, viral hepatitis status and alcohol consumption. Thus, it is premature to attribute the observed expression differences of CYP17A1 solely to biological sex without validation in a larger, multifactorial clinical cohort. Second, as indicated by the primary log-rank analysis (P=0.0726) and additional crossover-sensitive analyses (7-year truncated log-rank P=0.0933; weighted Gehan-Breslow P=0.0334; two-stage P=0.0908), the association between CYP17A1 expression and overall survival was not robust across methods, suggesting its role may be mechanistic or therapeutic rather than acting as a robust prognostic biomarker. Third, owing to the scarcity of comprehensive multi-omics HCC datasets with perfectly matched biological sex annotations, we were unable to construct a fully independent external validation cohort, which implies that the risk of overfitting cannot be entirely ruled out. Finally, the involvement of the androgen pathway is inferred from previous literature and our bioinformatics findings, rather than directly demonstrated through androgen level measurements or rescue experiments in this study. The in vivo efficacy and safety of Saikosaponin A also require systematic evaluation using humanized HCC mouse or patient-derived xenograft models in future studies.

In conclusion, the present study suggests that CYP17A1 is a key gene associated with HCC in men, potentially driving malignant phenotypes via the androgen signaling axis, while IRX3 emerges as a novel hypothesis-generating candidate. Furthermore, Saikosaponin A is highlighted as a putative therapeutic candidate. Future research efforts will focus on constructing a CYP17A1 conditional knockout mouse model to validate its necessity in HCC occurrence, performing target-dependency assays for Saikosaponin A and resolving the interaction network of IRX3 with sex-related signaling pathways. These studies are expected to provide a more solid theoretical basis for sex-specific strategies and personalized treatment in HCC.

Supplementary Material

Cell Counting Kit-8 assay of control and CYP17A1 knockout groups.

Inhibitory effects of Saikosaponin A on the proliferation of HepG2 cells assessed by Cell Counting Kit-8 assay.

Inhibitory effects of Saikosaponin A on the proliferation of Huh7 cells assessed by Cell Counting Kit-8 assay.

Acknowledgements

Not applicable.

Funding

Funding: The author(s) declare that financial support was received for the research, authorship and/or publication of this article. The present study was supported by the following projects: General Project of the Joint Special Project of Local Universities in Yunnan Province (grant nos. 202001BA070001-064 and 202101BA070001-102), Dali University Doctoral Research Start-up Fund Project (grant no. KYBS2018012), Open Project of Yunnan Provincial Key Laboratory of Entomological Biopharmaceutical R&D (grant no. AG2024002), Clinical Medicine Discipline Team Building Project of the First Affiliated Hospital of Dali University (grant no. DFYYB2024026), Open Project of Key Laboratory of Screening and Research of Resistant Plant Resources in West Yunnan, Yunnan Province (grant no. APKL2101) and Yunnan College Students Innovation and Entrepreneurship Training Project Fund (grant no. S202410679068).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author. The datasets analyzed for this study can be found in the GEO database (https://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE19665, GSE54236, GSE84402 and GSE121248, and in the TCGA-LIHC dataset via the TIMER2.0 platform.

Authors' contributions

ZW and JN contributed to conceptualization, data curation, formal analysis, investigation, methodology, project administration, supervision, CYP17A1 protein expression and functional validation, visualization, writing of the original draft, and review and editing. HZ and WL contributed to conceptualization, formal analysis, investigation, supervision, visualization, and review and editing. ZL and XY contributed to methodology, bioinformatics data analysis, visualization and writing of the original draft. YZhao and YW contributed to data curation, formal analysis and writing of the original draft. RQ contributed to data processing, conducting experiments, visualization and writing of the original draft. CL contributed to bioinformatics data analysis, visualization and writing of the original draft. QR contributed to the acquisition of clinical data, interpretation of data, and review of the manuscript. YL and YZhang contributed to conceptualization, data curation, funding acquisition, investigation, methodology, project administration, resources, supervision, visualization, writing of the original draft and review and editing.

Ethics approval and consent to participate

The study of human tumor cells was approved by the ethics committee of the First Affiliated Hospital of Dali University, Dali, China (approval no. OFY20221122001). The present study was conducted in accordance with the local legislation and institutional requirements.

Patient consent for publication

All patients provided written informed consent for the publication of their anonymized clinical data obtained and analyzed in the present study. The consent process ensured that participants were aware that their data would be used for research purposes and potentially included in scientific publications, with all personally identifiable information removed to protect their privacy.

Competing interests

The authors declare that they have no competing financial interests.

References