Figure 1 presents the research flowchart . Supplementary Table S1 was the List of Abbreviations. RNA-seq, mutation data, and clinical information were obtained from The Cancer Genome Atlas (TCGA) (424 HCC cases; https://portal.gdc.cancer.gov/v1) as the training set. We obtained 6908 exercise-related genes from GeneCards (search term: exercise; score: 1; www.genecards.org; Supplementary Table S2 ). HCC RNA-seq and clinical information obtained from the International Cancer Genome Consortium (ICGC) were used as the external validation set (233 HCC cases; https://docs.icgc-argo.org/docs/data-access/icgc-25k-data).

Weighted Gene Co-expression Network Analysis (WGCNA) is a method that constructs scale-free gene co-expression networks to identify functional modules significantly associated with clinical phenotypes and uncover key regulatory genes. In our study, WGCNA was performed to obtain genes most significantly associated with HCC. We obtained gene expression data from TCGA and removed genes with standard deviation of expression level < 0.5 to construct a scale-free co-expression network. The clustering dendrogram was cut at a height of 200 to detect outlier samples. Scale independence and mean connectivity were used to determine the optimal soft threshold, with a scale-free R² of 0.9. The adjacency matrix was converted into a topological overlap matrix (TOM; T O M i , j = ∑ u a i u a u j + a i j min k i , k j + 1 − a i j ) and calculated the corresponding dissimilarity (1-TOM). In the TOM formula, the numerator accounts for the direct connection weights and the indirect association strength mediated by all common neighbors, while the denominator acts as a normalization term. To convert TOM into a distance metric, the dissimilarity between two genes is defined as 1 − T O M x , y , where T O M x , y quantifies the similarity between gene x and the gene y based on their shared topological connectivity patterns. Subsequently, dynamic tree cutting and hierarchical clustering were used to identify modules, with the minimum module size of 50 genes, and the dendrogram was cut at a height of 0.1 to define module clusters.

We conducted differential gene expression analysis on the module genes obtained from WGCNA analysis with |log2 fold change (FC)| > 1 and FDRq < 0.05. Then, we used CIBERSORT (https://CIBERSORT.stanford.edu/) to assess the infiltration abundance of 22 immune cells in HCC patients. Finally, Spearman analysis was employed to assess the correlation between exercise-related differentially expressed genes (DEGs) and the abundance of 22 immune cells with |r| >0.4.

We performed univariate COX analyses to obtain EIGs that were significantly associated with survival in the training set for the construction of EIG-related prognostic signature (EIGPS).

We used 10 machine learning algorithms — LASSO, Ridge, StepCox, CoxBoost, survival support vector machine (survival-SVM), supervised principal components (SuperPC), random survival forest (RSF), generalized boosted regression modeling (GBM), elastic net (Enet), and partial least squares regression for Cox (plsRcox) — to screen variables based on 10-fold cross-validation to construct 101 models. Models containing less than 5 genes were excluded. The optimal model was selected based on the highest C-index. Using the linear combination of the optimal model, we calculated a risk score for each HCC patient, and subsequently divided them into high-risk group (HRG) and low-risk group (LRG). Furthermore, since CoxBoost+SuperPC was the subsequent optimal model, we briefly introduce these two algorithms here. CoxBoost integrates the Cox proportional hazards model with boosting algorithms. Boosting can iteratively optimize regression coefficients to yield a sparse model. SuperPC combines supervised learning with PCA. It first selects features significantly associated with survival outcomes, then applies PCA to reduce dimensionality and extract principal components (PCs) as new features in the subsequent model. This strategy can address high-dimensional collinearity and enhance model robustness.

We evaluated the model performance using the training set and the validation set. Principal component analysis (PCA) was used to reduce the dimensionality of gene expression data and project them onto the principal components to identify the distribution of the feature space of HRG and LRG after dimensionality reduction. The K-M survival curve, clinical ROC curve, temporal ROC curve, univariate and multivariate COX analysis, and risk curve were obtained using the “Survival”, “survminer”, and “timeROC” R package. The significance threshold was p < 0.05, and the area under the curve (AUC) of ROC >0.5 was considered significant.

We used the chi-square test to assess clinical features in the HRG and LRG, and further conducted survival analysis to examine the associations between clinical features and EIGPS.

Based on the expression of EIGPS genes, we performed clustering analysis to identify molecular subtypes using the “ConsensusClusterPlus” R package with Partitioning Around Medoids (PAM) algorithm and Euclidean distance; 80% of the samples were resampled for 10 repetitions. The optimal number of clusters (k) was determined by cumulative distribution function (CDF) plot and average cluster consensus.

To further investigate the potential functional and pathway differences between HRG and LRG, we used the “clusterprofiler” R package for GO and KEGG enrichment analysis of DEGs. In addition, gene set variation analysis (GSVA) was used to identify differentially enriched pathways.

For the sake of individualized prediction, we developed two nomograms for HCC patients. The first nomogram combined risk scores and clinical features, while the second utilized the signature genes. The C-index value was calculated to evaluate the consistency of the predicted values with the observed values.

We employed CIBERSORT to assess differences in the infiltration abundance of 22 immune cells between HRG and LRG. We used single-sample gene set enrichment analysis (ssGSEA) via “GSVA” R package to assess differences in enrichment levels for 29 immune traits (16 immune cells, 13 immune functions) between HRG and LRG. Additionally, TIMER, CIBERSORT−ABS, QUANTISEQ, MCPCOUNTER, XCELL and EPIC were used as additional immune infiltration algorithms to further validate the results.

Earlier research reported that all tumors could be divided into six immune subtypes, namely, wound healing (C1), IFN-γ dominant (C2), inflammatory (C3), lymphocyte depleted (C4), immunologically quiet (C5), and TGF-β dominant (C6) (7). Therefore, we performed survival analysis based on different subtypes.

To evaluate the association of EIGPS with immunotherapy, we used TISIDB (https://cis.hku.hk/TISIDB/index.php) to acquire immune checkpoint genes and assess the differences between HRG and LRG. The Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/) was used to score the immune evasion ability of HRG and LRG, with higher scores indicating worse response to immunotherapy. The immunophenoscore (IPS) of four immunotherapy regimens for HCC were obtained from the TCIA database (https://tcia.at/home), with higher IPS indicating better response to immunotherapy.

Macrophage-related genes were derived from Genecards ( Supplementary Table S3 ). TMB was quantified by the number of mutations per megabase. Using the optimal cutoff obtained by “survminer” algorithm, the samples were divided into high TMB group and low TMB group, and were then combined with the risk score for survival analysis.

The scRNA-seq data for HCC were obtained from the China National Gene Bank Nucleotide Sequence Archive (CNSA: CNP0000650; https://db.cngb.org/cnsa) and GEO (GSE162616), including a total of 15 HCC cases.

We used “Seurat” R package to process scRNA-seq data and obtained 34472 genes and 54982 cells. Subsequently, PercentageFeatureSet and FeatureScatter functions were used to calculate the percentage of mitochondrial genes and sequencing depth. We excluded cells with fewer than 20 mitochondria, cells with number of genes < 500 genes or > 10,000, and cells with UMI counts < 500 or > 20,000. The “Harmony” package and PCA were used for batch effect removal and dimension reduction, respectively. “FindNeighbors” and “FindClusters” functions were used for cell clustering and determining the resolution, and Uniform Manifold Approximation and Projection (UMAP) was then used for visualization. Cell types were manually annotated using marker genes. Finally, we performed cell-cell communication analysis using the “CellChat” package, with ligands and receptors pairs from the CellChat database (http://www.cellchat.org/).

Based on the GDSC database, we used the “oncoppredict” R package to evaluate the sensitivity of HRG and LRG to 198 FDA-approved drugs. Drug sensitivity was evaluated using half maximal inhibitory concentration (IC50), with lower values indicating higher drug sensitivity.

The cell lines (L-O2, HepG2) were purchased from the Institute of Cell Research, Chinese Academy of Sciences. Total RNA was extracted using Trizol (Invitrogen, 1596-026) according to the manufacturer’s protocols. cDNA was synthesized using reverse transcription kit (Fermentas, #K1622). qRT-PCR was performed using SYBR Green kit (Thermo, #K0223). Primers are shown in the Supplementary Table S8 .

We obtained the immunohistochemical results of the biomarkers in the normal population and HCC patients from the HPA database (https://www.proteinatlas.org/).

All statistical analyses were performed using R software (Version 4.3.2; the R Foundation, St. Louis, MO, USA). Chi-square test was used for categorical data, while t-test or Wilcox test was used for continuous data. Unless otherwise specified, p < 0.05 was the significance threshold.

After acquiring exercise-related genes from Genecards, we performed WGCNA using the training set to identify module genes most associated with HCC prognosis. We identified 8 co-expression modules using hierarchical clustering with an optimal power value of 10 ( Figures 2A, B, D ; Supplementary Table S4 ). Among these modules, the turquoise module (500 genes) had the highest correlation with HCC. And the gene significance (GS) and module membership (MM) for the turquoise module showed a significant correlation (r=0.66, P=7.5e-64; Figure 2C ). We then performed differential gene expression analysis on 500 genes in the turquoise module and identified 343 DEGs ( Figure 2E ; Supplementary Table S5 ).

Subsequently, we assessed the correlation between the identified DEGs and the infiltration abundance of 22 immune cells, and obtained 59 EIGs ( Supplementary Table S6 ).

We conducted univariate COX analysis to further identified 54 EIGs that were significantly associated with HCC prognosis ( Figure 2F ; Supplementary Table S7 ).

Among the 101 models constructed through machine learning, the CoxBoost +SuperPC model had the highest average C-index (0.742; Figure 3G ), and incorporated 7 EIGs (UPF3B, G6PD, ENO1, FARSB, CYP2C9, DLGAP5, SLC2A1). Low expression of CYP2C9 was associated with poor HCC prognosis, whereas the other genes were the opposite ( Supplementary Figures S1C–I ). qRT-PCR and immunohistochemical results also supported these findings ( Figure 3F ; Supplementary Table S8 ; Supplementary Figures S2A–G ).

To comprehensively assess the robustness of the constructed EIGPS, we calculated a risk score for each patient and divided all patients into HRG and LRG. Univariate (HR=2.229, 95%CI= 1.773−2.803, P < 0.001; Supplementary Figure S1B ) and multivariate (HR=1.978, 95%CI= 1.527−2.561, P < 0.001; Supplementary Figure S1A ) COX regression analyses showed that the risk score was associated with HCC survival independent of other clinical features.

PCA suggested that the risk score could significantly distinguish between HRG and LRG ( Figure 3E ). The K-M curves showed that the survival of LRG was significantly higher than that of HRG ( Figure 3B ). And the mortality gradually increased with the increase of risk scores ( Figure 3A ). The area under the ROC curves and decision curves suggested that the EIGPS had better predictive performance (AUC=0.804) than other clinical features ( Figure 3C ). The time-ROC indicated that EIGPS had good predictive performance in predicting survival at 1 (AUC=0.804), 3 (AUC=0.726), and 5 (AUC=0.696) years ( Figure 3D ).

We further validated the performance of EIGPS in the validation set, and observed that the results were all consistent with those in the training set ( Supplementary Figures S3A–I ).

Based on the expression of EIGPS, we divided HCC samples into 2 subtypes (optimal k=2; Supplementary Figure S4A ; Supplementary Figures S4A, B ). PCA showed a clear distinction between the two subtypes ( Figure 4C ). The expression patterns of EIGPS in the two subtypes were similar to risk scores ( Figure 4F ). The C2 subtype corresponded mainly to the LRG and had higher survival rates than C1 ( Figure 4E ). As for tumor immune microenvironment, C2 had significantly lower immune cell infiltration abundance ( Figure 4D ; Supplementary Figure S4C ), whereas C1 showed a large number of immune tumor-promoting phenotypes, especially higher macrophage (Mφ) abundance across 7 immune infiltration algorithms ( Figures 4B, G ).

We observed significant differences in T-stage, Stage, and Grade between HRG and LRG ( Supplementary Figure S5D ). We combined EIGPS with the significantly different features, and found that LRG combined with any clinical feature in early or late stages had better survival than HRG ( Supplementary Figures S5A–C ).

To optimize the clinical application of the risk model, we developed a nomogram based on clinical features and risk scores ( Figure 5H ), which could effectively predict 1-year, 3-year, and 5-year survival (C-index=0.754, 95%CI: 0.703-0.806; Figure 5G ). We also constructed a nomogram based on EIGPS gene expression (C-index=0.711, 95%CI: 0.664-0.759; Figures 5C, D ), which may benefit HCC patients who did not have the sequencing of all the 7 signature genes. Both nomograms showed independent predictive ability ( Figures 5A, E ). The mortality gradually increased with the increase of nomogram scores ( Figures 5B, F ).

Based on the immune subtypes identified by the previous study, we found 4 immune subtypes for HCC patients ( Supplementary Figure S6A ), with the C1 and C2 subtypes mainly found in the HRG ( Supplementary Figure S6B ). Patients with C1 and C2 had a significantly worse prognosis compared to other subtypes ( Figure 6F ).

Subsequently, 7 immune infiltration algorithms all suggested a positive correlation between risk scores and macrophage abundance ( Figure 6B ). CIBERSORT results showed that macrophages M0 (Mφ) had significantly higher abundance in HRG than in LRG, while monocytes were the opposite ( Figures 6A, D ). Additionally, the additional 6 immune infiltration algorithms also noted the significant differences in macrophage abundance between HRG and LRG ( Figure 6E ). Similar results were obtained by ssGSEA that the macrophage abundance as well as a number of immune tumor-promoting functions was significantly elevated in HRG ( Figure 6C ).

Based on PCA, the elbow method (selecting the top 10 principal components), and a resolution of 0.8, we obtained 25,664 genes and 50,016 cells, identifying 17 clusters. ( Supplementary Figures S7A, D ; Figure 7F ; Supplementary Figures S7B, C ). Subsequently, based on the expression of marker genes, we identified 5 cell types, including B cells, epithelial cells, macrophages, NK cells, plasma cells, and T cells ( Figures 7D, E ). EIGPS were mainly expressed in macrophages ( Figure 7A ). Furthermore, the interaction between macrophages and epithelial cells exhibited the highest interaction number and strength ( Figures 7B, C ). Macrophages mainly functioned as ligand-secreting cells to send signals to epithelial cells, a process potentially mediated by PPIA-BSG ligand-receptor pairs ( Supplementary Figure S8 ).

As expected, the expression of immune checkpoint-related genes and chemokine-related genes in the HRG was significantly higher than those in the LRG ( Figures 8A, B ). Among them, we selected immune checkpoint-related genes that were closely associated with clinical treatment, and found that all of them were positively correlated with risk scores ( Figure 8E ). The HRG had a significantly higher TIDE score than LRG, indicating its higher immune evasion ability and poorer response to immunotherapy ( Figure 8C ). The IPS results revealed that LRG had a higher IPS for CTLA4-/PD-1- treatment ( Figure 8D ). Drug sensitivity analysis showed that the HRG was more sensitive to Afatinib and Alpelisib ( Figure 8F ), whereas the LRG was more sensitive to Irinotecan and Oxaliplatin ( Figure 8G ), which provides reference for preferable choices of drugs for different patients.

We further revealed differences in potential biological functions and pathways between the HRG and LRG. GO-BP analysis revealed that DEGs of HRG and LRG were mainly enriched in leukocyte mediated immunity and positive regulation of leukocyte activation (top 2; Supplementary Table S9 ; Supplementary Figure S9A ). KEGG results showed that the top 2 most involved pathways were cytokine−cytokine receptor interaction and human T−cell leukemia virus 1 infection ( Supplementary Table S10 ; Supplementary Figure S9B ). GSVA enrichment analysis showed that the HRG and LRG had significant differences in complement and coagulation cascades, Fcγ receptor mediated phagocytosis and other immune-related pathways ( Supplementary Figure S9C ).

HCC patients were divided into high TMB and low TMB groups using the best cutoff value (cutoff = 1.316). The K-M curves showed that high TMB group had significantly lower survival rate than low TMB group ( Supplementary Figure S10A ). Furthermore, LRG combined with either low or high TMB had better survival than the HRG ( Supplementary Figure S10B ).

Since EIGPS was mainly expressed in macrophages, we further analyzed 1528 macrophage-related genes, and found that TP53 had highest mutation frequency in the HRG, while the mutation frequency of CTNNB1 was highest in the LRG ( Supplementary Figures S10C, D ).