Data collection and preprocessing
Gene expression data and clinical information for HCC were primarily sourced from TCGA data portal (https://portal.gdc.cancer.gov/), which includes 374 HCC samples and their corresponding clinical data. Additionally, we retrieved the GSE54236 HCC dataset from the GEO database, which contains gene expression data and clinical information for 80 HCC samples. To ensure consistency and comparability across datasets, batch effects between datasets were corrected using the ComBat function from the sva R package. Furthermore, we extracted eQTL GWAS data from the IEU OpenGWAS project (https://gwas.mrcieu.ac.uk/), covering 19,942 genes, and selected single nucleotide polymorphisms (SNPs) with a significance threshold of P < 5e−6 for further analysis. In addition, we integrated HCC-related GWAS data from the FinnGen R12 dataset (https://www.finngen.fi/fi) and the GWAS Catalog, specifically the GCST90013702 dataset. The FinnGen dataset includes data from 947 HCC patients and 378,749 healthy controls, while the GCST90013702 dataset includes data from 1,866 HCC patients and 195,745 healthy controls, providing valuable genetic association information for subsequent in-depth analysis.
MR analysis
In this study, a two-sample MR framework was implemented using the “TwoSampleMR” R package (v0.5.6) to infer the causal relationship between gene expression (eQTLs) and HCC risk. We selected independent instrumental variables (IVs) by applying linkage disequilibrium pruning (r2 = 0.1, kb = 10,000), and palindromic SNPs were excluded to avoid strand ambiguity. The inverse variance weighted (IVW) method was employed as the primary analysis strategy due to its superior statistical power under the assumption of valid instruments. To account for potential violations of the MR assumptions—such as horizontal pleiotropy—we conducted sensitivity analyses using multiple complementary MR methods, including MR-Egger regression, weighted median, and mode-based estimation. These methods provided robustness checks under different model assumptions and helped assess the consistency of causal estimates. The MR-Egger intercept test was used to detect directional pleiotropy, while heterogeneity was assessed via Cochran’s Q statistic. Additionally, we performed leave-one-out analysis to evaluate the influence of individual SNPs on the causal estimate. These sensitivity analyses confirmed the robustness of the MR results and minimized the likelihood of bias due to invalid instruments.
To address heterogeneity in the HCC GWAS data, we first integrated the FinnGen R12 (947 HCC cases/378,749 controls) and GCST90013702 (1866 HCC cases/195,745 controls) datasets using Metal software12, corrected for batch effects, and generated summary statistics, which were then used for the final MR analysis.
Feature selection and riskScore model construction
In this study, a prognostic prediction model for HCC was developed by integrating multiple analytical approaches. Univariate Cox regression analysis (coxph function) was performed using the “survival” R package (v3.5.7) and the “survminer” R package (v0.4.9) to identify candidate genes significantly associated with overall survival in HCC patients (selection criterion: P < 0.05). These candidate genes were cross-validated with eQTL genes identified through Mendelian randomization analysis. A Venn diagram was used to determine a feature gene set with both genetic association and clinical prognostic significance. Based on the expression profiles of these feature genes, the classification performance of 101 machine learning algorithms (including Lasso + StepCox, survival-SVM, etc.) was systematically evaluated. The 101 methods used by the machine learning model and their corresponding hyperparameters are presented in Table S9. The optimal algorithm for constructing the riskScore model was selected through cross-validation and C-index. Patients were stratified into high-risk and low-risk groups using the median riskScore as the threshold. Kaplan–Meier curves and Log-rank tests were employed to assess the prognostic discriminative power of this stratification.
To mitigate overfitting and ensure model generalizability, we implemented a rigorous tenfold cross-validation strategy for all machine learning algorithms during model training. Specifically, in the RunEnet function (used for Lasso, and Ridge), cross-validation was performed using the cv.glmnet function with nfolds = 10, where the dataset was randomly partitioned into 10 equally sized subsets. In each iteration, one subset was held out for validation while the remaining nine were used for training. This process was repeated across all folds, and the optimal penalty parameter (lambda) was selected based on the minimum cross-validated error. For other algorithms such as CoxBoost, GBM, and SuperPC, internal cross-validation procedures (cv.CoxBoost, cv.gbm, superpc.cv) with K = 10 folds were applied to select hyperparameters and assess performance stability. These consistent tenfold cross-validation frameworks were used to reduce the risk of overfitting, improve model robustness, and select the best-performing model across 101 algorithmic configurations.
Survival analysis and model validation
This study employed a multi-dimensional validation framework to assess the prognostic performance of the risk model. First, time-dependent ROC curves for the risk model and clinical indicators (such as age and stage) were generated using the “timeROC” R package (v0.4.3) to evaluate 1-year, 3-year, and 5-year survival rates. The “rms” R package (v6.7-0) and “pec” R package were then used to calculate C-index, quantifying the model’s predictive accuracy. The parameters for the C-index calculation were set as splitMethod = “bootcv” and B = 1000. Subsequently, univariate and multivariate Cox regression models (using the “survival” R package) were applied to evaluate the statistical significance of riskScore as an independent prognostic factor (hazard ratios and 95% confidence intervals). Finally, Kaplan–Meier survival curves were constructed based on the median riskScore to stratify patients into high-risk and low-risk groups. The Log-rank test was employed to assess the significance of survival differences between the two groups (P < 0.05).
Functional enrichment analysis
GSEA Pathway Enrichment Analysis: Differential gene expression analysis was performed using the “limma” R package (v3.58.1). Pathway enrichment analysis was conducted with the “clusterProfiler” R package (v4.12.0) on gene sets from KEGG (c2.cp.kegg_legacy.v2024.1.Hs.symbols.gmt), Gene Ontology (c5.go.v2024.1.Hs.symbols.gmt), and Hallmark (h.all.v2024.1.Hs.symbols.gmt), selecting significant pathways with an FDR < 0.05.
Immune-Related Pathway Analysis: The single-sample Gene Set Enrichment Analysis (ssGSEA) method was used to quantify the activity of immune-related pathways. Additionally, algorithms such as CIBERSORT were employed to calculate the abundance of 28 immune cell subtypes. Spearman correlation analysis was performed to evaluate the association between the riskScore and immune infiltration (P < 0.05).
Immune cell infiltration assessment
This study utilized multiple methods to estimate immune cell infiltration, including Quantiseq, Timer, Mcp_counter, Xcell, CIBERSORT, and EPIC algorithms, for a comprehensive assessment of immune cell distribution patterns. Tumor purity and the composition of stromal immune cells were precisely estimated using the “estimate” R package to quantify the immune infiltration status of the tumor microenvironment. Furthermore, the TIDE (Tumor Immune Dysfunction and Exclusion) online tool (http://tide.dfci.harvard.edu/) was employed to analyze the immune evasion potential13. Differences in immune therapy response between the high-risk and low-risk groups were compared to investigate the potential link between riskScore and immune therapy prognosis.
Genetic variation characteristic
This study employed a multi-dimensional approach to investigate the genetic variation features associated with HCC. First, somatic mutation annotation and mutation feature analysis were conducted on the TCGA-LIHC dataset using the “maftools” R package. SNV waterfall plots were created to systematically illustrate the heterogeneity of gene mutation profiles across samples. Next, a gene co-mutation heatmap was used to identify co-occurrence and mutually exclusive patterns of high-frequency mutated genes (e.g., TP53, CTNNB1). Fisher’s exact test was applied to identify significant co-mutation events (P < 0.05).
Drug sensitivity prediction
This study integrated genomic-based drug sensitivity prediction models to systematically evaluate potential targeted therapies for HCC. Drug sensitivity prediction matrices were constructed using the “oncopredict” R package and the “pRRophetic” R package14,15. Spearman correlation analysis was employed to quantify the association between the drug half-maximal inhibitory concentration (IC50) and the riskScore. Chemotherapeutic drugs significantly correlated with the riskScore were identified based on Spearman correlation analysis (P < 0.001).
Statistical analysis
All statistical analyses were performed using R software (https://www.r-project.org/). Statistical significance was defined as a P value < 0.05. Correlation analyses between variables were conducted using Pearson’s correlation coefficient (for normally distributed data) and Spearman’s correlation coefficient (for non-parametric data). For continuous variables, independent sample t-tests (for normally distributed data) or Mann–Whitney U tests (for non-normally distributed data) were applied for pairwise comparisons, based on the distribution characteristics. Differences among multiple groups were analyzed using the Kruskal–Wallis test (non-parametric) or ANOVA (for normally distributed data with homogeneous variance). Survival analysis was conducted using the Kaplan–Meier method to visualize the prognostic differences between high- and low-risk groups, and the Log-rank test was used to assess the significance of survival differences.
Reagents and antibodies
The siRNAs targeting SLC16A3 and STRBP were constructed by Shanghai Genechem Co., Ltd. (Shanghai, China). All antibodies, including anti-MCT4 (SLC16A3) (22,787-1-AP), anti-STRBP (17,362-1-AP), anti-α-Tubulin (11,224-1-AP), and anti-GAPDH (10,494-1-AP), were purchased from Proteintech Group (Wuhan, China).
Cell culture and transfection
The human normal liver cell line MIHC and HCC cell lines Hep3B, PLC/PRF/5, Huh7, and MHCC97H were obtained from the American Type Culture Collection (ATCC). Cells were cultured in a 37 °C incubator with 5% CO₂, and cell morphology and confluency were monitored every 8–24 h. When cell confluency reached 80–95%, cells were passaged using 0.25% EDTA trypsin. After centrifugation (1000 rpm, 5 min), cells were resuspended in fresh complete medium for continued culture. Transfection was performed in 6-well plates: when the cell confluency reached 70–80%, siRNA was mixed with opt-MEM medium and transfected using the Lipofectamine 3000 reagent (Thermo Fisher Scientific) according to the manufacturer’s instructions. After 6–8 h, the medium was replaced, and subsequent experiments were conducted 48 h post-transfection.
WB analysis
Sample Preparation: Cells were lysed with pre-chilled RIPA lysis buffer (containing protease and phosphatase inhibitors) and incubated on ice for 30 min. The lysates were then centrifuged at 12,000 rpm for 15 min at 4 °C to collect the supernatant. Protein concentrations were measured using the BCA assay. The samples were mixed with 5 × loading buffer and boiled for denaturation.
Electrophoresis and Transfer: A 10% separating gel (acrylamide/Bis 29:1, Tris–HCl pH 8.8, 0.1% SDS, 0.05% APS, 0.1% TEMED) and a 5% stacking gel (acrylamide/Bis 29:1, Tris–HCl pH 6.8, 0.1% SDS, 0.05% APS, 0.1% TEMED) were prepared. The electrophoresis was carried out with an initial voltage of 80 V (stacking gel), followed by 120 V (separating gel). Proteins were transferred to methanol-activated PVDF membranes (Millipore) using a wet transfer system (Bio-Rad) at a constant current of 250 mA for 90 min.
Antibody Incubation and Detection: The PVDF membranes were blocked with 5% non-fat milk for 1 h, then incubated overnight at 4 °C with the primary antibody, followed by a 1-h incubation at room temperature with the HRP-conjugated secondary antibody. Protein bands were visualized using ECL chemiluminescence reagents, and the results were captured using a chemiluminescence imaging system. The images were saved for further analysis. GAPDH and α-Tubulin were used as internal controls in WB assays to normalize protein loading across samples.
EdU-488 cell proliferation assay and transwell invasion assay
EdU-488 Cell Proliferation Assay: Cells were seeded in 6-well plates and incubated until they reached 30%-40% confluence. The medium was then replaced with fresh medium containing 50 μM EdU reagent, and cells were incubated at 37 °C for 2 h. After removing the medium, the cells were fixed with 4% paraformaldehyde for 15 min, permeabilized with 0.3% Triton X-100 for 15 min, and incubated with Click reaction solution (containing Alexa Fluor 488) for 30 min in the dark. The cells were subsequently counterstained with Hoechst for 10 min. Fluorescence images were acquired from random fields of view under a fluorescence microscope.
Transwell Invasion Assay: Transwell chambers (Corning, 8 μm pore size) were pre-chilled to 4 °C. Matrigel was diluted in serum-free medium at a 1:8 ratio and applied evenly to the upper chamber (40 μl per chamber), followed by incubation at 37 °C for 1 h to allow solidification. Log-phase cells were trypsinized, resuspended in serum-free medium (2 × 104 cells/200 μl), and seeded into the upper chamber. The lower chamber was filled with medium containing 20% fetal bovine serum. After 48 h of incubation at 37 °C, the cells were fixed with 4% paraformaldehyde for 30 min and stained with 0.1% crystal violet for 15 min. Non-invasive cells on the upper surface of the membrane were gently removed with a cotton swab, and images of random fields were captured under a microscope.
