News

Jul 27, 2023

Combi

Nature Communications volume

Nature Communications volume 13, Article number: 4450 (2022) Cite this article

9295 Accesses

5 Citations

127 Altmetric

Metrics details

Anti-cancer therapies often exhibit only short-term effects. Tumors typically develop drug resistance causing relapses that might be tackled with drug combinations. Identification of the right combination is challenging and would benefit from high-content, high-throughput combinatorial screens directly on patient biopsies. However, such screens require a large amount of material, normally not available from patients. To address these challenges, we present a scalable microfluidic workflow, called Combi-Seq, to screen hundreds of drug combinations in picoliter-size droplets using transcriptome changes as a readout for drug effects. We devise a deterministic combinatorial DNA barcoding approach to encode treatment conditions, enabling the gene expression-based readout of drug effects in a highly multiplexed fashion. We apply Combi-Seq to screen the effect of 420 drug combinations on the transcriptome of K562 cells using only ~250 single cell droplets per condition, to successfully predict synergistic and antagonistic drug pairs, as well as their pathway activities.

Despite major progress over the last decades, cancer remains a major cause of death. Our increased molecular understanding of the molecular basis of cancer has led to the development of targeted therapies. These therapies have so far provided limited efficacy and only in a small subset of patients1, despite major efforts to characterize patients genomically to find response biomarkers.

An approach that holds the promise to improve this situation is to complement large-genomic profiling in basal conditions with measurements after perturbing cancer cells with drugs2. While many approaches can be used to perform drug screenings, they are often low in throughput3, cost and time extensive4, and/or require large amounts of cells5, which together strongly restricts the number of potential drugs that can be screened per tumor biopsy. This limitation gets more pronounced when considering drug combinations due to the sheer number of potential combinations, which increases exponentially with the number of tested drugs.

Due to limited screening capacities, computational approaches to model drug–drug interactions have been developed6. While models on drug efficacies improved over the past years by an increase in available data resources, predictions on drug responses remain challenging and limited to well-characterized systems such as cell lines, thereby limiting their translatability into clinics. Among the different data types, gene expression states of cells were shown to be highly predictive of drug response7. Additionally, data repositories of drug-induced transcriptional changes, such as LINCS8, have proven to be a valuable resource. While there are already perturbation screening platforms available in plates for bulk9,10 and single-cell11,12 transcriptomics, they usually require large numbers of cells per tested condition, and they have not been used for screening drug combinations. Therefore, integrating transcriptomic readouts into a miniaturized combinatorial drug screening platform with the potential to screen tumor biopsies will enable more relevant predictions and increase our understanding of the mode of action of synergistic and antagonistic drug–drug interactions.

Droplet-based microfluidics, which uses picoliter to nanoliter-sized droplets as reaction vessels to perform cellular screens, provides a promising approach to achieve this goal. Due to the miniaturization over several orders of magnitude as compared to conventional plate-based screens, the number of drugs or drug combinations can be massively upscaled while working with low input cell numbers13. We previously demonstrated the first step in this direction by integrating Braille valves into a droplet microfluidic system to generate drug combinations in so-called plugs (~500 nl large droplets) stored sequentially in tubings14. Plugs were used to directly screen 56 combinatorial treatment options on pancreatic tumor biopsies to find the most potent drug pairs using a phenotypic apoptosis readout. While our previous approach provided the first proof of concept in directly screening patient material, the still relatively large volumes of 500 nl limited the number of drug pairs tested. Furthermore, an apoptosis assay provides only a single endpoint readout with limited insights into the drug pairs’ mode of action, which could significantly improve our understanding and the predictability of drug combinations to tackle resistance mechanisms.

To overcome these limitations, we present here a microfluidic platform that allows to perform highly multiplexed screens of hundreds of drug combinations in an emulsion of picoliter-sized droplets. By introducing a deterministic combinatorial barcoding approach, where sets of two barcodes encode drug pairs, we managed to screen all conditions in a highly multiplexed fashion, without the need to keep any spatial order (e.g., wells, plug-sequence). Since the DNA barcodes were designed for whole transcriptome analysis of cells after drug perturbation, we were additionally able to perform massively parallelized gene expression-based profiling of drug combinations. We demonstrated that the presented Combi-Seq approach can be applied to determine the impact of drugs on cell viability and cellular signaling, thus providing a high-throughput approach to discover synergistic drug pairs and decipher their mode of action.

Multiplexed combinatorial drug screens were performed in single-cell droplets by encapsulating drugs together with DNA barcode fragments (Fig. 1a). Each pairwise drug combination was encoded by a unique combination of two DNA barcode fragments, which together provided a priming site for reverse transcription (poly-dT) and PCR. After off-chip incubation of droplets, reagents for cell lysis, barcode fragment ligation, and reverse transcription were added to each droplet by picoinjection15. The ligation of two barcode fragments (BC-RT and BC-PCR) resulted in functional barcodes, encoding pairwise drug combinations (Fig. 1). Since the barcodes were used for the reverse transcription of mRNA released from lysed cells, transcriptomes were barcoded according to drug treatments (Fig. 1b). Subsequently, barcoded cDNA was extracted from the droplets to construct a sequencing library (Supplementary Fig. 1). Finally, sequencing was performed to demultiplex treatment conditions and to analyze their effects on gene expression.

a Overview Combi-Seq workflow: (1) Cells were encapsulated with drug combinations and pairs of barcode fragments encoding drugs. (2) After off-chip incubation, droplets were reinjected into a chip for picoinjection to add reagents for barcode ligation and reverse transcription (RT), enabling barcoding of the transcriptome according to the drug treatment. (3) Upon droplet breakage, pooled libraries for sequencing were generated in which cells from the same treatment group share the same barcode. This facilitated demultiplexing of drug treatments and gene expression-based readouts for pools of 250 cells. b Barcoding strategy applied to encode and decode drug combinations for low input cell pools. Pairs of barcoded PCR primers (BC-PCR) and barcoded poly-dT primers (BC-RT) were joined in a ligation reaction to form functional barcodes, which were used for reverse transcription, thereby encoding a combination of two drugs. By breaking droplets, barcoded cDNA can be recovered and amplified for sequencing. c Microfluidic pipeline used to generate drug combinations in droplets. (1) Braille valves were used to generate sequences of 20 drug-barcode (BC-RT) plugs within the delay tubing. The delay tubing was connected to a drop maker (2) into which cells and drugs-BC-PCR mixtures from multiwell-plates were injected. Finally, injecting the plugs into the drop-maker by opening two oil valves, droplets containing drug pairs with barcodes and cells were generated. M1–M3 indicate positions where fluorescence signals were acquired. d Generation of drug-barcode plugs in the delay tubing: (1) Plugs spaced out by oil were produced by sequentially opening the corresponding valves, (2) resulting in a sequence of 20 drug-barcode plugs. (3) By opening two oil valves, the plugs in the delay tubing were injected into the droplet-maker chip. e Droplet production from a cell suspension, drug-barcode plugs, and drug-barcode mixtures injected by the autosampler, resulting in the co-encapsulation of cells with drug-barcode combinations. f Scheme illustrating the generation of drug combinations from 20 drugs from the valve module with drugs from a 96-well plate. Each sequence of 20 drug-barcode (BC-RT) plugs was combined with drug-barcode (BC-PCR) mixes from one well.

In order to generate drug combinations in picoliter-sized droplets, we synchronized the Braille valve system and autosampler-based injection of drugs into a droplet-maker chip (Fig. 1c). In addition, cell suspensions were injected into the droplet-maker chip at a density of 0.1 cells per droplet volume, to obtain droplets containing single cells (Supplementary Movie 1). The autosampler (Dionex) was loaded with a 96-well plate, with each well containing a single drug together with the corresponding barcoded primer fragment (BC-PCR) and a marker dye enabling to monitor later mixing steps. Drugs were consecutively aspirated and injected into the droplet-maker. The time window between two samples from the autosampler (~3 min) was used to generate a sequence of 20 chemically-distinct plugs, each containing unique pairs of two drugs and two barcode fragments (BC-RT and BC-PCR), by injecting secondary drugs and barcodes (BC-RT) into a separate Braille valve chip (Fig. 1c and Supplementary Fig. 2a). In particular, each compound valve was opened sequentially and fluorinated oil was injected in between so that drug-barcode plugs spaced out by an immiscible oil phase could be injected into a delay tubing (Fig. 1d). Once the delay tubing was filled with a sequence of 20 plugs, two oil valves were opened to inject all plugs into the droplet-maker (~2 min, Supplementary Fig. 2b).

Thereby, drug-barcode plugs from the valve system were combined with the drug-barcode mixtures being injected from the autosampler and encapsulated together with single cells into droplets (Fig. 1e). By repeating this process, hundreds of combinations with specific pairs of barcode fragments were generated (Fig. 1f). It is important to note that scaling up the number of combinations can be achieved by increasing the number of drugs injected from the autosampler.

