Integrated use of bioinformatic resources reveals that co-targeting of histone deacetylases, IKBK and SRC inhibits epithelial-mesenchymal transition in cancer

Abstract

With the advent of high-throughput technologies leading to big data generation, increasing number of gene signatures are being published to predict various features of diseases such as prognosis and patient survival. However, to use these signatures for identifying therapeutic targets, use of additional bioinformatic tools is indispensible part of research. Here, we have generated a pipeline comprised of nearly 15 bioinformatic tools and enrichment statistical methods to propose and validate a drug combination strategy from already approved drugs and present our approach using published pan-cancer epithelial–mesenchymal transition (EMT) signatures as a case study. We observed that histone deacetylases were critical targets to tune expression of multiple epithelial versus mesenchymal genes. Moreover, SRC and IKBK were the principal intracellular kinases regulating multiple signaling pathways. To confirm the anti-EMT efficacy of the proposed target combination in silico, we validated expression of targets in mesenchymal versus epithelial subtypes of ovarian cancer. Additionally, we inhibited the pinpointed proteins in vitro using an invasive lung cancer cell line. We found that whereas low-dose mono-therapy failed to limit cell dispersion from collagen spheroids in a microfluidic device as a metric of EMT, the combination fully inhibited dissociation and invasion of cancer cells toward cocultured endothelial cells. Given the approval status and safety profiles of the suggested drugs, the proposed combination set can be considered in clinical trials.

cancer, combinatorial drug therapy, epithelial–mesenchymal transition (EMT), bioinformatics, systems pharmacology

Introduction

Generation of biological data has been biased and disconnected since the emergence of cellular biology approaches to identify the building blocks of cells and diseases [1]. This approach was expectedly reflected in pharmaceutical industry, which led to production of single-target drugs as magic bullets to inhibit a selected protein to cure complex diseases [2]. With the advent of high-throughput technologies, unbiased data generation is beginning to change the attitude of biologists to systems-level understanding of biological states and therapies [3].

The results of high-throughput experiments are usually conveyed as a list of genes that are differentially altered in that specific condition. From the biological point of view, if the genes within the list are validated to show that they can specifically and collectively discriminate between the conditions such as disease diagnosis, prognosis of patients and response to treatments, they can be regarded as a ‘gene signature’ [4]. Gene signatures can also be purposefully used to identify novel therapeutic target combinations for that disease, however; this approach may not always be so straightforward. First, correlation of expression of a set of genes observed in a disease does not necessarily show them to be ‘causative’ for that disease. Second, gene signatures are derived by different statistical methods such as regression models, neuronal networks and machine learning approaches [5]. So, they usually have redundancy and lack the required overlap [6], and such variations complicate selection of appropriate therapeutic targets from a single signature [7]. These considerations necessitate applying other complementary tools to identify therapeutic targets which converge interdisciplinary methods in statistics, mathematics and computation in shape of bioinformatic approaches that can be further translated into effective holistic treatments [8, 9].

With such outlook, we are facing significant advances in various bioinformatic resources such as repository databases, gene set libraries as well as analytical pipelines such as biologist-friendly enrichment tools and networks [10]. To manage these tools, efforts have been made to categorize these resources in repositories such as OmicTools [11], but the rational integration of these tools to answer specific biological questions is still lacking, which enables biological researchers to perceive principal mechanisms of complex processes as a whole [12, 13].

Here, we aim to show how integration of gene set-based analytical tools with network science can generate a predictive pipeline to identify combination of therapeutic targets from gene signatures. To show the relevance of our approach, we present a practical case study and apply selected gene signatures illustrating epithelial–mesenchymal transition (EMT), which is associated with aggressive behaviors in cancer such as metastasis, drug resistance and stemness. With commitment to simplicity, we follow sequential systems-based in silico experiments to hypothesize the essential targets that if hit together by co-drugging, will inhibit hallmarks of EMT in aggressive cancers (Figure 1). We also show how bioinformatic tools can be helpful in initial validation of the hypotheses. Finally, we proceed to direct experimental validation of drug combination efficacy in a three-dimensional (3D) coculture model of EMT and demonstrate that the anti-EMT efficacy of the suggested drug combination is significantly superior to single agents proposed so far in the literature for inhibition of EMT in lung cancer cells.

Figure 1

A systems-pharmacology workflow with integration of multiple bioinformatic tools sheds light on key players of EMT in cancer. (A) From published EMT signatures, genes associated with cell adhesion were extracted (BiNGO App in Cytoscape), converted to an up/down signature in various solid cancers (CancerMA) and extended by co-expression (cBioPortal) to identify drugs (L1000CDS²) that induce expression of epithelial-related genes while downregulate mesenchymal-associated genes. (B) EMT-associated kinome was identified by enrichment analysis (Expression2Kinases) and was used to create a kinase interaction network in EMT (STRING). Network analysis (Gephi) was performed to identify important kinases and their inhibitors (DRUGBANK). (C) The signatures of the identified drugs were extracted (iLINCS) and subjected to pathway enrichment analysis (Enrichr) to estimate drug combination efficacy. (D) Summary of drug classes in combination and their prospecting effects, which is hypothesized to inhibit multiple aspects in EMT.

It is noteworthy to mention that while the key findings of our straightforward approach were in line with previous literature derived from years of conventional wet-lab research; such approach is of paramount importance in transfiguring segmented data into a hypothesis-inspiring pipeline energizing another round of improved drug combination design for clinical research.

Results and discussion

Case study presentation and preparation of unified EMT gene set library

EMT is a complex reprogramming process through which carcinoma cells facilitate their dissemination and acquire stemness features and drug resistance [14, 15]. Despite recognition of EMT as the key machinery to evade environmental limitations in tumors, full understanding of the process and its clinical targeting is yet far from complete [16–18].

Cells in the EMT state are plastic and can easily change the expression profile of many genes according to the needs of their changing microenvironment. Thus, we are not faced with a single gene or even fixed list of genes that are upregulated or downregulated. That is why we did not choose low-throughput studies in our case presentation. For instance, in choosing a low-throughput study, Lou etal. [19] showed that murine 4T1 breast cancer cells represent mesenchymal behaviors indicative of undergoing EMT without necessarily expressing N-cadherin, which is well-known hallmark of EMT in other cell types. Thus, when trying to develop an effective drug combination to inhibit such complex and plastic process, it is important to hit ‘multiple players’ that are ‘robust’ in ‘various environmental conditions’ that is better reflected in high-throughput studies.

