Transcriptomics-Based Drug Repurposing Approach Identifies Novel Drugs against Sorafenib-Resistant Hepatocellular Carcinoma

Simple Summary Hepatocellular carcinoma (HCC), a type of liver cancer, remains a treatment challenge due to late detection and resistance to currently approved drugs. It takes 15–20 years for a single new drug to become FDA approved. The purpose of this study was to expedite identification of novel drugs against drug-resistant HCC. For this, we matched gene expression alterations in resistant HCC with gene expression changes caused by treatment of cancer cells with drugs already FDA approved for other diseases to find the drug that can reverse the resistance-related changes. Among the identified drugs, we validated the growth inhibitory effect of two drugs, identified their mechanism in HCC and, thus, provided proof of concept evidence for validity of this drug repurposing approach with potential for use in personalized medicine. Abstract Objective: Hepatocellular carcinoma (HCC) is frequently diagnosed in patients with late-stage disease who are ineligible for curative surgical therapies. The majority of patients become resistant to sorafenib, the only approved first-line therapy for advanced cancer, underscoring the need for newer, more effective drugs. The purpose of this study is to expedite identification of novel drugs against sorafenib resistant (SR)-HCC. Methods: We employed a transcriptomics-based drug repurposing method termed connectivity mapping using gene signatures from in vitro-derived SR Huh7 HCC cells. For proof of concept validation, we focused on drugs that were FDA-approved or under clinical investigation and prioritized two anti-neoplastic agents (dasatinib and fostamatinib) with targets associated with HCC. We also prospectively validated predicted gene expression changes in drug-treated SR Huh7 cells as well as identified and validated the targets of Fostamatinib in HCC. Results: Dasatinib specifically reduced the viability of SR-HCC cells that correlated with up-regulated activity of SRC family kinases, its targets, in our SR-HCC model. However, fostamatinib was able to inhibit both parental and SR HCC cells in vitro and in xenograft models. Ingenuity pathway analysis of fostamatinib gene expression signature from LINCS predicted JAK/STAT, PI3K/AKT, ERK/MAPK pathways as potential targets of fostamatinib that were validated by Western blot analysis. Fostamatinib treatment reversed the expression of genes that were deregulated in SR HCC. Conclusion: We provide proof of concept evidence for the validity of this drug repurposing approach for SR-HCC with implications for personalized medicine.


