Screening cell-cell communication in spatial transcriptomics via collective optimal transport

Abstract

Spatial transcriptomic technologies and spatially annotated single-cell RNA sequencing datasets provide unprecedented opportunities to dissect cell-cell communication (CCC). However, incorporation of the spatial information and complex biochemical processes required in the reconstruction of CCC remains a major challenge. Here, we present COMMOT (COMMunication analysis by Optimal Transport) to infer CCC in spatial transcriptomics, which accounts for the competition between different ligand and receptor species as well as spatial distances between cells. A collective optimal transport method is developed to handle complex molecular interactions and spatial constraints. Furthermore, we introduce downstream analysis tools to infer spatial signaling directionality and genes regulated by signaling using machine learning models. We apply COMMOT to simulation data and eight spatial datasets acquired with five different technologies to show its effectiveness and robustness in identifying spatial CCC in data with varying spatial resolutions and gene coverages. Finally, COMMOT identifies new CCCs during skin morphogenesis in a case study of human epidermal development.

Fig. 1. Overview of COMMOT.

a, COMMOT infers CCC in space while considering the competition between different ligand and receptor species. b, Collective optimal transport (COT) infers CCC in space by introducing multi-species distributions and enforcing limited spatial ranges. c, An example of inferring CCC for spatial distributions of ligand–receptor complexes from spatial distributions of the ligands and receptor where two ligand species (L1, L2) compete for one receptor species (R). d, Three applications of downstream analysis based on the inferred CCC network between cells or spots. DEG, differentially expressed gene; dir., direction; w.r.t., with respect to.

Fig. 2. Role of CCC in human skin development.

a, Predicted spatial origin of the skin subtypes of cells in intact tissue and the pseudotime projected to space. GRN, granular cell cluster; SPN, spinous cell cluster. b,c, The inferred amount of received signals of two example ligand–receptor pairs, GAS6-TYRO3 and PROS1-TYRO3 at the cell level (b) and cluster level (c). d, Immunostaining of proteins for GAS6, TYRO3 and PROS1. e, Fluorescent in situ hybridization against RNA molecules for predicted ligand–receptor interactions in human epidermis (solid white outline; regions of interest are marked by a white dashed square). The top row shows expression patterns of GAS6 (white) and TYRO3 (green); the bottom row shows expression patterns for PROS1 (white) and TYRO3 (green). In both cases, the middle and right panels show ligand–receptor signals, some of which colocalize to the stratum granulosum (white arrowheads). In merged images, the brightness of the GAS6 channel was increased to improve clarity against the prominent TYRO3 (green) signal. Experiments were repeated four times independently with consistent results. f, The signaling directions of four major signaling pathways. g, Heatmaps of selected signaling differentially expressed genes of the four signaling pathways, respectively. h, Immunofluorescence staining images of the identified signaling differentially expressed genes supporting the identified correlation between WNT signaling and the expression of these genes. Scale bars: d,e,h, 100 μm. The immunostaining experiments in d and h were repeated three times independently with consistent results.

Fig. 3. Inference of signaling direction in single-cell resolution spatial transcriptomics data.

a, MERFISH data of the mouse hypothalamic preoptic region with multiple slices across the anterior–posterior axis. b, Cluster-level summary of CCC through the OXT signaling pathway. c, Signaling directions of the OXT pathway. d, STARmap data of the mouse placenta. e, Signaling directions of the midkine, IGF, annexin and angiopoietin pathways.

Fig. 4. Downstream analysis of inferred CCC in single-cell resolution spatial transcriptomics data.

a–e, CCC analysis of seqFISH+ data of mouse secondary somatosensory cortex. a, Clustering of cell type based on gene expression. OPC, oligodendrocyte precursor cells. b, Enriched signaling in each cell type. c, Clustering based on inferred CCC. d, Enriched signaling in CCC-induced clusters. e, Differentially expressed genes in the CCC-induced clusters. f–h, CCC analysis of Slide-seq (v2) data of mouse hippocampus. f, Clustering of cell type based on gene expression. g, Clustering based on inferred CCC. h, Enriched signaling in CCC-induced clusters.

