{'fraction_non_alphanumeric': 0.033077519587132655, 'fraction_numerical': 0.020643367115009224, 'mean_word_length': 5.4468195138769175, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Prognostic comparative genes predict targets for sorafenib combination therapies in hepatocellular carcinoma

Graphic abstract


a b s t r a c t
With the increasing incidence and mortality of human hepatocellular carcinoma (HCC) worldwide, revealing innovative targets to improve therapeutic strategies is crucial for prolonging the lives of patients. To identify innovative targets, we conducted a comprehensive comparative transcriptome analysis of 5,410 human HCCs and 974 mouse liver cancers to identify concordantly expressed genes associated with patient survival. Among the 664 identified prognostic comparative HCC (pcHCC) genes, upregulated pcHCC genes were associated with prognostic clinical features, including large tumor size, vascular invasion and late HCC stages. Interestingly, after validating HCC patient prognoses in multiple independent datasets, we matched the 664 aberrant pcHCC genes with the sorafenib-altered genes in TCGA_LIHC patients and found these 664 pcHCC genes were enriched in sorafenib-related functions, such as downregulated xenobiotic and lipid metabolism and upregulated cell proliferation. Therapeutic agents targeting aberrant pcHCC genes presented divergent molecular mechanisms, including suppression of sorafenib-unrelated oncogenic pathways, induction of sorafenib-unrelated ferroptosis, and modulation of sorafenib transportation and metabolism, to potentiate sorafenib therapeutic effects in HCC combination therapy. Moreover, the pcHCC genes NCAPG and CENPW, which have not been targeted in combination with sorafenib treatment, were knocked down and combined with sorafenib treatment, which reduced HCC cell viability based on disruption to the p38/STAT3 axis, thereby hypersensitizing HCC cells. Together, our results provide important resources and reveal that 664 pcHCC genes represent innovative targets suitable for developing therapeutic strategies in combination with sorafenib based on the divergent synergistic mechanisms for HCC tumor suppression.

Introduction
Hepatocellular carcinoma (HCC) is the main subtype of primary liver cancer and the sixth most common cancer and the fourth leading cause of death of cancer patients worldwide [1]. Most HCC patients are found in developing countries in Asia and Sub-Saharan Africa. However, the incidence of HCC in developed countries, such as the USA, Japan and some European countries, is rising because of increasing cases of hepatitis C virus infection and nonalcoholic steatohepatitis (NASH) associated with metabolic syndrome or diabetes mellitus [2]. Hence, the World Health Organization estimates that more than 1 million patients will die from liver cancers in 2030 [3]. Depending on the tumor burden and liver function of HCC patients, treatments with hepatic surgical resection, local ablation and liver transplantation are the major curative interventions for early HCC patients [4]. Nevertheless, with the vast majority of HCC patients being first diagnosed at the non-curable advanced stage, patients will receive systemic therapies with antiangiogenic multiple kinase inhibitors and/or metronomic chemotherapeutic drugs as well as immune checkpoint inhibitors [4][5][6][7]. Sorafenib is the standard FDA-approved first-line therapeutic drug for the treatment of advanced-stage HCC, and it has a response rate of<10% and a therapeutic benefit of prolonging life by a median of 3 months for HCC patients [8,9]. Although therapeutic response rates exceeded 20% in HCC patients after treatment with the newly developed antiangiogenic multiple kinase inhibitors (e.g., lenvatinib) and/or immune checkpoint inhibitors (e.g., nivolumab and atezolizumab) [10][11][12], applying innovative genomic approaches to identify additional biomarkers for guiding suitable therapies and targets for developing novel therapeutic strategies remains critical for prolonging survival for the majority of HCC patients [4,13].
After applying comprehensive sequencing approaches to the HCC genomes, only a limited number of somatic mutations were revealed as potential therapeutic targets; moreover, frequently observed mutations are not currently applied for predicting therapeutic responses in clinical practice [14,15]. The genetic heterogeneity of HCC genomes could be due to the accumulation of somatic mutations during long-term exposure to divergent etiological risk factors, such as viral hepatitis, alcohol abuse, aflatoxin ingestion and metabolic syndrome. In particular, intratumoral heterogeneity with frequent HCC genomic alterations has been shown to confer selective advantages of cell proliferation, and it is commonly associated with tumor recurrence, drug resistance and metastasis, thus leading to poor outcomes in HCC patients [16][17][18].
Comparative and integrated cancer genomic approaches are commonly applied to identify cross-species aberrant cancer genes and in best-fit mouse models of cancer subtypes of several cancers to reveal the tumorigenic mechanisms and develop therapeutic approaches to improve therapeutic interventions [19][20][21][22][23][24]. Notably, limited efforts were applied to reveal the overall consensus of pathophysiological functions shared in human and mouse cancer genomes, reveal the common tumorigenic pathways and predict innovative targets for dissecting the molecular mechanisms underlying combination therapies for therapeutic improvement [22,25].
In this study, we applied the most comprehensive comparative transcriptome analysis to date (to the best of our knowledge) that included large-scale liver cancer transcriptome datasets including 5410 human HCC samples and 974 mouse liver cancers to reveal the major concordant pathophysiological functions in HCC progression. As enriched HCC differentially expressed genes (DEGs) are commonly associated with drug metabolism and cell proliferation [26], we hypothesized that comparative HCC genes associated with patient prognosis might include a large number of known and unknown theranostic targets suitable for sorafenib combination therapies for developing innovative strategies to prolong the lives of HCC patients. Indeed, we identified 664 prognostic comparative HCC (pcHCC) genes, $44% of which have not been thoroughly studied in HCC, uncovered targets for sorafenib combination therapies with divergent molecular mechanisms to synergize HCC tumor suppression, and revealed an innovative oncogenic p38/STAT3 axis that enhances sorafenib tumor suppression for potential improvement of therapeutic interventions in HCC patients.

Identification of 664 pcHCC genes associated with patient prognosis
We collected large-scale transcriptome datasets of human HCCs and mouse liver cancers from public databases, including 5410 human HCCs (52 datasets from 26 array platforms, Supplementary Table 1) and 974 mouse liver tumors (52 datasets from 23 array  platforms, Supplementary Table 2), for comparative analysis. All human transcriptomes were derived from HCC tumors, with the exception of the cell line transcriptomes. The majority of mouse transcriptomes were annotated as HCC ($70%) or liver cancer ($30%) in genetically modified models with and without carcinogen induction (Supplementary Table 2 for details). After inferring the gene expression intensity based on the log2 tumor/normal ratio and Z score normalization, we generated expression profiles of human HCCs that consisted of 3,339 tumor/normal ratios of 26,182 differentially expressed genes (DEGs) and expression profiles of mouse liver cancer that consisted of 633 tumor/normal ratios of 18,157 DEGs (analysis pipeline in graphic abstract). We then determined the concordant aberrantly expressed comparative DEGs and revealed that 1,904 (766 upregulated and 1,138 downregulated) DEGs participated in both human HCC and mouse liver cancer. To emphasize the clinical value, we matched the 1,904 comparative HCC DEGs with prognostic HCC DEGs in the TCGA_-LIHC dataset and identified 664 pcHCC genes for experimental and translational evaluation ( Fig. 1A and Supplementary Table 3). Notably, we found that few HCC studies have focused on 293 out of 664 ($44%) pcHCC genes (labeled ''unknown" as HCC genes) based on PubMed searches.

Validation of the aberrant known and unknown pcHCC genes
To validate aberrantly expressed pcHCC genes in HCC, we selected two independent HCC cohorts from two different platforms, namely, iCOD (integrated clinical omics database with transcriptomes conducted using Affymetrix HG-U133 plus 2.0 arrays) [27] and TCGA_LIHC (transcriptomes conducted by nextgeneration sequencing technology, RNA-Seq) [14], to confirm the altered gene expression in HCCs. Consistently, our results showed concordant aberrant pcHCC gene expression between two selected transcriptomic datasets (pcHCC genes from the meta-analysis results of this study, iCOD and TCGA_LIHC) with a very significant positive correlation (Pearson's correlation coefficient r > 0.61, P value < 2.2 Â 10 À16 ) ( Fig. 1B and Supplementary Fig. 1A and B, and Supplementary Table 4).
To detail the involvement of the 664 pcHCC genes in HCC tumor progression, we searched PubMed and found that a total of 354 upregulated (53.3% of 664) pcHCC genes, including 219 (62%, orange dots) that were previously reported to be related to HCC and 135 (38%, red dots) that were not previously associated with HCC ( Fig. 1A and Supplementary Table 3). The majority of the top 10 known upregulated pcHCC genes (orange dots), such as TOP2A, PRC1, and AKR1B10, were studied thoroughly, with multiple reports focusing on their roles in HCC tumor progression and poor prognosis [28][29][30]. In contrast, the top 10 unknown pcHCC genes (red dots) were rarely studied or even had zero HCC results (e.g., RRAGD) in PubMed searches. For the 310 (46.7% of 664) downregulated genes, 152 (49%, sky blue dots) were known to be downregulated, such as TAT, CYP2C9 and SLCO1B, and they were associated with better prognosis of HCC patients in some reports [31][32][33]. On the other hand, 158 unknown and downregulated pcHCC genes (51%, blue dots) have been seldom studied or even yielded no search results with HCC, such as CFHR3 and TTC36. Moreover, other known HCC genes that were upregulated (e.g., GPC3, AURKA and others) and downregulated (e.g., CLEC1B, CRHBP and others) were also identified in the majority of HCC samples (as percentages on the x-axis), with differential expression intensities between tumor and normal HCC samples (as percentages on the yaxis) ( Fig. 1A and Supplementary Table 3). Interestingly, several pcHCC genes were found to have very little or no representation among HCC studies, including CENPW and CENPU (upregulated) as well as OIT3, IDO2, CNDP1 and ANGPTL6 (downregulated), according to PubMed searches. Together, our results provide evidence that the 664 pcHCC genes not only participate in HCC tumorigenesis but are also useful for the development of innovative theranostic biomarkers for therapeutic interventions of HCC.

Validation of pcHCC genes associated with clinical features and patient survival
To further confirm the aberrant expression of the 664 pcHCC genes involved in altered liver function during tumor progression, we examined their aberrant expression associated with the clinical and demographic characteristics of HCC patients from the iCOD dataset ( Supplementary Fig. 2). Consistent with previous reports, aberrant expression of pcHCC genes in HCCs was associated with known independent predictors of patient prognosis, including vein invasion, tumor size, diabetes and late tumor stage [34][35][36]. Then, we conducted unsupervised clustering to divide the genes into two groups, namely, upregulated and downregulated pcHCC genes. As expected, according to the Kaplan-Meier plot, the cluster of HCC patients with high combined expression of upregulated pcHCC genes in the TCGA-LIHC dataset exhibited worse survival than the cluster of HCC patients with low combined expression, with the comparison of Cox coefficients showing a significant P value (log rank test, P < 0.0001) ( Fig. 2A). In contrast, the cluster of HCC patients with high combined expression of downregulated pcHCC genes showed a better prognosis than the cluster of HCC patients with low combined expression in the TCGA-LIHC dataset, with the comparison of Cox coefficients showing a significant P value (log rank test, P = 0.0021) (Fig. 2B).

Gene set enrichment analysis (GSEA) of 664 pcHCC genes
To reveal the core pathophysiological functions of the 664 pcHCC genes in HCC progression, we compared the functional hallmarks of the GSEA signatures of these genes with two HCC DEGs, including the HCC DEGs from the transcriptomes of human and mouse microarrays and from the TCGA-LIHC datasets. Based on the significance of the functional gene signatures, we found that genes associated with cell proliferation-related hallmarks, including genes for E2F targets, MYC targets, mitotic spindle, and G2M checkpoint, were upregulated, and genes associated with metabolism-related hallmarks, including genes for the metabolism of xenobiotics, bile acids, coagulation, fatty acids, and adipogenesis, were downregulated ( Fig. 3A and Supplementary Table 5). To further validate whether aberrant hallmarks of cell proliferation and metabolism play critical roles in the prognosis of HCC patients, we collected cell proliferation-related pcHCC genes (n = 116, Supplementary Table 6) and metabolism-related pcHCC genes (n = 88, Supplementary Table 7) to understand their relationship with the prognosis of HCC patients. Indeed, our results demonstrated that upregulation of cell proliferation-related pcHCC genes and metabolism-related genes was significantly associated with poor (log rank test, P = 0.004) and better (log rank test, P < 0.001) HCC prognosis, respectively (Fig. 3B). Together, our results showed that the core pathological features in the upregulation of cell proliferation and downregulation of certain metabolic pathways, especially lipid and xenobiotic (drug) metabolism-related pathways, presented significant involvement in tumor malignancy in human HCC and mouse liver cancer models.

Predicting targets for sorafenib combination therapies in pcHCC genes
Based on the 664 pcHCC genes exhibiting aberrant proliferation as known targets of the HCC first-line therapy drug sorafenib [26] and xenobiotic pathways known to have roles in drug metabolism, we hypothesized that pcHCC genes might involve multiple pathways that are correlated with sorafenib-altered HCC expression profiles in patients who received sorafenib therapy and thus may be useful for identifying targets and therapeutic agents for sorafenib combination therapy. In the TCGA-LIHC dataset, 28 patients who received sorafenib treatment and had RNA-Seq transcriptome data were downloaded for the analysis. After matching the sorafenib-altered HCC genes from the TCGA_LIHC dataset with the 664 pcHCC genes, we divided these genes into four possible quadrants, namely, Q1: up/up (n = 205), Q2: down/up (n = 116), Q3: down/down (n = 194), and Q4: up/down (n = 149), which represented the alterations in gene expression relative to the pcHCC genes/sorafenib-altered HCC genes, and the values were calculated as log2 ratios (Fig. 4A and Supplementary Tables 8 and 9). For instance, the representative top 10 known and unknown pcHCC genes were redistributed into four quadrants, as shown with the red lines/rectangles, for scenarios with potential further upregula-tion (Q1: up/up) and downregulation (Q3: down/down) after sorafenib treatment (Fig. 5A). In contrast, the green lines/rectangles indicate scenarios with opposite changes in expression (Q2: down/up; Q4: up/down) of the HCC genes after sorafenib treatment.
To further validate whether these sorafenib-altered pcHCC genes could be useful as targets in combination therapies, we conducted a GSEA and found that these genes participated in several signaling pathways, including the PI3K/AKT/mTOR, xenobiotic metabolism, angiogenesis, hypoxia and ferroptosis pathways, which were previously reported as sorafenib-related pathways and potential targets for sorafenib combination therapies ( Fig. 4B and Supplementary Table 10) [56]. For instance, AKR1B10 (aldoketo reductase family 1 member B10) in Q1 is known to convert carbonyl compounds on aldehydes or ketones to alcohols for detoxification in the liver, and its upregulation is associated with the prognosis of HCC patients [28]. Because AKR1B10 is upregulated in HCC after sorafenib treatment, combination therapy with an AKR1B10 inhibitor via suppression of the mTOR signaling pathway enhanced the inhibitory effect of sorafenib on HCC xenografts in a nude mouse model ( Fig. 4B and Table 1) [43]. In addition, SLC27A5 (in Q3), PCK1 (in Q2), SLC22A1 (in Q2) and SLCO1B3 (OATP1B3 in Q3) are known sorafenib transporters and/or partici- pate in sorafenib metabolism in the liver [45,46,50,57], and altered or aberrant expression of these transporters might modulate the therapeutic efficacy of sorafenib and/or activate the TXNRD1 and NRF2/KEAP1 pathways to increase oxidative stress and tumorigenicity in HCC. Therefore, combination therapy with sorafenib and NRF2/TXNRD1 inhibitors could be a promising strategy to improve HCC therapy [45,46,50,57]. Other known sorafenibaltered pcHCC genes, such as SLC7A11 (in Q4), are key modulators of ferroptosis, which has led to the application of another ferroptosis inducer, artesunate, in sorafenib combination therapy for HCC [58]. AURKB (Aurora kinase B, in Q4) is a known HCC oncogene, and its inhibitor was combined with sorafenib to further induce HCC tumor shrinkage [54] (Fig. 4B). Collectively, the results showing the different mechanisms to support sorafenib combination therapies (Table 1) provide a better understanding and rationale for exploring innovative mechanisms of sorafenib-altered pcHCC genes to further develop strategies for the improvement of sorafenib combination therapy.
2.6. Combination treatment with sorafenib and shRNAs of NCAPG or CENPW synergized suppression of cell viability by diminishing p38/ STAT3 signaling To explore the potential of innovative targets among these pcHCC genes for sorafenib combination therapies in HCC, we inves-tigated the mechanistic roles and testing therapeutic efficacy of two pcHCC genes, NCAPG (non-SMC condensin I complex subunit G, Fig. 3C) and CENPW (centromere protein W, Fig. 3D), which are involved in HCC cell proliferation but have not been studied in sorafenib combination therapy. In addition to the association of these genes with poor prognosis in HCC patients, aberrant upregulation of NCAPG and CENPW was significantly associated with extrahepatic metastasis (BCLC stage C) (Supplementary Fig. 3A and C) and tumor size but not with venous invasion (Supplementary Fig. 3B and D) in HCC patients. With these supportive results, we continued to evaluate the suppressive effects of knocking down the expression of upregulated NCAPG and CENPW with specific short hairpin RNAs (shRNAs) and then investigated the impact on downstream cell proliferation signaling in HCC cells to validate the sorafenib combination therapy (Fig. 5).
First, we confirmed the upregulated expression of NCAPG and CENPW in locally recruited HCC tumors compared to their adjacent normal tissues (n = 21 pairs) through RT-qPCR experiments ( Fig. 5A and B) and knocked down NCAPG and CENPW with their corresponding shRNAs to evaluate their roles in the predicted functions of HCC cell proliferation and migration. As expected, knockdown of NCAPG (shNCAPG #2 and #3) and CENPW (shCENPW #3 and #4) with two independent shRNAs in two HCC cell lines, Mahlavu and HCC36, significantly diminished cell proliferation and decreased cell migration, as shown by transwell migration assays (Fig. 5C-F). Our results suggested that upregulated expression of NCAPG and CENPW plays important roles in HCC tumor growth and that cell migration leads to tumor progression.
Second, to better understand the molecular pathways involved in NCAPG-and CENPW-driven HCC tumor progression, we examined the protein expression in the MAPK and STAT3 pathways, which are known for their involvement in HCC cell proliferation, survival and migration, through Western blotting analysis [59,60]. Our results not only confirmed the knockdown efficiency of shRNAs against CENPW and NCAPG but also validated the reduction in cell proliferation by decreasing the expression of the cell proliferation markers PCNA and MCM2 (Fig. 5G and 5H). Interestingly, knockdown of NCAPG and CENPW with shRNAs consistently diminished the expression of p-STAT3 (pSer727) and MAPK p-p38 (pThr180/pTyr182) signaling but not ERK/p-ERK (pThr202/Tyr204) and JNK/p-JNK (pThr183/Tyr185) MAPK signaling, as shown by the Western blotting analysis (Fig. 5G and H). Previous studies indicated that the p38/STAT3 axis (but not the ERK/JNK MAPK signaling pathway) is involved in the expression of proinflammatory cytokines (e.g., IL6) and the stress-induced acute phase response of HCC cells [61][62][63]. Nevertheless, the biological functions of NCAPG and CENPW in the p38/STAT3 axis in HCC might be different because knockdown of NCAPG diminished the expression of p-p38, p-STAT3 and STAT3 but not the overall protein level of p38, and knockdown of CENPW reduced the expression of p-p38, p38, and p-STAT3 but not the total protein level of STAT3.
Finally, we conducted combination treatments of sorafenib with shRNAs of NCAPG or CENPW to identify the viability of HCC cells. Our results demonstrated that knockdown of either NCAPG or CENPW in Mahlavu and HCC36 cells conferred significant hypersensitivity to sorafenib treatment in a dose-dependent manner (Fig. 6A-D). Moreover, sorafenib treatment of shRNA-knocked down NCAPG or CENPW Mahlavu cells further decreased the expression of p-p38 and p-STAT3 in a dose-dependent manner, as shown by the Western blotting analysis (Fig. 6E and F). Together, our results suggest that NCAPG and CENPW, which suppress p38/ STAT3 signaling, are potential targets for sorafenib combination therapies in HCC. In addition, our results suggest that the 664 pcHCC genes might contain innovative targets and pathways for developing therapeutic agents for HCC sorafenib combination therapies.

Discussion
In this study, we conducted a comprehensive comparative genomic analysis of large-scale transcriptome datasets (5410 human HCC and 974 mouse liver cancer) to identify concordantly expressed HCC genes associated with altered clinical features in HCC progression and patient prognosis (pcHCC genes) and reveal innovative targets for sorafenib combination therapy. Although HCC transcriptome analyses have been conducted in recent years with a focus on altered cell proliferation and/or metabolic reprogramming with different approaches [18,[64][65][66], to our knowledge, our results are the first to indicate that the enriched signatures harbored innovative targets ($44%) suitable for developing therapeutic agents/inhibitors in sorafenib combination therapy. We provided lines of evidence from in silico analyses ( Fig. 4 and Table 1) and experimental studies of the selected pcHCC genes NCAPG and CENPW (Figs. 5 and 6) to support our hypothesis. Experimentally, we revealed that the proinflammatory and stress-related p38/STAT3 axis, but not conventional ERK/JNK MAPK signaling, plays critical roles in HCC progression and is suitable as a target for sorafenib combination therapy. Moreover, we identified 664 pcHCC genes harboring numerous innovative HCC targets   ($44%) for developing HCC biomarkers associated with patient prognosis, thus revealing innovative signaling pathways in tumor progression and identifying theranostic targets to predict drugs/ agents for sorafenib combination therapies to potentially prolong the life of HCC patients. With the increasing incidence and mortality of HCC worldwide, rapid developments have been made in systemic therapies for advanced HCC and several therapeutic options have been FDA-approved as first-line therapies, including sorafenib [8,9], lenvatinib [10], and atezolizumab plus bevacizumab [12]. Other antiangiogenic multikinase inhibitors, such as regorafenib and cabozantinib, as well as the antibody drug ramucirumab, which targets anti-vascular endothelial growth factor-2 (VEGF-R2), and immune checkpoint inhibitors (ICIs) against PD-1 and cytotoxic T lymphocyte antigen 4 (CTLA-4), were approved for the secondline treatment of HCC [11,67,68]. Despite these therapeutic advances, only approximately 25% of HCC patients show a therapeutic response and improved overall survival, and the current therapeutic regimens remain ineffective for prolonging the life of most HCC patients [69]. Although sorafenib therapy in HCC offers some benefits, including a 2.8-month longer median overall survival than the placebo, endeavors to improve HCC sorafenib therapy, including the discovery of therapeutic biomarkers and combination therapy with other drugs, remain critical due to the limited availability, affordability and effectiveness of sorafenib as a single-drug HCC therapy.
For the discovery of therapeutic biomarkers of sorafenib, since frequent substitutional driver mutations with therapeutic options have not been identified after large-scale sequencing efforts in human HCC genomes [14], copy number alterations, especially gene amplification, leading to increasing expression of target genes have attracted specific attention to select subgroups of HCC patients for targeted sorafenib therapy. Indeed, the FGF3/FGF4 amplicon at chromosomal region 11q13 [70] and the VEGFA amplicon at chromosomal region 6p21.1 [71] were amplified in $2% and 7 $ 10% of human HCC genomes, respectively, and they encode the potential oncogenic driver proteins FGF4 and VEGF-A, which showed better responses to sorafenib treatment. Although copy number alterations in the HCC cancer genome commonly affect a small population of cancer patients, the identification of therapeutic targets/biomarkers could stratify subgroups of HCC patients to develop more precise HCC therapies. Notably, we also explored the correlations of the 664 aberrantly expressed pcHCC genes with copy number alterations (CNAs) in the HCCs of the TCGA_LIHC dataset ( Supplementary Fig. 4). The lack of statistical associations between CNAs and aberrant gene expression of the 664 pcHCC genes (r = 0.0058, P value = 0.1373) implied the low potential for developing pcHCC genes as CNA-based therapeutic biomarkers in HCC.
Although some comparative genomics approaches have been applied for decades to develop best-fit mouse models of cancer subtypes and identify driver genes in cancers, including HCC [21,23,72,73], a systematic and detailed explorations of comparative HCC genes implicated in the core HCC pathophysiological and clinical features and patient prognosis and even the targets of sorafenib combination therapy are lacking. Based on the aberrant enrichment of upregulated cell proliferation and downregulated lipid and xenobiotic metabolism-related genes of the 664 pcHCC genes, we hypothesized and systematically explored the underlying molecular mechanisms of sorafenib therapies in HCC. In addition to matching the divergently acting molecular mechanisms of sorafenib combination therapies reported previously (Table 1), we also selected two pcHCC genes, NCAPG and CENPW, because the aberrant upregulation of these genes is associated with the clinical features and prognosis of HCC patients (Supplementary Fig. 3A-D); moreover, these genes have not been previ- ously targeted in sorafenib combination therapies to explore the oncogenic mechanisms and develop therapeutic strategies. NCAPG encodes a non-SMC (structural maintenance of chromosome) condensin I complex subunit G protein that has been shown to be upregulated in HCC to activate the PI3K/AKT pathway and enhance cell proliferation, and it has been associated with poor prognosis in HCC patients [74][75][76][77][78][79]. Upregulated CENPW, the centromere protein W, which is sometimes also called cancer upregulated gene 2 (CUG2), was found to enhance cancer epithelial to mesenchymal transition (EMT) and stemness via activation of the canonical TGFb signaling pathway and tumor progression in other cancers [80,81]; however, its role in HCC is unknown. Neither the upregulation of NCAPG nor CENPW or related inhibitors has been shown to play a role in sorafenib combination therapy in HCC. Interestingly, we found that a rarely studied p38/STAT3 proinflammatory/stress-related axis (not JNK/ERK signaling) was required to sustain HCC oncogenic features. The combination of shRNAs of NCAPG or CENPW with sorafenib in the treatment of HCC cells showed synergistic suppression of oncogenic features and downstream signaling, further suggesting that pcHCC genes are an important resource as theranostic targets for early biomarker detection and therapeutic agents/inhibitors for sorafenib combination therapy to improve HCC therapy.
Because of the increasing need for cancer therapy worldwide, including HCC therapy, and the superior therapeutic benefits of combination therapy [12,[82][83][84], systematic exploration of sorafenib action mechanisms and innovative targets for combination therapy is crucial for improving HCC therapy. Our systematic prediction of potential targets among 664 pcHCC genes for the sorafenib combination therapy provides rationales and resources for further improving HCC therapy with the following features. First, we provide a comprehensive landscape of predictive sorafenib combination therapies with divergent mechanisms based on different scenarios of altered HCC gene expression ( Fig. 4 and Table 1), and the findings were validated using selected target genes (Fig. 6). Second, the innovative pcHCC targets identified in our study could also be biomarkers for predicting sorafenib therapeutic responses for tumor recurrence. Finally, a similar approach could be applied to other cancer types for the development of therapeutic agents/inhibitors to improve the combination therapy of other antiangiogenic multikinase inhibitors, and such approaches could be validated in comparative mouse cancer models. Our approach of applying concordantly expressed comparative cancer genes to forecast targets and therapeutic agents in combination cancer therapy could be difficult to implement due to the limited availability of transcriptomic datasets with clinical features derived from comparative cancer models and cancer patients treated with and without therapeutic agents. Nevertheless, with increasing cancer incidence and mortality predicted worldwide in the coming decade, the successes of combination cancer therapies reported in multiple cancers will encourage more clinical trials on combination therapies with the goal of reducing the toxicity and cost to prolong the survival of cancer patients.

Conclusions
To improve HCC therapies and meet emerging needs worldwide, we conducted a large-scale comparative transcriptome analysis and identified 664 pcHCC genes ($44% are innovative targets) associated with prognostic clinical features and enriched driver gene signatures suitable for developing sorafenib combination therapies. We provided lines of evidence to demonstrate the divergent mechanisms of sorafenib combination therapies in HCC and experimentally validated selected pcHCC genes by reducing their expression in HCC cells via a diminished innovative p38/STAT3 axis to hypersensitize HCC cells in sorafenib-combined treatments. We provided important insights on 664 pcHCC genes, and they include innovative targets and molecular signaling factors suitable for systematic targeting in sorafenib HCC combination therapy, which has the potential to prolong the life of HCC patients.

Data sources and preprocessing
The transcriptome datasets were downloaded from the NCBI Gene Expression Omnibus [85]. We selected transcriptomes derived from tumors of human hepatocellular carcinoma (HCC)only and mouse liver cancer models but excluded transcriptomes from cell lines ( Supplementary Tables 1 and 2). The majority of mouse transcriptomes were annotated as HCC ($70%) or liver cancer ($30%) in genetically modified models with and without carcinogen induction, as detailed in Supplementary Table 2. The ratio of tumor-versus-normal samples for each probe was calculated for each dataset. If there was no normal sample in a dataset, a pool of normal samples was constructed from normal samples from other datasets of the same platform. Two independent HCC datasets from iCOD and TCGA_LIHC (the level 3 RNA-Seq data of HCC tumors) were downloaded based on the original report [27] and The Cancer Genome Atlas [86] using the R package TCGAbiolinks (https://bioconductor.org/packages/release/bioc/html/ TCGAbiolinks.html) [87].

Integrated meta-analysis of HCC and mouse genes
All probes were mapped to the gene ID in the Ensembl database (https://www.ensembl.org, build 92), and mouse Ensembl gene IDs were then mapped to human Ensembl gene IDs [88]. An average ratio was calculated for the multiple probes mapped into one Ensembl gene ID. For all tumor/normal ratios, log2 transformations were performed, followed by the Z-score transformation of each sample to standardize data across samples in all datasets [89]. Student's t-test was performed for each Ensembl gene ID.

Determination of comparative differentially expressed genes (DEGs) and gene set enrichment analysis (GSEA) in HCC
To determine the comparative DEGs, we defined the upregulated probes when the processed log2 tumor/normal ratios were more than 1 and the downregulated probes when the ratios were less than À1. The percentage of DEGs was calculated from the sum of upregulated or downregulated labels versus the total sample numbers. An adjusted percentage for defining the upregulated and downregulated genes was calculated using a formula (percentage of upregulated labels minus percentage of downregulated labels) for each Ensembl gene ID. For ensembled gene IDs with <25% presented in the samples, those with a p value of more than 0.05 in Student's t-test, <5% total altered sample percentage, and 5% to À5% adjusted regulated percentage were filtered out. Only the Ensembl gene IDs having the same regulatory trend in both human and mouse samples were selected as comparative DEGs. Gene set enrichment analysis (GSEA) was applied to detect enriched pathways of the DEGs with FDR < 0.05 and P value <0.05 for consideration of significant enrichment [90].

Archived HCC samples
Archived HCC tissues and the corresponding normal tissues (n = 21) were collected from National Taiwan University Hospital with IRB approval [17]. Tissue was harvested after curative HCC surgery, frozen in liquid nitrogen, and stored in a À80°C freezer immediately after being sliced into small pieces.

RNA isolation, RT-PCR and quantification of RNA expression
Total RNA was extracted from HCC tissues using TRIzol (Invitrogen) according to the manufacturer's instructions. RNA purity was assessed with a NanoDrop spectrophotometer (Thermo Scientific). RNA was reverse-transcribed into cDNA with oligo-dT primers using a Maxima H Minus First Strand cDNA Synthesis Kit (Thermo Scientific). RT-PCR was performed with Taq DNA Polymerase Master Mix RED (Ampliqon). PCR cycling conditions were 95°C for 30 s, 60°C for 30 s and 72°C for 20 s for 40 cycles. The PCR products were separated through gel electrophoresis with DNA View (Tools Biotech) in a 1.5% agarose gel and identified through UV exposure. RNA expression quantification was performed with GelQuant. NET software provided by biochemlabsolutions.com. Target gene expression of certain HCC samples was normalized to the internal control GAPDH expression. The list of primers used for RT-PCR is included in Supplementary Table 11.

Transwell migration assay
Cell migration assays were conducted as described previously [84]. Each condition was performed in triplicate wells with seeding of a total of 3 Â 10 4 cells/well for Mahlavu or 8 Â 10 4 cells/well for HCC36 HCC cell lines. After 24 h, the migrated cells were stained with Giemsa solution (Merck) and counted in 100 Â microscopic fields with ImageJ (NIH, USA).

Cell viability assays
Cell viability was assessed in triplicate wells for each condition using the alamarBlue assay (BioSource International) [84]. A total of 1 Â 10 4 Mahlavu or HCC36 HCC cells/well were seeded in 6well plates and incubated in a 37°C incubator. Cell viability at days 0, 1, 3, 5, 7 and 9 post-seeding was analyzed by detecting the absorbance using a microplate reader. Sorafenib dissolved in dimethyl sulfoxide (DMSO) at various concentrations was added to the medium of HCC cell cultures.

Statistical analysis and data visualization
All in vitro conditions were performed in triplicate, and values are shown as the mean ± 95% confidence interval (CI). The differences between groups were determined using Student's t-test, and statistical significance was considered when P < 0.05. All analyses were performed in the R statistical programming environment. All plots were produced with the ggplot2 package (https:// ggplot2.tidyverse.org/). All scripts for data analysis and visualization are available upon request.

Declaration of Competing Interests
The authors declare no conflicts of interest.