Introduction
Liver cancer ranks seventh among the most common cancers in the world and, according to a recent report [1], it is the fourth leading cause of cancer-related death. As per the prediction of The American Cancer Society,~42,810 new individuals (30,170 in men and 12,640 in women) will be diagnosed with primary hepatocellular cancer and intrahepatic bile duct cancer and about 30,160 patients (20,020 men and 10,140 women) will die of these cancers in the United States in 2020 (https://www.cancer.org/cancer/liver-cancer.html). Since 1980, liver cancer incidence rates have more than tripled and death rates have more than doubled.
While a small proportion of hepatocellular carcinoma (HCC) patients diagnosed at an early stage can be treated by tumor resection, cryoablation or liver transplant, these treatments are not effective in the majority of HCC patients diagnosed at an advanced stage of the disease. Sorafenib, a multi-kinase inhibitor, is the most widely used drug since 2007 in treating such patients [2]. However, the median overall survival of sorafenib treated patients is only extended by 2.8 months compared to untreated patients [3]. This minimal therapeutic response is attributed to HCC tumors having an intrinsic resistance to the cytostatic effects of sorafenib [4]. For HCC patients who become resistant to sorafenib, FDA recently approved its congener regorafenib [5] and checkpoint blockade anti-PD-1 antibodies, nivolumab [6] and pembrolizumab [7], as second-line treatment. However, only a subset of such patients respond to this combination therapy. In recent clinical trials in the first-line setting, nivolumab [8] or pembrolizumab [9] could not significantly improve survival of HCC patients compared to sorafenib and best supportive care, respectively. Lenvatinib, another multi-kinase inhibitor, recently approved as first-line therapy for advanced HCCs, was not significantly superior to sorafenib in improving overall survival in clinical trials [10]. Thus, there is an urgent need to develop therapeutic strategies to overcome sorafenib resistance and discover new, more effective therapies.
Due to the heterogeneous molecular mechanisms underlying HCC tumor progression and sorafenib resistance, it is unlikely that targeting one major molecular mechanism will be sufficient to treat this disease [11]. Furthermore, collecting tissues through biopsy is not standard of care for advanced HCC patients that relapse on sorafenib. As a result, the lack of available biomolecular data of HCC patients treated with sorafenib severely limits the ability to study fundamental mechanisms of resistance and potential targets for combination therapies. Therefore, we sought to investigate sorafenib resistance in HCC in an unbiased way through global analysis of the transcriptome in experimental models of sorafenib resistance. We generated gene expression data from an in vitro model of HCC sorafenib resistance in the Huh7 cell line, and conducted a comprehensive analysis of other publicly available gene expression data from experimental models of sorafenib resistance (SR-HCC) and patient-derived HCC tumors. We evaluated these SR-HCC gene expression models for their coverage in human HCC tissue samples and their prognostic significance. Next, in order to discover drug-disease hypotheses, we utilized the aforementioned gene expression profiles as the basis for our computational drug repurposing analyses via connectivity mapping. Connectivity mapping uses pattern-matching algorithms to compare genome-wide gene expression changes observed in cultured human cells treated with drugs to those of biological states of interest: e.g., tumor vs. normal [12]. Connectivity scores quantify the drug-disease hypotheses through correlations between ranked gene lists of query gene signatures and drug reference gene signatures, commonly via the Kolgomorov-Smirnov statistic or the modified gene set enrichment analysis method [12,13]. For instance, drug-induced gene hypothesized to reverse or oppose the query gene signature characterizing a disease, and vice versa. Furthermore, the use of genome-wide expression profiles provides mechanistic insight into tumor biology and drug efficacy, which may be missed by other guilt-by-association approaches. In this study, we hypothesize that transcriptomics data from experimental models of sorafenib-resistant HCC (i) will enable validation of the in vitro models in the absence of tissue available from sorafenib resistant tumors, (ii) can be applied in connectivity mapping studies to predict novel therapies to curb resistance to sorafenib in HCC, and (iii) will reveal molecular mechanisms underlying sorafenib resistance in HCC tumors in the context of gene targets and enriched pathways.

Results
A workflow overview for this study is presented in Figure 1. Overview of computational drug repurposing workflow and drug target network. Gene expression signatures of experimental models of hepatocellular carcinoma (HCC) sorafenib resistance were (i) assessed for prognostic significance, and (ii) queried against gene expression signatures characterizing drug perturbations in the HepG2 cell line contained in the Library of Integrated Network-based Cellular Signatures (LINCS) database. Connectivity scores were calculated by the rank-based, non-parametric weighted Kolgomorov-Smirnov (KS) statistic. Drugs with negative connectivity scores (i.e., anti-correlated) represent those that are hypothesized to reverse HCC sorafenib resistance gene signature. Drug candidates were further prioritized based on FDA approval/clinical investigation status, known anti-neoplastic activity and literature evidence for drug target genes associated with HCC. Two drugs were subsequently selected for in vitro validation.

Generation and Microarray Analysis of Sorafenib-Resistant HCC Cell Line
Sorafenib-resistant HCC cell lines were generated from parental (sensitive) Huh7 cells (Huh7-S) following long-term exposure to sorafenib [14]. Viability assays demonstrated that the IC50 doses for resistant cells were approximately 5-fold higher than that of parental Huh-7 cells ( Figure S1). Microarray Overview of computational drug repurposing workflow and drug target network. Gene expression signatures of experimental models of hepatocellular carcinoma (HCC) sorafenib resistance were (i) assessed for prognostic significance, and (ii) queried against gene expression signatures characterizing drug perturbations in the HepG2 cell line contained in the Library of Integrated Network-based Cellular Signatures (LINCS) database. Connectivity scores were calculated by the rank-based, non-parametric weighted Kolgomorov-Smirnov (KS) statistic. Drugs with negative connectivity scores (i.e., anti-correlated) represent those that are hypothesized to reverse HCC sorafenib resistance gene signature. Drug candidates were further prioritized based on FDA approval/clinical investigation status, known anti-neoplastic activity and literature evidence for drug target genes associated with HCC. Two drugs were subsequently selected for in vitro validation.

Generation and Microarray Analysis of Sorafenib-Resistant HCC Cell Line
Sorafenib-resistant HCC cell lines were generated from parental (sensitive) Huh7 cells (Huh7-S) following long-term exposure to sorafenib [14]. Viability assays demonstrated that the IC 50 doses for Cancers 2020, 12, 2730 4 of 22 resistant cells were approximately 5-fold higher than that of parental Huh-7 cells ( Figure S1). Microarray analysis was performed on the parental cells (Huh7-S) and sorafenib-resistant cells (Huh7-R and Huh7-R-A7) (GSE94550). Comparing Huh7-R-A7 vs. Huh-S cells to define an HCC sorafenib resistance gene signature, we determined 368 genes to be differentially expressed (adjusted p < 0.001) ( Table S1). The top Gene Ontology (GO) molecular functions, biological processes and cellular components enriched in the sorafenib resistance gene signature (FDR < 0.05) are reported in Table S1.

Evaluation of Experimental Models of HCC Sorafenib Resistance against HCC Patient Datasets
We conducted a comprehensive analysis of publicly available gene expression datasets characterizing sorafenib resistance in HCC, including an in vitro HepG2 cell line (HepG2-R) [15], in vitro HCC patient culture (HCC-3sp-R) [16], and an in vivo xenograft model using transplanted Huh7 cells (Xeno-R) [17]. Sorafenib resistance (SR) dataset features are described in Table 1, and gene signature overlap is shown in Figure S2. Interestingly, the distinct gene signatures showed very little overlap among up-and down-regulated genes. To assess the clinical prognostic relevance of each of the SR gene signatures, we obtained additional publicly available gene expression datasets characterizing HCC patient tumors of diverse etiologies from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. First, the presence of the SR gene signatures within the HCC patient datasets was detected using the nearest template method (weighted cosine similarity metric, FDR < 0.05). The ability for the gene signatures to distinguish between HCC tumor and normal liver tissue was comparable across datasets. However, the average frequency of SR+ HCC tumors defined by the SR gene signatures among the GEO datasets were highest in the Huh7-R-A7 (52.7%) and Xeno-R (46.4%) signatures, while HepG2-R (26.3%) and HCC-3sp-R (27.7%) were lower (Table 1). Similarly, the Huh7-R-A7 and Xeno-R gene signatures detected a higher percentage of SR+ samples using RNAseq data from the TCGA liver hepatocellular carcinoma (LIHC) subset ( Figure 2A). Second, we assessed the relationship between the presence of a sorafenib resistance gene signature and survival using the TCGA dataset. Of the four SR gene signatures, HCC patients with tumors harboring the Huh7-R-A7 sorafenib resistance gene signature (SR+) exhibited significantly reduced survival as compared to those that showed the opposite gene signature pattern (SR−) (log-rank p = 0.0086; HR = 1.59, 95% CI: 1.13-2.45) ( Figure 2B). Due to the lack of survival outcome data for HCC patients with gene expression data from the GEO database, we assessed the sensitivity and specificity by which the four SR gene signatures could distinguish HCC tumor vs. normal liver tissue. The average sensitivity and specificity for the SR signatures are as follows: Huh7-R-A7 (0.50, 0.82), HepG2-R (0.53, 0.89), HCC-3sp-R (0.48, 0.89), Xeno-R (0.73, 0.80).    To identify drugs that can reverse sorafenib resistance in HCC, connectivity mapping analyses were conducted via the Library of Integrated Network-based Cellular Signatures (LINCS) L1000 system. This database consists of 476,251 gene expression profiles of drug and genetic perturbation conditions across 77 cellular contexts. We chose to focus on drug-treated profiles from the HepG2 cell line, as it represented the HCC cell line with the greatest number of unique drugs (n = 3740). To determine whether HCC cell line context can affect the distribution of drug repurposing hypotheses, we first used gene expression profiles for 18 HCC cell lines from the CellMinerHCC database to query HepG2 drug perturbation gene expression profiles from the LINCS database. We then performed hierarchical clustering analysis of connectivity scores for drug predictions across the 18 HCC cell lines. The heatmap shown in Figure 2C revealed two distinct clusters of drug connectivity scores. Interestingly, both HepG2 and Huh7 belonged to the same main cluster branch, suggesting that the LINCS HepG2 represents a suitable system to generate drug predictions from gene expression profiles originating in Huh7 cell line models. Next, we compared the connectivity scores across the four HCC sorafenib resistance models via hierarchical clustering analysis ( Figure 2D). We observed that drug prediction patterns derived from the two cell lines (Huh7-R-A7 and HepG2-R) and short-term culture (HCC-3sp-R) were nearly opposite those derived from the mouse model (Xeno-R). The highest degree of similarity of drug prediction patterns was between HCC-3sp-R and HepG2-R. These results are consistent with the observed gene signature overlap shown in Figure S2. For instance, the maximum observed gene overlap between HCC-3sp-R and HepG2-R was 88 up-regulated and 189 down-regulated genes, while Xeno-R showed the least overlap with any other model.

Drug Target Analysis and Candidate Prioritization
Due to the superiority of the Huh7-R-A7 sorafenib resistance gene signature to distinguish significant patient survival patterns using the TCGA LIHC data ( Figure 2B), we focused our prioritization and validation efforts for LINCS drug predictions from this gene signature. Amongst the drugs predicted to reverse sorafenib resistance from our LINCS analyses, we first determined the approval and investigational status from the DrugBank and Aggregate Analysis of ClinicalTrials.gov (AACT) databases to prioritize those drugs that would be most feasible for future preclinical and clinical testing (Table S2). For these drugs, we annotated drug activity, target pathway and gene information from KEGG and DrugBank databases. We then further selected only those drugs with known anti-neoplastic activity for initial validation. Finally, we conducted a systematic search of genes associated with hepatocellular carcinoma using a publicly available literature mining tool [18] to determine whether the genes targeted by drug candidates had known roles in HCC. The final prioritized drug list is shown in Table 2. We selected two drug candidates, one representative from each approval status, for subsequent validation: 1) dasatinib, SRC family of kinases inhibitor (FDA-approved) and 2) fostamatinib, SYK inhibitor (under clinical investigation). Additionally, we analyzed a protein-protein interaction (PPI) network of drug target genes, as shown in Figure S3. The initial drug targets were connected in a network of a total of 322 protein nodes through 1029 total interactions (edges). The average node degree (interaction partners) is 6.39, and the number of observed edges is significantly higher than expected (n = 475; PPI enriched p value < 0.05). PPI connections were recovered for five out of seven dasatinib gene targets, including SRC. A community detection network algorithm was applied to the network, and both SRC and SYK were found in the largest PPI module 5 (n = 46 total nodes). Notably, SRC was ranked with the third highest node degree (n = 33 connections), while SYK has n = 8 connections. SYK exhibited a clustering coefficient of 0.57, while SRC has a clustering coefficient of 0.14. Taken together, these results suggest SRC inhibition may have a broad impact on the overall network, while SYK appears to be involved with a more tightly regulated group of proteins. Finally, we selected the top 100 central genes based on the eigencentrality measure, which included all dasatinib and fostamatinib gene targets that were in the network. We found that these highly central drug target genes were altered (mutations, copy number variations, mRNA and protein expression levels) at a higher frequency in the SR+ TCGA patients (58%) as compared to the SRpatients (31%). The type and frequency of gene alterations of the top central genes in the PPI network in SR+ and SR-TCGA LIHC patient tumors are shown in Table S3. Dasatinib and fostamatinib were initially tested in vitro as single agents in HCC cell lines (parental Huh7-S, resistant pool Huh7-R and resistant clone Huh7-R-A7). Parental and sorafenib-resistant Huh7 cells were treated with increasing concentrations of dasatinib and fostamatinib independently, and cellular viability was assessed after 48 h. Sorafenib-resistant Huh7 cells were significantly more sensitive to dasatinib toxicity than parental cells ( Figure 3A). Parental cells displayed an IC 50 of > 60 µM, while the IC 50 of resistant cells was < 10 µM. On the contrary, all three cell lines displayed similar sensitivity to fostamatinib toxicity with IC 50 values between 20 and 35 µM ( Figure 3B), indicating that fostamatinib may be useful both before and after resistance to sorafenib occurs. We next hypothesized that combining either dasatinib or fostamatinib with sorafenib would synergistically inhibit cell viability.
Cancers 2020, 12, x; doi: www.mdpi.com/journal/cancers We thus assessed the effect of dasatinib and fostamatinib treatment, alone and in combination with sorafenib, on the reproductive ability of single cells using the clonogenic survival assay ( Figure 3C,D). At the low (2 µM) concentration tested, dasatinib completely inhibited colony formation in Huh7-R cells, while having little effect on parental cells. This observation is consistent with the higher sensitivity of resistant cells to dasatinib observed in viability assay ( Figure 3A). Furthermore, the addition of sorafenib did not further sensitize the Huh7-S cells to dasatinib. Similar to the viability assay, both parental and resistant cells were significantly sensitive to fostamatinib. The combination of fostamatinib and sorafenib appeared to have a greater ability to inhibit colony formation in Huh7-S, Huh7-R and Huh7-R-A7 cells compared to either drug alone.
Since sorafenib-resistant HCC cells appeared to be much more sensitive to long-term dasatinib toxicity than non-resistant HCC cells, we hypothesized that the activity of SRC kinases would be upregulated in the resistant cells. Kinase array analysis demonstrated that five out of seven SRC kinases (Src, Yes, Fyn, Lck and Lyn) were significantly hyper-phosphorylated in resistant cells as compared to parental cells ( Figures 4A and S4). On the contrary, SYK protein, the target of fostamatinb was not detected in the Huh7 cells. Therefore, in order to delineate possible mechanisms for the action of fostamatinib, we obtained the gene expression characterizing fostamatinib-treated HepG2 cells (10 µM, 6 h post-treatment) from the LINCS database, and performed Ingenuity pathway analysis. Several pathways were identified as significantly deregulated post-treatment: inhibition of oncogenic signaling through JAK/STAT (p < 1.0 × 10 −8 ), protein kinase A (p < 1.0 × 10 −7 ), PI3K/AKT (p < 1.0 × 10 −4 ), STAT3 (p < 1.0 × 10 −4 ) and ERK/MAPK (p < 1.0 × 10 −3 ); and activation of the tumor suppressor PTEN ( Figure 4B).  We thus assessed the effect of dasatinib and fostamatinib treatment, alone and in combination with sorafenib, on the reproductive ability of single cells using the clonogenic survival assay ( Figure 3C,D). At the low (2 µM) concentration tested, dasatinib completely inhibited colony formation in Huh7-R cells, while having little effect on parental cells. This observation is consistent with the higher sensitivity of resistant cells to dasatinib observed in viability assay ( Figure 3A). Furthermore, the addition of sorafenib did not further sensitize the Huh7-S cells to dasatinib. Similar to the viability assay, both parental and resistant cells were significantly sensitive to fostamatinib. The combination of fostamatinib and sorafenib appeared to have a greater ability to inhibit colony formation in Huh7-S, Huh7-R and Huh7-R-A7 cells compared to either drug alone.
Since sorafenib-resistant HCC cells appeared to be much more sensitive to long-term dasatinib toxicity than non-resistant HCC cells, we hypothesized that the activity of SRC kinases would be up-regulated in the resistant cells. Kinase array analysis demonstrated that five out of seven SRC kinases (Src, Yes, Fyn, Lck and Lyn) were significantly hyper-phosphorylated in resistant cells as compared to parental cells ( Figure 4A and Figure S4). On the contrary, SYK protein, the target of fostamatinb was not detected in the Huh7 cells. Therefore, in order to delineate possible mechanisms for the action of fostamatinib, we obtained the gene expression characterizing fostamatinib-treated HepG2 cells (10 µM, 6 h post-treatment) from the LINCS database, and performed Ingenuity pathway analysis. Several pathways were identified as significantly deregulated post-treatment: inhibition of oncogenic signaling through JAK/STAT (p < 1.0 × 10 −8 ), protein kinase A (p < 1.0 × 10 −7 ), PI3K/AKT (p < 1.0 × 10 −4 ), STAT3 (p < 1.0 × 10 −4 ) and ERK/MAPK (p < 1.0 × 10 −3 ); and activation of the tumor suppressor PTEN ( Figure 4B).

