Single-Cell Map of Diverse Immune Phenotypes in the Breast Tumor Microenvironment

Abstract

Knowledge of immune cell phenotypes in the tumor microenvironment is essential for understanding mechanisms of cancer progression and immunotherapy response. We profiled 45,000 immune cells from eight breast carcinomas, as well as matched normal breast tissue, blood, and lymph nodes, using single-cell RNA-seq. We developed a preprocessing pipeline, SEQC, and a Bayesian clustering and normalization method, Biscuit, to address computational challenges inherent to single-cell data. Despite significant similarity between normal and tumor tissue-resident immune cells, we observed continuous phenotypic expansions specific to the tumor microenvironment. Analysis of paired single-cell RNA and T cell receptor (TCR) sequencing data from 27,000 additional T cells revealed the combinatorial impact of TCR utilization on phenotypic diversity. Our results support a model of continuous activation in T cells and do not comport with the macrophage polarization model in cancer. Our results have important implications for characterizing tumor-infiltrating immune cells.

Keywords: Bayesian modeling; T cell activation; TCR utilization; breast cancer; single-cell RNA-seq; tumor microenvironment; tumor-infiltrating immune cells.

Figure 1:

Single-Cell RNA-Seq Experimental Design and Initial Data Exploration (A) Flow chart of experimental design and analysis. (B) Summary of samples and patient metadata; more details in S1A. (C) t-SNE of complete immune systems from two example breast cancer tumors. scRNA-seq data for each tumor is processed with SEQC (Figure S1B) and library size-normalized; each dot represents a cell colored by PhenoGraph clusters, labeled by inferred cell types. Additional tumors are shown in Figure S1C. (D) Pie charts of cell type fractions for each patient’s tumor-infiltrating immune cells, colored by cell type. (E) Left: Boxplots of expression of Hallmark Hypoxia signature (defined as the mean normalized expression of genes in the signature) across immune cells from each patient. Right: Heatmap of z-scored mean expression of genes in signature. Top: Barplot of total expression of each gene, across all patients. See Figure S1E–G for additional signatures.

Figure 2:

Unbiased Characterization of the Immune System Across Breast Cancer Patients (A) t-SNE of immune cells from 4 breast cancer tumors after library-size normalization (left) and Biscuit normalization (right). Cells are colored by tumor. Less mixing of tumors indicates either batch effects or patient-specific cell states. (B) Left: Boxplots of expression of CD8 T cell activation signature (defined as the normalized mean expression of genes in the activation signature in Table S4) across immune cells from each patient, in the same format as Figure 1E. Expression of T cell activation signature shows variability across patients. (C) t-SNE map of immune cells from all 8 breast tumors after Biscuit normalization and imputation showing rich structure and diverse cell types. Cells colored by Biscuit clusters and labeled with inferred cell types. (D) Histogram depicting entropy of the tumor (patient) distribution as a measure of sample mixing. Entropy is computed per cell, based on the distribution of patients in (30-NN) local cell neighborhoods after library-size normalization (left) as compared to Biscuit (right). (E) t-SNE of complete atlas of immune cells, post-Biscuit normalization, from all patients and tumor, blood, lymph, and contralateral normal tissues, labeled by inferred cell type (left) and normalized expression of 8 immune cell markers (right). Figure S2, S3A, Table S2 present further details on clusters. (F) Pearson correlations between cluster expression centroids and bulk RNA-seq data from purified immune populations from (Jeffrey et al., 2006; Novershtern et al., 2011). (G) Histogram of frequency of patients contributing to each cluster; 19 clusters are present in all 8 patients and 10 clusters are patient-specific. (H) Z-score normalized expression of canonical and cell type markers across clusters. (I,J) Differentially expressed genes (DEGs) in B cell (I) and NK cell (J) clusters, standardized by z-scores within cell type to highlight clusters with higher or lower expression of the marker compared to the average B or NK cluster (all DEGs are presented in Table S3).

Figure 3:

Impact of the Microenvironment on Breast Immune Cells (A) Breast immune cell atlas constructed from combining all patient samples (BC1–8) and tissues using Biscuit, projected with t-SNE. Each dot represents a cell, colored by cluster; major cell types are marked according to Figure 2F, H and Table S2, 3. (B) Subsets of immune atlas t-SNE in (A) showing cells from each tissue separately on the same coordinates as 3A to highlight the differences between tissue compartments. (C) Proportions of cell types across tissue types. (D) Distribution of variance of expression, computed for each gene across all immune cells (all patients), from tumor tissue compared to normal breast tissue. (E) Most significant Hallmark GSEA enrichment results on genes with highest difference in variance in tumor T cells vs normal tissue T cells. See Figure S3B, C for enrichment in monocytic and NK cells. Full lists of enrichments are presented in Table S5. (F) Phenotypic volume in log-scale (defined as pseudo-determinant of gene expression covariance matrix, detailed in STAR) of immune cell types in tumor compared to normal tissue, controlled for number of cells, showing siginificant (shown with asterisks) expansion of volume spanned by independent phenotypes in tumor compared to normal tissue.

Figure 4:

