Multi-omics data integration reveals novel drug targets in hepatocellular carcinoma

Background Genetic aberrations in hepatocellular carcinoma (HCC) are well known, but the functional consequences of such aberrations remain poorly understood. Results Here, we explored the effect of defined genetic changes on the transcriptome, proteome and phosphoproteome in twelve tumors from an mTOR-driven hepatocellular carcinoma mouse model. Using Network-based Integration of multi-omiCS data (NetICS), we detected 74 ‘mediators’ that relay via molecular interactions the effects of genetic and miRNA expression changes. The detected mediators account for the effects of oncogenic mTOR signaling on the transcriptome, proteome and phosphoproteome. We confirmed the dysregulation of the mediators YAP1, GRB2, SIRT1, HDAC4 and LIS1 in human HCC. Conclusions This study suggests that targeting pathways such as YAP1 or GRB2 signaling and pathways regulating global histone acetylation could be beneficial in treating HCC with hyperactive mTOR signaling. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07876-9.


Background
Liver cancer is the second leading cause of cancerrelated deaths worldwide, and hepatocellular carcinoma (HCC) accounts for approximately 90% of primary liver cancer cases [1]. Approximately 50% of HCC tumors exhibit loss of the tumor suppressors Pten, Tsc1, or Tsc2 leading to aberrant PI3K-AKT-mTOR signaling. However, the effector pathways via which mTOR promotes tumorgenicity are widely unknown. We generated an mTOR-driven HCC mouse model, by liver-specific deletion of the tumor suppressors Pten and Tsc1 [2,3], to investigate the molecular and cellular mechanisms of mTOR-driven tumorgenicity.
While DNA sequencing has enabled a comprehensive characterization of tumor genomes and stratification of patients, translating such information into treatment strategies has remained a major challenge. A limitation of relying entirely on genomic data to determine a therapeutic strategy is that it ignores functionally-related, non-mutated genes that could also encode potential drug targets. In addition, different mutations across cancer patients (genetic divergence) could result in the same pathways being activated (functional convergence) [4]. Apart from somatic mutations, tumorigenesis can be regulated by the levels of specific miRNAs, mRNAs, proteins and protein phosphorylation. miRNAs are key regulators of the transcriptome and can act as either oncogenes or tumor suppressors. Common mechanisms that can dysregulate miRNA expression in human cancers include amplification, deletion or epigenetic changes [5]. Transcriptomic and proteomic analyses have been performed to stratify HCC patients into clinically-relevant groups [6][7][8]. However, to further understand the effect of a genetic aberration or dysregulated gene expression (possibly due to aberrant miRNA expression) it is necessary to identify the mediators common to diverse alterations. Distinct genomic aberrations (in different tumors) are expected to converge functionally on the same downstream protein, referred to here as a 'mediator'. To identify such mediators, it is essential to integrate omics data, i.e., the genome, transcriptome, proteome and phosphoproteome (commonly referred to as multi-omics analysis), from diverse tumors.
Recently, multi-omics analysis has been informative in the characterization of tumors. For example, integration of DNA, RNA and phosphoproteomic data enabled stratification of prostate cancer patients and to identify individualized treatment options [9]. Computational methods that focus on the direct effect of genetic aberrations, i.e., the effect of a gene mutation on the encoded protein, have also been proposed [10]. However, a major drawback of these studies is that they rely solely on genomic analysis. New methods are necessary to integrate different types of omics data to identify dysregulated pathways.
In this study we use NetICS, a computational method to integrate multi-omic data (somatic mutations, miRNA differential expression, transcriptomics, proteomics and phospho proteomics) from an mTOR-driven mouse HCC tumor model [11], to understand the molecular mechanisms of mTOR-driven HCC. NetICS provides a comprehensive framework that reveals how specific genetic aberrations (i.e., deletion of the tumor suppressors Pten and Tsc1) and tumor-specific changes in miRNA expression can affect downstream mediators. NetICS employs a sample-specific network diffusion process that reveals the convergence of diverse changes in distinct tumors (mutations and differentially expressed miRNAs) on common downstream mediators. The identified mediators include transcription factors, kinases, phosphatases and deacetylases. While, some of the mediators are known oncogenes, others are novel oncogenic mediators. These mediators are potential, novel drug targets.

