Deconvolution of cell type-specific drug responses in human tumor tissue with single-cell RNA-seq.

BACKGROUND: Preclinical studies require models that recapitulate the cellular diversity of human tumors and provide insight into the drug sensitivities of specific cellular populations. The ideal platform would enable rapid screening of cell type-specific drug sensitivities directly in patient tumor tissue and reveal strategies to overcome intratumoral heterogeneity.
METHODS: We combine multiplexed drug perturbation in acute slice culture from freshly resected tumors with single-cell RNA sequencing (scRNA-seq) to profile transcriptome-wide drug responses in individual patients. We applied this approach to drug perturbations on slices derived from six glioblastoma (GBM) resections to identify conserved drug responses and to one additional GBM resection to identify patient-specific responses.
RESULTS: We used scRNA-seq to demonstrate that acute slice cultures recapitulate the cellular and molecular features of the originating tumor tissue and the feasibility of drug screening from an individual tumor. Detailed investigation of etoposide, a topoisomerase poison, and the histone deacetylase (HDAC) inhibitor panobinostat in acute slice cultures revealed cell type-specific responses across multiple patients. Etoposide has a conserved impact on proliferating tumor cells, while panobinostat treatment affects both tumor and non-tumor populations, including unexpected effects on the immune microenvironment.
CONCLUSIONS: Acute slice cultures recapitulate the major cellular and molecular features of GBM at the single-cell level. In combination with scRNA-seq, this approach enables cell type-specific analysis of sensitivity to multiple drugs in individual tumors. We anticipate that this approach will facilitate pre-clinical studies that identify effective therapies for solid tumors.

Background

Inter- and intra-tumoral heterogeneity present major challenges for cancer therapy. While scRNA-seq can determine the cellular composition of complex tumors and even reveal cell type-specific drug sensitivities, these measurements are ultimately limited by models of drug response. Acute slice cultures are an attractive approach to modeling drug response in solid tumors because multiple cultures can be rapidly generated from a single surgical specimen, and they do not require extensive culturing or manipulation, which leads to distortion of the native composition of the tissue, selection, and loss of heterogeneity by diluting populations that do not proliferate rapidly [1-3]. Furthermore, drug perturbation experiments in acute slice cultures can be carried out rapidly, beginning on the day of surgical resection, on timescales relevant for clinical decision-making. GBM is an ideal setting for testing this approach because it exhibits profound inter- and intratumoral heterogeneity and is the most common and deadly primary brain malignancy in adults. Surgical resection is part of the standard-of-care, and robust protocols for acute slice culture of human GBM have been established previously [4-6]. Furthermore, GBM has been extensively characterized by scRNA-seq, providing a detailed baseline for both the transformed populations that co-occur in individual patients and the microenvironment [7-11]. Indeed, single-cell characterization of GBM models has highlighted the importance of the tumor microenvironment in maintaining the phenotypic diversity of malignant cells [12]. We obtained GBM surgical specimens from seven patients, generated multiple 500 micron slices, and placed them in short-term culture for drug perturbation (Fig. 1a). Screens were completed within 24 h of surgery and analyzed immediately by scRNA-seq using our scalable microwell platform [13, 14] to deconvolve cell type-specific responses to multiple drugs (Fig. 1a).

Fig. 1

a Schematic illustration of experimental and analytical methods for slice culture drug perturbation and scRNA-seq. b UMAP embedding of scRNA-seq profiles from acutely isolated biopsies and slice cultures from different regions of the same tumor (PW032) colored by sample origin. c Same as b but colored by the log-ratio of Chr. 7 to Chr. 10 average expression where a high ratio (red) indicates malignant transformation. d Same as b but colored by cell type. e Heatmap of average expression of marker genes from cell types in the tumor microenvironment in each cell type and sample from PW032. f Fractional abundance of each major cell type in each biopsy and slice culture sample from PW032. g Two-dimensional model projecting each transformed cell from PW032 biopsies and slice into four major GBM transformed populations colored by sample origin

Methods

Preparation and culture of tissue slices

This work was approved by the Columbia University Irving Medical Center Institutional Review Board before commencing the study. All tumor specimens were procured from surgeries at Columbia University Irving Medical Center. Patient diagnosis information can be found in Additional file 1: Table S1. Tumor specimens were collected immediately after surgical removal and kept in ice-cold artificial cerebrospinal fluid (ACSF) solution containing 210 mM sucrose, 10 mM glucose, 2.5 mM KCl, 1.25 mM NaH2PO4, 0.5 mM CaCl2, 7 mM MgCl2, and 26 mM NaHCO3 for transportation. Preparation of ex vivo tissue slices was modified from methods described previously [5]. Briefly, the collected tumor specimens were placed in a drop of ice-cold ACSF and sliced using a tissue chopper (McIlwain) at a thickness of 500 μm under sterile conditions. The slices were immediately transferred to the ice-cold ACSF solution in 6-well plates using a sterile plastic Pasteur pipette half filled with ice-cold ACSF solution followed by a 15-min recovery in ACSF to reach room temperature. Intact and well-shaped slices (approximately 5–10-mm diameter) were then placed on top of a porous membrane insert (0.4 μm, Millipore). Then the membrane inserts were placed into 6-well plates containing 1.5 mL maintenance medium consisting of F12/DMEM (Gibco) supplemented with N-2 Supplement (Gibco) and 1% antibiotic-antimycotic (ThermoFisher). To ensure proper diffusion into the slice, culture medium was placed under the membrane insert without bubbles. A drop of 10 μl of culture medium was added directly on top of each slice to prevent the slice surface from drying. The slices were first rested for 6 h with the maintenance medium in a humidified incubator at 37°C and 5% CO2. Then, the medium was replaced with pre-warmed medium containing drugs with desired concentration (Additional file 1: Table S2) or corresponding volume of vehicle (DMSO). Drug dose was chosen as the estimated IC20 as measured in TS543 patient-derived glioma neurospheres (Additional file 1: Fig. S1) [15], which we have previously shown to harbor both proneural and mesenchymal GBM subpopulations [14]. Slices were then cultured with the treatment medium in a humidified incubator at 37°C and 5% CO2 for 18 h before being collected for dissociation.