Detailed Characterization of T Cells (A) Visualization of all T cells using first three diffusion components (two uninformative components denoting isolated clusters were removed). Each dot represents a cell colored by cluster, and by tissue type in insert. The main trajectories are indicated with arrows and annotated with the signature most correlated with each component. See Figure S4D for additional components. (B) Traceplot of CD8 T cell activation signature (defined as mean expression across genes in signature in Table S4) for all T cells along activation component. Cells are projected along the component (x-axis), and the blue line indicates the moving average of signature expression, using a sliding window of length equal to 5% of total number of T cells; shaded area displays standard error. (C) Heatmap showing expression of immune-related genes with the largest positive correlations with activation component, averaged per cluster and z-score standardized across clusters; columns (clusters) are ordered by mean projection along the component. See Figure S4 for additional components. (D) Violin plot showing the density of all T-cells (left), T cells in individual tissues (middle), and in individual clusters (right), along activation component. Number of dots inside each violin are proportional to number of cells. (E) Trace-plots (as in B) of (left) exhaustion/terminal differentiation signature along second component and (right) hypoxia signature along third component. Signatures are presented in Table S4. (F) Heatmap of cells projected on each diffusion component (rows) averaged by cluster (columns).

Figure 5.

Covariance Patterns Help Define Distinct T Cell Clusters (A, B, C) Heatmaps of mean expression for a curated set of transcriptomic signatures (Table S4) for (A) CD4 memory T cell, (B) CD8 memory T cell, and (C) Treg clusters. Only signatures with high expression in at least one T cell cluster are shown. Expression values are z-scored relative to all T cell clusters. (D) Cartoon illustration of two clusters of cells with similar mean expression for two example genes but opposite covariance between the same two genes. Distinct patterns in both mean and covariance of expression define clusters in Biscuit. (E) Scatter plot showing mean expression of GITR vs. CTLA-4 for each T cell cluster (represented by a dot). Treg clusters (red) show high mean expression of both genes. (F) Distribution of covariance values between GITR and CTLA-4 across all T cell clusters, with Treg clusters marked in red. Treg cluster covariance values exhibit differences despite sharing high mean expression levels. See Figure S5A for similar computation on the raw, un-normalized, data, verifying the result. See Figure S5C for similar results in CyTOF data. (G) Network visualization illustrating strength of covariance between pairs of checkpoint receptor genes in Treg clusters. Edge width denotes absolute magnitude and red and blue colors denote positive and negative signs of covariance respectively; the case of CTLA-4 and GITR is highlighted in yellow. Figure S5D shows networks for other T cell populations. (H) Heatmaps of covariance values for immune genes in two Treg clusters showing different modules of covarying genes. (I) Proportion of Treg clusters in each patient, indicating that differences in covariance patterns between clusters translate to patients.

Figure 6:

TCR Repertoire Shapes Diverse Phenotypic States (A) Pearson correlation between centroids of differentially expressed genes in T cell clusters inferred from BC1–8 patients using inDrop (rows) and clusters inferred from 27,000 T cells from BC9–11 tumors using 10× (columns) (see Figure S6A and STAR for details), showing near one-one-to-one mapping of clusters. (B) Same as in (A) computing Bhattacharyya similarity between pairs of clusters, accounting for both mean and covariance of clusters. Further details on clusters are presented in S6B, C. (C) Histogram of activation states of (top) all T cells from three tumors BC9–11 and (bottom) T cells separated by each of the top 20 most dominant TCR clonotypes in BC9, mapped using paired single-cell RNA and TCR sequencing. Frequencies of clonotypes are shown in S6D. Similar figures for BC10,11 are shown in Figure S6E. (D) Heatmap showing normalized mean expression levels for a curated set of transcriptomic signatures (rows, listed in Table S4) for T cell clusters in BC9–11. (E) Distributions of each CD8+ T or Treg cluster in BC9–11 across the 30 most frequent clonotypes from each tumor. Clusters (columns) are z-scored to highlight the combinatorial impact of clonotypes in shaping each phenotypic state, and sorted by activation level. Some clusters associated with the same clonotype have the same level of activation (seen as connected horizontal stretches in the heatmap), while others have similar environmental responses (Table S6; STAR). (F) t-SNE of normalized single-cell RNA-seq data for T cells from BC9–11 tumors colored by markers, T cell activation, and tumor (left); Biscuit clusters (middle, top); and examples of dominant clonotypes from each tumor identified with paired TCR sequencing, projected in the same coordinates (right). Separate projection of each dominant clonotype for each tumors is shown in Figure S6F.

Figure 7:

Detailed Characterization of Myeloid Cells (A) t-SNE map projecting myeloid cells from BC1–8 patients (all tissues). Cells are colored by Biscuit cluster and cell types are labeled based on bulk RNA-seq correlation-based annotations. (B through E) Projection of myeloid cells on macrophage activation, pDC, and monocyte activation diffusion components, colored by (B) cluster, (C) tissue, (D) cell type, and (E) expression of example lineage-demarcating genes. (F) Violin plots showing the density of cells along macrophage activation component and organized by overall density (left panel), tissue type (middle panel), and cluster (right panel). See Figure S7 for other components. (G) Scatter plot of normalized mean expression of M1 and M2 signatures per cell (dot); cells assigned to TAM clusters have been highlighted by cluster. (H) Scatterplot of mean expression of MARCO and CD276 in myeloid clusters; each dot represents a cluster; TAM clusters are marked in red, indicating high expression of both markers in macrophage clusters. (I) Distribution of covariance between MARCO and CD276 across all myeloid clusters. TAM clusters are marked in red and present substantial outliers. See Figure S7F for similar computation on the raw, un-normalized data, verifying the result. (J) Heatmaps showing covariance patterns of select macrophage marker genes in 3 TAM clusters.