Results
We isolated tumors from an HCC mouse model generated by liver-specific deletion of Pten and Tsc1 (Fig. 1A). We hereafter refer to this model as the liver-specific doubleknockout (L-dKO) mouse. We isolated twelve distinct liver tumors, three each from four 20 week old L-dKO mice. To detect somatic mutations, we compared exome sequence data from tumors and from matched muscle tissue (Fig. 1B). For all other analyses (RNA, miRNA, proteome and phosphoproteome), we compared the tumor nodules to healthy liver tissue from six control mice (cre-negative, age-and sex-matched littermates) (Fig. 1C). We detected a total of 157 point mutations and small insertions/deletions in the twelve tumors (Table 1). Except for the originally introduced Pten and Tsc1 deletions (Fig. S1), no specific mutation was found in more than a single tumor ( Fig. 2A-B). In contrast to somatic mutations, we detected mRNAs, miRNAs, proteins and phosphosites commonly dysregulated across multiple tumor samples ( Fig. 2C-F). On average, 4,348 mRNA, 108 miRNAs, 2,389 proteins and 906 phosphosites were dysregulated in the tumor samples ( Fig. 3A-D).
To identify the downstream mediators, we used NetICS, a network-based method that integrates multi-  Fig. 1 Experimental setting. A. Representative images of whole livers from 20-week-old L-dKO (tumors are indicated with arrowheads) and control mice. B. Three independent tumor samples were taken from each of four 20-week-old L-dKO mice. For each mouse, one muscle sample was used as the matched healthy tissue. This setting was used for exome sequencing. Each of the three tumor nodules was compared against the matched muscle tissue sample. C. Three independent tumor samples were taken from each of four 20-weekold L-dKO mice. Liver tissues from six cre-negative age-and sexmatched littermates were used as a control. This setting was used for mRNA sequencing, miRNA sequencing, proteome and phosphoproteome quantification omic data to prioritize cancer genes [11]. NetICS provides a framework to simulate how upstream events lead to the dysregulation of downstream genes and proteins. It detects how mediators are dysregulated in each sample, using sample-specific network diffusion. NetICS then systematically integrates the individual ranks to infer a global gene ranking across all tumor samples [11]. However, our NetICS framework failed to predict functional convergence among the 157 detected somatic mutations (except for Pten and Tsc1) after a random permutation test, suggesting that additional information is required to identify a common downstream mediator. Shown are the results of exome sequencing in the 12 tumor nodules. Relevant information are given such as the position in the chromosome, the reference and alternative alleles, the type of mutation, the amino acid substitution and variant allele frequency.
Hence, we integrated the Pten and Tsc1 deletions, the somatic mutations, and the differentially expressed miR-NAs in each sample as upstream events. As downstream events, we used the differentially expressed mRNAs, proteins and differentially regulated phosphosites per tumor. After systematic integration, NetICS analysis predicted 74 mediators that are functionally related to differentially expressed miRNA and somatic mutations (upstream) as well as to differentially expressed genes, proteins and phosphosites (downstream) ( Table 2). Pathway enrichment analysis of the mediators indicated a strong enrichment of cellular signaling pathways regulating cell cycle proteins and of TGF signaling, suggesting strong proliferation potential of malignant hepatocytes (Table 3). We also observed an upregulation of epithelial to mesenchymal transition (EMT) factors, suggesting increased metastatic potential. For example, many of the detected mediators are involved in Notch signaling, consistent with the observation that approximately 30% of human HCC displays active Notch signaling [12]. Furthermore, we observed upregulation of IL6 and leptin signaling which have been suggested to play crucial roles in the initiation and development of HCC [13,14]. Leptin is an important activator of cell proliferation and an inhibitor of cell death. Leptin signaling is also known to have angiogenic effects in multiple cancers including HCC [14]. The 74 mediators, consisting of various different types of regulatory proteins, include 14 kinases, 23 transcription factors, 2 deacetylates and 3 phosphatases. As expected, the mediators include known downstream targets of PTEN and TSC1, such as mTOR, AKT2 and AKT3. The mediators also include known HCC-related proteins, such as YWHAZ [15] and KLF4 [16]. Below, we discuss in detail five of the 74 detected mediators, namely YAP1, GRB2, HDAC4, SIRT1 and LIS1.
YAP1 is a transcription factor that activates genes involved in cell proliferation and suppresses apoptotic genes [17]. Directly upstream of YAP1 is the microRNA miR-375, the expression of which was significantly downregulated in L-dKO tumors compared to control liver tissues ( Fig. 4A and S2A). miR-375 has been shown to inhibit the expression of YAP1 [18]. Our data suggests that the reduced expression of miR-375 could in turn increase YAP1 protein levels that promote tumorigenesis. We also investigated multiple proteins whose expression is regulated by YAP1 (directly or indirectly). As expected, the transcript levels of YAP1 targets were significantly high in HCC tumors, including mRNAs of the direct YAP1 targets Ctgf, Birc5, Cyce1, Cyr61, Ki63 (Fig. S2B). Increased YAP1 levels upon immunoblot analysis in murine (4 out of 4 tumors) and human HCC (5 out of 5 patients) confirmed the dysregulation of YAP1 signaling in HCC (Fig. 5A-D).
The signaling adaptor protein GRB2 was found to be upregulated at both the mRNA and protein level in L-dKO tumors. According to the network (Fig. 4B), GRB2   Table 2 Genes are ranked based on the mediator score from the final ranked gene list across all tumor nodules. The top 5% ranked genes are shown, excluding the ones with a FDR adjusted P-value > 0.05 after the random permutation test. For each gene, the gene name and type are given. For each tumor nodule, +1 denotes upregulation at the RNA, proteome or phosphoproteome levels. Similarly, -1 denotes downregulation. If empty, the gene's, protein's or phosphosite's levels did not change significantly between tumor nodules and control samples Gene type RNA(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) PROTEOME(+1 significant upregulation, -1 significant downregulation (compared to normal samples))  Gene type RNA(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) PROTEOME(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) Mouse1-  Gene type RNA(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) PROTEOME(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) Mouse1- Gene type PROTEOME(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) PHOSPHOPROTEOME(number of phosphosites upregulated/number of phosphosites downregulated) Genes are ranked based on the mediator score from the final ranked gene list across all tumor nodules. The top 5% ranked genes are shown, excluding the ones with a FDR adjusted P-value > 0.05 after the random permutation test. For each gene, the gene name and type are given. For each tumor nodule, +1 denotes upregulation at the RNA, proteome or phosphoproteome levels. Similarly, -1 denotes downregulation. If empty, the gene's, protein's or phosphosite's levels did not change significantly between tumor nodules and control samples. Table 2 Genes are ranked based on the mediator score from the final ranked gene list across all tumor nodules. The top 5% ranked genes are shown, excluding the ones with a FDR adjusted P-value > 0.05 after the random permutation test. For each gene, the gene name and type are given. For each tumor nodule, +1 denotes upregulation at the RNA, proteome or phosphoproteome levels. Similarly, -1 denotes downregulation. If empty, the gene's, protein's or phosphosite's levels did not change significantly between tumor nodules and control samples

