The transcriptomic data of TCGA-LIHC dataset and clinical information of 371 HCC patients and 50 non-tumor samples were obtained from The Cancer Genome Atlas (TCGA) data portal up to 22 July 2022 (https://portal.gdc.cancer.gov/repository). The raw count data were normalized by the “limma” R package. The transcriptomic data and clinical information of ICGC-LIRI-JP dataset were downloaded from the HCCDB database (http://lifeome.net/database/hccdb/download.html), which included 212 HCC tumor samples and 177 non-tumor samples (Lian et al., 2018). Normalized read count values were used. Additionally, a third transcriptomic data and clinical information of GSE116174 dataset including 64 HCC patients were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/gds/). The TCGA cohort was used as the training set, and the ICGC and GSE116174 cohorts were as the validation sets. Another GEO cohort of GSE104580 was obtained to investigate the signature association with transarterial chemoembolization (TACE) treatment. A total of 1,639 TFs reported by Lambert et al. were included for signature establishment (Lambert et al., 2018). Since all the data were publicly available, it was not required for additional ethical approval.

Differentially expressed genes (DEGs) were identified by the “limma” R package, with |log2FC| > 1 and an adjusted p-value <0.05 in the TCGA and ICGC cohorts. Those overlapped with the TF list were DETFs, which was visualized by a Venn diagram drawn by an online tool (https://bioinformatics.psb.ugent.be/webtools/Venn/). The expression profile of DETFs in the TCGA cohort was plotted using the “pheatmap” R package.

The prognostic significance of individual DETFs was investigated by univariate Cox analysis based on the TCGA cohort, where those with follow-up time more than 0 day being included. The correlation between prognosis-related TFs was further analyzed by the “corrr” R package. Then, the Least Absolute Shrinkage and Selection Operator (LASSO) penalized Cox proportional hazards regression was applied to construct an optimal model via the “glmnet” R package (Friedman et al., 2010). TFs with non-zero coefficients were selected with 10-fold cross-validation. The risk score of each patient was determined by the following formula: risk score = (expression level of TF 1×coefficient) + (expression level of TF 2 × coefficient) + … + (expression level of TF n × coefficient). The patients were classified into high-risk and low-risk groups according to the optimal threshold of the risk score, which was determined by the “survminer” R package. To estimate the predictive accuracy of the TF signature, the area under the curve (AUC) value of receiver operating characteristic (ROC) curve was calculated via the “survivalROC” R package. Univariate followed by multivariate Cox regression analyses were performed to determine the independence of TF signature for the prognosis of patients using the R package “survival”, which was plotted by the R package of “forestplot”. Subsequently, the predictive value of TF signature was further confirmed in the validation cohorts of ICGC-LIRI-JP and GSE116174.

To infer the potential therapeutic agents, the “pRRophetic” package in R was used to predict the half-maximal inhibitory concentration (IC₅₀) of numerous drugs for HCC patients based on the TCGA and ICGC cohorts, by constructing the ridge regression model according to the expression profile of TF signature and the cell line expression spectrum in the Genomics of Drug Sensitivity in Cancer (GDSC) (www.cancerrxgene.org/) (Geeleher et al., 2014).

The DEGs between the high- and low-risk groups were identified by the “limma” R package with |log2FC| ≥ 1 and an adjusted p-value <0.05. Then, the DEGs underwent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses using the “clusterProfiler” R package, where adjusted p values with the BH method were used.

The common DEGs between groups in both the TCGA and ICGC cohorts were uploaded to the STRING online database (https://string-db.org/) to generate a protein-protein interaction (PPI) network. Those with an interaction score above 0.4 were selected. The PPI network was visualized by the Cytoscape software (version 3.9.1). The MCODE plug-in was used to screen functional modules. The top 10 hub genes were selected by the MCC method via the cytoHubba plug-in. Additionally, the interactions between the hub genes were validated by the “corrplot” R package.

Five patients with primary HCC were retrospectively collected at the Mengchao Hepatobiliary Hospital of Fujian Medical University, China. Formalin-fixed and paraffin-embedded (FFPE) tissues from tumor and adjacent non-tumorous specimens were obtained to determine the transcriptional expressions of TFs in the signature. The study protocol was approved by the Ethics Committee of our institution. All patients provided informed consent. Total RNA of the FFPE tissues was isolated by using the FFPE RNA Kit (AmoyDx, China), and subjected to reverse transcription using the FastKing-RT SuperMix (Tiangen, Beijing, China). Gene expression was then assessed by RT-qPCR using the GoTaq qPCR Master Mix (Promega, WI, United States) and the ABI 7500 Real-Time PCR system (Thermo Fisher Scientific Inc., MA, United States), with GAPDH used as a housekeeping gene. The temperature protocol consisted of 95°C for 1 min for initial denaturation followed by 45 cycles at 95°C for 10 s and at 60°C for 30 s. The relative expression levels were calculated using the 2^−ΔΔCT method and normalized by GAPDH. The primer sequences were as follows: HMGA1 forward: 5′-AAG​GGG​CAG​ACC​CAA​AAA-3′, reverse: 5′-TCC​AGT​CCC​AGA​AGG​AAG​C-3’; MAFG forward: 5′-TCT​AGG​GCT​TGG​GCT​GAT​CT-3′, reverse: 5′-TGT​AGC​CCT​TGT​CTG​CAC​TG-3’; GAPDH forward: 5′- CAG​GAG​GCA​TTG​CTG​ATG​AT-3′, reverse: 5′-GAA​GGC​TGG​GGC​TCA​TTT-3’.

The FFPE blocks of tumor and adjacent tissues were used for IHC staining by applying an ElivisionTM plus Polyer HP Kit (Maixin Biotechnologies, Inc., Fuzhou, China). Freshly cut slides (4-μm thick) underwent deparaffinage, hydration, blockage of endogenous peroxidase activity with 3% hydrogen peroxide, and heat mediated antigen retrieval. The sections were then incubated with anti-HMGA1 (ab129153; 1:100 dilution; Abcam, United States) and anti-MAFG (ab154318; 1:100 dilution; Abcam) antibodies at 4°C overnight, followed by incubation with a polymer helper reagent and poly-peroxidase-anti-mouse/rabbit IgG. Finally, the sections underwent diaminobenzidine (Maixin Biotechnologies, Inc.,) staining, dehydration, and mount. Additionally, the protein expressions of TFs in the risk score in HCC and normal hepatocellular tissues were obtained from the Human Protein Atlas (HPA, https://www.ptroteinatlas.org/) which provided the IHC data of human tissues (Uhlén et al., 2015). The staining intensity was scored as 0 (negative staining), 1 (weak staining), 2 (moderate staining) or 3 (strong staining). The extent of positive staining was scored as 1 (<25%), 2 (25%–75%), 3 (>75%). The final score was obtained by multiplying the intensity and extent scores.

All the statistic work was accomplished by the R software (version 4.1.2). Overall survival (OS) was calculated by the Kaplan-Meier method, with the difference between groups compared using the log-rank test. The paired t-test was used to compare the transcriptional and protein differences in HCC and paired non-tumorous tissues in our cohort. The comparison of other continuous data was performed by the Wilcoxon test between two groups and Kruskal–Wallis test among three or more groups. A p-value < 0.05 was considered statistically significant.

Based on the TCGA and ICGC cohorts, a total of 3,119 and 942 DEGs were identified respectively in tumor tissues compared with the respective non-tumorous tissues. Among the overlapping differentially expressed ones in the TCGA and ICGC cohorts, 25 were TFs (Figure 1A). According to the expression profile (Figure 1B), 14 DETFs were upregulated and 11 downregulated in tumor tissues in the TCGA cohort. The results of univariate analysis demonstrated that 16 of the 25 DETFs were associated with the prognosis of patients. ESR1, TCF21 and FOSB, which were downregulated in HCC, predicted a favorable outcome. The other two downregulated ones and 11 upregulated TFs in HCC indicated an adverse prognosis (Figure 1C). According to the correlation network, most of the TFs were positively related, except for TBX15, ESR1 and FOSB, which had negative correlations with some TFs (Figure 1D).

The 16 prognosis-related TFs were included to establish a prognostic model in the training set of TCGA cohort via LASSO Cox regression analysis. A 2-TF-based signature, comprising high mobility group AT-hook protein 1 (HMGA1) and MAF BZIP transcription factor G (MAFG), was identified based on the optimal value of λ (Supplementary Figure S1). The transcriptional (Figure 2A) and protein (Figures 2B, C) levels of HMGA1 and MAFG were validated to be elevated in HCC specimens compared to those in adjacent normal tissues in our retrospective cohort. High or medium staining intensity of protein expression of the two TFs by IHC was also observed in HCC tissues based on the HPA database (Supplementary Figures S2A, B). The insignificant difference of MAFG expression between tumor and normal tissues (Supplementary Figure S2B) may be attributed to the antibody specificity. The risk score = 0.0924 × HMGA1 + 0.0658 × MAFG. Then, the risk score for each patient was determined accordingly in the TCGA cohort as well as in the validation sets of ICGC and GSE116174 cohorts. The AUCs for the 2-TF signature were 0.73 for 1-year, 0.61 for 2-year, and 0.60 for 3-year OS in the TCGA cohort (Figure 3A). In the validation sets, the AUCs were 0.74, 0.70, 0.74 in the ICGC cohort (Figure 3B), and 0.59, 0.62, 0.60 in the GSE116174 cohort for 1-, 2- and 3-year OS, respectively (Figure 3C). A relative low predictive efficiency of the TF signature in the GSE116174 cohort may be attributed to a small sample size. The patients were then divided into high- or low-risk group based on the optimal cut-off. As the risk score increased, the number of deaths increased and the survival time decreased, with the expressions of the 2 TFs also in an increased trend (Figures 3A–C). According to the Kaplan-Meier curve, patients in the high-risk group suffered a significantly more unfavorable OS than their low-risk counterparts in the TCGA cohort as well as in the ICGC and GSE116174 cohorts (Figures 3A–C). Collectively, our results indicated that the 2-TF signature had a good performance for predicting the OS of HCC patients.

In the TCGA cohort, 245 patients with complete information including age, gender, tumor grade, TNM stage, vascular invasion and alpha-fetoprotein (AFP) were included for multivariate Cox regression analysis. The results demonstrated that the risk score was the only factor significantly associated with OS in the TCGA cohort (HR = 2.498, 95% CI = 1.280–4.876, p = 0.007) (Figure 4A). To validate our results, 203 patients in the ICGC cohort with age, gender, TNM stage, invasion (portal vein, vein or artery) and fibrosis underwent multivariate analysis. As shown in Figure 4B, the risk score was still capable of independently predicting OS (HR = 5.411, 95% CI = 2.467–11.868, p < 0.001). Therefore, the 2-TF signature could be an independent indicator for an adverse outcome in HCC patients.

To comprehensively evaluate the involvement of the 2-TF signature in the tumorigenesis of HCC, we further explored the associations of the risk score with clinicopathological features. In the TCGA cohort (Figure 5A), G3 or G4 tumors tended to have a higher risk score compared to G1 tumors. The risk score between G2 and G3 tumors was also significantly different. Regarding TNM stage, the risk score increased as the stage I increased to stage II or III. For patients with vascular invasion, either micro- or macro-invasion cases, the risk score was higher compared to those without invasion. In the ICGC cohort (Figure 5B), similar relationship of the risk score with TNM stage and invasion were observed. There was no significant correlation between the TF signature and age, gender, fibrosis or AFP level of HCC in the cohorts. Therefore, the 2-TF signature was involved in the tumor grade, TNM stage and tumor invasion in HCC.

According to the guideline, TACE is a standard treatment for intermediate-stage HCC (Chang et al., 2020). The clinical treatment cohort of GSE104580 was applied to investigate the association of the 2-TF signature with TACE treatment benefit. The results demonstrated that patients with a low-risk score had a desirable efficacy for TACE treatment (Figure 6A). Using the pRRophetic algorithm, we further evaluated the IC₅₀ of numerous agents to determine latent therapeutic drugs. There were 67 and 55 agents associated with the risk score in the TCGA and ICGC cohorts respectively, 38 of which were common between them (Figure 6B). The IC₅₀ difference of the 38 drugs was displayed in Figure 6C (all adjusted p < 0.05). Our results indicated that the high-risk patients tended to be more sensitive to numerous targeted agents, such as ABT-888 (veliparib, a PARP inhibitor), AZD8055 and rapamycin (both mTOR inhibitors), as well as the eight common chemotherapeutic agents, including ATRA, camptothecin, cisplatin, cytarabine, epothilone-B, etoposide, gemcitabine and methotrexate. Whereas patients in the low-risk group were prone to lapatinib, erlotinib and other targeted agents (Figure 6C). Based on our results, combined treatments of TACE with these drugs might be promising therapeutic regimens if administrated according to the classified subtypes of HCC patients.

To uncover the biological characteristics influenced by the risk score, the DEGs between the high- and low-risk groups in the TCGA and ICGC cohorts were identified and underwent GO enrichment and KEGG pathway analyses. According to the most significant GO terms as shown in Figure 7, upregulated genes in the high-risk group both in the TCGA (Figure 7A) and ICGC (Figure 7B) cohorts were mainly enriched in cell proliferation-related functions, such as nuclear division and chromosome segregation. Similarly, the top KEGG pathways related to the upregulated genes in the high-risk group were cell cycle, cellular senescence and oocyte meiosis in both cohorts (Figures 7A, B). The genes upregulated in the low-risk group were enriched in various biosynthetic and metabolic processes both in the TCGA and ICGC cohorts, such as fatty acid metabolic process, alpha-amino acid metabolic process, carboxylic acid biosynthetic and catabolic processes, and organic acid biosynthetic and catabolic processes (Figures 7A, B). Consistently, according to the KEGG pathways, upregulated genes in the low-risk group were preferentially related to retinol metabolism, cytochrome P450-related metabolism, complement/coagulation cascades and peroxisome proliferator-activated receptor (PPAR) signaling (Figures 7A, B). Therefore, the 2-TF signature divided the patients into proliferative and metabolic subtypes with different biological processes, which may explain the prognostic prediction performance.

There were 191 common DEGs between the high- and low-risk groups in the TCGA and ICGC cohorts (Figure 8A). Subsequently, a PPI network containing 164 nodes and 1,358 interactions was established based on the 191 DEGs. Additionally, two significant modules with a score ≥5, namely module 1 comprising 43 upregulated genes in the high-risk group and module 2 comprising 17 upregulated genes in the low-risk group, were identified via MCODE (Figure 8B). The top 10 hub genes were further screened out, including BUB1B, TYMS, PLK1, RRM2, NUF2, CDC20, BIRC5, CDC45, CDT1 and CCNB1, which were all included in module 1 (Figure 8C). All the hub genes were validated to be closely correlated in HCC based on the TCGA, ICGC and GSE116174 cohorts (Supplementary Figure S3A–C). Enrichment analysis demonstrated that the most significant pathway in module 1 were cell cycle, and module 2 were mainly enriched in complement and coagulation cascades pathway (Figure 8D).