Inferring pattern-driving intercellular flows from single-cell and spatial transcriptomics

Abstract

From single-cell RNA-sequencing (scRNA-seq) and spatial transcriptomics (ST), one can extract high-dimensional gene expression patterns that can be described by intercellular communication networks or decoupled gene modules. These two descriptions of information flow are often assumed to occur independently. However, intercellular communication drives directed flows of information that are mediated by intracellular gene modules, in turn triggering outflows of other signals. Methodologies to describe such intercellular flows are lacking. We present FlowSig, a method that infers communication-driven intercellular flows from scRNA-seq or ST data using graphical causal modeling and conditional independence. We benchmark FlowSig using newly generated experimental cortical organoid data and synthetic data generated from mathematical modeling. We demonstrate FlowSig's utility by applying it to various studies, showing that FlowSig can capture stimulation-induced changes to paracrine signaling in pancreatic islets, demonstrate shifts in intercellular flows due to increasing COVID-19 severity and reconstruct morphogen-driven activator-inhibitor patterns in mouse embryogenesis.

Fig. 1. Description of the FlowSig model.

a, We model intercellular flows to be directed from inflowing intercellular signals to GEMs that capture intracellular regulatory responses and that drive outflowing intercellular signals. FlowSig outputs an intercellular flow network describing directed edges from inflow signal variables (receptor gene expression weighted by the average expression of downstream TF gene set, R_i × TF_i), to GEMs (cell membership to latent GEM factors, GEM_i) to outflow variables (signal ligand gene expression, L_i). b, FlowSig uses additional perturbation data and pathway knowledge of immediate downstream TF targets to learn accurate intercellular flows resulting from cell–cell communication. c, From spatial transcriptomics data, we can infer the amount of inflowing signals received at each spatial location more accurately, enabling us to infer intercellular flows without additional perturbation data. FlowSig outputs an intercellular flow network describing directed edges from inflow signal variables (inferred amount of received signal ligand from COMMOT, rec. L_i) to spatially resolved GEMs (membership to GEMs, GEM_i), to outflow variable (ligand gene expression, L_i).

Fig. 2. Synthetic validation of FlowSig.

a–c, Causal diagrams representing activation (green arrow) or inhibition (red arrow) of unidirectional activation from inflowing SHH to outflowing BMP4 (a), SHH-inflow-driven patterning of NKX2.2, OLIG2, PAX6 and IRX3 (b) and competition between SHH inflow and BMP4 inflow to drive dorsoventral (dorsal (D), intermediate (I) and ventral (V)) patterning (c). d–f, The TPR and TNR of FlowSig output for a–c, respectively. We considered the effect of additional perturbation data and the effect of applying our biological flow model to constrain edges. In d–f, plots were generated over 500 simulations. Light blue boxes indicate the cases when total receptor expression (free plus bound receptor) was used as the inflow variable, while dark blue boxes indicate the cases when bound receptor expression was used as the inflow variable. Box plot whisker bounds are: d, minimum (TPR: control only 0.5, 0.5; control + perturbation: 0.5, 0.5; TNR: control only 0.56, 0.64; control + perturbation: 0.62, 0.72), maximum (TPR: control only 0.5, 0.5; control + perturbation: 0.5, 1.0; TNR: control only 0.79, 0.85; control + perturbation: 0.79, 0.85). Horizontal lines are defined by Q1 (TPR: control only 0.5, 0.5; control + perturbation: 0.5, 0.5; TNR: control only 0.56, 0.72; control + perturbation: 0.69, 0.77), median (TPR: control only 0.5, 0.5; control + perturbation: 0.5, 0.5; TNR: control only 0.69, 0.77; control + perturbation: 0.72, 0.77) and Q3 (TPR: control only 0.5, 0.5; control + perturbation: 0.5, 0.5; TNR: control only 0.77, 0.79; control + perturbation: 0.77, 0.79). e, Minimum (TPR: control only 0.18, 0.18; control + perturbation: 0.36, 0.45; TNR: control only 0.56, 0.60; control + perturbation: 0.56, 0.76), maximum (TPR: control only 0.73, 0.73; control + perturbation: 0.73, 0.73; TNR: control only 0.80, 0.92; control + perturbation: 0.76, 0.88). Horizontal lines are defined by Q1 (TPR: control only 0.55, 0.55; control + perturbation: 0.55, 0.55; TNR: control only 0.64, 0.76; control + perturbation: 0.76, 0.84), median (TPR: control only 0.55, 0.55; control + perturbation: 0.55, 0.55; TNR: control only 0.72, 0.84; control + perturbation: 0.76, 0.88) and Q3 (TPR: control only 0.55, 0.55; control + perturbation: 0.55, 0.55; TNR: control only 0.76, 0.88; control + perturbation: 0.76, 0.88). f, Minimum (TPR: control only 0.3, 0.3; control + perturbation: 0.4, 0.4; TNR: control only 0.56, 0.64; control + perturbation: 0.62, 0.72), maximum (TPR: control only 0.8, 0.8; control + perturbation: 0.8, 0.8; TNR: control only 0.79, 0.85; control + perturbation: 0.79, 0.85). Horizontal lines are defined by Q1 (TPR: control only 0.4, 0.4; control + perturbation: 0.4, 0.4; TNR: control only 0.67, 0.72; control + perturbation: 0.69, 0.77), median (TPR: control only 0.5, 0.5; control + perturbation: 0.5, 0.5; TNR: control only 0.69, 0.77; control + perturbation: 0.72, 0.77) and Q3 (TPR: control only 0.6, 0.6; control + perturbation: 0.6, 0.6; TNR: control only 0.77, 0.79; control + perturbation: 0.77, 0.79). Diamonds indicate outliers (less than Q1 – 1.5 × IQR or greater than Q3 + 1.5 × IQR).