Synchronization between the autosampler and the valve-based injection of drugs was crucial to ensure that combinations were only generated once the drug injected from the autosampler had reached its plateau concentration. Between each drug coming from the autosampler, a time window with decreasing and increasing concentrations was observed, as shown by the alternating injection of fluorescence dyes (Fig. 2a). This phenomenon is based on Taylor-Aris dispersion of solutes in the continuous, miscible carrier phase (PBS) of the autosampler16.No combinations were generated during that time window, which was rather used to produce compound plugs in the delay tubing of the Braille module. Once the plateau concentration was reached, as indicated by measuring a constant intensity of a fluorescent marker dye, the 20 compound plugs were injected into the droplet-maker and combined with the drug from the autosampler and cells into droplets. The injection of one such plug train took 2 min, resulting in an overall time (plug production and injection) of 15 seconds for generating ~2500 droplets containing cells and one barcoded combinatorial treatment condition. Once all 20 plugs were injected, the autosampler started aspirating the subsequent drug.

a Alternating injection of Alexa-488 or Cascade-Blue from a 96-well plate using the autosampler. Samples being injected in the time window of decreasing and increasing fluorescence signals (grey box) were discarded by ensuring no combinations were generated. Since no drug-barcode plugs from the Braille display were injected during time windows with unstable concentrations (grey box), droplets only contained non-functional barcodes (BC-PCR). Time windows with stable fluorescence signals (blue box) were used to generate drug combinations by the co-injection of 20 plugs generated on the Braille display chip resulting in functional barcodes. b Fluorescence intensities of a plug-sequence from the Braille valves supplemented with either Alexa-488 or Cascade-Blue measured at position M1 (=before combinatorial mixing). The blue overlay connecting (a) and (b) illustrates a time series during which one cycle of 20 plugs (second plug without reference dye) is combined with one sample from the autosampler. Alternating sequences of drugs supplemented with Cascade-Blue or Alexa-488 were used to quantify cross-contamination between specific (e.g., green) fluorescent-positive plugs into the subsequent negative (e.g., blue) sample for both injection modes, separately. c Cross-contamination from green positive plugs into green negative plugs from the Braille display (11 cycles as shown in b). d Cross contaminations of plugs for three different chips: The ratio between fluorescence intensities of a blue or green negative plug and the previous blue or green positive plug was analyzed to quantify the level of cross-contamination between sequential samples (n = 99, n = 80, and n = 171, respectively). e Fluorescence signals of droplets containing combinatorial mixtures. Scatterplot representing the fluorescence intensities measured for droplets generated from only Cascade-Blue labeled plugs and only Alexa-488 injected from the autosampler, measured at the droplet outlet (n = 91,899). f Fluorescence signals of droplets from (e) shown for 180 individual combinations. Colors represent cycles of 20 drug plugs combined with one drug from the autosampler. The boxplots show the median and first and third quartiles as a box, and the whiskers indicate the most extreme data points within 1.5 lengths of the box. M1–M3 indicate positions as shown in Fig. 1c. Source data are provided as a Source Data file.

To ensure droplet contents with marginal cross-contamination, we designed the geometry and delay tubing connectors of the Braille valve drop-maker chips such that no residual drug-barcode mixtures remained in the channels (Supplementary Fig. 3). Before each experiment, we measured a proxy for the level of contamination between plugs from the Braille display. This was done by using drugs supplemented alternatingly with Alexa-488 or Cascade-Blue, resulting in an alternating sequence of blue and green fluorescence peaks (Fig. 2b). Fluorescence intensities of plugs were measured on the droplet-maker and contaminations of drugs/dyes from one plug into the subsequent plug were detected by either green signals in UV peaks or UV signals in green peaks (Fig. 2c).

The ratio between the fluorescence signals from each negative peak (n) in either the green or the UV channel with the previous positive peak (n-1) was used as a proxy to quantify the level of cross contaminations between two drugs (Fig. 2d). Over three different chip setups (Braille valve and drop-maker), we found a mean of 1.5% of contamination in the UV channel and 0.7% in the green channel (Supplementary Table 1), indicating that the described systems can be applied to generate combinations with sufficiently high purity.

In the described microfluidic pipeline, drug combinations were generated by mixing drugs injected from an autosampler and a Braille valve module. To ensure precise and accurate drug concentrations within droplets, both drugs had to be encapsulated at a constant and predefined ratio. We validated the precise mixing of two drugs by supplementing all compounds on the Braille display with Cascade-Blue and all compounds from the autosampler with Alexa-488. We observed one highly dense main population of double blue and green double positive droplets, demonstrating that both compounds were co-encapsulated at a constant ratio (Fig. 2e). Furthermore, we confirmed the stable co-encapsulation of the two dyes for droplets over individual combinations (Fig. 2f). The median fluorescence intensities of individual combinations were highly stable with coefficients of variation (CV) over 180 combinations of 2.9 and 3% for blue and green intensities, respectively. The scattering of droplets around the main population can be explained by a short (<100 ms) flow equilibration phase at the beginning and end of plugs and fluctuation of the droplet trajectory within a focused laser beam (Supplementary Fig. 4 and Supplementary Movie 2). Consequently, we concluded that the injection modes (Braille valve or autosampler) can be robustly synchronized to generate drug combinations in droplets at high precision and purity.

To characterize the microfluidic pipeline and to demonstrate its applicability to perform gene expression-based combinatorial drug screens, we designed a small 4 × 4 drug screen (Table 1). First, we wanted to assess whether the injection mode of drugs from the Braille valves vs. autosampler cause any bias, and therefore loaded the same set of drugs on the Braille valves and autosampler. In case of an injection bias, we would expect to see differences between the same combination generated in reverse order (e.g., Imatinib and Trametinib vs Trametinib and Imatinib). Secondly, we aimed at assessing the impact of the barcoding mode on the gene expression readout. For this purpose, we first encoded treatment conditions such that drugs injected from the Braille valves were supplemented with barcoded BC-RT, whereas drugs from the autosampler were supplemented with BC-PCR. Then we repeated the experiment with BC-RT encoding drugs from autosampler and BC-PCR encoding drugs from the Braille valves, expecting comparable results if the barcoding mode is not impacting the readouts. We used the described pipeline to generate droplets, each containing single human leukemia K562 cells and all pairwise combinations of drugs and the corresponding barcodes, and incubated the emulsion for 12 h. After ligation, the two barcoded primer fragments formed one functional barcode encoding the pairwise drug combination. In order to obtain three replicates, the whole process was performed three times.

Each ligated barcode was used to reverse transcribe the transcriptomes from perturbed cells (Fig. 1a). After data preprocessing (Methods Gene-expression data preprocessing) and initial quality controls (Median read and gene count per sample of 3.47 × 105 and 3229, respectively, Supplementary Fig. 5) we performed dimension reduction using t-distributed stochastic neighbor embedding (t-SNE, Fig. 3a) on the demultiplexed count matrix. To analyze whether some systematic bias arises based on the injection source (autosampler or Braille valves), we analyzed samples according to the autosampler drug (Fig. 3a, top left panel), the Braille valves drug (Fig. 3a, top right panel), the "ordered" drug combination (where we made a distinction between e.g., Imatinib-Trametinib and Trametinib-Imatinib combination), and the "unordered" drug combination (where Imatinib-Trametinib and Trametinib-Imatinib samples were not distinguished). While the injection mode for single drugs from the autosampler (Fig. 3a, top left panel) or the Braille valves (Fig. 3a, top right panel) has only a moderate impact on the clustering of individual data points, their pairwise combinations (Fig. 3a, bottom panels) is the stronger determinant on the cohesion and separation between samples.

a TSNE plots of normalized gene expression data. Samples are color-coded based on Autosampler drug (top left panel), Braille valves drug (top right panel), ordered combination (bottom left panel), and unordered combination (bottom right panel). Color code is labeled for autosampler and Braille valves drugs (top panels). b Silhouette scores of sample clustering based on autosampler/Braille valves drugs and ordered / unordered combinations. Silhouettes scores are compared to random distributions (color code) created by permuting sample labels. c Pathway activity heatmap of samples. PROGENy pathway activities were calculated for each sample (z-scores of pathway activities, color code), and the pathway activity matrix was hierarchically clustered. Drugs of combinations are color-coded (light-green: autosampler drug, blue: braille valves drug, cyan: Same drug from autosampler and Braille valves). d Drug-induced pathway activity changes. A linear model (pathway_activity ~YM155 + Imatinib + Trametinib) was fitted for each pathway, and the linear model coefficients (color code) for each drug is plotted as a heatmap. e Drug-induced MAPK activity changes. MAPK activity (y-axis) grouped based on autosampler drug (x-axis) and Braille valves drug (color code), n = 3 biological independent experiments for all samples, except for YM155_Imatinib and DMSO_Imatinib n = 2. The boxplots show the median and first and third quartiles as a box, and the whiskers indicate the most extreme data points within 1.5 lengths of the box. f Correlation of MAPK pathway activities between two 4 × 4 screens, performed with swapped barcode assignments (r = 0.637, p = 0.01, and R2 = 0.406), the shaded area shows a 95% confidence interval of the regression estimate. Source data are provided as a Source Data file.

