DIALOGUE maps multicellular programs in tissue from single-cell or spatial transcriptomics data

Abstract

Deciphering the functional interactions of cells in tissues remains a major challenge. Here we describe DIALOGUE, a method to systematically uncover multicellular programs (MCPs)-combinations of coordinated cellular programs in different cell types that form higher-order functional units at the tissue level-from either spatial data or single-cell data obtained without spatial information. Tested on spatial datasets from the mouse hypothalamus, cerebellum, visual cortex and neocortex, DIALOGUE identified MCPs associated with animal behavior and recovered spatial properties when tested on unseen data while outperforming other methods and metrics. In spatial data from human lung cancer, DIALOGUE identified MCPs marking immune activation and tissue remodeling. Applied to single-cell RNA sequencing data across individuals or regions, DIALOGUE uncovered MCPs marking Alzheimer's disease, ulcerative colitis and resistance to cancer immunotherapy. These programs were predictive of disease outcome and predisposition in independent cohorts and included risk genes from genome-wide association studies. DIALOGUE enables the analysis of multicellular regulation in health and disease.

Extended Data Figure 1.. DIALOGUE identified MCPs in the mouse hypothalamus that are not recovered with other dimensionality reduction and clustering approaches.

(a)* Pearson correlation coefficient between genes, PCs, NMF, and DIALOGUE MCPs from either the training or the test set (x axis) across different pairs of cell types (panels) in spatial niches in the mouse hypothalamus. (b) Pearson correlation coefficient (red/blue, color bar) between the Overall Expression of the relevant MCP component when considering only defined subsets of the pertaining cell types (rows, columns), as previously identified by clustering. White: missing values (cell subtypes that cannot be compared). (c) MCPs are not merely driven by cell subtype composition in a niche. Fraction of cells from different clusters (as previously defined, y axis) among cells of a given type (label on top) that over- or under-express the relevant component of each pair-wise MCP1 (top or bottom 25%, respectively, x axis) involving that cell type. (d)* Similarity (y axis, Spearman’s r) between the gene loadings of MCPs identified in the microenvironment setting (x axis) and the gene loadings of matching MCPs identified in the micro-environment setting, when computed for different pairs of cell types using MERFISH data. *In both (a) and (d) middle line: median; box edges: 25^th and 75^th percentiles, whiskers: most extreme points that do not exceed ±IQR*1.5; further outliers are marked individually.

Extended Data Figure 2.. DIALOGUE captures spatial patterns.

(a) Average Overall Expression in a niche (dot, 15 cells on average) of the first MCP (MCP1) in the first (x axis) and second (y axis) cell type in that MCP. In red is the locally weighted polynomial (LOWESS) regression line. (b) As in (a), but depicting the Overall Expression residuals after regressing out impact of cell clusters, as previously defined. (a-b) Spearman correlation coefficient (R) and significance (P, one-sided). (c) Performance (AUROC, y axis) when predicting the expression of the corresponding DIALOGUE component in the neighboring cells located in the same macro-environment (dark blue, ~500 cells) or micro-environment (purple, and light blue, ~15 cells), when testing on unseen test set; the training data includes either spatial coordinates and single-cell profiles (light blue, “spatial data”) or only single cell profiles from ~500 cell aggregates, without spatial information (“dissociated”, Online Methods).

Extended Data Figure 3.. MCPs mark spatial patterns and phenotypes.

(a) Overall Expression of HMRF domain-specific programs in neighboring pairs of glutamatergic (y axis) and GABAergic (x axis) neurons from different regions (colors). (b) Overall Expression of the relevant components of MCPs 1–5 in glutamatergic (y axis) and their adjacent GABAergic (x axis) neurons from different regions (colors). (a-b) Spearman correlation coefficient (R) and significance (P, one-sided).

Extended Data Figure 4.. MCPs mark spatial patterns and phenotypes.