Characterization of Fostamatinib as Anti-HCC Drug
Fostamatinib, recently FDA approved for immune thrombocytopenia, has been clinically tested for autoimmune diseases and lymphoma. However, its anti-HCC efficacy has not been explored previously. Since fostamatinib could inhibit the growth of both sorafenib-sensitive and -resistant cells, we tested its inhibitory potential in several HCC cell lines with differing sensitivities to sorafenib. The IC50 of fostamatinib for most of these cell lines was 2-4 uM ( Figure S5). To confirm the in vivo efficacy of fostamatinib, we used a subcutaneous xenograft model of MHCCLM3 cells that inherently have reduced sensitivity to sorafenib. The growth of tumors was significantly inhibited in fostamatinib-treated mice compared to vehicle-treated mice ( Figure 5A). Comparable body weights of the tumor-bearing mice treated with vehicle or fostamatinib ( Figure 5A) indicate no systemic toxicity of the drug. Tumors harvested at the end of 28 days of treatment weighed an average of 3.2 g in vehicle-treated mice and 1.8 g in fostamatinib-treated mice ( Figure 5A). However, due to the small sample size, this difference was not statistically significant (p = 0.1031).
Based on the IPA analyses of the gene expression profile of fostamatinib-treated HepG2 cells from LINCS ( Figure 4B), we performed immunoblot analysis to validate the kinases that are inhibited by fostamatinib. Our data demonstrate that fostamatinib inhibits phosphorylation of ERK, AKT and STAT3 but has no effect on PTEN. Additionally, we observed inhibitory effect on EGFR and JNK ( Figure 5B).
Finally, we sought to validate predicted gene expression changes associated with fostamatinib treatment from LINCS in our Huh7 cell line model. We performed RNA-seq analysis of fostamatinib-vs. DMSO vehicle-treated sorafenib-resistant Huh7-R cells and determined the concordance of fold changes for differentially expressed genes (adjusted p < 0.05) with both the LINCS fostamatinib gene expression signature and Huh7-R sorafenib resistance gene expression signature ( Figure S6). Prioritized differentially expressed genes that were concordant with the LINCS fostamatinib signature and discordant with the Huh7-R sorafenib resistance signature are shown in Table S4. We observed 18 up-regulated and 9 downregulated genes significantly differentially expressed genes in the Huh7-R signature that were reversed in the fostamatinib treated Huh7 cells and LINCS fostamatinib HepG2 signature. Several examples of genes implicated in HCC tumorigenesis and/or sorafenib resistance that were up-regulated in the Huh7-R cells and down-regulated following fostamatinib treatment in both Huh7-R and HepG2 HCC cells include transforming growth factor beta 2 (TGFB2) [19,20], hypoxia inducible factor 1 alpha subunit (HIF1A) [21],

