Integrated Genomic Profiling and Drug Screening of Patient-Derived Cultures Identifies Individualized Copy Number-Dependent Susceptibilities Involving PI3K Pathway and 17q Genes in Neuroblastoma

Neuroblastoma is the commonest extracranial pediatric malignancy. With few recurrent single nucleotide variations (SNVs), mutation-based precision oncology approaches have limited utility, but its frequent and heterogenous copy number variations (CNVs) could represent genomic dependencies that may be exploited for personalized therapy. Patient-derived cell culture (PDC) models can facilitate rapid testing of multiple agents to determine such individualized drug-responses. Thus, to study the relationship between individual genomic aberrations and therapeutic susceptibilities, we integrated comprehensive genomic profiling of neuroblastoma tumors with drug screening of corresponding PDCs against 418 targeted inhibitors. We quantified the strength of association between copy number and cytotoxicity, and validated significantly correlated gene-drug pairs in public data and using machine learning models. Somatic mutations were infrequent (3.1 per case), but copy number losses in 1p (31%) and 11q (38%), and gains in 17q (69%) were prevalent. Critically, in-vitro cytotoxicity significantly correlated only with CNVs, but not SNVs. Among 1278 significantly correlated gene-drug pairs, copy number of GNA13 and DNA damage response genes CBL, DNMT3A, and PPM1D were most significantly correlated with cytotoxicity; the drugs most commonly associated with these genes were PI3K/mTOR inhibitor PIK-75, and CDK inhibitors P276-00, SNS-032, AT7519, flavopiridol and dinaciclib. Predictive Markov random field models constructed from CNVs alone recapitulated the true z-score-weighted associations, with the strongest gene-drug functional interactions in subnetworks involving PI3K and JAK-STAT pathways. Together, our data defined individualized dose-dependent relationships between copy number gains of PI3K and STAT family genes particularly on 17q and susceptibility to PI3K and cell cycle agents in neuroblastoma. Integration of genomic profiling and drug screening of patient-derived models of neuroblastoma can quantitatively define copy number-dependent sensitivities to targeted inhibitors, which can guide personalized therapy for such mutationally quiet cancers.


INTRODUCTION
Neuroblastoma is the most common pediatric extracranial malignant tumor and is responsible for a disproportionate 15% of all childhood cancer deaths. Despite current intensive multimodal therapy for patients with high-risk disease, 5-year survival remains at 30-50% even after decades of international multicenter trials. Due to its histologic and biologic heterogeneity, the classic stage-and risk-based clinical trial approach cannot adequately match treatments to the diverse individual susceptibilities of each patient's tumor, and is further limited by the small numbers of patients in each subgroup (1). Thus, therapeutic advancements for neuroblastoma may be better realized with a personalized approach.
Precision oncology approaches based on identification of targetable single nucleotide variants (SNVs) have limited usefulness in embryonal tumors because of their low mutational burden (2). Instead, copy number variations (CNVs) are more prevalent and are stronger prognostic factors in neuroblastoma (2)(3)(4). In particular, segmental chromosomal aberrations (SCAs) are associated with advanced disease stage and poorer prognosis. This corroborates with recent evidence that the pathogenicity of CNVs correlate with dosage sensitivity of involved genes, and are enriched for embryonal neurodevelopmental functions (5). This suggests that gene copy number could be used to predict and select targeted therapies, especially in embryonal tumors of childhood (6). While current sequencing-based panels may be an efficient and cost-effective manner to perform clinical genomic profiling, most are designed for adult cancers or hematological malignancies. Recently, the Oncomine Childhood Cancer Research Assay (OCCRA) was developed as a diagnostic-grade genomic profiling tool curated specifically for pediatric cancers (7,8). This is a promising new resource to detect significant and potentially actionable SNVs and CNVs in neuroblastoma.
Individualized preclinical tumor models can be a complementary means to rapidly test drug therapies ex vivo, and have been able to successfully uncover potential therapeutic leads (9,10). While this has been employed extensively in epithelial carcinomas, embryonal tumors like neuroblastoma have not been consistently engrafted as personalized in vitro models. Commercial cell lines are often significantly changed with multiple serial passages or have been contaminated over time, and do not capture individual patient genetic or phenotypic differences in treatment response (11)(12)(13)(14). We previously developed multi-lineage patient-derived cell cultures (PDCs) of neuroblastoma from pre-and post-treatment tumors and demonstrated their recapitulation of original tumors' chromosomal alterations, immunohistochemical and gene expression profiles, and ability to predict individualized responses to standard-of-care chemotherapy (15). However, the utility of patient-derived models for prediction of targeted therapies for neuroblastoma have not been well studied.
We hypothesize that in neuroblastoma, gene copy number and in vitro cytotoxicity to corresponding targeted inhibitors display distinct relationships in a dose-dependent manner. To study key genotype-phenotype correlations in neuroblastoma, we characterized SNVs and CNVs of neuroblastoma tumors using clinical genomic profiling and interrogated corresponding PDCs with a medium-throughput inhibitor screen, then correlated resulting genomic and phenotypic readouts. We then defined the quantitative associations between these genomic aberrations and the corresponding responses to targeted agents as a framework to potentially guide personalized targeted therapeutic options, particularly for patients with few or no targetable mutations.