(Continued)
Gene type PROTEOME(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) PHOSPHOPROTEOME(number of phosphosites upregulated/number of phosphosites downregulated) Mouse3-  Gene type PROTEOME(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) PHOSPHOPROTEOME(number of phosphosites upregulated/number of phosphosites downregulated) Mouse3-  Gene type PROTEOME(+1 significant upregulation, -1 significant downregulation (compared to normal samples)) PHOSPHOPROTEOME(number of phosphosites upregulated/number of phosphosites downregulated) Mouse3-      [19], loss of PTEN in the L-dKO tumors correlates with upregulation of PTK2-GRB signaling. GRB2 signaling activates several proteins including SHC1, K-RAS and H-RAS (Fig. 4B). H-RAS is a small GTPase that positively controls phosphorylation of the transcription factor JUN. Mass spectrometry analysis showed that phosphorylation of Ser63 and Ser73 (indicating active JUN) in JUN was significantly increased in L-dKO tumors. JUN in turn regulates transcription of the gene CDK4. CDK4 transcript levels and protein levels were upregulated in L-dKO tumors. CDK4 is a known oncoprotein that can be targeted by inhibitors [20]. SHC1 activates MAPK1 and MAPK3, two known protein-serine/threonine kinases that participate in the RAS-RAF-MEK-MAPK signal transduction cascade and are known to be involved in tumorigenesis [21]. NetICS also detected two deacetylases, namely HDAC4 (class II histone deacetylase) and SIRT1 (class III histone deacetylase) as mediators. HDAC4 is known to mediate tumorigenesis through chromatin structure remodeling and controlling protein access to DNA in colon cancer [22], glioblastoma [23], ovarian cancer [24], gastric cancer [18], and esophageal carcinoma [25]. Immunoblot analysis confirmed that HDAC4 protein levels were significantly increased in murine and human HCC tumor tissues (Fig. 5A-D). These observations also suggest that mechanisms similar to mouse L-dKO tumors (i.e., miRNAs) could be regulating HDAC4 protein levels in human HCC. NetICS analysis suggested that five significantly downregulated miRNAs (miR-10a, miR-140, miR-22, miR-29b and miR-29c) could lead to increased HDAC4 levels in tumors ( Fig. 4C and S3). SIRT1 is another histone deacetylase detected as a mediator. Immunoblot analysis revealed that SIRT1 protein levels were significantly reduced in murine and human (four of five patients) HCC (Fig. 5A-D). The full length images are shown in supplementary figure 5 ( Figure S5). However, unlike HDAC4, SIRT1 mRNA levels were reduced in four out of twelve L-dKO tumors. Upstream of SIRT1, NetICS detected four miRNAs that were significantly upregulated, namely miR-138, miR-146b, miR-34a and miR-9, that could contribute to reduced SIRT1 levels ( Fig. 4D and S4). The role of SIRT1 in tumorigenesis is debated due to conflicting reports on SIRT1 as a tumor promoter or suppressor. SIRT1 deacetylates and downregulates two well-known tumor suppressors, TP53 and General 74 Shown are pathway enrichment results generated by using Metacore version 6.35.69300 at 21.07.2018. The final list of predicted genes from NetICS were used (Table 2). Pathways are ranked based on the adjusted FDR P-value (column E) and the number and names of the pathway genes that are present in the final list of 74 predicted genes is given (column F-G).
E2F1, suggesting an oncogenic role [26]. Conversely, SIRT1 also deacetylates and represses the oncogenic transcription factor β-catenin, suggesting a role as a tumor suppressor [27]. Based on our analysis, we suggest that SIRT1 has a tumor suppressing role in mTORdriven HCC tumors. NetICS also detected proteins largely unexplored in cancer biology as mediators of tumorigenesis. For example, LIS1 (lissencephaly-1) is a conserved regulator of dynein. It binds to dynein's motor domain and induces a tight microtubule-dynein interaction [28]. A potential role of LIS1 in tumor progression is now being explored [29,30]. We examined TCGA transcriptome data for LIS1 expression. We found that 47.1% of HCC patients have reduced LIS1 expression, suggesting that LIS1 has a tumor suppressing role in mTOR-driven tumors (Fig. 5E).