To further quantify the extent of sample clustering based on injection source for ordered and unordered combinations, we performed silhouette analysis (Fig. 3b, and Methods Clustering-based quality control of gene expression data). As the distribution of silhouette scores are dependent on the number of clusters (4 for drugs, 16 for unordered, and 10 for ordered combinations), we compared the silhouette scores of clustering to random distributions created by permuting the sample labels. The silhouette scores for single drugs and combinations were significantly higher (p values <0.01) than the background distributions, showing that samples cluster together based on the used drugs and combinations. Consequently, also the barcoding mode for single drugs injected from the Braille valves or the autosampler encoded either with barcoded RT or PCR primers do not introduce a bias, since their impact on clustering and hence gene expression, is indistinguishable. Contrary to this, pairwise combinations were driven by clustering of the samples, showing that both drugs were detected together in an unbiased way. We also performed hierarchical clustering using the 100 most highly expressed genes across samples, which also showed drug and combination-based clustering of the samples (Supplementary Fig. 6). Additionally, treated samples were generally well separated from DMSO controls in a principal component analysis (PCA, Supplementary Fig. 7). To further demonstrate that our experimental pipeline does not introduce significant technical biases, we performed the small 4 × 4 screen with swapped barcodes (Braille valves drugs supplemented with BC-PCR and autosampler drugs supplemented with BC-RT). We observed similar quality (Supplementary Fig. 8) and clustering of samples based on tSNE and top 100 expressed genes (Supplementary Figs. 9a, b, 10).

To further analyze the gene expression signatures of cells treated with different combinations, we calculated pathway activity changes for each sample, using the PROGENy method17,18,19. PROGENy calculates pathway activities from gene expression data for 14 cancer-related pathways. Hierarchical clustering of samples based on pathway activities (Fig. 3c) also showed the drug and combination-based clustering. We observed two main clusters, one corresponding to combinations including YM155, while the other was dominated by Trametinib treated samples. Analyzing the associations between pathways activity changes and drugs (Fig. 3d), we found a decreased activity of the Hypoxia pathway in all YM155 treated samples, while all Trametinib (MAPK inhibitor) treated samples showed strong inactivation of MAPK (Fig. 3e, p value from a linear model: 0.03), and related EGFR pathways. Moreover, when correlating the MAPK pathway activity scores from the small 4 × 4 screen (Fig. 3e) and the swapped 4 × 4 screen (Supplementary Fig. 9e), we found a significant correlation (r = 0.637, p = 0.01), demonstrating reproducible detection of pathway activities (Fig. 3f). The pathway analysis suggests that the observed gene expression changes correspond to the known mechanism of action of the used drugs, which we further delineate in the discussion below. In summary, the results support the use of our screening method to analyze combination-induced gene expression changes in a high-throughput manner, enabling the characterization of drug responses in much greater detail as compared to phenotypic assays used previously14.

In order to assess the robustness of the Combi-Seq pipeline, we compared pathway activities from the 4 × 4 Combi-Seq drug screen in droplets to the activities from the same screen performed in a multiwell plate (median read count 1.07 × 106). Hierarchical clustering of correlations between pathway activities from bulk (Supplementary Fig. 11) and droplet (Fig. 3c) data resulted in two main clusters driven by positive correlations in pathway activities: For both formats, one cluster was formed by positive correlations between YM155 treated cells and the other cluster was driven by positive correlations of Imatinib treatments (Fig. 4a). Hence, we could demonstrate that low input droplet-based drug screens using shallow sequencing as a readout captures the main features of bulk high-input deep sequencing data at the functional level of pathway activities. Furthermore, we compared the Combi-Seq barcoding strategy with conventional polyA-based mRNA capturing (PolyA) in the bulk format. Gene expression between Combi-Seq and PolyA data showed correlation coefficients between 0.576 and 0.646 (Supplementary Fig. 12), comparable to the results of previous comparisons of different RNA-Seq sequencing approaches (e.g. Prime-Seq vs. True-Seq: R2 = 0.64)20, suggesting that the Combi-Seq approach does not introduce major biases during library preparation steps.

a Heatmap of correlations between pathway activities of bulk and droplet Combi-Seq data: PROGENy pathway activities were determined for each sample and correlations for matched samples using either data obtained in microtiter plates (rows) or droplets (columns) were calculated (color code blue to red) and used to perform hierarchical clustering. Drugs of combinations are color-coded (light green: Autosampler drug, blue: Braille valves drug, cyan: Same drug from autosampler and Braille valves). b Accuracy of spike-in detection as a correlation between ERCC input concentration and measured transcripts per million (TPM) for 250 cells per sample, shaded area shows 95% confidence interval of the regression estimate. (c) Summary of correlation coefficients for 125 (R = 0.63 and R2 = 0.40), 250 (R = 0.65 and R2 = 0.42) and 500 cells (R = 0.7 and R2 = 0.49) per sample, n = 5 biological independent experiments, data are presented as mean values ± 95% confidence interval. d Sensitivity as the detection probability for spike-in molecules for 250 cells per sample. A spike-in was considered as detected when the probability reached 50%. e Summary of sensitivities for different amounts (125, 250, and 500) of cells per sample, n = 5 biological independent experiments, data were presented as mean values ± 95% confidence interval. Source data are provided as a Source Data file.

