Dynamic gene regulatory network inference from single-cell data using optimal transport

Abstract

Motivation: Modelling gene expression is a central problem in systems biology. Single-cell technologies have revolutionized the field by enabling sequencing at the resolution of individual cells. This results in a much richer data compared to what is obtained by bulk technologies, offering new possibilities and challenges for gene regulatory network inference.

Results: In this work, we introduce GRIT (gene regulation inference by transport)-a method to fit a differential equation model and to infer gene regulatory networks from single-cell data using the theory of optimal transport. The idea consists in tracking the evolution of the cell distribution over time and finding the system whose temporal marginals minimize the transport cost with the observations. GRIT is finally used to identify genes and pathways affected by two Parkinson's disease associated mutations.

Availability and implementation: Matlab implementation of the method and code for data generation are at gitlab.com/uniluxembourg/lcsb/systems-control/grit together with a user guide. A snapshot of the code used for the results of this article is at doi: 10.5281/zenodo.15582432.

Figure 1.

Method pipeline: The model identification iteration is run until convergence. The iteration can be interpreted as a block coordinate descent for solving (3), alternating between minimization with respect to the transport plans $M_k$ and the model $(A,b)$. The results are then used in a separate variable selection step to obtain confidence scores for GRN link existence.

Figure 2.

(a) The squared model error $\parallel [A,b]-[\overline{A},\overline{b}]{\parallel }_F^2$ with data simulated from a linear discrete-time system with varying number of timepoints and number of cells per timepoint. The plot is in logarithmic scale and it shows the mean and 80th percentiles obtained from five replicates. The sloped gridlines correspond to a decay $m^{-1/2}$. (b) The AUROC and AUPR scores from 20 replicates with data simulated from linear discrete-time system (LD), linear continuous-time system (LC), nonlinear system (NL), or system with non-observed protein concentrations (Pr).

Figure 3.

Results summary on the BEELINE benchmarking pipeline for (a) the synthetic dataset, (b) curated dataset, and (c) the RNA-Seq data. The bars show the average rank of each method in the different tasks, and the error bars indicate the means of the bottom 50th percentile and the top 50th percentile.

Figure 4.

Histogram of perturbation target scores with high-confidence genes indicated for the LRRK2 mutation (a) and for the PINK1 mutation (b).