Among the high-throughput studies, although many research groups have proposed EMT signatures in various types of cancers, we used four pan-cancer EMT signatures [20–23] that provide a generic picture of EMT reflecting ‘key processes’ that are independent of tumor type and free from bias of the effect of stroma in tumor microenvironment [21]. Otherwise, the drug targets may be lost during the progression of EMT and metastasis to distant organs, and the cells may become resistant to drugs [18, 24]. In these four studies, two different statistical methods were applied for constructing the signatures. Byers group developed a 77-gene signature whose expressions were highly correlated with known EMT markers, including E-cadherin (epithelial marker) and N-cadherin, vimentin and fibronectin (mesenchymal markers) using The Cancer Genome Atlas (TCGA) database. Successively, Thiery’s research team proposed an EMT signature composed of 315 genes based on correlation analysis of EMT markers with cancer genes and developed a scoring system to quantify the extent of acquiring mesenchymal features in tumor samples. Two other studies in 2012 and 2016 also performed meta-analyses to integrate multiple gene expression datasets and proposed two gene lists manifesting two independent signatures, respectively. To propose a drug combination design, we first compiled EMT-associated gene signatures from four previously described studies to create a ‘meta-signature’ of pan-cancer 962 genes (for details, see Supplementary Methods) to be used as the input for two different drug target prediction pipelines in which we integrated gene set-based enrichment analyses and protein–protein interaction networks, respectively.

Case study part1: prediction of drugs to tune expression of cell adhesion-related genes

In the first pipeline to identify drugs that affect cell–cell interactions in EMT, we used enrichment-based tools such as ‘BinGO’ and ‘MSigDB’ as well as repository databases such as ‘CancerMA’, ‘cBioPortal’, ‘L1000CDS2’, ‘iLINCS’ and ‘Immgen’ to identify and validate drugs that affect cell interaction process in EMT, which are described in Box 1.

We thought that when the goal is to reverse or mimic a list of genes within a signature by use of drugs that mimic/inverse their expression, it would be more fruitful to extract subset of genes that together cooperate in a cellular function not just choosing all of the genes that are perturbed in a biological state. In our scenario, an important aspect of inhibiting EMT in invasive cells is to normalize the expression of epithelial versus mesenchymal set of genes and the intracellular machinery that regulates this balanced expression. Therefore, if we can reverse the expression of these sets of adhesion-related genes, we have attacked a principal mechanism of EMT initiation [25]. Taking into account the importance of cell adhesion-related genes in initiation of EMT, the biological process gene cluster under ‘Cell-adhesion, GO: 0007155’ term was dissected from the main library using ontology-based enrichment test in ‘BiNGO’, which resulted in a subset of 93 genes. The differential expression of the 93 subset genes was assessed by performing meta-analysis across 80 data sets and 13 cancer types in ‘CancerMA’ database. Complete list of differentially regulated cell adhesion genes across multiple cancers is provided in Supplementary File 1. This database also allowed us to convert the gene set into an upregulated (regarded as mesenchymal-related genes) and downregulated (regarded as epithelial-related genes) mode in a cancer-dependent manner.

The outputs of CancerMA are shown in Figure 2. We observed that 82 genes were deregulated in multiple solid cancers; however, no common cell adhesion-related gene was found to be deregulated in all cancer types of the study; e.g., a metalloproteinase (ADAM12) and fibronectin (FN1) were upregulated only in five cancer types (Figure 2A), while fibrillin (FBN1) was downregulated in seven cancer types and dermatopontin (DPT) and stromal cell-derived factor-1 (CXCL12) were downregulated in six and five cancer types, respectively (Figure 2B).

Figure 2

Meta-analysis from CancerMA database output shows inconsistent expression of cell adhesion-related genes across various solid cancers. (A) Upregulated cell adhesion-related genes dissected from collated EMT signatures. Expression status was determined by performing meta-analysis in CancerMA database for different cancer types. These upregulated genes were considered as mesenchymal genes (B) Downregulated cell adhesion-related genes from collated EMT signatures in various cancers obtained by meta-analysis. These downregulated genes were assumed to be associated with epithelial characteristics. The outer layer shows deregulated genes, and the inner layer represents cancer types. Edge thickness in the networks signifies fold change (FC) range, and ±1.5 FC was considered significant. Arrows indicate most similar genes between different cancers.

To identify drugs, which reverse cell adhesion-related genes in multiple cancers, the up/down lists for each cancer were submitted into ‘L1000CDS²’ search engine to explore LINCS expression data (Library of Integrated Network-based Cellular Signatures) of chemical perturbations. L1000CDS² prioritizes signature reversing drugs based on the overlap score calculated by characteristic direction method [26]. We observed that even though each cancer exhibited different set of expression patterns for cell adhesion-related genes, HDACIs, including trichostatin A, vorinostat, panobinostat and pracinostat reversed the expression of various cell adhesion-related genes in EMT meta-signature across 10 cancers compared with other classes of drugs (Figure 3A). Confirming these results, Tang et al. [27] observed that among multiple drug classes, HDAC inhibitors could restore E-cadherin expression with low cytotoxicity.

Figure 3

L1000CDS2 output and confirmations with MsigDB and Immgen databases show HDACIs are a notable class of drugs to reverse multiple EMT-associated processes by affecting cellular plasticity. (A) Illustration for the number of cancer types whose EMT induced cell-adhesion related genes are reversed by different classes of drugs. Each color represents a distinct class of drug, and the size of each colored box relatively shows the number of cancers whose EMT signature is affected by HDACIs. (B) Overlap between up/downregulated genes in embryonic plasticity and adult tissues, which are reversed by vorinostat. (C) GSEA of DEGs induced by vorinostat in multiple cell lines with stem cell signatures in MsigDB database version 6.0. (D) Heatmap and box-plot for mRNA expression of significantly downregulated genes by vorinostat in various mouse hematopoietic and stem cell lines in ImmGen project.

To further extend these deregulated genes with patient data in the in vivo conditions, the deregulated genes for each cancer were used as seed and were expanded with TCGA data sets (Supplementary File 2) using co-expression (Pearson’s r > 0.8) toolbox in ‘cBioPortal’ [28]. In resubmitting the expanded gene lists to L1000CDS², HDACIs were observed among the top ranked drugs to reverse the expression pattern of extended signatures in eight cancers, which are marked with asterisk in Table 1.