Characterization of Fostamatinib as Anti-HCC Drug
Fostamatinib, recently FDA approved for immune thrombocytopenia, has been clinically tested for autoimmune diseases and lymphoma. However, its anti-HCC efficacy has not been explored previously. Since fostamatinib could inhibit the growth of both sorafenib-sensitive and -resistant cells, we tested its inhibitory potential in several HCC cell lines with differing sensitivities to sorafenib. The IC50 of fostamatinib for most of these cell lines was 2-4 uM ( Figure S5). To confirm the in vivo efficacy of fostamatinib, we used a subcutaneous xenograft model of MHCCLM3 cells that inherently have reduced sensitivity to sorafenib. The growth of tumors was significantly inhibited in fostamatinib-treated mice compared to vehicle-treated mice ( Figure 5A). Comparable body weights of the tumor-bearing mice treated with vehicle or fostamatinib ( Figure 5A) indicate no systemic toxicity of the drug. Tumors harvested at the end of 28 days of treatment weighed an average of 3.2 g in vehicle-treated mice and 1.8 g in fostamatinib-treated mice ( Figure 5A). However, due to the small sample size, this difference was not statistically significant (p = 0.1031).
Based on the IPA analyses of the gene expression profile of fostamatinib-treated HepG2 cells from LINCS ( Figure 4B), we performed immunoblot analysis to validate the kinases that are inhibited by fostamatinib. Our data demonstrate that fostamatinib inhibits phosphorylation of ERK, AKT and STAT3 but has no effect on PTEN. Additionally, we observed inhibitory effect on EGFR and JNK ( Figure 5B).
Finally, we sought to validate predicted gene expression changes associated with fostamatinib treatment from LINCS in our Huh7 cell line model. We performed RNA-seq analysis of fostamatinibvs. DMSO vehicle-treated sorafenib-resistant Huh7-R cells and determined the concordance of fold changes for differentially expressed genes (adjusted p < 0.05) with both the LINCS fostamatinib gene expression signature and Huh7-R sorafenib resistance gene expression signature ( Figure S6). Prioritized differentially expressed genes that were concordant with the LINCS fostamatinib signature and discordant with the Huh7-R sorafenib resistance signature are shown in Table S4. We observed 18 up-regulated and 9 down-regulated genes significantly differentially expressed genes in the Huh7-R signature that were reversed in the fostamatinib treated Huh7 cells and LINCS fostamatinib HepG2 signature. Several examples of genes implicated in HCC tumorigenesis and/or sorafenib resistance that were up-regulated in the Huh7-R cells and down-regulated following fostamatinib treatment in both Huh7-R and HepG2 HCC cells include transforming growth factor beta 2 (TGFB2) [19,20], hypoxia inducible factor 1 alpha subunit (HIF1A) [21], intercellular adhesion molecule 1 (ICAM1) [22,23] and ETS proto-oncogene 1 transcription factor (ETS1) [24].  Figure S7. * indicates a statistically significant difference between the two groups at p = 0.05.

