Bioinformatics Data Mining Repurposes the JAK2 (Janus Kinase 2) Inhibitor Fedratinib for Treating Pancreatic Ductal Adenocarcinoma by Reversing the KRAS (Kirsten Rat Sarcoma 2 Viral Oncogene Homolog)-Driven Gene Signature

Pancreatic ductal adenocarcinoma (PDAC) is still one of the most aggressive and lethal cancer types due to the late diagnosis, high metastatic potential, and drug resistance. The development of novel therapeutic strategies is urgently needed. KRAS (Kirsten rat sarcoma 2 viral oncogene homolog) is the major driver mutation gene for PDAC tumorigenesis. In this study, we mined cancer genomics data and identified a common KRAS-driven gene signature in PDAC, which is related to cell–cell and cell–extracellular matrix (ECM) interactions. Higher expression of this gene signature was associated with poorer overall survival of PDAC patients. Connectivity Map (CMap) analysis and drug sensitivity profiling predicted that a clinically approved JAK2 (Janus kinase 2)-selective inhibitor, fedratinib (also known as TG-101348), could reverse the KRAS-driven gene signature and exhibit KRAS-dependent anticancer activity in PDAC cells. As an approved treatment for myelofibrosis, the pharmacological and toxicological profiles of fedratinib have been well characterized. It may be repurposed for treating KRAS-driven PDAC in the future.


Introduction
The occurrence of pancreatic cancer has significantly ascended throughout the past decade. Among them, pancreatic ductal adenocarcinoma (PDAC) accounts for most cases of pancreatic cancer, with its survival rate being lower than 8% [1,2]. The death rate is highly correlated with high incidence of metastasis, recurrence rate, and chemoresistance. Due to late diagnosis of most clinical cases, the aggressive type was often accompanied by angiogenesis and metastasis, resulting in high unresectable clinical

