Inferring signaling pathways with probabilistic programming

Abstract

Motivation: Cells regulate themselves via dizzyingly complex biochemical processes called signaling pathways. These are usually depicted as a network, where nodes represent proteins and edges indicate their influence on each other. In order to understand diseases and therapies at the cellular level, it is crucial to have an accurate understanding of the signaling pathways at work. Since signaling pathways can be modified by disease, the ability to infer signaling pathways from condition- or patient-specific data is highly valuable. A variety of techniques exist for inferring signaling pathways. We build on past works that formulate signaling pathway inference as a Dynamic Bayesian Network structure estimation problem on phosphoproteomic time course data. We take a Bayesian approach, using Markov Chain Monte Carlo to estimate a posterior distribution over possible Dynamic Bayesian Network structures. Our primary contributions are (i) a novel proposal distribution that efficiently samples sparse graphs and (ii) the relaxation of common restrictive modeling assumptions.

Results: We implement our method, named Sparse Signaling Pathway Sampling, in Julia using the Gen probabilistic programming language. Probabilistic programming is a powerful methodology for building statistical models. The resulting code is modular, extensible and legible. The Gen language, in particular, allows us to customize our inference procedure for biological graphs and ensure efficient sampling. We evaluate our algorithm on simulated data and the HPN-DREAM pathway reconstruction challenge, comparing our performance against a variety of baseline methods. Our results demonstrate the vast potential for probabilistic programming, and Gen specifically, for biological network inference.

Availability and implementation: Find the full codebase at https://github.com/gitter-lab/ssps.

Supplementary information: Supplementary data are available at Bioinformatics online.

${\lambda }_j\sim \text{Uniform} ({\lambda }_{\text{min}},{\lambda }_{\text{max}}) \forall j\in \{1\dots |V\left|\right\}$ $z_{ij} | c_{ij},{\lambda }_j\sim \text{Bernoulli}\left(\frac{e^{-{\lambda }_j}}{e^{-c_{ij}{\lambda }_j}+e^{-{\lambda }_j}}\right) \forall i,j\in \{1\dots \left|V\right|\}$ ${\sigma }_j^2\propto \frac{1}{{\sigma }_j^2} \forall j\in \{1\dots |V\left|\right\}$ ${\beta }_j | {\sigma }_j^2\sim \mathcal{N}(0,T{\sigma }_j^2{(B_j^{\mathrm{\top }}{B}_j)}^{-1}) \forall j\in \{1\dots \left|V\right|\}$ $X_{+,j} | B_j,{\beta }_j,{\sigma }_j^2\sim \mathcal{N}(B_j{\beta }_j,{\sigma }_j^{2}I) \forall j\in \{1\dots \left|V\right|\}$ Fig. 1. Our generative model. (Top) Plate notation. DBN parameters β_j and ${\sigma }_j^2$ have been marginalized out. (Bottom) Full probabilistic specification. We usually set ${\lambda }_{\text{min}}\simeq 3$ and ${\lambda }_{\text{max}}=15$. If ${\lambda }_{\text{min}}>0$ is too small, Markov chains will occasionally be initialized with very large numbers of edges, causing computational issues. The method is insensitive to ${\lambda }_{\text{max}}$ as long as it is sufficiently large. Notice the improper prior $1/{\sigma }_j^2$. In this specification, B_j denotes $X_{-,{\text{pa}}_Z(j)}$; i.e. the parents of vertex j depend on edge existence variables Z

Fig. 2.

Action probabilities as a function of parent set size. The reference size $\widehat{s}$ is determined from prior knowledge. It approximates the size of a ‘typical’ parent set. When $s<\widehat{s}$, add-parent is most probable; when $s>\widehat{s}$, remove-parent is most probable; and when $s=\widehat{s}$, all actions have equal probability

Fig. 3.

Heatmap of AUCPR values from the simulation study. Both DBN-based techniques (SSPS and the exact method) score well on this, since the data are generated by a DBN. On large problems the exact DBN method needs strict in-degree constraints, leading to poor prediction quality. LASSO and FunChisq both perform relatively weakly. See Supplementary Figure S2 for representative ROC and precision-recall curves

Fig. 4.

Heatmap of differential performance against the prior knowledge, measured by AUCPR paired t-statistics. SSPS consistently outperforms the prior knowledge across problem sizes and shows robustness to errors in the prior knowledge

Fig. 5.

Methods’ performances across contexts in the HPN-DREAM Challenge. MCMC is stochastic, so we run SSPS 5 times; the error bars show the range of AUCROC scores. The other methods are all deterministic and require no error bars. See Supplementary Figure S3 for example predicted networks, Supplementary Figure S4 for AUCPR scores and Supplementary Figure S5 for representative ROC and precision-recall curves