Table 1.

Drugs predicted by L1000CDS² to reverse EMT gene sets in various cancers

Cancer type	Drugs	Mechanism	Genes reversed by the drugs
Lung	Digoxigenin	Cardiovascular, inotropic	ADAM12, DSP, THBS2
	Vorinostat	HDACI*	ADAM12, COL4A3, FBLN5, COL4A3, MFAP4
	Sunitinib	Multi-kinase inhibitor*	DSP, THBS2, DLC1
Colorectal	Trametinib, saracatinib, linifanib	Multi-kinase inhibitor*	COL6A3, THBS2, AZGP1, COL6A3, COL4A3, CXCL12
	Vorinostat, pracinostat	HDACI*	AZGP1, COL6A3, THBS2, COL6A3, EZR, MFAP4, CEACAM1
	Albendazole	Anthelmintic*	THBS2, CEACAM1, EZR
Breast	Itraconazole	Antifungal*	COL14A1, MFAP4, SRPX
	Finasteride	5-alpha reductase inhibitor	COL14A1, DPT, PPAP2B
	Topotecan	Anticancer*	COL14A1, CXCL12, DPT
	Granisetron	Antiemetic	FN1, CXCL12, SRPX
	Regorafenib, dasatinib	Multi-kinase inhibitor*	FN1, PPAP2B, SRPX, CXCL12, MFAP4
	Mocetinostat	HDACI*	CXCL12, PPAP2B, SRPX
	Azacitidine	Anticancer	CXCL12, MFAP4, PPAP2B
Prostate	Vorinostat	HDACI	SPOCK1, CEACAM1
Pancreatic	Vorinostat	HDACI*	AEBP1, COL3A1, COL5A1, COL6A2, COL6A3, ITGB5, PERP, POSTN, THBS2, TNFAIP6, ADAM12, ITGB5
	Nintedanib	Multi-kinase inhibitor*	COL3A1, COL5A1, COL6A2, COL6A3, POSTN, THBS2, TNC, TNFAIP6
Head and neck	Vorinostat	HDACI*	ADAM12, CDH11, COL3A1, COL6A3, POSTN, TNFAIP6, VCAN, AEBP1
	Nintedanib	Multi-kinase inhibitor	CDH11, COL3A1, COL5A1, COL6A3, POSTN TNFAIP6, VCAN
	Cabergoline	Dopamine agonist	CDH11, COL5A1, POSTN, TNFAIP6, MPZL2
	Albendazole	Anthelmintic	AEBP1, CDH11, COL3A1, VCAN, CEACAM1
	Ruxolitinib	JAK 1 and 2 inhibitor	COL3A1, COL5A1, COL6A3, FN1, POSTN
Brain (GBM)	Vorinostat	HDACI*	COL13A1, COL5A1, COL6A1, COL6A2, VCAN, GPR56, AEBP1, PERP
	Saracatinib	Multi-kinase inhibitor*	ADAM12, COL5A1, FN1, GPR56, NEDD9, RHOB, FLRT3
	Cyclosporin A	Calcineurin inhibitor, immunosuppressant	AEBP1, COL14A1, COL5A1, NEDD9, TNC
	Parthenolide		AEBP1, COL3A1, COL5A1, FN1, NEDD9, RHOB
	Ivermectin	Anthelmintic	ADAM12, COL5A1, NEDD9, PLXNC1, RHOB
Adrenal	Vorinostat	HDACI	ADAM12, CXCL12, MFAP4
	Trifluoperazine	Antipsychotic	LOXL2, AZGP1, CXCL12, FBLN5
Renal	Vorinostat	HDACI*	COL4A3, CTGF, FBLN5, RHOB
	Topotecan	Topoisomerase inhibitor, anticancer	AZGP1, COL14A1, CXCL12, DPT
Ovarian	Vorinostat, Panobinostat	HDACI	ACTN1, FN1, POSTN, SRPX, THBS2
	Itraconazole	Antifungal	CDH11, MFAP4, SLIT2, SRPX, THBS2
	Doxorubicin*		NEDD9, PPAP2B, CDH11, COL3A1, FLRT2, NID2, SRPX, THBS2
Thyroid	Vorinostat	HDACI*	POSTN, CDH2, CTGF, FBLN5, SRPX, COL4A3, CXCL12, MFAP4
	Timolol	Antihypertensive	CDH2, FN1, DPT
	Doxepin	Tricyclic antidepressant	CDH2, DPT, FLRT2
	Digoxin	Cardiovascular, inotropic	FN1, POSTN, CTGF
	Trifluoperazine	Antipsychotic	FN1, CXCL12, FBLN5
	Sildenafil	Vasoactive agent	CDH2, FN1, FLRT2
	Azathioprine	Immunosuppressive	CDH2, CXCL12, DPT
	Itraconazole	Antifungal	CTGF, MFAP4, SRPX

*

Drugs that reverse co-expressed cell adhesion-related genes in EMT are shown with asterisk. GBM: Glioblastoma Multiforme.

We hypothesized that these drugs may modify the epigenetic mechanisms governing mesenchymal plasticity related to stemness features and pluripotency, which are conserved across many cancers [29, 30]. To confirm this hypothesis, we used data generated by Soundararajan et al. [31] developed a 160-gene signature of embryonic plasticity from epiblasts at Day 6.5, which represented the highest plasticity status and could mirror the capacity of tumors for distant metastasis and patient survival. The search for plasticity reversing drugs with L1000CDS² returned vorinostat, which reversed the expression of 13 genes related to plasticity (P-value in hypergeometric test = 0.011), while this drug did not affect the expression of a 664 adult tissue gene signature, which did not contain any of the embryonic-related genes (Figure 3B). Accordingly, downregulated genes by vorinostat in multiple cell lines were significantly enriched (FDR < 0.05) in stemness-related gene sets in ‘MSigDB’ (Figure 3C). Moreover, most significantly downregulated genes by vorinostat, obtained from ‘iLINCS’ database, were found to be expressed in stem cells and stromal fibroblasts compared with other mouse hematopoietic cell strains deposited in ‘ImmGen’ database (Figure 3D). These results pinpoint the role of epigenetic plasticity in aggressive behaviors of carcinoma cells, which would be inhibited by vorinostat.

Box 1. First bioinformatic pipeline showed histone deacetylase inhibitor (HDACI) drugs to be important in inhibiting EMT.