Cancer Genomics Analysis via the cBioPortal Website
The cBioPortal (http://www.cbioportal.org/) is a website to access, analyze, and visualize the large-scale TCGA (The Cancer Genome Atlas) cancer genomics data sets or other studies [36,37]. The "Pancreatic adenocarcinoma (TCGA, PanCancer Atlas)" dataset of 168 PDAC patients containing complete genetic status (mutation, copy number variation, and mRNA expression) was used in this study to compare the association between gene mutations and PDAC gene signature. In addition, a Kaplan-Meier survival plot was generated using the cBioPortal to investigate the impact of PDAC gene signature on patients' overall survival.

Connectivity Map Analysis
The Connectivity Map (CMap; https://clue.io/) database contains numerous gene signatures from cultured human cancer cell lines treated with drugs [38]. It is believed that a drug has the potential for treating a disease if this drug could reverse the disease-associated gene signature [39,40]. To identify the potential drugs to reverse PDAC gene signature, the commonly upregulated 53 genes were inputted to query the CMap database. The results were visualized as a heat map with a connectivity score between −100 and 100 corresponding to the magnitude of dissimilarity and similarity between queried and existing gene signatures.

Drug Sensitivity Profiling in Pancreatic Ductal Adenocarcinoma Cell Lines
The correlations between KRAS gene expression and drug sensitivity in PDAC cancer cell lines were obtained from the CellMinerCDB (https://discover.nci.nih.gov/cellminercdb/ [41]). The Cancer Therapeutics Response Portal (CTRP [42][43][44]) data from the Broad Institute of Massachusetts Institute of Technology (MIT) and Harvard (Cambridge, MA, USA) were used.

Identification of a Common Gene Signature in Human Pancreatic Ductal Adenocarcinoma
To identify the common gene signature associated with PDAC, five microarray data sets (Table 1) were obtained from the NCBI-GEO database [21]. Then, the DEGs were prepared using the R-based web application, GEO2R [21]. These DEGs (Supplementary Material, File S1) were analyzed using the InteractiVenn web-based tool [22]. As shown in the Venn diagrams ( Figure 1A), we identified 53 upregulated and 2 downregulated genes that were common in PDAC tissues when compared with the adjacent normal tissues ( Figure 1A). Their expression levels are listed in Table 2 and visualized in a heat map ( Figure 1B). To investigate the potential role of this common gene signature, pathway enrichment for the 53 upregulated genes was performed using the WebGestalt web-based tool [23] against GO biological processes [25,26], KEGG pathways [27], and cancer hallmarks [28]. We found that pathways related to cell-cell and cell-ECM interactions were significantly enriched, such as KEGG_ECM-receptor interaction, KEGG_Focal adhesion, HALLMARK_APICAL_JUNCTION, GO_Cell junction organization, GO_Integrin-mediated signaling pathway, and GO_Extracellular structure organization ( Figure 2A). The network for the 53 upregulated genes was further constructed and functional enrichment was performed for GO biological processes and KEGG pathways using the STRING database [24]. As shown in Figure 2B, ITGA2, ITGB4, LAMA3, LAMC2, LAMB3, and GPRC5A genes formed a major cluster, which participated in ECM-receptor interaction, focal adhesion, cell junction organization (together with CDH3 and ECT2 genes), and extracellular organization (together with SERPINB5 and MMP11 genes). Therefore, the alteration of genes related to cell-cell and cell-ECM interactions is a common gene signature in PDAC.     : a gradient color key shows the overlapped gene numbers in a pathway. In the volcano plot of (A), the purple circles for HALLMARK_ESTROGEN_RESPONSE_LATE/HALLMARK_KRAS_ SIGNALING_UP or KEGG_Focal adhesion/HALLMARK_MITOTIC_SPINDLE were overlapped. In the left part of (B), line colors indicate the types of interaction evidence. The cyan and pink lines indicate protein-protein interactions from curated and experimental data, respectively. The purple line indicates that two protein molecules share structural homology. Functional enrichment (gene ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways) in this network are shown in the right part of (B). Selected GO biological processes and KEGG pathways are highlighted with different colors. The term "count in gene set" indicates the overlapped genes (the first number) in a pathway (the second number). The term "false discovery rate" indicates the average rate of false coverage for the functional enrichment.

The Pancreatic Ductal Adenocarcinoma Gene Signature Was Associated with KRAS and TP53 Gene Mutations
To further confirm the role of PDAC gene signature, the TCGA-PAAD (pancreatic adenocarcinoma) data set with 168 PDAC cases was used to compare their mRNA levels. As shown : a gradient color key shows the overlapped gene numbers in a pathway. In the volcano plot of (A), the purple circles for HALLMARK_ESTROGEN_RESPONSE_LATE/HALLMARK_KRAS_ SIGNALING_UP or KEGG_Focal adhesion/HALLMARK_MITOTIC_SPINDLE were overlapped. In the left part of (B), line colors indicate the types of interaction evidence. The cyan and pink lines indicate protein-protein interactions from curated and experimental data, respectively. The purple line indicates that two protein molecules share structural homology. Functional enrichment (gene ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways) in this network are shown in the right part of (B). Selected GO biological processes and KEGG pathways are highlighted with different colors. The term "count in gene set" indicates the overlapped genes (the first number) in a pathway (the second number). The term "false discovery rate" indicates the average rate of false coverage for the functional enrichment.

The Pancreatic Ductal Adenocarcinoma Gene Signature Was Associated with KRAS and TP53 Gene Mutations
To further confirm the role of PDAC gene signature, the TCGA-PAAD (pancreatic adenocarcinoma) data set with 168 PDAC cases was used to compare their mRNA levels. As shown in Figure 3, most of them have higher mRNA expressions, especially for the 10 genes related to cell-cell and cell-ECM interactions: LAMB3 (10%), SERPINB5 (10%), CDH3 (8%), ECT2 (8%), LAMC2 (8%), ITGB4 (6%), ITGA2 (5%), LAMA3 (5%), GPRC5A (4%), and MMP11 (4%). PDAC patients with the higher gene signature (53 upregulated genes only) expression have poorer overall survival ( Figure 4A). In addition, we found that the upregulation of PDAC gene signature was significantly associated with KRAS and TP53 gene mutations ( Figure 4B,C). According to the TCGA-PAAD data set, 110 (65%) and 102 (61%) of the 168 PDAC cases harbored KRAS and TP53 gene mutations, respectively ( Figure 4B). Given the fact that KRAS and TP53 were the most two common mutated genes in PDAC [6,7], it is reasonable that the PDAC gene signature may be driven by KRAS and TP53 gene mutations during tumorigenesis. in Figure 3, most of them have higher mRNA expressions, especially for the 10 genes related to cellcell and cell-ECM interactions: , and MMP11 (4%). PDAC patients with the higher gene signature (53 upregulated genes only) expression have poorer overall survival ( Figure  4A). In addition, we found that the upregulation of PDAC gene signature was significantly associated with KRAS and TP53 gene mutations ( Figure 4B,C). According to the TCGA-PAAD data set, 110 (65%) and 102 (61%) of the 168 PDAC cases harbored KRAS and TP53 gene mutations, respectively ( Figure  4B). Given the fact that KRAS and TP53 were the most two common mutated genes in PDAC [6,7], it is reasonable that the PDAC gene signature may be driven by KRAS and TP53 gene mutations during tumorigenesis. The cases highlighted in red grids (labeled as "mRNA High") indicate those with mRNA expression higher than that of the average patient (z-score >2).

Figure 3.
A waterfall plot for the common gene signature expression in "The Cancer Genome Atlas-Pancreatic adenocarcinoma" data set. Genes highlighted in red or blue color indicate those commonly upregulated or downregulated in pancreatic ductal adenocarcinoma (PDAC) patients, respectively. The cases highlighted in red grids (labeled as "mRNA High") indicate those with mRNA expression higher than that of the average patient (z-score >2).

Gene Set Enrichment Analysis Revealed That the Pancreatic Ductal Adenocarcinoma Gene Signature Was Driven by KRAS Gene Mutation
To further investigate whether the PDAC gene signature was driven by KRAS and TP53 gene mutations, the effect of KRAS G12D (glycine 12 to aspartate) or TP53 R175H (arginine 175 to histidine; the human equivalent of mouse Trp53 R172H ) mutations on PDAC gene signature was analyzed by GSEA using the relevant microarray data sets (Table 3). In GSE58055 [30], a doxycycline (Dox)-inducible KRAS G12D mutation was introduced into the E6/E7-transformed human pancreatic ductal epithelial (HPDE) cells, in which the E6 and E7 proteins of the HPV16 virus inactivate p53 and RB, respectively [45,46]. We found that the PDAC gene signature was only enriched in HPDE cells with KRAS G12D induction by Dox, but not in cells with the induction of wild-type (WT) KRAS or green fluorescent protein (GFP) control vector ( Figure 5A), suggesting that the PDAC gene signature can be driven by KRAS G12D mutation. To confirm the above observation, another microarray data set GSE53659 [31] with the Kras G12D -driven PDAC in mice (Pdx1-Cre/Kras G12D/+ ; also known as KC mice) was used. As shown in Figure 5B, the PDAC gene signature was enriched in KRAS G12D -driven PDAC cells compared with that in the WT cells. Therefore, the PDAC gene signature is indeed driven by KRAS G12D mutation.
To investigate the role of TP53 (Trp53 in mice) gene mutation, two microarray data sets, GSE67358 [32] and GSE123646 [33], were employed. It has been shown that one-third of KC mice develop PDAC by 500 days [47], and the additional Trp53 mutation in KPC (Pdx1-Cre/Kras G12D/+ /Trp53 R172H/+ ) mice or Trp53 deletion in KP fl C (Pdx1-Cre/Kras G12D/+ /Trp53 -/+ ) mice accelerates the tumor development by 120-180 days [48]. However, only Trp53 mutation, but not deletion, can drive tumor metastasis in this model [49], suggesting a synergy between KRAS and TP53 mutations to promote PDAC progression. Because the incidence of tumor metastasis is about 65% in KPC mice [49], the gene expression profiles of both metastatic (meta) and non-metastatic (no meta) PDAC cells from KPC mice were compared with that from KP fl C mice ( Figure 5C, the left and middle parts). In addition, the gene expression profiles of metastatic and non-metastatic PDAC cells were also compared to each other ( Figure 5C, the right part). We found that the PDAC gene signature was not enriched in any group, suggesting that the PDAC gene signature was not associated with TP53 (Trp53) mutation and its metastasis-promoting effect. It was puzzling that inconsistent observation was found in the GSE123646 data set ( Figure 5D). The PDAC gene signature was enriched in KPC mice-derived PDAC cells (irrespective of their metastatic status) compared with that in KP fl C mice-derived PDAC cells. However, the PDAC gene signature was not enriched in KP fl C mice-derived PDAC cells transfected with the human equivalent of murine Trp53 R172H (TP53 R175H ). Such discrepancy may imply the minimal effect of TP53 gene mutation on PDAC gene signature expression, which warrants further investigation. Table 3. Microarray data sets from KRAS-mutant and TP53 (Trp53)-mutant cells.  The above results argue for the essential role of KRAS gene mutation in PDAC gene signature expression. To further investigate the expression of PDAC gene signature during KRAS G12D -driven pancreatic tumorigenesis, the gene expression profiles of normal pancreas, pancreatic intraepithelial neoplasia (PanIN) and PDAC in KC mice were obtained from the microarray data set GSE33323 [29]. GSEA showed that the PDAC gene signature is significantly correlated with PanIN and PDAC ( Figure 5E). The related expression of PDAC gene signature was visualized in a heat map showing that the PDAC gene signature was induced during KRAS G12D -driven PDAC development in KC mice ( Figure 5F). Taken together, we conclude that the PDAC gene signature is driven by KRAS, but not TP53, gene mutation.

Connectivity Map Analysis and Drug Sensitivity Profiling Identify TG-101348 (Fedratinib) as a Potential Drug Reversing KRAS-Driven Pancreatic Ductal Adenocarcinoma Gene Signature
To identify potential drugs that could reverse the KRAS-driven PDAC gene signature, we employed the CMap database that contains numerous gene signatures from cultured human cancer cell lines treated with drugs [38]. If a drug could reverse a disease-associated gene signature, this drug is believed to have the potential to cure the disease [39,40]. We queried the CMap database with the PDAC gene signature (53 upregulated genes) to connect the PDAC gene signature to drug-derived gene signatures. The CMap connectivity score ranging from −100 to 100 corresponds to the magnitude of dissimilarity and similarity between queried and existing gene signatures. Figure 6A showed the most dissimilar drugs (connectivity score < −95) representing the potential drugs that could reverse the queried PDAC gene signature. Interestingly, most of them belong to histone deacetylase (HDAC) inhibitors including trichostatin A (pan-HDAC), panobinostat (pan-HDAC), ISOX (HDAC6-specific), apicidin (pan-HDAC), and vorinostat (pan-HDAC). Therefore, inhibition of HDAC might have the potential to treat PDAC by reversing its KRAS-driven gene signature.
We hypothesized that a drug may exhibit KRAS-dependent cytotoxicity in cancer cells if this drug could reverse the KRAS-driven gene signature. To examine whether the predicted CMap drugs could exhibit KRAS-dependent cytotoxicity in PDAC cells, we employed the CellMinerCDB database that is a web-based tool enabling to explore and analyze pharmacological and genomic data of human cancer cell lines [41]. Due to the frequent KRAS mutation in PDAC, there was only one PDAC cell line (BxPC-3) harboring the wild-type KRAS gene ( Figure 6B). Thus, it is impossible to correlate the drug activity to KRAS gene mutation. According to the gene expression profiles from PDAC cell lines (CTRP-PAAD) and patients' tissues (TCGA-PAAD; Figure 6C), the mutant KRAS gene tended to be expressed higher compared to the wild-type KRAS gene. Thus, as an alternative, we correlated the drug activity with KRAS mRNA expression. The CTRP database only contained the drug sensitivity profiles for TG-101348, etoposide, ISOX, panobinostat, apicidin, and vorinostat. Surprisingly, pan-HDAC (panobinostat, apicidin, and vorinostat) and HDAC6 (ISOX) inhibitors, as well as etoposide, did not exhibit KRAS-dependent cytotoxicity ( Figure 6D). Only TG-101348 displayed significant association with KRAS expression ( Figure 6D). Therefore, TG-101348 may exhibit KRAS-dependent anticancer activity in PDAC cells via the reversion of KRAS-driven gene signature.
did not exhibit KRAS-dependent cytotoxicity ( Figure 6D). Only TG-101348 displayed significant association with KRAS expression ( Figure 6D). Therefore, TG-101348 may exhibit KRAS-dependent anticancer activity in PDAC cells via the reversion of KRAS-driven gene signature.

Discussion
Dynamic cell-cell and cell-ECM interactions maintain a tumor microenvironment that consists of acellular fibrous stroma and diverse populations of the non-neoplastic cancer-associated cells. Previous studies suggested that the tumor progression of PDAC as well as its deadly malignancy are

Discussion
Dynamic cell-cell and cell-ECM interactions maintain a tumor microenvironment that consists of acellular fibrous stroma and diverse populations of the non-neoplastic cancer-associated cells. Previous studies suggested that the tumor progression of PDAC as well as its deadly malignancy are highly associated with the tumor microenvironment. Thus, targeting the stromal compartment in PDAC may have anticancer effects and enhance chemo-/radio-sensitivity [50][51][52]. Our results imply that inhibition of PDAC gene signature (genes related to cell-cell and cell-ECM interactions) may be beneficial for treating PDAC via a remodeling of the tumor microenvironment.
TG-101348, also known as fedratinib, is a JAK2-selective inhibitor that has been approved for treating patients with myelofibrosis [53]. Myelofibrosis is a rare type of bone marrow cancer, which disrupts the normal production of blood cells. The discovery of JAK2 V617F (valine 617 to phenylalanine) mutation in myelofibrosis uncovers the activated JAK-STAT (signal transducer and activator of transcription) signaling as a primary driver for myelofibrosis, and supports the rationale for treating myelofibrosis by JAK2 inhibition [54]. Interestingly, a previous study has shown that STAT3 plays a critical role in KRAS-induced PDAC tumorigenesis. A large-scale cancer cell line screening identified a JAK2-selective inhibitor, AZ960, that blocks STAT3 activation and exhibits higher sensitivity against PDAC cell lines [55], which supports the utility of therapeutic targeting of JAK2-STAT3 signaling in PDAC.
HDAC inhibitors have been viewed as a prominent class of therapeutic agents for treating PDAC [56,57]. However, the impact of KRAS mutation on their anticancer activity was largely unclear. Our results suggest that the anticancer activity of HDAC inhibitors is unrelated to KRAS mutation status. Consistently, previous studies found that pan-HDAC inhibitors (vorinostat and AR-42) exhibit similar cytotoxicity in both KRAS WT and mutant PDAC cells [58,59]. However, it was also reported that a HDAC inhibitor, romidepsin, preferentially induces apoptosis in cancer cells harboring mutant KRAS [60]. More investigations are needed to clarify the exact role of KRAS mutation in the anticancer activity of HDAC inhibitors.

Conclusions
This study integrates bioinformatics resources to investigate the key driver mutation gene, KRAS, and the associated gene signature in PDAC. Our results demonstrate that the progression and prognosis of PDAC is highly associated with a KRAS-driven gene signature related to cell-cell and cell-ECM interactions. A FDA-approved JAK2-selective inhibitor, fedratinib (TG-101348), is predicted to reverse the KRAS-driven gene signature, thereby providing therapeutic benefit for KRAS-mutated PDAC patients.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2075-4426/10/3/130/s1, Table S1: The correlation between KRAS mutation types and PDAC gene signature, Figure S1: Role of KRAS mutation types in PDAC gene signature expression, Figure S1A: Visualization of the KRAS mutation burden and hotspots in 168 PDAC patients, Figure S1B: A heat map shows the correlation between KRAS mutation types and PDAC gene signature expression, Figure S1C: GSEA results for the role of KRAS mutation types in regulating PDAC gene signature, File S1: The differentially expressed genes (DEGs) prepared from microarray data sets.