Fig. 5. CCC inference using Visium spatial transcriptomics data.

a–e, Midkine (MK) signaling in human breast cancer tissue. a, Spatial signaling direction. b, Amount of received signal by each spot. c, Two examples of differentially expressed (DE) genes due to signaling. d, Identification of the differentially expressed genes due to the total amount of received signal in the MK signaling pathway. e, Unique impact on the identified differentially expressed genes by the individual ligand–receptor pairs. f,g, Signaling in mouse brain tissue. The signaling direction (left) and the level of received signal (right) are shown for PSAP signaling (f) and FGF signaling (g).

Extended Data Fig. 1. Validation using simulated data by partial differential equations (PDE) model.

The example PDE model where two ligand species can bind to the same receptor. The inference by COMMOT is compared to the simulation results in several 1-dimensional cases. b Comparison to simulated results in a 2-dimensional case with three ligand species and two receptor species. c An example of randomly generated 2-dimensional benchmark with two ligand species that binds to the same receptor. The simulated result, inference by COMMOT, and inference by pairwise method are shown. d Ten different cases of ligand–receptor binding and the performance of COMMOT and pairwise OT (with the same spatial limit as COMMOT but each LR pair examined separately) obtained by comparing to simulated results.

Extended Data Fig. 2. OXT CCC in MERFISH mouse hypothalamic preoptic region.

The inferred signaling directions and cluster-level CCC of OXT signaling in each of the slice of the MERFISH data.

Extended Data Fig. 3. AGT signaling pathway in mouse cortex.

1) Cell type plots, 2) spatial directions of CCC, and 3) heatmaps of cluster-level CCC of the AGT signaling pathway in a Visium, b STARmap, and c seqFISH+ mouse cortex data. Across these three datasets, AGT signaling was identified in neurons. Spatially, neurons in the L2-3 region were identified as strong receivers of AGT ligands across the three datasets. Interestingly, a striped signaling pattern was observed, wherein strong signals within individual layers form stripes, while weak signals form inter-stripe regions. Strong AGT signaling activity among oligodendrocytes was also identified in both STARmap and seqFISH+ datasets.

Extended Data Fig. 4. WNT signaling pathway in mouse cortex.

1) Cell type plots, 2) spatial directions of CCC, and 3) heatmaps of cluster-level CCC of the WNT signaling pathway in a Visium and b seqFISH+ mouse cortex data. In both Visium and seqFISH+ cortex datasets, we inferred WNT signaling to be active across different cortical layers. In both datasets, we identified WNT signaling to be relatively low in layer 5, compared to other layers.

Extended Data Fig. 5. TAC signaling pathway in mouse cortex.

1) Cell type plots, 2) spatial directions of CCC, and 3) heatmaps of cluster-level CCC of the TAC signaling pathway in a Visium and b STARmap mouse cortex data. TAC (tachykinin neuropeptide family) signaling activity was consistently found in both Visium and STARmap cortex datasets to be active in non-neuronal cells and in inhibitory neurons, especially in somatostatin-expressing neurons (Sst).

Extended Data Fig. 6. Robustness of CCC analysis on a well-studied drosophila embryo dataset.

a Spatial signaling direction and signaling among cell clusters for Dpp and Wg signaling pathways. b Robustness of inferred signaling direction evaluated by comparing the direction obtained from subsampled dataset to the one from the full dataset using cosine distance. Each point is an independent test and the line shows the average of the tests. c Robustness of inferred cluster-level communication evaluated by comparing random subsamples to the full dataset using the Jaccard distance. d Robustness of downstream gene identification. e Percentage of known downstream genes that are identified as differentially expressed gene due to signaling activity. f Examples of the identified positively, negatively, and partially differentially expressed genes associated to Dpp signaling. For panels b–e, the averages of 5 independent random subsampling are plotted.