Skip to main page content
U.S. flag

An official website of the United States government

Dot gov

The .gov means it’s official.
Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

Https

The site is secure.
The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

Access keys NCBI Homepage MyNCBI Homepage Main Content Main Navigation
. 2020 Aug 27;11(1):4296.
doi: 10.1038/s41467-020-17440-w.

Multiplexed single-cell transcriptional response profiling to define cancer vulnerabilities and therapeutic mechanism of action

James M McFarland #  1 , Brenton R Paolella #  1 , Allison Warren  1 , Kathryn Geiger-Schuller  1   2 , Tsukasa Shibue  1 , Michael Rothberg  1 , Olena Kuksenko  1   2 , William N Colgan  1 , Andrew Jones  1 , Emily Chambers  1 , Danielle Dionne  1   2 , Samantha Bender  1 , Brian M Wolpin  3   4   5 , Mahmoud Ghandi  1 , Itay Tirosh  2   6 , Orit Rozenblatt-Rosen  1   2 , Jennifer A Roth  1 , Todd R Golub  1   3   7   8 , Aviv Regev  1   2   8   9   10 , Andrew J Aguirre  11   12   13   14 , Francisca Vazquez  15 , Aviad Tsherniak  16
Affiliations

Multiplexed single-cell transcriptional response profiling to define cancer vulnerabilities and therapeutic mechanism of action

James M McFarland et al. Nat Commun. .

Abstract

Assays to study cancer cell responses to pharmacologic or genetic perturbations are typically restricted to using simple phenotypic readouts such as proliferation rate. Information-rich assays, such as gene-expression profiling, have generally not permitted efficient profiling of a given perturbation across multiple cellular contexts. Here, we develop MIX-Seq, a method for multiplexed transcriptional profiling of post-perturbation responses across a mixture of samples with single-cell resolution, using SNP-based computational demultiplexing of single-cell RNA-sequencing data. We show that MIX-Seq can be used to profile responses to chemical or genetic perturbations across pools of 100 or more cancer cell lines. We combine it with Cell Hashing to further multiplex additional experimental conditions, such as post-treatment time points or drug doses. Analyzing the high-content readout of scRNA-seq reveals both shared and context-specific transcriptional response components that can identify drug mechanism of action and enable prediction of long-term cell viability from short-term transcriptional responses to treatment.

PubMed Disclaimer

Conflict of interest statement

A.R. is a co-founder and equity holder of Celsius Therapeutics, equity holder in Immunitas, and an SAB Member of ThermoFisher Scientific, Neogene Therapeutics, Asimov, and Syros Pharmaceuticals. A.J.A. has consulted for Oncorus, Inc., Arrakis Therapeutics, and Merck & Co., Inc. and has research funding from Mirati Therapeutics and Deerfield Management that are unrelated to this work. A.T. is a consultant for Tango Therapeutics. T.R.G. is a consultant to GlaxoSmithKline, a founder of Sherlock Biosciences, and was formerly a consultant and equity holder in Foundation Medicine, acquired by Roche. T.R.G. also receives research funding unrelated to this project from Bayer Healthcare. B.M.W. receives research funding unrelated to this project from Celgene and Lilly and is a consultant for BioLineRx, Celgene, G1 Therapeutics, and GRAIL. The Broad Institute filed a United States Provisional Patent Application directed to work described in this article. The pending Application is entitled “Rapid Prediction of Drug Responsiveness” and was filed on September 2, 2019. The current inventor list includes A.T., A.J.A., F.V., B.R.P., and J.M.M. A.W., K.G.-S., T.S., M.R., O.K., W.N.C., A.J., E.C., D.D., S.B., M.G., I.T., O.R.-R., and J.A.R. declare no competing interests.

Figures

