We screened patients from The Cancer Genome Atlas (TCGA) as the discovery set and patients from the International Cancer Genome Consortium (ICGC) as the validation set. The following enrollment criteria were used: 1) available sequencing data for GMScore calculation, 2) pathological diagnosis of HCC, and 3) no prior history of radiation therapy, chemotherapy, target therapy, immunotherapy, or other anticancer medications (including neoadjuvant therapy). Sequencing data and corresponding clinical information, recorded up to August 1, 2021, were downloaded from TCGA (https://portal.gdc.cancer.gov/repository) and ICGC (https://dcc.icgc.org/projects/LIRI-JP) data portals. An additional melanoma cohort (GSE78220) treated with anti-PD-1 was used to verify the association of GMScore with immunotherapy outcomes, and data were downloaded from Gene Expression Omnibus (GEO) [7]. Genetic expression data presented as the fragments per kilobase per million (FPKM) values were transformed into transcripts per kilobase million (TPM) values to improve comparability between samples [8]. The American Joint Committee on Cancer criteria was used for clinical and clinicopathological classification and staging.

The 41 GM-related genes (Supplementary Table S1) were extracted from the Gene Ontology (GO) initiative and a published study [9]. Of them, differentially expressed genes between cancerous and the adjacent normal tissues were identified by the limma R package in HCC patients of TCGA, with a false discovery rate of <0.05. Subsequently, genes associated with overall survival (OS) were selected using univariate Cox regression models. The optimal cut-off values to define high and low expression subgroups were determined based on the association of genetic expression with OS using the Survminer R package. Subsequently, gene expression level was evaluated as 0 or 1; a value of 0 was assigned when the gene expression was less than the corresponding cut-off value, and a value of 1 otherwise. Furthermore, the least absolute shrinkage and selection operator (LASSO) Cox regression model was used to screen the most useful prognostic genes. A GMScore model was then constructed based on the fraction of selected genes using Cox regression coefficients. The formula was established as follows: GMScore = sum (each gene’s expression × corresponding coefficient).

The ESTIMATE algorithm was used to estimate the relative fraction of stromal and immune cells in the tumor microenvironment (TME) and was exhibited in the form of the StromalScore, ImmuneScore, and ESTIMATEScore. Of these, ESTIMATEScore correlated with DNA copy number-based tumor purity [10].

The total counts of somatic nonsynonymous variations, including missense, nonsense, splice-site, and frameshift mutations, in coding regions were defined as TMB [11].

Differentially expressed genes between the low and high GMScore subgroups were determined, with the criteria of adjusted p-value <0.05 and log2 (fold change) >1, and GSEA based on GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed using NetworkAnalyst 3.0 [12].

The CIBERSORT algorithm and the LM22 gene signature were utilized to quantify the abundance of infiltrating immune cells in the TME, based on transcriptome data [13].

TIDE, a computational method based on transcriptome to model the induction of T cell dysfunction and the prevention of T cell infiltration in tumor immune evasion, was used to calculate T cell dysfunction and exclusion scores, which were further merged as the TIDE score to predict tumor response to immune checkpoint inhibitors (ICIs) [14].

R software (version 4.0.3, http://www.r-project.org) or IBM SPSS Statistics ver.20 (IBM Corp., Armonk, NY, United States) were used to analyze the related data and plot graphs. The Student’s t-test, chi-square test, Fisher’s exact probability test, or Mann-Whitney U test were used to compare the differences between groups. The Kaplan-Meier method with the log-rank test was used to compare OS between different parts. Hazard ratios (HRs) and their 95% confidence intervals (CIs) for prognostic factors were calculated using univariate and multivariate Cox proportional hazard models. Patients were separated into high and low GMScore groups based on the optimal value associated with OS calculated using the Survminer R package. Receiver operating characteristic curve (ROC), time-dependent ROC curves, and the areas under the ROC curves (AUC) depicted or calculated by the timeROC R package, were used to evaluate the predictive power of the models. Statistical significance was set at p < 0.05, and all p-values were two-tailed.

A total of 365 patients from TCGA and 231 patients from ICGC were included (Supplementary Table S2). Patients in ICGC were older and had more stage III/IV diseases than those in TCGA (p < 0.05). Histology grade, alpha-fetoprotein (AFP), and vascular invasion data were only available in TCGA.

A total of 30 GM-related genes were differentially expressed between HCC and normal tissues (Figure 1A). Of these, 20 genes were significantly associated with OS (Figure 1B). After LASSO Cox regression analysis (Figures 1C,D), seven genes were selected to construct the GMScore of OS, as follows: GMScore = 0.374 ∗ expression level of SLC1A5 + 0.359 ∗ expression level of GAPDH +0.264 ∗ expression level of SLC38A1 + 0.112 ∗ expression level of SLC38A7 − 0.049 ∗ expression level of FTCD-0.113 ∗ expression level of MTHFS-0.157 ∗ expression level of GOT2.

Patients with high GMScores in both TCGA and ICGC included a significantly higher proportion of patients with stage III/IV disease (p < 0.001 and p = 0.009, respectively) than those with low GMScores. Histology grade 3/4 (p < 0.001), AFP larger than 200 ng/ml (p < 0.001), and vascular invasion (p = 0.001) were also more frequent in patients with high than in those with low GMScores (Table 1).

Patients were divided into high and low GMScore subgroups to establish a prognosis-predictive model in TCGA, which was validated using ICGC (Figure 2A). Kaplan-Meier survival analysis confirmed the survival discrepancy between high and low GMScore groups (Figure 2B), and the results of the ROC curve analysis verified the predictive value of the established risk model in both TCGA and ICGC (Figure 2C). After univariate selection for prognostic significance of variables (Supplementary Figure S1), multivariate analysis in both TCGA and ICGC showed that high GMScore was an independent predictor of poor OS (HR = 4.2, 95% CI 2.38–7.4, p < 0.001, and HR = 3.91, 95% CI 1.92–7.97, p < 0.001, respectively), even superior to staging in terms of HR (Figure 2D).

ImmuneScore, StromalScore, ESTIMATEScore, and TMB have been reported as prognostic predictors in many cancers [10, 11, 15]. We showed that they were also associated with or tended to be associated with OS of HCC (Supplementary Figure S2). However, prognostic stratification according to these factors could be further optimized using GMScore. The OS of both patients in the high and low subgroups of these biomarker values could be further stratified by the GMScore level (Figures 3A,B). Moreover, the ROC curve analysis revealed that GMScore was better than other biomarkers in terms of the prediction of 3-years OS (only two patients were followed up for over 5 years in ICGC) (Figures 3C,D).

Differentially expressed genes between the high and low GMScore subgroups were determined (Figure 4A and Supplementary Table S3). GSEA based on the GO biological process (BP) and KEGG pathways was performed (Supplementary Table S4). We revealed that genes associated with proliferation and cell cycle were significantly enriched in the high GMScore subgroup, indicating that tumors with high GMScores might have cell growth advantages (Figure 4B). Meanwhile, high GMScore also enriched genes associated with DNA repair and chromatin regulation, suggesting that tumors with high GMScores may be characterized by genomic stability (Figure 4C).

It is known that genomically stable or elevated DNA repair correlates with therapeutic resistance; thus, we examined the impact of high GMScore on outcomes of transcatheter arterial chemoembolization (TACE). In the TCGA HCC cohort, high GMScore significantly decreased OS in patients who received adjuvant TACE after surgery (p = 0.003; Figure 4D) or salvage TACE after recurrence (p < 0.001; Figure 4E) than those with low GMScore. A similar result was found in the ICGC HCC cohort (Figure 4F), although the timing of TACE was unknown.

TMB has been identified as a predictor of ICI efficacy in multiple cancers. In this study, we found that GMScore level did not affect TMB (Figure 5A). However, high GMScore subgroup had significantly higher gene expressions than low GMScore subgroup in terms of several immune checkpoints, including PD-1, CTLA-4, TIM-3, and TIGIT in both TCGA and ICGC (Figure 5B). For immune cell infiltration in tumors, high GMScore subgroup had a significantly higher abundance of regulatory T cells (Tregs) but a significantly lower abundance of M1 macrophages than low GMScore subgroup (Figure 5C). These findings indicate that GMScore may partly reflect TME.

Because no HCC cohorts treated by ICIs with transcriptome data were available, we used the TIDE algorithm to calculate T cell dysfunction and exclusion scores, and predict ICI response by the merged TIDE score. We observed that high GMScore subgroup had significantly higher exclusion score and TIDE score compared to the low GMScore subgroup (Figure 5D). Consequently, for the predicted response rate, high GMScore HCC was significantly inferior to low GMScore HCC in both TCGA and ICGC (Figure 5E). Furthermore, we used an ICI-treated melanoma cohort to verify the impact of GMScore on immunotherapy outcomes and found that the high GMScore subgroup had a significantly lower objective response rate than the low GMscore subgroup (14.3 vs. 66.7%, p = 0.016; Figure 5F). More importantly, high GMScore significantly reduced OS compared to low GMScore (p = 0.039; Figure 5G).