Drug and single-cell gene expression integration identifies sensitive and resistant glioblastoma cell populations

Domain: chemistry
Published: 2026-01-07
Editorial outcome: major_revision

Abstract
Glioblastoma (GBM) remains the most common and lethal adult malignant primary brain cancer with few treatment options. A significant issue hindering GBM therapeutic development is intratumor heterogeneity and plasticity. GBM tumors contain neoplastic cells within a fluid spectrum of diverse transcriptional states. Identifying effective therapeutics requires a platform that predicts the differential sensitivity and resistance of these states to various treatments. Here, we develop scFOCAL (Single-Cell Framework for-Omics Connectivity and Analysis via L1000), to quantify the cellular drug sensitivity and resistance landscape. Using single-cell RNA sequencing of newly diagnosed and recurrent GBM tumors, we identify compounds from the LINCS L1000 database with transcriptional response signatures selectively discordant with distinct GBM cell states, and leverage this capability to predict combination synergy. We validate the significance of these findings in vitro, ex vivo, and in vivo, and identify a combination of an OLIG2 inhibitor and Depatux-M for the treatment of GBM. Our studies suggest that scFOCAL identifies cell states that are sensitive and resistant to targeted therapies in GBM using a measure of cell and drug connectivity, which can be applied to identify new synergistic combinations.

Introduction
Glioblastoma (GBM) remains the most common malignant adult brain cancer, with a median overall survival of only 15 months1–3. No new targeted therapies have been approved for GBM since the introduction of the alkylating agent temozolomide (TMZ) in 20051–3. Thus, novel therapeutic options are direly needed for patients with GBM.

Intratumor heterogeneity poses a significant barrier to identifying effective targeted therapies for GBM4,5. Single GBM tumors contain diverse populations of cells with considerable genomic, transcriptomic, and proteomic differences6–9. GBM single-cell characterization has revealed a plastic transcriptional spectrum reminiscent of canonical neurodevelopmental cell types6. These transcriptional states encompass astrocyte-like (AC-like), neural-progenitor-like (NPC-like), oligodendrocyte-progenitor-like (OPC-like), and mesenchymal-like (MES-like) cells. The importance of their relative abundance is highlighted by the identification of transcriptional subtype (Classical, Proneural, Mesenchymal) plasticity and switching following standard-of-care treatment10. Therefore, it is necessary to develop novel means of targeting these dynamic tumor cell populations. However, there is no current method to predict which compounds target the different neoplastic cell states within GBM tumors. We and others have shown that cancer gene expression-based disease signature reversal is predictive of compound efficacy11–13. Using the NIH Library for Integrated Network-based Cellular Signatures (LINCS) L1000 assay dataset, we previously developed transcriptional consensus signatures (TCSs) for each small molecule, consisting of genes that are consistently up- or down-regulated irrespective of cell type. We showed that effective combination treatments can be identified through the integration of compound TCSs with TCGA-derived GBM disease signatures and that compound efficacy in a GBM subtype-specific manner can be predicted12–14. However, these analyses do not address the GBM intratumor heterogeneity that is evident after performing single-cell RNA sequencing (scRNAseq) to predict mechanisms of drug sensitivity or resistance.

Here, we describe a framework for the integration of scRNAseq data and L1000 TCSs to facilitate in silico perturbation analyses informing combination therapy design, scFOCAL (Single-Cell Framework for-Omics Connectivity and Analysis via L1000). Using scFOCAL, we demonstrate that the integration of GBM single-cell disease signatures with LINCS L1000-derived compound response signatures can identify drugs most likely to target transcriptionally distinct cell populations in individual GBM tumors. Furthermore, we demonstrate that this framework predicts the sensitive and resistant cell populations within GBM tumors in an orthotopic xenograft model. Using patient single-cell RNA sequencing data and a TCS for the aurora kinase inhibitor alisertib, we predict and confirm both the MES-like transcriptional identity of a distinct alisertib-resistant cell population and the depletion of transcriptionally NPC-like cells in vivo15. Moreover, we use publicly available datasets to demonstrate the utility of this approach to predict treatment-induced transcriptional response in discrete GBM cell populations and to identify synergistic small molecule combinations. Building on this, we leverage our framework to identify a combination of the OLIG2 inhibitor, CT-179, with an anti-EGFR antibody-drug conjugate, Depatuxizumab Mafodotin (Depatux-M, ABT-414), which synergizes to increase survival in an orthotopic xenograft model of GBM. Importantly, we introduce scFOCAL as a publicly available R package and Shiny application16. scFOCAL uses the methods defined here to generate and analyze cell-type-specific disease signatures, perform in silico drug connectivity analyses, characterize cell population sensitivity or resistance to specific treatments, and prioritize synergistic combinations.

Results
To identify the different cell types within GBM tumors, we performed single-cell RNA-sequencing of six tumors resected from three newly diagnosed and three recurrent GBM patients. We combine our sequencing results with those from published studies that use the same 10X Genomics platform to study GBM at the single cell level, thereby increasing sample size for our downstream analyses (Fig. 1)17. Our dataset (here to referred to as Suter et al.) and the Johnson et al. dataset are integrated and assessed using CytoTRACE, revealing a large tumor cell population marked by transcriptional diversity and dedifferentiation (Fig. 1d)18. This population contained tumor cells spanning the neurodevelopmentally-rooted GBM cell transcriptional states identified by Neftel et al. (Fig. 1e, f)6. In addition to neoplastic cells, both datasets contain populations of myeloid cells, T-cells, fibroblasts, endothelial cells, and oligodendrocytes. These non-neoplastic cell populations that comprise the tumor microenvironment (TME) are identified by the expression of specific marker genes, such as PTPRC for myeloid cells and CD3E for T-cells (Supplementary Fig. 1), and contain cells captured across both datasets (Fig. 1g). We then designate the transcriptional profiles of these TME cells as our normal cell reference to be able to compare the distinct GBM tumor cell states represented within our dataset. Utilizing this comparison, we were able to derive a pseudo-bulk disease signature for each of the cell states within GBM tumors (Supplementary Fig. 1h). The disease signatures unique to each transcriptional state reflected previously reported expression markers, such as CDK4 and SOX4 overexpression in NPC- and OPC-like cells, EGFR overexpression in AC-like cells, and CD44 overexpression in MES-like cells (Supplementary Fig. 1h).

Integration of L1000 TCSs with expression data of single GBM cells clusters, small molecules by mechanism of action
While distinct GBM single-cell transcriptional states have unique disease signatures relative to TME cell controls, this method relies on prior stratification and knowledge of tumor cell populations. We therefore seek a way to integrate L1000 TCSs with the expression of individual cells. We use Spearman’s ρ between a single cell’s normalized and scaled expression and TCS signatures to score compounds for their connectivity with each cell in the dataset, as done previously using bulk gene expression data (Fig. 2a)19,20. Through correlation analysis of each compound’s calculated connectivity across the single-cell landscape, we found that many compounds cluster by their mechanism of action based on their tumor-cell connectivities (Fig. 2b, c, and Supplementary Fig. 2a, b). Compound classes such as MEK, PI3K, and mTOR, HSP, HDAC, and BET inhibitors cluster individually based on their correlations with GBM tumor single-cell expression data (Pearson’s ρ > 0.7), suggesting that biologically relevant information is retained in this perspective (Fig. 2c). Interestingly, when looking only at 64 FDA-approved oncology drugs identified using annotations from the Drug Repurposing Hub21, distinct clusters of small molecules can be identified based on their calculated connectivity with all cells within the GBM scRNAseq atlas (Supplementary Fig. 2a, b).

The integration of L1000 TCSs of FDA-approved oncology drugs with the expression data of single GBM cells identifies distinct drug connectivity states within single tumors
To investigate whether drug connectivity patterns correlate with previously published GBM cell transcriptional states, we perform a clustering analysis of single cells from individual tumors based on their connectivity to 64 different FDA-approved drugs (Fig. 2d–g). Within individual tumors, clear patterns of distinct drug connectivity states are revealed, with some individuals containing at least 7 different states, where the number of unique states was determined using the ideal number of clusters k based on within-cluster sum of squares (Supplementary Fig. 3a). These patient-level drug connectivity clusters show patterns in their makeup by Neftel states, developmental potential quantified using CytoTRACE, and cell cycle phase (Fig. 2f, g, column annotations, and Supplementary Fig. 3a)18. Hierarchical clustering of all resulting patient clusters on their associated Neftel et al. state proportions, cell cycle proportions, and mean CytoTRACE score reveals an organization into distinct, high-level metaclusters—namely a proliferative cluster, an NPC-like cluster, an NPC/OPC-like cluster, an AC-like cluster, and a MES-like cluster (Supplementary Fig. 3b).

Integration of L1000 TCSs with single-cell-derived disease signatures predicts differential sensitivity to small molecules
We seek to determine whether the identified pseudo-bulk disease signatures for AC-like, MES-like, OPC-like, and NPC-like states within GBM tumors can be targeted with distinct small molecules or drugs. To accomplish this, we use a previously reported method of disease signature reversal scoring13. In brief, a disease-signature-specific discordance score was calculated for all compounds within the L1000 library as the number of genes perturbed by the compound TCS in the opposite direction of the disease signature expression, over the number of genes perturbed by compound TCS in the same direction as in the disease signature (Fig. 2b row annotations). Although a molecule can strongly reverse a bulk expression-based disease signature, reversal may be unique to a subpopulation of tumor cells. Many small molecules are predicted to reverse the disease signatures of MES-like, AC-like, NPC-like, or OPC-like cells, but rarely are they predicted to strongly reverse the disease signatures of all transcriptional states (Fig. 2b row annotations, and Supplementary Data 4). Utilizing the same FDA-approved oncology drugs used in tumor cell clustering, we find that many are differentially connected between the four transcriptional states (Fig. 2h).

In silico perturbation analysis using scFOCAL predicts the intratumor response to alisertib treatment in vivo
Utilizing our calculated cell-drug connectivity data, we find that GBM cell transcriptional states show differing mean connectivities to different clinical oncology drugs (Fig. 2h). Here, we find 54 out of 64 compounds to have significantly different connectivities between cell states, determined by generalized linear model analysis using limma (BH-adjusted p-value < 0.05, Fig. 2h), including the aurora kinase inhibitor alisertib. Prior studies reveal that alisertib increases survival in an orthotopic xenograft mouse model of GBM. However, tumors eventually grow back, suggesting an acquired mechanism of resistance13,22–24. Aside from their role in cell cycle progression, the aurora kinases also play an important role in neurodevelopment and have been extensively studied as targets not only in GBM but many other cancers15,24–50. For this reason, we seek to identify whether a specific cell state confers resistance to aurora kinase inhibition in GBM. Overall, the alisertib TCS consists of 139 genes, 54 that are upregulated relative to vehicle control, and 85 that are downregulated. Of these 139 genes, 49 were differentially expressed among AC-, MES-, NPC-, and OPC-like cells. Utilizing differential expression between cells belonging to each transcriptional state, we find that the alisertib TCS predicts that alisertib treatment promotes the expression of MES-like specific genes such as HMOX1 or SQSTM1, while reducing the expression of NPC-like specific genes such as CDK4 or EZH2 (Supplementary Fig. 4b). Additionally, the alisertib TCS discordance ratio, as utilized in SynergySeq, is highest against the NPC-like pseudo-bulk disease signature (Supplementary Fig. 4c). Using the calculated discordances of the alisertib TCS with the individual cells within the dataset, we performed an in silico drug connectivity analysis, where cells are predicted to be sensitive or resistant based on their expression concordance with the alisertib TCS (Fig. 3a, b). Predicted alisertib sensitive and resistant cell populations span source dataset, new diagnoses and recurrences, individual patient tumors, and cell cycle phases (Supplementary Fig. 5a–d). Labeling cells by their predominant meta-module enrichment, the predicted resistant cell populations across individual patients consistently show an increase in proportion of cells with MES-like states and depletion of NPC-like cells (Fig. 3c). The proportion shifts between NPC-like and MES-like states show the strongest differences, suggesting an NPC-like to MES-like predicted shift (Supplementary Fig. 4d). Differential expression testing between predicted alisertib-resistant and sensitive populations reveals a predicted resistance signature that was most strongly enriched in MES-like cells (Supplementary Fig. 5f). Further, expression of the NPC-like marker CDK4 is significantly decreased in the predicted resistant cell population (Supplementary Fig. 5e).

To test the predictions from our in silico drug connectivity analysis using the patient tumor dataset, we use scRNAseq of orthotopic xenografts treated with alisertib or vehicle (n = 3 per group) to assess alisertib-induced shifts in transcriptional state composition (Fig. 3d). As most tumor cells in our patient dataset were AC-like (Fig. 1e–g), we utilized GBM22 PDX cells, which are molecularly classical, for our in vivo studies. Single cell analysis of GBM22 tumors isolated from mice reveals the presence of tumor cells within each of the transcriptional GBM cell states (AC-, MES-, NPC-, and OPC-like) (Fig. 3g, h, and Supplementary Fig. 6a–c). Importantly, in vivo alisertib treatment induces a depletion of NPC-like expression enrichment and an increase of MES-like enrichment (Fig. 3i, j), which aligns with our predictions (Fig. 3c, Supplementary Fig. 4d). Likewise, relative to more modest changes in proportion for AC-like and OPC-like cells, the proportion of GBM22 cells in a predominantly NPC-like state is reduced, while the proportion of GBM22 cells in a predominantly MES-like state is increased (Fig. 3i, j).

In silico drug connectivity analysis predicts the intratumor compound discordance shifts observed in vivo
To assess whether resistant cell states predicted through in silico drug connectivity analysis of the patient tumor dataset represented those observed to persist through alisertib treatment in vivo, we directly compared the differences in global L1000 drug discordance between the predicted alisertib-sensitive and resistant cells in the patient data to those observed through comparing alisertib and DMSO-treated xenografts (Fig. 3k, and Supplementary Fig. 7a–h). We find that predictions of discordance shift from in silico drug connectivity analysis with alisertib in patient cells are predictive for certain classes of small molecules, including nucleoside reverse transcriptase inhibitors (Spearman’s ρ = 0.89, p = 0.033), tubulin polymerization inhibitors (Spearman’s ρ = 0.82, p = 0.0068), topoisomerase inhibitors (Spearman’s ρ = 0.82, p = 0.00011), serotonin receptor agonists (Spearman’s ρ = 0.82, p = 0.0068), MEK inhibitors (Spearman’s ρ = 0.77, p = 0.021), HCV inhibitors (Spearman’s ρ = 0.77, p = 0.021), HMGCR inhibitors (Spearman’s ρ = 0.83, p = 0.058), and mTOR and PI3K inhibitors (Spearman’s ρ = 0.94, p = 0.017) (Fig. 3k, and Supplementary Fig. 7a–h).

