A single-cell analysis of breast cancer cell lines to study tumour heterogeneity and drug response

Abstract

Cancer cells within a tumour have heterogeneous phenotypes and exhibit dynamic plasticity. How to evaluate such heterogeneity and its impact on outcome and drug response is still unclear. Here, we transcriptionally profile 35,276 individual cells from 32 breast cancer cell lines to yield a single cell atlas. We find high degree of heterogeneity in the expression of biomarkers. We then train a deconvolution algorithm on the atlas to determine cell line composition from bulk gene expression profiles of tumour biopsies, thus enabling cell line-based patient stratification. Finally, we link results from large-scale in vitro drug screening in cell lines to the single cell data to computationally predict drug responses starting from single-cell profiles. We find that transcriptional heterogeneity enables cells with differential drug sensitivity to co-exist in the same population. Our work provides a framework to determine tumour heterogeneity in terms of cell line composition and drug response.

Fig. 1. The breast cancer single-cell atlas.

A Representation of single-cell expression profiles of 35,276 cells from 32 cell lines color-coded according to cancer subtype (LA luminal A, LB luminal B, H Her2-enriched, TNA triple-negative type A, TNB triple-negative type B). B Expression levels of the indicated genes in the atlas, with red indicating expression, together with their C distribution within the cell lines, shown as a violin plot. D Dotplot of literature-based biomarker genes along the columns, for each of the 32 sequenced cell lines along the rows. Biomarker genes are grouped by type (Basal Epith. basal epithelial, Luminal Epith. luminal epithelial, L.P. luminal progenitor, EMT Epithelial to Mesenchymal Transition). E Graphical representation of 35,276 cells color-coded according to their cluster of origin. Clusters are numbered from 1 to 22. F For the indicated cluster, the corresponding pie chart represents the cluster composition in terms of cell lines. Cell lines in the same pie chart are distinguished by colour. Only the top 10 most heterogenous clusters are shown. In grey cell lines in the cluster contributing less than 5%, while the other colours represent a specific cell line. For example, Cluster 2 is the most heterogeneous cluster mainly composed of 3 cell line while cluster 19 is the most homogeneous since in its mainly composed by the cells coming from one cell line. G Expression levels in the atlas of the five luminal biomarkers identified as the most differentially expressed in each of the five luminal clusters (1, 2, 6, 8 and 18). H Expression of 22 atlas-derived biomarkers in the biopsies of 937 breast cancer patient from TCGA. I Accuracy in classifying tumour subtype for 937 patients from TCGA by using either the PAM50 gene signature or the 22 cluster-derived biomarker genes (scCCL) alone or augmented with HER2 gene (scCCL + HER2). Two-sided t-test is used to compare the performance of the different signatures. Source data are provided in a Source data file.

Fig. 2. Automatic classification of patients’ tumour cells.

A Cancer cells from triple-negative breast cancer (TNBC) biopsies of five patients were embedded in the BC atlas by means of the mapping algorithm in order to predict their tumour subtype. B For each patient, the pie chart shows cell-line composition obtained by mapping patient’s cells onto the atlas. C Tissue-slide of an oestrogen-positive breast tumour biopsy sequenced by means of the 10× Genomics Visium spatial transcriptomics (top-left) and the position of the mapped tissue tiles onto the atlas (top-left). Tiles are colour-coded according to the cell line (bottom-left) and to tumour subtype (bottom-right) as predicted by the mapping algorithm. D Cell-line composition of 937 BC patients from the TGCA database as estimated by the Bisque algorithm from their bulk RNA-seq data. For ease of interpretation, in the heatmap patients are clustered according to their cell-line composition. The bottom row reports the annotated cancer subtype in TGCA. E Predicted cell-line composition by the Bisque algorithm for four representative patients. F The distribution of the 937 BC patients across the 32 cell lines. For each cell line, the stacked bars report the percentage of patients of a given cancer subtype assigned by Bisque to that cell line. Since each patient is usually predicted by Bisque to be composed by multiple cell lines, the patient is associated to the cell line making up the largest fraction of the patient’s cell-line composition. G Performance of Bisque in classifying the tumour subtype of the 937 BC patients in TGCA from bulk gene expression profiles. Since each patient is usually predicted by Bisque to be composed by multiple cell lines, the patient is associated to the tumour type of the cell line making up the largest fraction of the cell-line composition. (PPV positive predictive value, AUC area under the curve). Source data are provided in a Source data file.

Fig. 3. Transcriptional heterogeneity in breast cancer cell lines and its impact on drug response.

A Percentage of cells expressing the indicated genes in each of the 32 cell lines. B Fluorescence cytometry of HCC38, MDA-MB-361 and AU565 cell lines stained with a fluorescent antibody against HER2. C Expression of HER2 protein in MDA-MB-361 cells is dynamic and re-established in less than 3 weeks. D Cell cycle phase for the HER2⁺ and HER2⁻ subpopulations of MDA-MB-361 cells. p-value refers to the Fisher’s exact test. E Enriched pathways (GSEA, FDR < 10%) across differentially expressed genes between the HER2⁺ (orange) and HER2⁻ (blue) MDA-MB-361 cells. F Gene expression versus drug potency for four anti-HER2 drugs. Each dot corresponds to a cell line with percentage of cells expressing ERBB2 or EGFR [y-axis] versus the experimental drug potency as Area Under the Curve (AUC) [x-axis]. PCC (Pearson correlation coefficient) and its p-value are also shown. G PCC values computed as in F for 66 drugs for which the cognate drug targets is known. The PCC distribution when choosing a random gene is also shown. Boxplots containing PCC distribution between a random gene and drug n = 1000, while n = 66 for boxplot containing PCC distribution between a drug and its cognate target gene. H Bioinformatics pipeline for the identification of drug sensitivity biomarkers for 450 drugs. I The top 250 most expressed genes in a single cell are used as input for a GSEA against the ranked list of genes correlated with drug potency for each one of the 450 drugs to predict its drug sensitivity. J Performance of DREEP in predicting drug sensitivity of 32 cell lines in the atlas to 450 drugs in terms of PPV (Positive Predicted Value) versus Recall. K Dose-response curve in terms of cell viability following treatment with either afatinib or etoposide at the indicated concentrations on sorted MDA-MB-361 cells (triplicate experiment). L Percentage of HER2⁺ cells in MDA-MB-361 after 72 h treatment with either afatinib or etoposide. (two-sided t-test), and M cell viability. N Percentage of HER2⁺ cells in MDA-MB-361 cell line at the indicated time-points either following 48 h of afatinib pre-treatment (red bars) or without any afatinib pre-treatment (black bars) and O the relative number of cells rescaled for the number of cells at the beginning of the experiment. Source data are provided in a Source data file. For K, L n = 3.