FlowSig

Python package to infer directed intercellular flows described by ligand-receptor interactions driving tissue-scale gene expression patterning.

FlowSig requires:

1. Single-cell RNA-sequencing (scRNA-seq) data with cell state annotations that compare a baseline control to one or more perturbed conditions, e.g. disease severity, OR spatial transcriptomics (ST) data.
2. Cell-cell communication (CCC) inferred for each condition of interest. For non-spatial data, we require input from CellChat. For ST data, we require input from COMMOT.

The code used to generate all of the results in Almet et al., "Inferring pattern-driving intercellular flows from single-cell and spatial transcriptomics", can be found in another GitHub repository here.

Update (March 3, 2026): I've edited the implementation of NSF provided by MOFA-FLEX to calculate the spatial gene expression models instead of the original spatial-factorization. This should also enable FlowSig to work for multiple spatial samples.

<details> <summary>Installation</summary>

Installation instructions

The easiest way to currently install FlowSig is to generate a Python virtual environment and then run pip.

N.B. make sure you're using at least Python 3.11!

# Create the virtual environment
python3 -m venv flowsigenv

# Activate the virtual environment
source flowsigenv/bin/activate

First approach: pip install

pip install flowsig

Backup approoach: install directly from GitHub repository

# Clone the repository
git clone https://github.com/axelalmet/flowsig.git
cd ./flowsig/

# Install
pip3 install .

</details> <details> <summary>Non-spatial scRNA-seq example </summary>

Application to non-spatial scRNA-seq of stimulated pancreatic islets

Here, we show how to apply FlowSig to an scRNA-seq dataset of wildtype and stimulated human pancreatic islets, as originally studied in Burkhardt et al. (2021). The processed data and cell-cell communication inference, which was obtained using CellChat, can be downloaded from the following Zenodo repository.

You can also look at the code in a Jupyter notebook found here.

Import packages

import flowsig as fs
import scanpy as sc
import pandas as pd

Load the data and cell-cell communication inference

Data is specified in the form of a Scanpy object, which is really just an annotated dataframe, i.e. AnnData object. All subsequent output generated from FlowSig is stored in the Scanpy object.

data_directory = '../data/'

# Load the scanpy object
adata = sc.read(data_directory + 'burkhardt21_merged.h5ad')
condition_key = 'Condition'

# Load the cell-cell communication inference
cellchat_Ctrl = pd.read_csv('../communication_inference/output/burkhardt21_leiden_communications_Ctrl.csv')
cellchat_IFNg = pd.read_csv('../communication_inference/output/burkhardt21_leiden_communications_IFNg.csv')

cellchat_output_key = 'cellchat_output'
# Make sure your keys for the cellchat output dictionary match the relevant condition labels
adata.uns[cellchat_output_key] = {'Ctrl': cellchat_Ctrl,
                                  'IFNg': cellchat_IFNg}

Construct GEMs

We now construct gene expression modules (GEMs) from the unnormalised count data. For non-spatial scRNA-seq where we have multiple conditions, we use the iNMF algorithm by pyliger.

fs.pp.construct_gems_using_pyliger(adata,
                                n_gems = 10,
                                layer_key = 'counts',
                                condition_key = condition_key)

Construct the flow expression matrices

We construct augmented flow expression matrices for each condition that measure three types of variables:

1. Intercellular signal inflow, i.e., how much of a signal did a cell receive. For non-spatial scRNA-seq, signal inflow is defined as receptor gene expression weighted by the average expression of immediate downstream transcription factors that indicate signal activation.
2. GEMs, which encapsulate intracellular information processing. We define these as cellwise membership to the GEM.
3. Intercellular signal outflow, i.e., how much of a signal did a cell send. These are simply ligand gene expression.

The kay assumption of flowsig is that all intercellular information flows are directed from signal inflows to GEMs, from one GEM to another GEM, and from GEMs to signal outflows.

For non-spatial scRNA-seq, we need to specify the model organism so that FlowSig knows which receptor-transcription factor targets list to look at.

fs.pp.construct_flow_expressions(adata,
                                cellchat_output_key=cellchat_output_key,
                                model_organism = 'human',
                                spatial = False,
                                method = 'cellchat'
                                )

To reduce the number of variables over which we have to infer intercellular flows—and thus computation time—and to prioritise 'informative variables', we only retain inflow and outflow signals that are sufficiently differentially flowing between the control and perturbed conditions. We determine differentially flowing signals using a Wilcoxon rank-sum test and retain variables only if they are below a specified adjusted p-value threshold (q-value) and above a specified log-fold-change threshold.

fs.pp.ddetermine_informative_variables(adata,  
                                    spatial = False,
                                    condition_key = 'Condition',
                                    control = 'Ctrl',
                                    qval_threshold = 0.05,
                                    logfc_threshold = 0.5)

If you wanted to visualise which variables remained, e.g., which are differentially outflowing, you can run the following code:

fig, ax = plt.subplots(figsize=(6, 4))
fs.pl.plot_differentially_flowing_signals(adata,
                                        condition_key = 'Condition',
                                        pert_key = 'IFNg',
                                        var_type = 'outflow',
                                        flowsig_expr_key = 'X_flow_orig',
                                        flowsig_network_key = 'flowsig_network_orig',
                                        qval_threshold = 0.05,
                                        logfc_threshold = 0.5,
                                        label_lowqval = True,
                                        ax=ax
                                        )
plt.show()

Setting pert_key='IFNg orients results so that those upregulated in the perturbed condition have positive log fold change. Note also that you also need to set flowsig_expr_key = 'X_flow_orig' and flowsig_network_key = 'flowsig_network_orig' to visualise both significant and non-significant variables.

Learn intercellular flows

We are now in a position to learn the intercellular flows. To increase reliability of objects, we bootstrap aggregate results over a number of realisations. For non-spatial data, we have to specify the condition label and the control condition.

This step uses UT-IGSP to learn what is called a completed partially directed acyclic graph (CPDAG), which encodes directed arcs and undirected edges that describe the Markov Equivalence Class of statistical dependence relations that were learned directly from the data using conditional independence testing (how do variables depend on one another) and conditional invariance testing (which variables changed significantly between conditions). For both tests, we use a parametric partial-correlation-based method. The main reason we used these tests were because they take the least time to run compared to nonparametric kernel-based tests. Any test like the Hilbert-Schmidt Independence Criterion takes way too long for even 10-20 variables. The big caveat is that partial correlation assumes the data is described by a linear Gaussian model, which obviously isn't true for scRNA-seq. It's a long-term goal to add different types of nonparametric conditional independence/invariance tests that can be run in a reasonable amount of time.

fs.tl.learn_intercellular_flows(adata,
                        condition_key = condition_key,
                        control = 'Ctrl', 
                        use_spatial = False,
                        n_jobs = 4,
                        n_bootstraps = 500)

Partially validate intercellular flow network

Finally, we will remove any "false positive" edges. Noting that the CPDAG contains directed arcs and undirected arcs we do two things.

First, we remove directed arcs that are not oriented from signal inflow to GEM, GEM to GEM, or from GEM to signal outflow and for undirected edges, we reorient them so that they obey the previous directionalities.

fs.tl.apply_biological_flow(adata,
                            adjacency_key = 'adjacency',
                            validated_key = 'validated')

Second, we will remove directed arcs whose bootstrapped frequencies are below a specified edge threshold as well as undirected edges whose total bootstrapped frequencies are below the same threshold.

edge_threshold = 0.7
fs.tl.filter_low_confidence_edges(adata,
                                edge_threshold = edge_t