Sensitivity to BUB1B Inhibition Defines an Alternative Classification of Glioblastoma

Glioblastoma multiforme (GBM) is the most aggressive common form of brain cancer in adults. GBM is characterized by poor survival (i.e., overall median survival time is one year; ref. 1), remarkably high tumor heterogeneity (2), and lack of effective therapies (3). Recent studies have revealed that GBM can be divided into distinct subtypes with prognostic relevance based on differences in gene expression (4, 5), somatic mutations (5, 6), DNA methylation (7), and DNA copy number (8). Although genomic analyses have defined GBM subclasses based on transcriptional and epigenetic regulations, the subclasses have failed to suggest effective subclass-specific treatment strategies, with the exception of the subclass based on promoter methylation status of the DNA repair gene MGMT. In this case, a subset of GBM tumors with MGMT promoter methylated and transcriptionally downregulated are more likely to benefit from the addition of temozolomide to radiotherapy (9). However, the majority of GBM patients show very little benefit from surgery, radiation, and temozolomide (i.e., standard-of-care therapy). Thus, to achieve better outcomes in the clinic, we need a better patient stratification and more in-depth understanding of the biology of these tumors.

Both adult and pediatric GBM tumors appear to be hierarchically organized, suggestive of a cancer stem cell origin (10–13). Consistent with this notion, tumor-initiating GBM stem cells (GSC) have recently been isolated that retain the development potential and specific genetic alterations found in the patient's tumor (10, 11, 14, 15). GSCs retain the specific genetic and epigenetic signatures found in the original tumor (10, 11, 15), and give rise to tumors with GBM patient-specific molecular signature and histologic features when implanted into the cortex of rodents (11, 15). Importantly, treatments targeting GSCs may be more effective because GSCs are radioresistant and chemoresistant due to their preferential activation of the DNA damage response, which eventually results in tumor recurrence (16). Therefore, the use of patient-derived GSC isolates can allow investigation of the molecular characteristics of subpopulation of tumors, and potentially develop more effective treatments.

Recently, we performed shRNA kinome screens in GSC isolates and non-neoplastic neural progenitor controls for genes required for GSC expansion (17). Combing these results with a GBM bionetwork created from patient molecular profiles, we identified BUB1B as the top GSC-specific hit. BUB1B encodes the highly conserved Bub1-like pseudokinase, BubR1, which possesses multiple functional domains implicated in mitotic checkpoint control, mitotic timing, and regulating KT-MT attachment (18). These include: N- and C-terminal KEN box domains required for Cdc20 binding and APC inhibition (19); a C-terminal kinase domain required for protein stability (20); and a GLEBS domain necessary for kinetochore localization during mitosis (21). Although BUB1B is essential for mammalian development (22), its essential function is contained solely within the N-terminal KEN box (23), which enables BubR1 to act as a pseudo-substrate inhibitor of APC/C^Cdc20 during G₂ and preanaphase mitosis, preventing premature anaphase onset (23). In contrast, we observed that in approximately 60% of GSCs, Ras-transformed cells, and HeLa cells, its GLEBs domain becomes essential for viability to promote kinetochore–microtubule attachment (17). Mechanistic experiments demonstrated that oncogenic Ras signaling triggers alterations in kinetochore regulation, resulting in added GLEB domain requirement and the primary reason we observe differentially sensitivity to BUB1B knockdown (17, 24).

BUB1B inhibition–sensitive (BUB1B^S) cells invariably have shorter metaphase interkinetochore distances (IKD), or shorter average distances between sister kinetochores during mitosis when stable end-on microtubule attachments have formed (17, 24). This serves as an indirect measure of the pulling forces generated by dynamic microtubules bound to kinetochores, such that stronger attachments lead to longer IKDs and weaker attachments produce shorter IKDs (24). Although IKDs are reliable predictors of BUB1B^R/S and in theory could be used to predict tumor sensitivity to BUB1B inhibition, in practice, taking IKD measurements is laborious and time consuming, requiring confocal microscopic z-sectioning of mitotic cells, and unlikely to be useful to “type” tumor samples. Here, we instead used gene expression signatures associated with BUB1B^S and BUB1B^R cells to create a computational framework for predicting BUB1B^S/R status for GSCs, which examines differential expression of 838 genes. We then applied our classifier to existing GBM patient tumor data, revealing that BUB1B^S GBM tumor patients have significantly worse prognoses, regardless of their molecular subtype. Examination of results from genome-wide RNAi screens in BUB1B^S/R isolates revealed differential reliance of genes enriched in mitotic cell cycle, and chromosome segregation, and suggested potential therapeutic treatments of BUB1B^S cells. Applying the classifier to drug sensitivity profiles confirmed BUB1B^S cells to be more sensitive to type I and II topoisomerase inhibitors as well as histone deacetylase (HDAC) inhibitors, and additionally uncovered BUB1B^S cells to be more sensitive to Raf inhibitors. Taken together, our results suggest that BUB1B^S/R classifier is useful for predicting tumor aggressiveness and therapeutic responses of GBM tumors in additional to their sensitivity to BUB1B-GLEBs domain inhibition.