Fig. 3. Experimental validation of FlowSig using a cortical organoid model.

a, Differentially inflowing signals between D18 and D35 timepoints. b, Differentially outflowing signals between D18 and D35 timepoints. c, Constructed gene expression modules capture time-specific and time-shared gene expression patterns. Some modules have been highlighted by the timepoint for which they are more enriched. d, Inferred intercellular flows due to FGF signaling. We speculated that EOMES is a key regulatory TF downstream of FGFR1 inflow. e, Inferred intercellular flows due to BMP signaling. We speculated that NR2F1 (CoupTF1) and PAX6 are downstream TF targets of BMP inflow. f, Measurement of the FC of EOMES gene expression using RT–qPCR following bath application of FGF, with four biological replicates and two technical replicates. Unpaired two-tailed t-tests were performed (t = 3.135, d.f. = 14, P = 0.0073). g, Measurement of the FC of PAX6 and NR2F1 (CoupTF1) gene expression using RT–qPCR following bath application of BMP, with two biological replicates and two technical replicates each. One-way ANOVA using Tukey’s multiple-comparisons test was used to calculate adjusted P values. For PAX6, control versus 10 ng ml^–1: mean diff. 0.50 with 95% confidence interval (CI) (0.40, 0.61), adjusted P < 1 × 10⁻⁴. For PAX6, control versus 50 ng ml^–1: mean diff. 0.57 with 95% CI (0.47, 0.68), adjusted P < 1 × 10⁻⁴. For NR2F1, control versus 10 ng ml^–1: mean diff. −0.91 with 95% CI (−1.88, 0.06), adjusted P = 0.066. For NR2F1 control versus 50 ng ml^–1: mean diff. −1.43 with 95% CI (−2.40, −0.45), adjusted P = 0.0068. In f, *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, ****P ≤ 0.0001. In g, *adjusted P ≤ 0.05, **adjusted P ≤ 0.01, ***adjusted P ≤ 0.001, ****adjusted P ≤ 0.0001. Error bars represent s.d.

Fig. 4. Application of FlowSig to perturbed non-spatial scRNA-seq of pancreatic islets.

a, Constructed GEM modules align primarily with cell types identified from clustering. b,c, Differentially inflowing (b) and outflowing (c) signals due to IFN-γ stimulation. d, Identified global intercellular flow network inferred by FlowSig, capturing condition-shared and condition-specific flows. e,f, Intercellular flows regulating outflowing signals upregulated (e) and downregulated (f) by IFN-γ stimulation.

Fig. 5. Application of FlowSig to scRNA-seq of human BALF sampled from people with moderate or severe COVID-19.

a,b, Constructed GEM modules align with both COVID-19 conditions (a) and cell types identified from the original study (b). NK, natural killer; mDC, myeloid dendritic cell; pDC, plasmacytoid dendritic cell. c, Differential expression analysis identifies distinct sets of outflowing signal ligands for each condition. d–f, Identified intercellular flows driving outflow signals differentially expressed in healthy controls (d) and individuals with moderate (e) or severe (f) COVID-19. g,h, Upset plots quantifying the inferred inflow signals (g) and mediating GEMs (h) that are shared across COVID-19 conditions. The vertical bars represent the sizes of the intersection sets (number of ligand-receptor interactions and GEMs, respectively), and the horizontal bars represent the number of inflow signals (g) and GEMs (h) in each condition. In (g), some intersection sets are annotated with the constituent ligand-receptor interactions, of the form L–R, where L is the ligand and R is the receptor. In the case where multiple ligands competitively bind to the same receptor, different ligands are separated by a slash (/). For example, the interaction L₁/L₂–R denotes that L₁ and L₂ are the ligands that can separately bind to the receptor, R.

Fig. 6. Application of FlowSig to spatial Stereo-seq data of E9.5 mouse embryo.

a, Twenty identified spatial GEMs. Spatial spots are labeled by the GEM to which the spot membership is highest. b, Inflowing signals that drive Shh outflow. c, Identified downstream outflowing signals that are driven by inflowing Shh (r-Shh). d, Top upstream TFs ranked by their regulatory effects on outflowing Shh, as measured by random forest feature (Gini) importance. Feature importances were calculated over 10 runs. Data are shown as mean ± s.d. e, Potential downstream TFs regulated by inflowing r-Shh, ranked by Spearman correlation. The heatmap shows the scaled (z-transformed) values of the fitted gene expression values as a function of r-Shh. f, The top upstream TFs of outflowing Wnt5a that are also downstream targets of inflowing r-Shh. Feature importance was averaged over 10 runs. Data are shown as mean ± s.d. g, Potential downstream TFs regulated by Wnt5a that are also upstream regulators of outflowing Shh. h, Suggested activator–inhibitor feedback between Shh and Wnt5a, as implicated by d–g.