From GWAS to drug screening: repurposing antipsychotics for glioblastoma

Glioblastoma is currently an incurable cancer. Genome-wide association studies have demonstrated that 41 genetic variants are associated with glioblastoma and may provide an option for drug development. We investigated FDA-approved antipsychotics for their potential treatment of glioblastoma based on genome-wide association studies data using a ‘pathway/gene-set analysis’ approach. The in-silico screening led to the discovery of 12 candidate drugs. DepMap portal revealed that 42 glioma cell lines show higher sensitivities to 12 candidate drugs than to Temozolomide, the current standard treatment for glioblastoma. In particular, cell lines showed significantly higher sensitivities to Norcyclobenzaprine and Protriptyline which were predicted to bind targets to disrupt a certain molecular function such as DNA repair, response to hormones, or DNA-templated transcription, and may lead to an effect on survival-related pathways including cell cycle arrest, response to ER stress, glucose transport, and regulation of autophagy. However, it is recommended that their mechanism of action and efficacy are further determined.


Background
Glioblastoma multiforme (GBM) classified by the WHO as a grade IV glioma is the most malignant and lethal tumor to occur in the brain with rapid de novo progression and limited survival rate [1,2]. The usual conventional therapies for GBM are surgical resection, radiation therapy, and chemotherapy. Concurrent radio-and chemotherapy will usually be enforced after surgery in case the cancer cells have invaded distal tissue or infiltrated the parenchyma. Currently, Temozolomide (TMZ) is the standard drug of choice used for the primary management against GBM. TMZ is a prodrug that is spontaneously hydrolyzed into an active form, 3-methyl-(triazen-1-yl) imidazole-4-carboxamide which alkylates DNA and triggers the death of tumor cells [3]. However, the median survival time (overall survival) of GBM patients is only extended a few months with TMZ treatment [4], and resistance to TMZ usually quickly develops through DNA repair mediated by methyl guanine methyl transferase, base excision repair or alkylpurine-DNA-N-glycosylase [3,5]. Although an FDA-approved anti-angiogenic therapy (bevacizumab) with a combination of TMZ has been applied in GBM management, it did not result in significant improvement in overall survival but only prolonged progressionfree survival [6,7]. Therefore, a therapeutic approach for GBM still falls far short of medical needs. Genome-wide association studies (GWASs) investigate the association between genetic variants and traits of interest, for example, the association between single nucleotide polymorphisms (SNPs) and diseases [8,9]. In rare cases, SNPs located in the promoter or the wellknown regulatory elements can be easily correlated to the regulation of genes. In most of cases (~ 88%), the SNPs of interest are located in intergenic or intronic regions, and are often associated with traits of interest with unknown reason [10]. Currently, four GWASs have reported that 41 risky SNPs are associated with GBM [11][12][13][14], only three of them are located in the untranslated regions or miRNA which is known to be implicated in oncogenesis [15][16][17]. Rs78378222 disrupts the polyadenylation of TP53 resulting in impaired mRNA processing and its reduced expression [15]. Rs10069690, providing an alternative splicing site and leading to a splice isoform of Telomerase reverse transcriptase [16], was found to have the most significant P value associated with GBM. Rs11558961 affects the secondary structure of Glial fibrillary acidic protein mRNA, which promotes the binding of miR-139, and thus decreases the susceptibility to chemotherapy [17]. The roles of other SNPs in GBM remain unclear to date.
Drug discovery and development are an expensive and highly time-consuming undertaking. It usually costs billions of dollars and takes decades to bring a compound to clinical application. Furthermore, newly developed drugs do not always confer a superior clinical benefit [18]. Scientists and pharmaceutical industries face growing demand for new drugs to fit unmet medical needs, and are under pressure to increase R&D productivity with limited resources [19]. Drug repositioning (or drug repurposing) is a process that resort to the approved drugs for new indications or answer of unmet medical needs. Clinically approved drugs have been proven to be safe, and the dosage range and formulation are already well studied. Taking advantage of abundant information about clinically approved drugs, drug development by drug repositioning may have a lower risk of failure, cost less and require less time to complete preclinical and phase I/II clinical trials [20,21]. Furthermore, the efficacy for unidentified targets of approved drugs within other diseases remains worthwhile pursuit.
Over the past decades, GWASs have uncovered a lot of genetic variants which may provide targets for drug repositioning. Lau and So summarized several approaches for drug repositioning using GWAS data [22]. The 'candidate gene approach' is the most straightforward method that maps the top risky SNPs to their related genes by functional annotation tools or eQTL, and then queries these genes in the gene-drug database. This approach is a clear-cut strategy with low computational cost, but meets difficulties such as limited druggable genes [23], uncertain effect size of risky SNPs [24], and the complexity of proper annotation of risky SNPs. This approach may also miss multi-target drugs which are considered to be more effective than single target drugs [25]. The 'pathway/gene-set analysis' organizes multiple function or biologically-related genes and further studies their perturbational expression under interruption of drugs. Although the 'pathway/gene-set analysis' seems to be able to overcome the defects of the 'candidate gene approach' , its challenge lies in the complexity of identifying the drug-mediated pathways or gene sets, because the mechanisms of many drugs are not fully understood. The computational tools Gene2drug [26] and Drug Set Enrichment Analysis (DSEA) [27] provide an opportunity to overcome the complexity of identifying drug-mediated pathways. Gene2drug analyses the genes perturbed by drug treatments and annotates those genes to biological/functional pathways; DSEA reversely annotates queried drugs to the pathways.
Previous review indicate that patients receiving drugs acting on the brain to sedate psychosis have a lower incidence of cancers than the general population [28]. Thus, antipsychotics have been considered as potential candidate agents against GBM [29][30][31]. This work aimed to discover whether there are current drugs that may be beneficial for GBM therapy using GWAS data via 'pathway/gene-set analysis' . Gene2drug, Gene Ontology (GO) Resource [32] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [33] were used to assist in in-silico screening. The antitumor activities of the candidate drugs and TMZ were extracted from DepMap for which data were identified by PRISM viability assay on 42 glioma cell lines [34].

