Inferring signalling dynamics by integrating interventional with observational data

Abstract

Motivation: In order to infer a cell signalling network, we generally need interventional data from perturbation experiments. If the perturbation experiments are time-resolved, then signal progression through the network can be inferred. However, such designs are infeasible for large signalling networks, where it is more common to have steady-state perturbation data on the one hand, and a non-interventional time series on the other. Such was the design in a recent experiment investigating the coordination of epithelial-mesenchymal transition (EMT) in murine mammary gland cells. We aimed to infer the underlying signalling network of transcription factors and microRNAs coordinating EMT, as well as the signal progression during EMT.

Results: In the context of nested effects models, we developed a method for integrating perturbation data with a non-interventional time series. We applied the model to RNA sequencing data obtained from an EMT experiment. Part of the network inferred from RNA interference was validated experimentally using luciferase reporter assays. Our model extension is formulated as an integer linear programme, which can be solved efficiently using heuristic algorithms. This extension allowed us to infer the signal progression through the network during an EMT time course, and thereby assess when each regulator is necessary for EMT to advance.

Availability and implementation: R package at https://github.com/cbg-ethz/timeseriesNEM. The RNA sequencing data and microscopy images can be explored through a Shiny app at https://emt.bsse.ethz.ch.

Supplementary information: Supplementary data are available at Bioinformatics online.

Fig. 1.

Synopsis. (a) Normal cells undergo EMT induced by TGF-β. (b) EMT is halted in cells perturbed by RNAi against any of the 59 regulators under investigation. Despite TGF-β treatment to induce EMT, the perturbed cells do not reach the mesenchymal phenotype, but are halted at intermediate stages. (c) Nested effects model to infer a non-transcriptional signalling network from downstream perturbation effects on genes measured using RNA sequencing. (d) Inferring the signal progression during an unperturbed EMT time course

Fig. 2.

Nested effects model. The S-nodes are modelled as hidden variables, and we aim to infer their causal graph Φ (solid arrows). In experiments separately perturbing each S-node, we observe the differential expression of all E-genes. Assuming that each E-gene is directly regulated by at most one S-node in Φ, we compute the maximum a posteriori attachment Θ (dashed arrows) of effect genes to S-nodes. We search for the signalling graph Φ which yields the most likely probabilistic nesting of effects

Fig. 3.

Adjacency matrix showing the average of bootstrap and jackknife supports for each edge in the graph. Parent nodes are topologically ordered (from root to leaf), and child nodes follow the same ordering. The entry in row i and column j records the average of the bootstrap and jackknife supports for the directed edge $S_i\to S_j$

Fig. 4.

Transitive reduction of the consensus network, which contains edges whose average of bootstrap and jackknife supports exceeds 50%. Nodes are topologically ordered so that all edges point downward. Transcription factors are coloured blue and miRNAs orange

Fig. 5.

Summary of luciferase reporter assays for Rbpj, Smad4 and Sox4. The volcano plots are based on Welch’s t-tests of perturbation versus control, after removing batch effects using a linear model. Experiments above the dotted lines have FDR-adjusted P-values <5%, and are labelled. Predictions by the NEM consensus network (Fig. 4) are indicated by green triangles for a predicted effect, and orange circles for a prediction of no effect. Part of this validation dataset was previously published by Meyer-Schaller et al. (2019)

Fig. 6.

Principal component analysis based on expression of 24 454 genes in biologically duplicated samples taken during an EMT time course. Top panel: PC1 versus PC2, bottom: PC2 versus PC3. The data was normalized using a regularized logarithmic transform (Love et al., 2014), and batch effects were removed using a linear model (Ritchie et al., 2015)

Fig. 7.

Inferred signal progression at $t\in \{0,12,24,36,48,60,72\}$ hours after TGF-β treatment. An S-node is uncoloured if it is inferred to be in a different signalling state at t versus 96 h; otherwise it is coloured green. The colour gradient comes from aggregating the inferred activity states over the networks resulting from 1000 bootstrap and 59 jackknife resamples. For illustration purposes, the signal progression is mapped onto the consensus network (Fig. 4)

Fig. 8.

Computational complexity of our method for mapping an observational effect profile onto a static NEM. The CPU time in milliseconds is shown on a logarithmic scale on the y-axis. The number of S-nodes n varied from 10, 20, …, 100 (left panel). The number of E-genes m varied from 2000, 4000, …, 20 000 (right panel). Simulations were timed using the R package microbenchmark (Mersmann, 2018)