Patients and Tumor Samples
Patients with neuroblastoma were prospectively recruited at KK Women's and Children's Hospital with Institutional Review Board approval (CIRB 2014/2079). Written consent from parents, and assent from children were obtained. Criteria for study inclusion were: patients under undergoing surgical biopsies or resections, with available excess tumor tissue; male and female patients with neuroblastoma aged 1-8 were included. No subject attrition was encountered, no randomization or power calculation was required, and investigators were blinded to patient identities. Excess tumor tissue from routine surgical procedures were transported on ice on the same day to the VIVA-KKH Paediatric Solid Tumour Laboratory for generation of PDCs. Corresponding tumor aliquots were snap frozen for molecular analysis.

OCCRA Comprehensive Genomic Profiling
Access to OCCRA (RRID : SCR_007834) (Thermo Fisher Scientific, Waltham, MA) was granted via an early access program. The DNA assay was utilized, which calls SNVs from hotspots of 86 genes and full exons of 44

Cell Culture
PDCs were generated as previously described by explant culture under growth conditions of 37°C with 5% CO 2 (15), and passaged at 1:2 split ratio upon reaching 80-100% confluence (15). For each PDC line, optimal seeding density to achieve 72hour log-phase growth was determined. Cells of each PDC line were seeded on clear bottom white 96-well plates in triplicate (Corning, Cat 3903) in a 1:2 split ratio and assessed for appropriate responses to positive and negative control agents 0.5% DMSO and 1mM Staurosporine over 72-hours using IncuCyte ® S3 Live-Cell Analysis System (RRID : SCR_019874) (Essen BioScience, Sartorius, Japan) and IncuCyte Base Software (v.2018B), using default analysis settings.

Short Tandem Repeat Fingerprinting
STR genotyping of all tumor-PDC pairs was performed using PowerPlex ® 21(Promega, Cat DC8902) using 1ng DNA. PCR products were resolved in an Applied Biosystems ® SeqStudio genetic analyzer and compared.

Medium-Throughput Screening Assay
Selleckchem Kinase Inhibitor Library (SelleckChem, Cat. L1200), consisting of 418 tyrosine kinase inhibitors, was reformatted into 96-well format and diluted with DMSO to achieve final concentrations of 1uM. Daughter plates were stored at -20°C and underwent no more than 15 freeze-thaw cycles.
Cells were seeded in a 1:2 dilution factor on 96-well white bottom plates (Corning, Cat 3917) using MultiFlo ™ FX (RRID : SCR_019746) (Biotek, Winooski VT). Plates were incubated for 24 hours at 37°C with 5% CO 2 , and 0.5mL of each drug was added to the corresponding wells using Bravo BenchCel Workstation (RRID : SCR_019468) (Agilent, Sata Clara CA) equipped with automated liquid and microplate handling. Treated plates were centrifuged for 1 minute at 1000rpm before and incubated for the next 72 hours. Cytotoxicity was determined using CellTiter-Glo ® Luminescent Cell Viability Assay (Promega, Cat G7572). Media was aspirated before adding 50uL of CellTiter-Glo reagent prepared as per manufacturer's protocol. Plates were shielded from light and agitated at 300rpm for 15 minutes; bioluminescence was measured using Infinite M1000 series microplate reader (Tecan, Mannerdorf, Switzerland) using a luminescence integration time of 250ms. Relative light units (RLUs) were normalized against DMSO negative controls of each corresponding plate to obtain normalized cytotoxicity values. Screening runs were validated only if high cytotoxicity to positive control staurosporine (1mM) and minimal DMSO effect were observed