Effect of Etiology with Drug Repurposing Hypothesis
We assessed whether HCC etiology could influence drug repurposing predictions. The GEO HCC patient gene expression datasets, described in Table 1, were filtered to include patient tumors specific to a given HCC etiology: hepatitis B virus (HBV), hepatitis C virus (HCV) and alcohol-induced (AI). Etiologyspecific gene expression signatures were used in connectivity mapping analysis. Using hierarchical clustering analysis of drug connectivity scores, we found that all three HBV patient datasets clustered together, and that two of the three HCV patient datasets clustered together with the one AI patient dataset ( Figure 6A).  Figure S7. * Indicates a statistically significant difference between the two groups at p = 0.05.

Effect of Etiology with Drug Repurposing Hypothesis
We assessed whether HCC etiology could influence drug repurposing predictions. The GEO HCC patient gene expression datasets, described in Table 1, were filtered to include patient tumors specific to a given HCC etiology: hepatitis B virus (HBV), hepatitis C virus (HCV) and alcohol-induced (AI). Etiology-specific gene expression signatures were used in connectivity mapping analysis. Using hierarchical clustering analysis of drug connectivity scores, we found that all three HBV patient datasets clustered together, and that two of the three HCV patient datasets clustered together with the one AI patient dataset ( Figure 6A).

Association of Clinical and Demographic Factors with Sorafenib Resistance Signature
We assessed whether clinical factors (etiology, clinical stage, Child-Pugh class) and patient demography (race, gender) influenced sorafenib resistance gene signature status. We observed that HCC etiology could affect the proportion of tumors harboring SR+ and SR− gene signatures in the TCGA LIHC dataset (chi-square test, p = 0.0279), as shown in Figure 6B. Two etiologies that exhibited significant differences in proportion of SR+/SR− HCC tumors included HBV (Fisher exact test, p = 0.0444) and nonalcoholic fatty liver disease (NAFLD) (Fisher exact text, p = 0.0180). We also examined the proportion of SR+ and SR− HCC tumors among different stages, Child-Pugh class, race and gender ( Figure 6C-F). HCC patient race exhibited a significant trend (chi-square test, p = 0.004), as shown in Figure 6E, where a higher proportion of SR+ and SR-HCC tumors were found in white and Asian patients, respectively. However, no significant trend was observed for clinical stage, Child-Pugh class or gender among SR+ vs. SR− HCC tumor status.

Discussion
In this paper, we address how different publicly available gene expression datasets derived from in vitro and in vivo models of sorafenib resistance in HCC may be reliably assessed in HCC patient tissue in

Association of Clinical and Demographic Factors with Sorafenib Resistance Signature
We assessed whether clinical factors (etiology, clinical stage, Child-Pugh class) and patient demography (race, gender) influenced sorafenib resistance gene signature status. We observed that HCC etiology could affect the proportion of tumors harboring SR+ and SR− gene signatures in the TCGA LIHC dataset (chi-square test, p = 0.0279), as shown in Figure 6B. Two etiologies that exhibited significant differences in proportion of SR+/SR− HCC tumors included HBV (Fisher exact test, p = 0.0444) and non-alcoholic fatty liver disease (NAFLD) (Fisher exact text, p = 0.0180). We also examined the proportion of SR+ and SR− HCC tumors among different stages, Child-Pugh class, race and gender ( Figure 6C-F). HCC patient race exhibited a significant trend (chi-square test, p = 0.004), as shown in Figure 6E, where a higher proportion of SR+ and SR-HCC tumors were found in white and Asian patients, respectively. However, no significant trend was observed for clinical stage, Child-Pugh class or gender among SR+ vs. SR− HCC tumor status.