Discussion
We have utilized NetICS, a multi-omics data integration method that predicts mediators, and an mTOR-driven HCC mouse model to detect novel drug targets in HCC. NetICS detected 74 mediators that were ranked in the top 5% among network proteins. These mediators were found to be significant after a random permutation test of the aberrant and differentially expressed genes and proteins. We described five of the mediators in detail, namely YAP1, GRB2, HDAC4, SIRT1, and LIS1, and suggest upstream causes of their dysregulation as well as their downstream effects.
Importantly, NetICS is able to predict 'silent' genes as mediators, i.e., genes not affected by mutation or differentially expressed ( Table 2). This could be because NetICS scans the neighborhood of the potential    mediator and detects aberrant expression and mutation patterns even if the gene itself is neither mutated nor aberrantly expressed. To demonstrate the power of NetICS approach, we tested the ability of multi-omic NetICS to detect mediators that would not be predicted in singleomic approaches, i.e. RNA, proteome or phosphoproteome data alone. Of the 74 top 5%-ranked mediators detected in the multi-omics NetICS approach, we detected 12 relying exclusively on the transcriptome, and none relying only on the proteome or phosphoproteome ( Table 4). Thus, NetICS has power in predicting silent genes that would not be detected by a single-omics approach.
Pathway enrichment on the detected mediator genes suggested multiple tumor-related pathways that could be potentially targeted to curb tumor growth. We focused on the mechanistic insights and pathways of 5 of these mediators that we picked manually. NetICS suggests that overexpression of HDAC4 -which is frequently dysregulated in human malignanciesdrives tumor growth in HCC. Inhibitors of HDAC4, such as LMK-235 [31], could be potentially useful in HCC with HDAC4 overexpression. Similarly, HCC with YAP1 overexpression could benefit from using inhibitors for YAP1 [32].