• BinGO (apps.cytoscape.org/apps/bingo): A Cytoscape plugin for performing gene ontology-based enrichment tests, including molecular function, biological processes and cellular components in Homo sapiens and a number of other biologically important species. The output is visualized as a graph and a table of significantly enriched terms. In our case study, we used BinGO to perform ontology-based enrichment analysis for extracting gene sets related to cell adhesion in EMT.

• CancerMA (cancerma.org.uk/information.html): A database with 80 manually curated cancer microarray data sets in 13 types of cancers and allows performing automatic meta-analysis on user-supplied gene lists. Owing to performing meta-analysis, the expression status of a gene can be more reliably obtained within a user-friendly interface. Here, we used CancerMA to assess differential expression status of cell adhesion-related genes and converting them to upregulation/downregulation mode.

• cBioPortal (www.cbioportal.org): A web resource to visualize, explore and analyze TCGA data through user-defined gene lists. It facilitates data integration because of containing various levels of data, including genomic and mutational data, epigenomic, transcriptomic and proteomic data and the clinical outcome. We used this pipeline to perform co-expression test, which helped us extend cell adhesion-related gene sets in EMT according to TCGA patient data.

• L1000CDS2 (amp.pharm.mssm.edu/L1000CDS2/): A web-based tool to prioritize drugs and ligands affecting user-provided genes based on chemical LINCS L1000 data. Similar to Connectivity map, L1000CDS2 nominates drugs that mimic or reverse a signature using characteristics direction method. Using this tool allowed us to identify Food and Drug Administration (FDA)-approved drugs that reverse the expression of upregulated mesenchymal genes and downregulated epithelial genes.

• iLINCS (ilincs.org/ilincs/): An integrative web-based platform to query and analyze genes, data sets and signatures against LINCS L1000 signatures for transcriptomic and proteomic data. We extracted L1000 data to obtain gene signature of selected drugs.

• Immgen (immgen.org/): A compendium of gene expression data sets of mouse hematopoiesis and immune cells with heatmap and box-plot visual output to demonstrate the expression of user-defined gene lists in various hematopoietic cell lines. Using this database, we could validate our hypothesis regarding the role of HDACIs in limiting cellular plasticity and stemness features.

• MSigDB (software.broadinstitute.org/gsea/msigdb): A collection of gene sets collected by Broad Institute, which allows performing gene set enrichment analysis (GSEA) to capture the biological function of large gene lists. By performing GSEA, we could validate the functional role of HDACIs gene signatures in differentiation and cellular plasticity.

Case study part2: elucidation of kinase interaction network in EMT and the drug inhibitors

As EMT is a stepwise process, we further analyzed the reversing drugs for a time-course data provided by Vetter et al. [32] who delineated early and late differentially expressed genes (DEGs) by SNAI upregulation as a model to induce EMT. We found trichostatin-A, dasatinib (SRC kinase inhibitor) and ketoprofen (a nonsteroid anti-inflammatory drug) as first-, fifth- and sixth-ranked drugs, respectively, to reverse early phase genes in EMT, and trichostatin-A as seventh-ranked drug to reverse late phase genes in EMT. These results not only underscore the significance of overcoming epigenetic barriers to switch from epithelial traits into mesenchymal features independent of the cancer type [33, 34] but also reinforce the contribution of multiple cellular kinases in accomplishment of EMT.

In our second bioinformatic pipeline (Box 2), we sought to extract the genes within the meta-signature based on their association with kinases. We thus determined EMT-associated kinome by performing kinase substrate enrichment embedded in ‘Expression2Kinases, X2K’ software. The identified kinases were connected using ‘STRING 10.0’ database to create a kinase interaction network (Figure 4A). Nodes and edges to construct the network can be found in Supplementary Files 3 and 4.

Figure 4

Integrated use of X2K, STRING and Gephi bioinformatic tools revealed SRC and IKBK as principal druggable kinases in EMT. (A) Illustrating the topology of SRC in the kinase network as a hub node and its first neighbor kinases in EMT-associated kinase interaction network based on centrality metrics degree, betweenness and closeness. (B) Identification of IKBK as another principal kinase in EMT functionally connected to three subcommunities in the network. The transparent nodes are the indirect neighbors of SRC and IKBK, while colored nodes are the directly connected nodes to IKBK and SRC, i.e. their first neighbors.

The network was consisted of five highly connected subcommunities (modules) regulating different pathways in EMT (Supplementary File 5). Calculation of main network centralities proper for our network [35], including degree, betweenness and closeness showed that the proto-oncogene SRC with 45 interactions was the hub node (a node with high degree of interactions) in the whole network. SRC also owned the highest betweenness (0.27) as a measure of shortest paths, which pass through this node and also highest closeness value (0.65) implying the role of SRC to be at the crossroad of multiple pathways (Figure 4B). Many of the identified kinases in the constructed network have already been proposed as druggable targets to inhibit EMT in separate parts of literature, including EGFR, PGFRR, EPHB receptors, YES1, LCK and AXL [36–42]. While it is not possible to inhibit all such redundant kinases in clinical settings, mathematical analysis of kinase interaction network showed that the addressed kinases are functionally connected to SRC family of kinases, which owned the highest degree and betweenness values in the network making them appealing druggable targets on which the information flow from the addressed upstream kinases is converged on.

Box 2. Second bioinformatic pipeline showed SRC and IKBK inhibitors to be important in EMT.

• X2K (maayanlab.net/X2K/): A multi-module tool that helps identifying additional data on regulatory patterns of transcription factors and kinases associated with genome-wide data. The tool performs enrichment tests on user-defined gene lists against libraries of transcription factors (ChipX enrichment analysis and position weight matrices methods) and kinase substrates database. It also allows for creating protein–protein interaction network from user-provided gene lists. Owing to embedding a wealth of gene libraries, X2K can also be used for GSEA of pathways and drugs. Using this tool helped us identify transcription factors and kinases associated with EMT gene signatures to create EMT kinase interaction network.

• STRING (string-db.org/): A comprehensive database of physical and functional protein–protein interactions. Queries of gene lists are returned as an interactive network based on experimental and computational predictions. Using STRING, we converted the EMT-associated kinases to protein–protein interaction network.

• Gephi (gephi.org/): A network visualization tool with powerful graphical capabilities to construct, manipulate and analyze graph-based biological networks. In the present case study, we used Gephi for analysis of network centrality metrics to identify the most influential nodes in the kinase interaction network.