Dissociation of tissue and slices

Collected tissue samples or tissue slices were dissociated using the Adult Brain Dissociation kit (Miltenyi Biotec) on gentleMACS Octo Dissociator with Heaters (Miltenyi Biotec) according to the manufacturer’s instructions.

Microwell scRNA-seq

Dissociated cells from each slice were profiled using microwell-based single-cell RNA-seq [14] as previously described [9, 16] with the following modifications: once the RNA-capture step was finished, sealing oil was flushed out of the devices by pipetting 1 mL of wash buffer supplemented with 0.04 U/μl RNase inhibitor (Thermo Fisher Scientific) and then beads were extracted from the device and resuspended in 200 μl of reverse transcription mixture. Bead-suspensions were divided into 50-μl aliquots and placed into PCR tubes (Corning) followed by incubation at 25°C for 30 min and at 42°C for 90 min in a thermocycler. Each cDNA library was barcoded with an Illumina sample index. Libraries with unique Illumina sample indices were pooled for sequencing on (1) an Illumina NextSeq 500 with an 8-base index read, a 21-base read 1 containing cell-identifying barcodes (CB) and unique molecular identifiers (UMIs), and a 63-base read 2 containing the transcript sequence, or (2) an Illumina NovaSeq 6000 with an 8-base index read, a 26-base read 1 containing CB and UMI, and a 91-base read 2 containing the transcript sequence.

scRNA-seq data preprocessing

Raw data obtained from the Illumina NextSeq 500 was trimmed and aligned as described previously [9]. For each read with a unique, strand-specific alignment to exonic sequence, we constructed an address comprised of the CB, UMI barcode, and gene identifier. Raw data obtained from the Illumina NovaSeq 6000 was first corrected for index swapping to avoid cross-talk between sample index sequences using the algorithm described by Griffiths et al. [17] before assigning read addresses for each sample. For samples had been sequenced on both Illumina NextSeq 500 and NovaSeq 6000, we combined the addresses from the NextSeq 500 and the corrected addresses from the NovaSeq 6000 for data processing as described previously [9, 16]. Briefly, reads with the same CB, UMI, and aligned gene were collapsed and sequencing errors in the CB and UMI were corrected to generate a preliminary matrix of molecular counts for each cell.

We applied the EmptyDrops algorithm to recover true cell-identifying barcodes in the digital gene expression matrix [18]. We then removed CBs that satisfied any of the following criteria: (1) fractional alignment to the mitochondrial genome greater than 1.96 standard deviations above the mean; (2) a ratio of molecules aligning to whole gene bodies (including introns) to molecules aligning exclusively to exons greater than 1.96 standard deviations above the mean; (3) average number of reads per molecule or average number of molecules per gene >2.5 standard deviations above the mean for a given sample; or (4) more than 40% of UMI bases are T or where the average number of T-bases per UMI is at least 4.

One important consideration, particularly with the drug-treated slice cultures where we expect increased cell death, is whether there are elevated levels of ambient mRNA (i.e., higher background) in the scRNA-seq profiles. We examined this by comparing the coverage of cell-identifying barcodes assigned to cells to those assigned to ambient RNA using EmptyDrops as described above and found that drug-treatment did not result in significantly increased background (Additional File 1: Fig. S2).

Unsupervised clustering, differential expression, and visualization

Clustering, visualization, and identification of cluster-specific genes were performed as described previously (www.github.com/simslab/cluster_diffex2018) [16, 19]. We used Louvain community detection as implemented in Phenograph for unsupervised clustering with k=20 for all k-nearest neighbor graphs [20]. For all clustering and visualization analyses of merged datasets, we first identified marker genes using the drop-out curve method described in Levitin et al. [16] (www.github.com/simslab/cluster_diffex2018) [19] for each individual sample and took the union of the resulting marker sets to cluster and embed the merged dataset. We projected drug-treated cells onto vehicle-treated cells with UMAP in Fig. 3 as described in Szabo et al. and Levitin and Sims (code available at www.github.com/simslab/umap_projection) [21, 22].

Fig. 3