Data Analysis
Clustering of copy number and z-score data was performed on ComplexHeatmap (v.2.0.0, RRID : SCR_017270) using R 3.6.1 with complete Euclidean clustering applied.
Pearson's correlation coefficients were calculated for each of 56,848 gene-drug pairs, comparing cytotoxicity and SNV, and cytotoxicity and copy number, using Fisher's z-transformation approach. Positive or negative drug-gene correlations were obtained for all 56,848 gene-drug combinations for both mutational variants and gene copy number. False discovery rates (FDR) for the Benjamini-Hochberg procedure were set at 0.1 as a threshold of statistical significance of gene-drug pairs (16). Data analysis was performed using R 4.0.3.

Probabilistic In Silico Model
Probabilistic Graphical Models (PGMs) were constructed with ReactomeFIViz (v. 7.2.3) using Reactome Functional Interaction (FI) network version 2018, which utilizes an adaptation of the PARADIGM approach (17)(18)(19) (RRID : SCR_009634). Briefly, this infers case-specific genetic variations, incorporating information from curated Reactome pathways, converting these reactions into factors then modeling the variations in each gene as constraint graphs which represent them as probability distribution functions. CNVs were input as continuous observation variables without discretizing and pairwise Markov random field models using empirical distributions were constructed over 100 permutations. Genes with mean differences in protein impact scores of <-0.01 and >0.01 between real and random samples were selected; genes with no CNV were removed. The PGM network was overlaid with associated Cancer Targetome drugs (20). A corresponding network was constructed from gene-drug pairs weighted according to z-scores, using Cytoscape (ver 3.8.0, RRID : SCR_003032). Both networks were clustered according to Reactome FI scores. Sub-networks (Reactome FI Modules) were annotated for significantly enriched pathways from CellMap (RRID : SCR_010642), Reactome, KEGG (RRID : SCR_012773), NCI Pathway Interaction Database, Protein ANalysis THrough Evolutionary Relationships (PANTHER), and BioCarta databases (RRID : SCR_006917).

Molecular Profiling of Neuroblastoma Identifies Few Actionable Pathogenic Mutations but Multiple Copy Number-Altered Genes Involved in Transcriptional Regulation
Fresh tumor specimens were obtained from surgical biopsies and resections of 13 patients with neuroblastoma following informed consent. Median age was 2.65 (range: 1.45-7.80). The study population appropriately demonstrated the expected spectrum of pathological and cytogenetic risk features in neuroblastoma: 9/13 (70%) had stage 4 disease and 3/13 (23%) had MYCN amplification (Supplementary Table S1). Tumor specimens were subjected to genomic profiling while corresponding PDCs were subjected to phenotypic drug-response testing ( Figure 1A).
Genomic profiling was performed using the OCCRA DNA panel. Among 130 genes interrogated for SNVs, 71 separate SNVs in 12 patients fulfilled criteria as sequence variants of IARC Class 3 and above according to American College of Medical Genetics standards (21); 27 of 71 SNVs were recurrent single nucleotide polymorphisms (SNPs) frequently encountered in our local population and excluded. Of the remaining 44 SNVs, 34 were Indels, amounting to a mean incidence of 3.1 mutations per case ( Figure 1B). The most frequent mutations were encountered in TSC2, TP53, ARID1A and ARID1B (Supplementary Table S2). Among these, 11 were predicted by Sift scores to be deleterious ( Figure 1B), and by PolyPhen scores, 5 as probably damaging and 4 as possibly damaging. TSC2 and TP53 mutations were observed in recognized hotspots (21) ( Figure 1C). Two patients had germline splice variants of SMARCB1, though neither exhibited clinical phenotype of the associated Coffin-Siris syndrome (22) ( Figure 1C).