We collected 18 GSC cell lines that have been previously published, one neuronal stem cell (NSC), the normal human astrocyte cell line (i.e., NHA, StemCell Technologies), and RasV12-transformed NHA (Russell Pieper, University of California San Francisco, San Francisco, CA; references in Supplementary Table S1). Three biological replicates for each cell line were sequenced using Illumina platform. Twenty million reads per sample on average were mapped to transcripts. We removed one GSC cell line (i.e., G25) from our further analysis due to its low number of mapped sequences for all three replicates (Supplementary Table S1). Sequencing data can be accessed at NCBI Gene Expression Omnibus (GEO) under GSE89623. In addition, we collected nine GBM tumor samples from Mount Sinai Hospital (September 2014–November 2015) and derived GSC isolates from each tumor sample. The tumor tissue samples and derived GSCs were sequenced using Illumina HiSeq2500. On average, 42 million reads per sample were mapped to transcripts. Sequencing data can be accessed at GEO under accession number GSE94874. The cell culture protocol and preprocessing of RNA-seq data are detailed in Supplementary Materials and Methods.

We first identified genes whose expression levels were significantly associated with BUB1B^S/R status. We designed a multiple linear model to explain the variation in mRNA expression level for each gene in terms of BUB1B^S/R status and cell types as covariate. The associated genes were defined by their significance level of the coefficients for BUB1B^S/R status. On the basis of these genes associated with BUB1B^S/R status, we trained an elastic net classifier, a regularization method to prevent overfitting problems in a generalized linear model (25). Then, the elastic net classifier was applied to gene expression profiles of new samples to predict BUB1B^S/R status. See Supplementary Materials and Methods for details.