a Experimental schematic for slice culture drug screening (6 drugs, 2 controls) from a single patient (PW030). b Heatmap showing the number of differentially expressed genes (FDR<0.01) in the tumor, myeloid, and oligodendrocyte populations between treated and control slices for each drug in the screen illustrated in a. c Same as b but showing only differentially expressed genes with FDR<0.01 and fold-change amplitude greater than two (both up- and downregulated genes). d UMAP embedding of scRNA-seq profiles of transformed cells from the control slices colored by expression of two proliferation markers (TOP2A, MKI67), two mesenchymal markers (CD44, VIM), and an astrocyte marker (GFAP). e Same as d but with UMAP projection density of scRNA-seq profiles of transformed cell from the treated slice cultures for each drug. Note that there is negligible projection density for the etoposide-treated cells onto the control cells for the small proliferative population expressing TOP2A and MKI67

Whole genome sequencing and analysis

Genomic DNA was extracted from a piece of frozen tissue from each tumor using the DNeasy Blood & Tissue Kits (Qiagen) according to the manufacturer’s instructions and was submitted to the Beijing Genomics Institute (BGI) for whole genome sequencing using their DNBseq technology. Raw sequencing data were aligned to the human genome using bwa mem and analyzed as described in Yuan et al. [9]. Briefly, we computed the number of de-duplicated reads that aligned to each chromosome for each patient and divided this by the number of de-duplicated reads that aligned to each chromosome for a diploid germline sample from patient PW034 (pooled blood mononuclear cells) after normalizing both by total reads. We then normalized this ratio by the median ratio across all chromosomes and multiplied by two to estimate the average copy number of each chromosome.

Identification of malignant glioma cells and non-tumor cells

Chr. 7 amplification and Chr. 10 deletion were observed from the whole-genome sequencing results for all patients in this cohort. Therefore, we identified the transformed cells and untransformed cells using a linear combination of normalized average chromosome 7 and 10 expression in each cell as follows. We first merged scRNA-seq data of all samples derived from the same patient for unsupervised clustering analysis and defined putative malignant cells and non-tumor cells using the genes most specific to each cluster. Putative tumor-myeloid doublet clusters were removed prior to malignant analysis. Next, we computed the average gene expression on each somatic chromosome as described in Yuan et al. [9]. We define the malignancy score to be the log-ratio of the average expression of Chr. 7 genes to that of Chr. 10 genes and plotted the distribution of malignancy score. We fit a double Gaussian to the malignancy score distribution and established a threshold at 1.96 standard deviations below the mean of the Gaussian with the higher mean (i.e., 95% confidence interval). Putative malignant cells with malignancy scores below this threshold and putative non-tumor cells with malignancy scores above this threshold were discarded as non-malignant or potential multiplets.

For the comparison of biopsy and acute slice cultures from PW032 shown in Fig. 1, we co-clustered all of the samples together using Phenograph and identified a cluster that was statistically enriched in genes associated with red blood cells (HBA1, HBA2, HBB), transformed glioma cells (SAA1, GFAP), and myeloid cells (CD14, C1QA). We discarded these cells as potential multiplets before completing our analysis. Similarly, for the drug screening analysis of PW030 in Fig. 3, we co-clustered all eight samples using Phenograph and identified a cluster that was statistically enriched in markers of transformed glioma cells (GFAP) and myeloid cells (CD14). We discarded these cells as potential multiplets before completing our analysis.

Single-cell hierarchical Poisson factorization (scHPF) analysis

For the scHPF model in Fig. 4, we combined scRNA-seq profiles from one vehicle-treated and one etoposide-treated slice from PW029; two vehicle-treated, one etoposide-treated, and one Panobinostat-treated slice from PW030, PW032, PW034, and PW036; and two vehicle-treated and one Panobinostat-treated slice from PW040 for a total of 21 samples (see Additional file 1: Table S1). To avoid dominant factors from any one sample, we randomly sub-sampled the scRNA-seq profiles such that each of the 21 samples contributed 803 cells to the model for a total of 16,863 cells. We then factorized the resulting merged count matrix using scHPF with default parameters and K = 17 (www.github.com/simslab/scHPF) [16, 23]. For all downstream analysis of the model, we removed two nuisance factors. The first was correlated with coverage and highly ranked housekeeping genes and ribosomal protein-encoding genes, and the second contained highly ranked genes associated with cell stress and heat shock, likely a result of dissociation artifacts in a subset of cells and samples (Additional file 1: Fig. S7a). This resulted in a scHPF model with 15 factors (Additional file 2: Table S4).

Fig. 4

a UMAP embedding of scRNA-seq profiles from slice cultures of six patients generated using the cell score matrix from joint scHPF analysis of the entire data set colored by patient. b Same as a but colored by treatment condition. c Same as a but colored by the scHPF-imputed log-ratio of Chr. 7 to Chr. 10 average expression where a high ratio (red) indicates malignant transformation. d Same as a but colored by expression of the oligodendrocyte marker PLP1. e Same as a but colored by expression of the myeloid marker CD14. f Same as a but colored by the total expression of the T cell receptor constant regions (TRAC, TRBC1, TRBC2). g Heatmap showing the log-ratio of the average expression of the top 100 genes in each eptoposide-treated to each control slice for each scHPF factor and each of three cell types—transformed (tumor), oligodendrocyte (oligo), and myeloid. h Same as g for panobinostat-treated slices. i Violin plots showing the distributions of the average expression of the top 100 genes in the Proliferation scHPF factor for each vehicle- and etoposide-treated slice for each patient in tumor cells. All within-patient, vehicle-treatment comparisons have p<0.05 (Mann-Whitney U-test) unless otherwise indicated (N.S. or not significant). j Same as i for the Panobinostat1/MT scHPF factor for each vehicle- and panobinostat-treated slice in tumor cells. k Same as j for the Panobinostat2/Chemokine scHPF factor in tumor cells. l Same as j for the Panobinostat3/Oligo scHPF factor in oligodendrocytes. m Same as j for the Myeloid2/Pro-Inflammatory scHPF factor in myeloid cells. n Same as j for the Myeloid3/CD163 scHPF factor in myeloid cells