(a) Spatial distribution of MCPs and HNRF programs. Overall Expression of MCPs identified by DIALOGUE and the HMRF domain programs in glutamatergic (circles) and GABAergic (dots) neurons in the mouse visual cortex. As shown, while many of the patterns follow either a more layered or salt and pepper pattern, MCP2 distinguished a more discrete region. While such boundaries sometimes reflect measurement artifacts, we did not find an association with number of genes/reads (typical quality measures) nor with simple alignment with Fields of View (FOV). (b,c) Shared and cell type specific components in DIALOGUE MCP1s in the mouse hypothalamus. (b) Fraction of genes (y axis) that are shared (yellow) or specific to one (A, dark blue) or another (B, light blue) of the cell types in each of the hypothalamus MCPs (x axis). (c) The two cell programs in each of the MCPs in (b) and their specific and shared (intersection) genes. P-values denote association with naïve animal behavior (multilevel mixed-effects models, two-sided test).

Extended Data Figure 5.. DIALOGUE identifies mis-localized cells and disease MCPs in single-cell data.

(a) ROC curves showing the true positive (y axis) and false positive (x axis) rate when predicting mis-localized cells of each major subset (panels) with different types of “contamination” with cells that are either from the same layer (LP/EPI) within control (black, from replicate biopsy) or UC (blue; from adjacent biopsy with a different clinical status: inflamed or non-inflamed); or from a different layer but same clinical status, when considering either all samples (green) or only samples from control (yellow) or UC patients (red). (b) UC multicellular program genes. Average expression (Z score residuals after regressing out the associations with the LP/EPI location, red/blue color bar) of top genes (columns) from the UC multicellular program, sorted by their pertaining cell type (top color bar), across samples (rows), sorted by Overall Expression (right, Online Methods), and labeled by clinical status, location and patient ID (left color bar). (c) Melanoma MCP1. Average expression (Z score, red/blue color bar) of top genes (columns) from MCP1 identified in four different cell types (top color bar), across melanoma tumor samples (rows), sorted by Overall Expression of MCP1 (right, Online Methods), and labeled by treatment status and ICB response (left color bar).

Figure 1.. DIALOGUE: a method for MCP identification.

(a) Multicellular programs (MCPs). Upper panel: MCPs can arise due to either shared cues (upper row) or direct cell-cell interactions (bottom row) and comprise both cell type specific (left column) or shared (right column) programs. Lower panel: MCPs can be observed also in dissociated tissues based on co-variation across samples. (b) DIALOGUE method. In Step I (left), DIALOGUE takes as input (left) mean gene/feature expression (columns) matrices for each cell type, across samples or physical niches (rows) and infers “sparse latent variates” across those samples/niches (middle matrices) and their activity (right matrices) for each cell type, such that the programs of each cell type are highly correlated with the corresponding programs of all the other cell types (in this case shown only for two). In Step II (right), it identifies a gene signature for each sparse canonical variate, by interrogating single-cell distributions, while accounting for confounding factors at different levels. This results in MCPs, each with a set of up- and down-regulated genes in each of its cell type compartments (right). Ligand-receptor interactions can then be used to identify potential mediators. (c-d) Spatial distribution of MCPs. Expression (color code) of MCP2 and MCP4 components in single cells (dots) identified for excitatory and inhibitory neurons (c) and astrocytes and inhibitory neurons (d); “high”: expression above the 80^th percentile; “moderate/low” otherwise. (e) Pearson correlation coefficient between genes, PCs, NMF, or DIALOGUE MCPs from either the training or the test set (x axis) across different pairs of cell types in spatial niches in the mouse hypothalamus (see Extended Data Fig. 1a for all other pairwise combinations). Middle line: median; box edges: 25^th and 75^th percentiles, whiskers: most extreme points that do not exceed ±IQR*1.5; further outliers are marked individually. (f) Receiver-Operator Curves (ROCs) for predictions of the expression of each MCP in unseen cells based on the expression profiles of their neighbors in the “micro-“ (radius of 37 μm, ~15 cells) or “macro-“ (radius of 220 μm, ~500 cells) environment, using spatial coordinates in the training data (black and purple) or only single-cell RNA profiles grouped to samples (red/blue, green/orange).

Figure 2.. DIALOGUE MCPs recover spatial information from different mouse brain regions and spatial genomics data types.

