DRUG AND SINGLE-CELL GENE EXPRESSION INTEGRATION IDENTIFIES SENSITIVE AND RESISTANT GLIOBLASTOMA CELL POPULATIONS

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. GBM tumors contain neoplastic cells within a spectrum of different transcriptional states. Identifying effective therapeutics requires a platform that predicts the differential sensitivity and resistance of these states to various treatments. Here, we developed a novel framework, ISOSCELES (Inferred cell Sensitivity Operating on the integration of Single-Cell Expression and L1000 Expression Signatures), to quantify the cellular drug sensitivity and resistance landscape. Using single-cell RNA sequencing of newly diagnosed and recurrent GBM tumors, we identified compounds from the LINCS L1000 database with transcriptional response signatures selectively discordant with distinct GBM cell states. We validated the significance of these findings in vitro, ex vivo, and in vivo, and identified a novel combination of an OLIG2 inhibitor and Depatux-M for GBM. Our studies suggest that ISOSCELES identifies cell states sensitive and resistant to targeted therapies in GBM and that it can be applied to identify new synergistic combinations.

Keywords: Drug Resistance; Glioblastoma; Pharmacology; Pharmacotranscriptomics; Single-Cell RNA Sequencing.

Figure 1:. Single-cell RNA sequencing reveals distinct transcriptional states present in newly diagnosed and recurrent GBM.

Single-cell RNA sequencing data of 6 patient glioblastoma tumors was integrated and harmonized with an external dataset obtained from Johnson et al. (2021) a-e. UMAPs of single-cell transcriptomes colored by (a) cell type, (b) source dataset, (c) whether the tumor was newly diagnosed or recurrent, (d) CytoTRACE score, and (e) assigned Neftel et al. (2019) GBM cell transcriptional state. f. Two-dimensional representation of relative enrichment of GBM cell transcriptional states in neoplastic cells. Cells are colored by assigned transcriptional state identity based on the most predominantly enriched signature of that cell. g. Sankey plot depicting proportions of cells grouped by source dataset, occurrence or recurrence, tumor ID, cell type, GBM cell transcriptional state, and expression-based cell cycle phase.

Figure 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 inverse 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).

Figure 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 (ANOVA with Games-Howell Post-hoc, AC-MES p.adj = 3.27e-08, AC-NPC p.adj < 1e-13, AC-OPC p.adj = 3.00e-13, MES-NPC p.adj = 2.88e-08, MES-OPC p.adj = 3.04e-08, NPC-OPC p.adj < 1e-13). j. Alluvial plot depicting shift in relative proportion of transcriptional state identities in alisertib and DMSO vehicle control treated xenografts. h. Bar plot of change in proportion of cell transcriptional identities in alisertib-treated xenografts relative to DMSO vehicle-treated xenografts. k. Scatterplot of L1000 small molecule TCS’s depicting predicted correlation shift (log2FC) vs. observed correlation shift (log2FC) in alisertib-treated xenografts. Differential small molecule correlations were calculated using limma.

Figure 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 ISOSCELES framework. 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 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.

Figure 5:. An ISOSCELES combination index predicts 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 connectivityand mean CT-179 resistant cell connectivity. Colors depict the calculated ISOSCELES combination index, the product of multiplying the differential connectivity log2FC values by mean CT-179 resistant cell connectivity values for each individual 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 Log-rank (Mantel-Cox) test.

Figure 6:

The ISOSCELES framework is available for use as a shiny web application or an R package.