• CREEDS (amp.pharm.mssm.edu/CREEDS/): This database archives crowdsourcing identification of DEGs in various biological conditions. The user interface allows for identifying similar/reverse gene signatures and conditions. CREEDS helped us validate the effect of manipulating selected nodes in kinase interaction network (SRC and IKBK) in inducing EMT status.

• DrugBank (drugbank.ca/): A comprehensive, well-curated source of FDA-approved drugs, experimental chemicals and their targets and associated pathways. Besides, DrugBank also contains metabolomic and pharmacogenomic and chemical structures of drugs. In present case study, we extracted FDA-approved drugs that inhibit selected kinases in EMT.

Closer inspection of the kinase interaction network revealed that by targeting SRC kinase alone, only four modules of the network (shown with blue, orange, yellow and green nodes) were manipulated, while an important subnetwork associated with cytokine–cytokine receptor interaction and TGF-beta signaling (red nodes) was unaffected. This finding may explain the inefficacy of SRC inhibitors in clinical trials [43]. From network analysis metrics, we observed that IKBKB owned the second rank for betweenness value (0.15) and the third rank for closeness (0.56) centrality confirming its important position in the network. IKBK was directly connected to SRC, TGFB-R, MAPK14 and CDK7, which were distributed among four other modules. The kinase interaction network elucidated that IKBK also interacts with ROCK1, AKT1, GSK and CHUK (IKK-α) kinases, which are dispersedly mentioned in the literature to fulfill EMT [44–47] (Figure 4C).

To further confirm the relevance of SRC and IKBK kinases in EMT, ‘CREEDS’ database of gene perturbation signatures [48] was queried, and DEGs by SRC and IKBK overexpression or constitutive activation were downloaded (for details of the queries, see ‘Methods’ section). Biological functions of DEGs in the SRC overexpression signatures were associated with cell adhesion and actin cytoskeleton organization. In IKBK, over-activation signatures, epithelial cell differentiation, extracellular matrix organization, cell junction organization and regulation of inflammatory response were observed (adjusted P-value < 0.05), all of which are evidently associated processes in EMT. Consistently, the overlap between DEGs in the SRC and IKBK perturbation signatures and the 962-EMT gene set was significant (P-value < 0.01) in hypergeometric statistical test (Supplementary File 6).

To identify FDA-approved drugs that inhibit SRC and IKBK targets, EMT-associated kinases were mapped onto drug–target network obtained from ‘DrugBank’ database [49, 50]. Four approved drugs, including dasatinib, bosutinib, nintedanib and ponatinib were found to inhibit SRC kinase family along with multiple other kinases such as PDGFR. Aspirin, sulfasalazine and mesalazine, currently used for treatment of inflammatory diseases were also found to inhibit IKBK. We thus chose dasatinib with 23 targets, including SRC kinase family as well as aspirin and sulfasalazine with 11 and 10 targets respectively, including IKBK, for further confirmation of their anti-EMT effects.

Case study part 3: computational and experimental validation of drug combination effect on inhibition of EMT

To elucidate the underlying pathways leading to significant inhibition of EMT by triple drug combination, we used different bioinformatic resources as described in Box 3. To this goal, gene expression changes of vorinostat, dasatinib and sulfasalazine were extracted from ‘iLINCS’ database for A549 and MCF-7 carcinoma cell lines. Expression changes for three drugs were also combined in each cell line to estimate the affected pathways following drug combination. We used ‘EnrichR’ tool to perform pathway enrichment analyses. Several pathways involved in EMT were enriched by combination of DEGs. In A549, a lung adenocarcinoma cell line with reversible mesenchymal traits, the combinatorial drug signature was associated with manipulation of several EMT-related pathways, including adherens junction and regulation of actin cytoskeleton, Foxo signaling, chemokine and NF-kB signaling (Figure 5A). Other pathways, including mismatch repair and cell cycle, were also affected by drug combinations, which are underlying mechanisms of EMT-induced drug resistance. However, in MCF-7, a hormone-responsive breast carcinoma cell line with epithelial traits only cell cycle, proteoglycans and Wnt signaling pathways were enriched in the combination signature (Figure 5B).

Figure 5

EnrichR tool output shows pathways involved in the anti-EMT effects of predicted drugs and their combination in two epithelial and mesenchymal cell lines. (A) KEGG pathway enrichment analysis with EnrichR in A549 and (B) MCF7 cells for DEGs following treatment by HDACIs, kinase inhibitors and their combination obtained from iLINCS database. Significantly enriched pathways are determined by FDR-adjusted q-value < 0.05). Unique pathways observed with drug combination are marked with dashed value. Color codes represent each drug and the combination.

Box 3. Third bioinformatic pipeline for computational validation of drug combination efficacy.

• iLINCS (ilincs.org/ilincs/): An integrative Web platform to query and analyze genes, data sets and signatures in LINCS L1000 data for transcriptomic and proteomic data. Using this database, we could access to L1000 data to extract signatures of our selected drugs.

• EnrichR (amp.pharm.mssm.edu/Enrichr/): Web-based software containing various biological libraries related to transcription factors, pathways, drug response and diseases. The tool allows fast enrichment analysis and visualization of provided gene lists against at least 60 libraries simultaneously. By EnrichR, we performed pathway enrichment analysis to understand affected pathways associated with selected drugs.

• CellMineR (discover.nci.nih.gov/cellminer/): A deposit of molecular profiles of NCI-60 panel of human cancer cell lines as the most widely used invitro models in cancer research. The profiles include DNA mutation, mRNA and protein expression as well response of cells to chemical and drug stimulation. User-defined gene lists are queried against various molecular and pharmacological data sets, and significant results are sent to the user. By CellMiner, we validated the role of selected drug target expression in epithelial versus mesenchymal ovarian cancer cell lines in NCI-60 panel.

• CSIOVDB (csibio.nus.edu.sg/CSIOVDB/CSIOVDB.html): Microarray database depositing expression data of 3431 human ovarian cancers with different molecular subtypes. Users can analyze the association of their queried genes to tumor stage and grade, molecular subtype, EMT score and patient survival in ovarian cancer. Using CSIOVDB, we confirmed the role of selected drug targets in stage and grade progression in ovarian cancer.