(a-b) DIALOGUE outperformed HMRF in extracting spatial signatures from Seq-FISH data from the mouse visual cortex. (a) Overall Expression of the HMRF O3 domain programs (left) and corresponding MCP1 components (right) in neighboring pairs of glutamatergic (y axis) and GABAergic (x axis) neurons from different regions (colors). Spearman correlation coefficient (R) and two-sided p-values (P) are shown. (b) HMRF-domain assignments (left) and Overall Expression of (middle) MCP1, and MCP2 (right) in glutamatergic (circles) and GABAergic (dots) neurons in their spatial context. (c-e) MCPs predict adjacent cells’s states on unseen test data. (c) Predictive accuracy (Area Under the Receiver Operator Curve, AUROC, y axis) for neighboring cells for different pairs of cell types (x axis), when using the Euclidean distance based on MCPs learned on training data (red), Principal Components from PCA (blue), or the original gene expression data (green). (d,e) True positive (y axis) and false positive (x axis) rates for predicting neighboring cells based on MCPs (red), PCs (blue) or gene expression (green) in mouse cerebellum Slide-Seq data (d) or scRNA-seq from spatially distinct regions in the mouse neocortex (e). AUROC values are shown in parentheses. (f) MCPs in the cerebellum distinguish coordinated expression between neurons and astrocytes in the inner and outer layers. Overall Expression (two left panels) in neurons and astrocytes and discretized expression (two right panels; (“high”: expression above the 80^th percentile; “moderate/low” otherwise)) of MCPs identified in the mouse cerebellum based on Slide-Seq data. (g) Overall Expression (y axis, Online Methods) of the excitatory-inhibitory MCP in male (left) or female (right) mice exhibiting different behaviors (legend). Middle line: median; box edges: 25^th and 75^th percentiles, whiskers: most extreme points that do not exceed ±IQR*1.5; further outliers are marked individually. (h) The two cell programs in each of two MCP1s and their specific and shared (intersection) genes. P-values denote association with naïve animal behavior (multilevel mixed-effects models, two-sided test).

Figure 3.. MCPs in NSCLC identify coordinated interferon responses in immune and stroma cells at the tumor edge.

Cell type assignments and MCPs in representative NSCLC tumors profiled by SMI data. (a) MCP2 in NSCLC tumor #13. Discretized expression (left, “High” > 80^th percentile; “moderate/low” otherwise), its inverse (middle, highlighting cells with low MCP expression), and Overall Expression (right) of MCP2 of cells (dots) in the tumor spatial context. (b) MCP1 in NSCLC tumor #13. Discretized expression (left) and Overall Expression (middle) of MCP1, and cell type annotations (right) of cells (dots) in the tumor spatial context. (c) MCP1 in in NSCLC tumor #9. Discretized (left) and Overall Expression (middle) of MCP1 and cell type annotations (right) of cells (dots) in the tumor spatial context.

Figure 4.. MCPs associated with ulcerative colitis and predictive of clinical outcomes.

(a) Deviation from multicellular patterns identifies mis-localized cells. ROC curves show the true positive (y axis) and false positive (x axis) rate when predicting mis-localized cells of each cell subset (color legend). (b) MCP1 in the human colon. Off diagonal panels: Comparison of Overall Expression scores (y and x axes) for each cell component of MCP1 (rows and columns, labels on diagonal) across the samples (dots, black: control, blue: non-inflamed IBD, red: inflamed IBD); the lines correspond to the linear fit. Pearson correlation (r) and significance (^***P < 1*10⁻³, one-sided) are shown in the panels above the diagonal. Diagonal panels: Distribution of Overall Expression scores for each cell type component, along with the kernel density estimates. (c) Most genes in five colon MCPs are specific to one cell component. Number of genes (y axis) in the up (left bars) and down (right bars) of each of five colon MCPs (x axis) that appear only in one cell type component or in multiple ones (color legend). (d) MCP1 is induced in UC samples. Distribution of Overall Expression scores (x axis) of MCP1 (“UC program”) across the different cell subtypes (y axis) for cells from UC patients (grey) and healthy individuals (light blue); black line denotes the mean of the distribution. P: One-sided p-values, mixed-effects models. (e) UC multicellular program genes. Average expression (Z score, red/blue color bar) of top genes (columns) from the UC multicellular program, sorted by their pertaining cell type (top color bar), across samples (rows), sorted by Overall Expression (right, Online Methods), and labeled by clinical status, location and patient ID (left color bar). (f,g) UC multicellular program predicts response to anti-TNF therapy. (f) Overall Expression (y axis) of the UC-program in bulk RNA-seq of colon biopsies from 24 UC patients pre- and post-infliximab infusion, stratified to responders (blue) and non-responders (grey) and in normal mucosa from 6 control patients (red). Middle line: median; box edges: 25^th and 75^th percentiles, whiskers: most extreme points that do not exceed ±IQR*1.5; further outliers are marked individually. P-value: linear regression model (Online Methods). (g) ROC curve obtained when using the UC multicellular program score in pre-treatment samples to predict the subsequent clinical responses to infliximab infusion in UC patients.

