Skip to main content

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

Abstract

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.

Background

Liver cancer is the second leading cause of cancer-related 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 double-knockout (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).

Fig. 1
figure 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-week-old L-dKO mice. Liver tissues from six cre-negative age- and sex- matched littermates were used as a control. This setting was used for mRNA sequencing, miRNA sequencing, proteome and phosphoproteome quantification

Table 1 Shown are the results of exome sequencing in the 12 tumor nodules. Relevant information are given such as the position in thechromosome, the reference and alternative alleles, the type of mutation, the amino acid substitution and variant allele frequency
Fig. 2
figure 2

A. Shown are Venn diagrams for each of the 4 mice demonstrating the common somatic mutations (point mutations and small insertions/deletions) between the tumor samples of the same mouse. We observe no common mutations between the tumor samples, implying that these are more likely passenger mutations. B. Graph showing the number of common mutated genes, differentially expressed mRNAs, C., miRNAs, D. proteins, E., and phosphoproteins, F., across the tumor samples as a cumulative histogram. In each of the subfigures, the number of genes or proteins being dysregulated in at least k samples is shown (y-axis), where k ranges from 1 to 12 (x-axis)

Fig. 3
figure 3

Shown are number of differentially expressed mRNAs, A., and miRNAs, B., proteins, C., and phosphosites, D., in each of the 12 tumor samples. Red color determines up regulation, blue color downregulation and grey color unchanged. mRNAs and miRNAs have been detected from RNA sequencing data and proteins and phosphosites have been detected from mass spectrometry data.

To identify the downstream mediators, we used NetICS, a network-based method that integrates multi-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. Hence, we integrated the Pten and Tsc1 deletions, the somatic mutations, and the differentially expressed miRNAs 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).

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

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.

Table 3 Pathway enrichment results by using Metacore tool

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).

Fig. 4
figure 4

A. Aberrant upstream and downstream interactors of YAP1. Upstream interactors include one significantly downregulated miRNAs (miR-375) and interaction with TSC1 and PTEN. Downstream, YAP1 interacts directly or indirectly with several genes the RNA expression of which has been found significantly upregulated. B. Aberrant upstream and downstream interactors of GRB2 are shown. Upstream interactors include an indirect interaction with PTEN (through PTK2). Downstream, GRB2 interacts with several genes like SHC1, K-RAS, JUN and CDK4. C. miRNAs upstream of HDAC4 significantly downregulated in HCC tumors compared to the control tissue. D. miRNAs upstream of SIRT1 significantly upregulated in HCC tumors compared to the control tissue. The expression counts have been normalized with respect to library sizes and have been transformed for variance stabilization. Shown are normalized expression counts of the miRNAs upstream of Sirt1. The expression counts have been normalized in the same way as in C. Blue color indicates downregulated or genomic deletion, grey color indicates not regulated and red color indicates upregulation

Figure 5
figure 5