To further confirm the relevance of identified drug targets in EMT in another type of aggressive cancer, expression of drug combination targets was analyzed in ovarian cell lines of the NCI60 panel using ‘CellMiner’ database.

To rank ovarian cell lines of NCI60 panel based on their EMT status, an EMT score was attributed to each cell line by subtracting the z-score expression of vimentin from CDH1, as these genes showed to be the most anticorrelated mesenchymal and epithelial markers in the ovarian cell lines, respectively (r = −0.7), so that cell lines with positive score were assumed as mesenchymal. The expression values of the predicted drug targets were then extracted from CellMiner and compared across mesenchymal ovarian cell lines. A positive correlation (r > 0.5) between increase in the EMT score and expression of FYN (a member of SRC family of kinases), PDGFRB, HDAC2, CHUK and IKBK was observed all of which are inhibited by our proposed drug combination. These results imply that by increasing mesenchymal properties in these cell lines, the expression of these drug targets increases as well. To see if the pattern was maintained in data sets of clinical samples, the association of these drug targets with various histopathological characteristics of ovarian cancers in ‘CSIOVDB’ database was assessed (Supplementary File 7). Among the drug targets, FYN and PDGFRB, CHUK and HDAC2 were correlated with mesenchymal subtype. Moreover, expression of PDGFRB was associated with transition from Stage 3 to 4 and was more enriched in mesenchymal, stem-A and Stem-B subtypes of ovarian cancer. Higher expression of CHUK and PDGFRB were also associated with decreased survival. From the indirect drug targets with decreased expression in stromal cells by vorinostat, expression of COL1A1, MMP2, TIMP2 and HTRA1 was enriched in mesenchymal subtype of ovarian cancer among them, expression of LOXL1 and MMP2 was also associated with decreased survival. These results indicate that the triple-drug combination that we proposed inhibits multiple aspects of EMT, especially tumor invasion, which is defined as the primary event allowing malignant tumor cells to penetrate through the stroma and metastasize to distant loci [51]. The drug combination can ultimately improve patient survival by regulation of (1) extracellular matrix remodeling; (2) cell adhesion and motility; (3) cell cycle-related pathways and (4) inflammatory pathways in mesenchymal-like cells and rewiring stromal cell reprogramming and (5) inhibition of tumor cell plasticity and stemness.

To experimentally validate the anti-EMT efficacy of the hypothesized drug combination, it is nontrivial to consider the transient nature of EMT as well as its dependency on signaling from surrounding microenvironment and basement membrane [52]. We thus implemented a 3D microfluidic device with collagen-embedded A549 lung carcinoma spheroids in coculture with HUVEC endothelial cells to induce EMT and monitored dispersion of carcinoma cells from the spheroids for 84 h following drug treatment. We used A549 cells because of their plastic and mesenchymal behavior in response to environmental stimulations. This measurement was formerly shown to correlate with upregulation of vimentin and downregulation of E-cadherin as hallmarks of EMT [53] and recapitulates the transient and dynamic nature of EMT following the real-time monitoring of drug effects in a physiologically relevant model.

To quantify the amount of cancer cell dispersion following drug treatment, we counted the distance of each green pixel (A549 cells) from green pixels (HUVECS) in the images using MATLAB and averaged the distance of red pixels for each time point. The resulting numbers were plotted and shown beside the figure of each microscopy. As shown, single drugs either did not inhibit cell dispersion or only posed a delay when the distance was compared to control (no-drug treatment group) and the density of cells reduced in the spheroid because of invasion of cells from their initial spheroids to toward the HUVECS. However, in the combination experiment, the A549 cells maintained a high-density condition within collagen spheres meaning that they were not able to digest the collagen and migrate toward the HUVEC cells.

Treatment of A549 cells with 1 µM vorinostat alone did not show a significant decrease in cell dispersion from the collagen spheroids compared to control after 84 h (Figure 6A and B). However, migration of HUVEC cells from the lateral channel of the devices toward A549 was less observed at the interface of disseminated carcinoma cells and the proliferated HUVECs (Figure 6C). Vorinostat affects cellular plasticity and given that this phenomenon in carcinoma cells is congruent with a permissive microenvironment, we analyzed the biological processes in A549 cells, which were affected by vorinostat using iLINCS data. Vorinostat-induced gene signature was enriched with processes involving cell migration, chemotaxis and angiogenesis. It could also decrease the expression of COL1A1, COL4A1, MMP2, TIMP2, LOXL1, CDKN1, HTRA1 and RAB31, which are associated with activated fibroblasts and stromal cells [54–56]. Accordingly, 24% of HDACI-induced gene products were located in extracellular vesicular exosomes and in extracellular matrix, which affected cell adhesion, cell migration, chemotaxis, epithelial differentiation, angiogenesis and regulation of NF-kB kinase expression, which are known to be associated with EMT (Figure 6D and E).

Figure 6

Experimental validation of the effect of vorinostat monotherapy on cell dispersion. (A) Confocal live cell microscopy images showing dispersion of A549 lung cancer cells (with red fluorescent) cocultured with HUVECs (green fluorescent) following treatment with 1 µM vorinostat versus control in microfluid device as a metric of EMT. Scale bar indicates 50 µm. (B) Time-resolved dispersion of A549 cells quantified by pixel distance measurement in vorinostat or no treatment control group normalized by dividing dispersion distance at each time-point to 0 h. Data are represented as mean ±SEM for three spheroids in two independent biological replicates. T-test was performed to assess statistical difference between vorinostat and control at 84 h, P-value < 0.05 was considered significant. (C) Migration of HUVEC cells toward dispersed A549 cells in presence and absence of vorinostat after 84 h. (D) Cellular location of vorinostat-induced gene products using Enrichr. FDR-adjusted P-value (q-value) < 0.05 was considered significant for cellular component enrichment. (E). Significantly enriched (q-value < 0.05) biological processes for shared genes between collated EMT signatures and vorinostat-induced DEGs from iLINCS database.