Data source
The risky SNPs of GBM (Trait NO.: EFO_0000519) and the annotated risk genes were collected from the GWAS Catalog, EMBL-EBI [35]. Gene2drug (https:// gene2 drug. tigem. it/) is a state-of-the-art online tool and easily accessible technique applying the Gene Expression Profiles to Pathway Expression Profiles (gep-2pep) algorithm [26] with Connectivity Map (CMap) database Molecular Signatures Database (MSigDB). The gene ontology of the annotated genes recoded by Gene Ontology Resource (GO) [32] or KEGG [33] were selected. The sensitivity of glioma cells to candidate drugs or TMZ were sorted from PRISM Repurposing Secondary Screen 19Q4 [34]. The gene mutation profiles (Mutation Public 21Q1) were extracted from DepMap portal (https:// depmap. org/ portal/). The gene expression profiles were downloaded from Cancer Cell Line Encyclopedia (CCLE_expression 21Q1) [36].

Gene to drug analysis
Candidate compounds were suggested by Gene2drug with the imported biological pathways in which the annotated genes were involved [26]. Briefly, the biological pathways of query genes were selected. A list of compounds will be generated by each pathway which suggests that the compounds in the list may interrupt the pathway based on data from MSigDB. Then, a merge process which selected compounds exist in all the lists and are ranked by P-value decides the output of Gene2drug. Nine annotated genes were found to be involved in a certain biological process. The significant compounds which had a P value lower than 1E-2 were included. The compounds co-existing in more than four lists were selected (six lists at most with candidate drugs included, a median number, four, was set as cut-off ), and then chosen within FDA approved drugs with known central nervous activities based on the treatment note of PRISM assay. The co-existence of compounds among the lists was analyzed with tool Venn Diagrams (Van der Peer Lab) and outputted in text format.

Prediction of pathways in which the candidate drugs are involved
The pathways which may be mediated by candidate drugs were predicted by using DSEA [27]. The downstream pathway shared by candidate drugs and the pathway based similarity were generated at the same time. The enrichment scores (E-score) from −1 to + 1 refer to down-regulation or up-regulation by drug treatments.

Validating the antitumor activity of candidate drugs
The sensitivity of candidate drugs and TMZ was extracted from DepMap for which data were identified by PRISM viability assay on more than 40 glioma cell lines. Briefly, cell lines with a barcode integrated in the genome were exposed to compounds for five days. Then, the mRNA was isolated from cells, amplified and detected. The viability of cells was calculated based on the level of barcode abundance and the sensitivity of cells to compounds was generated by comparing the treatment group and the control group [34]. The process sorted data from PRISM was performed by a python script (https:// github. com/ LinWZ-tw/ g2d-prism).

Gene expression profile
The expression profiles of cell lines were further analyzed with data downloaded from the CCLE via DepMap portal. The cluster and heatmaps of genes were generated by using ClustVis [37].

Prediction of drug target
The targets of candidate drugs were predicted by a webbased platform, GalaxySagittarius, which combines ligand similarity-based and receptor structure-based approaches with AUC up to 0.8 [38]. Briefly, a protein target database is generated by similarity-based screening, then the target protein is ranked by protein-ligand docking. The predictions were run with a search model of similarity combined prediction and 3D structure compound obtained from ZINC: ZINC000002040609 (Norcyclobenzaprine) and ZINC000001530764 (Protriptyline). The 3D binding poses and binding structures were downloaded from the web server; the 2D binding structures were generated by using BIOVIA Discovery Studio 2020. The gene ontologies of interest were analyzed by GO Ontology Resource (http:// geneo ntolo gy. org/) [32].

Limitations
For in-silico screening, we collected GBM-associated SNPs and their mapping genes from the GWAS Catalog which applies Ensembl mapping pipeline for genome annotation. Although the risky SNPs are likely associated with their mapping gene, the possibility that risky SNPs have an effect on distal genes via chromatin looping structures may be omitted [10]. The in-silico screening conducted by Gene2drug was limited by the size of MSigDB while it is one of the largest and most popular databases so far [26,39]. Validating work was relied on PRISM data which may contain divergent results in replicate tests.

In-silico screening of potential drugs
Four studies from the GWAS Catalog were used [11][12][13][14], which include 41 risky SNPs and 15 annotated genes (Additional file 1: Table S1). Only 9 of the annotated genes were found to be involved in biological processes or pathways collected in GO or KEGG (Additional file 1: Table S2). Then, lists of candidate compounds from each annotated gene were generated from Gene2drug [26]. Further analysis of compounds showed that one compound co-existed in six lists, eight compounds co-existed in five lists and 20 compounds co-existed in four lists (Table S3), no compounds could be found in more than six lists, and compounds found in fewer than four lists were discarded. Given that the aim was drug repositioning for GBM therapy, compounds had been approved by FDA for safety, and compounds known to affect the central nervous system were selected. A total 12 FDA-approved drugs were suggested (Table 1, Additional file 1: Figure S1). The protocol design and overall flowchart is shown (Fig. 1).

Prediction of pathways in which the candidate drugs are involved
The downstream pathways shared by candidate drugs were predicted by DSEA (see Methods). A list of 327 downstream pathways shared by 12 candidate drugs were collected (Additional file 1: Table S4). The much commonly shared pathways of significance (lowest P value) were ranked by DSEA utility (Additional file 1: Figure  S2A). Among them, several categories of pathway which may lead to cell death are summarized, including cell cycle arrest, response to ER stress, glucose transport, and regulation of autophagy ( Table 2).

The sensitivities of glioma cells to candidate drugs
Among the candidate drugs output from screening protocol, Chlorprothixene was excluded because two datasets of Chlorprothixene were found in which glioma cells show disparate responses. The sensitivities of 42 glioma cells to the rest 12 candidate drugs were extracted from PRISM (Additional file 1: Table S5). The sensitivities to all the candidate drugs were higher than that to TMZ based on mean sensitivity among cell types ( Fig. 2A, B). In particular, Norcyclobenzaprine and Protriptyline were interesting in that their mean sensitivities were significantly higher than that to TMZ (P < 1 E−4) ( Fig. 2A,  B). The sensitivities to Protriptyline have no significant correlation with well-known oncological genes such as expression level of EGFR or mutations in TP53 and PTEN (Additional file 1: Figure S3A-C and G). The sensitivities to Norcyclobenzaprine were correlated with mutations in TP53 (0.017 by Mann-Whitney U test) but were not correlated with PTEN or expression level of EGFR (Additional file 1: Figure S3 D-F and G). It suggests that these two candidate drugs may have potential application on tumors with a wide spectrum of genetic backgrounds, while the mutation in TP53 should be avoid when Norcyclobenzaprine is applied.
To investigate the difference between sensitive and not sensitive cells, genes showing a significantly (P < 1E−3) different level of expression were selected. A group of genes with a consistent pattern in cells not sensitive to Norcyclobenzaprine was disclosed (Fig. 2C). Another group of genes showed a similar pattern in cells not sensitive to Protriptyline (Fig. 2D). Subjects with similar expression level in those groups of gene may not respond to Norcyclobenzaprine and Protriptyline.

Prediction of drug target
Several proteins were predicted by GalaxySagittarius to be targets of candidate drugs, the top ten from each prediction were selected (Additional file 1: Table S6). The top ten for Norcyclobenzaprine were: Poly(ADP-Ribose)   Figure S4). The targets of both Norcyclobenzaprine and Protriptyline have similar gene ontology which are correlated with DNA-templated transcription, RNA synthesis, response to hormone, retinoic acid biosynthetic process, etc. (Table 3).

Direct screening protocol based on PRISM assay data
We built another screening protocol which directly ranks drugs according to PRISM assay data (Additional file 1: Figure S5A). In the field of neurology, there are 124 FDAapproved medications. Replicate tests were found in some of drugs result in a total of 140 treatments including TMZ as a control were employed. Glioma cells show higher mean sensitivities to 107 treatments (76.98%) than that to TMZ, but only 23 (16.55%) out of them show significantly higher sensitivities than TMZ (P < 1E−3). Among the candidate drugs with sensitivities higher than TMZ, the standard deviations were ranging from 0.255 to 0.805 (Additional file 1: Figure S5B). The top ten candidate drugs were Amitriptyline, Tranylcypromine, Mianserin, Triflupromazine, Doxylamine, Protriptyline, Metixene, Citalopram, Benserazide and Acepromazine (Additional file 1: Figure S5C-D). Moreover, nine of top ten candidate drugs were predicted to bind to PGR and eight were predicted to bind to Androgen Receptor (AR) (Additional file 1: Figure S5E). Both of PGR and AR have been reported to play important roles in glioma [40][41][42].

Discussion
A total 12 FDA-approved drugs were suggested as drugs of possible interest for GBM in our in-silico screening: Norcyclobenzaprine, Protriptyline, Iobenguane, Haloperidol, Alimemazine, Nortriptyline, Melatonin, Trifluoperazine, Perphenazine, Spiperone, Imipramine and Levomepromazine. Protriptyline, Nortriptyline and Imipramine are used for treatment of depression and classified as tricyclic antidepressants (TCAs). Norcyclobenzaprine is one of the major metabolites of Cyclobenzaprine which is usually used for muscle spasms, while both Norcyclobenzaprine and Cyclobenzaprine can act as TCAs to block serotonin receptors. A previous case-control study indicated that long-term use of TCAs may be associated with reduced risk of glioma [43]. TCAs are known to down-regulate β-adrenergic receptors which are involved in carcinogenesis and are considered as drug targets [44]. Melatonin is commonly used for sleep disorders. Haloperidol, Trifluoperazine, Perphenazine, Spiperone and Levomepromazine are clinically used for psychosis, especially schizophrenia or bipolar disorder. Alimemazine could be used as a sedative. Perphenazine and Levomepromazine are also used to control nausea and vomiting. Alimemazine, Trifluoperazine, Perphenazine and Levomepromazine are structurally similar compounds based on a phenothiazine. One of the candidate drugs, Iobenguane, is already an FDA-approved drug for low-grade glioma, indicating the credibility of our in-silico screening protocol. Among these candidate drugs, Norcyclobenzaprine, Iobenguane, Haloperidol, Melatonin, Trifluoperazine, Perphenazine, Spiperone and Imipramine have been studied for their effects on brain tumors. For example, Norcyclobenzaprine has been reported to inhibit the proliferation of GBM cells but with IC 50 values higher than 10 µM [45]. Iobenguane has been developed and approved by the FDA as a clinical agent for treating pheochromocytoma and paraganglioma [46]. The candidate drugs identified in our study have not only been studied for use on GBM but also many types of human malignancies  and many clinical trials have been completed [76][77][78] or are ongoing. Especially, there are 30 completed trials and 15 ongoing trials with intervention with Melatonin. Even though the antitumor activities of the candidate drugs have been reported, most studies exclusively report the antitumor activity of antipsychotic on a few representative cell lines. PRISM provided an overview of the sensitivity among 42 glioma cell lines to queried drugs. The cell lines which are outliers in sensitivity could be further studied as a foundation for selection of participants in future clinical trials. For example, Norcyclobenzaprine did not show significant antitumor activity on certain cell lines [45]; however, the overview of sensitivity on 42 glioma cell lines indicates promising potential. Protriptyline, Alimemazine, Nortriptyline and Levomepromazine have been studied for their effects on many types of human malignancies, but have not yet been reported to have been tested on brain cancers ( Table 1).
Most of candidate drugs are similar on correlation matrix according to GO biological process except iobenguane and melatonin (Additional file 1: Figure S2B). But the candidate drugs are not very similar on correlation matrix according to sensitivity (Additional file 1: Figure  S2C). Iobenguane, alimemazine, perphenazine, spiperone and imipramine show significant difference when comparing the correlation matrix according to GO biological processes and that according to sensitivities on glioma cells (Additional file 1: Figure S2D).
We further predicted the targets of candidate drugs by GalaxySagittarius. The binding poses of candidate drugs and predicted proteins show that the tricyclic structure  of Norcyclobenzaprine and Protriptyline contributes to π-interactions in the inner part of the binding pocket and the amine group at the 'tail' structure contributes to H-bonds with residues close to the mouth of the binding pocket (Additional file 1: Figure S4). Compounds with similar structures were found to have anti-tumor activity. For example, Triflupromazine which has a similar structure to trifluoromethyl-phenothiazine in the 'head' and a tertiary amine in the 'tail' shows potential against GBM (Additional file 1: Figure S5C). While most of the pathways participated in by each risk gene did not overlap (Additional file 1: Table S2), the candidate compounds suggested by Gene2drug co-existed in multiple lists predicted from different genes (Additional file 1: Table S3). It is common that a compound may have multiple off-target mechanisms. We hypothesize that the candidate compounds which disrupt several pathways associated with GBM may have greater potential than those with a single pathway associated with GBM. We further merged the compounds suggested by each input of a gene through ranking the candidate compounds by the number of lists in which they co-existed. Then, we selected a median number as cut-off criteria, and the candidate drugs co-existing in at least four lists of pathways were employed in further analysis. Using TMZ as a control, the sensitivity of cells to all the candidate drugs selected by this protocol was higher than that to TMZ. On the other hand, the top ten candidate drugs which were directly ranked according to mean sensitivity have lower P-value than the 12 candidate drugs correlated with multiple disease-risk genes (Fig. 2B, Additional file 1: S5D). We assume that a drug which only interrupts a few key survival-related genes still may trigger cell death and could lead to apparent antineoplastic outcome. For example, Triflupromazine was predicted to disrupt pathways which only correlate with TP53, but Triflupromazine shows an antineoplastic activity significantly higher than TMZ (P = 2.63E−05, Additional file 1: Figure S5D). While screening protocol directly based on PRISM cancer cell sensitivity data without prior search of inherent biological pathway and risk genes resulted in a parallel consistent list of candidate drugs, though more, it provided less information for further drug development other than sensitivity. Our current development regimen not only aims to work on cancer disease but also for other diseases of which potential drug targets and pathways are also involved. Furthermore, for some indications for which cell-or animal-models are difficult to construct, the approach that selects candidates correlated with several disease-related genes may reduce the number of candidates to a small figure, reducing the resources required for bench work.
Based the target prediction of candidate drugs by GalaxySagittarius and disrupted pathway predicted by DSEA, we suspect that Norcyclobenzaprine and Protriptyline may bind to targets of the "shell", such as PARP1, PARP2 or PGR, to interrupt a certain molecular function, for example, DNA-template transcription, response to hormone stimulation, etc. Then the influence of treatment is further extended to survival-related "core" pathways including cell cycle arrest, response to ER stress, glucose transport, and regulation of autophagy. However, the link between "shell" and "core" remains to be further investigated.

Conclusion
The current study presents a screening protocol following selection of candidates using Gene2drug. Selection from candidate compounds which correlate with multiple disease-risk genes or variants may reduce the number of candidates and decrease the burden of bench work to validate their therapeutic efficacy.
Our in-silico screening led to ten antipsychotics which show anti-tumor activity which is higher than TMZ in 42 cell lines. In particular, Norcyclobenzaprine and Protriptyline show significant potential against GBM; they are predicted to bind targets such as PARP1, PARP2, PRG, RBP1 to disrupt DNA repair pathways, respond to hormone and DNA-templated transcription and the retinoic acid signaling pathway, further effect survival-related pathways including cell cycle arrest, response to ER stress, glucose transport, and regulation of autophagy. Of these, the activity of Protriptyline against GBM has not yet been reported. The mechanism of action and therapeutic efficacy of Norcyclobenzaprine and Protriptyline are worth pursuing further.