A. Immunoblot analysis indicates increased Hdac4 and Yap1 and reduced Sirt1 protein levels in L-dKO tumors (n = 4) compared to age-matched littermate control (liver samples from control mice (n = 4)). B. Quantification of immunoblot (from Fig. 5A) indicates increased Hdac4 (****P = 0.000004) and increased Yap1 (**P = 0.0016) and reduced Sirt1 (**P = 0.0025) expression in tumors compared to age-matched control littermates (band intensities in each lane are normalized to intensity of corresponding total Akt protein levels). P values are from a two-sided unpaired t-test. Data is mean ± s.d. C. Immunoblot analysis indicates increased Hdac4 and Yap1 and reduced Sirt1 protein levels in liver tissue from patients with HCC compared to adjacent non- tumor liver tissue in a total of n = 5 HCC patients. D. Quantification of immunoblot (from Fig. 5C) indicates increased Hdac4 (*P = 0.040032) and increased Yap1 (*P = 0.013164) in 5 out of 5 patients and reduced Sirt1 (**P = 0.007155) in 4 out of 5 patients. 1 out of 5 patients did not show a reduction in Sirt1 protein levels (circles in red) and has been excluded in the significance analysis. (band intensities in each lane are normalized to intensity of corresponding total Akt protein levels). P values are from a two-sided paired t-test. Data are mean ± standard deviation. E. Lis1 mRNA is downregulated in liver cancer. Graphical representation of Lis1 mRNA regulation in n=373 liver cancer patients downloaded from TCGA (provisional) (accessed on 10.10.2018). Each circle represents a patient. 2 fold regulation of mRNA expression compared to control was used to define up regulation (log2 fold change ≥ 1) or down regulation (log2 fold change ≤ 1). Lis1 mRNA expression is downregulated in 47%, unchanged in 43% and upregulated in 10% of the liver cancer patients. The full length images are shown in supplementary figure 5 (Figure S5)

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 upregulation can be attributed to downregulation of PTEN. Consistent with the previous observation that PTEN inhibits PTK2 [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 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 mTOR-driven 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 single-omic 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.

Table 4 Tumor samples are compared against control samples at the RNA (column B), proteome (column C) and phosphoproteome (column D) levels for the 74 predicted mediators. The gene is indicated as “dysregulated”, if it is ranked at the top 5% of all genes based on P-value

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 malignancies - drives 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.

Methods

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.

Table 5 Statistics of whole-exome sequencing

Whole-exome sequencing data pre-processing was performed as described in Nuciforo et al, 2018 against the reference mouse genome GRCm38. In brief, paired-end reads in FASTQ format were aligned to the reference mouse genome GRCm38 using Burrows-Wheeler Aligner (v0.7.12) [33]. Local realignment was performed using the Genome Analysis Toolkit (GATK, v3.6) [34]. PCR duplicates were removed using Picard (v2.4.1, http://broadinstitute.github.io/picard/). Base quality adjustment was performed using GATK (v3.6) [34].

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. Primer sequences:

Table 6 Number of upregulated, downregulated and unchanged mRNA, miRNAs, proteins and phosphosites per tumor nodule

TBP F: ATGATGCCTTACGGCACAGG; R: GTTGCTGAGATGTTGATTGCTG;

CYR61 F: TAAGGTCTGCGCTAAACAACTC; R: CAGATCCCTTTCAGAGCGGT;

KI67 F: CGCAGGAAGACTCGCAGTTT; R: CTGAATCTGCTAATGTCGCCAA

CTGF F: GGCCTCTTCTGCGATTTCG; R: GCAGCTTGACCCTTCTCGG

BIRC5 F: GAGGCTGGCTTCATCCACTG; R: ATGCTCCTCTATCGGGTTGTC

CYCE1 F: CTCCGACCTTTCAGTCCGC; R: CACAGTCTTGTCAATCTTGGCA

miRNA sequencing

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 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.

Table 7 Combination rules for differentially expressed genes

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.

Antibodies

Yap1 ((G-6) sc-376830), HDAC4 (CST 7628), Sirt1 (CST 3931) and total AKT (CST, 9272) were obtained commercially. Horseradish peroxidase (HRP)-coupled anti-mouse (115-035- 774) and anti-rabbit (211-032-171) secondary antibodies were purchased from Jackson laboratories.

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.

Availability of data and materials

The DNA and RNA raw murine data analyzed in this study is available in the Sequence Read Archive repository, under the accession SRP156216 [https://www.ncbi.nlm.nih.gov/sra/?term=SRP156216]. The raw proteomic and phosphoproteomic murine data as well as the interaction network used for the analysis are available in our github repository [https://github.com/cbg-ethz/netics/tree/master/mouse_data]. The reference mouse genome GRCm38 was downloaded from [https://www.ncbi.nlm.nih.gov/assembly/GCF_000001635.20/].

Abbreviations

HCC:

Hepatocellular carcinoma

NetICS:

Network-based integration of multi-omics data

miRNA:

Micro RNA

mRNA:

Messenger RNA

L-dKO:

Liver-specific double-knockout

PCR:

Polymerase chain reaction

qPCR:

Quantitative polymerase chain reaction

FDR:

False discovery rate

SNV:

Single nucleotide variant

Indel:

Insertion/deletion

EMT:

Epithelial to mesenchymal transition

TCGA:

The cancer genome atlas

References

  1. Llovet JM, et al. Hepatocellular carcinoma. Nat Rev Dis Primers. 2016;2:16018–23.

    Article  Google Scholar 

  2. Hindupur SK, et al. The protein histidine phosphatase LHPP is a tumour suppressor. Nature. 2018;555(7698):678–82.

    CAS  Article  Google Scholar 

  3. Guri Y, et al. mTORC2 promotes tumorigenesis via lipid synthesis. Cancer Cell. 2017;32(6):807–823.e12.

    CAS  Article  Google Scholar 

  4. Gerlinger M, et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing. N Engl J Med. 2012;366(10):883–92.

    CAS  Article  Google Scholar 

  5. Peng Y, Croce CM. The role of MicroRNAs in human cancer. Signal Transduct Target Ther. 2016;1:15004.

  6. Yu CB, et al. Systemic transcriptome analysis of hepatocellular carcinoma. Tumor Biol. 2016;37(10):13323–31.

    CAS  Article  Google Scholar 

  7. Jiang Y, et al. Proteomics identifies new therapeutic targets of early-stage hepatocellular carcinoma. Nature. 2019;567(7747):257–61.

    CAS  Article  Google Scholar 

  8. Dazert E, et al. Quantitative proteomics and phosphoproteomics on serial tumor biopsies from a sorafenib-treated HCC patient. Proc Natl Acade Sci. 2016;113(5):1381–6.

    CAS  Article  Google Scholar 

  9. Ally A, et al. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell. 2017;169(7):1327–1341.e23.

    Article  Google Scholar 

  10. Mertins P, et al. Proteogenomics connects somatic mutations to signalling in breast cancer. Nature. 2016;534:55–62.

    CAS  Article  Google Scholar 

  11. Dimitrakopoulos C, et al. Network-based integration of multi-omics data for prioritizing cancer genes. Bioinformatics. 2018;34(14):2441–8.

    CAS  Article  Google Scholar 

  12. Villanueva A, et al. Notch signaling is activated in human hepatocellular carcinoma and induces tumor formation in mice. Gastroenterology. 2012;143(6):1660–1669.e7.

    CAS  Article  Google Scholar 

  13. He G, et al. Identification of liver cancer progenitors whose malignant progression depends on autocrine IL-6 signaling. Cell. 2013;155(2):384–96.

    CAS  Article  Google Scholar 

  14. Sun Q, Jiang N, Sun R. Leptin signaling molecular actions and drug target in hepatocellular carcinoma. Drug Design Dev Ther. 2014;8:2295–302.

    Article  Google Scholar 

  15. Zhao JF, et al. The ASH1-miR-375-YWHAZ signaling axis regulates tumor properties in hepatocellular carcinoma. Mol Ther Nucleic Acids. 2018;11:538–53.

    CAS  Article  Google Scholar 

  16. Sun H, et al. Loss of KLF4 and consequential downregulation of smad7 exacerbate oncogenic tgfbeta signaling in and promote progression of hepatocellular carcinoma. Oncogene. 2017;36(21):2957–68.

    CAS  Article  Google Scholar 

  17. Moroishi T, et al. The emerging roles of YAP and TAZ in cancer. Nat Rev Cancer. 2015;15(2):73–9.

    CAS  Article  Google Scholar 

  18. Kang ZH, et al. Histone deacetylase HDAC4 promotes gastric cancer SGC-7901 cells progression via p21 repression. PLoS One. 2014;9(6):e98894.

    Article  Google Scholar 

  19. Tamura M, Gu Jenny X, Yamada KM. Tumor suppressor pten inhibition of cell invasion, migration, and growth: differential involvement of focal adhesion kinase and p130cas. Cancer Res. 1999;59(2):442–9.

    CAS  PubMed  Google Scholar 

  20. Bollard J, et al. Palbociclib (PD-0332991), a selective CDK4/6 inhibitor, restricts tumour growth in preclinical models of hepatocellular carcinoma. Gut. 2016;66(7):1286–96.

    Article  Google Scholar 

  21. Chetram MA, Hinton CV. PTEN regulation of ERK1/2 signaling in cancer. J Recept Signal Transduct. 2012;32(4):190–5.

    CAS  Article  Google Scholar 

  22. Wilson AJ, et al. HDAC4 promotes growth of colon cancer cells via repression of p21. Mol Biol Cell. 2008;19(10):4062–75.

    CAS  Article  Google Scholar 

  23. Mottet D, et al. HDAC4 represses p21waf1/cip1 expression in human cancer cells through a sp1-dependent, p53-independent mechanism. Oncogene. 2008;28(2):243–56.

    Article  Google Scholar 

  24. Shen YF, Wei AM, Kou Q, Zhu QY, Zhang L. Histone deacetylase 4 increases progressive epithelial ovarian cancer cells via repression of p21 on fibrillar collagen matrices. Oncol Rep. 2015;35(2):948–54.

    Article  Google Scholar 

  25. Zeng LS, et al. Overexpressed HDAC4 is associated with poor survival and promotes tumor progression in esophageal carcinoma. Aging. 2016;8(6):1236–48.

    Article  Google Scholar 

  26. Li L, et al. Activation of p53 by SIRT1 inhibition enhances elimination of CML leukemia stem cells in combination with imatinib. Cancer Cell. 2012;21(2):266–81.

    CAS  Article  Google Scholar 

  27. Chalkiadaki A, Guarente L. The multifaceted functions of sirtuins in cancer. Nat Rev Cancer. 2015;15(10):608–24.

    CAS  Article  Google Scholar 

  28. Moon HM, et al. LIS1 controls mitosis and mitotic spindle organization via the LIS1–NDEL1–dynein complex. Hum Mol Genet. 2014;23(2):449–66.

    CAS  Article  Google Scholar 

  29. Li X, et al. Hepatic loss of Lissencephaly 1 (Lis1) induces fatty liver and accelerates liver tumorigenesis in mice. J Biol Chem. 2018;293(14):5160–71.

    CAS  Article  Google Scholar 

  30. Wang S, et al. Dysregulation of miR484-TUSC5 axis takes part in the progression of hepatocellular carcinoma. J Biochem. 2019;166(3):271–9.

    CAS  Article  Google Scholar 

  31. Wanke J, et al. Pharmacological Inhibition of Class IIA HDACs by LMK-235 in pancreatic neuroendocrine tumor cells. Int J Mol Sci. 2018;19(10):3128.

    Article  Google Scholar 

  32. Kandasamy S, et al. The YAP1 Signaling Inhibitors, Verteporfin and CA3, suppress the mesothelioma cancer stem cell phenotype. Mol Cancer Res. 2020;18(3):343–51.

    CAS  Article  Google Scholar 

  33. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–60.

    CAS  Article  Google Scholar 

  34. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.

    CAS  Article  Google Scholar 

  35. Cibulskis K, Lawrence MS, Carter SL, Sivachenko A, Jaffe D, Sougnez C, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol. 2013;31:213–9.

    CAS  Article  Google Scholar 

  36. Saunders CT, Wong WS, Swamy S, Becq J, Murray LJ, Cheetham RK. Strelka: accurate somatic small-variant calling from sequenced tumor-normal sample pairs. Bioinformatics. 2012;28:1811–7.

    CAS  Article  Google Scholar 

  37. Nuciforo S, et al. Organoid models of human liver cancers derived from tumor needle biopsies. Cell Reports. 2018;24(5):1363–76.

    CAS  Article  Google Scholar 

  38. Shen R, Seshan VE. FACETS: allele-specific copy number and clonal heterogeneity analysis tool for high-throughput DNA sequencing. Nucleic Acids Res. 2016;44:e131.

    Article  Google Scholar 

  39. Jean G, Kahles A, Sreedharan VT, De Bona F, Rätsch G. Rna-seq read alignments with palmapper. Curr Protoc Bioinformatics. 2010; Chapter 11:Unit 11.6. https://pubmed.ncbi.nlm.nih.gov/21154708/.

  40. Hafner M, et al. Identification of microRNAs and other small regulatory RNAs using cDNA library sequencing. Methods. 2008;44(1):3–12.

    CAS  Article  Google Scholar 

  41. Tyanova S, et al. The perseus computational platform for comprehensive analysis of (prote) omics data. Nat Methods. 2016;13(9):731–40.

    CAS  Article  Google Scholar 

  42. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  Google Scholar 

  43. Yekutieli D, et al. Resampling-based false discovery rate controlling multiple test procedures for correlated test statistics. J Stat Plan Inf. 1999;82:171–96.

    Article  Google Scholar 

Download references

Acknowledgements

Not applicable

Funding

The Swiss Institute of Bioinformatics (SIB) PhD Fellowship Programme, European Research Council (ERC) Synergy Grant 609883, SystemsX.ch Research, Technology and Development (RTD) Grant 2013/150 and the European Commission (EC) Horizon 2020 project 633974 SOUND, have supported this work. The Swiss Institute of Bioinformatics (SIB) PhD Fellowship Programme scholarship has funded the PhD of Christos Dimitrakopoulos. The rest of the funding sources funded HCC-specific studies. The funders had no role in the design, execution, or analysis of the study.

Author information

Authors and Affiliations

Authors

Contributions

CD, SKH, NB and MH conceived the idea and designed the analysis. CD and SKH analyzed the data, implemented the network analysis and interpreted the results. NB and MH have substantially revised the manuscript. SKH and MC performed the proteomics data analysis. DM performed the western blot analysis. LT and MM performed histological sectioning and analysis of the L-dKO murine liver tissues, obtained ethics approval for patient samples and provided patient material. DL performed the RNA rt PCR analysis. CKYN, SP, JB, ALM, JS and HLR performed the DNA sequencing analysis. JB performed the RNA sequencing analysis and DL performed the miRNA sequencing analysis. All authors have read and approved the manuscript.

Corresponding authors

Correspondence to Michael N. Hall or Niko Beerenwinkel.

Ethics declarations

Ethics approval and consent to participate

All relevant ethical regulations were followed in this study. An informed consent was obtained from the human subjects. The collection and analysis of patient samples was approved by the ethics committee of Northern Switzerland (EKNZ, Study ID: 310/12/PB_2019-00154/Genetische Analyse des Leberzellkarzinoms). Samples collected after 2016 were obtained under a written consent by the patient. Before 2016, samples were included if no active refusal of the patient existed, which was approved by the ethics committee. The animal license number was 2555. For two samples, we had no refusal from the patient to use the tissue for research purpose and because the tissue derived from 1999 and 2002 only an approval from the Ethics committee (EKNZ, Basel Switzerland) was needed at that time by the Swiss regulation to use the tissue for research purpose. For all remaining samples, a written informed consent of the patient was present to use the tissue for research purposes, as well as an approval from the Ethics committee (EKNZ, Basel Switzerland).

Consent for publication

Not applicable

Competing interests

There is no financial or non-financial competing interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Copy number profiles derived from whole-exome sequencing demonstrates the intragenic deletions of Tsc1 and Pten. For each nodule, segmented Log2 ratios (y-axis) were plotted according to their genomic positions (x-axis), for chromosomes 2 or 19. Red arrows indicate the loci of the intragenic deletions of Tsc1 and Pten.

Additional file 2: Figure S2. A.

Graph showing miRNA expression of miR-375 in L-dKO tumors (n=12) compared to control mice (n=6). B. mRNA expression analysis of indicated genes in 20-week-old L-dKO tumors compared to livers from age-matched control mice (n = 6). Expression for each gene is normalized to intensity of TBP gene expression (normalising control) in the corresponding mice. Two-sided unpaired t-test is used. Data is mean ± s.d.

Additional file 3: Figure S3.

Graph showing the expression of miRNAs (upstream of HDAC4) in L-dKO tumors (n=12) compared to control mice (n=6).

Additional file 4: Figure S4.

Graph showing the expression of miRNAs (upstream of SIRT1) in L-dKO tumors (n=12) compared to control mice (n=6).

Additional file 5: Figure S5.

Original western blots full length images of the images shown in Figures 5A and C.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Dimitrakopoulos, C., Hindupur, S.K., Colombi, M. et al. Multi-omics data integration reveals novel drug targets in hepatocellular carcinoma. BMC Genomics 22, 592 (2021). https://doi.org/10.1186/s12864-021-07876-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-021-07876-9

Keywords

  • HCC
  • mTOR signaling
  • NetICS
  • Omics