A total of 660 HCC samples datasets were procured from three publicly available datasets. Raw RNA sequencing data were standardized by variance-stabilizing transformation (VST) using the DESeq2 package in R, include 349 samples from the Cancer Genome Atlas (https://portal.gdc.cancer.gov/) TCGA-LIHC cohort (10) and 196 samples from the International Cancer Genome Consortium (https://dcc.icgc.org/) ICGC-LIRI-JP cohort (11). The microarray datasets, 115 samples of GSE76427, was downloaded from the Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/) (12). Genes associated with a favorable response to TTK in cancer immunotherapy were obtained from the TISIDB database (http://cis.hku.hk/TISIDB/) and used to established a gene set, referred to as GSTTKs (9).

GSTTKs differentially expressed between paracancerous and cancerous tissues were identified using the R package DESeq2 (13), with a false discovery rate < 0.05 and |Log fold change | > 1 as thresholds for significance. GSTTKs significantly associated with OS in HCC were identified by univariate Cox regression using the Survival package in R. A Venn diagram was generated using the VennDiagram package to identify the intersection of differentially expressed GSTTKs and prognostic GSTTKs. Somatic mutations in these genes in patients were described using the maftools R package (14). The copy number variation (CNV) status of each gene was retrieved from TCGA and delineated using GISTIC 2.0 to obtain chromosome information along with the gain or loss status, which was visualized in a circos plot (15). A principal component analysis (PCA) was performed using the PCAtools package in R to determine whether specific GSTTKs in the TCGA-LIHC dataset can distinguish between liver tumor samples and non-tumor samples.

The ConsensusClusterPlus package was employed for unsupervised clustering using the following parameter settings: partitioning around medoid (PAM) based on the center point, merge based on Ward’s distances using the minimum variance method (16). In addition, 1000 times repetitions were conducted for guaranteeing the stability of classification. The proportion of ambiguous clustering (PAC) was used to automatically select the optimal number of subtypes. PCA and tSNE analysis were performed to compare the transcriptional profiles between the different immune subtypes. For the clustering results for the TCGA-LIHC and ICGC-LIRI-JP cohorts, Kaplan–Meier survival curves were plotted and log-rank tests were performed using the survminer and survival packages in R.

Based on TCGA-LIHC dataset, a single sample gene set enrichment analysis (ssGSEA) was performed to quantitatively detect the relative levels of infiltration of 28 immune cells in the TME (17). The genetic signatures for these 28 immune cells were derived from Charoentong et al. (18). In the ssGSEA, differentially expressed marker genes were employed to evaluate the abundance of immune cells in individual samples. The relative abundance of each type of immune cell was represented as an enrichment score. To further explore the relationship between HCC subtypes and immune cell infiltration in HCC, the Wilcoxon rank sum test was used to analyze the differences in immune cell abundance between HCC subtypes. TIDE (Tumor Immune Dysfunction and Exclusion) algorithm developed by Liu can simulate the two main mechanisms of tumor immune escape: the induction of T cell dysfunction at high cytotoxic T lymphocyte (Cytotoxic T Lymphocytes, CTL) and the prevention of T cell infiltration at low CTL, and predict the response potential of tumor immunotherapy. This algorithm was used to evaluate the LIHC cohort, which is verified with the results of ssGSEA analysis to explore the difference of TME among different TTK patterns of HCC (19). The stromal and immune score was determined using the ESTIMATE package in R to assess the level of immune infiltration. These analyses were performed using the gene set variation analysis (GSVA) (20), ComplexHeatmap and estimate packages in R.

To evaluate the correlation between molecular subtypes and immune markers, the characteristic signatures related to differentially infiltrating immune cells in the HCC subtypes were collected from previous studies. Data for 148 immunomodulators and inhibitory immune checkpoints, including 41 chemokines, 21 major histocompatibility complex molecules, 18 receptor molecules, 44 immunostimulant molecules, and 24 inhibitory immune checkpoint molecules, were collected from previous studies (18, 21, 22). The Wilcoxon rank sum test was used to analyze the differential expression of these genes between the HCC subtypes. To determine the correlation between molecular subtypes and specific biological processes, annotated gene sets derived from the Kyoto Encyclopedia of Genes and Genomes (KEGG), Molecular Signatures Database (MSigDB), and a study by Mariathasan et al. were used for enrichment analysis and for comparing biological processes among subtypes (20, 23). For typing based on glycolipid metabolism, a glycolysis-cholesterol synthesis axis-related gene set was obtained from Schaeffer et al. (24) and modeling was completed using the R package ConsensusClusterPlus package. GSEA, GSVA, and over-representation analysis (ORA) were performed using the ClusterProfiler (25) and GSVA (20) packages in R.

Somatic mutation information in the mutect2 format for patients with HCC in TCGA-LIHC was converted to the mutation annotation format. The maftools package was used to generate waterfall diagrams to visually represent genes with high mutation frequencies. To investigate differences in the distribution of mutations among the HCC subtypes, differentially mutated genes were identified using p < 0.05 as the threshold for significance. Non-negative matrix factorization was carried out to reduce the dimensionality of the mutation matrix for the LIHC dataset, and the optimal number of mutation signatures in different HCC subtypes was identified (26). Thirty tumor mutational signatures that have been reported in COSMIC (https://cancer.sanger.ac.uk/cosmic) were downloaded for comparison with signatures identified by NMF, and mutational signature features of HCC were determined (27). A bar graph of 96 trinucleotide changes was generated for each sample to show the base change profile of each mutation feature. The whole process was performed using the NMF, BSgenome, and MutationalPatterns packages in R.

The index to represent the TTK level was establish based on the expression data for 18 GSTTKs including risk factors of CA9, SLC1A7, E2F1, RECQL4, AURKA, CENPF, RFPL4B, H2AFZ, KIF11, CDC7, TGIF2LX, MCM10, GRM4 and protective factors of SLC4A10, CAPN11, MYO1B, NR4A3, FGF12. The enrichment score (ES) of gene set that positively or negatively regulated TTK was calculated using single sample gene set enrichment analysis (ssGSEA) in the GSVA package (20), and the normalized differences between the ES of the risk factors minus protective factors was defined as the TTK potential index (TCscore) to computationally dissect the TTK trends of each sample:

The relationships between the TCscore and clinical characteristics, sensitivity to chemotherapeutics were evaluated. AJCC guidelines recommend the use of antineoplastic drugs such as doxorubicin, mitomycin, vincristine, cisplatin and sorafenib in the treatment of HCC. We predicted the chemotherapy response of each sample to these five drugs based on the GDSC database (the Genomics of Drug Sensitivity in Cancer, https://www.cancerrxgene.org/). The prediction process is realized by pRRophetic (28) packages in R.

All statistical analyses were conducted using R versions 3.6.3 and 4.0.2. For comparisons of continuous variables between two groups, normally distributed variables were evaluated using independent Student’s t-tests, and non-normally distributed data were analyzed using Mann–Whitney U tests (the Wilcoxon rank sum test). The chi-square test or Fisher’s exact test was used for comparisons of categorical variables between two groups. The relationships between gene expression levels were evaluated on the basis of Spearman correlation coefficients. Univariate and multivariate Cox analyses were used to identify independent prognostic factors. Receiver operating characteristic curves were plotted using the SurvivalROC package, and the area under the curve was used to evaluate the accuracy of the TCscore in predicting prognosis. The Rtsne package was used for a t-distributed stochastic neighbor embedding (tSNE) analysis. Two-sided p < 0.05 was the threshold for significance.

Comprehensive analysis of GSTTKs using multi-group data of TCGA-LIHC cohort. The result of difference analysis of transcriptome data shows that 92 of 641 GSTTKs were upregulated or downregulated in HCC, as shown in a volcano map in Figure 1A and a heatmap in Figure S1A. Univariate Cox regression analysis revealed that 125 out of 641 GSTTKs were related to prognosis in HCC. Taking the intersection of the two sets of genes, 37 GSTTKs simultaneously exhibited differential expression and prognostic value in HCC (Figure 1B). The univariate Cox analysis of 37 GSTTKs showed that 11 GSTTKs were protective factors with HR < 1 and 16 GSTTKs were risk factors with HR > 1 for HCC prognosis (Figure 1C). According to the genomic data of TCGA-LIHC, the top 10 oncogenic pathways and effects of HCC are shown in Figure S1B. The mutational landscape for the 37 GSTTKs is displayed in a waterfall plot in Figure 1D. Eighteen out of the 37 GSTTKs had a mutation frequency of >1% and were closely associated with progression or recurrence in HCC. Results of univariate cox regression analysis and differential analysis for 18 GSTTKs were shown in Table S1. As shown in Figure S1C, the co-occurrence of CA9 mutations and MCM10 mutations was significantly overrepresented in HCC. In addition, we detected widespread CNV in these 18 GSTTKs (Figure 1E). Copy number gains were most frequent, and RECQL4, CAPN11, and FGF12 showed extensive CNV amplification, whereas H2AFZ showed a copy number loss. The chromosomal locations of the 18 GSTTKs with CNV are shown in Figure 1F. HCC and non-tumor samples could be completely separated by the PCA (Figure 1G) based on these 18 GSTTKs with differential mRNA levels (Figure 1H), indicating high heterogeneity in the mutation status and expression of GSTTKs between normal and HCC tissues. Thus, GSTTK expression changes may play a crucial role in HCC occurrence and progression.

Based on RNA-seq data and clinical data for TCGA-LIHC, we identified four different patterns which show that comprehensive landscape of 18 GSTTKs interactions and their prognostic significance for HCC patients was depicted with the 18 GSTTKs network correlations (Figure 2A). The R package of ConsensusClusterPlus was used to classify patients with qualitatively different TTK patterns based on the expression of 18 GSTTKs, and two distinct modification patterns were eventually identified using unsupervised clustering, including 146 cases in Cluster1 and 203 cases in Cluster2 (Figure 2B). PCA algorithm and tSNE algorithm are used to evaluate the differences between the two TTK patterns, and it is found that there are significant differences in transcriptional profile among different TTK patterns (Figures 2C, D). To verify the stability and applicability of two TTK patterns in HCC, we repeated the unsupervised clustering analysis using LIRI-JP cohort from ICGC (Figure S2A) and GSE76427 cohort from GEO (Figure S2D); both populations could be well classified into two groups. The PCA (Figures S2B, E) and tSNE analysis (Figures S2C, F) results corroborated the two distinct patterns of TTK in HCC. Based on TCGA-LIHC expression profiling data, 16 out of 18 GSTTKs in the two clusters were significantly differentially expressed (Figure S2G). The clinical prognostic value of TTK patterns in patients with HCC was assessed through a survival analysis. Patients in the two clusters showed a significant difference in survival in both TCGA dataset (p = 0.0016, Figure 2E) and the ICGC dataset (p = 0.0025, Figure 2F).

By comparing the infiltrating immune cell composition in the TME of HCC between the two TTK subtypes (Figure 3A), we obtained the following key findings. 1) Samples in Cluster 1 mostly showed low immune cell infiltration, whereas samples in Cluster 2 mostly exhibited high immune cell infiltration. 2) The high immune infiltration zone in the heatmap contains immune cells that are established to mediate antitumor immune response (e.g., activated CD8+ T cells, type 1 T helper (Th1) cells, and dendritic cells) and multiple immunosuppressive cells (e.g., bone marrow-derived suppressor cells, regulatory T cells(Treg), immature dendritic cells, and neutrophils), suggesting that there may be a feedback mechanism, that is, TME may promote the recruitment or differentiation of immunosuppressive cells.

In order to determine the specific immune components that cause the difference of TME between Cluster 1 and Cluster 2, the differences of 28 immune cells among different subtypes were calculated. Combined with the results of the survival analysis (Figure 2E), samples in Cluster 2 corresponding to favorable survival outcome showed abundant infiltration by effector memory CD8⁺ T cells, Th1 cells, CD56 natural killer cells, eosinophils, natural killer T cells, neutrophils, and plasmacytoid dendritic cells, whereas those in Cluster 1 corresponding to an unfavorable clinical prognosis showed the infiltration of activated CD4⁺ T cells, effector memory CD4⁺ T cells, Th2 cells, and natural killer T cells (Figure 3B). A TIDE analysis based on RNA-sequencing data revealed that samples in Cluster 2 had higher scores for T cell dysfunction, microsatellite instability, and tumor-associated fibroblasts than those in Cluster 1, whereas samples in Cluster 1 scored higher for T cell exclusion, myeloid-derived suppressor cells, and tumor-associated M2 macrophages than those in Cluster 2; these findings were generally consistent with the ssGSEA results (Figure 3C).

Furthermore, we compared the two clusters with respect to biomarkers of infiltrating immune cells (Figure S3A) and molecular markers of the response to immunotherapy, including 41 chemokines, 21 major histocompatibility complex molecules, 18 receptor molecules, 44 immunostimulant molecules, and 24 inhibitory immune checkpoint molecules (Figure S3B). In addition, we also evaluated the correlations between the 18 GSTTKs and immune-infiltrating cells (Figure S4A) and identified significant correlations between NR4A3 and RECQL4 expression and most immune cells. Analysis of the differences in immune-infiltrating cells between the groups with high and low NR4A3 and RECQL4 expression were further analyzed (Figures S4B, C), and these results indicated that these genes may contribute to the difference between the TTK patterns.

Based on an enrichment analysis of TCGA-LIHC dataset by GSVA (Figure 4A), the two TTK clusters differed significantly with respect to metabolic pathways, suggesting that metabolic alterations as well as the TIME contributed to the distinct TTK patterns. Subsequent ORA (Figure S5A) and GSEA (Figure S5B) confirmed the difference in metabolic status between the two clusters. Next, we extracted glycolytic and cholesterogenic genes (Figure S6A) and used them to classify HCC into four metabolic subtypes: quiescent, glycolytic, cholesterogenic, and mixed (Figure 4B). PCA revealed a substantial separation among these four metabolic patterns (Figure S6B). The expression levels of genes involved in glycolipid metabolism are presented in Figure S6C. We detected significant differences in OS among the four metabolic clusters, with the quiescent and cholesterogenic subtypes being superior to the glycolytic and mixed subtypes (Figure 4C, p = 0.0032). This is consistent with the Warburg effect, in which aerobic glycolysis contributes to the aggressive cellular proliferation in malignant tumors. To investigate whether expression patterns across the glycolytic-cholesterogenic axis could underlie the differences between previously established immune subtypes (29), we determined the various HCC subtypes for each sample and investigated their degree of overlap with the metabolic phenotypes (Table S2). Quiescent and cholesterogenic subtypes could be classified into Cluster 2, whereas the glycolytic and mixed subtypes were mostly assigned to Cluster 1 or Lymphocyte Depleted Subtype (Figure 4D), suggesting that there is a relationship between TTK subtypes and the metabolic microenvironment in HCC. Analysis of the expression levels of the 18 GSTTKs (Figure S7A) and tumor-infiltrating cells (Figure S7B) according to the metabolic clusters uncovered the relationship between the immune and metabolic microenvironment in HCC.

We analyzed the distribution of somatic mutations in the two clusters using genomic data from TCGA-LIHC datasets (Figures 5A, B). Mutations in CTNNB1, a common therapy resistance gene in HCC, were predominant in Cluster 1, whereas mutations in TP53, a cardinal driver gene of HCC, were predominant in Cluster 2. Comparison of the mutant genes in the two clusters (Figure 5C) and revealed that both TP53 and RB1 showed the largest difference in mutation frequency between the two clusters. As somatic mutations are the result of multiple mutation processes, including DNA repair defects, and exposure to exogenous or endogenous mutagens, different mutation processes contribute to different combinations of mutation types or characteristics. To comprehensively characterize the landscape of genomic features, we identified five mutational signatures for the two HCC subtypes (Figure S8). C > A_DNA_Repair and DNA_MMR_Deficiency predominated in Cluster 1 (Figure 5D), whereas Smoking and DNA_MMR_Deficiency were the main patterns in Cluster 2 (Figure 5E).

We selected several markers of the stromal TME for a ssGSEA and found that the score for DNA damage response was significantly higher for Cluster 1 than for Cluster 2, whereas the scores for angiogenesis and epithelial interstitial transformation were significantly higher for Cluster 2 than for Cluster 1 (Figure 5F). To confirm these results, we calculated the stromal score as well as the ESTIMATE score for cases in TCGA-LIHC using the ESTIMATE algorithm and found significant differences between the two TTK types (Figures S9A, B). Moreover, there were significant differences in the stromal and ESTIMATE scores among the four metabolic subtypes (Figures S9C, D). Based on the differences in expression patterns and mutation frequencies between the two clusters, we selected genes involved in oncogenic pathways from the MsigDB and KEGG databases for a ssGSEA and found that only the receptor tyrosine kinase (RTK)-RAS pathway differed significantly between the groups (Figure 5G). Thus, RTK-RAS is the main pathway mediating the TTK patterns.

We developed a scoring system, referred to as the TCscore, to quantify TTK patterns based on the expression levels of the above 18 GSTTKs from TCGA-LIHC. A Spearman correlation analysis of the TCscore, stromal pathway score, and oncogenic pathway score revealed relationships between the TCscore and intertumoral or tumor microenvironment signaling pathways (Figure S9E). Most variables were negatively correlated with TCscore, but angiogenesis was most strongly related. TCscore were calculated for patients in the ICGC-LIRI-JP cohort using the formula applied to the TCGA-LIHC cohort to validate the prognostic ability of the GSTTKs signature. The sensitivity and specificity of the TCscore was assessed through time-dependent ROC analysis. The AUC values were 0.724 in TCGA-LIHC cohort (training set) and 0.729 in ICGC-LIRI-JP cohort (testing set), respectively (Figures 6A, B). A best threshold value of 0.738 was further selected for classification and Kaplan–Meier and cox regression analysis was performed after classification. Kaplan–Meier curves for OS were plotted according to the optimal cutoff value for TCscore for cases from TCGA-LIHC (p < 0.0001, Figure 6C) and ICGC-LIRI-JP (p < 0.0001, Figure 6E). Patients with high TCscore showed a relative shorter survival time than that of patients with low TCscore. The univariate (Figure S10) and multivariate cox regression analyses suggested that the TCscore is an independent prognostic factor for patients with HCC in TCGA-LIHC (p < 0.0001, hazard ratio (HR) = 1.910, 95% confidence interval (CI): 1.470–2.482, Figure 6D) and ICGC-LIRI-JP (p < 0.0001, HR = 2.136, 95% CI: 1.458–3.129, Figure 6F).

We combined the TCscore with other factors, including the mutation status of known oncogenes (e.g., TP53, ARID1A, AXIN1, CTNNB1, and TTN) and clinical characteristics (e.g., history of alcoholism, hepatitis B or C virus infection, expression of PDCD1, and tumor mutation burden), and plotted the Kaplan–Meier survival curves based on these parameters for cases from the TCGA database (Figure S11). When we also investigated the relationships between the TCscore and clinicopathological characteristics, we found significant correlations of the TCscore with survival, sex, T stage, grade, clinical stage, and TP53 and CTNNB1 mutation statuses (Figure S12). To assess the predictive value of the TCscore in immunotherapy, we used the TIDE algorithm to evaluate the associations with the treatment response and found that the differences in the TIDE score in view of the TCscore were similar to the TTK patterns (Figure 6G). A Spearman correlation analysis showed that the TCscore was negatively correlated with CAF cell (r = -0.4) and positively correlated with MDSC (r = 0.68, Figure 6H).

Using drug information from the GDSC database to calculate the half-maximal inhibitory concentration (IC50) values of common chemotherapeutics for HCC, we found that the IC50 value of sorafenib in the high TCscore group than in the low TCscore group (Figure 6I), whereas the IC50 values of four other drugs (doxorubicin, vinblastine, mitomycin, and cisplatin) showed the opposite pattern (Figures 6J–M), providing a basis for the selection of chemotherapy drugs when immunotherapy is combined with chemotherapy in clinical practice.