To visualize the scHPF model, we generated a UMAP embedding using a Pearson correlation matrix computed from the cell score matrix. To cluster the scRNA-seq profiles using the Phenograph implementation of Louvain community detection [20], we used the same Pearson correlation matrix and k=50 to construct a k-nearest neighbors graph. We conducted the aneuploidy analysis in Fig. 4c from the scHPF model by first computing the cell loading matrix Θ containing elements E[θ|x] for each cell-factor pair i,k and the gene sample weight matrix Β containing elements E[β|x] for each gene-factor pair g,k where x is the scRNA-seq count matrix. Next, we computed the diagonal cell scaling matrix Ξ containing elements E[ξ|x]*10,000 for each cell i and finally:

where G is the log-transformed scHPF-imputed expectation value matrix for the expression level of each gene in each cell. We colored the UMAP embedding in Fig. 4c by the difference in the average value of G for genes in chromosome 7 and that for chromosome 10. We scored each Phenograph cluster by the average of this value and took all cells in clusters with an above-average score to be malignantly transformed.

The fold-change values in the heatmaps in Fig. 4g, h were computed by dividing the average expression of the top 100 genes in each factor (rows) for the treated slice by that of each vehicle-treated slice (columns) and log-transforming. For select factors, the distribution of average expression of the top 100 genes across cells is shown for the tumor cells, oligodendrocytes, or myeloid cells for each slice in Fig. 4i–n.

Cell-type-specific differential expression analysis

To maximize our statistical power for the cell type-specific differential expression analysis shown in Additional file 1: Fig. S9, we added back all of the scRNA-seq profiles that we had subsampled out of the data set when we constructed the scHPF model, as described above. To project these additional data onto our existing scHPF model, we held the variational distributions for global, gene-specific variables fixed while updating the variational distributions for cell-specific local variables as when the model was originally trained. We used the “prep-like” command in scHPF to select the same genes that were used in the original model, and then projected the data with the “project” command in scHPF with the “—recalc-bp” option and default parameters. This results in variational approximations for the cell budgets ξ and weights θ for the additional, projected data, but does not alter the gene budgets η or weights β, nor does it alter the cell budgets or weights for the cells used to train the original model. To associate the projected cells with the originally defined Phenograph clusters, we used the “classify” command in Phenograph [20] with a Pearson correlation matrix derived from the cell score matrix computed by scHPF projection. This allowed us to assign the additional cells as transformed, myeloid, etc.

To identify differentially expressed genes for drug- vs. vehicle-treated tumor and myeloid cells, we first randomly sub-sampled the condition with a greater number of cells in each comparison to have the same number of cells as the condition with fewer cells. Next, we subsampled the count matrices for the two conditions such that they had the same average number of molecules per cell and normalized the resulting count matrix using scran [24]. We then conducted differential expression analysis for protein-coding genes using the two-sided Mann-Whitney U-test as implemented by the “mannwhitneyu” command in the Python module “scipy”. The resulting p-values were corrected for false discovery using the Benjamini-Hochberg procedure as implemented in the “mutlipletests” command in the Python module “statsmodels”. We note that this same approach was used for the differential expression analysis shown in Fig. 3b, c.

In situ hybridization, immunohistochemistry, and microscopy

To detect changes in situ of TOP2A, SOX2, CD163, MT2A, and CCL3 mRNA upon etoposide and panobinostat treatment, we performed RNAscope on vehicle- or drug-treated slices from three separate cohorts that were not processed for scRNA-seq. Treated slices were fixed in 4% PFA overnight at 4oC, paraffin-embedded, and cut into 5 mm sections. Probes against the abovementioned mRNAs were obtained from ACDBio (Additional file 1: Table S3). In situ hybridization (ISH) was performed according to the manufacturer’s protocol for the RNAScope® Multiplex Fluorescent V2 Assay (ACDBio). Briefly, serial sections were baked at 60 °C for 1 h before being deparaffinized in xylene and 100% ethanol. After drying the slides for 5 min at 60 °C, H2O2 was added for 10 min at RT. For antigen accessibility, slides were incubated in boiling 1X Target Retrieval reagents (~98 °C) for 15 min, washed in water, dehydrated in 100% ethanol, and finally treated with Protease Plus for 30 min at 40 °C. The C3 probes were diluted in C1 probes at a 1:50 ratio and incubated on the slides for 2 h at 40 °C. C1 probes were detected with TSA-Cy3 (Perkin Elmer, NEL744001KT) and C3 probes were detected with TSA-Cy5 (Perkin Elmer, NEL745001KT). DAPI was added to label the nuclei, and slides were mounted using Fluoromount. After drying at room temperature, the mounted slices were stored in the dark at 4 °C.