We then tested the effect of SRC and IKBK inhibition alone. As we had predicted, inhibition of A549 cells dissociation from the spheroids was not statistically significant following treatment with 1 µM Dasatinib alone for 84 h compared to control no-treatment group (Figure 7A–C). Independently, aspirin alone at 1 µM concentration did not significantly inhibit dispersion. Moreover, 1 µM sulfasalazine, another identified IKBK inhibitor, which also inhibits CHUK kinase, imposed a delay on A549 dispersion for about 72 h compared to negative control (Figure 7D and E). Based on the insights from kinase interaction network analysis, however, combination of 1 µM sulfasalazine with dasatinib significantly inhibited cell dispersion from the spheroids in the coculture (Figure 7G and H). Interestingly, addition of vorinostat to combination of dasatinib and aspirin, with a low concentration of 1 µM of each drug, fully inhibited cell dispersion (Figure 7I and J) and a distance of 200 µm (equal to the space between the inputs of the two channels of the device) was kept between the two cocultured cell types (Figure 7K). It was also grossly evident from green pixel count in the images that the triple-combination treatment maintained survival of cocultured HUVEC endothelial cells, which denotes tolerability of proposed combination regimen for normal cells (Supplementary File 8).

Figure 7

Effect of kinase inhibitors and their combination with vorinostat on cell dispersion. (A and B) Experimental results of A549 dispersion in spheroids before and after treatment with 1 µM dasatinib as a SRC kinase inhibitor for 84 h cocultured with HUVECs, respectively. (C). Measurement of dispersion distance of A549 cells from collagen spheroids in dasatinib and no treatment control normalized by dividing to distance at 0 time-point. (D and E) Treatment of A549 spheroids and monitoring cell dispersion with 1 µM aspirin or sulfasalazine for 84 h as IKBK inhibitors. (F). Normalized dispersion distance following treatment with aspirin or sulfasalazine for 84 h. (G and H) Significant (t-test, P-value < 0.05) inhibition of A549 dispersion distance normalized to 0 h time-point from the spheroids following co-treatment with dasatinib and sulfasalazine at 84 h. (I and J) Assessment of cell dispersion following combination treatment with 1 µM of each drug for 84 h. Significant (t-test, P-value < 0.05) dispersion distance between null control versus combinatorial triple treatment at 84 h was observed. Data represent mean± SEM, and the error bars indicate dispersion measurements in two independent spheroids. Scale bar indicates 50 µm. (K) Distance between dissociated A549 and HUVEC cells normalized to 0 h as a measure of endothelial cell migration toward carcinoma cells.

The proposed drug cocktail is composed of FDA-approved drugs that are safely prescribed in patients. To see if these drugs are able to affect viability of mesenchymal cancer cells, we further analyzed the effect of proposed drug combination on viability of two non-small cell lung cancer cell lines, which are closely similar to A549 in their invasive behavior. We used dual labeling via acridine orange (AO) and propidium iodide (PI) nuclear staining in 3D microfluidic device to quantify viable cells. Cells were treated with 1 µM of vorinostat, dasatinib and aspirin or sulfasalazine either alone or in combination. We observed that the drug combination was significantly more effective in decreasing the cell viability of HCC-44 and H23 cells (Supplementary File 9).

Conclusion and future perspectives

Here, we practically showed how integration of multiple bioinformatic resources can connect disparate data into knowledge of designing a drug combination to inhibit a complex process such as EMT. In light of bioinformatic analyses, we observed EMT inhibition when IKBK was co-targeted beside SRC, and their combination with HDAC inhibitors was effective in limiting EMT consequences. Addition of vorinostat to this combination also confined invasion of endothelial cell toward spheroids and maintained the inhibition of dispersion from the spheroids. Moreover, by implementing a 3D coculture system to simulate tumor microenvironment and use of clinically achievable doses of already prescribed drugs, we observed agreement between in silico and in vitro results to confirm the relevance of proposed combinatorial regimen in abrogating a process such as EMT. Besides, given the safety of the proposed drugs in this study, as well as their FDA approval for other indications, the outcomes presented here are of direct translatability to clinical trials.

Methods

Compilation of published EMT signatures

EMT-associated gene signatures were obtained from the literature sources, including Byers group [57] and Thiery research group [21]; also from two meta-analyses performed in 2012 [22] and 2016 [23]. The gene lists in the four aforementioned signatures were integrated to create an EMT seed gene library.

Enrichment-based and analytical tools

Gene ontology-based enrichment tests were performed in ‘BiNGO’ Cytoscape plugin [58], GSEA in ‘Molecular Signatures Database (MSigDB) version 6.0’ (software.broadinstitute.org/gsea/msigdb) [59] and also ‘Enrichr’ (amp.pharm.mssm.edu/Enrichr) Web platform as of Feb 2017 [60]. ‘CancerMA’ database (www.cancerma.org.uk/information.html) accessed on February 2017 was used to perform meta-analysis across 80 cancer data sets for 13 different cancer types [61]. ‘L1000CDS^2’ search engine (amp.pharm.mssm.edu/L1000CDS2) accessed on February 2017 was used to query L1000 data sets and identify EMT reversing drugs [62, 63]. ‘Expression2Kinases (X2K)’ downloaded in November 2016 (www.maayanlab.net/X2K) was used to identify EMT-associated kinome, which uses enrichment analysis of kinases and their substrates [64]. Initially, a bottom-up approach was taken starting from identification of transcription factors as described previously in [65] and [66], which mainly identified intracellular non-receptor kinases linking signal transduction pathways to gene expression changes in the nucleus. Individual EMT gene signatures were also directly submitted to KEA module of X2K to identify kinases regulating protein functions at the extracellular and cell membrane levels. ‘CellMiner’ analysis tool version 2.1 [67] (discover.nci.nih.gov/cellminer) and ‘CSIOVDB’ database of ovarian cancer microarray gene expression version 1.0 [68] (csibio.nus.edu.sg/CSIOVDB/) were applied to query drug targets expression analysis across NCI-60 cancer cell line panel and ovarian cancer data sets.

Databases, public resources and software