Figure 5.. MCP in the prefrontal cortex associated with AD pathology and aging.

(a) AD MCPs. Average expression (Z score, red/blue color bar) of top genes (columns) from the AD multicellular program, sorted by their pertaining cell type (top color bar), across samples (rows), sorted by Overall Expression (right, Online Methods), and labeled by clinical status (left color bar). (b) AD programs components across cell types. Off diagonal panels: Comparison of Overall Expression scores (y and x axes) for each cell component of MCP1 (rows and columns, labels on diagonal) across the samples (dots, black: control; red: AD); the lines correspond to the linear fit. Pearson correlation coefficient (r) and significance (^***P < 1*10⁻³, one-sided) are shown in the upper triangle. Diagonal panels: Distribution of Overall Expression scores for each cell type component, along with the kernel density estimates. In: inhibitory neurons; Ex: excitatory neurons, Oli: oligodendrocytes, Opc: oligodendrocyte-precursor cells, Ast: astrocytes, Mic: microglia. (c) Most genes in the AD MCPs are specific to one cell component. Number of genes (y axis) in the up-regulated (left bar) and down-regulated (right bar) compartment of the AD MCPs (x axis) that belong to one cell type component or multiple ones (color legend). (d) AD MCPs components are induced across cell types in AD. Distribution of Overall Expression scores (y axis) of AD multicellular program component across the different cell subtypes (y axis) for cells from AD patients (grey) and neurologically normal subjects (light blue); black line denotes the mean of the distribution. P: One-sided p-values, mixed-effects models. (e) The overall expression of AD MCPs in brain autopsies of AD (red) and non-AD (grey) individuals. (f) AD MCPs increase with age in the frontal cortex and cerebellum of neurologically normal subjects. Overall Expression (y axis) of the AD MCPs in bulk RNA-seq of the frontal cortex (left) and cerebellum (right) of neurologically normal subjects stratified by age (x axis). Middle line: median; box edges: 25^th and 75^th percentiles, whiskers: most extreme points that do not exceed ±IQR*1.5; further outliers are marked individually. (e-f) P: linear regression, one-sided p-value.

Figure 6.. An immunotherapy resistance MCP in melanoma tumors, linking T cell dysfunction to APOE repression in macrophages.

(a) Melanoma multicellular ICB resistance program (MCP2). Average expression (Z score, red/blue color bar) of top genes (columns) from the ICB resistance program, sorted by their pertaining cell type (top color bar), across samples (rows), with samples sorted by Overall Expression (right, Online Methods), and labeled by treatment status and ICB response (left color bar). (b) ICB program components across cell types. Off diagonal panels: Comparison of Overall Expression scores (y and x axes) for each cell component of MCP2 (rows and columns, labels on diagonal) across the samples (dots, black: responding samples, red: non-responding samples); the lines correspond to the linear fit. Pearson correlation coefficient (r) and significance (^***P < 1*10⁻³, one-sided) are shown in the upper triangle. Diagonal panels: Distribution of Overall Expression scores for each cell type component, along with the kernel density estimates. (c) Most genes in the ICB program are specific to one cell component. Number of genes (y axis) in the up (left bar) and down (right bar) compartment of the ICB program (x axis) that belong to one cell type component or multiple ones (color legend). (d) ICB multicellular program components are induced across cell types in ICB-resistant lesions. Distribution of Overall Expression scores (y axis) of ICB program components across the different cell types (y axis) for cells from non-responding (grey) and responding (light blue) lesions. P: One-sided p-values, mixed-effects models.