Comparison of high-throughput single-cell RNA-seq methods for ex vivo drug screening

Abstract

Functional precision medicine (FPM) aims to optimize patient-specific drug selection based on the unique characteristics of their cancer cells. Recent advancements in high throughput ex vivo drug profiling have accelerated interest in FPM. Here, we present a proof-of-concept study for an integrated experimental system that incorporates ex vivo treatment response with a single-cell gene expression output enabling barcoding of several drug conditions in one single-cell sequencing experiment. We demonstrate this through a proof-of-concept investigation focusing on the glucocorticoid-resistant acute lymphoblastic leukemia (ALL) E/R+ Reh cell line. Three different single-cell transcriptome sequencing (scRNA-seq) approaches were evaluated, each exhibiting high cell recovery and accurate tagging of distinct drug conditions. Notably, our comprehensive analysis revealed variations in library complexity, sensitivity (gene detection), and differential gene expression detection across the methods. Despite these differences, we identified a substantial transcriptional response to fludarabine, a highly relevant drug for treating high-risk ALL, which was consistently recapitulated by all three methods. These findings highlight the potential of our integrated approach for studying drug responses at the single-cell level and emphasize the importance of method selection in scRNA-seq studies. Finally, our data encompassing 27 327 cells are freely available to extend to future scRNA-seq methodological comparisons.

Introduction

Functional precision medicine (FPM) is a strategy whereby live cancer cells are perturbed with drugs, ex vivo, which has translational potential to provide personalized information for guiding therapy in vivo (1). The model of applying a drug to cancer cells and observing its behavior is not a new concept (2). Early initiatives towards clinical translation faced substantial problems, including poor growth of cells in ex vivo cultures, low throughput culture conditions, and limited availability of new drugs to analyze. The past decade yielded dramatic increases in new high throughput technologies and a dramatic increase in the number of drugs available, as well as widespread adoption of comprehensive genomic analysis of cancer samples (3–7). This has dramatically revived interest in ex vivo analysis as a method to guide clinical treatment (8–10).

In parallel, single-cell sequencing has greatly improved our understanding of molecular heterogeneity in cancer (11,12). In the context of FPM, identifying the genes and pathways underlying the heterogeneous behaviors of individual cells in response to a specific treatment is an important step to improve biological understanding and efficacy of treatment (13). Combining FPM with single-cell output for investigating cell-specific ex vivo treatment response holds great promise for understanding cellular responses to drugs that may be overlooked in bulk assays (14–16). While ex vivo drug screening has advanced our understanding of cellular responses to drugs, significant challenges remain with incorporating single-cell RNA-seq (scRNA-seq) as an output. One major difficulty is the high cost and technical complexity associated with single-cell technologies, including library preparation and the need for accurate cellular barcoding to ensure the scalability of these methods to handle large sample sizes and high-throughput screening. Additionally, there is a need for standardized protocols and benchmarks to ensure reproducibility and comparability of results across different scRNA-seq studies and platforms (17). To address this, we evaluate three different high throughput scRNA-seq methods that enable cellular barcoding to track cell-drug interactions in the context of the ex vivo drug screening.

Materials and methods

Cell culture