High Throughput Screening of Neuroblastoma PDCs Reveals Drug Hits
Targeting key Molecular Pathways Implicated in Neuroblastoma Concurrently, PDCs were generated from the same 13 tumors and their recapitulation of the corresponding original tumors were validated as previously described (15)  . The compound library was restricted to kinase inhibitors, to constrain target selectivity and limit offtarget effects. Drug classes associated with the highest mean normalized cytotoxicity targeted cytoskeletal signaling, cell cycle, angiogenesis and PI3K/Akt/mTOR pathway ( Figure 3A).
To determine the clinical significance of drug screen results, we correlated cytotoxicity and clinical variables using nonhierarchical unsupervised clustering, but no significant associations between PDC drug-response profiles and clinical variables were observed (Supplementary Figure S2) (28).

Correlation of Cytotoxicity and Copy Number Identify Mostly Novel Gene-Drug Associations
Having separately characterized the molecular aberrations of the patient tumors and the drug response profiles of their PDCs, we then sought to correlate both genomic and phenotypic readouts. Copy number and mutation calls of 131 genes and normalized cytotoxicity to 418 compounds were compared to determine significantly correlated gene-drug pairs. The z-score of Pearson correlation coefficients of 54,758 possible gene-drug combinations were normally distributed for both copy number and SNVs ( Figure 4A). Selecting for significantly-correlated gene-drug pairs at 0.1 false discovery rate (FDR) cutoff, copy number and cytotoxicity were significantly correlated in 1278 gene-drug  Figure S3). Among these were 49 hit compounds with >33% mean cytotoxicity, which involved 159 significantly correlated gene-drug pairs ( Table 1).