‘cBioportal’ (www.cbioportal.org) database version 1.X.0 was used to expand DEGs in each cancer with TCGA data sets by performing co-expression [28]. Random gene lists were obtained ‘Random Gene Set Generator’ embedded in ‘molbiotools’ (http://www.molbiotools.com/randomgenesetgenerator.html). ‘STRING 10.0’ (https://string-db.org) database [69] was used to convert the kinase lists into a kinase–kinase interaction network [70]. Analysis of proper network centrality metrics [71–73] was performed and visualized in ‘Gephi 0.8.0’ [74]. ‘CREEDs’ database version 1.0 (amp.pharm.mssm.edu/CREEDS/) was queried, and GEO IDs, including GSE37428, GSE15161, GSE17511 and GSE26410, were extracted, which represented SRC and IKBK overexpression [48]. To access drug-related signatures, the ‘iLINCS’ web platform (http://ilincs.org) was queried to extract data sets of vorinostat (LINCSCP_34467 and LINCSCP_66500), dasatinib (LINCSCP_36498 and LINCSCP_6177) and sulfasalazine (LINCSCP_34467 and LINCSCP_5717) in MCF-7 and A549 cell lines. Expression of drug signatures in murine hematopoietic cells was analyzed with ‘my gene set’ tool in ‘ImmGen’ project (immgen.org) [75].

Experimental validation of anti-EMT effects of predicted drugs

All procedures to coculture A549 spheroids with HUVECs in 3D microfluid devices were performed as previously described in [53, 76]. Briefly, stable histone H2B-mCherry expressing human lung carcinoma cell line A549 (ATCC, USA) was maintained in DMEM supplemented with 10% FBS. HUVEC cells were (Lonza, Basel, Switzerland) grown in microvascular endothelial growth media (Lonza EGM-2 MV, Basel, Switzerland). All procedures were performed as previously described in [77, 78]. Briefly, for spheroid generation, A549 cells were seeded in ultra-low attachment dishes (Corning Inc, NY, USA) for 10 days and were collected when forming spheroids of 40–70 µm in diameter using a cell strainer to sieve desired size of spheroids. Media was then removed following centrifugation. In total, 200 µl type I collagen solution (BD Biosciences) at a concentration of 2.5 mg/ml and pH 7.4 was mixed with cell suspension medium containing 30–50 tumor spheroids and pipetted into the central channel of the device and allowed to solidify at 37°C for 40 min in humidified incubator. DMEM for supplementation of A549 cells and endothelial growth medium containing HUVEC cells was then injected into related channels and was allowed to form a monolayer in the device with a nearly 200 µm distance from A549 spheroids, which permited diffusion of HUVEC condition medium to the spheroids. Drugs chosen based on predictions described above were dissolved in DMSO and reached to a final concentration of 1 µM with media and were injected into the channel inlets either alone or in combination. DMSO alone was used as control. Concentration of drugs was selected based on relevance to achievable clinical doses. High-content imaging was performed after 0, 12, 36, 60 and 84 h via FluoView 1000 confocal microscopy (Olympus, Japan). The extent of cell dispersion from the spheroids was quantified by measuring mean distance between red pixels in the images and was normalized by dividing to the mean distance in zero time point.

For live/dead assay, HCC-44 and H23 cells, which are closely similar to A549 cells in their properties, were cultured and embedded into collagen spheroids as stated previously in the methods and injected into 3D microfluid devices. Cells were then treated with either each drug alone or their combination and compared with DMSO control for 5 days. AO/PI staining solution (Nexcelom, CS2-0106) was added to the spheroids in the device and incubated for 20 min in dark, and images were then captured by Nikon Eclipse 80i fluorescence microscope (Roper Scientific). AO penetrates to the nucleus of live and dead cells and generates a green background in cells, while PI is permeable to dead cells and stains to produce red fluorescent. Images were captured, and total area occupied by each color was quantified using AutoQuant Module [79].

Statistical analysis

For enrichment tests FDR-adjusted P-value < 5%, for meta-analysis −1.5 < fold change <+ 1.5, and in co-expression analyses Pearson’s correlation co-efficient > 0.8 was considered significant. In cell-based experiments, two tailed t-test or analysis of variance with statistical significance of < 0.05 in GraphPad Prism 6 (GraphPad Software; La Jolla, CA) was used. Quantification of cell-based experiments was performed with MATLAB 2014 b (MathWorks; Natick, MA).

Key Points

EMT is known to be the underlying mechanism of aggressive behaviors in tumors. Here, we proposed and validated a combinatorial drug therapy to inhibit EMT.

• We integrated network biology and current bioinformatic tools to design a druggable target combination from published EMT signatures.

• We found that regardless of the tumor type, drugs that modify histone acetylation and epigenetic modifications, including vorinostat (HDACIs), are able to tune expression of multiple epithelial and mesenchymal genes.

• We also observed that kinases serve as important class of druggable targets in EMT, thus by constructing an EMT-associated kinase interaction network, identified SRC and IKBK as the main kinases to regulate EMT whose activity is inhibited by dasatinib and aspirin or sulfasalazine.

• We validated the efficacy of drug combinations in inhibiting cell dispersion as a metrics of EMT and observed that three-drug combination not only inhibited cell dispersion from 3D spheroids but also limited invasion of cocultured endothelial cells toward the spheroids.

• Taking into account the safety profile of the proposed combination of FDA-approved drugs, these results offer a translational value in treatment of invasive tumors.

Farnaz Barneh is a PhD candidate of Proteomics in Department of Basic Sciences at Faculty of Paramedical Sciences in Shahid Beheshti University of Medical Sciences. She is interested in systems pharmacology and drug repurposing in cancer treatment.

Mehdi Mirzaie is a faculty member in the Department of Applied Mathematics, Faculty of Mathematical Sciences, Tarbiat Modares University. His research interests include mathematical biology, network science and structural bioinformatics.

Payman Nickchi is a PhD student of Statistics and is interested to apply statistical methods for big data analysis.

Tuan Zea Tan is based at National University of Singapore and is an expert in Systems Biology and Bioinformatic Approaches toward EMT in cancer.

Jean Paul Thiery is based at National University of Singapore is an expert in Systems Biology and Bioinformatic Approaches toward EMT in cancer.

Mehran Piran is working in cooperation with a team of biologists for statistical and omics data analysis at Drug Design and Bioinformatics Unit at Pasteur Institute of Iran.

Mona Salimi is interested in anticancer drug development and characterization at Department of Physiology and Pharmacology in Pasteur Institute.

Fatemeh Goshadrou is interested in the functional role of adhesion molecules in the mammalian cells in Department of Basic Sciences at Faculty of Paramedical Sciences, Shahid Beheshti University of Medical Sciences.

Amir R. Aref is interested in using microfluidic devices to study cancer biology and immunotherapy.

Mohieddin Jafari is an assistant professor at Pasteur Institute and is interested in multidisciplinary team working to apply systems pharmacology and bioinformatic approaches in understanding diseases and their treatments.

Acknowledgement

Contributors to publicly available open-source tools are greatly appreciated. The authors also acknowledge Prof. Avi Ma’ayan for helpful comments.

Funding

This study was not performed with public or institutional grants.