The human glucocorticoid-resistant acute lymphoblastic leukemia (ALL) E/R+ Reh cell line (DSMZ, ACC 22) was maintained in RPMI-1640 medium (Sigma, R0883) supplemented with 2 mM l-glutamine (Sigma, G7513), 100 U/ml Penicillin, and 100 μg/ml Streptomycin (Sigma, P0781), along with 10% heat-inactivated fetal bovine serum (HIFBS; Sigma F9665). Cells were cultured at 37°C in humidified air (95%) and 5% CO₂. Cells were split twice a week, with the first two subcultures after thawing the cell line supplemented with an additional 10% HIFBS in the culture media to enhance cell growth. The Reh cell line was verified by karyotyping (https://doi.org/10.1101/2023.03.08.531483) and the Eurofins human STR profiling cell line authentication service (18).

Ex vivo drug treatment

The Fluorometric Microculture Cytotoxicity Assay (FMCA) was conducted following previously established protocols (9,19–21). Briefly, Reh cells were resuspended in cell culture media and seeded into a 384-well plate at a density of 7000 cells in 45 μl of media per well (0.156 × 10⁶ cells/ml). A total of 41 compounds, including positive and negative controls, were added to the plate using an Echo 550 liquid handler (Beckman Coulter, Brea, CA, USA), following a 1:3 dilution series in a five-step gradient. Subsequently, the plates were incubated for 72 h in a controlled environment (37 °C, 95% humidity, 5% CO₂). After the incubation period, the cellular response was assessed using a fully automated SCARA system (Beckman Coulter, Brea, CA, USA). The survival index (SI%) was determined by calculating the fraction of surviving cells relative to untreated cells. A comprehensive list of the compounds and their concentrations can be found in Supplementary Figure S1.

FMCA prior to scRNA-seq

Reh cells were seeded in a six-well plate at a density of 0.16 × 10⁶ cells per well, with 3 ml of cell culture media added to each well. Two six-well plates were prepared during the cell culture split, with each plate containing five wells for drug treatment and one well serving as a vehicle control treated with dimethyl sulfoxide (DMSO). The selected drugs and concentrations used in the experiment were as follows: 5 μM 5-azacytidine (A2385, Sigma), 10 μM dexamethasone (D1756, Sigma), 10 μM imatinib (I-5508, LC Lab), 1 μM XRP44X (X3129, Sigma-Aldrich), and 0.56 μM fludarabine (S1391, Selleckchem). The control well was treated with 0.9% (v/v) DMSO (D5879, Honeywell). Following a 72-h incubation (37°C, 95% humidity, 5% CO₂), the cells were harvested for analysis.

Cell harvesting

The FMCA and cell harvesting experiments were conducted in a consistent manner for each of the three scRNA-seq methods investigated. Following the treatment period, cells were harvested by transferring to 15 ml polypropylene tubes containing phosphate buffered saline (PBS). Prior to cell harvesting, the PBS was pre-warmed to 37°C. The harvested cells were centrifuged at 200 g for 10 min, and then resuspended in 2.0 ml of PBS. The concentration and viability of cells in each treatment condition were assessed using AO-DAPI staining and analyzed on a NucleoCounter-3000 (Chemometec).

Parse Biosciences Evercode library preparation (Parse)

The cells recovered from each treatment condition (range 616k-1720k) were fixed according to the manufacturer's protocol (Parse Biosciences) including 0.5% BSA in the fixation solution. After the fixation the cell concentration was measured for each suspension before being stored at –80°C. For barcoding, cell suspensions with equal number of cells from each treatment condition were thawed and split to two wells of the first barcoding plate of an Evercode whole transcriptome V1 mini kit (Parse Biosciences). Barcoding and sequencing library generation were performed according to the manufacturer's protocol, generating two sub-libraries, with a target of 5000 cells each (10000 cells in total). Prior to library generation the concentration of the pooled and barcoded cell suspension was evaluated by trypan blue (EVE cell counter) and propidium iodide label (Cellometer K2). The input cell concentration of the first sub-library (P1) was based on the trypan blue estimate and the second sub-library (P2) was based on the propidium iodide label estimate. The two final sequencing sub-libraries were evaluated on a TapeStation (Agilent) and quantitative real time PCR (qPCR) before being pooled together at equimolar concentrations.

MULTI-seq library preparation (MULTI-seq)

A volume equivalent to 500k cells from each treatment condition was transferred to six 15 ml tubes, and PBS was added to achieve a total volume of 5 ml. The cells were isolated by centrifugation at 200 g for 10 min and subsequently transferred to loBind 1.5 ml Eppendorf tubes. Barcode labeling was performed using lipid modified oligos (LMO, Sigma) as previously described (22). The processing was carried out in six parallel tubes, optimizing wash steps and volumes to minimize handling time and cell loss. The final resuspension of cells was performed in PBS supplemented with 1% BSA and 0.5 U/μl RNase inhibitor, using a low volume of 100 μl/sample. The LMO-barcoded samples were subsequently pooled equally by volume. Dead cells were removed by incubating with DAPI for 5 min, followed by negative selection using gentle sorting (Miltenyi MACS Quant Tyto). The concentration and viability of the purified cell suspension pool were evaluated using AO-DAPI staining (NucleoCounter-3000). The suspension pool with a volume equivalent to a target of 6000 cells was immediately processed on one lane of a Chromium next GEM Single cell 3′ gene expression v3.1 dual index kit (10x Genomics) following the manufacturer's protocol. During the initial cDNA purification step, the fraction containing the short cDNA, which contained the barcode libraries, was used for generating a barcode library according to the MULTI-seq protocol. The gene expression library and barcode library were evaluated using a TapeStation (Agilent) and qPCR before being pooled.

10x Genomics Single Cell Gene Expression Flex library preparation (10x Fixed)

Cells from each treatment condition (518k–1390k cells) were fixed with formaldehyde and additive according to the Chromium Next GEM Single Cell Fixed RNA Sample Preparation Kit (10x Genomics). Fixation was performed at +4°C overnight. After buffer exchange, the fixed cells were stored at –80°C until processing. A four barcode (BC) version of Chromium Fixed RNA Kit, Human Transcriptome (10x Genomics) was used for subsequent barcoding of the six different samples, i.e. three BC sets were combined in each 10x reaction (n = 2 libraries). The post hybridization washes were done individually for each barcoded sample, according to the manufacturer's protocol. Cell counting was performed using propidium iodide labeling (Cellometer K2). A target of 2500 cells per treatment condition were combined for a total of 7500 cells per library using a Chromium iX instrument. The two resulting libraries were generated according to the manufacturer's protocol. The libraries were evaluated using a TapeStation (Agilent) and qPCR before being pooled.

Sequencing

Sequencing was performed at the SciLifeLab National Genomics Infrastructure in Uppsala, Sweden. Each pool of libraries was sequenced independently by method on one lane of a NovaSeq6000 S Prime (SP) flowcell (Illumina), for three SP lanes in total. A paired-end read setup of 100–10–10–100 bp (Read1-index1-index2-read2) for Parse and MULTI-seq or 150–10–10–150 bp for 10x Fixed, was used with v1.5 sequencing chemistry. Base calling was performed with RTA v3.3.4 and bcl2fastq 2.20.0.422.

Data processing

The data from chromium-based methods (MULTI-seq and 10x Fixed) were demultiplexed and mapped with Cellranger v7.0.1 using a GRCh38 reference provided by 10x Genomics (refdata-gex-GRCh38-2020-A). For MULTI-seq, the gene expression analysis was run in parallel with feature barcode analysis. The feature barcode matrix was subsequently loaded into R where the two assays were merged and samples were demultiplexed using HTODemux from the Seurat package (23). Cells classified as negatives or doublets were removed. Parse data was demultiplexed and mapped with the ParseBiosciences-Pipeline v0.9.6p using the GRCh38 reference (refdata-gex-GRCh38-2020-A). The unfiltered count matrix was used for downstream analysis. Further analysis of the raw count matrices was performed in R using Seurat, harmony and ggplot2 packages in R. The inflection point of the total count versus barcode rank curve was determined for each library using the barcodeRanks function from DropletUtils (24). Cells with a total UMI count lower than the inflection point and cells with >10% of the total counts originating from mitochondrial genes were removed (Supplementary Figure S2). Following cell filtering, libraries from the Parse and 10x Fixed experiments were each merged to create a single dataset per method. Log-normalization and scaling was performed with Seurat v4 and variance stabilizing transformation was performed with SCTransform v2 without covariates (25,26). The data from the three scRNA-seq methods were combined using Seurat prior to clustering and visualization.

Data analysis

Down sampling. The libraries were downsampled using the downsampleReads function from DropletUtils (24). Down-sampling was performed on a cell-by-cell basis generating a downsampled datasets where all cells with read counts above the down sampling target had the same total read count and cells with read counts below the down sampling target were discarded. The downsampleReads function uses the hdf5 molecule information file generated by the Cellranger pipeline as input for MULTI-seq and 10x Fixed. For the Parse libraries, the tscp assignment.csv.gz was converted to hdf5 format with a custom a python script.

Dropout rate. Dropout was calculated in a similar approach to Ziegenhain et al. (27): First, we randomly subsampled 500 DMSO-treated cells per scRNA-seq method for a total of 1500 cells. Next, we selected the genes detected (non-zero count) in at least 25% of all cells by any scRNA-seq method. Finally, we calculated the dropout rate by scRNA method and by gene as the proportion of cells with zero counts.

Cell cycle scores. Cell cycle scores were computed using the CellCycleScoring function and marker genes classify cells into G1, G2M or S states in Seurat (23).

Differential gene expression. Differential gene expression was determined using MAST (28). Differentially expressed genes were determined by comparing each drug treatment separately to control (DMSO).

Pathway analyses. The R-package clusterProfiler was used for over-representation analysis for GO biological process, WikiPathways, KEGG and CMAP on the genes determined to be differentially expressed (fold change > 0.50, adjusted P< 0.01) as compared to DMSO controls (29,30). Significant overlaps (FDR q-values < 0.05) between differentially expressed genes and Hallmark gene sets were computed for each sequencing method by Gene Set Enrichment Analysis (31).

Cell cycle analysis by DNA stain

Reh cells were seeded in a six-well plate at a density of 0.16 × 10⁶ cells per well, with 3 ml of cell culture media added to each well and treated with 5 μM 5-azacytidine (A2385, Sigma), 10 μM dexamethasone (D1756, Sigma), 10 μM imatinib (I-5508, LC Lab), 0.56 μM fludarabine (S1391, Selleckchem) or 0.9% (v/v) DMSO (D5879, Honeywell). Following a 72-h incubation (37°C, 95% humidity, 5% CO₂), the cells were harvested (see Cell Harvesting above). For each treatment 0.5 × 10⁶ cells were lysed and stained with 10 μg/ml DAPI according to the ‘Two-step cell cycle analysis’ application note (Chemometec). The fluorescence distribution was quantified in a NC-3000 (Chemometec) image cytometer.

Results

Functional drug resistance testing by FMCA

We performed a systematic methods comparison study of three approaches for multiplexing scRNA-seq experiments in an ex vivo drug screening model system (FMCA). To identify drugs suitable for the method comparison, we initially assessed the ex vivo cellular drug resistance of the glucocorticoid-resistant E/R+ ALL Reh cell line using FMCA. A panel of 41 drugs from diverse classes, including targeted drugs, glucocorticoids, and chemotherapeutic agents, was employed at varying concentrations (Supplementary Figure S1). Based on the dose-response curves after 72 h of treatment, five drugs were selected for further investigation of their transcriptional effects (Figure 1). The chosen drugs were as follows: 5-azacytidine, a nucleoside analogue that inhibits DNA methyltransferase, resulting in DNA hypomethylation and cell death at high doses (32,33). Dexamethasone, a glucocorticoid, was included as a negative control since the Reh cell line is known to exhibit glucocorticoid resistance (34,35). Imatinib, a tyrosine kinase inhibitor used in the treatment of ALL patients with ABL-class fusions, although the Reh cell line does not possess such fusions; its inclusion was based on the responsive profile (Figure 1). XRP44X, an inhibitor of Ras-activated Elk3 transcription-factor activity, was selected due to a previous observation in the Reh cell line (36). Fludarabine, a purine analog, is currently employed in the treatment of high-risk ALL patients as part of conditioning before stem cell transplantation and lymphodepletion before CAR-T therapy (37). For each drug, a concentration was determined for subsequent single-cell analysis, aiming to achieve a SI% of ∼50% to balance between inducing a cellular response, while maintaining a sufficient proportion of surviving cells for robust assessment of transcriptional effects.

Figure 1.

Fluorometric Microculture Cytotoxicity Assay (FMCA) survival indexes (SI%) for five selected drugs. The SI% (y-axis) after exposed to 72 h treatment at a range of concentrations (x-axis). Each dot represents the mean SI% measured from two replicates. The dotted vertical line denotes the concentration selected for single-cell transcriptomic analyses (scRNA_seq). The results from the FMCA experiment used to generate each of the scRNA-seq datasets are indicated according to the legend in the right side of the panel.

Multiplexed scRNA-seq library preparation

We evaluated the ex vivo transcription profiles of the Reh cells after treatment with the five aforementioned drugs and DMSO controls. Three library preparation protocols were tested: (i) the Parse Biosciences Evercode Whole Transcriptome mini v1 (Parse) (38), (ii) the cellular LMO-barcoding approach in combination with scRNA-seq using the 10x Genomics Chromium 3′ GEX V3 protocol (MULTI-seq) (22) and (iii) the Chromium Single Cell Gene Expression Flex kit from 10x Genomics (10x Fixed).

To capture subtle differences in cell states by each drug condition, we aimed to profile 6000–15000 cells with each scRNA-seq protocol (1000–2500 cells per drug condition), based on each kit's ability to incorporate a drug-specific cellular barcode to mark surviving cells according to which drug they were exposed. Due to protocol differences and restraints in number of sample barcodes available for each kit, the library set-up for each method was slightly different. For Parse we split each sample into two initial barcoding wells (technical replicates) and in the final step before library amplification the pooled cells were split into two sublibraries containing all six conditions (two sub-libraries: P-1 and P-2). For MULTI-seq we performed a single experiment (M-1). For 10x Fixed we used three barcodes (corresponding to 5-azacytidine, XRP44X, dexamethasone) in library F-1 and three barcodes (corresponding to imatinib, fludarabine and DMSO control) in library F-2. Differences in sequencing depth, UMIs and number of genes detected prior to cell filtering and library merging were observed between the 10x Fixed (F-1 and F-2) libraries and the Parse (P-1 and P-2) sublibraries (Supplementary Figure S2). The differences in sequencing depth per cell of P-1 and P-2 sublibraries was correlated with cell concentration measurement as described below. The characteristic inflection points in the barcode rank plot (knee plot) that separate cells from background noise demonstrate superior signal-to-noise in the 10x Fixed and MULTI-seq approaches in comparison to Parse (Supplementary Figure S3). After removing non-cell UMIs, cells with > 10% mitochondrial gene reads, and merging the sub-libraries we characterized 27327 cells in total: 7225 cells with the Parse method, 10516 with MULTI-seq, and 9586 with 10x Fixed (Table 1).

Table 1.

Single-cell RNA-sequencing libraries and statistics

Method	Parse		MULTI-seq	10x Fixed
Sub-library name	P-1	P-2	M-1	F-1	F-2
Number single cells targeted	5000	5000	6000	7500	7500
Number single cells identified	2296	4929	10516	4257	5329
Total reads sequenced/cells identified	111k	42k	43k	41k	62k
Median reads per cell	26771	9605	27483	32241	49279
Median genes per cell	3772	2698	4167	6262	6927

Evaluation of library complexity

Since several approaches were used for cell counting prior to library preparation, we estimated how well each method performed. In our hands, the propidium iodide labeling used to count cells in P-2, F-1 and F-2 more accurately estimated the cell concentration, while quantification by trypan blue overestimated the cell concentration, resulting in underloading of library P-1, while AO-DAPI underestimated the cells loaded in M-1, resulting in overloading. We therefore conclude that accurate cell counting is important for downstream data quality and propidium iodide labeling was the most reliable method for determining cell concentration prior to scRNA-seq library preparation for fixed cells.

The proportion of raw reads retained for downstream analysis after demultiplexing and quality filtering varied by scRNA-seq method (Figure 2A). The most sequencing reads were retained in 10x Fixed (85.5 %), followed by MULTI-seq (62.2 %) and Parse (51.6 %). We further examined read retainment through the different steps of the pipeline. Incomplete or non-matching sample or cell barcode sequences were the primary contributors to read loss in the Parse and MULTI-seq libraries (Supplementary Figure S4). The proportion of filtered reads that mapped to exonic, intronic or intergenic regions also varied between the libraries generated by the three methods, which likely reflects differences in the method design (Figure 2B). The 10x Fixed approach resulted the highest proportion of exonic reads (98.7 %), consistent with the assay using probes targeting exons. MULTI-seq, based on the 10x Genomics 3′ assay that uses oligo (dT) primers, produced 75% reads mapped to exons in line with what others have shown previously (17). Parse uses random primers in addition to oligo(dT), resulting in the majority of reads mapping to intronic regions (60.4 %), which is in line with what other reports using this method have observed (https://doi.org/10.1101/2023.06.28.546827).

Figure 2.

Quality metrics of single cell RNA sequencing data generated from three different library preparation methods (Parse, MULTI-seq, 10x Fixed). (A) Stacked bar plot showing the percentage of sequencing reads lost across each step of data processing. (B) Proportion of mapped reads that map to exonic, intronic and intergenic regions by scRNA-seq method. (C) Violin plots of the of mitochondrial (left) and ribosomal protein genes (right) read fractions separated by scRNA-seq method. (D) The number of genes detected (y-axis) plotted against the number of UMIs (x-axis) for each scRNA-seq library preparation method. The numbers are visualized as hexagonal binning plots, the color indicates the number of cells contributing to each hexagonal bin according to the figure legend scale to the right of the panel. The total number of cells captured by each method is indicated in brackets. (E) Violin plot showing the distribution of genes detected per cell for each scRNA-seq method. (F) Pairwise scatterplot of the pseudobulk gene expression level per gene and method. The Pearson's correlation coefficient (ρ) is indicated in each panel.

The presence of a high mitochondrial gene read fraction in scRNA-seq data can indicate cellular stress or cell death, while the read fraction of ribosomal protein genes can vary across cell state and may be an indication of RNA degradation (39,40). To assess these parameters, we compared the proportion of reads mapping to mitochondrial and ribosomal protein genes across all cells by scRNA-seq method. Notably, the Parse method exhibited the lowest mitochondrial (median 0.39% Figure 2C) and ribosomal protein genes (median 0.42%) read fraction in the filtered data as well as prior to the initial UMI filtering step (Supplementary Figure S5). This may be attributed to the multiple washing steps incorporated in the Parse workflow, which likely minimizes ambient RNA and the impact of mitochondrial/ribosomal protein genes expression. The MULTI-seq library contained a low read fraction of mitochondrial genes (median 1.55%), but higher read fraction of ribosomal protein genes (median 22.37%). The low read fraction of mitochondrial genes likely is a result of the dead cell removal step used for the MULTI-seq library, where the cells were >99% viable prior to cell capture. The Chromium 3′ chemistry, used in MULTI-seq captures a higher proportion of ribosomal protein gene reads compared to other scRNA-seq methods (https://doi.org/10.1101/2023.06.28.546827) (41). The 10x Fixed libraries contained the highest read fraction of mitochondrial genes (median 3.53%), even though the method had the shortest time between cell isolation and fixation. Unlike the other protocols, the10x Fixed pre-processing steps did not include dead cell removal nor extensive cell washes, which may attribute to this difference. No reads mapping to ribosomal protein genes were detected in the10x Fixed libraries as the there are no probes for ribosomal protein genes included in the design. Despite our experimental design involving the use of five drug treatments aimed at inducing cell death, we did not observe an increased ratio of mitochondrial genes when comparing drug-treated cells to control (DMSO) cells within each method (Supplementary Figure S5).

The relationship between the number of UMIs sequenced and genes detected per cell differed by method (Figure 2D), with Parse achieving the cells with the greatest total number of genes detected (maximum 15594 genes in a single cell). However, taking the average number of genes detected per cell across the all the sequenced cells as a measure of sensitivity, the 10x Fixed approach achieved the highest sensitivity (mean ± SD = 6572 ± 1064 genes/cell), followed by MULTI-seq (mean ± SD = 4199 ± 881), while Parse showed the lowest sensitivity in detecting genes, as well as greatest variability between the cells (mean ± SD = 3690 ± 2140) (Figure 2E). Yet, the correlation of average expression values was high (ρ = 0.991–0.997) between sample replicates for Parse (Supplementary Figure S6). To further examine the agreement between methods, we performed pairwise comparisons of the average expression of each gene by method for the DMSO control cells (Figure 2F). Overall, moderate correlations were observed between the methods across all of the treatments (ρ range = 0.60–0.80), with MULTI-seq and 10x Fixed achieving the highest concordances (Supplementary Table S1).

Down sampling

To fairly compare the efficiency of mRNA capture between protocols, we downsampled the sequencing reads per cell to a common depth and stepwise-reduced fractions (80000–2500 reads per cell) for each of the three scRNA-seq methods, revealing large differences in sensitivity (Figure 3A, B). The 10x Fixed approach consistently detected more genes per cell than the other methods, especially at lower sequencing depths. Once the sequencing depth exceeded 60000 reads per cell the number of genes detected per cell was similar between the 10x Fixed and Parse libraries, albeit in a significantly smaller proportion of cells in the Parse libraries. This was further visualized in UMAP plots where data the original dataset was compared to data downsampled to 20000 and 5000 reads per cell. The MULTI-seq and 10x Fixed approaches displayed a more even distribution of genes detected per cell and the 10x Fixed approach detected the most genes across the cell population (Figure 3C–E). To further compare between the libraries, we sub-sampled the DMSO cells from the three scRNA-seq datasets (20k and 5k reads/cell) to equal number of cells per library to avoid bias in differing number of cells and subsequently calculated the average expression, variance, and dropout rate for the geneset present in 25% of the cell subset (Figure 3F). In this comparison, 10x Fixed displayed highest average expression and standard deviation, as well as lowest dropout rates across the downsampled datasets (27.0–65.4%), while MULTI-seq had the lowest overall expression levels as well as the highest median dropout rates (53.0–80.4%). In summary, 10x Fixed detects the common sets of genes in more cells than the other two methods.

Figure 3.

Gene detection and drop out in downsampled data. (A) The median number of genes detected per cell (y-axis) by downsampled sequencing depth in reads per cell (x-axis). The single-cell RNA-sequencing (scRNA-seq) library preparation method is indicated by color according the legend. (B) The proportion of cells with sufficient sequencing reads to be compared at each downsampled read depth, color-coded by scRNA-seq method. (C-E) UMAPs by scRNA-seq. Each cell is color coded by the number of genes detected in each cell, according to the legend on the right side of each panel. The number of cells in each UMAP panel is indicated in the brackets. The UMAPs in panel C are generated from non-downsampled data. The UMAPs in panel D are generated from downsampled data to 20000 reads per cell. The UMAPs in panel E are generated from downsampled data to 5000 reads per cell. (F, G) The average expression (F) and standard deviation (G) calculated using across 7766 genes detected in 25% of 1500 sub-sampled DMSO-treated cells (500 cells per scRNA-seq method) downsampled to 20k reads/cell and the 3627 genes detected in the same dataset, downsampled to 5k reads per cell. (H) The dropout rate is calculated by method and gene as the proportion of the 500 cells in which the gene is not detected (zero counts).

Drug-specific transcriptional profiles

Next, we examined the influence of specific drug conditions on gene expression profiles of the single cells recovered by the three methods. Unsupervised clustering of co-normalized data revealed that these sequencing approaches captured orthogonal aspects of the data in a highly complementary way (Figure 4A). The cells exposed to the drugs revealed a largely overlapping pattern with DMSO treated control cells, except for the cells exposed to fludarabine, where the majority of the fludarabine treated cells displayed a different transcriptional state and were confined to one cluster (Figure 4B). This pattern was evident in each of the three scRNA-seq methods (Supplementary Figures S7-9). Each cell was further assigned into G1, S or G2M phases of the cell cycle based on cell cycle scoring in Seurat (23) (Figure 4C, Supplementary Figures S7-9), which revealed that the majority of fludarabine treated cells displayed a transcriptional profile related to the S-phase of the cell cycle (Figure 4D). We examined the cell cycle with orthogonal analysis with DNA stain intensity, which confirmed that a large proportion of the cells appear to be in S-phase, with few cells in G2/M, further supporting that the fludarabine treated cells were in cell cycle arrest (Supplementary Figure S10).

Figure 4.

Global gene expression in the dataset. (A–C) UMAPs of the merged dataset of 27631 cells from the Parse, MULTI-seq and 10x Fixed libraries. Each dot in the graph represents a single cell. The UMAPs are color coded by library preparation method (A), drug condition (B), and by cell cycle phase marker gene expression (C). The color key is indicated at the bottom of each respective panel. (D) Stacked barplots indicating the proportions of cells annotated to G1-, S- or G2/M-phases of the cell cycle by scRNA-seq method.

Ability of scRNA-seq methods to detect drug-specific differential gene expression

To further examine the underlying transcriptional profiles and how they differed between methods and treatment conditions, differential gene expression analyses were performed using MAST (28), by scRNA-seq method, comparing cells between drug-treated condition and DMSO controls. We used a cut-off (absolute log fold change > 0.50, adjusted P< 0.01) for determining differentially expressed genes (DEGs). The resulting DEGs of were compared for each of the treatments (Figure 5A). As expected due to the glucocorticoid resistance of the Reh cell line (35), few or no DEGs were identified in the dexamethasone treated cells (Supplementary Table S2). A small number of DEGs were detected in the 5-azacytadine, Imatinib, or XRP44X treated cells and only a few of these were detected by more than one method (Supplementary Table S3-S5). The largest number of DEGs were found in the fludarabine treated cells: 1245 by 10x Fixed, 331 by MULTI-seq, and 156 by Parse (Supplementary Table S6), with 42 DEGs detected across all methods. Gene Set Enrichment Analysis (31) of fludarabine DEGs suggested a high concordance of Hallmark gene sets enriched in differentially expressed genes called by the scRNA-seq methods (Table 2, Supplementary Table S7). Examples among the top enriched gene sets detected across the three scRNA-seq methods include upregulation of p53 pathway signaling members (e.g. CDKN1A, DDB2, ZMAT3; FDR q-value = Parse 1.08 × 10⁻⁴; MULTI-seq 3.87 × 10⁻¹⁸; 10x Fixed 9.00 × 10⁻²¹) and deregulation of genes encoding cell cycle related targets of E2F transcription factors (e.g. CDKN1A, ASF1B; FDR q-value = Parse 8.92 × 10⁻⁴; MULTI-seq 3.58 × 10⁻¹¹; 10x Fixed 8.78 × 10⁻⁸, Figure 5B). Among the genes involved in these gene sets, we particularly detected p53-dependent mediators of apoptosis, as well as genes involved in DNA damage response and oxidative stress. We additionally observed a significant downregulation of MYC targets, but only in the MULTI-seq and 10x Fixed gene sets (e.g. MYC, PA2G4, NME1, SRSF2/3, HSPD1; FDR q-value = Parse NS; MULTI-seq 7.98 × 10⁻³⁰; 10x Fixed 2.28 × 10⁻⁴⁴). Altogether, these results support the ability of the tested scRNA-seq methods to detect biologically relevant pathways associated with drug perturbations in ALL, with 10x Fixed and MULTI-seq achieving a higher level of concordance in DEGs and pathways.

Figure 5.

Differential gene expression analyses. (A) Upset plots show the number of genes with significant (P< 0.01) altered expression (absolute log fold change > 0.50) in the cells surviving each of the five treatment conditions compared against the DMSO-treated control cells. To the left of each plot the numbers of significantly differentially expressed genes are indicated and the upper bars indicate the numbers of genes common (or unique) between the methods. (B) Heatmap of differentially expressed genes in the enriched pathways in the fludarabine-treated cells detected by the scRNA-seq methods. The cells are ordered along the x-axis by scRNA-seq method and then by control (DMSO, red) and fludarabine treated (pink). The differential expressed genes are ordered along the y-axis and grouped according to the pathway. The color key for the gene expression levels is indicated at the bottom of the panel. Purple indicates low expression, yellow high expression, and grey is missing data.

Table 2.

Gene set enrichment analysis for hallmark gene sets enriched in at least two scRNA-seq methods

			Parse				MULTI-seq				10x Fixed
Hallmark Gene Set Name	Description	# Genes in Gene Set (K)	# Genes in Overlap (k)	k/K	P-value	FDR q-value	# Genes in Overlap (k)	k/K	P-value	FDR q-value	# Genes in Overlap (k)	k/K	P-value	FDR q-value
P53 Pathway	Genes involved in p53 pathways and networks	200	7	0.04	6.49e-6	1.08e-4	22	0.11	2.32e-19	3.87e-18	28	0.14	3.60e-22	9.00e-21
E2F Targets	Genes encoding cell cycle related targets of E2F transcription factors	200	8	0.04	5.12e-7	2.56e-5	25	0.13	3.12e-23	7.81e-22	27	0.14	5.34e-21	8.90e-20
G2M Checkpoint	Genes involved in the G2/M checkpoint, as in progression through the cell division cycle	200	6	0.03	7.14e-5	8.92e-4	16	0.08	2.86e-12	3.58e-11	15	0.08	1.23e-08	8.78e-08
Allograft rejection	Genes up-regulated during transplant rejection	200	4	0.02	5.17e-3	1.99e-2	13	0.07	4.13e-9	2.95e-8	12	0.06	3.72e-06	1.86e-05
MYC Targets V1	A subgroup of genes regulated by MYC - version 1 (v1)	200	NA	NA	NA	NA	31	0.16	1.59e-31	7.96e-30	46	0.23	4.56e-46	2.28e-44
MYC Targets V2	A subgroup of genes regulated by MYC - version 2 (v2)	58	NA	NA	NA	NA	9	0.16	4.71e-10	3.93e-9	16	0.28	3.28e-18	4.10e-17
MTORC1 signaling	Genes up-regulated through activation of mTORC1 complex	200	NA	NA	NA	NA	14	0.07	3.94e-10	3.93e-9	20	0.10	2.28e-13	2.28e-12
PI3K, AKT, MTOR. signaling	Genes up-regulated by activation of the PI3K/AKT/mTOR pathway	105	6	0.06	1.82e-6	4.55e-5	6	0.06	1.38e-4	8.65e-4	NA	NA	NA	NA
Apoptosis	Genes mediating programmed cell death (apoptosis) by activation of caspases	161	5	0.03	2.48e-4	2.07e-3	7	0.04	2.13e-4	1.07e-3	NA	NA	NA	NA

Discussion

In this study, we conducted a comprehensive assessment of three multiplexed experimental approaches to investigate transcriptional changes in individual cells during cytotoxic ex vivo drug screening (FMCA). Our aim was to compare the performance of three single-cell scRNA-seq methods using reference leukemia (Reh) cells treated with various cytotoxic drugs. By evaluating a wide range of performance metrics, we were able to identify technical strengths and limitations intrinsic to each approach. Our findings demonstrate, through the example of fludarabine-treated cells, that any of the scRNA-seq methodologies could robustly discern mechanistic underpinnings of the drug in ALL cells.

Each of the three scRNA-seq methods employed in our study generated high-quality data for single-cell gene expression profiling. Nevertheless, it is essential to recognize that these methods come with their inherent advantages and limitations. The Parse and 10x Fixed methods offered the advantage of cell fixation, allowing for potential sample storage and collection at multiple time points. This feature may be particularly beneficial for studies involving delicate or easily disrupted cell types (42). On the other hand, the MULTI-seq protocol has the greatest potential for achieving ultra-high multiplexing (22) and thus lower costs, although it requires a larger number of input cells and is more time-consuming compared to the other methods. It is worth noting that minimizing the inclusion of ambient RNA from dying cells during multiplexing is crucial for obtaining optimal results when using LMO or similar antibody-based cell hashing protocols, as suggested previously (43,44). In our hands, the additional washing steps and the need for performing a dead cell removal step after labeling with LMOs contributed to the longer processing time until mRNA capture, which made the MULTI-seq approach the most laborious of the three methods tested. At present, the 10x Fixed protocol is most limited in terms of number of barcodes available (n = 16), however it is feasible to process several pooled reactions in parallel. The probe-based design of the 10x Fixed method exhibited superior performance in the quality metrics we assessed in this study. However, it is crucial to note that this method does not generate cDNA during the library preparation, which restricts its utility when aiming to detect isoforms in conjunction with long-read sequencing (45), or for analysis of SNVs for genotyping or mutation detection to deconvolute pooled samples (46). Future endeavors focused on i) increasing the multiplexing capabilities of scRNA-seq experiments, ii) streamlining the laboratory procedure, and iii) efficiently fixing cells, will be needed to bring down the cost of fully implementing single-cell read outs into high-throughput ex vivo drug screening.