To test the sensitivity and accuracy of the Combi-Seq approach, we furthermore made use of a library of 92 different DNA sequences from the External RNA Controls Consortium (ERCC). We spiked in ERCC molecules into droplets and analyzed the accuracy (Fig. 4b, correlation between ERCC molecule concentration and measured expression) and sensitivity (Fig. 4d, the threshold for ERCC molecule detection). We observed an increasing accuracy (Fig. 4c, 0.63, 0.65, 0.7 Pearson's correlation for 125, 250, and 500 input cells, respectively) and sensitivity (Fig. 4e, 2.55, 2.06, and 1.48 log10(attomoles/ul) detection limit for 125, 250, and 500 input cells, respectively) of Combi-Seq method with increased input material (number of cells per droplet, where the amount of added ERCC spike-ins was increased proportionally to the number of cells). These results are within the range of other low input RNA-Seq approaches21.

Based on the promising results of the 4 × 4 drug combination experiment, we performed a high-throughput screen, using a total of 420 different combinatorial treatment conditions. For this study, we purposely selected drugs with low logD values (Supplementary Fig. 13 and Supplementary Data). This way, undesired exchange of typically hydrophobic compounds can be minimized22. Additionally, we performed a series of systematic experiments to test if drug exchange impacts our readout. In particular, we mixed droplets containing only one drug and cells (treated) with droplets containing only DMSO and cells (DMSO controls) and incubated them for 24 h at 37 °C (Supplementary Fig. 14a). Additionally, droplets with DMSO and cells were incubated and processed separately as a positive control for an untreated phenotype (separate DMSO controls). Subsequently, cells were lysed and barcodes ligated to perform reverse transcription for whole transcriptome sequencing (Supplementary Fig. 14b). Projections of the data using UMAP, separated the treated from the DMSO control samples, with the majority of DMSO control samples clustering together (Supplementary Fig. 14c). This indicates that over an incubation period of 24 h, the effect of drugs can only be observed for droplets that contained the drug in the first place, but not in droplets that contained only DMSO. Convincingly, the untreated DMSO samples that were incubated with droplets containing drugs, clustered well with the DMSO controls that were handled separately, even for drugs with positive LogD values such as Imatinib. In order to estimate optimal drug concentrations for large-scale combinatorial screens, we generated dose-response curves using K562 cells to determine their 35% growth inhibition values (GR) for each single drug. Compared to the traditional metric, the IC50 value, the GR50 (or, in our case, the GR35) is not sensitive to the number of divisions that occur during the drug treatments, which might vary depending on cell lines and conditions (Supplementary Data)23. Drugs were assigned to the Braille valves or autosampler to achieve a balanced distribution of drugs with high and low GR35 values. Drugs loaded on the Braille valves or autosampler were supplemented with BC-RT or BC-PCR, respectively. We aimed at generating 250 droplets containing a single-cell for each of the 420 treatment conditions and incubated droplets for 12 h at 37 °C before performing picoinjection for cell lysis, barcode ligations, and RT (Fig. 1a). This process was performed three times to obtain replicates.

After initial preprocessing and quality control (median reads and genes per sample of 32892 and 547, respectively, Supplementary Fig. 15), we performed the same dimensionality reduction (Fig. 5a) and silhouette analysis (Fig. 5b), as for the 4 × 4 screen. Again, samples clustered significantly better based on the used combination, than randomly expected (p values of silhouette scores vs. random distribution: <0.01, 0.47, and <0.01 for autosampler drug, Braille valves drug, and combination, respectively). We also found that the low cell numbers did not impact the ability to distinguish treated from non-treated cells based on gene expression (Supplementary Fig. 16a).

a TSNE plots of normalized gene expression data. YM155 (left panel), Blebbistatin (middle panel), and YM155-Blebbistatin combination treated samples are labeled as representative examples. b Silhouette scores of sample clustering based on autosampler/Braille valves drugs and combinations. Silhouettes scores are compared to random distributions (color code) created by permuting sample labels. (c) ROC analysis of drug signature similarity to LINCS-L1000 data. For each drug, a consensus signature was calculated and similarity (Spearman's rank correlation) to corresponding LINCS-L1000 signatures was calculated. The similarity values were used as predicted values for ROC analysis, while true positives were the matched drug pairs between the high-throughput screen and LINCS-L1000. Source data are provided as a Source Data file.

To further investigate whether the gene expression values of the high-throughput screen are biologically meaningful, we compared the obtained gene expression signatures to those available for the same drugs in the public LINCS-L1000 dataset8. As LINCS-L1000 contains only expression signatures of monotherapy drug treatments, we calculated consensus signatures for each drug of our high-throughput screen (Methods Functional genomic analysis of gene expression signatures) and compared these to consensus signatures generated from the LINCS-L1000 database across all available cell lines and concentration doses (note that LINCS-L1000 does not include data obtained directly from K562 cells). For 32 drugs used in our microfluidic screen, corresponding data on LINCS-L1000 was available. To compare signature similarities of these, we calculated Spearman's correlation coefficients for all pairs of drugs across the two datasets. Our ROC analysis showed that signatures of the same drugs from the two screens (true positives) are more similar, than signatures of unrelated drug pairs (Fig. 5c, ROC AUC: 0.59), and this area under the ROC curve is statistically significant compared to a random distribution created by permuting drug labels (Supplementary Fig. 17, p = 0.019).

Taken together, we demonstrated that for the high-throughput combinatorial drug screen, the injection (Braille valve or autosampler) and barcoding mode did not bias the data and that drug signatures correlation between our data and LINCS-L1000 showed significant similarity. Consequently, we concluded that our Combi-Seq approach can be applied to perform combinatorial drug screens at large scale using low input material and shallow RNA-seq as a readout.

As all used drug concentrations were GR35 values, we expected that synergistic combinations could lead to decreased cell viability, while in the case of antagonistic combinations, we expected increased cell viability values. While we did not measure cell viability directly, the CEVIChE algorithm (Methods Functional genomic analysis of gene expression signatures)18 allowed us to infer cell viability changes for all used drug combinations from gene expression data (Supplementary Fig. 18), which worked robustly for low cell numbers and was independent of the percentage of detected mitochondrial genes (Supplementary Figs. 16b, 19). By comparing predicted decreased viabilities between drug combinations and single drug treatments, we determined synergy scores for all 420 combinations (Fig. 6a). We found several clusters of potential synergistic and antagonistic combinations (e.g.: Triciribine-Dacarbazine and Razoxane-Trametinib, respectively). To experimentally validate the predictions on the synergy of the high-throughput drug screen, we performed 5 × 5 dose matrix combinatorial cell viability screens with all possible combinations of Triciribine, YM155, Razoxane, and Doxorubicin with Dacarbazine, Imatinib, and Trametinib in a microtiter plate format.

a Heatmap showing predicted synergy scores that were determined by comparing viabilities from all 420 drug combinations to corresponding single drug treatments. b Correlation between experimental and predicted cell viability (Pearson correlation r = 0.66, p = 0.018, R2 = 0.44). Drug synergy (y-axis) was measured for 12 combinations in a microtiter plate format (color code) and plotted against the predicted synergy value (x-axis), the shaded area shows a 95% confidence interval of the regression estimate. Source data are provided as a Source Data file.

We calculated synergy scores (positive: synergistic, negative: antagonistic) for the tested 12 combinations using the Bliss independence synergy model (Methods Plate-based viability measurements to validate hits from the microfluidic screen)24. Among the 12 validated combinations, antagonism was more common, which is comparable to previously published data (Supplementary Fig. 20). Our 12 measured (plate experiment) and predicted (from gene expression data obtained in the microfluidic system) synergies showed a significant correlation (Fig. 6b, Pearson correlation r = 0.66, p = 0.018, see Supplementary Table 2 for values without Razoxane-Trametinib). The comparison confirmed that the combinations of Trametinib and topoisomerase 2 inhibitors show antagonistic effects, while BCR-ABL inhibition (Imatinib) synergizes with increased induction of apoptosis by inhibiting survivin with YM155, further confirming the discovery potential of our Combi-Seq platform.

In summary, using our combinatorial microfluidic gene-expression platform, we showed that (i) the measured gene expression values cluster based on the chemical perturbation, (ii) the resulting data is in good agreement with public monotherapy perturbation profiles, and (iii) predicted cell viabilities and drug synergies could be validated in a multiwell plate format for selected hits. Taken together, this illustrates how comprehensive information can be gained from gene expression profiles obtained in a highly multiplexed microfluidic format, sequencing only about 250 cells per drug treatment. This should make the workflow particularly interesting for use with very limited material, such as patient samples.

Cancer patient stratification for personalizing treatments with chemotherapeutics and targeted drugs have shown to increase the success of cancer therapies25,26,27. These efforts are largely driven by dissecting the genomic and transcriptional landscape of tumors or cell lines in order to identify traits that explain drug sensitivities28,29. While a variety of identified genomic and/or transcriptomic markers are successfully used in clinics, they are available for only a small subset of tumor types and patients1. Furthermore, many patients often suffer from tumor relapse30, which is largely rooted in intra-tumor heterogeneity31. The relapse is often driven by the surge of a resistance mechanism to the drug that renders the efficacy of single drugs short-lived32. While treatments with drug combinations offer the potential to reduce the risk of drug resistance, their prediction and empirical evaluation remain challenging.

To advance in solving these challenges, we present here a microfluidic pipeline enabling highly multiplexed combinatorial drug screens in single-cell droplets using global transcriptomics as a readout. By integrating deterministic barcoding of treatment conditions, we were able to assess the efficacy of drug combinations by changes in gene expression and gained comprehensive readouts from whole transcriptome sequencing. We applied our pipeline to screen 420 combinatorial treatment conditions in a single tube, illustrating the high level of multiplexing. Based on assay miniaturization in a droplet format, only about 250 cells were needed per tested condition, hence opening a way for personalized screens directly on patient material and drastically increasing the scale at which combinatorial screens can be performed on patient-derived cell lines or organoids and spheroids.

We have designed the microfluidic platform as a modular system in which the Braille display valves allow us to quickly change between injected drugs overcoming the limitations of a slow autosampler-based injection. Since both are combined on the droplet generator, fast and efficient generation of drug combinations becomes feasible and allows the encapsulation of single cells into droplets of high chemical diversity. Since the autosampler used here injects drugs from up to three 96 or even 384-well plates, the number of drug combinations can be further scaled up to a theoretical maximum of 3 × 384 × 20 = 23,040 combinations in a single experiment. What becomes most limiting at that scale are sequencing costs and available material (when e.g., using primary cells) rather than instrument throughput.

We see significant added potential by the possibility to screen such large numbers of drug combinations at the single-cell level: Integrating fluorescence-based droplet sorting upstream of the picoinjection (cell lysis) step could, for example, be used to physically separate and sequence resistant clones for all 420 treatment options in a single experiment (e.g., implementing the phenotypic Caspase-3 assay we used previously)14. This way, one could analyze the difference in their transcriptomic signature as compared to responding cells, opening the way for highly multiplexed studies to reveal new biomarkers for resistance and (chemical) sensitizers to overcome these. As recently demonstrated, single-cell readouts of drug perturbation provide great insights into heterogeneous drug response11. Performing such screens on patient biopsies will allow us to dissect the impact of tumor heterogeneity on drug responses and thereby to define more efficacious drug combinations and additionally gain insights into their potential resistance mechanisms. In order to enable the encapsulation of primary tumor cells (or microscopic fragments) and/or their expansion into spheroids, we aim at integrating support structures for cell adherence and growth in droplets as previously shown33,34.

All generated datasets demonstrate that neither the barcoding approach nor the injection mode biased the gene expression-based readout. Monotherapies from both injection and barcoding modes had similar impacts, while their combinations had the strongest impact on gene expression and was the main driver of the observed clustering. This confirms the highly precise and accurate operation of the presented microfluidic workflow and the specificity of the deterministic barcoding approach. Additionally, we found significant similarities between consensus gene expression signatures of monotherapies from our large-scale screen with drug signatures from the LINCS-L1000 dataset, illustrating a high level of reproducibility. In the pathway activity analysis, we found that the hierarchical clustering was largely driven by the three drugs YM155, Imatinib, and Trametinib, which further supports the detection of drug-specific effects. The two main clusters were driven by Trametinib treatments inhibiting MAPK and EGFR pathway activities, and the opposing effects of YM155 treatments, inducing the upregulations in MAPK and EGFR pathway activities while inhibiting hypoxia and TRAIL-related pathways. While the effects of Trametinib on MAPK are expected35, the effects of survivin inhibition by YM155 are less well understood, due to complex signaling and incomplete knowledge on survinin36. The increased MAPK pathway activity is likely to reflect a counteractive mechanism since survivin expression was described to be regulated by Sp1 and c-Myc activation through the MAPK pathway37, and a higher concentration of the drug target will reduce the drug effect. As survinin expression has been linked to drug resistance in leukemia, combinatorial treatment with YM155 and Trametinib could potentially have a beneficial effect on decreasing the chances for relapse, due to the inhibition of survivin and putative compensatory expression induced by the MAPK pathway. Taken together, these findings show that the described microfluidic pipeline can be applied to disentangle the effects of drug combinations on pathway activities. Such information will be of great impact when screening patient biopsies to identify potential resistance mechanisms and to predict efficacious drug pairs. Analyzing the pathway activities upon perturbations was limited by the number of detected genes, and, therefore, to the small 4 × 4 screens, since these samples were sequenced at higher depth. To detect a comparable number of genes for all 420 drug combinations, a ten times higher coverage would have been necessary. We instead used the large screen to show that shallow sequencing data is sufficient to determine synergistic drug pairs. We mined the dataset with 420 treatment conditions for drug pairs with synergistic or antagonistic effects. We found that combining Trametinib with topoisomerase 2 (Top2) inhibitors (Doxorubicin or Razoxane) has antagonistic effects (Fig. 6b). We hypothesize that the inhibition of the MAPK pathway by Trametinib and the resulting G1 cell cycle arrest, counteracts the DNA damage normally caused by Top2 inhibition when cells enter the subsequent S-phase. Among the predicted and validated synergistic combinations, we found the combination of Triciribine with Dacarbazine and YM155 with Imatinib among the top hits (Fig. 6b). BCR-ABL (the target of Imatinib) has previously been linked to upregulating the expression of survivin (the target of YM155)38. Therefore, survivin antisense silencing was found to reduce the viability and to increase the efficacy of Imatinib treatment in chronic myeloid leukemia (CML) cell lines, as well as in myeloid progenitors from CML patients39,40. By confirming a previously described efficient treatment combination for CML, that has the potential to prevent the onset of Imatinib resistance, we further illustrate the potential of our pipeline in identifying clinically relevant drug pairs. Furthermore, we validated that there is a good correlation between viability scores from the gene expression data and the experimentally validated synergy scores obtained from the plate-based drug screen. These results show that it is possible to use cost-effective low sequencing depth in large transcriptomics screens to discover synergistic drug pairs.

Together with the inferred pathway activities under perturbation, this should not only allow for the identification of synergistic combinations but also gain insights into their mechanisms of action. Compared to our previous single-measurement phenotypic assay platform14, the global transcriptomic readout provides orders of magnitude more data points per sample, while the cell consumption could be reduced further by a factor of about sixfold. The higher content readouts should enable more robust predictions on the best combinatorial treatments and the discovery of new drug sensitizers and biomarkers, and the even smaller needs of material further facilitate the application in the clinic for patient stratifications and treatment prioritization.

All devices used for the valve module were replicated from molds prepared using soft-lithography with AZ-40XT positive photoresist (Microchemicals) according to the manufacturer's instructions. Structures from 25400 dpi photomasks (Selba) were patterned on 4-inch silicon wafers (Siltronix) in a mask aligner (Suess MicroTec MJB3) using light with a wavelength of 375 nm. Structures were covered with a ~1 cm thick layer of PDMS mixed with a curing agent at a 1:10 ratio (Sylgard 186 elastomer kit, Dow Corning Inc) and cured overnight at 65 °C. In addition, we prepared PDMS membranes by mixing PDMS with a curing agent at a 1:10 ratio and distributing it over a transparent sheet using a spin coater at 500 rpm (Laurell WS 650), which were cured overnight at 65 °C. The drug inlet and waste ports of the valve chip were punched using 0.75 mm biopsy punches (Harris Unicore), whereas the plug outlet port was punched horizontally to the outlet channels using a 0.5 mm biopsy punch (Harris Unicore). Chips were bonded to a PDMS membrane using a plasma oven (Diener Femto). We inserted PTFE tubings with an inner diameter of 0.4 mm (Adtech) into the horizontally punched outlet port until the tubing reached the funnel-like structure of the outlet channel. Subsequently, chips were bound to a glass slide to support chip structures with inlets and outlets. In order to prevent surface wetting, channels were treated with Aquapel (PGG Industries) before use. The valve structures of Braille chips were aligned (Supplementary Fig. 2a) on top of the pins of a Braille display (KGS Corporation, Supplementary Fig. 21a) and mounted using a plexiglas holder (Supplementary Fig. 21b). Using our "CombinatorialPlugFluidics" LabVIEWTM 2013 software (all required software can be downloaded from www.epfl.ch/labs/lbmm/downloads). the movement of the pins was controlled, so that opening the collection channel resulted in the closing of the waste channel and vice versa. Closing a channel was achieved by a pin pushing into the elastic PDMS membrane. This could be opened again by moving the pin down. For all experiments, we used 20 syringes (Becton Dickinson) filled with 5 ml drug-barcode solutions and four syringes filled with 5 ml HFE oil (Novec™ 7500, 3 M). These were connected to the inlet ports of the valve chip with PTFE tubing, and fluids were injected at 500 µl h−1 using syringe pumps (Harvard Apparatus). The waste outlets were connected with a piece of PTFE tubing to direct fluids to a waste container.

Drop-maker molds were manufactured from negative photoresist SU-8 2075 as described by the manufacturer (Microchemicals). PDMS containing a 10% (w/w) curing agent was poured over the molds and cured overnight. Inlet ports for cells and HFE were punched vertically using 0.75 or 1 mm biopsy punches for PEEK tubing coming from connecting the autosampler. Inlet ports for compound plugs from the Braille valves were punched horizontally using 0.5 mm punches. Chips were first plasma bound to a PDMS membrane and subsequently to a glass slide. The channel wall hydrophobicity was increased by injecting Aquapel (PGG Industries) into the channels.

To facilitate the injection of drugs from microtiter plates into the drop-maker device, we used a Dionex 3000SL Autosampler aspirating drugs from a 96-well plate. The autosampler was programmed to sequentially aspirate 310 µl of the compound from wells into a 125 µl sample loop. The large excess of aspirated volume accounted for a needle volume of 60 µl and a loop overfill factor of 2. By overfilling the sample loop by twice its volume, we ensured that the remaining compound mixture was not diluted by the carrier fluids that remain in the sample loop after each cycle due to washing. The injection of compounds from the sample loop into the droplet-maker was driven by a syringe pump (Harvard Apparatus) injecting PBS (Thermo Fisher) at 500 µl h−1. After the aspiration of one drug, a delay time started to ensure that each drug got injected and combined with all drugs from the Braille display before the next drug was aspirated.

Random 10 nt long DNA sequences with balanced base distributions were generated using the bgen tool (gear.embl.de). Barcoded PCR primers were functionalized with a 5′-end biotin for purification, followed by a spacer sequence, a common primer sequence, a unique barcode sequence, and a ligation site (Table 2). The reverse complements (RC) were functionalized with a free 5′-end phosphate to enable ligation. Barcoded RT primers comprised a dT(20)-VN sequence, a unique barcode sequence, and a phosphate group at the 5′-end. The RC for this had a ligation site complementary to the ligation site of the PCR primers. A list of all barcode sequences can be found in the supplementary materials. Complementary sequences were annealed at equimolar concentrations by heating mixtures to 95 °C for 10 min in a thermal block (Eppendorf) followed by their cooling to room temperature (RT) for 1 h.

Barcode-drug mixtures for the valve module were prepared by diluting barcoded RT primers in FreeStyle media (Thermo Fisher) to 1 µM. Drugs dissolved in DMSO were added to their corresponding barcode at 2x the final concentrations (see supplementary materials). The barcode-drug mixtures were supplemented with either Cascade-Blue (Thermo Fisher) or Alexa-488 (Thermo Fisher) at 10 µM for monitoring purposes and subsequently aspirated into 5 ml luer-lock syringes (BD) connected with PTFE tubings using 27 G ¾ needles (BD). Barcode-drug mixtures for the autosampler-based injection were prepared in round bottom 96-well plates by diluting barcoded PCR primers to 4 µM in FreeStyle media and the corresponding drugs to 4x the final concentrations. Mixtures were supplemented with Alexa-488 at 10 µM and plates were sealed with adhesive qPCR seals (Thermo Fisher).

K562 cells (ATCC, CCL-243) were cultured in IMDM media (Thermo Fisher) supplemented with 10% FBS (Thermo Fisher) and 1% Penicillin-Streptomycin (Thermo Fisher). On the day of the experiments, cells were washed twice in PBS and resuspended in FreeStyle Media (Thermo Fisher) supplemented with 4% FBS. The concentration of the cell suspension was adjusted to 2 × 106 cells ml−1 and subsequently aspirated into a 3 ml luer-lock syringe (BD).

Syringes containing drug-barcode mixtures and HFE oil were connected to the Braille valve chip as shown in Fig. 1c and all injected at 500 µl h−1. The default mode for all compound valves was to direct fluids to the waste outlets and two HFE oil valves to direct fluids to the outlet tubing. The length of outlet tubing of the Braille valve was adjusted to harbor all 20 compound plugs spaced out with HFE oil and then connected to the Braille inlet on the drop marker chip (Supplementary Fig. 21b). A syringe containing cells were mounted on a pump and connected to the drop-maker and injected at 500 µl h−1. Cell sedimentation was prevented by low-speed rotations of the magnetic disc using the Multi Stirrus™ system (VP Scientific). Autosampler output tubing and HFE carrier phase supplemented with 1% Pico-Surf1 (Sphere Fluidics) were connected via the respective inlets and injected at 500 and 6000 µl h−1, respectively. The droplet-maker chip was mounted on a microscope (Nikon, Eclipse Ti2-E) with an optical setup for measuring fluorescence intensities of compound plugs or droplets14. Lasers with a wavelength of 375 or 488 nm were used to excite dyes and emitted light (450 or 520 nm) was measured using photomultiplier tubes. Lasers were focused at one of the three positions (M1–M3) shown in Fig. 1c to measure fluorescent dyes that were injected together with the drug-barcode mixtures into the drop-maker. At positions M1 and M2, fluorescence signals were recorded using the PlugAcquisition LabViewTM 2014 software. Fluorescence signals from droplets were acquired at position M3 using the DropletAcquisition LabViewTM 2014 software (all required software can be downloaded from www.epfl.ch/labs/lbmm/downloads).

A CSV file containing the Braille valve opening sequence (Table 3) was loaded into the sample on-demand software. The number of cycles was set to 21 since we were combining all 20 drugs from the Braille valves with 21 drugs from a 96-well plate. Once all tubings were connected and fluids were injected into the drop marker (carrier fluid from the autosampler and HFE oil from the Braille valves), plug production and autosampler-based injection were started simultaneously. 20 plugs were produced into the delay tubing (180 s) and subsequently injected into the drop-maker, by opening two HFE oil valves (tot. flow rate of 1000 µl h−1). This ensured a continuous and stable flow rate for plug injections, and, therefore, resulted in a laminar flow of compounds from the Braille valves, autosampler compounds, and cells from which droplets were generated at the flow focusing junction (Movie S1). Droplets of ~800 pl were collected in an Eppendorf tube which was kept on ice. Once 420 combinations were generated (~100 min) the Eppendorf tube was placed in a humidified incubator at 37 °C and a 5% CO2 atmosphere to incubate cells for 12 h.

Chips for picoinjection were produced by replicating SU-8 molds using PDMS with a curing agent as described above. Casts were plasma bound to glass slides and treated with Aquapel. Chips were heated to 95 °C, and first low melting solder and second cables were inserted into the ports for the two electrodes (Supplementary Fig. 22). The chip was mounted on the microscope of a microfluidic station and the power electrode was connected to a high voltage amplifier, while the chip was grounded over the grounding electrode.

After an incubation of 12 h, droplets containing cells, drug combinations, and corresponding DNA barcodes were transferred into a 3 ml syringe and injected through the droplet inlet port into a picoinjection chip (Supplementary Fig. 2). Droplets were flushed at 180 µl h−1 and individual droplets were spaced out by injecting HFE oil with 1% PicoSurf surfactant at 700 to 1000 µl h−1 over the oil inlet. For cell lysis, barcode ligation, and reverse transcription of mRNA released from lysed cells, we pico-injected a one-pot reaction mix containing 0.9% Igepal (Sigma Aldrich), 3x ligation buffer (NEB), 60,000 U/µl T4 Ligase (NEB), 1.5 mM dNTPs (Thermo Fisher), 7.5 µM Template Switching Oligonucleotide (IDT), 12 U/µl Maxima -H reverse transcriptase (Thermo Fisher) and 6000 U/µl NxGen RNase Inhibitor (Lucigen). Flow rates for the reagent mix were adjusted according to the droplet frequency and size in order to inject the equivalent of 1/3 of the final droplet volume (Movie S3). In order to achieve the injection of reagents into the droplets passing by the injector nozzle, we applied a continuous electrical field of 0.1 V using a function generator (Rigol). Picoinjection was performed over ~1 h during which all droplets (injection and collection) were kept on ice. Subsequently, the emulsion was kept at RT for 30 min and then incubated for 90 min at 42 °C.

Drug exchange experiments were performed like the described combinatorial drug screens with small adaptations: The different drugs were only injected from the autosampler while only media supplemented with BC-15 or BC-O and DMSO were injected continuously from the braille inlet, using a normal syringe pump instead of the braille display. Emulsions of all samples were collected in individual wells of 96-well plates, incubated for 24 h at 37 °C in a 5% CO2 atmosphere, pooled, and then processed as described for the combinatorial drug screens.

First, 10 µl of BC-PCR and 5 µl BC-RT, both 20 µM, were added to the respective wells of a 96-well plate, followed by 4 µl from one of the corresponding 50-fold drug stock of Imatinib, Trametinib, or YM155 (Table S4). K562 cells were resuspended in FS-media with 1% (wt/vol) FBS to a concentration of 1050 cells/µl. For each well containing the mixture of drugs and barcodes, 181 µl of this cell suspension were added and then incubated at 37 °C in a 5% CO2 atmosphere for 12 h. To increase the efficiency, cell lysis and ligation were done first, followed by purification and washing steps, before carrying out reverse transcription. In particular, 20 µl of cells were transferred to a well containing 10 µl of a 0.9% Igepal (Sigma Aldrich), 3x ligation buffer (NEB), 60,000 U/µl T4 Ligase (NEB), and 6000 U/µl NxGen RNase Inhibitor (Lucigen) for cell lysis (10 min on ice) and barcode ligation (30 min at RT). About 200 µl of 6x SSC were added to the mixtures before transferring them to 1.5 ml Eppendorf tubes for centrifugation at 2000×g for 5 min. Supernatants were transferred to fresh Eppendorf tubes supplemented with 50 µl of C1 dynabeads at 10 µg µl−1 in 6x SSC buffer and incubated for 15 min at RT on a nutator (VWR). Beads were washed three times in 6x SSC buffer, briefly centrifuged in a benchtop centrifuge, and then resuspended in 80 µl of reverse transcription mix, containing 1x Maxima RT buffer (Thermo Fisher), 4% Ficoll-PM 400 (Sigma Aldrich), 1 mM dNTPs, 2000 U/µl NxGen RNase Inhibitor (Lucigen), 2.5 µM Template Switching Oligonucleotide (IDT) and 4 U/µl Maxima -H reverse transcriptase (Thermo Fisher). Mixtures were incubated under nutation (VWR) for 30 min at RT and for 90 min at 42 °C. Beads were washed twice in TE-SDS (10 mM Tris pH 8.0, 1 mM EDTA, and 0.5% SDS) and twice in nuclease-free water (Thermo Fisher) before cDNA from different conditions was pooled and processed for library preparations as described in the section "Library preparation and sequencing".

K562 cells were cultured as described in section "Preparation of cell suspensions". On the day of the experiments, cells were washed in PBS and resuspended in FreeStyle Media (Thermo Fisher)", supplemented with 1% FBS. For comparing Combi-seq and PolyA-based RNA capture strategies, 250 cells/μl were transferred into an Eppendorf tube and mixed with a 5 μl solution, containing the corresponding drugs or DMSO, together with Combi-Seq or PolyA-Seq barcodes suspended in FreeStyle Media, 1% FBS. After 16 h incubation at 37 °C and 5% CO2 atmosphere, cells were mixed with 15 μl of the solution containing lysis buffer, ligase, and reverse transcriptase as in the droplet experiment (see section "Picoinjection for cell lysis, barcode ligation and reverse transcription" for details).

ERCC spike-ins (Thermo Fisher, Cat No. 4456739) were diluted at 1:100 in RNase-free water. Then 30 μl of the diluted ERCC spike-ins were added to a total of 600 μl of the one-pot reaction mix. This reaction mix was pico-injected at a ratio of 1:3 to a droplet as described in section "Picoinjection for cell lysis, barcode ligation and reverse transcription".

Upon reverse transcription of mRNA, all cDNA was barcoded according to drug treatments, and therefore, we broke the emulsion by adding 0.5 to 1 ml of 1H,1H,2H,2H-Perfluorooctanol (Abcr). The supernatant was transferred into a fresh Eppendorf tube, supplemented with 1x the volume of C1 dynabeads (Thermo Fisher) at 2.5 µg µl−1 in 6x SSC buffer (Thermo Fisher), and incubated at RT for 20 min. Beads were washed 2x in TE-SDS (10 mM Tris pH 8.0, 1 mM EDTA, and 0.5% SDS) and 2x in nuclease-free water (Thermo Fisher). Beads were resuspended at 5 µg µl−1 followed by MseI (NEB) digestion according to the manufacturer's instructions. The supernatant was purified 2x using SPRIselect beads (BD), first at 0.6x and then at 0.8x the volume of cDNA, and finally amplified in KAPA HiFi ready mix (Roche) with 0.8 µM of SMART primer (Supplementary Table 3) over a total of 13 cycles (PCR program in Supplementary Table 4). Products were purified on 0.6x the volume SPRIselect beads and then analyzed using high sensitivity DNA chips on a 2100 Bioanalyzer (Agilent). Fragmentation of cDNA was performed to shorten the fragments and to introduce linker sequences. This was achieved using a Tn5-based tagmentation protocol for 3′ end libraries developed in house41. Fragments were amplified using a P5-SMART primer (Supplementary Table 3) and an i7 indexed P7 adapter primers (Illumina, Supplementary Table 3) at 0.75 µM in KAPA HiFi ready mix (PCR program, Supplementary Table 4). Fragments were purified on 1x the volume of SPRIselect beads, and size distributions were determined using a Bioanalyzer. All replicates were pooled at equimolar ratios and sequenced on a NextSeq 500 (Illumina) machine together with 10% PhiX spike-ins. Paired-end sequencing was performed by sequencing the barcode combination (Read 1, 26 bp) using the sequencing custom primer (Supplementary Table 3) and the mRNA (Read 2, 59 bp).

In boxplots, the center line represents the measured median, and the upper box and lower box hinges correspond to the first and third quartiles. The whiskers extending from the lower and upper hinges of the box represent the 1.5-fold interquartile range. The dots with the lines shown in the violin plot in Fig. 2d correspond to the mean with standard deviation for each of the measured cross contaminations.

We used the SCANPY pipeline42 for gene expression preprocessing and quality control. Samples with low gene count and a high ratio of mitochondrial genes (>15 percent) and genes with a high dropout rate were filtered out. Read counts were normalized based on sequencing depth and z-score transformed. The batch effect (replicates) was removed by using the combat function of SCANPY. For dimension reduction, we used Principal Component Analysis, followed by t-distributed stochastic neighbor embedding (TSNE)43. Additional data analysis was performed in custom Python 3.7 scripts using NumPy44, and pandas as statsmodels libraries.

To analyse the clustering of samples based on the different factors (Autosampler Drug, Braille Valves Drug, Combination) we used silhouette score analysis. Silhouette coefficient (b − a)/max(a,b) was calculated for each sample, where a was the mean intra- and b was the mean nearest-cluster distance. For each clustering factor, the average of Silhouette Coefficients were calculated (scikit-learn Python library). As silhouette score is dependent on the number of clusters, we created random clusters by permuting sample labels, thus cluster membership.

Pathway activities were calculated using PROGENy method17,18,19. Individual drug-specific pathway activities were calculated by fitting a linear model (pathway_activity ~ YM155 + Imatinib + Trametinib).

To compare the similarity between gene expression signatures from the high-throughput screen and LINCS-L1000 dataset8, we calculated consensus signatures for each drug of the high-throughput screen. To calculate consensus signatures, we fitted a linear model (gene_expression ~ Drug1 + Drug2 + … + Drugn) for each gene of the expression matrix and used the linear model coefficients as drug-specific signatures. To compare signature similarities, we calculated Spearman's rank correlation coefficient between the drug signatures of the high-throughput screen and LINCS-L1000 signatures. The similarity values were used as predicted values for ROC analysis, while true positives were the matched drug pairs between the high-throughput screen and LINCS-L1000.

For cell viability predictions, we used the CEVIChE method. CEVIChE predicts cell viability from gene expression changes based on a linear model, trained on a large compedion of matched cell viability and gene expression dataset18. As the measured genes of the high-throughput screen showed low overlap with the genes used by the original CEVIChE model, we retrained CEVIChE using only the genes measured in the high-throughput screen. This retrained CEVIChE model showed comparable performance (Pearson correlation between predicted and observed cell viability: 0.31) to the original method.

Drug-plates were prepared in advance in a 4 × 4 checkerboard for each combination, such that after the addition of cells, each drug was present at its GR35 concentration and a fourfold dilution series thereof. Each plate also contained media and DMSO negative controls and monotherapies for each drug. K562 cells were passaged the day before each experiment. On the day of the experiment, cells were washed once with PBS, then resuspended in FreeStyle 293 media (Thermo Fisher), containing 1% FBS. Cells were added using the multistep function of a multichannel pipette to each pre-prepared drug plate, such that each well had a 200 µL final volume and ~2 × 104 cells. The reservoir from which to aspirate cells was frequently refilled with a freshly resuspended stock solution to ensure that cells remained in suspension. Plates were sealed with a gas-permeable foil (Sigma) and incubated for 48 h. To prevent evaporation, plates were kept in the incubator within a box with ~1 cm water, within the box, the plates rested on tip boxes. After incubation, 22 µL PrestoBlue (Thermo Fisher) cell viability reagent was added to each well, and plates were resealed and returned to the incubator for 1 h. Plates were then read using a Tecan microplate reader with excitation/emission wavelengths of 535/615 nm (20 and 10 nm wavelength bandwidth, respectively). Based on the measured cell viability for monotherapies, we calculated expected cell viabilities using the Bliss independence model24 for each combination, for each concentration pair. The difference between expected and measured cell viability for combinations was averaged across all concentrations and was given as a synergy score.

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

The RNA-Sequencing data generated in this study have been deposited in the GEO database under accession code GSE174696. The fluorescence spectrometric data acquired on the microfluidic stations and viability measurements of cells in plates generated in this study are provided in the Source Data file. The logD values of drugs used in this study are available in the ChEMBL database45. Source data are provided with this paper.

Microfluidic control software can be downloaded from www.epfl.ch/labs/lbmm/downloads. The CAD files for the microfluidics chips are available from Zenodo: https://zenodo.org/record/6845607#.YtWZlMFBz1K. The code used to analyze the transcriptomic data can be downloaded from https://github.com/saezlab/Combi-Seq-analysis46.

Marquart, J., Chen, E. Y. & Prasad, V. Estimation of the percentage of US patients with cancer who benefit from genome-driven oncology. JAMA Oncol. 4, 1093–1098 (2018).

Article PubMed PubMed Central Google Scholar

Letai, A. Functional precision cancer medicine-moving beyond pure genomics. Nat. Med. 23, 1028–1035 (2017).

Article CAS PubMed Google Scholar

Cetin, A. E. et al. Determining therapeutic susceptibility in multiple myeloma by single-cell mass accumulation. Nat. Commun. 8, 1613 (2017).

Article ADS PubMed PubMed Central CAS Google Scholar

Gao, H. et al. High-throughput screening using patient-derived tumor xenografts to predict clinical trial drug response. Nat. Med. 21, 1318–1325 (2015).

Article CAS PubMed Google Scholar

Montero, J. et al. Drug-induced death signaling strategy rapidly predicts cancer response to chemotherapy. Cell 160, 977–989 (2015).

Article CAS PubMed PubMed Central Google Scholar

Menden, M. P. et al. Community assessment to advance computational prediction of cancer drug combinations in a pharmacogenomic screen. Nat. Commun. 10, 2674 (2019).

Article ADS PubMed PubMed Central CAS Google Scholar

Iorio, F. et al. A landscape of pharmacogenomic interactions in cancer. Cell 166, 740–754 (2016).

Article CAS PubMed PubMed Central Google Scholar

Subramanian, A. et al. A next generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell 171, 1437–1452 e1417 (2017).

Article CAS PubMed PubMed Central Google Scholar

Bush, E. C. et al. PLATE-Seq for genome-wide regulatory network analysis of high-throughput screens. Nat. Commun. 8, 105 (2017).

Article ADS PubMed PubMed Central CAS Google Scholar

Ye, C. et al. DRUG-seq for miniaturized high-throughput transcriptome profiling in drug discovery. Nat. Commun. 9, 4307 (2018).

Article ADS PubMed PubMed Central CAS Google Scholar

Srivatsan, S. R. et al. Massively multiplex chemical transcriptomics at single-cell resolution. Science 367, 45–51 (2020).

Article ADS CAS PubMed Google Scholar

McFarland, J. M. et al. Multiplexed single-cell transcriptional response profiling to define cancer vulnerabilities and therapeutic mechanism of action. Nat. Commun. 11, 4296 (2020).

Article ADS CAS PubMed PubMed Central Google Scholar

Clausell-Tormos, J. et al. Droplet-based microfluidic platforms for the encapsulation and screening of Mammalian cells and multicellular organisms. Chem. Biol. 15, 427–437 (2008).

Article CAS PubMed Google Scholar

Eduati, F. et al. A microfluidics platform for combinatorial drug screening on cancer biopsies. Nat. Commun. 9, 2434 (2018).

Article ADS PubMed PubMed Central CAS Google Scholar

Abate, A. R., Hung, T., Mary, P., Agresti, J. J. & Weitz, D. A. High-throughput injection with microfluidics using picoinjectors. Proc. Natl Acad. Sci. USA 107, 19163–19166 (2010).

Article ADS CAS PubMed PubMed Central Google Scholar

Miller, O. J. et al. High-resolution dose-response screening using droplet-based microfluidics. Proc. Natl Acad. Sci. USA 109, 378–383 (2012).

Article ADS CAS PubMed Google Scholar

Schubert, M. et al. Perturbation-response genes reveal signaling footprints in cancer gene expression. Nat. Commun. 9, 20 (2018).

Article ADS PubMed PubMed Central CAS Google Scholar

Szalai, B. et al. Signatures of cell death and proliferation in perturbation transcriptomics data-from confounding factor to effective prediction. Nucleic Acids Res. 47, 10010–10026 (2019).

Article CAS PubMed PubMed Central Google Scholar

Holland, C. H. et al. Robustness and applicability of transcription factor and pathway analysis tools on single-cell RNA-seq data. Genome Biol. 21, 36 (2020).

Article CAS PubMed PubMed Central Google Scholar

Janjic, A. et al. Prime-seq, efficient and powerful bulk RNA sequencing. Genome Biol. 23, 88 (2022).

Article CAS PubMed PubMed Central Google Scholar

Svensson, V. et al. Power analysis of single-cell RNA-sequencing experiments. Nat. Methods 14, 381–387 (2017).

Article CAS PubMed PubMed Central Google Scholar

Gruner, P. et al. Controlling molecular transport in minimal emulsions. Nat. Commun. 7, 10392 (2016).

Article ADS CAS PubMed PubMed Central Google Scholar

Hafner, M., Niepel, M., Chung, M. & Sorger, P. K. Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat. Methods 13, 521–527 (2016).

Article CAS PubMed PubMed Central Google Scholar

Holbeck, S. L. et al. The National Cancer Institute ALMANAC: a comprehensive screening resource for the detection of anticancer drug pairs with enhanced therapeutic activity. Cancer Res. 77, 3564–3576 (2017).

Article CAS PubMed PubMed Central Google Scholar

Paez, J. G. et al. EGFR mutations in lung cancer: correlation with clinical response to gefitinib therapy. Science 304, 1497–1500 (2004).

Article ADS CAS PubMed Google Scholar

Lynch, T. J. et al. Activating mutations in the epidermal growth factor receptor underlying responsiveness of non-small-cell lung cancer to gefitinib. N. Engl. J. Med. 350, 2129–2139 (2004).

Article CAS PubMed Google Scholar

Iyer, G. et al. Genome sequencing identifies a basis for everolimus sensitivity. Science 338, 221 (2012).

Article ADS CAS PubMed PubMed Central Google Scholar

Dietrich, S. et al. Drug-perturbation-based stratification of blood cancer. J. Clin. Invest 128, 427–445 (2018).

Article PubMed Google Scholar

Garnett, M. J. et al. Systematic identification of genomic markers of drug sensitivity in cancer cells. Nature 483, 570–575 (2012).

Article ADS CAS PubMed PubMed Central Google Scholar

Van Allen, E. M. et al. The genetic landscape of clinical resistance to RAF inhibition in metastatic melanoma. Cancer Disco. 4, 94–109 (2014).

Article CAS Google Scholar

Dagogo-Jack, I. & Shaw, A. T. Tumour heterogeneity and resistance to cancer therapies. Nat. Rev. Clin. Oncol. 15, 81–94 (2018).

Article CAS PubMed Google Scholar

Al-Lazikani, B., Banerji, U. & Workman, P. Combinatorial drug therapy for cancer in the post-genomic era. Nat. Biotechnol. 30, 679–692 (2012).

Article CAS PubMed Google Scholar

Sart, S., Tomasi, R. F., Amselem, G. & Baroud, C. N. Multiscale cytometry and regulation of 3D cell cultures on a chip. Nat. Commun. 8, 469 (2017).

Article ADS PubMed PubMed Central CAS Google Scholar

Headen, D. M., García, J. R. & García, A. J. Parallel droplet microfluidics for high throughput cell encapsulation and synthetic microgel generation. Microsyst. Nanoeng. 4, 17076 (2018).

Article CAS Google Scholar

Abe, H. et al. Discovery of a highly potent and selective MEK inhibitor: GSK1120212 (JTP-74057 DMSO Solvate). ACS Med. Chem. Lett. 2, 320–324 (2011).

Article CAS PubMed PubMed Central Google Scholar

Wheatley, S. P. & Altieri, D. C. Survivin at a glance. J. Cell Sci. https://doi.org/10.1242/jcs.223826 (2019).

Zhang, Y. et al. Sp1 and c-Myc modulate drug resistance of leukemia stem cells by regulating survivin expression through the ERK-MSK MAPK signaling pathway. Mol. Cancer 14, 56 (2015).

Article PubMed PubMed Central CAS Google Scholar

Wang, Z., Sampath, J., Fukuda, S. & Pelus, L. M. Disruption of the inhibitor of apoptosis protein survivin sensitizes Bcr-abl-positive cells to STI571-induced apoptosis. Cancer Res. 65, 8224–8232 (2005).

Article CAS PubMed Google Scholar

Stella, S. et al. Suppression of survivin induced by a BCR-ABL/JAK2/STAT3 pathway sensitizes imatinib-resistant CML cells to different cytotoxic drugs. Mol. Cancer Ther. 12, 1085–1098 (2013).

Article CAS PubMed Google Scholar

Carter, B. Z. et al. Regulation of survivin expression through Bcr-Abl/MAPK cascade: targeting survivin overcomes imatinib resistance and increases imatinib sensitivity in imatinib-responsive CML cells. Blood 107, 1555–1563 (2006).

Article CAS PubMed PubMed Central Google Scholar

Hennig, B. P. et al. Large-scale low-cost NGS library preparation using a robust Tn5 purification and tagmentation protocol. G3 8, 79–89 (2018).

Article CAS PubMed Google Scholar

Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

Article PubMed PubMed Central Google Scholar

Kobak, D. & Berens, P. The art of using t-SNE for single-cell transcriptomics. Nat. Commun. 10, 5416 (2019).

Article ADS PubMed PubMed Central CAS Google Scholar

Harris, C. R. et al. Array programming with NumPy. Nature 585, 357–362 (2020).

Article ADS CAS PubMed PubMed Central Google Scholar

Davies, M. et al. ChEMBL web services: streamlining access to drug discovery data and utilities. Nucleic Acids Res. 43, W612–W620 (2015).

Article CAS PubMed PubMed Central Google Scholar

Szalai, B. Combi-Seq for multiplexed transcriptome-based profiling of drug combinations using deterministic barcoding in single-cell droplets. GitHub https://doi.org/10.5281/zenodo.6761914 (2022)

Download references

We thank all members of the Genomics Core Facility at EMBL for their help in preparing libraries and sequencing our samples, the EMBL electronic and mechanical workshop for building controllers and holders for the Braille valves, the EMBL SIM for maintaining the mask aligner, Bart Deplancke and Cathrin Brisken for their feedback on the manuscript, and all current and previous members of the Merten lab for discussions and feedback. B.S. thanks the Premium Fellowship Program of the Hungarian Academy of Sciences (460044).

B. Szalai

Present address: Turbine Simulated Cell Technologies Ltd, Budapest, Hungary

These authors contributed equally: L. Mathur, B. Szalai.

Genome Biology Unit, European Molecular Biology Laboratory (EMBL), Heidelberg, Germany

L. Mathur, R. Utharala, M. Ballinger, J. J. M. Landry, V. Benes & C. A. Merten

Collaboration for joint PhD degree between EMBL and Heidelberg University, Faculty of Biosciences, Heidelberg, Germany

L. Mathur

Department of Physiology, Faculty of Medicine, Semmelweis University, Budapest, Hungary

B. Szalai

Institute of Enzymology, Research Centre for Natural Sciences, Budapest, Hungary

B. Szalai

Institute of Bioengineering, School of Engineering, École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland

N. H. Du & C. A. Merten

Université de Strasbourg, CNRS, Architecture et Réactivité de l’ARN, UPR, 9002, Strasbourg, France

M. Ryckelynck

Faculty of Medicine and Heidelberg University Hospital, Institute of Computational Biomedicine, Heidelberg University, Heidelberg, Germany

J. Saez-Rodriguez

Faculty of Medicine, Joint Research Centre for Computational Biomedicine (JRC‐COMBINE), RWTH Aachen University, Aachen, Germany

J. Saez-Rodriguez

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

L.M. conceptualized and developed the microfluidic platform and barcoding approach, designed and performed experiments, and analyzed the data under the supervision of C.A.M. B.S. devised and performed the computational analysis of the data under the supervision of J.S.-R. N.H.D. designed and performed experiments and analyzed the data. R.U. wrote the LabVIEW software used for the microfluidics platform. M.B. performed plate-based validation experiments. J.J.M.L. processed the sequencing data under the supervision of V.B. M.R. supported the development of the picoinjection procedures. L.M. and B.S. interpreted the results and wrote the manuscript with help from N.H.D., C.A.M., and J.S.-R. reviewed and edited the manuscript.

Correspondence to J. Saez-Rodriguez or C. A. Merten.

L.M., R.U., and C.A.M. are inventors of patent applications covering the barcoding strategy and microfluidic methods described here (WO2016207441A1). J.S.-R. has received funding from GSK and Sanofi and consultant fees from Travere Therapeutics and Astex Pharmaceuticals. B.S. is a full time employee of Turbine Ltd., Budapest, Hungary. The remaining authors declare no competing interests.

Nature Communications thanks Jeremy Jenkins, Angela Wu and the other, anonymous, reviewer(s) for their contribution to the peer review of this work

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Mathur, L., Szalai, B., Du, N.H. et al. Combi-seq for multiplexed transcriptome-based profiling of drug combinations using deterministic barcoding in single-cell droplets. Nat Commun 13, 4450 (2022). https://doi.org/10.1038/s41467-022-32197-0

Download citation

Received: 26 May 2021

Accepted: 21 July 2022

Published: 01 August 2022

DOI: https://doi.org/10.1038/s41467-022-32197-0

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

Nature Protocols (2022)

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.