Connectivity analysis with an L1000-derived panobinostat TCS separates vehicle-treated from panobinostat-treated cells in both neoplastic and myeloid cell populations
To assess the application of our prediction model in ex vivo settings, we utilized a publicly available scRNAseq dataset of GBM acute slice cultures treated with the pan-HDAC inhibitor panobinostat (Zhao et al., 2019) (Supplementary Fig. 8a, b)51. A panobinostat TCS is derived from the 2020 L1000 data release (Supplementary Fig. 8c) and used to score panobinostat connectivity of DMSO or panobinostat-treated neoplastic and myeloid cells. In both populations, panobinostat-treated cell populations showed significantly higher panobinostat connectivity values (Supplementary Fig. 8d, e). These findings are in agreement with those reported within the original manuscript in which this dataset was reported, wherein panobinostat elicits alterations to the transcriptional states of both tumor and myeloid cell populations51.

Connectivity analysis predicts drug combinations that synergize in vitro and in vivo
Having demonstrated the predictive power of our connectivity analysis framework to identify tumor cells resistant or sensitive to different perturbations, we seek to leverage our framework to identify synergistic combinations for glioblastoma (Fig. 4a). First, we utilized a recently published synergy screen performed across 24 independent glioma stem cell lines in spheroid culture to assess the ability of scFOCAL to predict relative synergy of small molecule combinations (Supplementary Fig. 9)52. Using erlotinib, lapatinib, pazopanib and sunitinib as reference compounds (Supplementary Fig. 9a), we ranked 13 small molecule combinations with these inhibitors using an scFOCAL based combination score (Supplementary Fig. 9b, top heatmap annotation). This score is representative of the specificity of a partner molecule in targeting the cell population resistant to the reference molecule and thus reflects the breadth of coverage a combination would have against the GBM cell transcriptional states globally represented within our patient GBM scRNAseq atlas. We find that this scFOCAL-based combination score strongly correlates with the mean observed BLISS synergy across cell lines of the combinations present within this screen (Spearman ρ = 0.7, p = 0.01), indicative of the potential that a representative GBM single-cell atlas can be used as a representative model for identifying more effective combinations (Supplementary Fig. 9c). Next, we focus on the OLIG2 inhibitor CT-179 due to its presumed ability to target OPC-like cells53–57. Using bulk RNA sequencing of molecularly classical GBM8 cells treated with 200 nM CT-179 or vehicle, we generate a dose-dependent CT-179 transcriptional signature comprising 841 differentially expressed genes for analysis of cell-CT-179 connectivity in our patient tumor atlas (Fig. 4b, and Supplementary Fig. 10a–c). By gene ontology analysis, the CT-179 response signature is enriched for processes including microtubule and tubulin binding, motor activity, and processes related to mitotic cell division (Supplementary Fig. 10b). Using the Spearman correlation of the directional CT-179 response signature with each tumor cell’s scaled expression, cells are binned into predictive sensitive or resistant populations, which spanned all patients within our dataset (Fig. 4c). As expected, the CT-179 predicted sensitive population is primarily OPC-like (Fig. 4d), as evidenced by low CT-179 connectivity in OPC-like cells and high CT-179 connectivity in AC-like cells (Fig. 4e). Further, the relative proportions of cells in an AC-like state is predicted to increase with CT-179 treatment, relative to the predicted depletion of cells in both NPC- and OPC-like states (Supplementary Fig. 10d). Using our established workflow for predicting the differential connectivities of L1000 small molecules between predicted sensitive and resistant cell populations, we identify 721 L1000 small molecules that exhibit a significantly reduced connectivity in CT-179 resistant cells compared to CT-179 sensitive cells (logFC < 0, adj. p < 0.05) (Fig. 5a). As above, we also identified 706 discordant L1000 small molecules by their mean connectivity to the CT-179 resistant population (mean resistant cell connectivity, mRCC < 0) (Fig. 5b). 411 L1000 small molecules exhibited both decreased connectivity and overall discordance with the CT-179 resistant population. By weighting the connectivity differentials of these molecules with their mean discordance to the CT-179 resistant cell population (mean Resistant Cell Connectivity, RCC), we established a combination index score for each molecule within this intersect (Supplementary Data 6). Top hits by this metric included the EGFR inhibitor varlitinib, as well as docetaxel and indibulin, which elicit anti-neoplastic effects through alteration of tubulin and microtubule assembly and dynamics (Fig. 5c)58–61. Indeed, EGFR and OLIG2 are expressed by distinct GBM cells in AC-like and OPC- and NPC-like states, respectively (Supplementary Fig. 10d). These prioritizations lead us towards the use of the anti-EGFR monomethyl auristatin F (MMAF) antibody-drug conjugate Depatux-M (ABT-414), which would target EGFR-expressing cells with MMAF, a tubulin poison, in combination with CT-179. Mice bearing GBM6 orthotopic xenograft tumors were treated with vehicle, CT-179 alone, Depatux-M alone, or CT-179 and Depatux-M in combination. Treatment of mice bearing GBM6 patient-derived xenografts with the CT-179/Depatux-M combination increases survival over vehicle treatment or either monotherapy alone (Fig. 5e). In corroboration of the increased survival, bioluminescence imaging (BLI) reveals a significant reduction in tumor size following combination treatment with CT-179 and Depatux-M relative to CT-179 monotherapy alone (Fig. 5d). Collectively, these studies suggest that targeting OPC-like and AC-like cells simultaneously may be advantageous in vivo, and that maximizing drug discordance across the tumor cell landscape identifies more effective combinations.

Release of the scFOCAL framework as an R package and a web application
We package the methodology defined and utilized here into an analytical framework termed scFOCAL (Single-Cell Framework for-Omics Connectivity and Analysis via L1000), which we have made available as a Shiny web application and as an R package. As we have demonstrated, scFOCAL integrates bulk drug-response TCSs derived from the LINCS L1000 dataset with multi-subject scRNAseq data and facilitates the analysis of drug and cell discordance from multiple perspectives. First, scFOCAL provides a simple approach to single-cell-derived disease signature generation and permits the analysis of disease signature heterogeneity across tumors relative to non-transformed cell types. Next, scFOCAL scores small molecules for their reversal potential of the generated disease signatures. Compounds predicted to reverse single-cell-derived disease signatures can then be utilized in scFOCAL’s in silico drug connectivity approach. Here, as demonstrated with alisertib and panobinostat, L1000 TCSs are scored for discordance against the scaled expression of single-cells in a dataset. A custom transcriptional response signature can also be used, as we show in the case of CT-179. Users can then select a sensitivity cut-off (i.e., pseudo-dose), and the single-cell model will be split into predicted sensitive and resistant populations. Characteristics and visualizations, including proportions of annotated identities and predicted drug sensitivities, can then be explored. Additionally, users can identify predicted synergistic compounds through calculation of an scFOCAL combination index score for each L1000 compound with their drug of interest. All data generated by users of scFOCAL are available to be downloaded for ease of access. Lastly, the scFOCAL shiny application integrates seamlessly with the R package OrthologAL and the OrthologAL shiny app, which can prepare Seurat objects for dual-species xenografts and non-human expression data for analysis in scFOCAL62. Single-cell connectivity scores for all compounds against the patient GBM scRNAseq atlas presented in this study are downloadable using the scFOCAL Shiny web application, which can be found at https://robert-k-suter.shinyapps.io/scFOCAL alongside instructions for R package download at https://github.com/AyadLab/scFOCAL (Fig. 6).

Methods
Collection of newly diagnosed and recurrent GBM tumors
Six (6) patients with GBM permitted tissue collection before tumor resection under approved IRB study 20170887 (MODCR00002196). Of the 6 enrolled patients, 3 were newly diagnosed with GBM and 3 had tumor recurrence after previously receiving standard of care treatment with temozolomide chemotherapy and radiotherapy. The 3 newly diagnosed patients had not received any treatment before resection. After surgical resection, University of Miami Hospital pathologists separated the resected tissue, which was then transferred from the operating suite to the laboratory on ice. Tissue from the enhancing edge was then used for dissociation and single-cell RNA sequencing (see below). Patient samples were obtained with written, informed consent.

GBM cell xenograft models
Human patient-derived GBM cells (GBM6, GBM8, GBM22) were obtained from the Mayo Clinic Brain Tumor Patient-Derived Xenograft National Resource. GBM22 cells were transduced to express GFP (pSMAL-CellTag-V1), while GBM6 cells were transduced to express eGFP-FLUC2 as previously described76,81–83. Briefly, lentiviral particles were packaged in HEK293T cells, and transduction to primary GBM22 and GBM6 cells was performed in the presence of 5 μg/ml polybrene. GBM8 was derived from a newly diagnosed tumor in a female patient and is molecularly IDHwt, Classical, EGFR amplified, and MGMT methylated. GBM22 was derived from a newly diagnosed tumor in a male patient and is molecularly IDHwt, Classical, MGMT methylated. GBM6 was derived from a newly diagnosed tumor in a male patient and is molecularly IDHwt, with EGFRVIII mutation, MGMT unmethylated.

Mouse models
All research was conducted in compliance with ethical regulations. Athymic nude mice (Crl:NU(NCr)-Foxn1nu) from Charles River Laboratories (stock 553) were utilized for in vivo experiments. Orthotopic implantations took place when mice were 6–8 weeks old. All procedures were approved by the University of Miami IACUC, protocol number 18-014, or by the Mayo Clinic IACUC, protocol number A4595-19. All animal procedures were performed per guidelines and regulations, including ARRIVE guidelines13,76,84,85. Mice were housed in the respective Department of Veterinary Resources on corncob bedding in individual HEPA ventilated cages (Innocage® IVC, Innovive USA) on a 12-h light-dark cycle at 68–79°F (20–26.1°C) and 30–70% humidity. Animals were fed water (reverse osmosis, 2ppm Cl2) and an irradiated rodent diet (LabDiet 5053 PicoLab Rodent Diet 20) ad libitum. Animals were monitored for signs of lethargy or head-tilting. For IVIS studies, maximum tumor volume did not exceed 60 mm3, in accordance with the approved protocol. Experimental groups for these experiments were assigned using a random number generator.

Preparation of newly diagnosed and recurrent GBM tumors for scRNAseq
We performed single-cell RNA sequencing on the 6 patient GBM tissue samples collected at the time of surgical resection (n = 3 newly diagnosed, n = 3 recurrent, Fig. 1) After resection, University of Miami hospital pathologists separated the GBM tissue, which was transferred from the operating suite to the laboratory on ice. Tumor samples were dissociated using the Worthington Biochemical Papain Dissociation System86. Tumor tissue was finely minced with a sterile surgical blade and incubated in sterile Earl’s Balanced Salt Solution (EBSS) containing 20 units of papain per mL in 1 mM L-cysteine with 0.5 mM EDTA and DNase for 45 min at 37°C while shaking. Following papain incubation, the tissue was gently triturated with a 5 mL serological pipette until it was adequately dissociated. The resulting suspensions were pelleted and resuspended in EBSS with ovomucoid protease inhibitor, bovine serum albumin (BSA), and DNase. This suspension was layered on top of ovomucoid protease inhibitor solution with BSA to form a gradient, which was centrifuged at 70×g for 6 min. The supernatant, which contained cellular debris, was removed. The pellets were resuspended, and the gradient centrifugation was repeated. The resulting pellet was then washed 2 times with ice cold PBS with 0.1% BSA, passed through a 0.45 μM cell strainer, and resuspended in ice-cold PBS with 0.1% BSA. This resulting suspension was then checked for cell viability by staining an aliquot with acridine orange and propidium iodide (AO/PI) and counted on a Nexcelom K2 Cellometer. All samples had cell viability greater than 90% prior to loading of the single-cell chips. Single cells were captured using a 10X Genomics Chromium Controller, and patient single-cell gene expression libraries were prepared using 10X Genomics 5′ Gene Expression chemistry. All libraries were sequenced on an Illumina NextSeq 500.

Pre-processing and quality control of patient-derived scRNAseq data
Raw scRNAseq dataset was aligned and demultiplexed using Cell-Ranger (10X Genomics), and the resulting filtered unique molecular identifier (UMI) matrices were used for downstream analysis in the R package Seurat (v4)87–91. Within Seurat, poor quality cells were further filtered by removing outliers based on the number of detected UMIs and unique transcripts, as well as percent expression of mitochondrial transcripts. Pass-filter cells were those with greater than 200 and less than 4500 detected features, and with less than 12.5% mitochondrial expression. Data from 42,839 pass-filter single-cell transcriptomes across 33,469 features were normalized and variance stabilized using a regularized negative binomial regression implemented within the R package scTransform, and individual patient datasets were integrated using an anchor-based technique utilizing the Pearson residuals87–91.