Immunohistochemistry was performed on 4% paraformaldehyde-fixed paraffin-embedded tissue sections (5 μm thick) to analyze γH2AX induction in etoposide-treated slice cultures. Sections were de-paraffinized in xylene (3 × 5 min), followed by rehydration in 100% ethanol (2 × 5 min), 95% ethanol (2 × 5 min), and 75% ethanol (1 × 5 min). Slides were washed in water and antigen retrieval was performed in 10 mM Citrate buffer (pH 6) in a pressure cooker for 10 min. After cooling for 30 min, slides were washed in phosphate-buffered saline (PBS; pH 7.4) and blocked in 10 % normal goat serum for 30 min. Primary antibody incubation was performed overnight at 4 °C using a mouse monoclonal anti-SOX2 antibody (clone 20G5; Thermo Scientific, MA1-014, 10 μg/mL) and a rabbit monoclonal anti-phospho-Histone H2A.X (Ser139) (clone 20E3; Cell Signaling Technology, #9718, 1:100 dilution). Following three washes in PBS, slides were incubated with goat anti-rabbit Alexa Fluor 647 and goat anti-mouse Alexa Fluor 568 secondary antibodies (Thermo Scientific; A21244 and A11031 respectively, 1:500 dilution) for 1 h. Slides were washed three times in PBS, incubated for 10 min with DAPI (Thermo Scientific, D1306, 0.5 μg/mL), and mounted using Fluoro-Gel with TES buffer (Electron Microscopy Sciences)

Images were acquired on a Zeiss LSM 800 confocal microscope with a 40×/1.3 NA oil immersion objective, using 405-nm, 561-nm, and 639-nm excitation. Five to six fields per probe were selected based on high SOX2 expression in serial sections, a pervasive marker of transformed glioma cells [9]. Confocal stacks were acquired with a 1 Airy pinhole and at 0.58-μm steps. Images were exported to ImageJ for further analysis.

To estimate cell counts, confocal images were first segmented into individual nuclei based on DAPI staining in ImageJ. First, we converted the DAPI channel into a binary image using the Huang intensity threshold followed by a watershed filter. Next, nuclear masks were generated using the “Analyze Particles” function. For each RNAscope or immunofluorescence image, we applied rolling ball background subtraction (50 pixel radius) followed by quantification of the average fluorescence intensity in each nuclear mask. To identify positive and negative nuclei for a given probe, we fit the resulting intensity distribution for all masks with an area greater than 6 square-microns to a Gaussian mixture model and used a threshold of two standard deviations above the mean intensity of the lowest intensity mode (representing the fluorescence background).

Results

Acute slice culture preserves tumor cell states and microenvironment in GBM at the single-cell level

We first demonstrated that acute slice culture preserves the cellular heterogeneity of GBM using scRNA-seq data from three uncultured biopsy specimens and three cultured slices obtained from the same patient (PW032). To identify subpopulations, we performed unsupervised clustering as previously reported [16, 20] on the entire data set containing 10,480 cells (4358 from uncultured biopsies; 6122 from cultured slices) and embedded the profiles in two dimensions using Uniform Manifold Approximation and Projection (UMAP, Fig. 1b–d) [25]. Transformed and untransformed subpopulations or clusters were distinguished by chromosome 7 amplification and chromosome 10 deletion, which were supported by both the scRNA-seq and whole genome sequencing (WGS) data (Fig. 1c, Additional file 1: Fig. S3a; see the “Methods” section). By identifying highly enriched marker genes for each cluster, the non-malignant cells were further classified into myeloid cells (CD14, AIF1, TYROBP, CD163), oligodendrocytes (PLP1, MBP, MAG, SOX10), T cells (TRAC, TRBC1, TRBC2, CD3D), endothelial cells (ESM1, ITM2A, VWF, CLDN5), and pericytes (PDGFRB, DCN, COL3A1, RGS5) (Fig. 1d, e, Additional file 1: Fig. S3b-d). We observed transformed cells and all untransformed cell types with similar fractional composition (Fig. 1f) along with expression of their marker genes (Fig. 1e) in both uncultured biopsy samples and the cultured slices. Although the biopsies and slice cultures are similar, they are not identical, and systematic shifts in gene expression are evident in the embedding in Fig. 1b and Additional file 1: Fig. S3e. Cell type-specific analysis in Additional file 1: Fig. S3b-d highlights the major gene expression differences between the biopsies and slice cultures which could result from both spatial heterogeneity and culture conditions.

In previous studies, we used scRNA-seq to show that transformed cells in high-grade gliomas resemble oligodendrocyte lineage cells (including progenitors or OPCs), astrocytes, neuronal precursors, and mesenchymal cells at the level of gene expression [9], consistent with earlier work using bulk analysis [26-28]. More recently, Neftel et al. developed an elegant model based on gene signatures of these four major states to classify scRNA-seq profiles of glioma cells [11]. We used this model to examine the transformed cells in the biopsies and slice cultures in detail and found that all four states were well-represented in both the biopsy and slice cultures and that most cells classified as astrocyte-like or mesenchymal in this particular tumor (Fig. 1g, h, Additional file 1: Fig. S3f). However, the slice cultures contained more mesenchymal cells whereas the biopsy cells were more astrocyte-like (Fig. 1g, h, Additional file 1: Fig. S3f). While this could be due, in part, to culture conditions, we expect this level of variation based on previous studies of spatial heterogeneity since the slice cultures and biopsies were obtained from different regions of the tumor [8, 16, 27, 29].

We conducted a similar analysis of slice cultures from six patients (including PW032) and found representation of transformed cells (Additional file 1: Fig. S4) and the same major cell types with similar relative abundances after 24 h of culture (Fig. 2a). We also analyzed their transformed populations using the Neftel et al. model (Fig. 2b, Additional file 1: Fig. S5) [11] and found good representation of all four major GBM states with some tumors appearing more proneural (OPC/NPC – PW034) and others more astrocytic (PW029), mesenchymal (PW030, PW040), or both (PW032, PW036) as quantified in Additional file 1: Fig. S4d. Finally, to analyze spatial effects across slice cultures within a single resection, we profiled five slice cultures such that each was 500-μm thick and the interval between the two most adjacent slices was also 500 μm (maximum spatial distance of 3.5 mm). scRNA-seq profiles of the five slices co-clustered well based on UMAP embedding (Fig. 2c, Additional file 1: Fig. S6a-e) and the four-state model of the transformed cells (Fig. 2d, Additional file 1: Fig. S6f) and showed good representation of major cell types (Additional file 1: Fig. S6). Taken together, these data suggest that cultured slices preserve the major cellular and molecular features of the tumor microenvironment and represent the well-established inter- and intra-tumoral heterogeneity observed in gliomas.

Fig. 2

a Fractional abundance of each major cell type in all untreated slice culture scRNA-seq data sets from the six patients in the study. b Two-dimensional model projecting each transformed cell from all untreated slice culture scRNA-seq data sets from the six patients in the study. c UMAP embedding of scRNA-seq profiles from five untreated slice cultures taken within 3.5 mm of each other from PW040 colored by sample of origin. d Same as b but for the transformed cells from the five untreated slice cultures from PW040

Acute slice culture and scRNA-seq for personalized drug screening

To test the feasibility of drug screening with patient-derived slice cultures and scRNA-seq, we perturbed slices derived from one GBM resection (PW030) with six different drugs chosen for diverse mechanisms of action and included two vehicle controls (Fig. 3a and Additional file 1: Table S1, S2). We profiled 48,404 cells from eight slices and identified transformed and untransformed populations as described above. To identify drug-induced transcriptional changes, we performed differential expression analysis for the tumor, myeloid, and oligodendrocyte populations (Fig. 3b, c). Treatment with the histone deacetylase (HDAC) inhibitor panobinostat resulted in the strongest response with 9632, 4228, and 3183 significantly differentially expressed genes (p<0.01) in the tumor, myeloid, and oligodendrocyte populations (Fig. 3b), respectively, with similar results when we restricted our analysis to fold-changes >2 (Fig. 3c). To identify drugs with highly specific effects on subpopulations of tumor cells, we first computed a UMAP embedding for the transformed cells from the control slices (Fig. 3d), which we use as a reference for comparison to treated cells. The majority of control tumor cells appear mesenchymal (Fig. 2b) with pervasive expression of CD44 and VIM and an astrocytic subpopulation expressing GFAP at high levels (Fig. 3d). However, there is a small subpopulation of proliferating cells marked by TOP2A and MKI67 (Fig. 3d). Next, we projected the profiles of transformed cells from each treated slice into this embedding (Fig. 3e). Consistent with the differential expression analysis, we observed that panobinostat had a dramatic compositional effect on the transformed cells. We also noticed that etoposide selectively eliminated the small, proliferative subpopulation, consistent with its mechanism-of-action as a topoisomerase poison [30]. Given the disparate effects of these two drugs in PW030, we made them the focus of our subsequent analysis.

Conserved, cell type-specific responses to etoposide and panobinostat in GBM across six patients

To identify cell type-specific responses to etoposide and panobinostat that are conserved across patients, we conducted slice culture drug-perturbations across six GBM patients followed by scRNA-seq. Additional file 1: Table S1 contains a summary of all the slice culture samples used for this analysis. After subsampling the scRNA-seq profiles from each vehicle- and drug-treated slice culture from all six patients to the same number of cells, we generated a low-dimensional representation of the merged data using single-cell hierarchical Poisson factorization (scHPF) [16]. This Bayesian algorithm operates directly on the count matrix and identifies latent factors corresponding to the major gene expression programs that define the population. We identified 15 factors associated with canonical markers of neural cell types, GBM subpopulations, biological processes (e.g., proliferation), and drug response (Additional file 1: Fig. S7, Additional file 2: Table S4). Two nuisance factors were associated with coverage (enriched in ribosomal and other housekeeping genes) and cell stress (e.g.. heat shock—likely a dissociation artifact) and removed from the model (Additional file 1: Fig. S7a, Additional file 2: Table S4).

To visualize the model, we created a UMAP embedding of the scHPF cell score matrix (Fig. 4a–f, Additional file 1: Fig. S8; see the “Methods” section). Based on aneuploidy analysis of chromosomes 7 and 10, the transformed cells from each patient separate into essentially non-overlapping clusters (Fig. 4a, c). In contrast, untransformed oligodendrocytes, myeloid cells, and T cells overlap significantly across the six patients (Fig. 4a, d–f), consistent with previous studies of fresh resections [9]. We also observed that panobinostat- and vehicle-treated cells generally showed little overlap across all cell types, while etoposide-treated cells tended to co-cluster with the controls (Fig. 4b). This is consistent with our screening results above which suggest that panobinostat significantly alters gene expression, whereas etoposide primarily impacts genes involved in proliferation.

To identify conserved drug responses, we compared the expression of top genes in each factor between the drug- and vehicle-treated slices from each patient. As expected, the most conserved response to etoposide was a decrease in expression of the proliferation factor in the tumor compartment (Fig. 4g, i). This occurred in all but one patient, PW034, despite its high levels of cells expressing proliferation markers (Fig. 4i). Etoposide did not show consistent effects on other factors and had limited impact overall on oligodendrocytes and myeloid cells (Fig. 4g). We validated the loss of TOP2A+ tumor cells, which we also observed by conventional differential expression analysis (Additional file 1: Fig. S9a), using in situ hybridization analysis of etoposide-treated slice cultures from a separate cohort (Additional file 1: Fig. S10). These results suggest that the alterations in transformed glioma cell-specific gene expression results from an etoposide-mediated change in cellular composition. We also show the expected DNA damage response of transformed glioma cells to etoposide treatment using double immunofluorescence staining of slice cultures with γH2AX and SOX2 (Additional file 1: Fig. S11). γH2AX is an established marker of DNA damage and repair that is known to respond to etoposide [31]. We previously showed the SOX2 is a pervasively expressed marker of transformed glioma cells in GBM using scRNA-seq and immunohistochemistry [9]. In contrast to etoposide, panobinostat affected multiple factors for both tumor and non-tumor cells (Fig. 4h). We observed a modest decrease in expression for the proliferation factor across all patients except for PW040 (Fig. 4h). Interestingly, panobinostat induced expression of LEFTY1, BEX5, and SAXO2 as part of the Panobinostat3/Oligo factor, which was predominantly oligodendrocyte-specific (Fig. 4h, Additional file 1: Fig. S7b). However, the most notable effect was upregulation of metallothionein family genes (Panobinostat1/MT factor) across all cell types (Fig. 4h, j), consistent with previous reports that HDAC inhibitors can perturb this highly inducible gene cluster [32, 33]. Interestingly, cell type-specific differential expression analysis not only confirmed metallothionein induction but also revealed upregulation of mature neuronal genes (e.g., SNAP25, SLC17A7, KCNB1, RAB3A), a component of the same factor, specifically in tumor cells (Additional file 1: Fig. S9c). Because we observed metallothionein induction in all six patients and the three major cell populations analyzed here, it is a potentially useful biomarker of panobinostat response.

Panobinostat treatment significantly impacted gene expression in myeloid cells. In slice cultures from 3/5 patients, we observed a modest decrease in a factor marked by pro-inflammatory cytokines (Fig. 4h, m), which have been shown to be predominantly expressed by microglia in the glioma microenvironment [9, 34]. We observed a more consistent effect on a myeloid factor marked by CD163, which is likely expressed by macrophages which are thought to be immunosuppressive in GBM (Fig. 4h, n) [34]. We verified that CD163 exhibited significant, myeloid-specific differential expression (Additional file 1: Fig. S9d) and the loss of CD163+ macrophages in general and relative to CCL3+ pro-inflammatory myeloid cells by in situ hybridization analysis of vehicle and panobinostat-treated slice cultures from a separate group of patients (Additional file 1: Fig. S12). In addition, we validated the widespread induction of metallothionein by panobinostat using in situ hybridization of MT2A (Additional file 1: Fig. S13).

Cell type-specific responses to etoposide and panobinostat within an individual patient

In the experimental design above, the majority of replicates occur across patients. While this enables analysis of conserved drug responses, we cannot assess the consistency of a cell type-specific response within an individual patient. Figure 5 shows the results from an alternative experimental design in which we perturb multiple slice cultures for each drug with spatially adjacent vehicle control replicates so that we can assess the reproducibility of cell type-specific drug response within an individual resection. We tested this approach on etoposide- and panobinostat-treated slice cultures generated from a single GBM patient with three replicate slices for each drug and four vehicle-treated slices. We observed a significant reduction in cell viability in the etoposide- and panobinostat-treated slice cultures compared to vehicle controls (Additional file 1: Fig. S14). Figure 5a shows a UMAP embedding of the corresponding scRNA-seq data. As in the above experiments, we can use aneuploidies in chromosomes 7 and 10 to identify malignantly transformed glioma cells (Fig. 5b) and other markers to identify myeloid cells, T cells, and oligodendrocytes in the microenvironment (Fig. 5c, Additional file 1: Fig. S15). Unlike the primary tumors described above, these slice cultures were generated from a surgical resection of a recurrent tumor and showed typical features of recurrent GBM including high levels of myeloid infiltration [9, 27]. We first performed differential expression analysis between the transformed tumor cells from all of the etoposide- and vehicle-treated slices (Fig. 5d). Consistent with our previous analysis of conserved drug responses, we find the etoposide significantly reduces the expression of proliferation markers (Fig. 5d), suggesting a robust cell type-specific response in this specific patient sample. We repeated this analysis for the transformed GBM cells in the panobinostat-treated slices and found consistent induction of metallothioneins and neuronal markers as expected (Fig. 5e). We also analyzed the effects of panobinostat on myeloid cells, and just as in the above analysis, we found expression changes that are consistent with a loss of CD163+ macrophages (Fig. 5f). Finally, to assess reproducibility across replicates, we performed gene set enrichment analysis (GSEA) with gene sets taken from the top 100 scoring genes in each of the scHPF factors from our analysis in Fig. 4 on the differentially expressed genes from independent comparisons of etoposide- (Fig. 5g) and panobinostat-treated (Fig. 5h) slice cultures to adjacent control slices. For the etoposide-treated slices, we found that the proliferation signature in the transformed glioma cells was the only factor with significant reduction across all three replicates (FDR<0.05). For panobinostat, the CD163+ macrophage signature was reduced in myeloid cells in all three replicates with FDR<0.05 in two out of three, while the metallothionein signature was increased in both myeloid and tumor cells for all three replicates. with FDR<0.05 in two out of three. These results suggest that our conserved findings from primary GBM are applicable to recurrent GBM and demonstrate that we can use our approach to identify robust, cell type-specific drug responses within an individual patient.

Fig. 5

Ten slices from a single patient (TB6393) were treated with panobinostat (3 slices), etoposide (3 slices), or vehicle (DMSO, 4 slices adjacent to drug-treated slices). a UMAP embedding of scRNA-seq profiles of ten slices colored by treatment condition. b Same as a but colored by the log-ratio of Chr. 7 to Chr. 10 average expression where a high ratio (red) indicates malignant transformation. c Same as b but colored by cell type. d Volcano plot of differential expression analysis between all transformed tumor cells from etoposide-treated slices and all adjacent vehicle-treated slices. Genes highlighted in red and blue have fold-increase or decrease, respectively, greater than two and FDR<0.05. A large set of cell cycle control markers highly downregulated in the etoposide-treated cells are labeled and highlighted in cyan. e Same as d but for panobinostat-treated transformed tumor cells showing strong induction of metallothioneins and several mature neuronal markers labeled and highlighted in orange. f Same as e but for the myeloid cells showing downregulation of the macrophage markers highlighted in cyan and strong induction of metallothioneins highlighted in orange. g Heatmap showing the normalized enrichment score (NES) from gene set enrichment analysis (GSEA) analysis. GSEA was performed using gene sets from the top 100 genes of each scHPF factor from Fig. 4 to analyze the ranked differentially expression genes between tumor or myeloid cells from each etoposide-treated slice and that of its adjacent vehicle-treated slice. scHPF factor with consistent enrichment and FDR<0.05 in at least 2 treated vs. untreated comparisons are marked with asterisk. h Same as g but showing NES from GSEA analysis for differentially expression genes between tumor or myeloid cells from each panobinostat-treated slice and that of its adjacent vehicle-treated slice

Discussion

Collectively, this work establishes a multiplexed experimental and analytical pipeline for deconvolving cell type-specific drug responses in GBM tissue from individual patients. Acute slices generated from fresh tumor tissues preserve the key molecular and cellular features of the original tissue and provide a setting for drug response to be evaluated on multiple tumor cell subpopulations and cell types in the tumor microenvironment. We further demonstrated the feasibility of conducting drug screens using this approach with a turnaround time of less than 1 week after surgery. Focused analysis of etoposide and panobinostat across six patients (five for each drug) revealed drug-induced responses in specific populations of transformed and microenvironmental cells, patient-specific drug sensitivities, and drug effects conserved across patients. Etoposide consistently downregulated cell cycle genes in proliferating tumor cells with minimal conserved effects on untransformed or less proliferative transformed cells. The HDAC inhibitor panobinostat induced the expression of metallothionein family genes and mature neuronal genes in tumor cells and significantly re-modeled the myeloid population in the tumor microenvironment.

We note that drug delivery across the blood brain barrier is a major challenge in the treatment of GBM, and that local delivery may be necessary for many drugs to be effective. Our slice culture experiments are only informative of how a given therapy might impact malignantly transformed GBM cells and cells in the brain microenvironment once it has been effectively delivered to the tumor.

Conclusions

Acute slice cultures recapitulate the major cellular and molecular features of both transformed and untransformed cells in GBM. In combination with scRNA-seq, we can conduct drug screens directly on intact human surgical specimens and deconvolve cell type-specific drug responses. Overall, we hope that this approach will find broad utility for pre-clinical studies and the development of cellular and molecular enrollment criteria for clinical trials.

Additional file 1: Supplementary Materials For: Deconvolution of Cell Type-Specific Drug Responses in Human Tumor Tissue with Single-Cell RNA-seq. This file includes supplementary figures: Fig. S1-S15. and supplementary tables: Table S1-S3.

Additional file 2: Table S4. Gene score matrix for each scHPF factor in the model shown in Fig. 4.