Discussion
In this paper, we address how different publicly available gene expression datasets derived from in vitro and in vivo models of sorafenib resistance in HCC may be reliably assessed in HCC patient tissue in terms of their prognostic significance and ability to derive drug repurposing predictions. We also provide proof of concept evidence for the successful use of the computational drug repurposing approach in the context of sorafenib-resistant HCC. Similarly, we report on differences in drug repurposing hypotheses with respect to HCC etiology and cell line source. Although connectivity mapping has been applied across diverse domains, this is the first study that has sought to explicitly identify drugs for sorafenib-resistant HCC using drug-induced gene expression signatures from the Library of Integrated Network-based Cellular Signatures (LINCS) database.
Several previous studies have applied HCC gene signatures to the Connectivity Map (CMap) database, which consists of a collection of gene expression profiles of five human cancer (non-HCC) cell lines treated with 1309 compounds [12]. Two studies used HCC gene expression signatures to query against the CMap database, and validated several drug candidates in vitro and in vivo [25,26]. Another group used a combination of CMap and LINCS to discover novel HCC drugs, and validated three anthelminitics in primary hepatocytes and two mouse models [27]. Lv et al. queried CMap using the HCC-3sp-R gene expression data evaluated in this study, and generated six drug predictions to reverse resistance to sorafenib in HCC; however, these predictions were not validated in any context in this publication and no information regarding drug mechanisms was presented [11]. While these previous studies have demonstrated the feasibility of validating connectivity mapping predictions in HCC, all studies to date have been mostly limited to the use of the CMap drug reference database. Two distinct advantages of this LINCS-based study include the ability to gauge the effect of drug perturbations on liver cancer cells (HepG2), as well as the increased number of perturbagens tested in the LINCS systems (n = 3740 total in HepG2 cell line).
In this study, fostmatinib was discovered to be a potentially useful first-line therapy or following resistance to sorafenib, either alone or in combination with sorafenib. Fostamatinib is an inhibitor of spleen tyrosine kinase (SYK), and is currently under investigation for the treatment of several autoimmune diseases [28]. Fostamatinib has been shown to have anti-cancer properties for hematological malignancies [29,30], and this is the first study investigating its use in HCC and sorafenib resistance. SYK is a non-receptor cytoplasmic tyrosine kinase involved in signal transduction in cells of hematopoietic origin, and more recently, implicated both as a tumor suppressor and promoter of cell survival in various hematopoietic and epithelial cancers [31,32]. Reduction in SYK expression has been described as a potential prognostic biomarker in several cancers, including HCC [33][34][35][36]. Although SYK mRNA has prognostic significance in HCC, the lack of its expression at protein level in some HCC cell lines that are sensitive to fostamatinib suggested that the drug functions through other targets that remain to be discovered. We have demonstrated that several kinases (ARK, AKT, STAT3, EGFR, and JNK) are inhibited by fostamatinib but future follow-up studies will be needed to identify direct target(s) of fostamatinib in HCC.
Dasatinib was confirmed to have a unique role in inhibiting cell growth of sorafenib-resistant HCC cells. The Src family kinase inhibitor dasatinib is approved for the treatment of Ph+ chronic myeloid leukemia (CML) in chronic phase and imatinib-resistant disease, Ph+ acute lymphoblastic leukemia with resistance to prior therapy, and is under clinical investigation for solid cancers. Src family kinase activity has been implicated in several oncogenic processes, including cellular proliferation, survival, migration and angiogenesis, and increased activity has been demonstrated in HCC in vitro [37][38][39]. Our experiments demonstrated that dasatinib, a Src family kinase inhibitor, was effective in reducing HCC cell viability and colony formation alone and in combination with sorafenib in sorafenib-resistant HCC cells. Additionally, we showed that Src family kinases were significantly activated in the sorafenib-resistant HCC cells as compared to sorafenib-sensitive HCC cells, consistent with the known mechanism of action of dasatinib.
Recently, dasatinib was shown to be successful in reducing HCC cell proliferation, adhesion, migration and invasion in vitro via inhibiting Src and several downstream signaling pathways, including PI3K/PTEN/Akt and SFK/FAK [40]. Another study found that phosphorylation of Src was inhibited in a panel of HCC cells that were sensitive and resistant to dasatinib, and that cell proliferation was not affected by knocking down Src and p-Src in dasatinib-sensitive cells. The authors concluded that dasatinib-mediated inhibition of Src alone is not sufficient to induce its anti-proliferative or pro-apoptotic effects, and that dasatinib may mediate its effects via other targets in addition to Src [41]. Finally, dasatinib was tested in patients with advanced HCC in a recent phase II clinical trial (NCT00459108), but was terminated early due to futility. The primary objectives were to determine the progression-free survival (PFS) rate and response rate at 4 months in patients with unresectable advanced HCC treated with dasatinib. Several factors that may have influenced the results of this clinical trial include compromised liver status in advanced HCC patients and the use of RECIST criteria to determine response rate, which is known to be ineffective in evaluating cytostatic agents, including sorafenib [42]. Furthermore, our results suggest that dasatinib may be most useful for an enriched HCC patient population with transcriptomic biomarkers characteristic of sorafenib resistance.
We also found that HCC etiology may influence sorafenib resistance and drug repurposing hypothesis generation using the transcriptomics-based LINCS system. Interestingly, sorafenib was previously observed to be more effective for HCC patients with an underlying HCV infection compared to HBV infection or alcoholic cirrhosis [43]. In another study, dasatinib was shown to be most effective in a group of HCC patients with a "progenitor molecular subtype", as assessed by gene expression profiling [41]. Taken together, these findings highlight the importance of considering HCC patient etiology and other molecular features in drug prediction studies.

Cell Culture and Acquired Sorafenib Resistance
All cells were maintained in minimum essential media supplemented with L-glutamine (2 mM), 10% FBS, sodium pyruvate (0.11 g/L) and penicillin/streptomycin (100 U/mL). Cell media for sorafenib resistant cell lines were also supplemented with sorafenib (6 µM and 0.1% DMSO). Sorafenib was withdrawn from the cell media of resistant Huh7 cells for 5-7 days prior to performing all experiments. The HCC cell line Huh7 cells were generously provided by Dr. James Taylor (Fox Chase Center, PA, USA). Sorafenib resistant cells (Huh7-R) and several resistant clones ( Figure S1) were generated as described previously [14]. Huh-R-A7 was randomly selected to represent sorafenib resistant clones for all subsequent experiments.