Conclusions
To conclude, application of NetICS to multi-omics data from an mTOR-driven HCC mouse model detected new potential drug targets. This approach could be used to identify drug targets in other tumor types.

Animal experiments
Liver-specific Tsc1 and Pten double knockout mice were generated as described in [2] and [3] at Biozentrum, University of Basel. In short, tumors from 20 week-old L-dKO mice and whole liver from control mice were snap-frozen and pulverized. This powder was used for subsequent exome sequencing, total RNA sequencing (including miRNA and mRNA), proteomics and phosphoproteomics. For exome sequencing, muscle tissue from the quadriceps of 4 L-dKO mice was used as a control. The mice were on mixed genetic background (C57BL/6J, 129/SvJae, BALB/cJ). Age and sex matched littermate mice without the Cre gene were used as controls. Only male mice were used in all experiments. Mice were fasted overnight before euthanasia by CO2 inhalation. The total number of mice used were 6 control mice and 4 L-dKO mice.

Exome sequencing
DNA extracted from three tumor nodules and a muscle tissue sample each from four mice were subjected to whole-exome capture using the SureSelect Mouse All Exon (Agilent) capture system and to massively parallel sequencing on an Illumina HiSeq 2000 at the Genomics Facility Basel, ETH Zurich, Switzerland. A median of 141 and 97 million 101-bp paired-end reads were generated from DNA extracted from tumor nodules and the muscle, respectively, equivalent to median depths of 78x (tumor nodules, range 34x-124x) and 57x (germline, range 35x-129x; Table 5). Exome sequencing data have been deposited in the Sequence Read Archive under the accession SRP156216.
Somatic single-nucleotide variants (SNVs) were identified using MuTect (v1.1.4) [35] and somatic small insertions and deletions (indels) were identified using Strelka (v1.0.15) [36]. To remove false mutation calls resulting from sequencing and/or alignment artifacts, a panel of normal was created from the four normal samples in this cohort using the artifact detection mode of MuTect2 (packaged in GATK, v3.6). Variants present in at least two of the four samples in the panel of normal were disregarded. Variants outside the target regions, covered by <10 reads in the tumor or <5 reads in the germline were disregarded. Variants supported by <3 reads in the tumor or for which the tumor variant allele fraction was <5 times than that of the normal variant allele fraction were disregarded [37]. 157 putative somatic mutations passed the filters (Table 1).
FACETS [38] was used to define copy number alterations. Specifically, read counts for positions within the target regions with dbSNP (Build 142) entries were generated for each matched tumor nodule and normal samples as input to FACETS, which performs a joint segmentation of the total and allelic copy ratio and infers allele-specific copy number states. To enable detection of the intragenic deletions of Tsc1 and Pten, 15-20 evenly-spaced positions per deleted exon were tiled within the regions of the deletions (Fig. S1).

Transcriptome sequencing and quantitative PCR (qPCR) analysis
Raw fastq files were aligned to the reference genome Mus_musculus.GRCm38.72 using PALMAPPER with default parameters [39]. The length of the seeds of the PALMAPPER index was set to 15. Then we computed read counts using htseq-count against the reference genome annotation (Mus_musculus.GRCm38.72.gtf). Based on these counts we performed the differential gene expression analysis using DESeq2 where we compared each tumor sample individually against all six control samples. Exact numbers of detected dysregulated mRNA per tumor sample are given in Table  6. For quantitative PCR analysis, RNA was prepared as shown above, 500ng RNA was used to make cDNA using Superscript III (Invitrogen) as per manufacturers instructions. ABI Step One (Applied Biosystems) machine was used together with Syber Green PCR Kit (Invitrogen) and the primers below (100pM) to perform qPCR as per manufacturer instructions. TBP was used as a normalizer and ddCT method was used for analysis. miRNA sequencing libraries were generated using a modified protocol from [40]. Briefly, RNA from tissues was isolated using the Qiagen miRNAeasy kit as described above (section Animal experiments). 10 microgram of total RNA was run in a 15% polyacrylamide gel, the part containing small RNAs was cut and subjected to nucleotide extraction using overnight 0.4M NaCl and ethanol precipitation. Isolated small RNA mix was subjected to Illumina TrueSeq Small library preparation kit used as per manufacturer's instructions. Afterwards the small-RNA libraries were run in a 10% Polyacrylamide gel to clean up the adaptor-adaptor fraction. The gel part containing the small-RNA libraries was cut and libraries were extracted using overnight 0.4M NaCl and ethanol precipitation. They were run using Illumina NextSeq500 sequencer as per manufacturer's instructions. Exact numbers of detected dysregulated miRNA per tumor sample are given in Table 6.

