Uncovering drug repurposing candidates for head and neck cancers: insights from systematic pharmacogenomics data analysis

Effective treatment options for head and neck squamous cell carcinoma (HNSCC) are currently lacking. We exploited the drug response and genomic data of the 28 HNSCC cell lines, screened with 4,518 compounds, from the PRISM repurposing dataset to uncover repurposing drug candidates for HNSCC. A total of 886 active compounds, comprising of 418 targeted cancer, 404 non-oncology, and 64 chemotherapy compounds were identified for HNSCC. Top classes of mechanism of action amongst targeted cancer compounds included PI3K/AKT/MTOR, EGFR, and HDAC inhibitors. We have shortlisted 36 compounds with enriched killing activities for repurposing in HNSCC. The integrative analysis confirmed that the average expression of EGFR ligands (AREG, EREG, HBEGF, TGFA, and EPGN) is associated with osimertinib sensitivity. Novel putative biomarkers of response including those involved in immune signalling and cell cycle were found to be associated with sensitivity and resistance to MEK inhibitors respectively. We have also developed an RShiny webpage facilitating interactive visualization to fuel further hypothesis generation for drug repurposing in HNSCC. Our study provides a rich reference database of HNSCC drug sensitivity profiles, affording an opportunity to explore potential biomarkers of response in prioritized drug candidates. Our approach could also reveal insights for drug repurposing in other cancers.