Cell Viability Assay
Cells were seeded into 96-well plates (~2000 cells/well) and allowed to incubate overnight. The next day, cell media were replaced with MEM media containing specified concentration of sorafenib, fostamatinib or dasatinib (with 0.1% final DMSO concentration). After 48 h of treatment, the CellTiter-Glo ® luminescent viability assay was utilized following the manufacturer's instructions. Prior to measuring viability with a luminometer, the luminescent supernatant was transferred to an opaque luminometer 96-well plate.

Colony Formation Assay
Cells were seeded into 6-well plates (~2000-5000 cells/well) and allowed to incubate for 24-48 h. Cells were then treated with a continuous dose of therapeutics (with 0.1% final DMSO concentration) for 14-18 days. Media were replaced every 2-3 days. After colonies grew to a sufficient size, cells were fixed with 3.7% paraformaldehyde (in PBS) and stained with a 0.05% crystal violet solution.

Mouse Strains, Animal Husbandry and Treatment
Male NSG (NOD scid gamma) mice were purchased from Target Validation Shared Resource (TVSR) core facility at the Ohio State University. All animals were housed in a temperature-controlled room under a 12 h light/12 h dark cycle and under helicobacter-free conditions and fed normal chow diet. All animal studies were reviewed and approved by the Ohio State University Institutional Laboratory Animal Care and Use Committee (Protocol # 2008A0236). For drug treatment, fostamatinib disodium (Cat# DC1013, DC-Chemicals, Shanghai, China) was dissolved with 30% polyethylene glycol 400 (PEG 400) (Cat# 91893, Sigma-Aldrich, St. Louis, MO, USA), 5% propylene glycol (PPG) (Cat# P4347, Sigma-Aldrich), and 0.5% Tween 80 (Cat# P4780, Sigma-Aldrich) to the final concentration of 20 mg/mL immediately before delivery through oral gavage. For subcutaneous xenografts, 10-12 weeks old NSG mice were injected subcutaneously with HCCLM3 (2.5 × 10 6 ) cells into the right flank. When tumor volume reached 100 mm 3 , mice were randomized into 2 groups for treatment with either vehicle (30% PEG 400, 5% PPG and 0.5% Tween 80) or fostamatinib disodium (80 mg/kg) administrated daily through oral gavage. Tumor volumes based on digital caliper measurements were calculated by the ellipsoidal formula (1/2(length × width 2 )). After 28 days of treatment, mice were euthanized and tumor tissues were collected, weighed and photographed.

Human Phospho-Protein Array
Cell lysates of Huh7 parental and Huh7 sorafenib-resistant pool cells were subjected to phosphoprotein analysis using the Proteome Profiler Human Phospho-Kinase Array Kit (#ARY003B, R & D Systems, Minneapolis, MN, USA) following the manufacturer's protocol to quantify phosphorylation levels of 43 proteins phosphorylated at tyrosine, serine, or threonine residues.

Sorafenib Resistance Gene Expression Signatures
Microarray analysis was performed on the parental cells (Huh7-S), a pool of sorafenib-resistant cells (Huh7-R) and sorafenib-resistant clone A7 (Huh7-R-A7) using the Affymetrix GeneChip Human Transcriptome Array 2.0 platform. Three biological replicates of each cell type were used for the microarray analysis. A differential gene expression signature was defined by comparing microarray data from Huh7-R-A7 vs. Huh7-S cells. Signal intensities were analyzed by Affymetrix Expression Console software. Gene expression levels were RMA-normalized and log-transformed [44]. A filtering method based on the percentage of arrays (85%) below a noise cutoff of 6 (log2 scale) was applied to filter out low expression genes, and a linear model was employed to detect differentially expressed genes. In order to improve the estimates of variability and statistical tests for differential expression, a variance smoothing method with fully moderated t-statistic was employed for this study [45]. The significance level was adjusted for multiple hypothesis testing by controlling the mean number of false positives, and a threshold p-value < 0.0001 was maintained to determine statistical significance [46]. We averaged gene expression values for multiple probeset ID's mapping to the same gene. A fold-change cutoff of >3 for up-regulated genes and <0.25 for down-regulated genes was imposed. Raw and normalized data were deposited in the Gene Expression Omnibus (GEO) database (accession: GSE94550).

HCC Tumor Gene Expression Data
A systematic search of the Gene Expression Omnibus (GEO) database was conducted for datasets containing human HCC and normal liver tissue. We obtained raw (CEL files) gene expression data from six GEO microarray datasets: GSE14323 (GPL571), GSE14520 (GPL571), GSE14520 (GPL3921), GSE45267 (GPL570), GSE62232 (GPL570) and GSE6764 (GPL570). Signal intensity values were RMA-normalized and gene expression values were log-transformed. Differential gene expression analysis (tumor vs. normal) was conducted via the limma (Linear Models for Microarray Analysis) R package [47], which employs an empirical Bayes method to moderate the standard errors of estimated log-fold changes. We averaged gene expression values for multiple probeset ID's mapping to the same gene. Benjamini-Hochberg False Discovery Rate (FDR) correction was applied to adjust for multiple hypotheses testing, and significance cutoff was set at adjusted p < 0.0001 [48]. Gene overlap Venn diagrams were generated using Venny 2.1.0 (available at https://bioinfogp.cnb.csic.es/tools/venny/ index.html).

RNA-Seq Analysis of Fostamatinib-Treated Huh7-R cells
RNA-Seq reads of DMSO-or fostamatinib-treated Huh7-R cells (n = 3 each) were first mapped to the human genome hg19 using Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT) [49]. Raw read counts for each gene were quantified by using the featureCounts program [50]. Then RNA-seq counts were obtained by using GENCODE v.22 Gene Transfer Format (GTF) file as a transcript reference (GENCODE annotation). Genes with read counts below 5 for at least 2 samples out of 3 within each group were filtered out. Then the read counts were normalized with the TMM method [51]. To identify genes differentially expressed between samples, the limma R package was used to calculate p-values for group comparisons under a linear model [47]. The p-value cut-offs were determined by controlling the mean number of false positives [52]. Raw and normalized data were deposited in the Gene Expression Omnibus (GEO) database (accession: GSE113005).

Library of Integrated Network-Based Cellular Signatures (LINCS) Analyses
Connectivity mapping analyses were conducted via the Library of Integrated Network-based Cellular Signatures (LINCS) system (database version A2) using the web-based platform (http://www.lincscloud.org/). Connectivity scores were calculated using the weighted Kolgomorov-Smirnov (KS) statistic to rank predictions from the LINCS database, as previously described [12]. We selected LINCS compound perturbations tested exclusively in the HepG2 liver cancer cell line. Connectivity scores were averaged for individual LINCS compound perturbations tested at different concentrations and time points in the HepG2 cell line. Individual LINCS gene signatures were obtained from the Broad LINCS Cmap C3 Cloud Compute platform using the slice_slice_tool command.

Database Access: DrugBank, AACT, KEGG, Cell Miner, Beegle
We downloaded the "Approved" and "Investigational" external drug link files from the DrugBank database (v 4.5.0, released date 20 April 2016, available at https://go.drugbank.com/releases/4-5-0) and conducted searches in the DrugBank web interface to categorize drug approval status [53]. We downloaded the Aggregate Analysis of ClinicalTrials.gov (AACT) database, a publicly available resource produced by the Clinical Trials Transformation Initiative (CTTI), for drugs under clinical investigation (AACT accessed: 17 June 2016). We used KEGG DRUG database ID's mapped to drugs in DrugBank to extract drug activity and target pathway information in KEGG Drug (KEGG DRUG accessed: 20 June 2016) [54]. Drug target genes were obtained from DrugBank and KEGG Drug databases, and assessed for association with HCC using the Beegle literature-mining tool [18]. Gene expression signatures for 18 HCC cell lines vs. a pool of 19 normal liver samples were obtained via the CellMinerHCC database [55].

Bioinformatics Analyses
Hierarchical clustering using the Euclidean distance of gene expression and drug connectivity scores was performed using the heatmap.2 function from the "gplots" R package. The nearest-template prediction (NTP) method was applied to normalized, log-transformed gene expression data to classify tumor samples as SR+ or SR− (FDR < 0.05) [56]. Gene mutation, copy-number and mRNA expression data for SR+/SR− liver HCC (LIHC) patients in the TCGA analysis were obtained via the cBioPortal (v 1.2.4) tool [57]. The drug target protein-protein interaction (PPI) network was generated using the STRING (v 10.0) database [58]. Only high confidence interaction scores (0.700 and above) from experiments, databases, neighborhood and gene fusion sources were included in the final network. Network analysis of the drug target PPI network was performed using Gephi 0.9.1 software [59], including algorithms to detect the following network features: degree, modularity class, eigencentrality and clustering coefficient. Enriched gene ontology functions in the sorafenib resistance gene signature were obtained via functional enrichment analysis using the ToppGene Suite [60]. Fostamatinib-treated HepG2 gene expression data from LINCS was analyzed via QIAGEN's Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, CA, USA, www.qiagen.com/ingenuity).

Survival Analysis
Clinical data and RNAseq data (Level 3, v2, RSEM-normalized) from 377 liver HCC patients contained in the TCGA database were obtained via the Broad Institute Firebrowse tool (http:// firebrowse.org/; TCGA data version 2016_01_28). Survival analysis comparing SR+ and SR− patients was conducted using Prism 7 software. Survival analysis of HCC patients for SYK based on gene expression data was conducted using the cBioPortal (v 1.2.4) tool for TCGA data [57]. Death was selected as the survival measure, and the median was chosen as the bifurcation point to define "high" vs. "low" gene expression. Survival curves were generated using the Kaplan-Meier method, and the log-rank test was used for statistical comparison.

Ethics Information for Publicly Available Datasets
Patient gene expression data were obtained from publicly available databases, including Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA). We confirm that the publicly available data were collected with institutional approval and informed consent, and this information can be found in the original study publications: GSE14323 [61], GSE14520 [62], GSE45267 [63], GSE62232 [64], GSE6764 [65] and GSE26391 [16]. Data collection policies for the TCGA can be found at the following website: https://cancergenome.nih.gov/abouttcga/policies.

Conclusions
We show the feasibility of the drug repurposing workflow by validating novel drugs to use alone and in combination with sorafenib in HCC. Our analysis of publicly available gene expression datasets of sorafenib resistance models and HCC patient datasets to determine prioritized drug candidates may inform future additional validation studies. Future studies utilizing data from HCC patient samples with resistant disease will be needed as it becomes available. Importantly, future studies of novel HCC drugs must carefully consider safety and toxicity profiles, as one of the greatest barriers to clinical trial success for approval in this patient population is liver function [4].
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-6694/12/10/2730/s1, Figure S1: Sorafenib dose response curve for sorafenib sensitive and resistant Huh7 cells; Figure S2: Comparison of HCC sorafenib resistance gene signatures; Figure S3: Drug target network. Resultant drug target protein-protein interaction (PPI) network from prioritized drug candidates using PPI interaction data from the STRING database (v 10.0); Figure S4: Dot blots of Proteome Profiler Human Phospho-Kinase Array; Figure S5: Fostamatinib dose response curve for HCC cell lines; Figure S6: Correlation plot of gene expression log2 fold change between Huh7 sorafenib resistance (Huh7-R) signature and fostamatinib-treated Huh7 cells; Figure S7: Uncropped Western blot images; Table S1: Gene expression signature of sorafenib resistance in HCC tumors cells; Table S2: LINCS drug predictions to reverse A7 clone sorafenib resistance gene signature (Huh7-R-A7 vs. Huh7-S); Table S3: Drug target gene PPI network alterations in SR+ and SR-HCC primary tumors for patients in the TCGA LIHC dataset; Table S4: Concordance of fostamatinib-induced gene expression signatures with HCC sorafenib resistance gene expression signature. Funding: This work was supported in part by NIH National Library of Medicine T15LM0112750 to KRF and TM, NIH National Cancer Institute R01 CA086978 to STJ, Pelotonia Idea Grant to STJ, and NIH National Cancer Institute R01 CA193244 to KG.