Mass spectrometry for proteomics and phosphoproteomics
Liver tissues from L-dKO tumors and control mice were obtained as detailed above (Animal experiments). Label free mass spectrometry was performed on the tumor nodules. Tumor proteome was always compared to the proteome obtained from livers of six control mice pooled together. A detailed description about the Statistics about whole-exome sequencing are given. These include the total number of reads, the mean target coverage and the percent of target bases covered at least at 10X, 20X, 50X and 100X for each tumor and muscle tissue sample.
proteomics method used to analyse L-dKO tumor nodules and the softwares used for data analysis can be found in [2]. For the phosphoproteome, the desalted peptides were enriched for phosphopeptides using TiO2 beads. Detailed protocol is available in [3]. After data processing, the protein groups datasets and phospho peptide datasets were exported into a FileMaker Pro-12 databank. For statistical analysis, an R-based program -Perseus, version 1.4.0.2, was used [41]. ANOVA-based two-sample t-test was performed by adjusting S0 to 1 and the number of randomizations to 250 (default). The 5% FDR was used for analysis. Exact numbers of detected dysregulated proteins and phosphosites per tumor sample are given in Table 6.

Detection of differentially expressed mRNA and miRNA
The DESeq2 tool [42] with default settings was used to detect differentially expressed genes and miRNA between tumor and normal tissue. Every tumor sample was compared against the six control samples from healthy liver tissue. We considered as significant the genes detected with an FDR adjusted P-value lower than 0.05.

Interaction network
In order to construct a directed functional network, we downloaded functional interactions for the species Mus Musculus from three different databases including Kegg, Signor and miRTarBase. From miRTarBase, we only kept the interactions supported by strong experimental evidence (either reporter assay or western blot). The interactions cover a variety of types at different cellular levels, including (de) phosphorylation (phosphoproteome), expression/repression (RNA) and activation/ inhibition (proteome). Interactions characterized as "binding" or "complex" were treated as undirected edges. The network contained 5,546 genes and 44,423 interactions in total. In order for network diffusion to converge to a unique solution (steady state), we only used the largest strongly connected component of the network, which contains 2,484 genes and miRNAs and 32,954 interactions. We excluded self-interactions.

NetICS
We employed NetICS [11] for data integration and network gene ranking. We used differentially expressed miRNAs and somatic mutations (SNV, indels) as upstream causal events. miRNA differential expression was computed in a sample-specific manner by comparing each tumor nodule to the 6 control samples from healthy liver tissue. As downstream events, we used differentially expressed genes/proteins at the RNA, proteome and phosphoproteome levels.
At the phosphoproteome level, one gene was included if there was at least one differentially expressed phosphosite in its protein when tested between tumor and normal tissue. Data at the downstream level were integrated by using the rule described at Table 7.
After we run NetICS, we kept the top 5% of the genes in the ranked list. We performed a random permutation test by permuting the labels of differentially expressed genes, miRNAs and mutated genes for each sample. We then recomputed the gene list and computed an empirical p-value for each gene by counting how many times the score given by NetICS was higher than the original score. We repeated the random permutation procedure 10,000 times and adjusted the p-value by FDR correction [43]. We ended up with 74 genes in total.

Immunoblotting
Both human and murine liver tissue was homogenized in T-PER (ThermoFisher scientific, 78510) supplemented with 1 mM PMSF, 1× Complete Mini Protease Inhibitors (Roche), 1× PhosSTOP (Roche) using a Polytron (PT 10-35 GT) at 500g for 2 min. Equal amounts of homogenate were SDS-PAGE fractionated and transferred onto a nitrocellulose membrane that was incubated, after blocking (5% BSA in TBST), with appropriate antibodies.

Source of human samples
All human samples used in this study were obtained after following the relevant ethical regulations. An informed consent was obtained from the human subjects.