combinations (Supplementary
Due to the small number of mutations, no significant genedrug combinations could be identified from a similar correlation of SNVs and cytotoxicity ( Figure 4A). Thus, we subsequently focused on evaluating the gene-drug pairs with significantly correlations of copy number and cytotoxicity.
To further study these observations, we compared the expected gene-drug associations (according to drug library specifications) with the significantly correlated gene-drug pairs (according to z-scores of copy number-cytotoxicity correlations) identified ( Figure 4B). A total of 416 gene-drug pairs were expected (Supplementary Data 1). Among these, copy number and cytotoxicity of 4 gene-drug pairs showed significant correlations at 0.1 FDR cutoff (true positive pairs): JAK3 with AT9283, JAK1 with CYT387, FGFR1 with tyrphostin (AG1296), and FGFR3 with brivanib (BMS540215) ( Figure 4C). Only the former 2 drugs demonstrated >33% cytotoxicity and qualified as hits ( Figure 4C).
In 1274 gene-drug pairs, correlations were not expected but were observed. To identify pairs with the strongest treatment effect and correlation, we ranked them according to the product of their cytotoxicity and z-score (see Supplementary Data 1). Among the top-and bottom-ranked pairs, were PI3K and cell cycle agents that correlated with copy number losses of transcriptional regulators CBL and DNMT3A (log 2 ratios -0.294, and 0.039, respectively), and copy number gains of GNA13 and PPM1D ( Figure 4D). The drugs most commonly associated with these pairs were PI3K/mTOR inhibitor PIK-75, and CDK inhibitors P276-00, SNS-032, AT7519, flavopiridol and dinaciclib (Supplementary Figure S4). Intertumor heterogeneity of copy-number dependent drug sensitivities were observed. For example, MYCN non-amplified tumors showed increased drug sensitivity to these agents with increased PPM1D and GNA13 copy number, while MYCN amplified tumors did not show similar trends (Supplementary Figure S4). As this was a substantial number of novel gene-drug associations that were not expected, we next verified them against available data sources and via in-silico modelling.

Significantly Correlated Gene-Drug Pairs Are Verified in Public Datasets
To verify if the significantly correlated gene-drug pairs were consistent with published literature, we compared them against known gene-drug associations in 6 public datasets. In each dataset, median 224 (range:46-965) pairs were identified as having an established association of drug and gene target (Supplementary Figure S5 and Supplementary Data 1). In all, 1305 of 54,758 (2.4%) pairs were identified in at least 1 dataset. Verified and unverified gene-drug pairs did not significantly differ in the distribution of cytotoxicity of the drug involved, copy number of the gene involved, and z-score of the pair across the 6 datasets (Supplementary Figure S5).
Among the 1278 significantly correlated gene-drug pairs, 19 pairs (1.5%) were verified in at least 1 data set (Supplementary Table 4) ( Figure 5A), though only 1 pair (5.3%) qualified as a hit with >33% mean cytotoxicity [JAK3 with AT9283 (mean cytotoxicity 64.3 ± 16.8%)]. Among the remaining unverified 1259 gene-drug pairs, a higher percentage of pairs [12.9% (n=162)] were predicted to be hits with >33% mean cytotoxicity. Among these pairs, cytotoxicity was positively   Since published data on this majority of the significantly correlated gene-drug relationships was lacking, we performed an in-silico prediction of how the CNV profile of our neuroblastoma tumors would be expected to influence drug response, and compared that with the corresponding z-scorebased relationships. We constructed a Probabilistic Graphical Model (PGM) network by weighting our CNV data with Reactome pathway-based interactions to infer variations in inter-gene relationships (17,19), then overlaying networks with associated Cancer Targetome drugs (20) (Figure 5B). A corresponding network was constructed from the 1278 significantly correlated gene-drug pairs and weighted according to z-scores ( Figure 5C). Both networks were clustered according to Reactome Functional Interaction (FI) scores, and subnetworks of commonly represented drugs were examined. Since the curated source included only FDA-approved cancer drugs, as expected only verified gene-drug pairs were found in the PGM network (Supplementary Data 1). Among them, genedrug pairs involving sunitinib and vandetanib were modeled by the PGM networkboth receptor tyrosine kinase inhibitors with known effects against neuroblastoma (29)(30)(31). Examples in Figure 5D show how the PGM, derived only from CNV data, predicted associations between sunitinib malate and ALK and MAP2K2, and between vandetanib and ACVR1 and FGFR1which were similarly demonstrated in the z-score-weighted network. The similarities observed between the CNV-based PGM predictions and the z-score-weighted network supported the integrity of the gene-drug associations identified by our experimental system.
Since FI-based clustering might reveal associations between drugs and signaling pathways with potential functional significance, we isolated FI-clustered sub-networks of genedrug pairs and annotated them for enriched Reactome pathways and Gene Ontology features (Supplementary Data 2). Interestingly, 2 sub-networks coincided with the unverified significantly associated gene-drug pairs that were earlier identified (Supplementary Figures S6A, B). In the former, copy number of ALK, DNMT3A and ASXL2 was inversely associated with efficacy of multiple kinase inhibitors like saracatanib (Module 1). In the latter, copy number of 17q genes CBL, PPM1D and GNA13 were positively associated with cytoxicity to dinacilib, flavopiridol and ponatinib (Module 8). Other subnetworks shared common involvement of the PI3K-mTOR (Modules 3, 6, 11), JAK-STAT (Modules 0, 3, 5, 6) pathways. Notably the z-score network also illustrated the lack of inhibitors associating with MYCN despite its substantial degree of copy number amplification.

DISCUSSION
CNVs are more prevalent in mutationally quiet embryonal cancers like neuroblastoma and may reflect fundamental oncogenic dependencies that can be potentially exploited for cancer therapeuticsparticularly since this phenomenon has already been utilized to target oncogene amplifications in adult cancers. In this study we hypothesized that CNVs may modulate gene dosage and correspondingly, drug response to respective targeted inhibitors. Thus, we integrated diagnostic-grade genomic profiling with phenotypic readouts of personalized preclinical tumor models to identify copy number-dependent  therapeutic susceptibilities in patients with neuroblastoma. We used correlational coefficients to quantify the strength of association between tumor copy number and cytotoxicity of PDCs in corresponding gene-drug pairs, and validated this against public data and in silico. The findings defined a prioritized set of targetable copy number-altered genes and their ideal corresponding inhibitors for copy-number dependent drug targeting. To provide a statistically robust context to apply these in clinical practice, we generated prediction models of the expected therapeutic response for given degrees of copy number change for each gene-drug pair. These potentially form the basic elements of a decision-making framework that could be used to inform personalized therapeutic decisions by predicting potential drug hits based on CNVs. While the identified SNVs were known genomic aberrations in neuroblastoma (32)(33)(34), they were mostly unactionable. This underscored the limitations of targeted mutational profiling for identifying personalized treatment options for neuroblastoma and other mutationally quiet embryonal cancers (2,35). Among the commonest variants observed in our study cohort, TP53 mutations are still largely considered untargetable (36,37). Similarly, options for TSC2-deficient tumors are limited and largely aimed at downstream mTOR pathway targets (38). Frequency of these SNVs differed from Western series where ALK, TERT, PTPN11 and NRAS mutations are more frequent, likely reflecting known inter-ethnic differences in mutational incidence. For example, de novo somatic TP53 mutations are uncommon in Western series of neuroblastoma (33), but polymorphisms such as R72P R337H are more frequent in non-Caucasian populations and are associated with varying incidence of neuroblastoma (39,40). Similarly, recurrent mutations in TSC, ARID1A and ARID1B have been reported at a lower incidence (1.7%, 11% and 6.9%, respectively) in Western populations (41)(42)(43).
Conversely, CNVs were more prevalent and importantly, demonstrated significant and biologically consistent correlations with cytotoxicity to corresponding targeted inhibitors. PPM1D, GNA13, and STAT family genesall on 17q and implicated in up to 70% of neuroblastomasshowed the strongest correlations of copy number and cytotoxicity to targeted inhibitors (44)(45)(46)(47)(48)(49)(50). This could be due to SCAs, which are frequent and associated with worse outcome in neuroblastoma (4,51), and regulate sensitivity to drug response in various cancers (52)(53)(54). Hence, with the known prevalence of SCAs in neuroblastoma, these associations between CNVs of key transcription factors and sensitivity to corresponding inhibitors suggests that differential copy number may be a potential quantitative predictor of drug response in neuroblastoma.
Significant associations between PI3K/Akt/mTOR pathway genes and their corresponding inhibitors were most frequently observed, accounting for half of the top 20 drug hits. PI3K/Akt/ mTOR pathway upregulation, with cytoplasmic phosphorylated Akt, occurs in 62% of neuroblastomas and is a promising therapeutic target (55,56), due to the central role played by PI3K activation in stabilizing MYCN, rendering neuroblastoma tumors sensitive to PI3K inhibition via downregulation of MYCN and corresponding tumor proliferation (27,57). BGT226 was the most potent PI3K inhibitor identified (83.6% cytotoxicity), but its effects on neuroblastoma remain understudied, and development has been discontinued following phase I/II trials in advanced adult solid tumors for lack of efficacy, inconsistent target inhibition and dose-limiting toxicities (58,59). PI3Ka inhibitor PIK-75 (78.9% cytotoxicity) stimulates antineuroblastoma activity by destabilizing MYCN and sensitizing tumor cells to anthracyclines (27,60), and indeed increased cytotoxicity with PIK-75 was associated with PPM1D and GNA13 gain, and DNMT3A, CBL, EED and ASXL2 loss. Next most potent was PI-103 (68.3% cytotoxicity), which sensitizes neuroblastoma cells to TRAIL-induced apoptosis (61), and primes them to cytotoxicity from anthracyclines (62).
Cyclin dependent kinases (CDKs) critically affect cell growth, proliferation and transcriptional regulation in neuroblastoma (63). Pan-CDK inhibitors flavopiridol and dinaciclibboth among our top hitsare in now early phase trials for neuroblastoma and advanced solid tumors (23)(24)(25). AT7519, a CDK2 inhibitor (70.0% cytotoxicity), has shown excellent preclinical anti-tumor responses against MYCN-amplified neuroblastoma (64), and is currently in phase II clinical trials for adults (NCT00390117) and children (NCT01627054). CDK4/6 expression is upregulated in 30% of neuroblastomas resulting in E2F overexpression and S and G2/M phase progression (32), and is another attractive therapeutic target. Palbociclib, a CDK4/6 inhibitor, is currently approved for advanced breast cancer and in phase III clinical trials (65); P276-00, a selective CDK4 inhibitor (79.7% cytotoxicity), has shown anti-proliferative effects against other cancers (66), but too has not been assessed in neuroblastoma.
Mixed patterns of gene-drug associations were observed in JAK/ STAT pathway genes, MAPK and ALK. JAK3 showed significant correlation with its corresponding inhibitor AT9283 and qualified as a drug hit, but not with other JAK inhibitors like ruxolitinib. While RAS-MAPK mutations are frequent in relapsed neuroblastoma, response to MAPK inhibitors (67), and BET and MEK inhibitors (68), have been inconsistent in vivo, possibly due to more complex relationships between receptor tyrosine kinase inhibitors and MAPK family genes. Similarly, we also did not find significant correlation between ALK copy number and cytotoxicity with crizotinib, but instead observed associations between ALK copy number and response to other kinase inhibitors. Despite the well-established prognostic significance of ALK mutations in neuroblastoma (69), a phase I/II clinical trial for pediatric patients with relapsed and refractory solid tumors yielded complete responses in only 1 out of 11 patients with known ALK mutations (NCT00939770) (70). Results are awaited from more recent trials for patients with ALK mutations investigating loratinib, a third generation ALK inhibitor (NCT03107988), ceritinib and ribociclib (NCT02780128), and ensartinib (NCT03155620). These findings may indicate variations in gene dosage-related sensitivity to response within drug classes.
This study was limited by the small number of cases, which may limit the robustness of the gene-drug correlations identified, and the inherent inability to test every drug hit in individual patients. Few of our shortlisted gene-drug associations could ultimately be verified against available datasets, limiting the confidence of our observations. However, our low incidence of verified gene-drug pairs is mirrored in the other drug repurposing studies mapping gene-drug associations in human disease (71). In future, corroboration of our findings in larger cohorts, and further functional validation will be required to verify the true clinical efficacy of our identified gene-drug pairs.
Associations between gene copy number and drug response could have been influenced by prior systemic treatment, or activation of alternate signaling pathways, though these are frequently encountered in relapse clinical cases where personalized therapeutic solutions are most often pursued. For example, the higher frequency of TSC2 mutants in our series could have contributed to the high rates of sensitivity of the PDCs to PI3K/ mTOR-targeted agents (72). The frequently observed massive amplifications of MYCN may have also affected responses to cell cycle agents. While MYCN is similarly recognized as an undruggable target due to the lack of appropriate drug binding surfaces on its DNA binding domain, BH3 mimetics and pan-BCl-2 inhibitors can exploit the Bcl-2 over-expression frequently accompanying MYCN amplification (73)(74)(75). Alternative treatment strategies have also been suggested that capitalize on its interactions with ALK and AURKA (76,77).
Treatment susceptibility prediction via genomic profiling and phenotypic assays also have their respective limitations: genomic profiling makes fundamental assumptions on correlations between genotype and phenotypic response to targeted inhibitors, while in vitro models may not fully recapitulate the original tumor drugresponse phenotype. Also, while OCCRA offers a more efficient and focused means to survey for meaningful genomic aberrations, other novel gene targets may be missed. Future studies could complement mutational profiling with gene expression results. Also, in future, drug sensitivity could be evaluated more accurately with three-dimensional organoid models, which may better represent the true tumor microenvironment and expression profile that may influence response to test agents (41,78).
In conclusion, integration of individualized genomic profiling and drug screening of patient-derived models of neuroblastoma showed that susceptibility to PI3K and cell cycle agents was significantly associated with copy number gains of PI3K and STAT family genes, particularly on 17q. This approach can facilitate CNV-based quantitative treatment response prediction in other mutationally quiet embryonal cancers like neuroblastoma.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by SingHealth Duke NUS Central Institutional Review Board (number 2014/2079). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.  MYCN amplification clustered together, of which 3 cases had concurrent copy number losses in 1p genes MPL, CSF3R, ARID1A, and ID3, and 19p genes MAP2K2, GNA11, CIC, CALR, JAK3, CRLF1, and SMARCA4 (NBL170120, NBL061218, NBL031219B) (solid box) and 1 case had multiple losses involving chromosome 9, 10, 15 and 16p genes (NBL061218) (dotted box).

AUTHOR CONTRIBUTIONS
Supplementary Figure S3 | Heatmap of z-scores of 418 drugs (columns) and 136 genes (rows) using unsupervised Euclidean clustering. Blue-red color scale highlights gene-drug pairs with 2.5% top-and bottom-most z-scores. Significant clusters of genes and drugs are annotated, and corresponding OCCRA panels for respective genes are indicated in grey and black on the right panel.
Supplementary Figure S4