In terms of sequencing efficiency and gene detection sensitivity the 10x Fixed approach outperformed the other two methods. A higher number of genes per cell (mean = 6.5k) were consistently detected with 10x Fixed, which is in line with a previous study utilizing the same method (https://doi.org/10.1101/2023.04.25.538273). Our differential expression analysis further underscored findings revealing more differential expressed genes using the 10x Fixed approach. Despite that the Parse approach achieved the cells with the highest number of UMIs sequenced and genes detected, this was observed in only a small fraction of individual cells.

Our study revealed that treatment with fludarabine elicited the most pronounced transcriptional changes, exhibiting a large number of differentially expressed genes that converge on p53 pathway signaling and Myc transcription factor activity. Fludarabine treatment is known to induce a p53-dependent transcriptional responses in chronic lymphocytic leukemia (45,46), and our study highlights a similar mechanism of action in an additional hematological malignancy model. Although we observed only minimal overlap in differentially expressed genes between the scRNA-seq method, of the genes identified, DDB2, was consistently detected by all three methods and is known to be induced by p53. Another gene, CDKN1A, which encodes the protein p21, was also upregulated and is involved in cell cycle arrest (47) and potentially inhibits apoptosis (48). CDKN1A is highly expressed in primary leukemic cells at relapse and, when overexpressed in Reh cells, has been associated to an increased resistance to cytarabine (12), an antimetabolite similar to fludarabine. We additionally detected a consistent downregulation of Myc transcriptional targets in fludarabine-treated cells with the MULTI-seq and 10x Fixed approaches. A recent drug repurposing effort found fludarabine phosphate to suppress MYCN signaling in neuroendocrine prostate cancers (49). These findings motivate further exploration of MYC transcriptional activity as a biomarker of response to fludarabine and other DNA synthesis inhibitors.

Contrary to our expectations, our study did not reveal analogous transcriptional perturbations across the remaining drugs (36). Several factors may contribute to this observation. Firstly, the combination of a single high-dose treatment and 72-h incubation period may have surpassed the therapeutic window for some drugs, potentially causing us to miss time-sensitive transcriptional insights. Secondly, the administration of drugs at high concentrations may have led to polypharmacological responses rather than a targeted effect on the primary mechanism of action. Thirdly, our choice of stringent statistical thresholds, aimed at enhancing comparability between scRNA-seq methods, may have inadvertently masked subtle drug-induced transcriptional responses. To address these potential limitations, future experiments could consider implementing lower drug doses to capture more nuanced responses, adopting time series analyses with shorter incubation durations to better capture dynamic transcriptional changes, and employing a more flexible approach to statistical thresholds. These modifications may allow for a more comprehensive exploration of drug-induced transcriptional dynamics.

While our study comprehensively assessed three multiplexed scRNA-seq methods for the purposed of ex vivo drug screening, there are limitations to be discussed. In conducting this comparative analysis, we employed a conventional commercially available B-cell leukemia model (34). Although cell line biology is distinct from primary samples in many aspects, these data can potentially serve as a valuable reference and resource, facilitating future integration and comparison of results to make biological inferences by the wider scientific community (50). Primary leukemia specimens might conceivably result in higher variability, particularly in terms of increased cellular heterogeneity, mRNA expression levels, and viability, which could have added important information on the performance of the different scRNA-seq methods. However, our cell line strategy not only served to mitigate potential sources of experimental variation across the scRNA-assays, but also permitted open data sharing and dissemination. We also acknowledge the presence of other techniques that facilitate the multiplexing of cells originating from multiple conditions in a single RNA-seq library. Notable examples include transient transfection of short barcoding oligos (SBOs) (51), antibody-based cell hashing (52), the 10x Genomics 3′ Cellplex approach, among others. Consequently, our study does not cover all possible approaches in this field. Moreover, we recognize that our analyses were conducted leveraging the first version of the Parse Biosciences Evercode Whole Transcriptome method. It is plausible that the improved second iteration (V2) of the kit may yield different metrics and performance outcomes (https://doi.org/10.1101/2022.08.27.505512). Finally, our primary aim was to identify drugs featuring dose-response curves where ∼50% of cells persisted, enabling the exploration of drug-resistant cellular states via scRNA-seq. This selection aimed to achieve balance between capturing substantial drug effects while ensuring a sufficient number of viable cells. Although we recognize the strengths of this approach in examining cellular phenotypes post-extended drug exposure, it is worth noting that the initial FMCA drug screening might have overlooked numerous additional interesting drug-responses because the high drug concentrations resulted in complete cell clearance.

Overall, the outcomes of this study underscore both the technical disparities inherent to scRNA-seq methods and the remarkable concurrence in the biological processes underpinning drug response. This level of concurrence in the scRNA-seq data is important for extending future possibilities to explore the mechanism of action for other therapeutics. In conclusion, we show that combining ex vivo drug screening with the high resolution of scRNA-seq is a feasible approach to uncovering mechanisms of action for drugs, bridging struggles associated with cellular heterogeneity in bulk samples, and opens an exciting opportunity for generating new biological insight in future FPM studies.

Data availability

The scRNA-seq data was deposited at the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) under accession number GSE229617.

Supplementary data

Supplementary Data are available at NARGAB Online.

Acknowledgements

Sequencing was performed by the SciLifeLab National Genomics Infrastructure (NGI) SNP&SEQ unit at Uppsala University. NGI is funded by SciLifeLab, the Knut and Alice Wallenberg Foundation, and the Swedish Research Council through grant agreement no. 2019-0222. The computation analyses were enabled by resources in project SNIC 2022/22-63 provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at UPPMAX, funded by the Swedish Research Council through grant agreement no. 2022-06725. We thank Katarina Tegnér, Marie Söderberg, Anna Haukkala, Magnus Lindell, Elin Övernäs and Sara Ekberg from NGI for assistance with sequencing and data processing.

Author contributions: Henrik Gezelius: Conceptualization, Methodology, Investigation, Formal analysis, Validation, Visualization, Writing-original draft. Anna Pia Enblad: Conceptualization, Formal Analysis, Writing-review & editing. Anders Lundmark: Formal analysis, Data curation, Visualization. Martin Åberg: Investigation, Validation, Writing-review & editing. Kristin Blom: Methodology. Jakob Rudfeldt: Investigation, Validation. Arja Harila: Conceptualization. Amanda Raine: Conceptualization. Verónica Rendo: Formal Analysis, Writing-review & editing. Merja Heinäniemi: Writing-review & editing. Claes Andersson: Methodology, Recourses, Formal analysis, Writing-review & editing. Jessica Nordlund: Conceptualization, Recourses, Funding acquisition, Visualization, Writing-original draft.

Funding

The Swedish Childhood Cancer Fund [PR2019-0046, TJ2020-0039, PR2022-0082]; Swedish Research Council [2019-01976]; Göran Gustafsson Foundation; European Union's Horizon 2020 research and innovation program [824110 EASI-Genomics]. Funding for open access charge: Uppsala University.

Conflict of interest statement. None declared.