Harmonization and integration of in-house and external scRNAseq datasets
Previously published single-cell RNA sequencing data from Johnson et al. (2021) were downloaded from Synapse (https://synapse.org/singlecellglioma)17. The five IDHwt GBM samples out of the eleven total glioma samples available in the dataset were used17. To integrate the Johnson et al. dataset with our internal dataset, each was first split into individual datasets based on the patient source. Each individual patient’s expression matrix was then independently prepared for integration using scTransform, and integrated together using the anchor-based technique in Seurat (Fig. 1b). 28,130 single-cell transcriptomes from Johnson et al. were integrated with our dataset for a total of 70,969 single cells. Integrated expression values were used for initial shared-nearest-neighbor (SNN) clustering, dimensionality reduction, and cell type identification (Fig. 1a). Downstream analysis, including differential expression and enrichment scoring, was performed using log-normalized counts.

Generation of GBM cell transcriptional state-specific disease signatures
SNN clustering on the patient-derived, variance-stabilized, integrated scRNAseq data identified distinct clusters of neoplastic and non-neoplastic cells. To confirm the identification of these neoplastic cells, we assessed transcriptional diversity with the R package CytoTRACE using raw counts as input. CytoTRACE scores cells based on their developmental potential, reflecting the plasticity of GBM tumor cells relative to the TME (Supplementary Fig. 1)18. Non-neoplastic cell identities were assigned by expression of canonical cell type markers identified through differential expression with MAST (Model-based Analysis of Single-cell Transcriptomics) (Supplementary Fig. 1, and Supplementary Data 1)92. The neoplastic cells were further subdivided and assigned GBM cell state identities using GBM cell transcriptional state signatures previously reported in Neftel et al. (2019) and the R package singscore (Supplementary Fig. 1)6,93. Briefly, the raw count matrices across all samples were log-normalized and scaled, regressing out S-Phase and G2/M-Phase expression module scores as calculated in Kowalczyk & Tirosh et al. (2015). These scaled data were then used to generate gene rank lists for each individual cell. Using the R package Singscore, gene set enrichment analysis was performed on each cell independently for the six Neftel et al. GBM cell state signatures (meta-modules). The means of NPC1- and NPC2-like, and MES1- and MES2-like state enrichment scores were calculated and represented as NPC-like and MES-like enrichment, respectively. Each cell’s assigned transcriptional identity reflects its predominant, most-enriched GBM cell state expression program (AC, MES, NPC, OPC). Thus, the identities assigned within this manuscript indicate the predominant transcriptional program of that cell. Hierarchy plots were generated using the R package scrabble6. Differential expression testing using MAST was performed between each of the state-specific clusters and all non-neoplastic cells to generate disease signatures (Fig. 2a, and Supplementary Data 2)92. MAST was used to account for the technical pitfall of stochastic dropout and the resulting bimodal expression distributions inherent in scRNAseq. Individual patient signatures were generated by comparing each individual’s tumor cells to all non-neoplastic cells from all patients within the dataset (Supplementary Data 3). Heatmaps were generated using the R package pheatmap94.

In silico drug connectivity analysis with the aurora kinase inhibitor, alisertib
The overlap of the alisertib TCS and results from differential expression testing using MAST between patient neoplastic cells grouped by predominant transcriptional state enrichment was calculated to identify potential state-specific expression perturbed by alisertib (Fig. 4a). Single-cell Spearman’ ρ values for the alisertib TCS were extracted from the L1000 small molecules versus single-cell correlation matrix. Neoplastic patient cells were divided into predicted alisertib resistant and alisertib sensitive populations wherein cells concordant (Spearman’s ρ > 0) with the alisertib TCS were identified as a predicted resistant population, and those discordant (Spearman’s ρ < 0) with the alisertib TCS were identified as a predicted sensitive population (Fig. 3a, b). The predicted shift in cell states was then calculated, demonstrating an increase in MES-like cells and a decrease in NPC-like cells (Fig. 3c). Proportions across the sourcedataset, initial occurrence or recurrence, source patient, and cell cycle phase were also compared (Supplementary Fig. 4a–d). Proportion shifts of cells predominantly in each transcriptional state within each individual were compared using ANOVA with a Games-Howell post-hoc test. Resulting adjusted p-values were used to visualize the relative probability of predicted directional state abundance shifts (Supplementary Fig. 4d). Differential expression testing between predicted resistant and sensitive populations was performed using MAST, and cells were scored for their enrichment of overexpressed genes using Seurat’s AddModuleScore function (Supplementary Fig. 5f).

GBM cell orthotopic xenograft, alisertib treatment, and scRNAseq library creation
Human patient-derived GBM cells (GBM22) expressing GFP obtained from the Mayo Clinic Brain Tumor Patient-Derived Xenograft National Resource were implanted in the brains of nude mice (Fig. 3d)22,76,85. Cells were visualized to confirm GFP expression prior to intracranial implantation into mice. Briefly, 6–8 weeks old Nu/Nu mice (Charles River Laboratory) were anesthetized with ketamine (100 mg/kg) and xylazine (10 mg/kg), and a hole was drilled through the skull (1 mm anterior and 2 mm lateral to bregma). Using a Hamilton syringe, 150,000 GBM22 cells suspended in PBS were injected stereotactically at a depth of 2 mm. The skin incision was closed using surgical glue with sutures or wound clips and mice were given buprenorphine as analgesic. Procedures were approved by the University of Miami IACUC, protocol number 18-014. All procedures were performed in accordance with guidelines and regulations, including ARRIVE guidelines.

Following tumor implantation, mice recovered for a period of 7 days. GBM22 orthotopic tumor-bearing mice were treated intraperitoneally with either alisertib (25 mg/kg) (drug formulation: 5% DMSO/ 30% PEG400/ 5% Tween-80/ 60% PBS) or vehicle control 5 days per week and monitored daily for signs of cognitive decline and body weight (BW). Tumors were collected at the time of initial signs of cognitive decline or 15% of BW loss.

Tumor cells from three individual mice within each treatment group were captured and pooled for single cell RNA-sequencing. Tumor-bearing mice were perfused with PBS, and their brains were removed. Tumor was isolated and dissociated using the Worthington Biochem Papain dissociation kit. Cells were washed two times using ice-cold PBS with 0.1% BSA (w/v), then checked for viability using AO/PI dye and a Nexcelom K2 Cellometer. Cells from the dissociation suspension were loaded into a 10X Genomics’ Chromium chip using 10X Genomics’ NextGEM 3′ chemistry targeting the capture of ~10,000 cells per sample. With each individual animal experiment, paired-end library sequencing of DMSO and alisertib-treated single-cell libraries was performed in tandem on a single S1 NovaSeq flow cell, with asymmetric read lengths as recommended by 10X Genomics.

Pre-processing and quality control of in vivo scRNAseq data
Sequencing base calls were input into Cell Ranger V3.0.2 for demultiplexing and fastq generation. Alignment and gene quantification were first performed in tandem to a dual-species reference transcriptome for both hg19 (human) and mm10 (mouse), obtained from 10X Genomics (refdata-cellranger-hg19-and-mm10-3.0.0.). Using the R package Seurat (V4.0.1), the percentage of raw UMI counts within each cell aligning to either mm10 or hg19 transcriptomes was calculated and used in tandem with SNN clustering to separate out human GBM22 cells from contaminating mouse cells (Fig. 3e). GBM22 cell barcodes were then used to isolate these human cells from the data aligned to the GrCH38 (human) transcriptome for downstream gene expression analysis (Fig. 3f). Outlier cells were removed based on UMI count, detected gene count, and percent of mitochondrial transcript counts. Cells with less than 200 or more than 6000 detected genes, or with more than 10% of mitochondrial genome alignment were removed for downstream analysis. Pass-filter data was comprised of 16,102 GBM22 single-cell transcriptomes with 10,394 measured features.

Generation of a CT-179 GBM cell transcriptional response signature
GBM8 cells (Primary, IDHwt, Female, Classical, EGFR amp, MGMT methylated) were obtained from the Mayo Clinic Brain Tumor Patient-Derived Xenograft National Resource. GBM8 cells were treated with 100 nM CT-179, 200 nM CT-179, or vehicle control, and RNA was extracted following a 24h treatment period. RNA sequencing was performed on the NovaSeq 6000. FastQ files were aligned to the hg19 human genome to determine raw gene counts95. Differential expression between 200 nM CT-179 and vehicle-treated GBM8 cells was performed using DESeq2 (Supplementary Fig. 10a). Differentially expressed genes retained within the signature showed a logFC greater than 0.425, calculated as two standard deviations from the mean of all logFC values, and an adjusted p-value less than 0.05 (Fig. 4b). Gene ontology analysis was performed on this filtered signature using clusterProfiler (Supplementary Fig. 10b)96–100. Enrichment of the signature in 100 nM CT-179 treated, 200 nM CT-179 treated, or vehicle-treated GBM8 samples was calculated using the R package Singscore and compared to confirm the dose-dependence of our signature (Supplementary Fig. 10c).

In vivo evaluation of CT-179 and Depatux-M in combination
Orthotopic implantation of GBM6-eGFP-FLUC2 was performed using a previously established protocol76. Mice were anesthetized with ketamine and xylazine (100 mg/kg and 10 mg/kg, respectively), then placed in the Stoeling #51615 frame for sterotactic guidance. The incision site was disinfected with Nolvasan and the skull was exposed. Using a Dremel with #8 bit, the skull was drilled at a point 1 mm anterior and 2 mm lateral to bregma. 150,000 cells (in a suspension of 50,000 cells/mL) were injected sterotaxically over 1 min using a 26G Hamilton syringe at a depth of 2 mm. One week following surgery, mice were treated with vehicle, CT-179, Depatux-M, or a combination of CT-179 and Depatux-M (n = 10 per group). CT-179 was formulated in 0.5% methylcellulose, 0.5% Tween 80 in water, and stored at 4ºC. Fresh batches were made on the first day of dosing, and on day 18 of dosing. CT-179 formulation was stirred at RT for 30 min prior to use, and was delivered at a dose of 128 mg/kg QD PO from days 8 to 35. Depatux-M (ABT-414) was formulated in PBS, and was delivered at a dose of 5 mg/kg IP on day 8 and day 22. Tumor formation was confirmed, and tumor volume was assessed using BLI (Fig. 5d). BLI was performed using an IVIS 2.50 cooled charged-coupled device camera (Xenogen 200 series) and Living Image software76. Mice were dosed with 10 mg/kg Cycluc 1 by intraperitoneal injection and imaged 10 min later under isoflurane anesthesia76. A moribund state was defined by outward neurological deficits (altered gait, hunched posture, circling, poor body condition score). Survival differences were assessed using the Log-Rank test and Kaplan–Meier analysis (Fig. 5e).

Quantification and statistical analysis
Ranking of L1000 compounds against GBM cell transcriptional state disease signatures. The L1000 assay contains 978 landmark genes that can be measured to algorithmically infer the expression of 11,300 additional genes77,78,101. Using a previously published algorithm to rank discordance of L1000 perturbations, small molecules’ response signature discordance with scRNAseq-derived disease signatures were calculated for each compound in the dataset13,14,77,101. This disease-specific discordance score was calculated for all compounds within the L1000 library as a ratio of the number of genes with discordant expression (opposite gene expression direction/sign) between the TCS and the cell state disease signature over the number of genes with concordant expression (same gene expression direction/sign) between the TCS and the disease signature (Annotation Bar Fig. 2b, and Supplementary Data 4)13.

Ranking of single compounds against individual cell expression programs. The filtered single-cell counts matrix output by Cell Ranger was log-normalized and scaled using the R package Seurat. Using the normalized consensus signatures from the LINCS L1000 2017 dataset, predicted relative sensitivity of each single-cell to each L1000 perturbation (connectivity) was measured as the Spearman’s ρ between a compound’s consensus signature and each individual cell’s scaled gene expression of the L1000 genes present in each compound’s consensus signature (Fig. 2a). Pearson correlation analysis was applied to the compound versus cell coefficient matrix to visualize small molecule clustering based on integration with the single-cell data (Fig. 2b). Mechanism of action annotations were obtained from the Broad Institute Drug Repurposing Hub (Supplementary Fig. 2)21. Single-cell drug connectivities were compared between cell states using a generalized linear model implemented through limma. The mean connectivities of significantly different FDA-approved oncology compounds in addition to alisertib (p < 0.05) for each cell state were then scaled and visualized by heatmap (Fig. 2h).

Tumor cell drug connectivity clustering analysis. Using the Drug Repurposing Hub, 63 FDA-approved oncology drugs with L1000-derived TCSs from the 2017 data release were identified. A pair-wise Spearman correlation matrix was generated on the calculated drug connectivity values of all tumor cells from each patient in the single-cell atlas for these drugs. Euclidean distance was calculated using the base R function dist(). Ward’s hierarchical clustering of the distance matrix was performed using the base R function hclust(). Within cluster variance (sum of squares) was visualized for each number of resulting clusters K using the R package factoextra102,16. K of 7 was selected for representative visualizations based on the resulting elbow plot.

Xenograft tumor cell stratification and state shift analysis. Xenograft tumor cells were assigned GBM cell state identities using GBM cell transcriptional state signatures (Neftel et al., 2019), as described above for patient data (Fig. 3g, h, and Supplementary Fig. 5b, c). Singscore-derived enrichment scores and proportions of cells predominantly within each transcriptional state were compared between DMSO vehicle and alisertib-treated xenografts (Fig. 3i, j). The analysis pipeline is available at https://github.com/AyadLab/scFOCAL-dataProcessing.

Patient and xenograft pharmacotranscriptomic comparisons. To investigate the relationship between cells in GBM22 xenografts and those captured directly from patient tumors, GBM22 cells from DMSO-treated xenografts were input into scFOCAL and were scored for correlation with all small-molecule TCSs. The mean compound-cell connectivity values (Spearman’s ρ) were calculated across each transcriptional state identity. Using Spearman’s ρ once more, these connectivity values were compared to the same mean connectivity values obtained using patient neoplastic cells (Supplementary Fig. 7a). To assess whether scFOCAL predictions of alisertib-resistant and sensitive cell identities from patient data adequately modeled the pharmacotranscriptomic changes occurring in vivo, we utilized a linear model to calculate global L1000 small molecule discordance shift across both patient and xenograft datasets. A pseudo-count of 1 was added to all values in the compound-cell connectivity matrices, which were then fit to a linear model using limma and voom. Patient-based compound-cell connectivity values were grouped by predicted resistance or sensitivity, and differential discordance was calculated across populations within each individual patient tumor. The xenograft cell discordance differential was calculated in the same manner. Resulting logFC values for each compound were compared across in silico and in vivo analyses using Spearman’s ρ (Supplementary Fig. 6b) (Supplementary Data 5).

Panobinostat TCS connectivity analysis in ex vivo GBM acute slice culture. Single-cell RNA sequencing data of 5 paired acute slice culture samples derived from patient glioblastoma resections treated with either DMSO vehicle or the HDAC inhibitor Panobinostat (0.2 μM) was obtained from GEO accession GSE14884251. Single-cell transcriptomes with fewer than 200 detected genes or with greater than the 75th percentile of reads detected, calculated respectively per sample, were kept for downstream analysis. The pass-filter dataset was comprised of 62,250 single-cell transcriptomes with 58,828 features. Individual samples were then integrated using reciprocal PCA within Seurat. SNN clustering and differential expression analysis were used to identify discrete cell types within the integrated dataset. A TCS was generated for panobinostat utilizing L1000 data from within the 2020 release as described above (Supplementary Fig. 7c). Cell connectivity to the panobinostat TCS was calculated across all cells within the dataset. A Wilcoxon rank sum test with continuity correction was utilized to compare panobinostat connectivity across treatment groups within the determined neoplastic (37,682 cells) and myeloid cell (13,017 cells) populations respectively, and was visualized by violin and area-normalized kernel density estimate plots created with ggplot2 (Supplementary Fig. 7d, e). Analysis of proportion change in Neftel et al. states was performed utilizing Q1 and Q3 quantile cut-offs for panobinostat TCS connectivity, and across treatment groups. Comparisons of proportion were performed using a paired Wilcoxon test. Shifts in proportion across these two groupings were compared using the Spearman correlation.

Combination scoring analysis. We determined the potential to apply scFOCAL connectivity analysis in ranking synergistic small molecule combinations. An scFOCAL combination score was implemented as the differential logFC of each small molecule’s connectivity between cells predicted to be sensitive or resistant to a reference small molecule, multiplied by that small molecule’s mean connectivity to cell population predicted to be resistant to the reference small molecule. To assess the utility of this approach, synergy screen data was obtained from Houweling et al. (2023) supplementary data52. Raw sum synergy BLISS scores from short-term treatment assays of 24 different cell lines in spheroid culture were used. Following filtering of the synergy screen data for compounds present within the L1000, 4 potential reference small molecules with more than two associated combinations with molecules also in the L1000 data were used for analysis. These reference molecules were erlotinib, lapatinib, pazopanib, and sunitinib. Neoplastic cells from the patient single-cell atlas were scored for connectivity of these and all partner molecules within the screen. Connectivity cut-offs for sensitive or resistant binning of the patient neoplastic cells was calculated as the mean of all cells’ connectivity for each reference molecule, respectively. Using limma, the differential connectivity across reference molecule sensitive and resistant cell populations, for each of the partner compounds present in the screen, were calculated for each individual patient within the atlas. The resulting logFC values from each patient for each partner molecule were averaged across all individual patients. These aggregated differential connectivities were then multiplied by the mean connectivity score of each respective small molecule to the resistant cell population to the respective reference small molecule across all patients. These scores were then scaled between 0 and 1. Calculated combination score and observed BLISS raw sum synergy for each combination was compared using Spearman correlation in each of the 24 cell lines tested, and across the aggregate mean BLISS scores for all cell lines for each combination.

Scoring of L1000 small molecules for predicted synergy with CT-179. CT-179 single-cell connectivity was calculated as the Spearman correlation between the normalized and scaled expression cells within our dataset against the logFC calculated for 200 nM CT-179 treatment of GBM8 cells. As above, sensitive and resistant cells were defined as those having CT179 connectivities less than and greater than 0 respectively. L1000 drug-cell connectivities were transformed using the Fisher’s Z transformation, and limma was used to calculate the differential connectivities to L1000 molecules between predicted CT-179 sensitive and resistant cells, using patient ID as a latent variable. Compounds with logFC < 0.5 and FDR < 0.05 were retained. The logFC values for these compounds were then weighted by multiplying by their respective mean resistant cell-drug connectivities to derive our scFOCAL-based combination index. Thus, a more positive combination index value indicates both decreased connectivity to that drug in the CT-179 reference population and an overall discordance of the drug with CT-179 predicted resistant cells. This combination index can be defined as follows:

Combination Index = logFC * ρ * I(logFC < 0 and ρ < 0)

Where a combination index is only calculated when both the differential connectivity (logFC) and the mean resistant cell connectivity (RCC) (ρ) are negative values. This combination index is a feature implemented within the scFOCAL R package and associated shinyGUI.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, Suter et al. developed a framework to identify cell type specific effective drugs for glioblastoma (GBM) using patient derived scRNAseq and drug induced transcriptional profiles from LINCS L1000 dataset. Their method, named ISOSCELES, uses cell type specific single cell transcriptomics data, and is based on the “drug induced signature reversal” principle. Using signature reversal methods with single cell data is an important step for the applicability of these methods, so the manuscript can be interesting for the systems biology community. The authors also provide an R package and a Shiny app for ISOSCELES. I have mainly methodological comments, and also some regarding the released tool.\nMajor:\nIn a recent paper (Koudijs et al., Validation of transcriptome signature reversion for drug repurposing in oncology, Briefings in Bioinformatics 2023) a potential issue regarding drug repurposing approaches in oncology has been raised. While tumor cells have a proliferation signature, oncology drugs have an anti-proliferation signature, so the observed anti-similarity can be aspecific, only a consequence of trivial difference between proliferative and anti-proliferative signature. The authors should discuss the potential effect of this fact on their results.\nMore specifically, I assume that alisertib strongly inhibits cell cycle. Based on previous studies (e.g.: the already cited Nefel et al. Cell 2019 paper), NPC-like cells have the largest proportion of cycling cells, while MES-like cells have the lowest proportion. Couldn’t explain the observed changes in cell proportions after alisertib changes simply the fact that alisertib has anti-similarity to cycling cell signature (more abundant in NPC-like cells), and inhibits this cycling cell population, thus decreasing their proportion? A more detailed analysis of cell cycle score in different populations could have help to discard this possibility.\nIn the “In silico drug connectivity analysis predicts intratumor compound discordance shifts seen in vivo” section it is hard to understand what is the consequence of these results. Could the authors define Observed / Predicted Corr. Coefficient FC more precisely (SFigure 5)? If I understand correctly, here the authors compare predicted cell proportion changes, for a given compound, and observed proportion changes after alisertib treatment? What is the reasoning behind this comparison? Also, the drugs with high correlation in SFigure 5 are strong inhibitors of cell cycle (tubulin polymerization / topoisomerase inhibitors), which can be connected to the two points above.\nMinor:\nThe GitHub repository for the manuscript is currently empty, I would suggest sharing it with the reviewers during the review process.\nFor ISOSCELES Shiny application some tutorial dataset would be helpful.", "reviewer_2": "The manuscript describes application of disease signature reversal to the transcriptional profiles of single cells as a prediction of drug sensitive and resistant populations in glioblastoma tumors. Based on scRNAseq from 11 GBM cases (a mix of primary and recurrent cases), the authors derive disease signatures for each of four GBM transcriptional cell states. The reversal of these disease signatures in L1000 small molecule transcriptomic dataset is then calculated to predict which subpopulation of GBM cells could be more sensitive or resistant to specific drugs.\nConceptually, this is a very interesting approach, which could help identify drug combinations or sequence in which drugs should be combined to more effectively target all subpopulations in GBM. However, the current manuscript provides only a single example of a drug for which the computational analysis was performed and only a single PDX for experimental validation. It is thus unclear how robust is the ISOSCELES platform.\nSpecific comments:\n1. In the example provided, it is mentioned that several drugs are predicted to target the Alisertib resistant population. It would be important to demonstrate if these predictions are correct by testing these drugs on resistant populations in vitro or combination treatment in vivo.\n2. The ISOSCELES figure (Fig. 5, Results, Proportion shift tab) suggests that there are two population expanding as resistant – MES1 and AC. How can this be reconciled with results presented in Fig. 4? Experiments performed on another cell line or two could help determine whether these expansions depend on intial fraction of the cell state in untreated tumor.\n3. The current ISOSCELES model is built on 11 datasets, of which 4 are recurrent GBM cases. Would the predictions be better if only primary cases were used? Would they look different for only recurrent cases? There is several publicly available datasets based on 10xGenomics platform that could be used to test this: ex. 16 patients in PMID: 32641768 , 10 patients in PMID: 35122077, 11 patients in PMID: 31901251 and many others.\n4. The overall robustness of the study would be significantly strengthened by additional PDXs tested for their Alisertib response. In a heterogeneous disease, such as GBM, a single cell line/PDX might not be representative.\nMinor comments:\n1. Please include a summary of the number of cells and gene count per cell to describe the single-cell datasets used in this study.\n2. CytoTRACE-based assignment of malignant and non-malignant state should be briefly explained in the text to clarify what is presented in Fig. 1d.\n3. Reference to Fig.3a in the text is missing.\n4. Fig. 3 – would supervised clustering based on disease signature reversal specific to each cell state show significant differences between states? In other words, is L1000 dataset integration with expression data helpful as a discovery tool?\n5. Fig. 4g – is the change in percent composition significantly different after Alisertib treatment? Please provide statistical test p value.", "reviewer_3": "In this manuscript, the authors present an interesting general approach: the development of an algorithm that allows to predict sensitivity of glioblastoma cell subpopulations (i.e., NPC-, OPC-, MES- and AC-like cells) to certain drugs and drug classes. Tumor heterogeneity is a formidable challenge in the treatment of glioblastoma, and a better idea of which transcriptional tumor cell subpopulations repond to which drugs/classes can help to develop tailored multi- or oligo-drug strategies; potentially even identifying one drug that is able to target all, or at least the most important, tumor cell subpopulations.\nHaving said that, this manuscripts falls short of providing evidence that the computational in silico approach presented here is of any practical (biological / therapeutic / mechanistic) value. The investigation of one single drug in one single cell line in vivo demonstrating one major change - a shift towards MES-like tumor cell subpopulation - is no evidence for the utility of this entire approach. A shift towards MES-like states is well known for many therapeutic interventions now, and far away from any proof of a specific alteration induced.\nThe authors would need to demonstrate that their predictions regarding multiple drug classes is really true, at least relevantly overlapping with experimental data. They need to study the response of all 4 subpopulations to a panel of drugs experimentally; subpopulations can be studied by enrichment/depletion analyses, better individually after FACS sorting or, ideally, and something this reviewer would strongly recommend, using in vivo reporter systems or bar code systems that allow to distinguish all 4 subpopulations. This should be done in state-of-the-art in vitro assays. After that, the authors need to demonstrate in vitro and in vivo that the most promising compound (mix) is indeed providing a therapeutic advantage compared to other approaches. The results need to be clear and convincing to make this study to a relevant advance in the field. The authors need to provide evidence that ISOSCELES is of any de facto meaning and usability for the community. This most important part is simply missing from this study.\nMore minor point: the presentation of the results and figure legends is often confusing and should be significantly improved. All figure panels should be cited in the order of appearance; it should be explained in the results section what they actually show; and the reader needs to get more explanatory informations there, too."}, "author_response": "Response to Reviewers:\nReviewer 1: Computational pharmacology, omics (Remarks to the Author): In this manuscript, Suter et al. developed a framework to identify cell type specific effective drugs for glioblastoma (GBM) using patient derived scRNAseq and drug induced transcriptional profiles from LINCS L1000 dataset. Their method, named ISOSCELES, uses cell type specific single cell transcriptomics data, and is based on the “drug induced signature reversal” principle. Using signature reversal methods with single cell data is an important step for the applicability of these methods, so the manuscript can be interesting for the systems biology community. The authors also provide an R package and a Shiny app for ISOSCELES. I have mainly methodological comments, and also some regarding the released tool.\nIn a recent paper (Koudijs et al., Validation of transcriptome signature reversion for drug repurposing in oncology, Briefings in Bioinformatics 2023) a potential issue regarding drug repurposing approaches in oncology has been raised. While tumor cells have a proliferation signature, oncology drugs have an anti-proliferation signature, so the observed anti-similarity can be aspecific, only a consequence of trivial difference between proliferative and anti-proliferative signature. The authors should discuss the potential effect of this fact on their results. More specifically, I assume that alisertib strongly inhibits cell cycle. Based on previous studies (e.g.: the already cited Nefel et al. Cell 2019 paper), NPC-like cells have the largest proportion of cycling cells, while MES-like cells have the lowest proportion. Couldn’t explain the observed changes in cell proportions after alisertib changes simply the fact that alisertib has anti-similarity to cycling cell signature (more abundant in NPC-like cells), and inhibits this cycling cell population, thus decreasing their proportion? A more detailed analysis of cell cycle score in different populations could have help to discard this possibility.\nWe thank the reviewer for their important perspective on the consideration of the cell cycle in the use of disease signature reversal approaches. We have made our best effort to take cell-cycle-related gene expression differences into consideration throughout our different analyses. First, enrichment scores for Neftel states in both the presented patient GBM scRNAseq atlas and in vivo studies were calculated using cell cycle-regressed and scaled data. Regression was performed within the ScaleData() function of Seurat on G2M and S phase module scores calculated as in Tirosh et al., 2015, as implemented within the CellCycleScoring() function of Seurat, using updated G2M and S phase gene sets as provided within the Seurat package. This function also bins cells into their cell cycle phase identities based on the G2/M and S phase module scores, and the proportions of cells in each phase within predicted alisertib-sensitive or resistant cell populations are depicted in Supplementary Figure S4d. Further, as described in the methods section, scFOCAL cell-drug connectivity is calculated against cell-cycle-regressed single-cell expression.\nSupplementary Figure S4d: Sankey diagrams depicting proportions of predicted sensitive and resistant cell populations across cell cycle phase.\nComparing the cell-cycle differences between the predicted sensitive and resistant and sensitive cell populations, it can be seen that an abundance of cells determined to be in G1 are predicted to be alisertib-sensitive, also that an abundance of cells determined to be in S-phase are predicted to be resistant. Importantly, alisertib has been reported to arrest cells in S and G2/M phases. In both the predicted resistant population from the patient GBM scRNAseq atlas and the alisertib-treated xenografts, the proportions of cells in both G2/M and S are increased relative to the predicted alisertib-sensitive population or the DMSO-treated xenografts, respectively. Taken together, our prediction and result in vivo that NPC-like cells are preferentially depleted by alisertib does not seem to be an artifact related to cell-cycle-related gene expression or a decrease in proliferative cells.\nFurther, to demonstrate more effectively that our platform is identifying biology beyond the cell cycle, we now include a second in vivo experiment to demonstrate that our approach is useful for predicting drugs targeting different GBM cell transcriptional states in tumors. We included a new in vivo study with the Olig2 inhibitor CT-179, as we expected this compound to preferentially target the OPC-like GBM cell state. Our new bulk RNA sequencing data also indicates an increase in proliferation and cell-cycle-related pathways in GBM8 Mayo PDX cells after treatment with CT-179 (Supplementary Figure S9a-b).\nSupplementary Figure S9: Generation of a GBM CT-179 transcriptional response signature.\na. Volcano plot of differential expression between 200nM CT-179 and vehicle-treated GBM8 cells calculated via DESeq2. b. Barplot depicting gene ontology analysis of the filtered CT-179 transcriptional response signature. Bars are colored by the respective ontology used. Selected GO terms highlighted include ‘microtubule’, ‘microtubule binding’, ‘tubulin binding’, and ‘microtubule motor activity’.\nVia scFOCAL, as expected, CT-179 is predicted to target OPC-like cells. Interestingly, though, CT-179 is predicted to yield an AC-like resistant cell population. Through the prioritization of EGFR and tubulin targeting molecules by our combination prioritization strategy, which we previously included in our original submission as Supplementary Figure S6, we validate our predictions from a CT-179 response signature by demonstrating synergy in vivo with CT-179 and Depatux-M, an anti-EGFR MMAF (tubulin inhibitor) antibody-drug conjugate (Fig. 4, Fig. 5, Fig. S9). While further mechanistic investigation of the mechanism behind this synergy is required, we hypothesize that Depatux-M is targeting a more proliferative AC-like population of cells that express EGFR. We sincerely hope that this new use-case demonstration and in vivo validation of our scFOCAL platform’s performance ease caution on cell-cycle related bias in its application.\nImportantly, our alisertib experiment still fits the altered narrative of this manuscript through demonstration that in silico perturbation of our patient GBM scRNAseq atlas accurately predicts differential connectivities to other drugs between the vehicle and alisertib-treated PDX tumors (Fig. 3k, below), which led us to investigate the potential for our framework to prioritize synergistic combinations. (Also see added in vitro validation of our scFOCAL combination score performed using data obtained from Houweling et al, 2023 (Neuro Oncology Advances)).\nFigure 3k. Scatterplot of L1000 small molecule TCS’s depicting predicted correlation shift (log FC) vs. observed correlation shift (log FC) in alisertib-treated xenografts. Differential small molecule correlations were calculated using limma.\nIn the “In silico drug connectivity analysis predicts intratumor compound discordance shifts seen in vivo” section it is hard to understand what is the consequence of these results. Could the authors define Observed / Predicted Corr. Coefficient FC more precisely (SFigure 5)? If I understand correctly, here the authors compare predicted cell proportion changes, for a given compound, and observed proportion changes after alisertib treatment? What is the reasoning behind this comparison? Also, the drugs with high correlation in SFigure 5 are strong inhibitors of cell cycle (tubulin polymerization / topoisomerase inhibitors), which can be connected to the two points above.\nWe thank the reviewer for their important comments on the presentation of the drug connectivity analysis. For the original Supplementary Figure S5, we sought to characterize the pharmacotranscriptomic differences between the cells predicted to be sensitive or resistant to alisertib treatment using scFOCAL. Importantly, these predicted sensitive and resistant populations should be considered analogous to the DMSO-treated or alisertib-treated cells, respectively, from tumors in our in vivo experiment.\nIn our revisions, we have made this point a focus, as herein lies the potential to utilize differential drug connectivity to identify more effective drug combinations. As such, we have moved this figure from the supplement to a main figure (Fig. 3k). Here, we identified small molecules with a differential drug-cell connectivity (cell expression vs drug TCS) between predicted alisertib-sensitive (Alisertib TCS Connectivity < 0) and -resistant (Alisertib TCS Connectivity > 0) cell populations using limma, within each patient sample in the dataset. These differentials represent our predictions of alisertib-induced increases or decreases in sensitivities to other L1000 molecules. Calculated results were aggregated from each patient to a consensus of which small molecules had conserved shifts across all patients’ predicted sensitive and resistant cell populations. We then repeated this analysis using the data from our in vivo alisertib experiments to identify small molecules with differential drug-cell connectivity scores across actual treatments (DMSO vs alisertib). In our new Fig 3k (previously supplemental), we demonstrate a significant correlation across all L1000 molecules of their predicted connectivity differentials determined in silico vs those differentials calculated from actual treatment data. This analysis serves to show that the cells which we are predicting to be sensitive or resistant from our patient dataset are transcriptionally (in the context of L1000 signatures) similar to the tumor cells that are sensitive or resistant (persisting after treatment) respectively in our xenograft experiments, and serves as a basis for our new scFOCAL combination index algorithm. Importantly, this differential drug connectivity measure is utilized within our combination scoring approach, as depicted below in Fig. 4a.\nFigure 3: In silico perturbation of GBM tumor cell scRNAseq data using an L1000-derived alisertib TCS predicts an NPC-like to MES-like tumor response confirmed in vivo\na. Histogram of single-cell alisertib TCS correlations (connectivities). Cells with negative correlation coefficients (blue) are predicted to be sensitive to alisertib, while cells with positive correlation coefficients (red) are predicted to be resistant to alisertib. b. Hierarchy plot of patient GBM tumor cells colored by predicted sensitivity (ρ < 0) or resistance (ρ > 0) to alisertib. c. Bar plot depicting mean shift in proportions of cells in resistant vs. sensitive populations within individual patient tumors. Error bars depict the standard error of the mean percentage differences of individual patient tumors. d. Schematic of in vivo experiments. e. UMAP plot of pre-filter single-cell transcriptomes colored by treatment with either alisertib or DMSO vehicle control. f. UMAP of captured cells from orthotopic xenografts colored by percent alignment to the human transcriptome (hg19). g. Dot plot of pass-filter, human xenograft single-cell transcriptome expression of Neftel et al. signatures, grouped by assigned transcriptional state identity as assigned based on predominant transcriptional state module expression. h. Two-dimensional hierarchical representation of GBM22 xenograft cells’ relative enrichment scores for GBM cell transcriptional state modules. Cells are colored by assigned transcriptional state identity. i. Bar plot of mean enrichment shift of transcriptional state signatures in alisertib-treated xenograft cells normalized to DMSO vehicle-treated xenograft cells. Error bars represent 95% confidence interval. (Wilcoxon with Benjamini-Hochberg correction: *p-adjusted = 8e-9; **p-adjusted = 1.1e-60;***p-adjusted = 3.2e-238;NS p-adjusted = 0.79) j. Alluvial plot depicting shift in relative proportion of transcriptional state identities in alisertib and DMSO vehicle control treated xenografts. k. Scatterplot of L1000 small molecule TCS’s depicting predicted correlation shift (log FC) vs. observed correlation shift (log FC) in alisertib-treated xenografts. Differential small molecule correlations were calculated using limma.\nThe thought behind this analysis is now more complete and leads to subsequent predictions of potential synergistic combinations based on a reference small molecule, where the reference molecule connectivity divides single-cell input into resistant and sensitive populations, used in the calculation of an scFOCAL-based combination score discussed in more detail further down in this rebuttal. Importantly, we have revised our figures to reflect the reasoning behind this analysis more clearly, displaying our in silico prediction vs in vivo results in the same figure as our predictions and experimental data. Further, we have made a strong effort to better explain the significance of this approach and the figure in the text. Also importantly, additional analysis and figures demonstrate the utility of this scoring by effectively ranking the synergy of combinations as seen in vitro, and by identifying a novel combination therapy which is more effective than either monotherapy alone in vivo.\nFigure 4a. Diagram of workflow for identifying CT-179 combinations using scFOCAL and cell-drug connectivity.\nMinor: The GitHub repository for the manuscript is currently empty, I would suggest sharing it with the reviewers during the review process. For ISOSCELES Shiny application some tutorial dataset would be helpful.\nWe thank the reviewer for pointing this out. The GitHub repository for the scFOCAL application itself contains the scFOCAL R package source code and is available at https://github.com/AyadLab/scFOCAL. Within the readme file for this git, there are links to download tutorial datasets for use with the application. We have also since deposited all processing pipelines, including those for analyses performed during revision, within https://github.com/AyadLab/scFOCAL-dataProcessing. While the application is available on shinyapps.io, we suggest installing it as an R package locally as per the readme file in the GitHub repository, as analysis with the tutorial dataset will run more efficiently, and shinyapps.io limits the needed memory usage.\nScreenshot of the scFOCAL GitHub README\nReviewer 2: The manuscript describes application of disease signature reversal to the transcriptional profiles of single cells as a prediction of drug-sensitive and resistant populations in glioblastoma tumors. Based on scRNAseq from 11 GBM cases (a mix of primary and recurrent cases), the authors derive disease signatures for each of four GBM transcriptional cell states. The reversal of these disease signatures in L1000 small molecule transcriptomic dataset is then calculated to predict which subpopulation of GBM cells could be more sensitive or resistant to specific drugs. Conceptually, this is a very interesting approach, which could help identify drug combinations or sequences in which drugs should be combined to more effectively target all subpopulations in GBM. However, the current manuscript provides only a single example of a drug for which the computational analysis was performed and only a single PDX for experimental validation. It is thus unclear how robust is the ISOSCELES platform.\nSpecific comments: Minor comments:\nFig. 3 – would supervised clustering based on disease signature reversal specific to each cell state show significant differences between states? In other words, is L1000 dataset integration with expression data helpful as a discovery tool?\nRegarding the reviewer’s comment on supervised clustering, we now show in new Figure 2d-h that clustering cells on their scFOCAL calculated drug connectivities is, in fact, very interesting. Importantly, we show that even in the context of only 63 FDA-approved oncology drugs, we can cluster and stratify cells in a way that reflects not only Neftel et al transcriptional state, but also other features such as differentiation state, and not surprisingly also by cell cycle. Importantly, we also demonstrate that in the context of these FDA-approved molecules, they are clustered by their connectivities to GBM tumor cells in such a way that reflects their respective mechanisms of action. (Supplemental Figure S2). Lastly, we demonstrate in Figure 2h that alisertib represents an FDA-approved (Orphan approval) compound which has significantly lower connectivity to NPC-like tumor cells in our dataset relative to the other states, as determined using limma.\nFigure 2: Integration of single-cell expression and small molecule L1000 TCS signatures permits clustering of both compounds and cells by reversal of GBM cell transcriptional state-specific disease signatures. a. Schematic of single-cell sensitivity and resistance scoring. b. Correlation matrix depicting similarities of L1000 small molecule TCSs by their connectivity to all individual cells within our single-cell atlas. Row annotations depict compound TCS scores for the reversal of AC-, MES-, NPC-, and OPC-like disease signatures. c. Network plot of select L1000 small molecules colored by mechanism of action. Connections indicate a Pearson’s ρ > 0.7 between small molecules by their calculated discordances against all single-cells in the GBM dataset as represented in (b). d. Elbow plot depicting the within-cluster sum of squares by number of clusters k for GBM tumor cells from patient sample GBM21. e. Elbow plot depicting the within-cluster sum of squares by number of clusters k for GBM tumor cells from patient sample SM006 (Johnson et al. dataset). f-g. Correlation matrices depicting pairwise Spearman correlations of single GBM tumor cells from individual patient tumor samples (GBM21, in-house) (f), SM006 (Johnson et al.) (g) by their connectivity values to 63 FDA-approved oncology drug TCSs. Row annotation bars depict hierarchical clustering identities (k = 7). Column annotations depict CytoTRACE scores, cell cycle phase, and assigned Neftel et al. state identity. h. Heatmap showing the scaled single-cell drug connectivities for FDA-approved oncology compounds plus alisertib that were significantly different between Neftel et al. cell states (linear model, BH-adjusted p < 0.05). Heatmap color indicates predicted cell state sensitivity to each compound (red = relative sensitivity, blue = relative resistance). Alisertib, used in an in vivo experiment in the next figure, is highlighted by a red arrow.\nSupplementary Figure S2: FDA-approved oncology drugs can be clustered on their L1000 TCS connectivity to GBM tumor cell expression. a. Elbow plot depicting the within-cluster sum of squares by number of clusters k for connectivity values of 64 clinical oncology drug TCSs with all GBM tumor cells from 11 patients in the scRNAseq atlas. b. Correlation matrix heatmap depicting pairwise Spearman correlations of 64 clinical oncology drugs based on their connectivity values to all GBM tumor cells from the 11 patients in the scRNAseq atlas.\nIn addition to the use of our framework as a discovery tool to stratify cell states by their relative drug connectivities, we also demonstrate its utility in the drug discovery pipeline – leveraging scFOCAL to identify an effective combination with a new promising clinical candidate.\nFig. 4g – Is the change in percent composition significantly different after Alisertib treatment? Please provide the statistical test p-value.\nWe thank the reviewer for this important comment. As individual xenograft samples with shared treatment were pooled together, we are unable to perform a statistical significance assessment in regards to the percent composition before and after treatment. However, at the single-cell level, we can assess the shifts in relative Neftel et al state enrichment by comparing the enrichment scores for alisertib-treated tumor cells normalized to the mean enrichment scores of DMSO-treated tumor cells. Using a Wilcoxon signed-rank test and Benjamini-Hochberg multiple testing correction, we demonstrate significant enrichment shifts in AC-like, MES-like, and NPC-like module expression. Notably, OPC-like module enrichment was not significantly altered by alisertib treatment in vivo. These statistics have been added to Figure 3i as pictured below, and the prior panel depicting percent composition change as a bar plot has been removed in favor of the alluvial barplot in Fig 3j.\nFigure 3i: Bar plot of the mean enrichment shift of transcriptional state signatures in alisertib-treated xenograft cells relative to the mean enrichments of DMSO vehicle-treated xenograft cells. Error bars represent the 95% confidence interval. (Wilcoxon with Benjamini-Hochberg correction performed for each state’s normalized shift vs 0 (steady-enrichment): *p-adjusted = 8e-9; **p-adjusted = 1.1e-60;***p-adjusted = 3.2e-238; NS p-adjusted = 0.79) Figure 3j: Alluvial plot depicting shift in relative proportion of transcriptional state identities in alisertib and DMSO vehicle control-treated xenografts.\nIn the example provided, it is mentioned that several drugs are predicted to target the Alisertib-resistant population. It would be important to demonstrate if these predictions are correct by testing these drugs on resistant populations in vitro or combination treatment in vivo.\nWe thank the reviewer for this important perspective on the utility of the scFOCAL platform to identify combinations of small molecules that are predicted by the platform to target unique cell populations present within the same GBM tumor.\nIn our revisions, we’ve used the alisertib experiment to show that our platform predicts the relative sensitivites of alisertib resistant cells to different drugs (the compounds that we identify to be discordant with the predicted alisertib-resistant population from our patient dataset, also show increased discordance to cells in alisertib-treated tumors relative to vehicle treated tumors) to frame the potential to use our framework to design combinations. In summary, scFOCAL predicts differential drug connectivities occurring in vivo. We then go on to validate the significance of this demonstrating prediction and observed synergy correlation in vitro and in vivo using different small molecules. As our in vivo experiments with alisertib showed the OPC-like expression module was not significantly altered by alisertib, we hypothesized that the novel Olig2 inhibitor CT-179 would be predicted to target the OPC-like state using scFOCAL. Eventually, we demonstrate scFOCAL’s ability to identify novel synergistic combinations in vivo using this inhibitor.\nTo this end, during revisions, we first refined the algorithm within the scFOCAL platform for identifying potentially synergistic combinations and visualization of its output. We applied this algorithm to the patient scGBM atlas used to generate combination scores for compounds used within the high-throughput synergy screen published in Houweling et al, 2023. This synergy screen was performed by the original authors on 25 different glioma stem cell lines in spheroid culture. Raw sum synergy measured by BLISS was compared to our scFOCAL combination index for tested molecules in combination with either erlotinib, lapatinib, pazopanib, or sunitinib. In both cases, overall, the combination index was significantly correlated with observed BLISS raw sum synergy. For the majority of cell lines tested, this combination index was very predictive of the relative synergy of the combinations tested. The most potent combinations in these instances were of those made of compounds with distinct TCS-connected cell populations.\nSupplementary Figure S8: scFOCAL combination scoring predicts relative synergy of small molecule combinations in vitro. a. Stacked bar plots depicting the proportion of cells predicted to be sensitive or resistant across individual patient GBM tumors to reference small molecules erlotinib, lapatinib, pazopanib, and sunitinib, respectively. b. Heatmap of BLISS raw sum synergy of small molecule combinations obtained from Houweling et al., 2023. The top bar plot annotation depicts scFOCAL calculated combination score using labeled reference molecules depicted on the bottom annotation bar. Sidebar annotations depict Spearman’s ρ and p-value between scFOCAL combination score and observed BLISS raw sum synergy for each cell line. c. Scatterplot of scFOCAL calculated combination score and mean BLISS raw sum synergy across all cell lines tested for small molecule combinations screened in Houweling et al., 2023 (Spearman's ρ = 0.7, p-value = 0.01). Points are colored by the reference molecule used for combination scoring. The shaded gray area around the line of best fit depicts the standard error.\nBuilding on this analysis of prior in vitro screens, we moved to demonstrate an effective combination prioritized by our platform in a new in vivo experiment. Using a recently developed OLIG2 inhibitor, CT-179, a GBM cell line response signature was generated using bulk RNA sequencing (Figure 4a-b, Supplementary Fig. S9a-c). We demonstrate that this signature is enriched by CT-179 treatment in a dose-dependent manner, and then utilize this signature to predict relatively sensitive or resistant cell populations using our patient scRNAseq dataset (Figure 4c). Our predictions suggest that OPC-like cells would be sensitive, as would be expected, and the AC-like population to be relatively more resistant (Figure 4e). Using our combination algorithm, we identified that compounds targeting EGFR and tubulin scored highly for synergy with CT-179 (Figure 5a-c). With this information, we demonstrate the combination of CT-179 and Depatux-M (a-EGFR /MMAF), an antibody drug conjugate combining the targeting of EGFR (over-expressed in AC-like cells) expressing cells with a tubulin poison, to increase survival in vivo relative to each monotherapy alone (Figure 5d-e).\nFigure 4: The integration of a bulk-derived transcriptional response signature for the OLIG2 inhibitor CT-179 predicts targeting of OPC-like GBM cells. a. Diagram of workflow for identifying CT-179 combinations using SCFOCAL and cell-drug connectivity. b. Heatmap of GBM8 cells treated with vehicle or 200nM CT-179 for 24 hours. Columns represent biological replicates. c. Bar plot depicting the proportion of tumor cells within each patient tumor predicted to be sensitive (CT-179 response ρ < 0) or resistant (CT-179 response ρ > 0) to CT-179 treatment. d. Hierarchy plot of GBM tumor cells arranged by their relative expression of Neftel et al. states and colored by predicted sensitivity or resistance to CT-179 treatment. e. Violin and box plot of GBM cell CT-179 connectivity (ρ) grouped by assigned GBM cell transcriptional state.\nFigure 5: An scFOCAL combination index predicts a synergistic combination of an OLIG2 inhibitor CT-179 with Depatux-M (ABT-414), an anti-EGFR antibody MMAF drug conjugate. a. Volcano plot depicting the results from limma-based differential drug connectivity between predicted CT-179 sensitive and resistant GBM cells. Patient ID was used as a covariate in the model design. b. Barplot of L1000 small molecule mean resistant cell connectivity (RCC) to predicted CT-179 resistant cells. Color depicts this same value. c. Scatterplot of L1000 small molecules plotted by resistant vs. sensitive differential connectivity and mean CT-179 resistant cell connectivity. Colors depict the calculated scFOCAL combination index, the product of multiplying the differential connectivity log FC values by the mean CT-179 resistant cell connectivity values for each molecule. Compounds highlighted in (a), (b), and (c) are varlitinib, an EGFR inhibitor, and indibulin and docetaxel, which act through inhibition of tubulin. d. Bioluminescence signal quantification of GBM6-eGFP-FLUC2 orthotopic xenograft tumors (means SD, Combo vs Vehicle: p.adj < 0.005 from Day 14, Combo vs. CT-179 monotherapy: p.adj < 0.05 from Day 11, Adjusted p-values from multiple t-tests with Holm–Šídák correction to control the family-wise error rate) e. Kaplan-Meier survival curves of mice bearing GBM6-eGFP-FLUC2 orthotopic xenografts treated with indicated therapies (n = 10 per group). MS: median survival. P-values determined using the Log-rank (Mantel-Cox) test.\nSupplementary Figure S9: Generation of a GBM CT-179 transcriptional response signature. a. Volcano plot of differential expression between 200nM CT-179 and vehicle-treated GBM8 cells calculated via DESeq2. b. Barplot depicting gene ontology analysis of the filtered CT-179 transcriptional response signature. Bars are colored by the respective ontology used. Selected GO terms highlighted include ‘microtubule’, ‘microtubule binding’, ‘tubulin binding’, and ‘microtubule motor activity’. c. Barplot of singscore-calculated enrichment for the directional CT-179 transcriptional response signature in 200nM CT-179-treated, 100nM CT-179-treated, and vehicle-treated GBM8 cells after 24 hours. Significance bars depict p-value following pairwise t-tests (n = 3 per treatment group). d. Box and violin plots depicting the expression of EGFR and OLIG2 in patient GBM tumor cells grouped by assigned Neftel et al. GBM cell transcriptional state.\nThe scFOCAL figure (Fig. 5, Results, Proportion shift tab) suggests that there are two populations expanding as resistant – MES1 and AC. How can this be reconciled with results presented in Fig. 4? Experiments performed on another cell line or two could help determine whether these expansions depend on the initial fraction of the cell state in untreated tumor.\nThank you for bringing this to our attention. This figure was meant to be illustrative and was based on an outdated and sub-sampled testing dataset, prior to calculation of overall MES and NPC scores (i.e MES1/2, NPC1/2). We have updated Figure 6 to illustrate scFOCAL analysis of the sample data included with the GEO release.\nFigure 6: The scFOCAL framework is available for use as a shiny web application or an R package.\nThe current scFOCAL model is built on 11 datasets, of which 4 are recurrent GBM cases. Would the predictions be better if only primary cases were used? Would they look different for only recurrent cases? There is several publicly available datasets based on 10xGenomics platform that could be used to test this: ex. 16 patients in PMID: 32641768 , 10 patients in PMID: 35122077, 11 patients in PMID: 31901251 and many others.\nThank you for your comment. We sought to develop predictions regardless of newly diagnosed or recurrent tumor status, as many patients seek out trials upon tumor recurrence. As such, we had not tried to separate newly diagnosed from recurrent tumors since we did not have enough recurrent GBM patient tumors for these predictions. A single-cell atlas with an adequate number of paired longitudinal primary and recurrent samples would be required to adequately examine this further. However, such a dataset is not currently available. Our statistical analysis suggests we need at least 10 samples per group to identify differences, which is not available in our current dataset. Therefore, we will perform this analysis when more samples are available.\nThe overall robustness of the study would be significantly strengthened by additional PDXs tested for their Alisertib response. In a heterogeneous disease, such as GBM, a single cell line/PDX might not be representative.\nWe thank the reviewer for this important comment, and have since shifted the focus of our manuscript in response to this comment as well as others. To address this comment specifically, we have included additional analyses demonstrating our platform's performance in the context of different small molecules and perspectives.\nTo highlight the utility of our platform, we have first applied scFOCAL to publicly available scRNAseq data of patient GBM acute slice cultures treated with either DMSO or the HDAC inhibitor panobinostat, and showed that scFOCAL calculated panobinostat connectivity of single-cells correlates with panobinostat treatment ex vivo (Panobinostat TCS derived from the L1000 (Supp. Fig 7c)). Here we demonstrate that L1000-derived panobinostat connectivity calculated with scFOCAL is a predictor of whether a cell has been treated with panobinostat or not. Interestingly, in the original paper where this data originated, the authors identify that panobinostat treatment dramatically affected the tumor-associated myeloid cell population as well. scFOCAL cell-drug connectivity identifies this difference as well (Supp. Fig 7de).\nSupplementary Figure S7: scFOCAL connectivity analysis with an L1000-derived panobinostat TCS separates vehicle-treated from panobinostat-treated cells in both neoplastic and myeloid cell populations. a. UMAP of 62,250 single-cell transcriptomes from patient-derived acute slice culture samples obtained from Zhao et al., 2021, treated with DMSO or the HDAC inhibitor panobinostat (0.2 μM). b. UMAP of single-cell transcriptomes colored by discrete cell type. c. Heatmap showing z-scores of panobinostat-specific gene expression (vs. treatment naïve) for genes retained in the Panobinostat TCS across cell lines tested in the L1000 dataset. d. Violin plots showing panobinostat TCS connectivity of Panobinostat-treated vs. DMSO-treated cells in neoplastic and myeloid populations, respectively (Wilcoxon rank sum test with continuity correction, p-value < 2.2e-16). e. Area-normalized kernel-density estimate (KDE) plots of scFOCAL calculated panobinostat TCS connectivity colored by treatment with DMSO or panobinostat, across neoplastic and myeloid cell populations, respectively.\nWe too recognize that this extremely heterogeneous cancer likely exhibits varying responses to different drugs depending on the cell line utilized or the patient from whom a tumor originated. Thus, we have focused more on the aspect of combination design, where input data will change and thus will affect predictions, and that scFOCAL represents a step towards personalized medicine, rather than towards the identification of a one-size-fits-all treatment. As such, we have included a separate analysis of a different small molecule predicted to target a different GBM cell state, with a different predicted resistant GBM cell state. Leveraging this, we identify a novel combination and demonstrate its increased efficacy relative to each component as monotherapy alone.\nTaken together, even though we have not profiled additional PDXs with alisertib, we have included a variety of different use-cases and perspectives in support of our framework’s capability that we believe best suit its intended purpose.\nPlease include a summary of the number of cells and gene count per cell to describe the single-cell datasets used in this study.\nThank you for addressing this important oversight. We have since added cell and detected gene counts within the manuscript for patient, xenograft, and newly analyzed acute slice culture scRNAseq datasets.\nCytoTRACE-based assignment of malignant and non-malignant state should be briefly explained in the text to clarify what is presented in Fig. 1d.\nThank you for your comment. We have updated the methods with a brief explanation of CytoTRACE to clarify the importance of Fig. 1d, that the tumor cell population we identify within our dataset is de-differentiated relative to known non-tumor cells, and further that the tumor cell population is heterogeneous in the context of transcriptional diversity and thus developmental potential.\nReference to Fig.3a in the text is missing.\nThank you for bringing this to our attention. We have made sure that every figure, including supplementary figures, is mentioned and referenced within the text.\nReviewer 3: GBM therapy, preclinical models (Remarks to the Author): In this manuscript, the authors present an interesting general approach: the development of an algorithm that allows to predict sensitivity of glioblastoma cell subpopulations (i.e., NPC-, OPC-, MES- and AC-like cells) to certain drugs and drug classes. Tumor heterogeneity is a formidable challenge in the treatment of glioblastoma, and a better idea of which transcriptional tumor cell subpopulations repond to which drugs/classes can help to develop tailored multi- or oligo-drug strategies; potentially even identifying one drug that is able to target all, or at least the most important, tumor cell subpopulations. Having said that, this manuscripts falls short of providing evidence that the computational in silico approach presented here is of any practical (biological / therapeutic / mechanistic) value.\nThe investigation of one single drug in one single cell line in vivo demonstrating one major change - a shift towards MES-like tumor cell subpopulation - is no evidence for the utility of this entire approach. A shift towards MES-like states is well known for many therapeutic interventions now, and far away from any proof of a specific alteration induced.\nWe thank the reviewer for their important comment. We agree that as this tool was developed, the understanding of a common proneural-to-mesenchymal transcriptional response in glioma has been increasingly characterized. In line with comments from another reviewer, we have addressed this issue by identifying different small molecules which are predicted to have varying responses to that predicted and observed in response to alisertib. Importantly, not all proneural to mesenchymal programs are the same, and differing NPC-like and MES-like states have already been characterized. While not within the scope of this paper, more mechanistic studies on the proneural to mesenchymal response in GBM need to be performed.\nHowever, to address this comment, we have dramatically shifted the focus of the paper more towards our platform’s ability to effectively identify more effective combinations by identifying compounds that target distinct GBM cell transcriptional states present within the same tumor. Additionally, we have included an additional analysis demonstrating the predictive value of scFOCAL cell-drug connectivity utilizing different small molecules.\nFirst, to highlight the utility of scFOCAL in predicting heterogeneous cell population response outside of a proneural to mesenchymal shift, we demonstrate that cell-drug connectivity is predictive in the context of the HDAC inhibitor panobinostat using publicly available data of panobinostat-treated patient GBM acute slice cultures. Here, we also demonstrate the utility of scFOCAL to predict dynamics of the tumor-associated macrophage population within these same tumors, as described by the authors of the data’s originating manuscript (Zhao et al, 2021). Panobinostat TCS connectivity very effectively separates DMSO and panobinostat-treated tumor cells and macrophages on a probabilistic distribution (Supp. Fig 7de).\nSupplementary Figure S7: scFOCAL connectivity analysis with an L1000-derived panobinostat TCS separates vehicle-treated from panobinostat-treated cells in both neoplastic and myeloid cell populations. a. UMAP of 62,250 single-cell transcriptomes from patient-derived acute slice culture samples obtained from Zhao et al., 2021, treated with DMSO or the HDAC inhibitor panobinostat (0.2 μM). b. UMAP of single-cell transcriptomes colored by discrete cell type. c. Heatmap showing z-scores of panobinostat-specific gene expression (vs. treatment naïve) for genes retained in the Panobinostat TCS across cell lines tested in the L1000 dataset. d. Violin plots showing panobinostat TCS connectivity of Panobinostat-treated vs. DMSO-treated cells in neoplastic and myeloid populations, respectively (Wilcoxon rank sum test with continuity correction, p-value < 2.2e-16). e. Area-normalized kernel-density estimate (KDE) plots of scFOCAL calculated panobinostat TCS connectivity colored by treatment with DMSO or panobinostat, across neoplastic and myeloid cell populations, respectively.\nNext, we have developed a combination scoring approach within our platform, and demonstrate its ability to predict the relative synergy of small molecule combinations using publicly available synergy screen data obtained from Houweling et al. By effectively performing in silico perturbation studies using our patient scRNAseq dataset and performing differential connectivity analysis between predicted reference-drug-sensitive and -resistant populations, we derive a combination score that positively correlates with BLISS synergy in vitro in a large panel of GBM cell cultures.\nSupplementary Figure S8: scFOCAL combination scoring predicts relative synergy of small molecule combinations in vitro. a. Stacked bar plots depicting the proportion of cells predicted to be sensitive or resistant across individual patient GBM tumors to reference small molecules erlotinib, lapatinib, pazopanib, and sunitinib, respectively. b. Heatmap of BLISS raw sum synergy of small molecule combinations obtained from Houweling et al., 2023. The top bar plot annotation depicts scFOCAL calculated combination score using labeled reference molecules depicted on the bottom annotation bar. Sidebar annotations depict Spearman’s ρ and p-value between scFOCAL combination score and observed BLISS raw sum synergy for each cell line. c. Scatterplot of scFOCAL calculated combination score and mean BLISS raw sum synergy across all cell lines tested for small molecule combinations screened in Houweling et al., 2023 (Spearman's ρ = 0.7, p-value = 0.01). Points are colored by the reference molecule used for combination scoring. The shaded gray area around the line of best fit depicts the standard error.\nWe wanted to demonstrate that our approach is useful for predicting drugs targeting different cell populations in GBM tumors. For this reason, we included in vivo studies with the Olig 2 inhibitor CT179. CT-179 is predicted to target OPC-like cells and yield an AC-like resistant cell population. Through the prioritization of EGFR and tubulin targeting molecules by our combination prioritization strategy, which we previously included in our original submission as Supplementary Figure S6, we validate our predictions from a CT-179 response signature by demonstrating synergy in vivo with CT-179 and DepatuxM, an anti-EGFR MMAF (tubulin inhibitor) antibody-drug conjugate (Fig. 4, Fig 5, Fig S9). Importantly, our alisertib experiment still fits the altered narrative of this manuscript through demonstration that in silico perturbation of our patient GBM atlas accurately predicts differential connectivities to other drugs between the vehicle and alisertib treated PDX tumors (Fig 3k above, Supplementary Fig. S6), which led us to investigate the potential for our framework to prioritize synergistic combinations in our original submission.\nSupplementary Figure S6: In silico perturbation and predictive differential connectivity analysis models the drug-connectivity perturbation response to alisertib in vivo. a-h. Scatterplots of L1000 small molecule TCSs within individual compound classes depicting scFOCAL-predicted correlation shift (log FC) vs. observed correlation shift (log FC) in alisertib-treated xenografts. Differential small molecule correlations were calculated using limma.\nThey need to study the response of all 4 subpopulations to a panel of drugs experimentally; subpopulations can be studied by enrichment/depletion analyses, better individually after FACS sorting or, ideally, and something this reviewer would strongly recommend, using in vivo reporter systems or bar code systems that allow to distinguish all 4 subpopulations. This should be done in state-of-the-art in vitro assays.\nWe thank the reviewer for these suggestions. We have steered the focus of our manuscript towards our platform’s ability to effectively rank more effective combinations. To do this, we have focused on analyzing existing datasets that utilize L1000-tested compounds, which supported our novel platform scFOCAL, and performed new experiments using the OLIG2 inhibitor CT-179, generating a CT-179 transcriptional response signature, which we used to demonstrate in vivo that scFOCAL can identify novel synergistic combinations.\nAfter that, the authors need to demonstrate in vitro and in vivo that the most promising compound (mix) is indeed providing a therapeutic advantage compared to other approaches.\nWe have mined previously published data for panobinostat in slice cultures, and we demonstrate that the scFOCAL predictions are consistent with what was observed by single-cell analysis ex vivo. These studies further support that our pipeline predicts the cells that are targeted pharmacologically in GBM, and that drug-cell connectivity is a predictive value of how cells respond to different treatments. As suggested, we have also extended our findings to develop a scFOCAL combination index to score combinations for potential synergy. Using mined synergy screen data from Houwelling et al, 2023, we score small molecule combinations using scFOCAL, and found that our predictions correlate with relative synergy results across 24 GBM cell lines in spheroid culture. Having validated our combination index using in vitro data, we moved to identify a novel combination using the OLIG2 inhibitor CT-179, and demonstrated the efficacy of our prioritized combination in vivo.\nMore minor point: the presentation of the results and figure legends is often confusing and should be significantly improved. All figure panels should be cited in the order of appearance; it should be explained in the results section what they actually show; and the reader needs to get more explanatory information there, too.\nWe thank the reviewer for these suggestions and have now modified the order of the figures to be consistent with the order of appearance, in addition to the major overhaul of all of the figures included. We have also explained in depth what the results show, as suggested in the figure legends and results."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank the authors for answering my questions, and also for the rich additional analysis. I think the synergy prediction strongly increased the importance of this paper.", "reviewer_2": "This is a revised manuscript describing a new computational tool for predicting drug sensitivity at the single cell level. The main point raised in the prior round of review was the robustness of the prediction, as limited experimental data were shown, and the question about the ability to predict combinations of drugs to target the subpopulations.\nThe presented approach to identify drug combinations is highly warranted. However, the scFOCAL wouldn’t have been needed to come up with a combination of OLIG2i and EGFR-MMAF, since both targets are what distinguishes the transcriptional cell states in GBM. Docetaxel and varlitinb don’t seem to be the top hits in the connectivity analysis. It would be important for the readers to understand how to prioritize hits that may be selected when they aply scFOCAL to their own dataset.\nSince this manuscript has been significantly altered from the previous version, there are also several other points that should be addressed:\n1. CytoTRACE should be better described in the text, as it is unclear what in represents in Fig. 2f&g.\n2. Line 168: “Within each individual, cycling cells seem to separate into a unique drug connectivity state”. Is this true for all patients in this cohort or only for the two tumors presented in Fig. 2?\n3. The authors focus on alisertib, based on previous connection to MES cell state. There are a number of other candidates, including alpelisib and afatinib, which have highly differential connectivity between cell states. This should be clarified.\n4. Fig.3a. Predicted sensitive vs resistant populations are divided right at 0 point for TCS correlation coefficient. Would the differences be more striking if top and bottom quartiles were used? With the current split, in Fig.3b the distribution looks very similar.\n5. Fig. 3j. Is the increase in MES and decrease in NPC after alpelisib treatment statistically significant? The changes are not massive, so statistical test would be helpful.\n6. Fig. 3k. The correlation between observed and predicted shifts is weak overall, yet for some classes of compounds they work well (Fig. S6). It might be worth to put S6 panels in the main figure, since they represent a more meaningful result. Is there any connection between these drug classes and action or transcriptional consequences of alisertib treatment?\n7. The Panobinostat response analysis seems a bit disconnected from the rest of the manuscript. Since this is an scRNAseq dataset, are the authors able to see changes in Neftel et al. cell states in the cultured slices?\n8. Docetaxel and varlitinb don’t seem to come up as top hits in the connectivity analysis. Could the authors comment on the top hits and why they were not pursued for experimental validation?\n9. Fig. 2b – scale for disease signature reversal is missing", "reviewer_3": "The authors have responded well to many of my comments and have now provided a better proof-of-principle that their novel cell state - based approach can indeed help to identify new therapeutic combinations in GBM."}, "author_response": "We thank the reviewers for their positive assessment of our revised manuscript. We have now addressed Reviewer 2’s concerns point-by-point in a second revision. We hope our manuscript is now acceptable for publication at Nature Communications.\nREVIEWER COMMENTS:\n● Reviewer #1:\n○ “Thank the authors for answering my questions, and also for the rich additional analysis. I think the synergy prediction strongly increased the importance of this paper.”\n● Reviewer #3:\n○ “The authors have responded well to many of my comments and have now provided a better proof-of-principle that their novel cell state-based approach can indeed help to identify new therapeutic combinations.”\n● Reviewer #2 (Comments Addressed)\n○ Comment #1: “This is a revised manuscript describing a new computational tool for predicting drug sensitivity at the single-cell level. The main point raised in the prior round of review was the robustness of the prediction, as limited experimental data were shown, and the question about the ability to predict combinations of drugs to target the subpopulations. The presented approach to identify drug combinations is highly warranted. However, the scFOCAL wouldn’t have been needed to come up with a combination of OLIG2i and EGFR-MMAF, since both targets are what distinguish the transcriptional cell states in GBM. Docetaxel and varlitinb don’t seem to be the top hits in the connectivity analysis. It would be important for the readers to understand how to prioritize hits that may be selected when they apply scFOCAL to their own dataset.”\n○ Comment #1 Response: We thank the reviewer for their positive assessment of the need for our platform to predict combinations in an unbiased manner at the single-cell level. Regarding whether scFOCAL would have been needed to come up with a combination of OLIG2i and ɑ-EGFR-MMAF, we would like to point out that prior to our findings, no other group has predicted and demonstrated that the combination of an OLIG2i and ɑ-EGFR-MMAF produces a survival benefit in preclinical models of GBM. Our work does exactly that and puts forth a novel clinical candidate, OLIG2i CT-179, as a possible combination therapy in GBM. These findings are essential to clinical trial designs that will emanate from this work. Therefore, our findings are significant and important for the medical and scientific GBM community.\nA neuro-oncology expert may indeed identify that this combination makes sense in the context of GBM cell transcriptional states. Uniquely, scFOCAL identifies that OLIG2 inhibition may sensitize cells to drugs acting on tubulin, simultaneously identifying that these sensitized cells likely also express EGFR in a single quantitative measure, our combination score. The results of our in silico combination screening, taken together in context with our identification of EGFR-expressing AC-like GBM cells as the most concordant with our CT-179 response signature, motivated us to validate the combination of CT-179 and Depatux-M in vivo over exploring alternative targets.\nImportantly, Depatux-M does not necessarily inhibit EGFR signaling, but utilizes EGFR expression to guide a tubulin poison to the target cells. While varlitinib, an EGFR inhibitor, scores highly for our combination index, this is mainly contributed to by its corresponding logFC differential in connectivity, suggesting that the connectivity between CT-179 sensitive and resistant cells differs, and that varlitinib is predicted to target different cells than CT-179. By contrast, indibulin’s high combination index is mainly contributed to by the magnitude of its mean discordance with CT-179-resistant cells. Combined with the interpretation of the scFOCAL prediction of AC-like cell state resistance, Depatux-M, leveraging this approach to kill the cell, while homing to EGFR expressing cells (regardless of EGFR activation), makes sense in the context of the individual factors comprising the scFOCAL combination index.\nIn reference to the observation that docetaxel and varlitinib are not the top hits in our connectivity analysis, we would like to point out that they are all within in the top 31 candidates from 1674 molecules we initially began with (Top 1.8%), from the 1382 annotated as clinically tested by the repurposing hub (Top 2.25%), and from the 411 passing filter for being sensitized by CT-179 in silico and for being discordant with CT-179-predicted-resistant cells (Top 7.5%). Depatux-M is a more attractive therapeutic than either alone due to its status as an antibody drug conjugate targeting EGFR. Further, Depatux-M has been utilized in a Phase III clinical trial for glioblastoma, which demonstrated extension of progression-free survival (PMID: 35849035), indicating its safety and BBB penetration. Lastly, our alisertib studies demonstrated the accuracy in predictions of tubulin polymerization inhibitors, lending to further confidence in prioritizing indibulin and docetaxel.\nRegarding how an individual reader can prioritize hits that may be selected when they apply scFOCAL to their own dataset, we provide a combination score that can rank various compounds. This combination score must be considered in the context of their disease of interest and in the biology underlying heterogeneity within that specific disease. Further studies will be required to determine the possible cut-offs for scores to predict combinations that are validated through in vitro and in vivo combination studies.\nWe now add to the limitations of the study this important point, that although our combination index provides a predicted measure of synergy of compounds, as a function of how broadly the combination targets a heterogeneous tumor cell landscape, the identification of effective combinations further requires an iterative process, where data and information from the scientific literature must be considered.\nImportantly, our manuscript demonstrates the presence and utility of a biologically relevant cell-and-drug transcriptional connectivity signal, of which further developments can be applied to further strengthen its application.\n○ Comment #2: “CytoTRACE should be better described in the text, as it is unclear what it represents in Fig. 2f&g.”\n○ Comment #2 Response: In the discussion section, we have further explained our interpretation of drug connectivity clustering of cells exhibiting patterns with tumor cell CytoTRACE scores, elaborating on the interpretation that cells clustering on drug connectivity correlating with different relative CytoTRACE scores indicates that tumor cells in different differentiation states exhibit unique sensitivities from more differentiated cell states. Cells with a higher CytoTRACE score have a higher transcriptional diversity, and thus represent more stem-like de-differentiated tumor cell states. In our new Supplementary Fig. S3b, it can be better seen that these CytoTRACE scores are correlated with NPC-like GBM cell transcriptional states.\n○ Comment #3: “Line 168: 'Within each individual, cycling cells seem to separate into a unique drug connectivity state”. Is this true for all patients in this cohort or only for the two tumors presented in Fig. 2?”\n○ Comment #3 Response: We thank the reviewer for this important point. In response, we have added a new supplementary figure (Supplementary Fig. S3). This figure first presents the individual clustering heatmaps for each of the 11 patients. We then performed hierarchical clustering on the resulting cell groups using three key features: 1) proportions of Neftel GBM states, 2) cell cycle phase proportions, and 3) the mean CytoTRACE score. This analysis reveals that the patient clusters organize into distinct, high-level metaclusters. We have added further comments on this within our revised manuscript and have removed the statement in question. While a proliferative meta-cluster is apparent, there are also variations in cycling proportion within metaclusters.\nSupplementary Figure S3: Cell-drug connectivity using FDA-approved oncology drugs clusters patient GBM tumor cells reflecting patterns of GBM cell transcriptional state. a) Correlation matrices depicting pairwise Spearman correlations of single GBM tumor cells from each patient’s tumor cells. Column annotations depict each cell’s CytoTRACE score, cell cycle phase, and dominant Neftel et al. transcriptional state. b) Heatmap depicting hierarchical clustering of all clusters k from all 11 patients on the relative proportions of cells within each GBM cell transcriptional state, cell cycle phase, and mean CytoTRACE score.\n○ Comment #4: “The authors focus on alisertib, based on previous connection to MES cell state. There are a number of other candidates, including alpelisib and afatinib, which have highly differential connectivity between cell states. This should be clarified.”\n○ Comment #4 Response: We thank the reviewer for this important request, as other small molecules may be of interest to other researchers. Our prior unpublished work shows that alisertib monotherapy elicits a significant increase in survival time in orthotopic xenografts. This is consistent with what was previously reported about alisertib (PMID:27816996). Thus, we were interested in characterizing the eventual development of resistance and tumor growth despite treatment with alisertib. Further, our prior publications indicate that alisertib is well tolerated in vivo and does get into the brain. While we agree that alpelisib and afatinib are of interest, we do not have in vivo data to support their efficacy or brain penetrance in our hands. We can pursue these studies in the future and comment in the Discussion that similar molecules should be analyzed further.\nRebuttal Figure 1: Kaplan-Meier survival curve for GBM22-bearing mice treated with DMSO or alisertib. (unpublished)\n○ Comment #5: Fig.3a. Predicted sensitive vs resistant populations are divided right at 0 point for TCS correlation coefficient. Would the differences be more striking if top and bottom quartiles were used? With the current split, in Fig.3b the distribution looks very similar.\n○ Comment #5 Response: We thank the reviewer for this important question, and have added additional analysis to Supplementary Figure S4 to address this. To demonstrate that the scFOCAL prediction for relative sensitivity or resistance to alisertib in the context of Neftel et al. states is robust, we have iteratively re-created Figure 3c across different quantiles of alisertib connectivity (i.e., top 5% vs bottom 5%, top 10% vs bottom 10%, through to the top 45% vs the bottom 45%. Within each patient, the proportion change of each cell state was compared using Wilcoxon rank-sum tests and Holm multiple comparisons correction. At each different quantile comparison, the differences in proportion are significant between MES-like and NPC-like populations. At the 5% comparison, MES and OPC-like proportion shifts are also statistically significant.\nSupplementary Figure S4: Interrogation of an alisertib TCS using single-cell RNA sequencing indicates the predicted sensitivity of an NPC-like GBM cell state. a. Diagram describing interrogation of the L1000-derived alisertib TCS using GBM cell state markers as in (b) or by discordance ratio with GBM cell state disease signatures as in (c). b. Heatmap depicting differential expression of cell state markers that overlap with alisertib TCS genes. Cell color indicates the relative log2FC of each gene, and black heatmap cells indicate a gene was not differentially expressed by a given cell state. Row annotation depicts the magnitude and direction of gene expression changes of the alisertib TCS. c. Box plot of cell state disease signature discordance to the alisertib TCS by individual patient tumor. Cell state disease signatures were calculated using MAST to compare cells within each GBM cell state to non-neoplastic cell types within each individual. d. Bar plot depicting the pairwise post-hoc -log(p.adj) values of proportion change for each state combination as calculated from main Figure 3c. e-m. Histograms of patient GBM tumor cells colored by assigned alisertib TCS quantile-based sensitivity or resistance, paired with corresponding box plots of predicted proportion change for each Neftel et al. GBM cell transcriptional state. For each quantile cut-off, the change in proportion was calculated by subtracting the proportion in the sensitive group from that of the resistant group. Post-hoc pairwise comparisons of relative state shifts were performed with the rstatix package using Wilcoxon rank-sum tests, and the resulting p values were adjusted for multiple comparisons using the Holm method (*p.adj < 0.05, **p.adj <0.005).\n○ Comment #6: Fig. 3j. Is the increase in MES and decrease in NPC after alpelisib treatment statistically significant? The changes are not massive, so statistical test would be helpful.\n○ Comment #6 Response: We thank the reviewer for this important comment. The Neftel et al. state identities of the tumor cells are assigned based on the state most enriched in that cell by singscore enrichment using the original 50 gene signatures. Because for this experiment, individual tumors from 3 mice per treatment were pooled into single GEM captures (10X), we are unable to perform statistical tests on proportions of cells. However, in Figure 3i, depicting the normalized shift in each state’s enrichment utilizing individual cells as datapoints, we are able to perform statistics and demonstrate that the NPC-like signature depletion and MES-like signature enrichment are statistically significant, while the change for OPC-like enrichment does not reach significance. Panel 3j is another perspective of 3i, depicting changes in the predominantly expressed signatures of tumor cells, calculated from the significantly changed enrichment values depicted in Fig. 3i.\nFig. 3i. Bar plot of mean enrichment shift of transcriptional state signatures in alisertib-treated xenograft cells normalized to DMSO vehicle-treated xenograft cells. Error bars represent 95% confidence interval (Wilcoxon with Benjamini-Hochberg correction: *p-adjusted = 8e-9; **p-adjusted = 1.1e-60;***p-adjusted = 3.2e-238; NS p-adjusted = 0.79). Fig. 3j. Alluvial plot depicting shift in relative proportion of transcriptional state identities in alisertib and DMSO vehicle control-treated xenografts.\n○ Comment #7: Fig. 3k. The correlation between observed and predicted shifts is weak overall, yet for some classes of compounds they work well (Fig. S6). It might be worth to put S6 panels in the main figure, since they represent a more meaningful result. Is there any connection between these drug classes and action or transcriptional consequences of alisertib treatment?\n○ Comment #7 Response: We thank the reviewer for this important suggestion. We have now moved the panels for topoisomerase inhibitors, NRT (nucleoside reverse transcriptase), and serotonin receptor agonists from Supplementary Fig. S6 to the main figure, Fig. 3k, as suggested. We chose to further highlight these compound classes to demonstrate the predictive power to predict compounds with decreased cell discordance following alisertib treatment, such as the majority of topoisomerase inhibitors, but also to show the power to predict compounds with increased discordance against tumor cells, such as NRT and select serotonin receptor agonists. In the manuscript, we have added a discussion of the relationship between topoisomerase inhibition and reduced expression of Aurora Kinase A, the target of alisertib. Further, we have added a comment on the potential for future investigation of NRT and serotonin receptor agonists, as we have previously published on the potential of repurposing psychiatric drugs for brain cancers, particularly due to their known blood-brain barrier penetrance.\nFigure 3: In silico perturbation of GBM tumor cell scRNAseq data using an L1000-derived alisertib TCS predicts an NPC-like to MES-like tumor response confirmed in vivo. a. Histogram of single-cell alisertib TCS correlations (connectivities). Cells with negative correlation coefficients (blue) are predicted to be sensitive to alisertib, while cells with positive correlation coefficients (red) are predicted to be resistant to alisertib. b. Hierarchy plot of patient GBM tumor cells colored by predicted sensitivity (ρ < 0) or resistance (ρ > 0) to alisertib. c. Bar plot depicting mean shift in proportions of cells in resistant vs. sensitive populations within individual patient tumors. Error bars depict the standard error of the mean percentage differences of individual patient tumors. d. Schematic of in vivo experiments. e. UMAP plot of pre-filter single-cell transcriptomes colored by treatment with either alisertib or DMSO vehicle control. f. UMAP of captured cells from orthotopic xenografts colored by percent alignment to the human transcriptome (hg19). g. Dot plot of pass-filter, human xenograft single-cell transcriptome expression of Neftel et al. signatures, grouped by assigned transcriptional state identity, as assigned based on predominant transcriptional state module expression. h. Two-dimensional hierarchical representation of GBM22 xenograft cells’ relative enrichment scores for GBM cell transcriptional state modules. Cells are colored by assigned transcriptional state identity. i. Bar plot of mean enrichment shift of transcriptional state signatures in alisertib-treated xenograft cells normalized to DMSO vehicle-treated xenograft cells. Error bars represent 95% confidence interval (Wilcoxon with Benjamini-Hochberg correction: *p-adjusted = 8e-9; **p-adjusted = 1.1e-60;***p-adjusted = 3.2e-238; NS p-adjusted = 0.79). j. Alluvial plot depicting shift in relative proportion of transcriptional state identities in alisertib and DMSO vehicle control-treated xenografts. k. Scatterplots of L1000 small molecules and select subsets of compound classes’ TCSs depicting predicted correlation shift (log2FC) vs. observed correlation shift (log2FC) in alisertib-treated xenografts. Differential small molecule correlations were calculated using limma. R (ρ) values determined using the Pearson correlation.\n○ Comment #8: The Panobinostat response analysis seems a bit disconnected from the rest of the manuscript. Since this is an scRNAseq dataset, are the authors able to see changes in Neftel et al. cell states in the cultured slices?\n○ Comment #8 Response: We thank the reviewer for this important request to increase the relevance of the meta-analysis of the acute slice culture treatment data included from Zhao et al. To do this, we have further analyzed the dataset in the context of Neftel et al. transcriptional states, and show that the pattern of shifts in transcriptional states induced by panobinostat treatment is also correlated with increases in panobinostat TCS connectivity. To do this, we have split the data on the panobinostat TCS connectivity quartiles, and analyzed the shift in proportions of cells predominantly within each of the 4 transcriptional states. (Supplementary Fig. S4). We directly compare these predicted shifts to those observed in the actual treatment groups across all 5 individual patient sources. While the proportion shifts between treatment groups do not reach statistical significance, a clear trend is apparent, which panobinostat TCS connectivity also identifies and directly correlates with. Further, the inclusion of this dataset and meta-analysis within the supplementary material is warranted due to the perceived targeting of AC-like cell states by panobinostat, unique from both alisertib and CT-179, which were found to be most transcriptionally discordant with NPC-like and OPC-like cells, respectively.\nSupplementary Figure S7: scFOCAL connectivity analysis with an L1000-derived panobinostat TCS separates vehicle-treated from panobinostat-treated cells in both neoplastic and myeloid cell populations. a. UMAP of 62,250 single-cell transcriptomes from patient-derived acute slice culture samples obtained from Zhao et al., 2021, treated with DMSO or the HDAC inhibitor panobinostat (0.2 μM). b. UMAP of single-cell transcriptomes colored by discrete cell type. c. Heatmap showing z-scores of panobinostat-specific gene expression (vs. treatment naïve) for genes retained in the Panobinostat TCS across cell lines tested in the L1000 dataset. d. Violin plots showing panobinostat TCS connectivity of Panobinostat-treated vs. DMSO-treated cells in neoplastic and myeloid populations, respectively (Wilcoxon rank sum test with continuity correction, p-value < 2.2e-16). e. Area-normalized kernel-density estimate (KDE) plots of scFOCAL calculated panobinostat TCS connectivity colored by treatment with DMSO or panobinostat, across neoplastic and myeloid cell populations, respectively. f. Barplots of cell proportion for each GBM cell transcriptional state, with cells split based on panobinostat TCS connectivity quartiles (high vs low). g. Barplots of cell proportion for each GBM cell transcriptional state, comparing DMSO vehicle and panobinostat-treated tumors h. Scatterplot of proportion shifts predicted from Q1-Q3 comparisons, plotted against proportion shift between treatment groups. R value determined using Spearman correlation.\n●\n○ Comment #9: “Docetaxel and varlitinb don’t seem to come up as top hits in the connectivity analysis. Could the authors comment on the top hits and why they were not pursued for experimental validation?”\n○ Comment #9 Response: We thank the reviewer for this important comment. In response, we would first like to clarify that these 3 compounds highlighted behind the interest of Depatux-M, namely indibulin, docetaxel, and varlitinib, are prioritized within the top 7% of the 411 molecules passing both the negative-logFC and negative mean resistant cell connectivity filters (Indibulin: 16/411 (top 3.9%); Docetaxel: 21/411 (top 5.1%); Varlitinib: 31/411 (top 7.6%)). Prior to this filtering, these 411 compounds are drawn from a pool of 1382 small molecules selected for use in clinical trials as annotated in the repurposing hub (Indibulin: 16/1382 (Top 1.2%); Docetaxel 21/1382 (Top 1.5%); Varlitinib 31/1382 (Top 2.3%)).\nImportantly, the coloring in Figure 5c depicts the scFOCAL combination index, which we leverage for prioritizing. While other top hits are definitely interesting, we sought to prioritize a combination that would demonstrate that targeting distinct GBM cell states simultaneously improves survival. Depatux-M, harboring qualities of multiple scFOCAL prioritized hits, has been utilized in phase III trials, is safe, and extended PFS in glioblastoma (PMID 35849035). To this end, we have added a further detailed discussion of our focus on Depatux-M. Lastly, we have made sure to include a table of CT-179-resistant vs sensitive logFC, mean CT-179-resistant cell connectivity (MRC), and combination index values in the supplementary data, within Supplementary Data 6.\n○ Comment #10: Fig. 2b – scale for disease signature reversal is missing\n○ Comment #10 Response: We thank the reviewer for catching this important omission. We have since added the scale for each individual Neftel et al. state discordance ratio for the annotation of the heatmap in Fig. 2b.\nUpdated Fig 2b. Correlation matrix depicting similarities of L1000 small molecule TCSs by their connectivity to all individual cells within our single-cell atlas. Row annotations depict compound TCS discordance ratios for the reversal of AC-, MES-, NPC-, and OPC-like disease signatures calculated against non-neoplastic cells in the dataset."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors clarified all of the points raised in review. The manuscript will be of high interest to Nature Communications readers hoping to streamline synergistic drug combination testing."}, "author_response": ""}]

Decision letter