We expected that molecular profiles of GSCs and GBM tumor cells were similar. Thus, we extracted molecular profile of GBM tumor cells from expression profiles of bulk tumor samples by following three steps. First, we determined gene expression profile of each normal brain cell type and GSC. The single-cell expression levels of six different brain cell types (26) and our GSC data were merged, normalized (27), and averaged to define the profile of each cell type. Second, we deconvoluted bulk tissue expression data by using CIBERSORT (https://cibersort.stanford.edu/) with the expression levels of normal brain cells and GSC, resulting in relative proportion of six normal cell types and GSC. Finally, we extracted molecular profile of GBM tumor cells from the measurements of bulk tumor samples based on relative proportion and gene expression profile of each cell type. The extracted molecular profiles of GBM tumor cells were used to predict BUB1B^S/R status. See details in Supplementary Materials and Methods.

We characterized BUB1B^S/R subtypes in term of the clinical survival rate for three datasets including The Cancer Genome Atlas (TCGA; ref. 28), Gravendeel and colleagues dataset (29), and Joo and colleagues dataset (see Supplementary Materials and Methods for download and preprocessing of datasets; ref. 30). First, the overall survival was associated with the predictive values based on the elastic net classifier by fitting Cox proportional hazards model (31). In addition, the Kaplan–Meier curves for BUB1B^S and BUB1B^R subtypes that are defined as samples with the highest and lowest 25% predicted values were compared (32). The test of equality for survival distributions was performed using the log-rank test.

We performed genome-wide shRNA screen and Barcode array analysis for three GSC cells and one NSC cell (CB660) as described previously (33). We detected each isolate's candidate targets as shRNAs underrepresented in GSCs relative to NSCs. On the basis of these candidate targets, we identified genes only required for the expansion of BUB1B^S and BUB1B^R GSCs. These genes were used to define BUB1B^S/R lethal subnetworks. See Supplementary Materials and Methods for details.

With genes in the two large BUB1B^S-specific lethal subnetworks (see Supplementary Materials and Methods for details), we determined whether these genes were upregulated or downregulated in BUB1B^S GSCs by comparing their gene expression levels between BUB1B^S GSCs and BUB1B^R GSCs (Supplementary Fig. S7). Then, we investigated whether these genes were perturbed by drug treatments by using the tool (http://www.broadinstitute.org/cmap/) from the connectivity map (34). A negative enrichment score from the output of the tool represents the treatment of the drug might inhibit the BUB1B^S lethal genes (therapeutic to BUB1B^S group), resulting in selectively killing of BUB1B^S GSCs.

We predicted BUB1B^S/R status of glioma cell lines from Cancer Cell Line Encyclopedia (CCLE) and Cancer Genome Project (CGP) database (35, 36). We collected gene expression levels of 59 and 43 glioma cell lines from CCLE and CGP database, respectively, and applied our classifier on gene expression profiles to predict BUB1B^S/R status. Next, we explored the drug sensitivity measurements of cell lines to identify drug that has different sensitivity on BUB1B^S/R subtypes. We used the concentration at which the drug response reached an absolute inhibition of 50% (IC₅₀) and an “activity area” to capture simultaneously the efficacy and potency of a drug (36). See details in Supplementary Materials and Methods. For each measurement, we calculated the correlation between drug sensitivity measurements and predictive values that determine BUB1B^S/R status across glioma cell lines to identify drugs that have different sensitivities between BUB1B^S and BUB1B^R cell lines.

The effects of etoposide and irinotecan on GSC isolates were determined by treating them with various concentrations of each drug for incubation period 72 hours at concentrations ranging from 1 μmol/L to 4 mmol/L for etoposide and 0.01 μmol/L to 1 mmol/L for irinotecan. See details in Supplementary Materials and Methods.

On the basis of our previous results, where we demonstrated that GBM tumor isolates are either BUB1B inhibition–sensitive (BUB1B^S) or resistant (BUB1B^R Fig. 1A; ref. 17), we hypothesized that BUB1B^S/R could be captured by genome-wide gene expression analysis. To this end, we analyzed the RNA-seq profiles of 18 GSCs as well as astrocyte and neural progenitor cells, classifying them by BUB1B^R/S status. BUB1B^R/S status was determined by sensitivity to LV-shBUB1B and/or measurement of inter-kinetochore distance at metaphase (Supplementary Table S1 and Supplementary Fig. S1; refs. 17, 24). This analysis revealed nine BUB1B^S isolates, four BUB1B^R isolates, and five isolates that had not been experimentally tested for BUB1B^R/S (Fig. 1B, one isolate was removed from further analysis due to low sequencing quality). In addition, the neural progenitors and the normal human astrocyte cell line (i.e., NHA) were BUB1B^R, on the other hand, the RAS-treated astrocyte was BUB1B^S (Fig. 1B). To identify a BUB1B^R/S gene expression differences, we used a linear model to fit the variation in mRNA expression level for each gene in terms of BUB1B^R/S. We identified 838 genes whose expression is associated with BUB1B^S/R at FDR <1% corresponding to P < 10^−3.9. Gene set enrichment analysis (37) of the 838 BUB1B^S/R-associated genes revealed that they were significantly enriched for cell-cycle–related pathways such as the mitotic cell cycle (Supplementary Fig. S2).

We next developed an elastic net–based classifier (see Materials and Methods for details; ref. 25) to predict sensitivity to BUB1B inhibition using the 838 BUB1B^S/R–associated genes (Fig. 1C and D). The cross-validation results indicate that the elastic net classifier can achieve low mean squared error (Supplementary Fig. S3). We applied the classifier to four GSC isolates of unknown BUB1B^R/S status, including G144, 1406, 448T, and 559T (i.e., white in Fig. 1B). The predictive values of new samples obtained by fitting elastic net classifier on its gene expression profiles were used to determine BUB1B^S/R (Fig. 1C and D). We predicted that G144 and 1406 would be BUB1B^R, whereas 448T and 559T would be BUB1B^S. Consistent with these predictions, IKD measurements of metaphase cells confirmed that both G144 and 1406 have long IKDs (i.e., BUB1B^R), whereas 448T and 559T have short IKDs (i.e., BUB1B^S; Fig. 1E). Thus, these results demonstrate that our classifier can accurately predict BUB1B^R/S status.

We next tried to apply our classifier to GBM patient tumors to investigate the relationship between BUB1B^R/S status and GBM patient prognosis. However, tumor samples are complex mixtures of malignant and normal tissue cells; molecular signal might be averaged out between normal and cancer cells for tumor tissue. To investigate the relationship between tumor tissues and GSCs, we collected nine pairs of GSCs and the corresponding tumor tissues from which they were derived, and measured gene expression levels. When considering the gene expression, GSCs, which were purely neoplastic, were clearly separated from the corresponding tumor tissues (Supplementary Fig. S4A). This indicates that molecular signature of tumor tissues likely contained features of both GBM tumor cells and nontumor cells. The predicted value of the BUB1B^R/S status based on gene expression levels of tumor tissues and corresponding derived GSCs were not consistent (Supplementary Fig. S4B). A reasonable explanation is that the normal, nontransformed cells are predicted to be BUB1B^R (17), so that the prediction of the BUB1B^R/S status based on the mixed molecular signals of tumor tissue expression is confounded by proportion of normal cells in the tumor tissues. These results indicate needs of decomposition of the molecular signals of GBM tumor cells from normal brain cells in bulk tumor tissues to accurately predict the BUB1B^R/S status of tumor cells.

To extract the molecular signal of GBM tumor cells, we decomposed tumor tissue expression into expression of diverse cell types in the human brain (i.e., astrocytes, endothelial, microglia, neurons, oligodendrocytes, and oligodendrocyte progenitor cells; ref. 26) and our GSCs using CIBERSORT (38), and estimated relative fraction of diverse brain normal cell types and GBM tumor cell in each tissue sample (Supplementary Table S2; see Methods and Materials for details). The relative fraction of GBM tumor cells varied from 19% to 55% (Supplementary Fig. S4C), and indeed, the predicted value of the BUB1B^R/S status based on tissue expression levels was highly associated with the predicted proportion of GBM tumor cells in the mixture (Supplementary Fig. S4C). The predicted cancer cell fractions were consistent with histologic results. For example, the tumor tissue sample with the lowest GBM tumor cell fractions (i.e., 9217A), the histologic staining of tumor tissue indicated lower percentage of cancer cells than others (Supplementary Fig. S4D). On the basis of the relative fraction and molecular profile of each cell type, we extracted the molecular signals of GBM tumor cells from bulk tumor tissue expression levels (see Materials and Methods for details). The deconvoluted expression profiles of GBM tumor cells were used to predict BUB1B^S/R, and resulted in the similar order of predictive values with those based on GSC expression levels (Supplementary Fig. S4E), Taken together, the deconvolution of bulk tumor tissue expression levels to extract molecular signatures of GBM tumor cells is the essential step to predict accurate BUB1B^S/R based on tumor tissue expression levels.

We characterized BUB1B^S/R subtypes in term of the clinical survival rate for three datasets for GBM: TCGA (28, 39), Gravendeel and colleagues dataset (29), and Joo and colleagues dataset (30). We applied the deconvolution procedure to these tissue gene expression datasets, and predicted BUB1B^S/R status based on extracted molecular profiles of GBM tumor cells. The overall survival was significantly associated with the predicted sensitivity value based on the elastic net classifier by fitting Cox proportional hazards model (31) for TCGA dataset–based all three platforms (e.g., Agilent platform P < 0.040 based on Cox regression; Table 1). The BUB1B^S subtype showed worse survival rate than the BUB1B^R subtype (Fig. 2A), and the difference assessed by log-rank test (32) was significant. The application to the Gravendeel and colleagues dataset and Joo and colleagues dataset yielded similar results that BUB1B^S/R status was associated with survival rate (e.g., Gravendeel and colleagues P < 0.0097 and Joo and colleagues P < 0.003 based on Cox regression; Fig. 2A; Table 1) with poor prognosis for patients in the BUB1B^S subtype. Note that the predicted sensitivity scores based on tissue expression levels without deconvolution were not consistently associated with survival rate (Table 1), underscoring the importance of extracting expression signatures of GBM tumor cells before applying our classification method. To be able to perform more robust exploration of the relationship of BUB1B^S/R subtype to survival rate, we incorporated the other clinical information into predicted sensitivity values. Both gender and the methylation status of MGMT promoter that is prognostic to temozolomide treatment (9) were found not to be significantly associated with predicted sensitivity values (Supplementary Fig. S5A) for any datasets, on the other hand, age at diagnosis was found to be significantly associated with survival rates for some of datasets (Supplementary Fig. S5B). Even though the inclusion of age as a covariate decreased the significance level of predicted sensitivity value to BUB1B inhibition in Cox regression model for most datasets, the BUB1B^S/R subtypes were still significantly associated with survival rate for most datasets (Supplementary Table S3).

We further examined the other clinical phenotypes associated with BUB1B^S. In particular, we analyzed orthotopic glioblastoma xenograft data from 58 GBM patients (30) with expression profiles from GBM patients' primary tumors and parallel in vivo xenograft tumors. We examined the invasiveness of GBM patients as well as xenograft inherited from these GBM patients, and found suggestive association between our predicted sensitivity values and invasiveness of patient's tumor as well as the corresponding xenograft (Fig. 2B), indicating that patients with BUB1B^S tumors were likely to have more invasive disease, resulting in worse overall survival.

We compared the existing molecular GBM tumor classifications (5, 7) to our classification based on BUB1B^S/R. First, we compared the 838 genes associated with BUB1B^S/R to 840 genes that were used to determine the conventional subtypes, including classical, mesenchymal, neural, and proneural subtypes (5). Interestingly, only 35 genes were in common (P > 0.2), representing the uniqueness of our classification of GBM samples based on BUB1B^S/R.

In addition, we compared membership of our BUB1B^S/R subtypes to the previously identified molecular subtypes based on gene expression and CpG island methylation status including G-CIMP, classical, mesenchymal, neural, and proneural subtypes (7). Interestingly, for TCGA data, mesenchymal subtype was enriched for the samples predicted as a BUB1B^R subtype (Fig. 2C; P values from Fisher exact test were 1.6 × 10⁻¹², 0.01, 0.31 for Affymetrix, Agilent, and RNAseq platform, respectively); however, Neuronal subtype was enriched for the samples predicted as a BUB1B^S subtype (Fig. 2C; P values from Fisher exact test were 5.9 × 10⁻⁹, 5.7 × 10⁻⁹, 0.34 for Affymetrix, Agilent, and RNAseq platform, respectively). Other subtypes consisted of similar numbers of BUB1B^R/S-classified samples.

We also compared our classification to the previously identified molecular subtypes for GBM patients of Gravendeel and colleagues data (29) and Joo and colleagues data (30). As there was no reported molecular subtype information associated with the dataset, we classified the GBM samples into four subtypes (i.e., classical, mesenchymal, neural, and proneural) based on the expression signatures previously defined (5). The result was similar to those from TCGA data: only neural subtype (P from Fisher exact test are 7 × 10⁻⁴) and mesenchymal subtype (P from Fisher exact test are 1.0 × 10⁻⁷) were enriched for samples classified as a BUB1B^S subtype and BUB1B^R subtype, respectively, for Gravendeel and colleagues dataset (Fig. 2C). Thus, the previously used molecular subtypes cannot be fully reflected in our classification based on BUB1B^S/R, demonstrating the uniqueness of our classification based on sensitivity to BUB1B inhibition.

To further understand the relationship between the BUB1B^S versus BUB1B^R state in tumor cells, we analyzed the results from genome-wide RNAi screens in BUB1B^S and BUB1B^R GSC isolates, and neural progenitor isolates (33). This screening approach revealed 576 and 122 candidate targets required for expansion of BUB1B^R and BUB1B^S GSCs, respectively (see Materials and Methods for details). We prioritized these candidate therapeutic targets by identifying the biological pathways enriched for screen hits in each subtype, assuming that there are different subnetworks selectively lethal to each subtype. We projected these subtype-specific screen hits onto a GBM network, constructed de novo from GBM samples of TCGA (4) by integrating gene expression and DNA copy number variation data (40, 41). Examination of subnetworks in the GBM network revealed BUB1B^S- and BUB1B^R-specific subnetworks, which were enriched for candidate lethal hits of BUB1B^S or BUB1B^R GSC isolates (Supplementary Fig. S6A and S6B), respectively. The projection of candidate lethal hits of BUB1B^R GSC isolates on GBM network failed to detect significantly large BUB1B^R-specific subnetwork (Supplementary Fig. S6B; see Supplementary Materials and Methods for details); however, we detected two significantly large BUB1B^S-specific subnetwork (Supplementary Fig. S6A and Fig. 3A). The two subnetworks that are specific to BUB1B^S cells (Fig. 3A), containing 11 screen hits (i.e., RACGAP1, POLE, CCNA2, MYBL2, RFWD3, ANXA2, KRR1, PRRP1R12A, NCL, and CPSF6), were significantly enriched for several biological processes related to mitotic cell cycle (P < 1.3 × 10⁻¹¹), microtubule cytoskeleton organization (P < 8.9 × 10⁻⁶), and chromosome segregation (P < 3.4 × 10⁻⁶; Fig. 3B). In addition, the expression levels of these genes were mostly upregulated for BUB1B^S GSCs compared with BUB1B^R GSCs (Supplementary Fig. S7), suggesting that BUB1B^S cells might experience the heightened mitotic stress and render the cell hypersensitive to perturbation of the mitotic machinery. Together, inhibition of these subnetworks will lead to selectively impair the viability of BUB1B^S GSCs. We focus on these selectively lethal subnetworks rather than the individual lethal gene, and further used the subnetworks to identify potential drugs that may be effective to BUB1B^S subtypes in the next section, because the perturbation of the subnetwork tightly connected with the selectively lethal genes will be more efficient to kill cancer stem cells.

We next wished to determine whether we could use our BUB1B^R/S classifier as a predictor of sensitivity to existing drugs, because there are no existing therapeutics strategies for inhibiting BUB1B-GLEBS domain activity. To this end, we endeavored to systematically identify drugs that would selectively impair the viability of BUB1B^S cells using the Connectivity Map (CMAP), encompassing 7,000 genome-wide expression responses to 1,309 different compounds (34). We analyzed the expression changes of genes in the subnetwork specific to BUB1B^S cells (Fig. 3A) after each drug's treatment to identify the drug that inhibits genes within this subnetwork. A negative enrichment score from connectivity map represents the treatment of the drug might inhibit the genes specifically lethal to BUB1B^S cells, resulting in selectively killing of BUB1B inhibition–sensitive GSCs. We identified several potential drugs, including MS-275, etoposide, camptothecin, irinotecan, thioridazine, and azacitidine (Table 2). Interestingly, etoposide targets DNA topoisomerase II whose disruption perturbs centromere organization leading to intense Bub1 on kinetochores (42), and both camptothecin and irinotecan target DNA topoisomerase I, whose treatment in GBM was beneficial (43). Also, thioridazine, an antipsychotic drug, has been identified as an antiglioblastoma and anticancer stem cell agent (44, 45).

To validate our candidate drugs that are specifically cytotoxic to BUB1B^S tumor cells, we analyzed cell line datasets, CCLE and CGP, with gene expression profiles as well as drug sensitivity measurements (35, 36). Application of our classifier to CCLE data yielded 12 and 14 BUB1B^R or BUB1B^S cell lines, respectively, out of 59 glioma cell lines (Supplementary Fig. S8A). For CGP data, 13 and 8 out of 43 glioma cell lines were classified into BUB1B^R or BUB1B^S subtypes, respectively (Supplementary Fig. S8B). Furthermore, we found 27 common cell lines within both CCLE and CGP datasets, and the predictive values of these common cell lines were similar (Fig. 4A; correlation coefficient = 0.69 and P < 5.7 × 10⁻⁵), resulting in consistent classification. In summary, we found that five cell lines are consistently classified into BUB1B^R subtype, and four are classified into BUB1B^S subtype in both datasets.

We further investigated the drug sensitivity of cell lines from both CCLE and CGP dataset. On the basis of the drug sensitivity measurements (i.e., IC₅₀ and activity area) for each glioma cell line and predicted BUB1B^S, we tested whether the two groups showed significantly distinct drug sensitivity. We measured the correlation between drug sensitivity measurements and the predictive value from the elastic net classifier. Encouragingly, 4 of 13 potential drugs (Table 2) were tested in at least one of two datasets (CCLE or CGP) and the sensitivity of three of four drugs was significantly associated with BUB1B^S status. We validated the BUB1B^S subtype showed more sensitive viability to etoposide (P < 0.009) targeting topoisomerase II, MS-275 (P < 0.008) targeting HDAC from CGP dataset, and irinotecan (P < 0.10) targeting topoisomerase I from CCLE dataset (Fig. 4B and C). Even though the sensitivity to camptothecin targeting topoisomerase I was not significantly associated with BUB1B^S from CGP data, sensitivity to topotecan and irinotecan, both targeting topoisomerase I, was different between BUB1B^S and BUB1B^R subtype from CCLE data (Fig. 4B and C), suggesting the drug targeting topoisomerase I as the BUB1B^S subtype-specific drugs. This discrepancy can be explained by two potential reasons: (i) the different drug dose usages of camptothecin between CGP data and CMAP. CGP data uses 0.0003–0.1 μmol/L concentration and CMAP uses 11.4 μmol/L concentration; (ii) tissue-specific effect. Only MCF7 and PC3 that are not glioma cell lines were tested from CMAP; on the other hand, we considered glioma cell lines for CGP data. Interestingly, most glioma cell lines tested showed sensitivity to camptothecin compared with other drugs (Supplementary Fig. S9). Furthermore, we detected PLX4720 targeting BRAF as a potential drug that is more effective to BUB1B^S subtypes from both CCLE and CGP datasets. This suggests the potential importance of RTK–Ras signaling affecting KT function (46).

In addition, we collected drug sensitivities in several glioma cell lines reported in published studies (Fig. 4D). We compared predictive scores and IC₅₀ values from each article. Interestingly, even though CGP data did not confirm the significant sensitivity difference of camptothecin, the previous studies showed that BUB1B^R cell line (i.e., U87) is less sensitive to camptothecin than BUB1B^S cell line (i.e., DBTRG-05; ref. 47) and two cell lines (i.e., U118 and SW1783) that have higher predictive scores (48), consistent with our prediction. Also, Metz and colleagues (49) confirmed that BUB1B^R cell line (i.e., U87) is less sensitive to irinotecan than BUB1B^S cell line (i.e., U251), and the study by Ciesielski and colleagues (48) indicates different sensitivity to etoposide, consistent with our prediction.

Among the predicted drugs that are selectively lethal to BUB1B^S GSCs, we selected two drugs (i.e., etoposide and irinotecan) due to their previous treatment in GBM patients (50, 51). To validate the drugs that show different cytotoxicity to subtypes, we collected nine additional GSC isolates. We measured genome-wide expression levels for these GSCs and applied our classifier to gene expression profiles of these GSC isolates to predict the sensitivity to BUB1B inhibition, resulting in four BUB1B^S and four BUB1B^R samples and one sample ranked in between (Fig. 5A).

These GSC isolates were treated with etoposide and irinotecan to determine their cytotoxicity (Fig. 5B and C). As shown in Fig. 5B (listed in Supplementary Table S4), the exposure to etoposide reduced cell growth of 9217B compared with that of 10647B: the IC₅₀s of etoposide for 9217B and 10647B were 144 and 4.857 mmol/L, respectively. We detected the same pattern for irinotecan as etoposide. The BUB1B inhibition–sensitive GSC–9217B and the GSC ranked between sensitive and resistant 9260B showed more reduced cell growth compared with BUB1B inhibition–resistant GSC–10647B (Fig. 5C). These results indicate that etoposide and irinotecan have selective cytotoxicity to BUB1B^S GSCs, consistent with our predictions.

One of the promises of precision medicine is the improved ability to predict which treatments will work best for specific patients. Our results suggested that GSCs can be grouped into BUB1B^S or BUB1B^R subtypes. Our BUB1B^S/R-associated gene signatures are not associated with previously identified gene lists of GBM classification, suggesting the uniqueness of our classification. Our signature was highly enriched for mitotic cell cycle, supporting the association with short IKDs and lethal KT-MT instability. We extracted GBM tumor cell profiles deconvoluted from bulk GBM tissue expression levels to identify new prognostic subtypes of GBM that may have been obscured previously by confounding normal and GBM tumor cells. Our BUB1B^S/R subtype based on sensitivity to BUB1B inhibition was associated with overall survival, and BUB1B^S subtype had poor survival rate presumably because patients with these tumors tended to have more invasive disease. We identified lethal genes to GSCs of each subtype, and further suggested potential treatments that would be more effective to BUB1B^S subtype by connectivity map database accompanied with siRNA screens and GBM causal networks.

The gene signature associated with the BUB1B^S/R subtype was able to classify GSCs and GBM tumors into two subtypes with clinical relevance. Previous studies showed that patient cohorts treated with etoposide had significant improvement in median overall survival (50), and the treatment of bevacizumab with irinotecan was effective in recurrent GBM patients (43). We predicted and validated the effective responses to etoposide and irinotecan of BUB1B^S GSCs. We expect that classification of BUB1B^S/R subtype would improve the effectiveness of these drugs to GBM patients. Whether the patients with BUB1B^S tumors responding better to these treatments needs to be studied further. Also, in the future, high-throughput drug screening experiments, we may better focus on screening compounds for only BUB1B^S cells instead of screening compounds for all GBM tumor cells. However, we identified several drugs, such as MS-275, which show differential sensitivity depending on BUB1B^R/S classification in in vitro cell lines. However, drugs including MS-275 show marginal to moderate antiglioma effects in clinical trials (52, 53). Because of the blood–brain barrier (BBB) that protects tumors from exposure to many active drugs, drugs that are very effective in solid tumors frequently fail to demonstrate activity in brain tumors (54). Alternative strategies targeting GBM tumors are needed. Recently, FDA approved electric fields therapy referred to as tumor-treating fields (TTFields) as a treatment of recurrent glioblastoma. TTFields target rapidly dividing tumor cells, resulting in aberrant mitotic exit of tumor cells and disruption of cells during mitosis (55). The molecular mechanism targeted by TTFields is similar to the one of classifying BUB1B^R/S cells. It will be useful to investigate the molecular targets of TTFields by genome-wide transcriptomes and compare them with our genes used to classify BUB1B^R/S as well as to investigate whether BUB1B^S GBM patients have better prognosis with TTFields treatment.

Even though targeting GSCs is a promising therapeutic option to treat GBM, most GBM tumor expression profiles measured in bulk tissues resulted from a mixture of tumor as well as normal cells. The classifier based on GSC expression levels was not suitable to apply to tissue expression levels directly because normal cells were BUB1B^R (17). Our proposed procedure was not only applicable to classify GSC isolates, but also to classify GBM tissue into BUB1B^S and BUB1B^R subtypes, based on the extracted GBM tumor cell profiles from the tissue expression decomposition. Our analysis of tissue and GSC expression indicated that the proportion of tumor cells was variable across samples.

The heterogeneous nature of the tumors from the same patient needs to be taken into account to develop therapies. One challenge will be how to treat the patients with a mixture of tumor cells from distinct BUB1B^S/R subtypes. Indeed, our results indicate that some of tumors from the same patients have potentially distinct BUB1B^S/R subtypes (e.g., both 9260B and 10647A were derived from two different tumors occurred at different time in the same patient). Also, when considering tumor tissue expression of GBM patients, we detected nonclearly classified samples whose predictive scores were close to zero, indicating the potential subpopulations with distinct BUB1B^S/R subtypes in the tumor tissues.

To assess tumor cell heterogeneity with regard to BUB1B^S/R subtypes, we applied our approach to a publicly available single-cell sequencing dataset (2), which consists of single-cell transcriptomic profiles of stem-like tumor-propagating glioblastoma cells (GBM6 and GBM8) and single cells isolated from five tumors as well as population transcriptomic profiles of five tumor tissues and stem-like cells from three tumors as controls. For each transcriptomic profile of either single-cell or population controls, we applied our BUB1B^R/S classifier and calculated predicted values from elastic net. First, we analyzed single-cell sequencing data for GSCs (GBM6 and GBM8), and expected homogeneous transcriptional profiles as well predicted BUB1B^R/S scores. However, we found that the predicted BUB1B^R/S score was associated with the quality of single-cell data with the predicted score anticorrelated with the number of nonexpressed genes in each single-cell sequencing data (Supplementary Fig. S10). If the cells of low-sequencing quality data were removed, the predicted scores for GSCs were indeed homogeneous as expected (Supplementary Fig. S10). The result suggests that we need to filter out low-quality single-cell sequencing data that contains a large number of nonexpressed genes in further data analyses. Second, we examined intratumor heterogeneity of single cells from five individual tumor tissue samples. Again, we found that the predicted BUB1B^R/S score was associated with the quality of single-cell data with the predicted score anticorrelated with the number of nonexpressed genes in each single-cell sequencing data (left panels in Supplementary Fig. S11A–S11E), which was consistent with the trend observed in GSCs (Supplementary Fig. S10). After removing cells of low-quality sequencing data, there was still transcriptional heterogeneity within cells from the same individual tumor tissue sample (left panels in Supplementary Fig. S11A–S11E). Third, we used the dataset to further test our deconvolution method for predicting transcriptional feature of extracted GBM tumor cells from bulk tumor samples. The predicted BUB1B^R/S values for GSCs were more similar to the values for the extracted GBM tumor cells than to the values for bulk tumor samples (right panels in Supplementary Figs. S11A–S11E and S12), and the predicted scores of the extracted GBM tumor cells were close to the average of scores of single cells. For example, the tumor MGH31 was predicted as BUB1B^R based on bulk tumor tissue transcriptional profile, whereas it was predicted as BUB1B^S based on both GSC and extracted GBM tumor cell transcriptional profiles. The results of single-cell data indicate intratumoral heterogeneity of tumors. Thus, for future personalizing treatments, it is necessary to explore intratumoral heterogeneity with respect to the BUB1B^S/R subtypes.

In conclusion, our results suggest a new way to classify GBM tumors into two molecular subtypes, which have different prognosis and response to potential drug treatments. Further understanding of these subtypes may provide critical information to develop precision medicine for the deadliest cancer, GBM.

No potential conflicts of interest were disclosed.

Conception and design: E. Lee, R.L. Yong, P. Paddison, J. Zhu

Development of methodology: E. Lee, M. Pain, J.G. DeLuca, R.L. Yong, J. Zhu

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): H. Wang, J.A. Herman, C.M. Toledo, R.L. Yong, P. Paddison

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): E. Lee, J.A. Herman, R.L. Yong, P. Paddison, J. Zhu

Writing, review, and/or revision of the manuscript: E. Lee, M. Pain, R.L. Yong, P. Paddison, J. Zhu

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): E. Lee, J.G. DeLuca, J. Zhu

Study supervision: R.L. Yong, J. Zhu

We thank members of Zhu laboratory for discussions.

This work was partially supported by the NIH (R21-CA170722, R01-AG046170, and U01-HG008451 to J. Zhu) and (R21-CA170722 and R01 CA190957 to P. Paddison).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.