Abstract
Osteosarcoma is a highly aggressive cancer for which treatment has remained essentially unchanged for more than 30 years. Osteosarcoma is characterized by widespread and recurrent somatic copy-number alterations (SCNA) and structural rearrangements. In contrast, few recurrent point mutations in protein-coding genes have been identified, suggesting that genes within SCNAs are key oncogenic drivers in this disease. SCNAs and structural rearrangements are highly heterogeneous across osteosarcoma cases, suggesting the need for a genome-informed approach to targeted therapy. To identify patient-specific candidate drivers, we used a simple heuristic based on degree and rank order of copy-number amplification (identified by whole-genome sequencing) and changes in gene expression as identified by RNA sequencing. Using patient-derived tumor xenografts, we demonstrate that targeting of patient-specific SCNAs leads to significant decrease in tumor burden, providing a road map for genome-informed treatment of osteosarcoma.
Osteosarcoma is treated with a chemotherapy regimen established 30 years ago. Although osteosarcoma is genomically complex, we hypothesized that tumor-specific dependencies could be identified within SCNAs. Using patient-derived tumor xenografts, we found a high degree of response for “genome-matched” therapies, demonstrating the utility of a targeted genome-informed approach.
This article is highlighted in the In This Issue feature, p. 1
Introduction
Osteosarcoma is the most common form of bone cancer in children and young adults. For patients who have metastatic disease at diagnosis or who relapse, the 5-year survival rate is below 30% (1–3). Established treatment regimens consisting of intensive multidrug therapy and surgical resection have significant short- and long-term toxicities and morbidities. A recent multinational effort to improve outcome in osteosarcoma by intensifying chemotherapy for high-risk groups failed to demonstrate improved survival (4), underscoring a critical need for new treatment strategies.
Osteosarcoma is characterized by significant somatic copy-number alteration (SCNA) and structural variation (SV) with few recurrent point mutations in protein-coding genes, with the exception of tumor suppressors RB1 and TP53 (5–7). This genomic landscape suggests that copy number (CN)–amplified genes within SCNAs may be critical drivers of disease progression and maintenance. Indeed, recent analysis of pan-cancer sequencing data suggests that the distinction between “C-class” (CN-driven) versus “M-class” (mutation-driven) cancers may be due to fundamentally different mechanisms of oncogenesis (8). Despite the prevalence of SCNAs in many cancers, efforts to develop targeted therapies and companion biomarkers have focused primarily on point mutations in protein-coding cancer genes. Comparatively less attention has been placed on identifying key tumor-specific vulnerabilities contained within SCNAs as a strategy for genome-informed targeted cancer therapy.
Previous studies have analyzed CN alterations in an attempt to identify candidate osteosarcoma driver oncogenes or tumor suppressors (9–15). These studies have focused on the identification of recurrent changes. However, genome-wide sequencing studies highlight the complexity of osteosarcoma genomes and underscore the likelihood that there are few common alterations across osteosarcoma samples (5, 6). Therefore, a strategy that seeks to identify key vulnerabilities in subclasses of osteosarcoma tumors may be more likely to succeed. Whether specific alterations present in SCNAs could harbor significant tumor-specific vulnerabilities is unknown, although prior work has established that CN changes can identify key vulnerabilities in other cancers (16, 17).
A major limitation in advancing genome-informed therapy for osteosarcoma is the lack of adequate models. Patient-derived tumor xenografts (PDTX) are increasingly utilized as model systems to test novel therapeutic approaches (18, 19). PDTX models for adult and pediatric cancers have been shown to closely recapitulate the genomic alterations present in the tumor of origin (20, 21). To nominate potential therapeutic opportunities for osteosarcoma, targetable cancer genes with high levels of amplification were identified by analyzing whole-genome sequencing (WGS) data from primary patient tumor samples. For a subset of these samples, PDTX models were also established, and a close correlation in SCNAs between the primary tumor and the corresponding PDTX was verified. The therapeutic efficacy of the subclass-specific targets was then tested in the PDTX models. Using this approach, significant responses were seen in all PDTXs tested, providing proof of principle for a genome-informed strategy to target a cancer that currently lacks adequate therapeutic options.
Results
Identification of Recurrent CN Changes in Targetable Cancer Genes by WGS of Osteosarcoma
WGS was performed on 30 tumor samples and corresponding germline DNA obtained from 23 patients. Samples were obtained from prechemotherapy diagnostic biopsies of primary tumors (n = 8), postchemotherapy resections (n = 8), or metastasis (n = 14), therefore representing the full spectrum of disease progression (Supplementary Table S1). All samples sequenced were reviewed by a pathologist and confirmed to be osteosarcoma. Regions with >70% tumor were microdissected for DNA isolation and sequencing. Average sequencing depth was 65X for tumors and 34X for germline samples (Supplementary Table S2). This dataset was expanded by merging with a previously published nonoverlapping WGS dataset of 33 samples from 31 patients (28 biopsies, 5 metastases; ref. 6). Consistent with previous reports, these osteosarcoma genomes were characterized by significant CN change and multiple structural rearrangements (Fig. 1A). When comparing two samples from the same patient, whole-genome CN profiles were largely consistent (Fig. 1B and C). For all samples, SVs and SCNAs were calculated. For samples with a matched germline sample, somatic nucleotide variants (SNV) were also calculated. Tumor purity was estimated using read-depth density and used to correct CN measurements (Supplementary Fig. S1A–S1D). Purity-adjusted CN values were classified as >4, >8, or >12 copies. Two hundred fifty-two genes annotated in public databases as druggable and clinically actionable were selected for further analysis (see Methods, Supplementary Fig. S1E and S1F, and Supplementary Tables S3 and S4). Genes with at least 2 samples indicating SCNA of greater than 8 copies were selected (Fig. 1D; Supplementary Table S4). Among the most commonly amplified genes (by patient) were MYC (39%) and CCNE1 (33%). Recurrent amplifications were also common in RAD21 (38%), VEGFA (23%), AURKB (13%), and CDK4 (11%) and others. Some of these alterations are likely co-occurring as they are in adjacent genomic locations. All but 2 patients had an amplification in at least one of the actionable genes. Deletions and mutations in a subset of known tumor suppressors were also evaluated with recurrent deletions, SVs, and SNVs found in TP53 (74%), RB1 (64%), PTEN (56%), and others. A landscape of somatic gains and losses across the evaluated osteosarcoma genomes indicates widespread alterations (Fig. 1E) but suggests that some of the most consistent SCNAs are those in actionable genes. Furthermore, although osteosarcoma genomes are highly complex, these results suggest the possibility of classifying osteosarcoma tumors based on the presence of specific SCNA in potential driver genes. In addition, comparison of the SCNA profile of patients with multiple samples suggests that SCNAs are highly stable between samples from the same patient (Supplementary Fig. S2). Although SCNAs were common, non-CN variations were found less often, resulting in low tumor mutational burden (Fig. 1D; Supplementary Table S1). Moderate- or high-impact SNVs and short indels were found in only 69 of the 252 druggable genes, of which only 20 were found recurrently altered. The most commonly altered genes were TP53 (12), ATRX (7), RB1 (6), and PRKDC (4), with the remaining alterations found in 3 or fewer patients. And although osteosarcoma tumors show a great deal of SV, SVs were found within only 81 of the druggable genes (30 recurrently). The most common genes with an SV were TP53 (29), LRP1B (14), RB1 (8), and FHIT (8). No potentially in-frame gene fusions were detected in more than 2 patients. Taking these observations into account, we then sought to determine whether analysis of SCNAs could have therapeutic relevance.
Genomic analysis of osteosarcoma (OS) and identification of recurrent SCNAs in primary tumors. A, Circos plot for indicated sample. CN gain and losses (outermost circle; red and blue, respectively), loss of heterozygosity (intermediate circle), and structural alterations (inner arcs) are shown. B and C, Genome-wide SCNA plot for diagnostic and resection samples from the same patient. D, Analysis of alterations across a cohort of 63 samples from 54 patients in actionable and druggable genes. CN gain was classified as >4 copies, >8 copies, or >12 copies. Actionable genes where at least 2 patients have gains of >8 copies were included in the top plot. Selected genes of interest (AKT1 and FOXM1) were also included. Losses for selected tumor-suppressor genes (TSG) were calculated and classified as <1.2 copies (minor) or <0.8 copies (major) and included in the bottom plot. SV truncations, gene fusions, and SNVs were calculated and included as indicated. Genes contained in segments with different CN states are annotated as a breakpoint (CNA). Top bar plot summarizes the CN gain/loss per sample. Right bar plot summarizes the CN gain/loss for a gene. Each column represents a single sample. Numbers to the left indicate percentage of alteration across patients (samples derived from the same patient were aggregated). In the top plot, only gains were included in the alteration percentages. For tumor suppressors, losses, SV, and SNVs were included. Purity estimates for each sample were calculated, and all CN gains/losses were adjusted accordingly. The genome-wide tumor mutation burden (TMB) was also calculated for each sample with a matched germline (number of variants across the genome per megabase). Samples collected from the same patient are labeled with the same letter in the “multisample” row. E, Combined genome-wide SCNA across all patients in the cohort (samples derived from multiple samples were combined). Percentages are of patients with gains and/or losses in 10 kb bins tiled across the genome. Gain/loss calculated and annotated as above. The loci of selected genes of interest are shown.
Genomic analysis of osteosarcoma (OS) and identification of recurrent SCNAs in primary tumors. A, Circos plot for indicated sample. CN gain and losses (outermost circle; red and blue, respectively), loss of heterozygosity (intermediate circle), and structural alterations (inner arcs) are shown. B and C, Genome-wide SCNA plot for diagnostic and resection samples from the same patient. D, Analysis of alterations across a cohort of 63 samples from 54 patients in actionable and druggable genes. CN gain was classified as >4 copies, >8 copies, or >12 copies. Actionable genes where at least 2 patients have gains of >8 copies were included in the top plot. Selected genes of interest (AKT1 and FOXM1) were also included. Losses for selected tumor-suppressor genes (TSG) were calculated and classified as <1.2 copies (minor) or <0.8 copies (major) and included in the bottom plot. SV truncations, gene fusions, and SNVs were calculated and included as indicated. Genes contained in segments with different CN states are annotated as a breakpoint (CNA). Top bar plot summarizes the CN gain/loss per sample. Right bar plot summarizes the CN gain/loss for a gene. Each column represents a single sample. Numbers to the left indicate percentage of alteration across patients (samples derived from the same patient were aggregated). In the top plot, only gains were included in the alteration percentages. For tumor suppressors, losses, SV, and SNVs were included. Purity estimates for each sample were calculated, and all CN gains/losses were adjusted accordingly. The genome-wide tumor mutation burden (TMB) was also calculated for each sample with a matched germline (number of variants across the genome per megabase). Samples collected from the same patient are labeled with the same letter in the “multisample” row. E, Combined genome-wide SCNA across all patients in the cohort (samples derived from multiple samples were combined). Percentages are of patients with gains and/or losses in 10 kb bins tiled across the genome. Gain/loss calculated and annotated as above. The loci of selected genes of interest are shown.
Generation and Analysis of Osteosarcoma PDTXs
The intertumor heterogeneity of osteosarcoma suggests that no single model system will be effective to test the therapeutic potential of specific drugs across all osteosarcomas. We therefore developed a panel of PDTXs by directly grafting tumor samples into immunocompromised mice (NSG) with the goal of developing tumor models that reflect the diversity of the human disease. Fifteen PDTX models were established (6 from pretreatment biopsies, 3 from postneoadjuvant therapy surgical resections, and 6 from metastatic sites; Supplementary Table S5). The histology of these PDTX models resembled the original tumor (Supplementary Fig. S3A). WGS analysis of each PDTX was performed (average 34X sequencing depth; Supplementary Table S6, and Supplementary Fig. S3B and S3C). Ten PDTXs had matching WGS analysis from the corresponding primary tumor, and one PDTX had matching WGS from a related primary tumor from the same patient. As with the patients with multiple samples, CN changes in these PDTX models were highly correlated with those in the corresponding primary tumor, and where available these correlations extend across multiple PDTX passages (Fig. 2A and B; Supplementary Fig. S2, and Supplementary Tables S4 and S7). For the primary tumors and PDTXs with available RNA-sequencing (RNA-seq) data, we compared the expression levels for key amplified genes of interest with the corresponding CN changes for the same genes using z-scores (Supplementary Tables S8 and S9; Supplementary Fig. S4A). We found a positive correlation in SCNAs and gene expression between the PDTX and their matched primary tumors, confirming that increased CN for a gene generally results in increased levels of gene expression (Supplementary Fig. S4B). This effect was particularly evident when genes showed extreme gains in CN (>12). These results suggest that osteosarcoma PDTX models reflect the genomic state of their primary tumor and could therefore be used to functionally test patient-specific targets.
Genomic analysis of osteosarcoma PDTX models and comparison with primary tumors. A, Scatter plots comparing SCNA changes in representative primary tumor versus PDTX pairs. CN represented as the normalized log2 ratio between somatic and germline samples. All protein-coding genes shown with actionable and druggable genes in red. B, Correlation matrix of CN across primary tumor versus exact matched PDTX pairs (all protein-coding genes). For all samples, the PDTX sample correlated best with the matched primary tumor. C, Genomic alterations across all PDTX samples for recurrent genes shown in Fig. 1D and other gene targets tested in PDTX models. All samples annotated and alterations calculated as in Fig. 1D. SCNAs targeted and tested in this study indicated with a white diamond. PDTXs derived from two separate samples from the same patient are indicated in the multisample row (letter matched with Fig. 1D). The relationship of the PDTX and the related primary tissue is shown. The companion primary tissue was from the exact tissue used for PDTX generation in all cases but one where a different stage primary tissue was used. In four instances, a comparable primary tissue was unavailable. D, CN for genes tested in this study for PDTX and exact matched primary samples. Additional primary and PDTX passages for selected samples are also shown. PDTX samples indicated with *, and PDTX samples from additional passages indicated with **. E, Schema for proposed treatment subclasses.
Genomic analysis of osteosarcoma PDTX models and comparison with primary tumors. A, Scatter plots comparing SCNA changes in representative primary tumor versus PDTX pairs. CN represented as the normalized log2 ratio between somatic and germline samples. All protein-coding genes shown with actionable and druggable genes in red. B, Correlation matrix of CN across primary tumor versus exact matched PDTX pairs (all protein-coding genes). For all samples, the PDTX sample correlated best with the matched primary tumor. C, Genomic alterations across all PDTX samples for recurrent genes shown in Fig. 1D and other gene targets tested in PDTX models. All samples annotated and alterations calculated as in Fig. 1D. SCNAs targeted and tested in this study indicated with a white diamond. PDTXs derived from two separate samples from the same patient are indicated in the multisample row (letter matched with Fig. 1D). The relationship of the PDTX and the related primary tissue is shown. The companion primary tissue was from the exact tissue used for PDTX generation in all cases but one where a different stage primary tissue was used. In four instances, a comparable primary tissue was unavailable. D, CN for genes tested in this study for PDTX and exact matched primary samples. Additional primary and PDTX passages for selected samples are also shown. PDTX samples indicated with *, and PDTX samples from additional passages indicated with **. E, Schema for proposed treatment subclasses.
To identify potential therapeutic targets in the PDTX models, we focused on SCNAs found in the same list of actionable genes recurrently amplified in primary tumors (Fig. 2C; Supplementary Table S3). In addition to reflecting the SCNA profile of the primary tumors from a whole-genome perspective (Supplementary Fig. S2), PDTXs were also consistent from a locus-specific perspective where gene targets with SCNA present in the primary tumors were maintained in the PDTXs (Fig. 2D; Supplementary Fig. S4A and S4C). The sole exception to this is OS128 where only a related tumor sample (metastasis) was available, whereas the PDTX was derived from a resection sample. The PDTX for this patient was noted to have a CDK4 amplification. For any individual cancer, multiple cancer-related genes could show SCNA and are thus potential drivers. We hypothesized that the degree of amplification would be indicative of selective pressure, reflecting its importance for tumor development. By prioritizing gene targets based on the degree of amplification and increases in gene expression, we identified several potentially targetable candidate driver pathways (Fig. 2E) and tested them using the PDTX models. The six candidate driver pathways tested were MYC gain, CDK4/FOXM1 gain, PTEN loss/AKT gain, AURKB gain, VEGFA gain, and CCNE1 gain, covering 93% (14 of 15) of PDTX models.
Genome-Informed Targeting of MYC
MYC was amplified in 39% of patients having at least >4 copies of MYC (Fig. 1D) and 3 patients having at least 12 copies of MYC. Two PDTX models with high levels of MYC amplification were identified (Fig. 3A). Both of these PDTXs were derived from posttreatment metastatic tumor resections, thus representing aggressive forms of osteosarcoma, and both PDTXs exhibited increased levels of MYC protein compared with PDTXs without MYC SCNA (Fig. 3B). Several approaches to targeting MYC-driven cancers have been described (22). Two of the most well-studied are inhibition of BRD4 and inhibition of CDK9 (23). Although neither approach is entirely specific to MYC, both have been shown in at least some tumors to preferentially affect tumors carrying MYC amplifications (24, 25). Treatment of both MYC-amplified PDTXs with the CDK inhibitor AT7519 (26) resulted in decreased tumor growth (Fig. 3C–H). AT7519 is a multi-CDK inhibitor that targets CDK1, 2, 4, 6, and 9 and has shown tolerability in phase I clinical trials (27, 28). Notably, in one PDTX model, single AT7519 treatment resulted in 3 of 8 tumors demonstrating greater than 30% tumor shrinkage (Fig. 3E). Similar results were obtained using two other multi-CDK inhibitors that also target CDK9 (ref. 29; Supplementary Fig. S5A).
MYC-amplified patient xenografts respond to CDK inhibition. A, Rank-ordered list of SCNAs identified 2 PDTXs with MYC amplification. B, Western blot across PDTXs with varying CNs of MYC. C, Growth curve for MYC-amplified PDTX (OS152) treated with AT7519 compared with vehicle control. D, Individual tumor volume at last time point in C. ***, P < 0.0001. Error bars, SEM. E, Waterfall plot of individual tumors in C. F, MYC-amplified PDTX (OS186) treated with AT7519 compared with vehicle control. G, Individual tumor volume at last time point in F. ***, P < 0.0001. Error bars, SEM. H, Waterfall plot of individual tumors in F. I, Western blot of short-term treated tumors OS152 PDTX, 4 doses of drug and sacrifice 4 hours after last dose. J, Representative IHC and quantitation of OS152 [10X fields of view (FOV)] after short-term treatment for CC3 as a measure of apoptosis, ***, P = 0.0006, and proliferation at end of study as measured by pH3 and quantitation, **, P = 0.003. Statistics calculated by the Student t test. Error bars, SD. Scale bar, 100 μm.
MYC-amplified patient xenografts respond to CDK inhibition. A, Rank-ordered list of SCNAs identified 2 PDTXs with MYC amplification. B, Western blot across PDTXs with varying CNs of MYC. C, Growth curve for MYC-amplified PDTX (OS152) treated with AT7519 compared with vehicle control. D, Individual tumor volume at last time point in C. ***, P < 0.0001. Error bars, SEM. E, Waterfall plot of individual tumors in C. F, MYC-amplified PDTX (OS186) treated with AT7519 compared with vehicle control. G, Individual tumor volume at last time point in F. ***, P < 0.0001. Error bars, SEM. H, Waterfall plot of individual tumors in F. I, Western blot of short-term treated tumors OS152 PDTX, 4 doses of drug and sacrifice 4 hours after last dose. J, Representative IHC and quantitation of OS152 [10X fields of view (FOV)] after short-term treatment for CC3 as a measure of apoptosis, ***, P = 0.0006, and proliferation at end of study as measured by pH3 and quantitation, **, P = 0.003. Statistics calculated by the Student t test. Error bars, SD. Scale bar, 100 μm.
To determine whether MYC was a potential target of AT7519 in osteosarcoma, we analyzed protein expression in tumors after short-term treatment. Levels of MYC protein as well as the MYC target MCL1 were significantly reduced (Fig. 3I). AT7519 also decreased phosphorylation of RNAPII (Fig. 3I; Supplementary Fig. S5B) and increased levels of cleaved PARP (Fig. 3I). Histologic analysis also showed increased cleaved caspase-3 (CC3) staining after short-term treatment, indicating that AT7519 induces apoptosis in some osteosarcoma tumor cells (Fig. 3J). At the end of study, there was a significant reduction in proliferation as marked by pH3 staining in treated tumors (Fig. 3J). Overall, these results indicate that AT7519 results in decreased proliferation and increased apoptosis in MYC-amplified osteosarcoma. Given prior literature on the mechanism of action of AT7519 and based on our own biochemical studies, we postulate that these results are most likely due to inhibition of MYC via CDK9, although we cannot exclude the possibility of a combined effect on several targets in addition to CDK9.
We then tested whether the BRD4 inhibitor JQ1, which has been shown previously to target some MYC-driven tumors (30), would be similarly effective in osteosarcoma tumors carrying MYC SCNAs. Surprisingly and in contrast to the effect of AT7519, JQ1 had only a marginal effect on tumor growth inhibition in a MYC-amplified PDTX (Supplementary Fig. S5C). Furthermore, JQ1 treatment did not cause reductions in MYC or MCL1 levels at the transcript or protein level in vivo (Supplementary Fig. S5D and S5E). To further assess differences between BRD4 inhibition and CDK inhibition in MYC-amplified osteosarcoma, we used a cell line generated from a PDTX (OS186; see Methods). We assessed viability after treatment with two different BRD4 inhibitors and compared this with three CDK inhibitors that have been shown to target CDK9. CDK inhibitors were more effective at decreasing viability compared with BRD4 inhibitors (Supplementary Fig. S5F). We also noted a decrease in pRNAPII (S2) with CDK inhibitor treatment but no decrease after JQ1 or iBET151 treatment. The decrease in pRNAPII (S2) was correlated with a decrease in MYC and canonical target MCL1 (Supplementary Fig. S5G). Notably, apoptosis was also increased after CDK inhibition but not after BRD4 inhibition (Supplementary Fig. S5H). These observations are consistent with previous reports that JQ1 acts independent of MYC inhibition in osteosarcoma (31).
To assess whether other “nonmatched” therapies (i.e., not matched to the SCNA in this PDTX) would also inhibit tumor growth, we treated a MYC-amplified PDTX with drugs targeting other pathways and observed no significant reduction in tumor growth with any of these agents (Supplementary Fig. S5C; see below for summary of all “matched” vs. “nonmatched” treatments). These results suggest that MYC SCNA analysis could be used to identify a subset of patients with osteosarcoma sensitive to MYC-directed therapy.
Genome-Informed Targeting of Cyclin E
Cyclin E (CCNE1) amplification is common in many cancers and is associated with poor prognosis and chemotherapy resistance (32–34). CCNE1 was amplified in 33% of assessed patients with osteosarcoma (Fig. 1D). Five PDTX models also carried CCNE1 amplification, and four had overexpression of Cyclin E protein compared with PDTXs without CCNE1 SCNA (Supplementary Fig. S6A and S6B; Fig. 3A). CDK2 inhibitors have been proposed as targeted therapy for Cyclin E–amplified tumors (35, 36). The CDK inhibitor dinaciclib (SCH 727965), which targets CDKs 1, 2, 5, and 9 (29), was therefore used to determine the effect of CDK inhibition in the context of CCNE1 amplification in osteosarcoma. Treatment of three different CCNE1-amplified PDTX models resulted in significant inhibition of tumor growth (Supplementary Fig. S6C–S6L). This result was confirmed in one PDTX on a subsequent passage (Supplementary Fig. S6F). Analysis after short-term treatment confirmed a modest but statistically significant increase in apoptosis as marked by CC3 staining (Supplementary Fig. S6M). At the end of study, we observed a decrease in proliferation as measured by pH3 staining (Supplementary Fig. S6M). Treatment with two “nonmatched” targeted agents (AZD1152 and palbociclib) led to only limited effects on tumor growth (Supplementary Fig. S6N). Thus, osteosarcoma tumor SCNAs with CCNE1 amplification may also be susceptible to therapy with multi-CDK inhibitors.
Genome-Informed Targeting of CDK4
CDK4 is a cyclin-dependent kinase that regulates cell-cycle progression during G1–S and is amplified in a variety of cancers including breast, head and neck, and lung (37). Palbociclib is a specific inhibitor of CDK4/6 and has been used successfully to treat breast cancer and other cancers (38, 39). CDK4 amplification was observed in 11% of patients, with 5 patients having gains of >12 copies (Fig. 1D). Two PDTXs with CDK4 amplifications were identified by rank-order analysis (Fig. 4A), and both demonstrated increased CDK4 gene and protein expression (Fig. 4B). When treated with palbociclib, both PDTXs exhibited significant growth arrest (Fig. 4C–H). To determine the early effects of drug treatment, tumors were analyzed after short-term treatment, and decreases in phospho-RB1, total RB1, and total FOXM1 were observed (Fig. 4I), consistent with on-target effects. Treatment with palbociclib led to a modest but statistically significant increase in apoptosis after short-term treatment as determined by CC3 staining (Fig. 4J). At end of study, pH3 was significantly decreased compared with vehicle, indicating a decrease in proliferation (Fig. 4J). These results suggest that CDK4 inhibitors could be effective as a targeted therapy in osteosarcoma tumors with CDK4 amplification.
CDK4- and FOXM1-amplified PDTXs respond to palbociclib treatment. A, Rank order of SCNA gains and losses in PDTX OS156 (left) and OS128 (right). B, Western blot of PDTXs with various CNs of CDK4. C, CDK4-amplified PDTX (OS156) growth curve treated with palbociclib compared with vehicle. D, Individual tumor volume at end of study in C,***, P < 0.0001. Error bars, SEM. E, Waterfall plot of individual tumors in C. F, PDTX (OS128) growth curve treated with palbociclib compared with vehicle. G, Individual tumor volume at end of study in F, **, P = 0.0035. Error bars, SEM. H, Waterfall plots of individual tumors in F. I, Western blot of short-term treated tumors from OS156 with palbociclib for 3 days. J, IHC of PDTX OS156 short-term treated tumors for CC3 and quantitation (per 10X FOV), *, P = 0.02, and at end of study for pH3 and quantitation (per 10X FOV), **, P = 0.002. Error bars, SD. Scale bar, 100 μm. K, Rank-ordered gene list of SCNA gains and losses for PDTX OS107. L, Western blot of PDTXs with various CNs of FOXM1. M, OS107 PDTX growth curve treated with palbociclib, tumor volume at treatment day 14, ***, P < 0.0001. N, Tumor volume at day 14 of treatment in M, ***, P < 0.0001. Error bars, SEM. O, Waterfall plot at 14 and 28 days of treatment from M. P, Western blot of PDTX OS107 short-term treated tumors with palbociclib compared with vehicle.
CDK4- and FOXM1-amplified PDTXs respond to palbociclib treatment. A, Rank order of SCNA gains and losses in PDTX OS156 (left) and OS128 (right). B, Western blot of PDTXs with various CNs of CDK4. C, CDK4-amplified PDTX (OS156) growth curve treated with palbociclib compared with vehicle. D, Individual tumor volume at end of study in C,***, P < 0.0001. Error bars, SEM. E, Waterfall plot of individual tumors in C. F, PDTX (OS128) growth curve treated with palbociclib compared with vehicle. G, Individual tumor volume at end of study in F, **, P = 0.0035. Error bars, SEM. H, Waterfall plots of individual tumors in F. I, Western blot of short-term treated tumors from OS156 with palbociclib for 3 days. J, IHC of PDTX OS156 short-term treated tumors for CC3 and quantitation (per 10X FOV), *, P = 0.02, and at end of study for pH3 and quantitation (per 10X FOV), **, P = 0.002. Error bars, SD. Scale bar, 100 μm. K, Rank-ordered gene list of SCNA gains and losses for PDTX OS107. L, Western blot of PDTXs with various CNs of FOXM1. M, OS107 PDTX growth curve treated with palbociclib, tumor volume at treatment day 14, ***, P < 0.0001. N, Tumor volume at day 14 of treatment in M, ***, P < 0.0001. Error bars, SEM. O, Waterfall plot at 14 and 28 days of treatment from M. P, Western blot of PDTX OS107 short-term treated tumors with palbociclib compared with vehicle.
In reviewing the SCNA data, we identified a PDTX with FOXM1 amplification (Fig. 4K). Three primary tumors also had amplification of FOXM1 (Fig. 1D). FOXM1 is involved in cell-cycle control and is amplified and overexpressed in several tumor types including basal-type breast cancer and diffuse large B-cell lymphoma (40). In addition, FOXM1 overexpression predicts poor outcome in a variety of cancers (41). FOXM1 is phosphorylated by CDK4/6, and this phosphorylation is required for FOXM1 activity and stabilization (42). Thus, CDK4/6 inhibitors could be predicted to inhibit FOXM1-amplified tumors in addition to tumors with CDK4 amplification. We confirmed overexpression of FOXM1 in the PDTX with FOXM1 amplification compared with other PDTXs by western blot (Fig. 4L). Treatment with palbociclib resulted in a decrease in tumor volume compared with initiation of drug treatment (Fig. 4M–O). Short-term palbociclib treatment led to decreased levels of pFOXM1 (Thr600), total FOXM1, pRB1 (807/811), and total RB1 (Fig. 4P). Palbociclib caused tumor shrinkage compared with slowed tumor growth or stasis as observed with CDK4-amplified PDTXs. Next, we assessed the effect of AT7519 (also targets CDK4/6) and dinaciclib (low specificity for CDK4/6). Only AT7519 had a similar effect on tumor growth as palbociclib (Supplementary Fig. S7A), suggesting the importance of specifically targeting CDK4 in these PDTXs. To directly assess the consequence of FOXM1 loss in our PDTX model, we used two shRNAs directed toward FOXM1 (Supplementary Fig. S7B) and implanted spin-infected cells to determine tumor growth in vivo. We observed a significant delay in tumor formation and growth with both shRNAs compared with control (Supplementary Fig. S7C). At end of study, tumors were excised and real-time analysis was preformed. FOXM1 transcript levels were higher at the end of study compared with cells used at the beginning to implant, suggesting that loss of FOXM1 repression by shRNA is required for tumors to grow and indicating an important role of FOXM1 expression in these tumors (Supplementary Fig. S7D).
CDK4 monophosphorylates RB1, and this interaction is critical for the mechanism of action of CDK4 inhibitors, such that tumors with loss of RB1 are generally insensitive to CDK4 inhibitors. However, FOXM1 is also directly phosphorylated by CDK4 (42), and thus the effect of CDK4 inhibitors on FOXM1-amplified tumors may be independent of RB1 status. To test this hypothesis, we knocked down expression of RB1 in a FOXM1-amplified PDTX (see Methods and Supplementary Fig. S7E and S7F). As expected, RB1 knockdown led to increased tumor growth (Supplementary Fig. S7G). Palbociclib treatment decreased tumor growth in both shRB1 and shGFP tumors, suggesting that in the context of FOXM1 amplification, palbociclib is effective even in the absence of RB1 (Supplementary Fig. S7G). In contrast, knockdown of RB1 in a CDK4-amplified PDTX cell line led to reduced effect of palbociclib treatment (Supplementary Fig. S7H and S7I). These data suggest that the mechanism of palbociclib response in CDK4-amplified osteosarcoma PDTX is through the canonical RB1 inhibition, whereas in the context of FOXM1 amplification, wild-type RB1 may not be required for response.
Genome-Informed Targeting of AURKB
AURKB is part of the chromosomal passenger complex and is a key regulator of mitosis (43). Prior studies have shown that some osteosarcoma tumors contain SCNAs involving AURKB (44) or overexpress AURKB and that treatment with inhibitors can lead to hyperploidy and apoptosis in vitro (45). Of the osteosarcoma samples in the cohort described above, 13% had a gain of AURKB (Fig. 1D). We observed higher protein levels of AURKB in an amplified PDTX by western blot and IHC (Fig. 4K; Supplementary Fig. S8A and S8B). To test whether inhibition of AURKB would be effective for this subclass of osteosarcoma, we used AZD1152, an AURKB-specific inhibitor (46). Treatment of these PDTXs resulted in a significant decrease in growth rate compared with vehicle (Supplementary Fig. S8C–S8E). Short-term AZD1152 treatment was associated with a decrease in pH3 staining and an increase in apoptosis as measured by CC3 (Supplementary Fig. S8F). Notably, 3 other PDTXs that did not contain SCNAs for AURKB and thus would be “nonmatched” to this drug exhibited no difference in growth rate compared with vehicle (Supplementary Figs. S5C, S6N, and S8G). These results suggest that AURKB inhibition may be specifically effective in osteosarcoma tumors with AURKB amplification.
Genome-Informed Targeting of PI3K–AKT–mTOR
Alteration of the PI3K–AKT–mTOR signaling pathway has also been reported in a subset of osteosarcoma (47). Loss of PTEN was seen in 56% of patients assessed by WGS. AKT1 amplification (>4 copies) was seen in 3 patients (Fig. 1D), and a single >12 copy gain of AKT2 was seen in 1 patient (data not shown). Taken together, PTEN loss or AKT1/2 gain may represent another targetable subclass of osteosarcoma. We identified two PDTX models with alterations in this pathway, one with biallelic PTEN loss and one with amplified AKT1 (Fig. 5A). PTEN protein levels were low in one PDTX with PTEN loss, and AKT1 protein was high in the PDTX with amplification of AKT1. In addition, western blot analysis indicated high levels of phosphorylation of ribosomal protein S6 in both of these samples, consistent with increased signaling downstream of AKT (Fig. 5B). To target this pathway, we used MK2206, a pan-AKT inhibitor (48). Significant reduction in tumor growth was seen with treatment in both PDTXs compared with control vehicle (Fig. 5C–H). Similar results were seen in a subsequent passage of one of these PDTXs (Supplementary Fig. S9A–S9C). Short-term treatment with MK2206 revealed a near-complete inhibition of pAKT1 (S473) and a decrease in the downstream target pS6RP (Fig. 5I). These early changes were not a result of a decrease in proliferation, as levels of PCNA remained similar between treatment and control groups (Fig. 5I). At this time point, we observed an increase in apoptosis as marked by CC3 staining (Fig. 5J), whereas at the end of study, a decrease in proliferation as indicated by pH3 staining (Fig. 5J) was observed. We also tested one PDTX with the mTOR inhibitor rapamycin and observed a similar decrease in tumor growth compared with vehicle as with MK2206 treatment, suggesting that either AKT or mTOR inhibition may be sufficient to target this pathway (Supplementary Fig. S9D–S9F).
AKT/PTEN pathway alterations respond to MK2206. A, Rank-ordered SCNA of gains and losses for OS525 (left) and OS052 (right). B, Western blot for PDTX with altered CN for PTEN and AKT. C, PTEN loss PDTX (OS052) treated with MK2206. D, Individual tumor volume at last time point in C, **, P = 0.005. Error bars, SEM. E, Waterfall plot of individual tumors in C. F, AKT1 gain PDTX (OS525) treated with MK2206. G, Individual tumor volume at end of study in F, **, P = 0.004. Error bars, SEM. H, Waterfall plots of individual tumors in F. I, Western blot of PDTX OS525 (AKT gain PDTX) short-term treatment of MK2206, 2 doses (M,W), and sacrificed 12 hours after last dose. J, IHC of CC3 short-term treatment and quantitation of OS525 (per 10X FOV), *, P = 0.015, and pH3 at end of study and quantitation (10X FOV), *, P = 0.023. Error bars, SD. Scale bar, 100 μm.
AKT/PTEN pathway alterations respond to MK2206. A, Rank-ordered SCNA of gains and losses for OS525 (left) and OS052 (right). B, Western blot for PDTX with altered CN for PTEN and AKT. C, PTEN loss PDTX (OS052) treated with MK2206. D, Individual tumor volume at last time point in C, **, P = 0.005. Error bars, SEM. E, Waterfall plot of individual tumors in C. F, AKT1 gain PDTX (OS525) treated with MK2206. G, Individual tumor volume at end of study in F, **, P = 0.004. Error bars, SEM. H, Waterfall plots of individual tumors in F. I, Western blot of PDTX OS525 (AKT gain PDTX) short-term treatment of MK2206, 2 doses (M,W), and sacrificed 12 hours after last dose. J, IHC of CC3 short-term treatment and quantitation of OS525 (per 10X FOV), *, P = 0.015, and pH3 at end of study and quantitation (10X FOV), *, P = 0.023. Error bars, SD. Scale bar, 100 μm.
To assess whether other targeted agents not directed to AKT or mTOR and therefore not matched to the genome of these tumors would have an effect, we treated one PDTX with either palbociclib (CDK4/6 inhibitor) or AZD1152 (AURKB inhibitor) and observed no statistically significant decrease in growth compared with vehicle (Supplementary Fig. S8G). These results suggest that SCNA can be used to identify a subset of osteosarcoma samples susceptible to treatment with AKT or mTOR inhibition.
Genome-Informed Targeting of the VEGF Pathway
Alterations in VEGFA CN and possible therapeutic implications in osteosarcoma have been previously reported (49). We observed VEGFA amplifications in 23% of patients (Fig. 1D). However, VEGFA protein across several PDTXs with varying CNs at this gene were similar. Nevertheless, a striking increase in VEGFR2 protein (also referred to as KDR), which is the main receptor for VEGFA, was observed in one VEGFA-amplified PDTX (Supplementary Figs. S6A and S10A). We hypothesize that VEGFR2 may be a better biomarker of VEGFA activation for VEGFA-amplified tumors given the secreted nature of the latter and its autocrine signaling mechanism of action (50). Sorafenib is a multikinase inhibitor active against several receptors including VEGFR2, RAF1, and BRAF which has been tested in several pediatric malignancies including osteosarcoma (51–53). We treated a VEGFA-amplified PDTX with sorafenib and observed a significant decrease in tumor volume compared with vehicle (Supplementary Fig. S10B–S10D). This study was repeated on a subsequent passage of the same PDTX, and a similar result was confirmed (Supplementary Fig. S10E–S10G). IHC at the end of study revealed a decrease in VEGFA protein and pERK, an established downstream target of VEGFR2 signaling in osteosarcoma (54). A significant decrease in proliferation as measured by pH3 was also observed (Supplementary Fig. S10H). Next, we sought to determine the specificity of sorafenib for VEGFA-amplified PDTX by testing in “nonmatched” PDTXs. We treated OS128 with sorafenib and observed no difference in growth between treatment and vehicle (Supplementary Fig. S10I). Together, these data suggest that inhibition of the VEGFA–VEGFR2 pathway could be used as a targeted agent for the subset of osteosarcoma tumors with VEGFA amplification.
Specificity of Targeted Therapies Based on SCNA
In this study, multiple PDTX models were treated with both “matched” and “nonmatched” therapy. Therefore, we sought to develop a general statistical framework to compare the efficacy of therapy targeted specifically to the genomic alterations in individual tumor samples. First, the % tumor growth inhibition index (TGI) was calculated for matched genome-informed therapy for each of the 10 PDTX models used and 12 matched therapies tested, taking into account that multiple genome-matched therapies could be indicated and were tested in a single PDTX (Fig. 6A; Supplementary Table S10). For example, CDK inhibitors can target multiple CDKs and thus would be considered matches for several PDTXs. Overall, CDK inhibition for MYC-amplified tumors was the most effective therapy, with TGIs of 86% to 97% for two different PDTXs using three different “matched” drugs. CDK4 inhibition in CDK4- or FOXM1-amplified PDTXs was the second most effective targeted therapy, with TGIs ranging from 61% to 111%. We observed a 79% TGI using sorafenib in one PDTX with VEGFA amplification and a 57% TGI with AZD1152 targeting an AURKB amplification. CDK inhibition by dinaciclib for CCNE1-amplified PDTXs had TGIs from 54% to 94%. Inhibition of AKT–PTEN pathway by the AKT inhibitor MK2206 was 61% to 67% for either PDTX tested (Fig. 6A).
Specificity of targeted therapies based on SCNA. A, Calculations of %TGI per PDTX and genome-matched targeted therapy tested for 10 PDTX and 12 targeted therapies. Asterisks indicate PDTXs that have multiple “matched” targeted therapies tested. B, Forest plot of mixed-effects model and pooled analysis calculation of growth rate of matched versus vehicle (P = 0.0058) of 10 PDTXs tested. C, Forest plot of mixed-effects model and pooled analysis calculation of growth rate of matched versus nonmatched (P = 0.0456) of 5 PDTXs.
Specificity of targeted therapies based on SCNA. A, Calculations of %TGI per PDTX and genome-matched targeted therapy tested for 10 PDTX and 12 targeted therapies. Asterisks indicate PDTXs that have multiple “matched” targeted therapies tested. B, Forest plot of mixed-effects model and pooled analysis calculation of growth rate of matched versus vehicle (P = 0.0058) of 10 PDTXs tested. C, Forest plot of mixed-effects model and pooled analysis calculation of growth rate of matched versus nonmatched (P = 0.0456) of 5 PDTXs.
To more rigorously compare the effectiveness of “matched” therapies, we used a meta-analytic framework to compare all the matched therapies across the PDTXs with “nonmatched” therapies. As tumors grew approximately linearly over time within each mouse on the log-transformed scale, we used a mixed-effects model and approximated the per mouse linear tumor growth trajectories allowing different rates of growth per mouse, per drug, and per PDTX, pooling the combination of multiple PDTX studies carried out over varying time periods (see Methods). Importantly, the mixed-effects model accommodates varying degrees of drug effect comparisons across PDTXs, including the possibility that the comparison may be significant in some PDTX models but not in others. First, we analyzed the growth effect of all matched therapies tested compared with vehicle control. Tumors from mice treated with matched drugs, on average, grew by 1.052-fold per day (corresponding to rate of growth of 0.0509 per day log scale), whereas tumors grew significantly faster in vehicle-treated mice, by 1.127-fold per day (corresponding to rate of growth of 0.1194 per day log scale; P = 0.0058; Fig. 6B; Supplementary Table S11). To evaluate the validity of our SCNA-targeted therapies, a similar meta-analysis comparing “matched” versus “nonmatched” therapies in the five PDTX models that were treated with both was performed. Overall, mice treated with nonmatched drugs grew by 1.1144 per day (corresponding to rate of growth of 0.1083 per day log scale), whereas mice that were treated with matched drugs grew significantly more slowly at a rate of 1.052-fold per day (corresponding to rate of growth of 0.0514 per day log scale; P = 0.0456; Fig. 6C; Supplementary Table S11). This daily change corresponds to a decrease in tumor size of 38.09% for the matched group over the course of 1 week and a decrease of 61.67% over 2 weeks. Moreover, the growth observed for the “nonmatched” drugs was not significantly slower compared with the vehicle treatment (P = 0.6884), underscoring the value of matched targeted therapy in these models.
Taken together, these results indicate that although osteosarcoma is highly heterogeneous, SCNAs represent a potentially novel avenue to define targeted, patient-specific therapies for this disease. Importantly, no single drug was universally beneficial to all samples, highlighting the importance for matched targeted therapies based on SCNA in the treatment of osteosarcoma.
Discussion
The last major advance in the treatment of osteosarcoma was made more than 30 years ago with the demonstration of clinical efficacy of a combined regimen including doxorubin, cisplatin, and methotrexate (55). Despite intensive efforts, no new therapeutic regimens have been found to improve survival for metastatic patients. In addition, no biomarkers to stratify patients to distinct therapeutic options currently exist. The complex genomic landscape of osteosarcoma suggests the need to address the heterogeneity of this disease in the design of future clinical trials. We analyzed genomic alterations in CN across a large cohort of osteosarcoma samples, combining previously unpublished WGS with a published cohort. By focusing on potentially actionable genes, we were able to reduce the complexity of the genomic landscape of osteosarcoma to identify alterations most likely to be of direct clinical relevance. Notably, this analysis also highlighted the extreme diversity among patients with osteosarcoma, as most genes were amplified in a subset of cases. Very few patients had no significant SCNA gains in druggable, clinically actionable genes, even though they exhibit the characteristic osteosarcoma pattern of genome-wide SCNA. However, each patient showed various degrees of loss of at least one canonical tumor suppressor, with TP53 being the most commonly altered (74%). Thus, it is likely that several distinct oncogenic drivers are responsible for the aggressive nature of this disease in individual patients, and a patient-specific approach to treating osteosarcoma is therefore likely to be warranted.
To develop models to assess the role of specific drivers within subsets of osteosarcoma, we assembled a collection of PDTX models and characterized them using WGS and RNA-seq. Genomic alterations were found to be relatively stable between primary tumors and their corresponding PDTXs, with only small variations between primary tumors and their established PDTXs. The genomic alterations were also found to be highly consistent across multiple passages in the PDTXs and derived cell lines. Together, these results suggest that these PDTX models may serve as faithful preclinical models to evaluate patient-specific therapies in osteosarcoma. However, a much larger collection of PDTX models will be needed to fully capture the full heterogeneity of osteosarcoma seen in the human disease, justifying continued efforts to generate such models. In one instance where a PDTX had a different SCNA in a targetable gene (CDK4 in OS128), the PDTX was derived from a different tissue (resection) than what was analyzed with WGS (metastasis). This suggests that although in most cases SCNAs present in the primary tumor are sustained in the PDTX, amplification of different drivers can occur during cancer progression. Using these PDTX models, we tested the potential effectiveness of genome-informed “matched” targeted therapies directed at putative driver genes with SCNAs. Taken together, these studies indicate a strong predictive value for presence of an SCNA in a PDTX model and drug response. Indeed, all of the “genome-matched” drugs had greater than or near 60% tumor growth inhibition in vivo, suggesting the potential for a positive effect in patients (56).
The ability to rapidly and inexpensively sequence tumor genomes has raised the possibility of “personalized” approaches for cancer therapy. The primary focus of most efforts has been on the targeting of point mutations in key oncogenes. However, many cancers and in particular many pediatric cancers have a low frequency of recurrent mutations in protein-coding genes (57, 58). Osteosarcoma in particular is characterized by significant SVs but few recurrent point mutations. Thus, we hypothesized that CN alterations rather than point mutations may be the dominant oncogenic mechanism. As CN changes are likely to select for amplicons highly supportive of oncogenesis, we reasoned that ranking CN changes across the genome would help identify key driver genes in individual samples and prioritize their therapeutic potential. In addition, using the rank order of gene targets helps to avoid issues of tumor purity, ploidy, and subclonality by selecting the targets most likely to be drivers.
In our studies we used only single drugs rather than combinations, as we wished to identify the target-specific effect of each drug in each PDTX. Nonetheless, it is well established that single-agent treatment most often leads to rapid development of resistance and is thus not an effective approach in most cancers. We expect that the next stage to advance genome-informed therapy in osteosarcoma will be the rational design of combination therapies, potentially with multiple targeted agents or a single targeted agent in concert with cytotoxic chemotherapy. For example, in some of the PDTXs used in this study, the same PDTX exhibited amplification of more than one of the targets tested. Combination therapies would therefore be expected to lead to increased response, provided this was tolerable in the preclinical models and ultimately in patients.
Many preclinical studies have evaluated the role of targeted therapies in osteosarcoma, with variable results. In the majority of cases, these studies have been carried out in established cell lines or cell line xenografts without specific attention to the genomic characteristics of the cell lines being used. Although PDTX models of osteosarcoma have been described (59–62), to our knowledge none have carried out a comprehensive analysis of PDTX models to test genome-informed therapy for this disease. A recent report described using a genomic approach similar to that proposed here to identify targeted therapies in osteosarcoma (63). In 2 patients evaluated, targeting of genes identified as altered using a DNA panel failed to lead to a clinical response. One difference between the integrative genomic analysis utilized here and in the prior study is that we used WGS and matched RNA-seq to identify SCNAs containing the most highly amplified genes, thereby potentially selecting the most likely drivers for that tumor. We propose that such a comprehensive approach may be particularly important for genomically complex diseases, such as osteosarcoma, and may help identify the most significant vulnerabilities for a given patient.
In summary, we report a comprehensive analysis of CN alterations and their therapeutic relevance for osteosarcoma, the most common bone malignancy in children and young adults and a disease for which traditional approaches to advance therapeutic discovery have been mostly unfruitful in the last 30 years. Our approach provides a blueprint for studies directed at genome-informed therapy and underscores the potential utility of precision medicine trials for this disease.
Methods
Sample Preparation and Sequencing
Samples.
Written informed consent was obtained from each patient (or from a parent in the case of patients <18 years of age) with recognized guidelines (Belmont Report) and institutional review board approval at each participating institution. All samples analyzed were reviewed by a pathologist at diagnosis and confirmed to be osteosarcoma. Samples were received fresh and a representative sample was reserved for PDTX generation, with remaining samples snap-frozen for DNA/RNA extraction.
Extraction of Nucleic Acids.
Snap-frozen samples were embedded in Tissue-Tek O.C.T. compound. Frozen section slides were cut on a cryostat with a section depth of 5 μm and stained with hematoxylin and eosin. Diagnosis was confirmed, and all samples were reevaluated for cellular content and quality by a pathologist (F.K. Hazard or S.-J. Cho). After visual inspection, samples with >70% tumor were macrodissected from the O.C.T. block at a depth of up to 5 mm for extraction. Samples were disrupted with a mortar and pestle under liquid nitrogen. DNA and RNA were co-extracted using the AllPrep Kit (Qiagen, 80204) with QIAshredder homogenization (Qiagen, 79654). Germline DNA (peripheral blood) was extracted using the DNeasy Blood and Tissue kit (Qiagen, 69504) according to the manufacturer's instructions. DNA was quantified using the Nanodrop 2000c (Thermo Fisher) and the QuBit High Sensitivity dsDNA assay (Thermo Fisher, Q32851). DNA integrity was quantified using the Genomic DNA Analysis ScreenTape (Agilent, 5067-5365) on the TapeStation 4200 (Agilent). RNA was quantified using Nanodrop 200c (Thermo Fisher) and QuBit High Sensitivity RNA assay (Thermo Fisher, Q32852). RNA was quantified using High Sensitivity RNA ScreenTape (Agilent, 5067-5579) on a TapeStation 4200 (Agilent).
WGS.
Libraries were made using the TruSeq Nano kit (Illumina, FC-121-4001) with a 350-bp insert as per the manufacturer's instructions. Libraries were made using 200 ng of input genomic DNA and sequenced to a depth of 30X (germline/PDTX) or an input of 400 ng DNA for sequencing depth of 60X (somatic). The majority of samples were sequenced on an Illumina HiSeq system with paired-end 2 × 150 bp reads by Macrogen, Inc. A subset of samples were sequenced to a depth of 50X (germline and tumor) on an Illumina HiSeq system with paired-end 2 × 100 bp reads by Illumina, Inc.
Raw DNA FASTQ data were preprocessed using NGSUtils (64) and aligned to the sex-specific GRCh38 reference genomes using bwa-mem (65). PDTX FASTQ data were aligned to both GRCh38 and GRCm38 human and mouse reference genomes separately using the same process. Reads were then classified as either human or mouse based upon the MAPQ and AS alignment scores using NGSUtils. Ambiguous reads with the same score for both organisms were discarded. Postprocessing of the alignment data included marking duplicate reads, indel realignment, and base-quality score recalibration using Picard (http://broadinstitute.github.io/picard) and GATK (66). For samples with a matched germline, somatic variants were called using MuTect2 (67). SVs were determined using Delly(V0.7.8) (68). SNV and SVs were annotated using VEP (69) and NGSUtils. Only SNVs rated as “moderate” or “high” impact by VEP were used for further analysis.
Purity Estimation.
The purity of samples was estimated by examining the density of the log2 ratios across bins. Purity-adjusted CN was then calculated for the range of 20% to 100% purity. Each purity level was evaluated to determine the difference between the purity-adjusted CN density peaks. The highest purity percentage resulting in a minimal distance between the density peaks and integer values was chosen as the best estimate. Ploidy and subclonality were not evaluated.
CN Analysis.
The number of reads aligned to the genome across 10-kb bins was determined for somatic and germline samples using NGSUtils and a normalized log2 ratio was calculated. For samples without a matched germline, a surrogate sex-matched germline sample was used instead. Bins where the germline counts were outliers (Tukey method) were excluded. Using the log2 ratio, bins were assembled into segments using DNACopy (70). The number of reads within segments was then recounted, and a log2 ratio was calculated for each segment. For genes that were not contained in a segment, the read counts across the whole gene region were determined and used. A purity-adjusted CN was then calculated as: ; where Pref is the reference ploidy, Pexp is the expected ploidy, v is the (normalized) log2 somatic/germline ratio, and p is the estimated purity. Purity-adjusted CN gain was classified as >4 copies, >8 copies, or >12 copies. Losses for selected tumor-suppressor genes were calculated and classified as <1.2 copies (minor) or <0.8 copies (major).
Clinically Actionable Genes.
A list of clinically actionable genes was determined by combining the gene lists from different cancer-related gene panels: FoundationOne and FoundationOne Heme (combined, Foundation Medicine, list accessed May 31, 2017), MSK-IMPACT (71), Mi-oncoseq (72), and UCSF 500 cancer gene panel (http://cancer.ucsf.edu/research/molecular-oncology/ucsf500). Genes present in more than one of these panels were included in the “actionable” gene list. Actionable genes that also had a drug interaction listed in the the Drug Gene Interaction Database (DGIdb; ref. 73; accessed June 13, 2018) were included in the actionable/druggable gene list, and this gene list was used to identify recurrent SCNA.
RNA-seq.
RNA-seq libraries were made using the TruSeq Stranded mRNA Kit (Illumina, RS-122-2101) with an input of 200 ng in accordance with the manufacturer's instructions. All manufacturer controls were used in preparation of the libraries. Libraries were quantified using the High Sensitivity DNA Kit (Agilent, 5067-4626) on the BioAnalyzer 2100 (Agilent). Sequencing was performed on an Illumina HiSeq system using chemistry for at least 2 × 75 bp reads at the Stanford Functional Genomics Facility.
RNA-seq FASTQ data were preprocessed using NGSUtils and aligned to the same sex-specific GRCh38 reference genomes as above using STAR (74). As above, PDTX RNA-seq data were aligned separately to human and mouse reference genomes and assigned to an organism based on alignment scores. Gene-level counts were obtained using the GENCODE v24 gene annotations (75) and NGSUtils. Gene expression was calculated as log2(CPM + 1).
Generation of PDTXs
All animal studies were done in accordance with the Institutional Animal Care and Use Committee at Stanford University or the University of California, San Francisco. Fresh patient samples were cut into 1 × 1 mm fragments and implanted either fresh or frozen in 90% FBS or 10% DMSO for later use. Prior to implantation, tumor fragments were dipped in Matrigel (Corning Matrigel #356234) and implanted in the subrenal capsule of NSG mice (Jackson Laboratory Strain #005557). Tumor growth was monitored for up to 1 year after implantation. Successfully implanted tumors were harvested at approximately 1 to 2 cm. A fragment was kept for histology, and the remainder was digested with collagenase and filtered through a 70 μm filter. For RNA/DNA isolation, cells were depleted of mouse stroma (Ter119, CD45, CD31, Mouse MHC class I) on a MACS column (Miltyni Biotech), followed by positive selection using Human HLA A,B,C (eBioscience) and sorting on a FACS Aria II. For subsequent passages and drug studies, cells were implanted subcutaneously in flanks of NSG mice (5 × 105 cells per flank) in 30 μL alpha MEM and 20 μL Matrigel (Corning).
Generation of PDTX Cell Lines
After successful generation of PDTXs, we generated a single-cell suspension and removed mouse stroma as above by depletion on a MACS column. We plated cells using standard DMEM supplemented with 10% FBS and 1% PSG. Cells were allowed to expand and sorted for human HLA positive to enrich for human osteosarcoma tumor cells, and this was performed twice to generate a pure population. Cell lines were submitted for karyotyping to confirm they were derived from osteosarcoma. We submitted osteosarcoma PDTX cell lines for low-pass WGS to confirm that the cell lines were derived from the patient and corresponding PDTX as listed.
Treatment of Mice with Targeted Compounds
When tumor cohorts of mice reached an average size of 100 mm3 per tumor, mice were stratified into treatment arms based on average tumor size per group. Mice were then dosed with drug or vehicle for 2 to 3 weeks. Mice were weighed at the beginning of study and periodically throughout drug treatment. Tumor volume was measured with digital calipers 3 to 4 times per week using the formula (length × width × width)/2 in mm3. Statistical significance at the end of study was calculated using a two-tailed, unpaired t test using Prism 6 software. Mice were dosed as follows: MK2206 (Selleckchem) 120 mg/kg MWF, p.o. in 30% captisol, rapamycin (Selleckchem) 4 mg/kg daily i.p., palbociclib (Pfizer CTP) 100 mg/kg daily p.o. in 50 mmol/L sodium lactate buffer, AT7519 (Selleckchem) 15 mg/kg daily i.p., flavopiridol (Selleckchem) 7.5 mg daily i.p., AZD1152 (Selleckchem) 25 mg/kg 4 consecutive days per week i.p., dinaciclib (Selleckchem) 30 mg/kg daily i.p., JQ1 (gift of Dr. James Bradner) 50 mg/kg daily i.p., and sorafenib (Selleckchem) 60 mg/kg 6 days per week p.o. All drugs delivered intraperitoneally were dissolved in DMSO and dosed with 10% B-hydroxycyclodextran (Sigma) unless otherwise stated.
Knockdown Studies in PDTX
pLKO shRNA constructs were purchased from Thermo. Lentivirus for each construct was generated by transfecting 293 cells (Invitrogen) with lipofectamine (Invitrogen), and viral supernatants were collected on days 2 and 3 after transfection, pooled, and stored at −80°C until use. Viral supernatants were then thawed and filtered through 0.45 μm PES filters and concentrated by ultracentrifugation for 2 hours at 24,000 rpm at 4°C. Viral pellets were resuspended on a platform rocker for 2 hours with approximately 500 μL fresh media. Cells were dissociated into a single-cell suspension using collagenase digestion buffer and filtered through a 70-μm filter and depleted for lineage (as above) on a MACS column. The resulting cell suspension was then plated at approximately 5 million cells per well of a 6-well plate and spin-infected with polybrene and virus in media at 1,500 rpm at room temperature for 30 minutes (Sorvall XRT centrifuge) then placed in the incubator. Media were changed the following day and 48 hours after infection, and puromycin (Invitrogen) was added (2 μg/mL) for 48 hours. After selection, the cells were gently trypsinized, filtered, and counted for viable cells. Cells were then implanted (as above) keeping the cell numbers consistent between study groups. Remaining cells were kept for confirmation of gene knockdown.
Western Blotting and IHC
Xenograft tumor fragments were stored at −80°C until use for western blot or fixed in formalin and embedded in paraffin for histology. Frozen tumors were thawed and mechanically disrupted using a 1.5 mL tube plastic homogenizer on ice and RIPA buffer. Protein quantitation was determined by BCA assay (ThermoFisher) and run on Bio-Rad 4%–15% gradient gels transferred to PVDF. Western blot analyses were performed using the LI-COR Odessey system and LI-COR blocking buffer or Bio-Rad Chemi Doc touch using LI-COR blocking for primary antibodies and 5% NFMP for secondaries in TBST. Primary antibodies for western blot and IHC were as follows: FOXM1, RB, pRB (807/811), PTEN, AKT, pAKT(S473), cleaved PARP, Cyclin E, CC3, and pH3 (Cell Signaling Technology); cMYC, MCL1, and AURKB (Abcam); RNAPII and RNAPII (S2) (Bethyl Labs); PCNA (SCBT and Cell Signaling Technology); and B-actin (Sigma). For both western blot and IHC analyses, short-term drug treatment consisted of 3 to 4 days of treatment, and mice were sacrificed 4 hours after the last dose. MK2206-treated mice were treated for 2 doses and sacrificed 12 hours after the last dose. Western quantitation was determined using LI-COR secondary antibodies and Image J quantitation software. For long-term drug treatment for pH3 staining, mice were sacrificed 24 to 48 hours after the last dose. CC3 and pH3 IHC staining were performed on two tumor specimens per PDTX treatment condition. Representative 10X fields of view (FOV) were imaged, avoiding areas of necrosis and low cellularity. Five to 8 FOVs were analyzed per condition. Data are displayed as mean per 10X FOV, and error bars are SD. All quantitation was done in a blinded fashion.
Statistical Analysis of PDTX Cohorts
With linear growth trajectories well approximating within-mouse tumor growth on the log-transformed scale, we log-transformed the tumor volume data and used a linear mixed-effects model to compare rates of growth of different treatment drugs. The per-PDTX analysis included random intercepts and random slopes to account for the within-mouse correlation among the longitudinal tumor volume measurements and differential rates of growth per mouse, respectively. In contrast to the per-PDTX analysis where different drug effects within each PDTX were included as fixed slope effects, the pooled mixed-effects analysis included random per-drug, per- PDTX slopes to accommodate varying degrees of drug effects within PDTXs and differential rates of growth per PDTX. The additionally included random slopes allowed the possibility that the comparison of “matched” versus “nonmatched” drugs might be significant in some PDTX models but not in others. All computation was conducted using SAS 9.3.
Accession Codes
Whole-genome and RNA-seq data that support the findings of this study have been deposited in the European Genome-phenome Archive (EGA; EGAS00001003201—whole-genome and RNA-seq analysis of pediatric osteosarcoma). Previously published whole-genome sequences that support the findings in this study are available from the EGA (EGAD00001000159, EGAD00001001053, doi:10.1016/j.celrep.2014.03.003).
Disclosure of Potential Conflicts of Interest
N. Marina is Director at Five Prime Therapeutics, Inc. D.G. Mohler has served as an expert witness for the University of California, San Francisco. D.S. Hawkins is a consultant/advisory board member for Loxo Oncology, Bayer, Bristol-Myers Squibb, and Celgene. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: L.C. Sayles, M.R. Breese, S.L. Spunt, E.A. Sweet-Cordero
Development of methodology: L.C. Sayles, M.R. Breese, A.L. Koehne, E.A. Sweet-Cordero
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L.C. Sayles, A.L. Koehne, S.G. Leung, H.-Y. Liu, A. Spillinger, A.T. Shah, F.K. Hazard, S.L. Spunt, N. Marina, G.E. Kim, S.-J. Cho, R.S. Avedian, D.G. Mohler, S.G. DuBois, D.S. Hawkins, E.A. Sweet-Cordero
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L.C. Sayles, M.R. Breese, A.L. Koehne, A.G. Lee, B. Tanasa, K. Straessler, N. Marina, M.-O. Kim, E.A. Sweet-Cordero
Writing, review, and/or revision of the manuscript: L.C. Sayles, M.R. Breese, A.L. Koehne, S.G. Leung, F.K. Hazard, S.L. Spunt, R.S. Avedian, D.G. Mohler, M.-O. Kim, S.G. DuBois, D.S. Hawkins, E.A. Sweet-Cordero
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L.C. Sayles, M.R. Breese, A.L. Koehne, S.G. Leung, A.G. Lee, H.-Y. Liu, A. Spillinger, E.A. Sweet-Cordero
Study supervision: L.C. Sayles, E.A. Sweet-Cordero
Other (worked with E.A. Sweet-Cordero on the initial concept and participated in enrolling patients and looking at the data): N. Marina
Other (tumor board review of the cases): D.G. Mohler
Acknowledgments
This work was funded in part by the Lucile Packard Children's Hospital (LPCH) of Stanford University. We wish to especially thank the leadership of LPCH (Chris Dawes) and the Department of Pediatrics (Hugh O'Brodovich) for their generous support of the Stanford Pediatric Cancer Genomics Initiative. Further support for this effort was provided by generous funds from the Helen Diller Comprehensive Cancer Center and Dr. Alan Ashworth. E.A. Sweet-Cordero is a Benioff Professor of Children's Health, which provided him with funds to carry out this work. We thank all members of the Sweet-Cordero lab for helpful discussion. We thank also Dr. Julien Sage (Stanford University) and Dr. Erin Breese (Cincinnati Children's Hospital Department of Pediatrics) for thoughtful discussion regarding experimental and study design. We also thank the Stanford Tissue Bank and UCSF tissue bank for help in collecting patient samples. Most importantly, we would like to express the deepest appreciation to all the patients and their families who generously donated tissue for this research.
This work used resources from the Stanford Genetics Bioinformatics Service Center and the Stanford Functional Genomics Facility, including support from an NIH S10 Shared Instrumentation Grant (S10OD018220). E.A. Sweet-Cordero was supported by a Morgridge Scholar Grant and by the Children's Health Research Institute (Stanford). Funding was also provided to E.A. Sweet-Cordero by the Hyundai Hope on Wheels Foundation, the St. Baldrick's Team Clarkie Award, the CureSearch Foundation, and the Sarcoma Foundation of America. This study was funded in part by the HDFCCC Laboratory for Cell Analysis Shared Resource Facility through a grant from the NIH (P30CA082103). A.L. Koehne was funded by training grant 2T32OD011121-11 from the NIH. K. Straessler was funded by training award R25 CA180993-02 from the NIH (Cancer Systems Biology Scholars Program, Stanford University).