Comprehensive transcriptomic analysis of molecularly targeted drugs in cancer for target pathway evaluation.

Targeted therapy is a rational and promising strategy for the treatment of advanced cancer. For the development of clinical agents targeting oncogenic signaling pathways, it is important to define the specificity of compounds to the target molecular pathway. Genome-wide transcriptomic analysis is an unbiased approach to evaluate the compound mode of action, but it is still unknown whether the analysis could be widely applicable to classify molecularly targeted anticancer agents. We comprehensively obtained and analyzed 129 transcriptomic datasets of cancer cells treated with 83 anticancer drugs or related agents, covering most clinically used, molecularly targeted drugs alongside promising inhibitors of molecular cancer targets. Hierarchical clustering and principal component analysis revealed that compounds targeting similar target molecules or pathways were clustered together. These results confirmed that the gene signatures of these drugs reflected their modes of action. Of note, inhibitors of oncogenic kinase pathways formed a large unique cluster, showing that these agents affect a shared molecular pathway distinct from classical antitumor agents and other classes of agents. The gene signature analysis further classified kinome-targeting agents depending on their target signaling pathways, and we identified target pathway-selective signature gene sets. The gene expression analysis was also valuable in uncovering unexpected target pathways of some anticancer agents. These results indicate that comprehensive transcriptomic analysis with our database (http://scads.jfcr.or.jp/db/cs/) is a powerful strategy to validate and re-evaluate the target pathways of anticancer compounds.

Many cancer cells are addicted to driver oncogenes or to cancer-selective survival factors, and their proliferation and survival is highly dependent on oncogenic signaling pathways.1,2 Therefore, molecularly targeted drugs that selectively inhibit these pathways are critically important for the pharmacological treatment of advanced cancer.3 Presently, various inhibitors of oncogenic kinase pathways are available for the clinical treatment of cancer, such as inhibitors of oncogenic tyrosine kinases (for example, EGFR, HER2, BCR-ABL, and ALK), the RAF/MEK/ERK pathway, the PI3K/AKT/mTOR pathway, and multikinases.4 However, after treatment with each agent, cancer cells soon acquire drug-resistant phenotypes by several mechanisms including gatekeeper mutations in the target kinases and bypassing of signaling pathways.5,6 To improve treatment outcomes, additional next-generation inhibitors that possess better activity or overcome drug resistance to the primary agent should be further developed.

Target validation of agents is critically important for the development of new compounds as clinical antitumor agents. In the initial stages of drug development, high-throughput screens are usually carried out based on enzyme inhibition assays. As a result, candidate agents that have the potential to inhibit target enzymes are screened out. In some cases, however, the agents are found to affect additional target molecules in cancer cells and cause unexpected cytotoxicity during drug development or in clinical trials,7,8 which may mislead the selection of proper cancer subtypes for the agents and cause delay or failure in clinical trials. To ensure rational targeted therapy, target validation of compounds should be carried out with multiple reliable and unbiased methods.

Genome-wide gene expression analysis is an unbiased method to evaluate the mode of action of chemical compounds.9 We previously analyzed gene expression data of cancer cells that were mainly treated with classical antitumor agents, including DNA topoisomerase inhibitors, anti-metabolites, and tubulin-binding agents. We showed that the gene signature data reflected the modes of action of the respective agents.10 However, it is still not clear whether this signature-based analysis could widely be applied to classify the target pathways of molecularly targeted agents in cancer. To address these questions, in this study, we comprehensively obtained and analyzed gene expression data of cancer cells treated with 83 anticancer drugs or related agents covering most clinical (small molecule) anticancer drugs, such as oncogenic receptor tyrosine kinase inhibitors and other kinase inhibitors as well as inhibitors of promising molecular cancer targets. Our data indicated that this gene expression-based analysis efficiently classified the oncogenic kinase inhibitors as well as other classes of agents in a target pathway-dependent manner. Our data provide a platform to evaluate molecular pathways or primary cellular targets of compounds for further development of antitumor agents.

Materials and Methods

Cell lines and compounds

Human colon cancer HT-29 cells, ovarian cancer SKOV3 cells, leukemia K562 cells, and prostate cancer PC3 cells were obtained and cultured as described previously.10–12 Human lung cancer H2228 cells were obtained from ATCC (Manassas, VA, USA). Human lung cancer PC-9 cells were a kind gift from Dr. Kazuto Nishio (Department of Genome Biology, Kinki University Faculty of Medicine, Osaka, Japan).13 These cells were cultured in RPMI-1640 medium supplemented with 10% heat-inactivated FBS and 100 μg/mL kanamycin. The anticancer drugs or compounds used in our analysis are listed in Table1. The agents were obtained as described in Table S1. Stock solutions of the compounds were prepared using dimethyl sulfoxide as a solvent or as described previously.10 We examined the growth inhibitory effect of each agent (Fig. S1) and determined the GI50 values (Table S1). Growth inhibition assays were carried out and the GI50 values for each agent was determined as described previously.10

Table 1

Cancer cell line–anticancer drug combinations used in this study

Cell	Compound	Criteria	Target/Mode of action
K562	Imatinib	BCR-ABL inhibitor	BCR-ABL/KIT
	Dasatinib	BCR-ABL inhibitor	BCR-ABL/Src
	Nilotinib	BCR-ABL inhibitor	BCR-ABL
	Bosutinib	BCR-ABL inhibitor	BCR-ABL/Src
	Ponatinib	BCR-ABL inhibitor	BCR-ABL (T315I)
	SN-38	DNA damaging agent	Topoisomerase I
	Doxorubicin	DNA damaging agent	DNA intercalator/Topoisomoerase II
PC-9	Gefitinib	EGFR/HER2 inhibitor	EGFR
	Erlotinib	EGFR/HER2 inhibitor	EGFR
	Afatinib	EGFR/HER2 inhibitor	EGFR/HER2
	Trametinib	RAF/MEK/ERK inhibitor	MEK
	SN-38	DNA damaging agent	Topoisomerase I
	Doxorubicin	DNA damaging agent	DNA intercalator/Topoisomerase II
H2228	Crizotinib	ALK inhibitor	ALK
	Alectinib	ALK inhibitor	ALK
	SN38	DNA damaging agent	Topoisomerase I
	Doxorubicin	DNA damaging agent	DNA intercalator/Topoisomerase II
SKOV3	Lapatinib	EGFR/HER2 inhibtor	EGFR/HER2
	SN-38	DNA damaging agent	Topoisomerase I
	Doxorubicin	DNA damaging agent	DNA intercalator/Topoisomerase II
HT-29	Vemurafenib	RAF/MEK/ERK inhibitor	BRAF (V600E)
	Dabrafenib	RAF/MEK/ERK inhibitor	BRAF (V600E)
	Trametinib	RAF/MEK/ERK inhibitor	MEK
	U-0126	RAF/MEK/ERK inhibitor	MEK
	Everolimus†	PI3K/AKT/mTOR inhibitor	mTOR
	Temsirolimus†	PI3K/AKT/mTOR inhibitor	mTOR
	PP242†	PI3K/AKT/mTOR inhibitor	mTOR
	BKM120	PI3K/AKT/mTOR inhibitor	PI3K
	BEZ235	PI3K/AKT/mTOR inhibitor	PI3K/mTOR
	AKT Inhibitor VIII	PI3K/AKT/mTOR inhibitor	AKT 1/2
	Regorafenib	Multikinase inhibitor	VEGFR, RAF, KIT, RET etc
	Sorafenib	Multikinase inhibitor	VEGFR, RAF etc
	Pazopanib	Multikinase inhibitor	VEGFR, PDGFR,, KIT, FGFR etc
	Sunitinib	Multikinase inhibitor	VEGFR, PDGFR, KIT etc
	Cabozantinib	Multikinase inhibitor	VEGFR, MET,RET,KIT,FLT1/3/4 etc
	Vandetanib	Multikinase inhibitor	VEGFR, EGFR etc
	Axitinib	Multikinase inhibitor	VEGFR, KIT, PDGFR etc
	Gefitinib	EGFR/HER2 inhibitor	EGFR
	Erlotinib	EGFR/HER2 inhibitor	EGFR
	Afatinib	EGFR/HER2 inhibitor	EGFR/HER2
	Lapatinib	EGFR/HER2 inhibitor	EGFR/HER2
	Crizotinib	ALK inhibitor	ALK
	Alectinib	ALK inhibitor	ALK
	SU11274	MET inhibitor	MET
	AG1024	IGFR inhibitor	IGF1R
	PDGFR inhibitor V	PDGFR inhibitor	PDGFR
	Dasatinib	BCR-ABL/Src inhibitor	BCR-ABL/Src
	CDK4 inhibitor	Cell cycle inhibitor	CDK4
	NU6102	Cell cycle inhibitor	CDK1/Cyclin B
	ATM/ATR kinase inhibitor	DNA damage check point inhibitor	ATM,ATR
	SB218078	DNA damage check point inhibitor	CHK1
	CHK2 inhibitor II	DNA damage check point inhibitor	CHK2
	GSK-3 inhibitor IX	GSK-3 inhibitor	GSK-3
	FH535	β-catenin/TCF inhibitor	β-catenin/TCF
	Celecoxib	COX2 inhibitor	COX2
	BI 2536	Mitosis inhibitor	Polo-like kinase
	Aurora kinase inhibitor III	Mitosis inhibitor	Aurora kinase
	Docetaxel†	Mitosis inhibitor	Tubulin
	Paclitaxel†	Mitosis inhibitor	Tubulin
	Vincristine†	Mitosis inhibitor	Tubulin
	Trichostatin A†	HDAC inhibitor	HDAC
	Vorinostat†	HDAC inhibitor	HDAC
	Romidepsin	HDAC inhibitor	HDAC
HT-29	5-Aza-2′-deoxycytidine	DNA methyltransferase inhibitor	DNA methyltransferase
	Decitabine	DNA methyltransferase inhibitor	DNA methyltransferase
	Bortezomib†	Proteasome inhibitor	Proteasome
	Carfilzomib	Proteasome inhibitor	Proteasome
	MG-132†	Proteasome inhibitor	Proteasome
	MLN-4924	Nedd8 conjugation inhibitor	Nedd8 activating enzyme
	17-AAG†	Hsp90 inhibitor	Hsp90
	Geldanamycin†	Hsp90 inhibitor	Hsp90
	PKR inhibitor	RNA-dependent protein kinase inhibitor	RNA-dependent protein kinase (PKR)
	Ruxolitinib	JAK inhibitor	JAK
	TX-1918	Eukaryotic elongation factor-2 kinase inhibitor	Eukaryotic elongation factor-2 kinase (eEF2K)
	Vismodegib	Hedgehog pathway inhibitor	SMO
	SN-38†	DNA damaging agent	Topoisomerase I
	Doxorubicin†	DNA damaging agent	DNA intercalator/Topoisomerase II
	Camptothecin†	DNA damaging agent	Topoisomerase I inhibitor
	Topotecan†	DNA damaging agent	Topoisomerase I inhibitor
	Mitoxantrone†	DNA damaging agent	DNA intercalator/Topoisomerase II
	Etoposide†	DNA damaging agent	Topoisomerase II inhibitor
	Amrubicin	DNA damaging agent	Topoisomerase II inhibitor
	Cisplatin†	DNA damaging agent	DNA cross-linker
	Melphalan†	DNA damaging agent	DNA cross-linker
	Oxaliplatin†	DNA damaging agent	DNA cross-linker
	Neocarzinostatin†	DNA damaging agent	DNA cleavage
	Bleomycin†	DNA damaging agent	DNA cleavage
	Nimustine†	DNA damaging agent	DNA alkylator
	Mitomycin C†	DNA damaging agent	DNA alkylator
	5-FU†	DNA damaging agent	Pyrimidine
	Gemicitabine†	DNA damaging agent	Pyrimidine
	Methotrexate†	DNA damaging agent	DHFR
	6-Mercaptopurine†	DNA damaging agent	Purine
	Actinomycin D†	DNA damaging agent	DNA replication/RNA synthesis
	Pemetrexed†	DNA damaging agent	DNA/RNA synthesis
	2-Deoxyglucose†	ER stress inducer	Glycolysis
	Tunicamycin†	ER stress inducer	N-glycosylation
	Thapsigargin†	ER stress inducer	SERCA
	A23187†	ER stress inducer	Ca2+ ionophore

Gene expression data of these compounds were reported previously.10 17-AAG, 17-N-allylamino-17-demethoxygeldanamycin; AKT, protein kinase B; ALK: anaplastic lymphoma kinase; ATM, ataxia telangiectasia mutated; ATR, ataxia telangiectasia and Rad3-related protein; BCR-ABL, fusion gene of breakpoint cluster region protein (BCR) and Abelson murine leukemia viral oncogene homolog (ABL); CDK4, cyclin-dependent kinase 4; CHK, checkpoint kinase; DHFR, dihydrofolate reductase; EGFR, epidermal growth factor receptor; ER, endoplasmic reticulum; FGFR, fibroblast growth factor receptor; 5-FU, 5-fluorouracil; GSK3, glycogen synthase kinase 3; HDAC, histone deacetylase; HER2, human EGFR-related 2; Hsp90, heat shock protein 90; IGF1R, insulin-like growth factor 1 receptor; KIT, mast/stem cell growth factor receptor; MET, hepatocyte growth factor receptor; mTOR, mammalian target of rapamycin; PDGFR, platelet-derived growth factor receptor; PI3K, phosphoinositide 3-kinase; PKR, protein kinase RNA-activated; SERCA, sarco/endoplasmic reticulum Ca2+-ATPase; SMO, smoothened; T-cell factor (TCF); VEGFR, vascular endothelial growth factor receptor.

Drug treatment and GeneChip analysis

For gene expression analysis, we chose a concentration of drugs that were 3- to 10-fold greater than the GI50 value and caused >80% growth inhibition after 48 h of treatment, and gene expression data were obtained after 6 h of treatment10 Drug treatment concentrations and treatment duration for each agent are summarized in Table S1. Total RNA was extracted using an RNeasy Mini kit (Qiagen, Hilden, Germany). Microarray analysis was carried out as described previously with the GeneChip Human Genome U133 Plus 2.0 array (Affymetrix, Santa Clara, CA, USA).10 The signature data will be released on our website (http://scads.jfcr.or.jp/db/cs/).

Statistical analysis

All analyses were carried out using the statistical programming language R version 2.15.0 (http://www.r-project.org/) and Bioconductor version 2.10 (http://bioconductor.org/).

Data preprocessing

The R package software of Affymetrix Microarray Suite 5.0 was used to generate signal intensities for each of the HG-U133 Plus 2.0 arrays in the study. Expression values were normalized to a mean target level of 100.

Identifying gene signatures

Gene sets were extracted and classified as up- or downregulated after exposure to the drug. For each treatment sample, we calculated treatment-to-control ratio statistics, where, if any intensity value was <50, the value was replaced as 50. For the hierarchical clustering and the principal component analyses, we selected probe sets if the treatment-to-control ratio was >3 for upregulated genes or less than one-third for downregulated genes and the intensity of at least the treatment or control was >300 (Table S2).

To identify the signature gene sets characteristic of some drug subsets, we extracted the probe sets whose expression changes after drug treatment were statistically significantly different between the drug subsets and other agents. Statistical evaluations were carried out using Student’s t-test. Probes with more than a twofold differential expression and a P-value of <0.05 were extracted.

Hierarchical clustering

Probe sets for hierarchical clustering comprised the collection of all gene signatures. We carried out hierarchical clustering using the logarithm of the sample and probe set ratio statistics. We used Ward’s method for linkage and Pearson’s correlation for distance metric.

Principal component analysis

We carried out a principal component analysis based on the cancer cell gene expression data to examine 3-D clustering patterns of subclasses of anticancer drugs (oncogenic kinase inhibitors, HDAC inhibitors, proteasome inhibitors, mitosis inhibitors, and DNA damaging agents). We plotted cancer cells treated with anticancer drugs in a 3-D space consisting of three principal components. We used the same probe sets used for hierarchical clustering in this analysis.

Gene ontology analysis

To interpret the extracted gene signatures, we used gene ontology analyses using the DAVID analytical tool;14,15 this analysis is a method of highlighting relevant gene ontology terms associated with a given gene signature.

Analysis with the C-map algorithm (connectivity scoring analysis)

To investigate the relationship between gene signature and compound, we adopted the connectivity score based on the Kolmogorov–Smirnov statistic as developed by Lamb et al.9 For each treatment sample, all probe sets were ranked based on the treatment-to-control ratio and the rank matrix was configured using a similar method to Lamb et al.9 We modified our program and calculated the connectivity scores for all compounds, as described previously.10

Western blot analysis

Cells were lysed in TNE buffer (150 mM NaCl, 1.0% NP-40, 1 mM EDTA, and 10 mM Tris–HCl, pH 8.0) supplemented with 1× protease inhibitor cocktail (Nacalai Tesque, Kyoto, Japan) and PhosSTOP phosphatase inhibitor cocktail (Roche, Mannheim, Germany). Western blot analysis was carried out as described previously,16 using the following primary antibodies: anti-phospho-p70S6 kinase (p70S6K), anti-p70S6K, anti-phospho-AKT, anti-AKT, anti-phospho-ERK (Cell Signaling Technology, Danvers, MA, USA), anti-ERK (Santa Cruz Biotechnology, Santa Cruz, CA, USA), and anti-β-actin (Sigma, St. Louis, MO, USA).

Results

Comprehensive collection of gene expression data related to molecularly targeted, anticancer drug effects

In our previous analysis, we obtained gene expression data from human colon cancer HT-29 cells treated with 35 compounds mainly consisting of classical antitumor agents.10 For comprehensive transcriptomic analysis, we further obtained gene expression data from cancer cells treated with the most commonly used clinical molecularly targeted anticancer drugs, such as inhibitors of driver oncogenes (EGFR, HER2, BCR-ABL, ALK), RAF/MEK/ERK pathway inhibitors, PI3K/AKT/mTOR pathway inhibitors, multikinase inhibitors, HDAC inhibitors, and proteasome inhibitors (Table1). Alongside the anticancer compounds that are presently in clinical trials, we also included “promising” next-generation targeted inhibitors in our analysis, such as inhibitors of several receptor tyrosine kinases (MET, IGF1R, PDGFR), regulators of the cell cycle/check point (CDK4, ATM/ATR, CHK1, CHK2, Aurora kinase, and Polo-like kinase), β-catenin/TCF, COX2, and NAE.17–20 We used HT-29 cells because it is a commonly used, solid tumor cell line and we have used it in our previous analyses.10 We obtained transcriptomic data for all the agents in HT-29 cells with the exception of the BCR-ABL inhibitors that did not suppress HT-29 cell proliferation. Moreover, in the cases of drugs whose primary targets preferentially exist in specific types of cancer cell lines, we also used additional cell lines such as BCR-ABL-positive K562 cells21 for the BCR-ABL inhibitors, mutant EGFR-expressing PC-9 cells13 for the EGFR inhibitors, EML4-ALK fusion-positive H2228 cells for the ALK inhibitors,22 and HER2-overexpressing SKOV3 cells23 for the HER2 inhibitor. To estimate the effect of cell type difference on the gene expression analysis, we treated the cell lines with SN38 and doxorubicin, and obtained gene expression data as reference data (the gene expression data will be released on our website, http://scads.jfcr.or.jp/db/cs/).

Gene signatures reflect the target pathways of molecularly targeted drugs

As summarized in Table S2, we extracted genes whose expression was up- or downregulated by the analyzed agents. To compare the gene expression data of the compounds, we carried out a hierarchical clustering analysis with the acquired 129 gene expression datasets for cancer cells treated with 83 agents (4869 probe sets whose expression was up- or downregulated more than threefold in at least one of the datasets). As shown in Figure1, we observed that the compounds targeting similar molecules or molecular pathways were clustered together, such as DNA damaging agents, HDAC inhibitors, proteasome inhibitors, and inhibitors of mitosis-related molecules. These results indicate that the gene expression signatures reflect the primary target pathways of the drugs, as shown in our previous study.10 Moreover, in this study, we found that most of the inhibitors of oncogenic kinase pathways formed a large cluster distinct from classical antitumor agents or from other classes of agents. The data for the oncogenic kinase inhibitors in K562, PC-9, and SKOV3 cells were also clustered together with those of the oncogenic kinase inhibitors in HT-29 cells, whereas the data for the DNA damaging agents, SN-38 and doxorubicin, in multiple cancer cell lines were clustered together. Principal component analysis confirmed that the kinase inhibitors were clustered together and that this cluster was distinct from those of other classes of agents (Fig.2). These data indicate that the kinase inhibitors affect a shared molecular pathway in cancer cells distinct from other classes of antitumor agents. We further extracted signature genes whose expression was commonly modified by oncogenic kinase inhibitors (Table S3). Subsequent gene ontology analysis with the DAVID bioinformatics database revealed that several categories of genes, such as those involved in transcriptional regulation or apoptosis, were enriched in the signature genes (Table2).

Fig 1

Hierarchical clustering analysis based on 129 gene expression datasets of cancer cells treated with 83 anticancer drugs or related agents. For the analysis, we selected and used 4869 probe sets as gene signatures if the treatment-to-control ratio was greater than 3 for upregulated genes or less than one-third for downregulated genes and the intensity of at least the treatment or control was greater than 300 in at least one of the datasets. The values in the heat map are the logarithm values of the sample-to-control ratio of intensity values. Orange bars indicate 16 h of treatment samples. For agents with two treatment dosages, the samples of higher dosage are shown with asterisks. ER, endoplasmic reticulum; HDAC, histone deacetylase.

Fig 2

Principal component analysis based on gene expression data of cancer cells treated with subclasses of anticancer drugs. The subclasses contained a total of 73 datasets for oncogenic kinase inhibitors, HDAC inhibitors, proteasome inhibitors, tubulin-binding agents, and DNA damaging agents. In the principal component analysis, we plotted the data in a 3-D space consisting of three principal components.

Table 2

Gene ontology (GO) analysis of oncogenic kinase inhibitor signature genes

GO term	P-value	FDR
GO:0009952 anterior/posterior pattern formation	0.0004	0.0052
GO:0003002 regionalization	0.0013	0.0185
GO:0048806 genitalia development	0.0014	0.0204
GO:0045944 positive regulation of transcription from RNA polymerase II promoter	0.0019	0.0274
GO:0006355 regulation of transcription, DNA-dependent	0.0023	0.0332
GO:0007242 intracellular signaling cascade	0.0025	0.0353
GO:0042127 regulation of cell proliferation	0.0025	0.0355
GO:0051252 regulation of RNA metabolic process	0.0028	0.0397
GO:0042981 regulation of apoptosis	0.0028	0.0400
GO:0043067 regulation of programmed cell death	0.0030	0.0422
GO:0010941 regulation of cell death	0.0030	0.0431
GO:0043065 positive regulation of apoptosis	0.0036	0.0513
GO:0043068 positive regulation of programmed cell death	0.0037	0.0528
GO:0010942 positive regulation of cell death	0.0038	0.0538
GO:0007389 pattern specification process	0.0039	0.0549
GO:0010557 positive regulation of macromolecule biosynthetic process	0.0045	0.0638
GO:0045893 positive regulation of transcription, DNA-dependent	0.0056	0.0785
GO:0031328 positive regulation of cellular biosynthetic process	0.0057	0.0793
GO:0051254 positive regulation of RNA metabolic process	0.0058	0.0812
GO:0007548 sex differentiation	0.0058	0.0815
GO:0009891 positive regulation of biosynthetic process	0.0061	0.0848

Signature probe sets whose expression changes after drug treatment were significantly different between the oncogenic kinase inhibitors and other agents were extracted based on the Student’s t-test (fold-change values of more than 2 and the P-value of less than 0.05). We carried out GO analyses using the DAVID analytical tool to extract relevant GO terms associated with the gene signature. FDR, false discovery rate.

Classification of oncogenic kinase inhibitors based on gene expression signature

To examine whether the gene signature analysis could further distinguish the kinase inhibitors depending on their modes of action, we next focused on the gene signatures in HT-29 cells. As shown in Figure3, within the kinome-targeted agents, drugs with similar target pathways were clustered together, such as: (i) RAF/MEK/ERK pathway inhibitors; (ii) PI3K/AKT/mTOR pathway inhibitors; (iii) EGFR/HER2 inhibitors; (iv) multikinase inhibitors targeting VEGFR and PDGFR (shown as “multikinase inhibitors (1)” in Fig.3); and (v) multikinase inhibitors targeting VEGFR and RAF (shown as “multikinase inhibitors (2)” in Fig.3). Analyses of the gene expression signatures of BEZ235, vemurafenib, and gefitinib with the C-map algorithms further confirmed that the signatures of these agents were significantly similar to those of other drugs targeting the same or similar pathways in HT-29 cells (Table3A–C). These results indicated that the gene signature analysis could classify the kinome-targeted agents in a target pathway-dependent manner.

Fig 3

Hierarchical clustering analysis of the gene signatures of HT29 cells treated with 38 kinome-targeted drugs. For the analysis, we selected 2458 probe sets as gene signatures if the treatment-to-control ratio was greater than 3 for upregulated genes or less than one-third for downregulated genes and the intensity of at least the treatment or control was greater than 300 in at least one of the datasets. The values in the heat map are the logarithm values of the sample-to-control ratio of intensity values. Orange bar indicates 16 h of treatment sample. For the agents with two treatment dosages, the samples of higher dosage are shown with asterisks. AKT, protein kinase B; ALK, anaplastic lymphoma kinase; EGFR, epidermal growth factor receptor; HER2, human EGFR-related 2; mTOR, mammalian target of rapamycin; PI3K, phosphoinositide 3-kinase; RAF.

Table 3

Compounds similar to (A) BEZ235, (B) vemurafenib, (C) gefitinib (10 μM in HT29 cells) and (D) gefitinib (0.6 μM in PC-9 cells) with regards to gene expression changes after treatment

Rank	Cell	Compound	Concentration	Unit	Score	Up_score	Down_score
(A)
1	HT-29	BEZ235	1.00E-06	M	1.00000	0.9979	−0.99954
2	HT-29	BKM120	3.00E-06	M	0.97648	0.97167	−0.97878
3	HT-29	AKT Inhibitor VIII	1.00E-05	M	0.89995	0.96704	−0.83055
4	HT-29	Temsirolimus	1.00E-05	M	0.87040	0.80446	−0.93412
5	HT-29	PP242	1.00E-05	M	0.85915	0.82501	−0.89109
6	HT-29	6-Mercaptopurine	1.00E-04	M	0.84061	0.70322	−0.97586
7	HT-29	Cabozantinib	3.00E-05	M	0.81859	0.77389	−0.86118
8	HT-29	Crizotinib	1.00E-05	M	0.80053	0.80713	−0.79189
9	HT-29	Lapatinib (10 μM)	1.00E-05	M	0.79923	0.64398	−0.95245
10	HT-29	ATM&# x002F;ATR kinase inhibitor	1.00E-05	M	0.79442	0.7299	−0.85690
11	HT-29	Methotrexate	1.00E-06	M	0.78992	0.70036	−0.87746
12	HT-29	Sorafenib	1.00E-05	M	0.77319	0.56009	−0.98431
13	HT-29	Everolimus	1.00E-05	M	0.75731	0.78196	−0.73073
14	HT-29	Vandetanib	1.00E-05	M	0.74860	0.72206	−0.77323
15	PC-9	Gefitinib (30 μM)	3.00E-05	M	0.74129	0.63593	−0.84476
(B)
1	HT-29	Vemurafenib	3.00E-05	M	1.00000	0.99762	−0.99770
2	HT-29	Cabozantinib	3.00E-05	M	0.95797	0.96799	−0.94347
3	HT-29	U-0126	3.00E-05	M	0.93570	0.89185	−0.97516
4	HT-29	Dabrafenib	1.00E-05	M	0.87493	0.80785	−0.93791
5	HT-29	Vandetanib	1.00E-05	M	0.86775	0.86776	−0.86367
6	HT-29	Sunitinib	1.00E-05	M	0.85795	0.82050	−0.89138
7	HT-29	Sorafenib	1.00E-05	M	0.84555	0.83440	−0.85274
8	HT-29	Regorafenib	3.00E-05	M	0.81791	0.74200	−0.89000
9	HT-29	PDGF inhibitor V	1.00E-05	M	0.77796	0.83640	−0.71588
10	HT-29	Gefitinib (30 μM)	3.00E-05	M	0.77393	0.75086	−0.79339
11	HT-29	Pazopanib	3.00E-05	M	0.74890	0.69189	−0.80240
12	HT-29	Gefitinib (10 μM)	1.00E-05	M	0.74553	0.75449	−0.73308
13	HT-29	PP242	1.00E-05	M	0.73347	0.65308	−0.81042
14	HT-29	AKT inhibitor VIII	1.00E-05	M	0.72928	0.76986	−0.68529
15	HT-29	Erlotinib	3.00E-05	M	0.72384	0.70626	−0.73802
(C)
1	HT-29	Gefitinib (10 μM)	1.00E-05	M	1.00000	0.99927	−0.99945
2	HT-29	Gefitinib (30 μM)	3.00E-05	M	0.96045	0.96838	−0.95129
3	HT-29	Erlotinib	3.00E-05	M	0.94112	0.99669	−0.88435
4	HT-29	Sunitinib	1.00E-05	M	0.93169	0.99170	−0.87050
5	HT-29	Sorafenib	1.00E-05	M	0.91256	0.94111	−0.88283
6	HT-29	Pazopanib	3.00E-05	M	0.90385	0.8882	−0.91834
7	HT-29	Lapatinib (10 μM)	1.00E-05	M	0.89179	0.95223	−0.83022
8	HT-29	PDGF inhibitor V	1.00E-05	M	0.80332	0.83498	−0.77063
9	HT-29	Dasatinib	1.00E-07	M	0.76608	0.58031	−0.95086
10	HT-29	Thapsigargin	1.00E-08	M	0.74753	0.95102	−0.54308
11	HT-29	Vandetanib	1.00E-05	M	0.74082	0.89791	−0.58278
12	HT-29	AG1024	3.00E-05	M	0.73856	0.93070	−0.54548
13	HT-29	Vemurafenib	3.00E-05	M	0.72601	0.89795	−0.55314
14	PC-9	Erlotinib (30 μM)	3.00E-05	M	0.70436	0.75877	−0.64905
15	HT-29	Tunicamycin	3.00E-06	g/mL	0.68796	0.88138	−0.49367
(D)
1	PC-9	Gefitinib (0.6 μM)	6.00E-07	M	1.00000	0.99652	−0.99634
2	PC-9	Erlotinib (0.6 μM)	6.00E-07	M	0.98035	0.96886	−0.98486
3	PC-9	Erlotinib (30 μM)	3.00E-05	M	0.93176	0.93554	−0.92133
4	PC-9	Gefitinib (30 μM)	3.00E-05	M	0.92112	0.92387	−0.91180
5	PC-9	Afatinib	3.00E-08	M	0.86916	0.82167	−0.91045
6	PC-9	Trametinib	1.00E-06	M	0.60445	0.45342	−0.75116
7	HT-29	U-0126	3.00E-05	M	0.60392	0.58254	−0.62100
8	HT-29	Cabozantinib	3.00E-05	M	0.59445	0.56836	−0.61630
9	HT-29	Vemurafenib	3.00E-05	M	0.58051	0.53379	−0.62308
10	HT-29	PP242	1.00E-05	M	0.54759	0.55305	−0.53822
11	HT-29	Vandetanib	1.00E-05	M	0.52336	0.54008	−0.50290
12	HT-29	Dabrafenib	1.00E-05	M	0.51811	0.40252	−0.63000
13	HT-29	Sunitinib	1.00E-05	M	0.51786	0.45561	−0.57641
14	HT-29	Gefitinib (30 μM)	3.00E-05	M	0.50888	0.5259	−0.48823
15	HT-29	PP242	1.00E-05	M	0.50778	0.5421	−0.46985

AKT, protein kinase B; ATM, ataxia telangiectasia mutated; ATR, ataxia telangiectasia and Rad3-related protein; PDGF, platelet-derived growth factor. Compounds in our data that showed high similarity in their gene signatures to the given compounds were extracted using C-map algorithms. Top 15 data among the acquired 129 datasets are shown.

Moreover, we also evaluated the signature of gefitinib obtained in the mutant EGFR-expressing PC-9 cells using the C-map algorithms. The gefitinib signature of PC-9 cells showed significant similarity to those of oncogenic kinase inhibitors of HT-29 cells, including the gefitinib signature in HT-29 cells, while the top hits were other EGFR inhibitors of PC-9 cells (Table3D). These data indicated that, for the agents whose targets are selectively expressed in certain subtypes of cancer, use of data obtained in specific cancer cell lines could aid accurate evaluation of the drug target pathways based on the signature analysis.

To observe biological differences in the signature genes between subclasses of the kinase inhibitors, we further extracted genes that showed significantly selective expression in cells treated with RAF/MEK/ERK and PI3K/AKT/mTOR pathway inhibitors (Table S4). Gene ontology analysis revealed characteristic features of each gene set (Table4). Namely, the gene set specific for the RAF/MEK/ERK pathway inhibitors not only contained genes related to cell proliferation, protein kinase cascades, and cell death, but also genes involved in phosphate metabolic processes. In contrast, the gene set specific for the PI3K/AKT/mTOR pathway inhibitors was characteristically related to erythrocyte homeostasis, response to hypoxia, and angiogenesis, as well as cell proliferation and protein kinase cascades.

Table 4

Gene ontology (GO) analysis of signature genes of (A) RAF/MEK/ERK inhibitors and (B) phosphoinositide 3-kinase/protein kinase B/mammalian target of rapamycin (PI3K/AKT/mTOR) inhibitors

GO term	P-value	FDR
(A)
GO:0042127 regulation of cell proliferation	<0.0001	<0.0001
GO:0008285 negative regulation of cell proliferation	<0.0001	0.0002
GO:0006469 negative regulation of protein kinase activity	0.0001	0.0011
GO:0033673 negative regulation of kinase activity	0.0001	0.0013
GO:0007243 protein kinase cascade	0.0001	0.0013
GO:0051348 negative regulation of transferase activity	0.0001	0.0017
GO:0043407 negative regulation of MAP kinase activity	0.0007	0.0115
GO:0008219 cell death	0.0007	0.0116
GO:0016265 death	0.0008	0.0122
GO:0006793 phosphorus metabolic process	0.0008	0.0131
GO:0006796 phosphate metabolic process	0.0008	0.0131
GO:0007242 intracellular signaling cascade	0.0010	0.0158
GO:0045321 leukocyte activation	0.0012	0.0191
GO:0044092 negative regulation of molecular function	0.0013	0.0198
GO:0010557 positive regulation of macromolecule biosynthetic process	0.0013	0.0206
GO:0045859 regulation of protein kinase activity	0.0015	0.0238
GO:0043549 regulation of kinase activity	0.0019	0.0289
GO:0031328 positive regulation of cellular biosynthetic process	0.0019	0.0290
GO:0009891 positive regulation of biosynthetic process	0.0021	0.0322
GO:0051338 regulation of transferase activity	0.0024	0.0363
GO:0040012 regulation of locomotion	0.0026	0.0397
GO:0019220 regulation of phosphate metabolic process	0.0026	0.0402
GO:0051174 regulation of phosphorus metabolic process	0.0026	0.0402
GO:0051270 regulation of cell motion	0.0026	0.0406
GO:0001775 cell activation	0.0029	0.0446
GO:0002521 leukocyte differentiation	0.0040	0.0609
GO:0000188 inactivation of MAPK activity	0.0043	0.0655
GO:0043405 regulation of MAP kinase activity	0.0052	0.0784
GO:0045449 regulation of transcription	0.0057	0.0849
GO:0006366 transcription from RNA polymerase II promoter	0.0060	0.0896
GO:0051252 regulation of RNA metabolic process	0.0061	0.0913
GO:0030097 hemopoiesis	0.0062	0.0927
GO:0042113 B cell activation	0.0063	0.0940
GO:0045941 positive regulation of transcription	0.0065	0.0968
(B)
GO:0042127 regulation of cell proliferation	0.0001	0.0015
GO:0034101 erythrocyte homeostasis	0.0001	0.0015
GO:0007169 transmembrane receptor protein tyrosine kinase signaling pathway	0.0008	0.0122
GO:0048872 homeostasis of number of cells	0.0014	0.0227
GO:0007243 protein kinase cascade	0.0021	0.0333
GO:0048514 blood vessel morphogenesis	0.0036	0.0569
GO:0008284 positive regulation of cell proliferation	0.0039	0.0614
GO:0001666 response to hypoxia	0.0041	0.0643
GO:0070482 response to oxygen levels	0.0049	0.0766
GO:0001525 angiogenesis	0.0059	0.0902
GO:0007167 enzyme-linked receptor protein signaling pathway	0.0063	0.0973

Signature probe sets whose expression changes after drug treatment were significantly different between the RAF/MEK/ERK inhibitors (or PI3K/AKT/mTOR inhibitors) and other agents in HT29 cells were extracted using the Student’s t-test (fold-change values of more than 2 and the P-value of less than 0.05). We carried out GO analyses using the DAVID analytical tool. FDR, false discovery rate. Characteristic GOs for each signature were indicated as bold letters (phosphate metabolic process-related GOs for the RAF?MEK?ERK inhibitors and the GOs related to erythrocyte homeostasis,response to hypoxia, and angiogenesis for the PI3K ?AKT?mTOR inhibitors).

Gene signature analyses revealed potential new target pathways of some anticancer drugs

As described above, the anticancer drugs were basically clustered in a target pathway-dependent manner. However, we also found several agents that were clustered in unexpected positions. As shown in Figure1, CDK4 inhibitor, AG1024 (IGF1R inhibitor), and FH535 (β-catenin/TCF inhibitor) unexpectedly showed similar gene expression signatures with the ER stress inducers. Amrubicin is an anthracycline drug that is supposed to target DNA topoisomerase II.24 However, the agent was not clustered together with other topoisomerase II inhibitors but instead with the proteasome inhibitors (Fig.1). These data suggest potential novel modes of action for these agents. Among these drugs with unexpected gene signatures, we focused on vismodegib, a Hedgehog pathway inhibitor,25 because our clustering analysis suggested its possible similarity with the oncogenic kinase inhibitors (Fig.1). To validate whether vismodegib could affect kinase signaling pathways, we examined its effect on the phosphorylation of components in the MEK/ERK and AKT/mTOR pathways. As shown in Figure4, vismodegib clearly suppressed the phosphorylation of p70S6K, a molecule downstream of mTOR, in HT-29 cells as well as in PC3 cells in which the AKT/mTOR pathways are strongly activated. As a positive control, we also observed inhibition of p70S6K phosphorylation by temsirolimus, a clinically used mTOR inhibitor. In contrast, ERK and AKT phosphorylation was not significantly affected by vismodegib treatment, although we observed a marginal inhibition of ERK phosphorylation in PC3 cells (Fig.4). These data indicated that our gene signature analysis successfully revealed a novel action of vismodegib on the mTOR pathway.

Fig 4

Effect of vismodegib on the ERK and protein kinase B (AKT)/mammalian target of rapamycin (mTOR) signaling pathways. HT-29 and PC3 cells were treated with vismodegib or temsirolimus at the indicated concentrations for 2 h. The phosphorylation and expression of ERK, AKT, and p70S6 kinase were analyzed by Western blotting. Actin expression was also examined as a loading control.

Discussion

During anticancer drug development, the molecular target of each candidate compound should be strictly determined with reliable methods. In the present study, we showed that gene signature-based analysis can classify oncogenic pathway inhibitors in a target pathway-dependent manner and is a powerful tool to evaluate the molecular targets of compounds. We prepared subsets of genes for the signature analysis and showed that the signature reflected the modes of action of the agents (Figs3). These data indicated that the analysis worked well to validate target pathways of the agents.

Overall, most inhibitors of oncogenic kinase pathways formed a unique cluster in the hierarchical clustering analysis as well as in the principal component analysis. These data suggested that this signature-based analysis could predict the potential of compounds to affect oncogenic signaling pathways. Our analysis further revealed that, of the kinome-targeted drugs, agents with similar molecular targets showed similar gene expression signatures. These data indicate the gene signature analysis is effective in validating target molecules or pathways of kinome-targeted compounds.

In addition to the kinase inhibitors, other compounds that target similar molecular pathways were also clustered together. For instance, NAE is a component of the NEDD8 conjugation pathway that regulates protein turnover upstream of the proteasome.20 MLN4924, a specific inhibitor of NAE, was clustered with proteasome inhibitors (Fig.1). Aurora kinase and Polo-like kinase are involved in the process of mitosis.18 Inhibitors of these molecules (Aurora kinase inhibitor III and BI2536) were clustered with tubulin-binding agents, the classical inhibitors of mitosis. These results clearly indicate that the gene expression signatures reflect the primary target pathways of the agents.

However, we observed some kinase inhibitors that were clustered in an unexpected way in our signature analysis, such as ALK inhibitors that were not clustered with the majority of the oncogenic kinase inhibitors (Fig.1). These data potentially suggest that these agents could affect unique downstream pathways; however, we should also take into account off-target effects, because in this study we used these agents at higher doses than the clinically relevant concentrations.

Our gene expression analysis also assigned some antitumor drugs to unexpected modes of action. One agent was vismodegib, an inhibitor of the Hedgehog pathway,25 whose gene expression pattern showed significant similarity with those of the oncogenic kinase inhibitors. Our “wet” experiments confirmed that vismodegib actually inhibits the mTOR pathway. These data indicate that the signature-based analysis was effective in identifying novel target pathways of the drugs.

Endoplasmic reticulum stress is involved in the mode of action of some anticancer drugs.26,27 In this study, celecoxib, a selective inhibitor of COX2, showed a similar gene signature to that of ER stress inducers (Fig.1). This result is consistent with previous reports showing that the cytotoxic effect of celecoxib correlates with ER stress.7,28 Additionally, we also found several agents that were clustered together with ER stress inducers, such as a CDK4 inhibitor, AG1024 (IGF1R inhibitor), and FH535 (β-catenin/TCF inhibitor). These results suggest that these agents could affect ER stress pathways.

In the signature-based analysis, careful interpretation of results was required. First, we needed to administer relatively high doses of agents that were less cytotoxic to cancer cells. As for the EGFR inhibitors, gefitinib and erlotinib, we tested how the drug concentration would affect the result of the signature analysis and found that high-dose treatment (30 μM) still showed significant similarity in gene signature to that of low-dose treatment (0.6 μM) (Table3D). Nevertheless, for such high-dose treatment data, we should be careful to confirm whether the signature reflects the physiological mode of action of the agents. Second, as we mentioned above, gene signatures of the agents could depend on cell context in some cases. As we have shown, the signatures of oncogenic kinase inhibitors in different cancer cell lines showed significant similarity (Fig.1). However, the target pathway-based classification was more accurately achieved using the data of a single cancer cell line (Fig.3). These data would be valuable to examine the cell context effect on the signature analysis. We further showed that, for the agents whose targets are selectively expressed in a certain subtype of cancer, use of data obtained in specific cancer cells could help accurate evaluation of drug target pathways (Table3D). In this aspect, our data would be valuable because we obtained the gene signature data using multiple specific cancer cell lines. Finally, the signature analysis could reveal target “pathways” of each agent, but the analysis would not be enough to completely define target “molecules” of the agent (for example, inhibitors of mitotic pathways showed similar gene signatures despite the direct target of each agent being different). Considering these points, integrated approaches with signature analysis and other methods would be important for accurate evaluation of the molecular targets of antitumor compounds.

There are several other publicly available databases related to compounds’ transcriptomic data. Connectivity map (C-map) (https://www.broadinstitute.org/cmap/) is a pioneering database that contains genome-wide transcriptome data for more than 1000 compounds.9 In addition, several other databases containing drug-related gene expression data have recently been established, such as the Library of Integrated Cellular Signatures (http://www.lincsproject.org/) and the Cancer Cell Line Encyclopedia (https://www.broadinstitute.org/ccle/home). These are huge databases, but they do not focus on anticancer drugs, nor do they cover all antitumor agents. Our database is unique in that it is a compact database focusing on anticancer drugs and it covers genome-wide gene expression data of most clinically available anticancer compounds as well as promising inhibitors of molecular cancer targets. Moreover, we are updating the database by adding newly approved agents’ data. Our website (http://scads.jfcr.or.jp/db/cs/) also provides an online analysis tool for users to easily compare the gene signature of query compounds to those in our database. These aspects make our database more updated and user-friendly, particularly for oncologists, than other public databases providing gene expression data. It should also be noted that our data were obtained using HT-29 cells as well as the specific driver oncogene-expressing cell lines, whereas the C-map and the other databases used different types of cells. Therefore, we believe that the combination of our database and others would provide more robust information to estimate modes of action of anticancer compounds.

In conclusion, we obtained and analyzed gene expression data for a wide variety of molecularly targeted agents. This is a unique, comprehensive analysis of gene expression related to the pathways of molecularly targeted anticancer drugs. Our data will not only be beneficial in classifying antitumor agents but could also be valuable as a reference database to evaluate the modes of action of new candidate compounds in drug development.