Fig. 1
Fig. 1. Multiplexed transcriptional profiling across pools of cell lines.
a Schematic diagram illustrating the MIX-Seq platform. b Heatmap showing likelihoods assigned by the SNP classification model for each cell coming from each parental cell line. The model consistently picks out a single cell line, from among the 24 “in-pool” cell lines, with high confidence. c UMAP representation of cells treated with DMSO control (gray) or nutlin (red) across a pool of 24 cell lines. Arrows indicate the shift in the population median coordinates for each cell line. d Heatmap showing average log fold-change estimates for each cell line for top differentially expressed genes. Nutlin sensitivity is given by 1 − area under dose response curve (AUC, see “Methods”). e Volcano plot showing strong gene expression changes in response to nutlin treatment across TP53 WT cell lines (n = 7). Effect size estimates and p values (not corrected for multiple comparisons) for this and subsequent differential expression analyses are estimated using the limma-trend pipeline, (“Methods”). Vertical lines indicate a logFC threshold of 1. f Same as e for TP53 mutant cell lines (n = 17), showing little gene expression change in response to nutlin treatment. g Gene set analysis identifies gene sets that are upregulated (right) and downregulated (left) by nutlin treatment in the TP53 WT cell lines.
Fig. 2
Fig. 2. Viability-related and -independent transcriptional response components.
a UMAP representation of single-cell expression profiles in a 99 cell line pool treated with vehicle control (gray) or trametinib (red). Arrows indicate trametinib-induced shift of population median coordinates for each cell line. b Histogram showing the number of cells recovered in each cell line and condition. c Volcano plot showing the viability-independent response for each gene, representing the “y-intercept” of a linear fit of expression change to drug sensitivity (1-AUC; “Methods”). Inset at left shows this relationship for an example gene: EGR1. The blue line shows the linear regression trend line (with the 95% CI interval shown in shaded gray). Inset below shows top upregulated (right, red) and downregulated (left, blue) gene sets. d Same as c for the viability-related response.
Fig. 3
Fig. 3. Machine learning analysis powered by large-scale transcriptional profiling.
a Accuracy of models trained to predict sensitivity for each individual drug, using either measured transcriptional responses or baseline “omics” features of the cell lines (using the same set of cell lines). Predictions based on transcriptional profiling at 6 and 24 h posttreatment are indicated by the gold and blue dots, respectively. b Plot of the Pearson correlation between each gene’s transcriptional response (24 h posttreatment) and drug sensitivity (1 − AUC) across cell lines and drugs, compared with its feature importance for the random forest predictive models. Gene set enrichment analysis of highlighted genes shown above. c Matrix of measured transcriptional responses across cell lines 24 h after trametinib treatment. d Eigenvalue spectrum of PCA applied to the matrix from c. e The projection of each cell line’s response onto PC1 is plotted against its measured trametinib sensitivity (red and green points indicate cell lines with activating BRAF or KRAS mutations respectively). Linear regression trend line (with the shaded gray region showing the 95% CI interval) is shown in blue. f Comparison of PC2 scores with expression of the melanoma-specific transcription factor SOX10 across cell lines. g UMAP embedding of transcriptional response profiles across drugs, cell lines, and posttreatment time points. Points are colored by treatment condition (drug and time point). Inset below shows zoomed view of the region indicated by the rectangle, with the larger green dots representing responses of BRAF mutant melanoma lines to the BRAFi dabrafenib.
Fig. 4
Fig. 4. Population heterogeneity in pre- and post-perturbation transcriptional programs.
a UMAP representation of gene expression profiles of DMSO (smaller points) and nutlin-treated (larger points) cell populations for a pool of 24 cell lines (as in Fig. 1). Cells are colored by their inferred cell cycle phase. TP53 WT cell lines (red arrows) show predominance of G0/G1-phase cells after nutlin treatment not observed in TP53 mutant cell lines (black arrows). b Quantification of the change in proportion of cells in G0/G1 in each cell line (y-axis) shows that nutlin treatment elicits G1 arrest selectively among the TP53 WT cell lines. Error bars show the 95% CI for estimates of the change in G0/G1 cell proportions. c Average change in the fraction of cells in each cell cycle phase for each drug treatment (averages are weighted by measured drug sensitivity, all for 24 h time points). d UMAP plot showing two subpopulations of RCC10RGB cells emerging 24 h after treatment with bortezomib. Dot color depicts inferred cell cycle phase and dot size represents treatment condition. e Heatmap showing LFC relative to control cells of top differentially expressed genes between the two clusters of bortezomib-treated cells for each of the 10/24 cell lines that showed a bimodal response pattern. In nearly all cases, one of the clusters (labeled cluster 2 by convention) was characterized by a predominance of S-phase cells, compared with mostly G0/G1 cells in cluster 1. f Volcano plot showing differentially expressed genes in cluster 2 vs cluster 1 averaged across the nine cell lines depicted in e. Genes that are part of the S-phase signature are highlighted in red.
Fig. 5
Fig. 5. Dual-multiplexed transcriptional profiling across cell lines and time points.
a Schematic diagram illustrating experiment using Cell Hashing to multiplex scRNA-seq of cell line pools sampled at different time points following drug treatment. b UMAP plot showing 13,713 cells across a pool of 24 cell lines at different times following treatment with trametinib (shades of blue) or DMSO control (pink). c Single-cell expression levels of EGR1 at different time points following trametinib treatment for an example insensitive/sensitive cell line (left/right). Red dots depict the mean expression levels at each time point, and error bars show the interval +/− s.e.m. d Same as c, for MCM7. e (Top) Time course of the viability-independent response for top downregulated genes. (Bottom) Enrichment of HALLMARK_KRAS_SIGNALING_UP genes in the downregulated viability-independent response at each time point. f (Top) Same as e, showing time course of the viability-related response for top downregulated genes. (Bottom) Enrichment of HALLMARK_G2M_CHECKPOINT genes in the viability-related response at each time posttreatment. g Average time course of G0/G1 arrest across cell lines (n = 24 cell lines). Error bars indicate interval +/− s.e.m.
See this image and copyright information in PMC

Similar articles

See all similar articles

Cited by

See all "Cited by" articles

References

    1. Tsherniak A, et al. Defining a cancer dependency map. Cell. 2017;170:564–576.e16. - PMC - PubMed
    1. Meyers RM, et al. Computational correction of copy number effect improves specificity of CRISPR-Cas9 essentiality screens in cancer cells. Nat. Genet. 2017;49:1779–1784. - PMC - PubMed
    1. McDonald ER, et al. Project DRIVE: a compendium of cancer dependencies and synthetic lethal relationships uncovered by large-scale, deep RNAi screening. Cell. 2017;170:577.e10–592.e10. - PubMed
    1. Iorio F, et al. A landscape of pharmacogenomic interactions in cancer. Cell. 2016;166:740–754. - PMC - PubMed
    1. Behan FM, et al. Prioritization of cancer therapeutic targets using CRISPR-Cas9 screens. Nature. 2019;568:511–516. - PubMed

Publication types