Hierarchical clustering and heatmap visualization of the data matrix of active compounds using Morpheus. A data matrix of 418 active targeted cancer compounds, 404 active non-oncology compounds and 64 active chemotherapy compounds for the 28 HNSCC cell lines were uploaded onto the Morpheus tool, an online platform for data analysis and visualization (https:// softw are. broad insti tute. org/ morph eus/). To determine if there is any pattern of drug sensitivity that is associated with the subsite of the HNSCC, annotations for cancer subtype (oral squamous cell carcinoma (OSCC) or non-OSCC) were included. Further, the MOA was annotated for each compound. Hierarchical clustering analysis was performed on Morpheus, using the metric of "one minus spearman correlation" and "complete" linkage method, on both rows (compounds) and columns (cell lines). Interactive heatmaps can be accessed and explored by uploading the json files [Files S1-S3], onto the Morpheus tool [ Fig. S3A]. These files are downloadable from figshare-https:// doi. org/ 10. 6084/ m9. figsh are. 16443 759. v1.

t-distributed stochastic neighbour embedding (tsne) analysis and interactive visualization
on RShiny app. The data matrix of 886 active compounds in 28 HNSCC cell lines consist of 24,808 data points, and the Morpheus tool was utilized to perform t-distributed stochastic neighbour embedding (tsne) analysis [Metric: one minus spearman rank correlation; Run on: Rows (compounds); Epsilon (learning rate): 100, perplexity : 30]. Other combinations of the epsilon and perplexity were also tested (300/30, 500/30), where compounds including EGFR inhibitors, MEK inhibitors (MEKi) were shown to cluster together 7 , were also recapitulated in our analysis. Representative T1 and T2 values from one of the ten runs were used for tsne plotting. Non-oncology compounds that consistently fall into proximity with (or consistently clustered close-by) the groups of targeted cancer compounds in all ten runs of tsne analysis were highlighted. To enable interactive visualization, we have created a RShiny app webpage: https:// annie chai. shiny apps. io/ tsne/, where details on each of the compounds can be obtained (Fig. 2D). Further, searches based on compounds or MOA or targets of interest can be conducted where results of the search query would be highlighted as enlarged dot(s) in the tsne plot as shown in Fig. S3B.
Identification of compounds with preferential activity in HNSCC. The "Data explorer" tool from DepMap was utilized to perform the "two-class comparison" analysis (https:// depmap. org/ portal/ inter active/), to compare the means of drug sensitivity data between two groups and generate the estimates of effect size and corresponding p-value. The "drug sensitivity (PRISM repurposing Primary screen 19Q4)" was selected as input dataset, with the 28 HNSCC cell lines being placed in the "in" group cell lines, while all other cell lines tested in PRISM were used as "out" group cell lines. A total of 179 unique compounds (ruxolitinib and calcitriol had duplicated screen results) were found to have preferential activity in HNSCC (P < 0.01 or − log10 (P value > 2) Analysis workflow in identifying potential gene expression-based biomarkers for MEKi. The logFC data for all 20 MEKi were visualized using a heatmap in Morpheus [File S4, downloadable from figsharehttps:// doi. org/ 10. 6084/ m9. figsh are. 16443 759. v1]. Hierarchical clustering was performed using a metric of "one minus spearman correlation" and "complete" linkage method. A distinct cluster of MEKi-resistant cell lines was observed, which included six lines: FADU, SNU1041, BICR18, HSC2, SNU46 and YD8. For the compounds, a few sub-branches consisting of eight active MEKi were found to share largely similar drug sensitivity profiles. This cluster consists of the three MEKi that were differentially enriched for killing activity in HNSCC (trametinib, AZD8330 and cobimetinib), as well as another five compounds (MEK162, Ro-4987655, AS-703026, PD-0325901 and TAK-733). The clinical indication and drug development status of these eight MEKi were manually curated [ Table S6].
Differentially expressed genes (DEGs) were identified using the limma package (Bioconductor), by comparing the six resistant lines versus the six most sensitive lines (BICR6, CAL27, SCC25, PECAPJ49, SNU1076, YD38), based on their average logFC against the eight MEKi. The list of DEGs is available from Table S7. The 136 significantly upregulated DEGs in MEKi-sensitive lines were used as input query into the STRING portal (https:// string-db. org, Version 11) network. Pathway enrichment analysis was performed using the portal and the most enriched REACTOME pathway-"Cytokines signalling in immune system" was highlighted in red (n = 24). The gene expressions of these 24 cytokines were used to correlate with MEKi average logFC, those with significant Pearson's correlation (P-value < 0.05) were shortlisted as potential biomarkers of response to MEKi. GSEA comparison across multiple datasets. GCT file with gene expression matrix of the six MEKisensitive and six MEKi-resistant HNSCC was used as input for GSEA using GenePattern 11 . Gene set database "h.all.v7.4.symbol.gmt [Hallmarks]" was selected, with the rest following default parameters. Similar GSEA was also performed using the two following datasets as described.
For the drug screen by Lepikhova et al. 12 , six HNSCC lines were contained in each MEKi low and high sensitivity groups. High MEKi sensitivity group-UT_SCC_47, UT_SCC_42A, UT_SCC_106A, UT_SCC_40, UT_ SCC_24A; Low MEKi sensitivity group-UT_SCC_54a, UT_SCC_6A, UT_SCC_21, UT_SCC_29, UT_SCC_28, UT_SCC_44. The microarray gene expression dataset of these 12 cell lines was downloaded from Gene Expression Omnibus (GEO) with the accession number GSE108062. GSEA was performed as described above, with Illumina HumanHT-12 V4.0 expression bead chip selected as the chip platform.
Hallmark gene sets enriched in MEKi-resistant groups of respective datasets (false discovery rate < 25%) were compared for similarity. Hallmarks commonly identified from all three datasets were highlighted in the Venn diagram. Full results of the GSEA are available in Table S8.
Statistical analysis and Graph Pad Prism. The differences in the mean of logFC between the two groups were evaluated for their statistical difference using the two-tailed, unpaired Mann-Whitney test in the GraphPad Prism software 8.0.2. Pearson correlation between gene expression and drug sensitivity was also computed using GraphPad Prism, P-value < 0.05 is considered as statistically significant.

Results
Systematic analysis of high-throughput drug screen data on the HNSCC subset. Several highthroughput drug screens had include HNSCC cell lines, including the CCLE 2 , CTRP 13 , GDSC 4,14 and PRISM 7 . Additionally, other independent research groups have also published their HNSCC screen outcomes, including Ghasemi et al. 15 and Lephikova et al. 12 . These drug screens vary diversely in terms of screening techniques, platforms, cell lines and compounds. Amongst these, the "PRISM Repurposing dataset" conducted by the Broad Institute 7 , represents the largest resource for drug repurposing, and is the only dataset that includes non-oncology compounds published to date [ Fig. 1A]. The primary screen of PRISM tested 4,518 compounds, including targeted cancer compounds (n = 956, 21%), chemotherapy (n = 96, 2%) and a vast variety of non-oncology  Fig. 1C, with limited selectivity. Of the 3,062 non-oncology compounds, 404 (12%) were found to be active in HNSCC lines. The top MOAs included glucocorticoid receptor agonist (n = 13), cyclooxygenase inhibitor (n = 11) and adrenergic receptor agonist (n = 10) [ Fig. 1D]. In the group of glucocorticoid receptor agonist, steroids such as hydrocortisone, prednisolone and methylprednisolone are among those active compounds. These non-oncology compounds had been widely used as palliative/supportive treatment or anti-inflammation agent for cancer patients 16 , and the current data suggest that they may also pose an anti-cancer effect that warrants further investigations. More than half (n = 64, 66%) of the screened chemotherapy drugs are active in HNSCC, with the topoisomerase inhibitors (n = 20), tubulin polymerization inhibitor (n = 13) being the top MOAs [ Fig. 1E]. The proportion of drugs being active in > 90% of HNSCC lines are also highest for these two MOAs, since they target essential cellular processes and are broadly toxic across all lines. Topoisomerase inhibitors that belong to anthracyclines drug such as doxorubicin, epirubicin, daunorubicin, nemorubicin, idarubicin all showed active killing in all tested lines. By contrast, the platinum-based chemotherapy including cisplatin and carboplatin that are widely used for HNSCC appears not as potent in vitro.

Reference database and visualization of drug response profile of HNSCC
To investigate the relationship between the drug response profile and cancer subtype (OSCC or non-OSCC) of HNSCC cell lines, we separated the PRISM screening compounds into three categories-targeted cancer compounds [ Fig Fig. S2]. For targeted compounds, the variability in response is seen across the majority of the drugs, suggestive of response based on genotype/ targets status of each cell line. As expected, the non-oncology compounds are generally less potent across all HNSCC, but instead, their activities are rather selective in a small subset of cell lines [ Fig. 2B]. On the contrary, most of the chemotherapy compounds are rather unselective and show potency across many of the HNSCC cell lines. Notably, the BICR18 line with a mismatch repair defect (MSH2 and MLH1-mutated), appears to be sensitive to most of the chemotherapeutic agents [ Fig. 2C].
To enable easier exploration and analysis to the interactive, customizable and searchable data matrix via the online Morpheus tool, we provide the data matrix for the 4,518 compounds × 28 HNSCC cell lines [Table S1] and the json files [Files S1-S3] for users to explore this data via Morpheus. www.nature.com/scientificreports/ Exploring the opportunity for repurposing of non-oncology drugs. Next, we investigated 404 active non-oncology compounds in HNSCC for drug repurposing for this disease. From the heatmap with hierarchical clustering, we noticed that active targeted compounds with clear and similar mechanisms of action often tend to cluster together, suggestive of consistent effect on cell lines and their putative activities were recapitulated in the screen. It has been previously shown that compounds with similar target profile and structure (and presumably similar MOA) clustered closely together when their drug sensitivity profiles are dimensionally reduced using tsne analysis 17 . As such, we performed tsne analysis [ Fig. 2D] for unsupervised clustering of the 886 active compounds to identify any non-oncology compounds that show a similar killing pattern with those targeted compounds of distinct MOAs. Consistent with the Uniform manifold approximation and projection (UMAP) done by Corsello et al. 7 on the pan-cancer cell lines, some of the targeted compounds consistently show strong clustering-this includes the EGFR inhibitors, MEKi, PI3K/AKT/MTOR inhibitors, multi-tyrosine kinase inhibitors (mTKIs) and Aurora kinase inhibitors [ Fig. 2D]. Active compounds that kill > 90% of the HNSCC cell lines tend to cluster together, albeit not tightly and contains the majority of the less selective chemotherapy drugs. Non-oncology compounds were very scattered across the plot, reflecting the vast diversity of MOAs and likely polypharmacology activities. Seven non-oncology compounds were clustered with the EGFR inhibitors (piroxicam, FK-3311, spironolactone, fenofibrate, palovarotene, levocetirizine) [ Fig. 2E]. Among these, spironolactone and fenofibrate were also found in the cluster of EGFR inhibitors in Corsello's pan-cancer UMAP 7 Table S2]. To strengthen our analysis, we also analyzed the PRISM secondary screen data. Consistently, the top three most significant compounds are also found in the secondary screen, together with 86 other compounds [ Fig. 3B] [ Table S3]. In total, 36 compounds are consistently found to be selectively effective in HNSCC in both the primary and secondary screens [ Fig. 3C]. Analysis of the MOAs of these compounds revealed that nearly half (n = 16, 44.4%) are EGFR inhibitors, followed by nine mTKIs and three MEKi [Fig. 3D]. The observation that EGFR inhibitors are significantly enriched and more effective in HNSCC compared to other non-HNSCC is consistent with the observation that EGFR gene dependency is also selectively enriched among HNSCC [ Fig. S4A]. A previous HNSCCspecific analysis of the GDSC also found EGFR inhibitors (afatinib and gefitinib) among the four drugs (50%) with significantly higher sensitivity in HNSCC 18 .  Table S4]. Of the 22 (61%) non-FDA approved drugs, 19 were considered "novel" as they have not been tested in HNSCC before [ Fig. 3E] [ Table S4]. There are three compounds (osimertinib, neratinib, and bosutinib) that represent promising novel candidates for drug repurposing, as these are FDA-approved, but not tested in clinical trials on HNSCC patients. Osimertinib is approved for EGFR-mutated NSCLC (e.g. exon 19 deletions, p.L858R and p.T790M). However, these EGFR mutations are rare in HNSCC 19 . Interestingly, when we compare the drug response data of unselected HNSCC versus that of NSCLC with or without the actionable mutations, HNSCC sensitivity towards osimertinib was comparable with that of EGFR-mutated NSCLC (P = 0.4180), and also significantly more sensitive than EGFR-wildtype NSCLC (P < 0.0001) [Fig. 4A]. Some of the HNSCC cell lines are even more responsive than the NSCLC with EGFR mutations. This suggests that despite the lack of actionable mutation in EGFR, HNSCC are still very dependent on EGFR signalling. In fact, osimertinib treatment resulted in a significant antitumor effect in an in vivo study using two EGFR-wildtype HNSCC models (FADU and CAL27) 45 . Therefore, the repurposing of osimertinib to inhibit this essential EGFR signalling in HNSCC is promising. Unlike in NSCLC, where activating mutation, high EGFR expression (Pearson's R = − 0.3675, P = 0.0120) and dependency (Pearson's R = 0.4757, P = 0.105) were associated with osimertinib sensitivity, EGFR expression (Pearson's R = 0.2957, P = − 0.1341) and dependency (Pearsons' R = 0.2723, P = 0.1980) were not associated with osimertinib sensitivity in HNSCC [ Fig. S4B-C]. However, further mechanistic investigation is needed to delineate the mode of action and to identify any relevant biomarker in the HNSCC context.

Identification of HNSCC tissue-specific biomarker of response for osimertinib.
A recent proteogenomics study of HNSCC has proposed two models of EGFR activation: ligand-dependent and ligandindependent pathways 20 . In the context of HNSCC, the ligand-dependent mode is often observed, where the availability of EGFR ligands is the rate-limiting factor where the abundance of EGFR ligands (AREG, EREG, HBEGF, EPGN and TGFA), are associated with EGFR pathway activation and response towards cetuximab in the patient-derived xenograft (PDX) models 20 . Hence, we examined if the sensitivity towards osimertinib can also be better predicted by the EGFR ligands abundance. We correlated the expression levels of the five EGFR ligands    www.nature.com/scientificreports/ individually, as well as their average expression, with osimertinib sensitivity [Fig. 4B]. We found that the average expression of all EGFR ligands, instead of EGFR expression, is significantly correlated with the osimertinib sensitivity [ Fig. 4C], corroborating the recent proteogeonomic findings. In contrast, in NSCLC (ligand-independent mode), the EGFR expression is significantly correlated with osimertinib sensitivity [Fig. S4D]. This further highlight the tissue-specific difference in the mechanism underlying drug sensitivity and biomarker selection. Development of resistance towards tyrosine kinase inhibitor is common, but an earlier understanding of intrinsic resistance could shed light on how to manage unresponsiveness or acquired resistance. We performed GSEA on the four most sensitive and four most resistant lines to identify activated pathways that might govern the response towards osimertinib [ Table S5]. The TGF-beta signaling pathway is significantly enriched (normalized enrichment score = 1.4347; nominal p-value = 0.0051) among the resistant lines [ Fig. 4D]. Significantly upregulated genes included the ligands of TGF-beta receptors (TGFB1, TGFB2, and TGFB3), as well as the effector proteins (SMAD4, SMAD6, and SMAD9) and some of its target genes (ID1, ID2 and ID3) [Fig. 4E].
Uncovering biomarker of response/resistance for MEKi. We found three MEKi (trametinib, cobimetinib and AZD8330) among the compounds that were enriched for killing activity in HNSCC when compared to non-HNSCC [ Fig. 2D]. Trametinib and cobimetinib (in combination with vemurafenib) have been approved to treat unresectable or metastatic melanoma with BRAF V600E or V600K mutations 21 . Although BRAF mutations are rare in HNSCC, mutations in the MAPK pathway are found in 18.6% of HNSCC 22 , suggesting that targeting MAPK pathway could be a promising strategy in HNSCC. Trametinib has only been tested in a window-of-opportunity trial [NCT01553851], where it resulted in a median of 46% reduction in tumour size among 17 HNSCC patients 23 . However, the biomarkers that were being investigated (p-ERK1/2 and CD44) did not correlate well with response 23 .
We next sought to explore potential biomarkers of MEKi in HNSCC. Of the 20 MEKi in PRISM, 14 (70%) were active compounds in HNSCC. Unsupervised hierarchical clustering showed a distinct cluster consisting of eight active MEKi (red branches) [ Fig. 5A] [ Table S6]. Six cell lines were distinctively non-responsive towards all the MEKi, while the other 22 cell lines show variable response. To identify potential gene expression signature or single gene marker that may predict response, we performed DEG analysis between the six most sensitive and six most resistant cell lines. We identified 136 significantly up-regulated genes among the MEKi-sensitive HNSCC, and 428 significantly down-regulated genes [ Fig. 5B] [ Table S7]. Among the 136 genes, significant enrichment of genes in the Reactome pathway of "Cytokine Signalling in Immune System" (HSA-1280215) was observed. As shown in the STRING functional protein association analysis [ Fig. 5C], 24 of the genes in this pathway are significantly up-regulated and interconnected. Among these, we further shortlisted six cytokines (IL1A, SAA1, LCN2, CSF2, IL1B and CXCL1), where their gene expression significantly correlated with MEKi sensitivity across the 28 HNSCC cell lines [ Fig. 5D]. To complement the DEG analysis, we also performed GSEA and showed that immune-related pathways such as the interferon gamma/alpha responses and inflammatory response hallmark pathways were significantly enriched among the MEKi-sensitive group [ Fig. 5E]. On the other hand, proliferation or cell-cycle related hallmark pathways such as the G2M checkpoint, MYC targets (v2) and E2F targets are significantly enriched among the MEKi-resistant group [ Fig. 5E]. To confirm our findings, we performed GSEA on MEKi-resistant HNSCC from the GDSC (v2) and Lepikhova et al. dataset [Table S8]. G2M checkpoint, MYC targets (v2), E2F targets and spermatogenesis were commonly enriched hallmark pathways [ Fig. 5F]. This further substantiates that aberrantly activated cell cycle and pro-proliferative pathways are associated with resistance towards MEK inhibition. Up-regulation of genes such as CCND1, a critical cell-cycle regulator is observed among the MEKi-resistant HNSCC. CCND1 overexpression is known to promote tumorigenesis in various cancers, and CCND1 amplification is a common event in HNSCC resulting in the dysregulation of cell cycle pathways. The CCND1-CDK4/6 complex phosphorylates Rb and releases E2F into the nucleus, leading to transcriptional activation of downstream E2F targets 24 . Whilst clinical trials utilizing CDK4/6 inhibitors are underway, this data suggests that a combination of CDK4/6 inhibitors such as palbociclib and abemaciclib with MEKi could be an effective strategy to overcome intrinsic resistance towards MEKi in HNSCC. www.nature.com/scientificreports/

Discussion
The increased availability of high-throughput pharmacogenomics data have provide an opportunity to perform drug repurposing in specific cancer type. Here, we have performed a systematic in silico analysis on the PRISM drug screening data to identify compounds for HNSCC. We demonstrated that some of the non-oncology compounds could be exploited for this devasting disease. More importantly, we identified and validated the findings using independent data sets to gain insights the MOA of these compounds and their potential biomarkers for HNSCC.
Our HNSCC-specific analysis of the PRISM dataset has uncovered a rich resource of targeted compounds that holds the potential to be repurposed, and also revealed a large number of active non-oncology compounds that could have activity in HNSCC. Despite new anti-cancer drugs being tested on HNSCC patients, most fail to show clear clinical benefit. This could be due in part to the absence of biomarkers for patients stratification, as the use of biomarkers was associated with a four-fold increase in the success rates for oncology clinical trials 25 . One of the unique propositions of this dataset is that all 28 HNSCC lines have deep molecular characterization, where genomics, transcriptomics, epigenomics, proteomics and essentialomics (essential genes profile from CRISPR/ Cas9 screens) data are available 7 . Integration of multi-omics data could enable the identification of biomarker of response and improve our understanding of drug MOAs and potential resistance mechanism. Among these features, the transcriptomics (mRNA expression) was found to be most predictive of drug sensitivities 7,26 . Our findings shed light on the development of gene-expression based biomarker for HNSCC, especially where mutations in actionable genes are relatively rare in this cancer type.
A limitation while interpreting this dataset is that this PRISM primary screen is done at a single dose of 2.5 µM. Some highly potent compounds might have saturated anti-cancer effect at this dose and cannot be differentiated from those that are less potent. Further, the use of logFC relative to the DMSO control as a surrogate marker of drug sensitivity, instead of IC50 or EC50 (half-maximal effective concentration), limits our ability to determine the potency of each of the compounds. While this could be overcome by examining the secondary PRISM screen, there was a substantial reduction in the number of compounds tested (4,518 in primary screen vs 1,207 in secondary screen) and only 22 HNSCC lines were included. Furthermore, there were missing data points in the secondary screen where, of the expected 1,207 × 22 data matrix, only 60% of the data points were available [ Fig. S5]. Hence, the secondary screen data was used only for validating the data from the primary screen. Furthermore, due to the pooling of a mixture of molecular-barcoded cancer cell lines, the drug response might be confounded by paracrine-mediated modulation 6,7 .
Whilst an inherent limitation of high-throughput drug screening precludes the role of the stromal and immune compartments, known clinical responses were well recapitulated in this PRISM dataset, demonstrating the utility and relevance of cancer cell line-derived drug response data. Large-scale analysis found that oncogenic alterations in tumour tissue can be faithfully recapitulated by cancer cell lines 5 . Besides, a recent pan-cancer study called the Celligner, revealed that HNSCC cell lines (which included the 28 lines tested in PRISM) rank second after Ewing sarcoma, where the transcriptomics profile of cell lines and patient tumours are highly correlated 27 . In the HNSCC dataset, we identified possible FDA-approved drug candidates including osimertinib and trametinib that could be evaluated for the treatment of HNSCC. Further, we also offer insights into the potential candidates of non-oncology compounds to be repurposed for HNSCC, which can be explored further with our RShiny app or the interactive heatmap.
Using tsne projection of the drug sensitivity profile [ Fig. 2D], some non-oncology drugs showed similar activities with subsets of targeted compounds with defined MOAs. Piroxicam and FK-3111 which are nonsteroidal anti-inflammatory drugs (NSAIDs) approved for treating pain and inflammation, are two examples of cyclooxygenase inhibitors that clustered closely with EGFR inhibitors. Epidemiological studies have reported associations of NSAIDs use with reduced cancers risk, thus prompting the investigation of their anti-cancer properties 28 . Piroxicam is used widely to treat OSCC in dogs 29 , its anti-cancer properties have also been demonstrated in several human cancers [30][31][32] . In pre-malignant and malignant human oral cancer lines, piroxicam selectively inhibited malignant cell growth via cell cycle arrest in the S phase 33 .
Another interesting non-oncology drug is salinomycin, that clustered with PI3K/AKT/MTOR inhibitors. Salinomycin is widely used as an anticoccidial drug in poultry but was later found to be a potent and selective inhibitor of epithelial cancer stem cell (CSC) growth from a high-throughput screen 34 . Salinomycin preferentially inhibits breast CSCs at a potency of more than 100-fold that of paclitaxel and resulted in potent suppression of tumour growth in vivo 34 . Subsequently, numerous studies showed that salinomycin is effective as an anti-cancer drug in various cancers [35][36][37][38] . In HNSCC, salinomycin treatment resulted in a 71.5% increase in apoptosis and significantly reduced the sphere-forming ability of HNSCC CSCs 39 . Interestingly, salinomycin was shown to induce phosphorylation of AKT, this is consistent with the effect of specific AKT inhibitors when tested in HNSCC 40,41 . This further suggests that salinomycin may share similar MOA with the AKT inhibitors, as indicated by their clustering in the tsne. Notably, salinomycin (specifically HSB-1216, delivered via QUATRAMER technology) has been granted orphan drug designation by the FDA in 2020, and a Phase I trial is planned for small cell lung cancer [42][43][44] . Further investigation would be of interest to uncover the potential utility of salinomycin for HNSCC treatment, and our discovery here may shed light into uncovering its MOA and enabling selection of patients with biomarker of response. To facilitate more discovery research, our interactive heatmap and the tsne plot on RShiny app (https:// annie chai. shiny apps. io/ tsne/) would serve as a useful starting point.
Importantly, from the thousand of compounds screened, we have shortlisted 36 promising compounds [ Fig. 3C-E] that should be prioritized for further investigation, due to their enriched killing activities among HNSCC lines. Although EGFR is an enriched dependency in HNSCC and EGFR inhibitors are most active in killing HNSCC relative to other cancers, Cetuximab is the only approved targeted therapy for use in HNSCC. In this study, we focused on osimertinib, a third-generation EGFR inhibitor approved for NSCLC patients with www.nature.com/scientificreports/ exon 19 deletions, L858R and T790M point mutations in EGFR. Despite not having these mutations, osimertinib was found to be equally effective in HNSCC, with no significant difference in logFC when compared to the EGFR mutation-positive NSCLC. In HNSCC, EGFR expression/amplification is not predictive of response to EGFR inhibitors 20,45,46 . Consistent with more recent studies on cetuximab and panitumumab, we have demonstrated that the expression of the EGFR ligands but not EGFR is more predictive of osimertinib response. This supports the notion that in HNSCC, the rate-limiting factor governing EGFR pathway activation is the expression of EGFR ligands (ligand-dependent pathway) 20 ; while in another context such as NSCLC, the EGFR ligand-independent pathway might be prevailing 47 . Using GSEA, we found that activation of TGF-beta signalling might play a role in mediating resistance towards osimertinib. The use of TGF-beta blocking antibody diminishes the emergence of cetuximab-resistant HNSCC 48 and TGF-beta activated CAFs were shown to suppress cetuximab activity in HNSCC in vitro and in vivo PDX models 49 demonstrating that integrative analysis of drug response and genomics features in cancer models could identify mechanisms underlying response to cancer therapeutics. Although a preclinical study has shown the efficacy of osimertinib as a single agent in inhibiting in vivo tumour growth 50 , no clinical trial yet has tested osimertinib in HNSCC. Our data suggests that osimertinib could be efficacious in HNSCC expressing high levels of EGFR ligands, with low expression of members of the TGF-beta pathway, and clinical trial designs that examine and validate these biomarkers is warranted. Another class of targeted cancer compounds with promising efficacy in HNSCC is MEKi. Clear selectivity of MEKi are seen where subsets of HNSCC cell lines demonstrated distinct responses to MEKi suggesting that predictive markers of response are crucial for patient selection. Unlike in melanoma, where BRAF mutation status could help to identify responders to vemurafenib, other mechanisms of response must be at play in HNSCC as BRAF mutations are uncommon in HNSCC 51 . Shared mutations among the six most resistant lines were not found, highlighting that the lack of response could be heterogeneous and likely not driven by mutations. Emerging studies have revealed that transcriptomic profiles could be more powerful in predicting drug responses compared to mutational profiles in some context 7,26 . From our DEG analysis, we revealed that up-regulation of several cytokines and immune-related pathways are associated with sensitivity towards MEKi. This is consistent with the increasingly recognized immunomodulatory roles of the MAPK pathway in cancers 52 . Among these are interleukin genes such as the IL1A and IL1B, which are pro-inflammatory cytokines that could induce activation of the MAPK pathway 53 Besides, SAA1, CXCL1 and CSF2 have also been implicated to regulate the MAPK pathway 54 . HNSCC has been shown to be one of the most immune-enriched solid tumors with the majority belonging to the IFNγ-dominant immune subtype 55 . The upregulation of pro-inflammatory cytokines that confers sensitivity towards MEK inhibition could therefore explain why a majority of HNSCC cell lines are responding to MEKi despite the lack of activating mutations in the Ras/Raf/MAPK pathways. This concept is currently being tested in an on-going clinical trial [NCT03264066] of combining cobimetinib with atezolimumab (anti-PDL1 mAb) in HNSCC and other cancers. Some HNSCC however, were intrinsically resistant towards MEK inhibition. In this regard, we found cell cycle and proliferation-related pathways to be consistently associated with resistance towards MEKi. This supports the combination of CDK4/6 inhibitors with MEKi for better synergy, which has been investigated in other cancers, such as melanoma, colorectal and NSCLC 56 . However, these were investigated in the context of the presence of mutation(s) in the Ras/Raf/MAPK pathway which is lacking in HNSCC. Instead, in HNSCC, we showed that the immune-cytokines enriched HNSCC could promote responsiveness to MEKi. The potential utility of those shortlisted biomarkers of response, as well as the efficacy of combining MEKi with CDK4/6 warrant further investigation in HNSCC.

Conclusions
In summary, our HNSCC-specific analysis of the PRISM repurposing dataset has helped to reveal the targetable MOAs within HNSCC and shortlisted 36 drug candidates to prioritize for repurposing in HNSCC. Our processed HNSCC drug response dataset that is made available in the form of interactive heatmaps or RShiny app webpage is an important resource that could accelerate drug development in the context of HNSCC. We have also identified several HNSCC-specific biomarkers for EGFR inhibitor and MEKi that warrant further validation. Furthermore, using transcriptomic-based pathway-level enrichment analysis, we demonstrated how rational drug combination may be proposed, such as in the case of MEKi and CDK4/6 inhibitor. These exemplify the potential of our study in enabling new discoveries in identifying promising new drugs and biomarkers for further investigation, to expand the treatment options for HNSCC. In addition, our in silico analysis pipeline can also be applied to other cancers covered in this PRISM dataset, for their corresponding drug repurposing exploration.