CARNIVAL#
CARNIVAL (CAusal Reasoning for Network identification using Integer VALue programming) is a method for the identification of upstream reguatory signalling pathways from downstream gene expression (GEX). Applications of CARNIVAL include the identification of drug’s modes of action and of deregulated processes in diseases (even if the molecular targets remain unknown) by deciphering the alterations of main signalling pathways as well as alternative pathways and off-target effects.
Figure 1: Liu A., Trairatphisan P., Gjerga E. et al. From expression footprints to causal pathways: contextualizing large signaling networks with CARNIVAL npj Systems Biology and Applications volume 5, Article number: 40 (2019) (equal contributions).
The aim of the CARNIVAL pipeline is to identify a subset of interactions from a prior knowledge network that represent potential regulated pathways linking known or potential targets of perturbation towards active transcription factors derived from GEX data. The pipeline includes a number improved functionalities comparing to the original version and consists of the following processes:
Transcription factors’ (TFs) activities and pathway scores from gene expressions can be inferred with our in-house tools (Dorothea, CollecTRI).
TFs’ activities and signed directed protein-protein interaction networks with or without the provided target of perturbations and pathway scores are then used to construct an optimization problem with CORNETO.
CORNETO is used to solve the optimization problem with any of the supported solvers (CPLEX, GUROBI, SCIPY, etc), which identifies the sub-network topology with minimised fitting error and model size.
The original version of CARNIVAL was implemented in R and CPLEX. The new re-implementationo of CARNIVAL in CORNETO support a wide variety of solvers thanks to the support of both CVXPY and PICOS. It also has more flexibility since the problem is symbolically defined, and can be modified through the CORNETO API after creating the CARNIVAL problem. This gives user extra flexibility to modify the problem or to use CORNETO as a building block for other optimization problems.
import corneto as cn
cn.info()
|
|
A toy example#
G = cn.Graph.from_sif_tuples(
[
("I1", 1, "N1"), # I1 activates N1
("N1", 1, "M1"), # N1 activates M1
("N1", 1, "M2"), # N1 activaes M2
("I2", -1, "N2"), # I2 inhibits N2
("N2", -1, "M2"), # N2 inhibits M2
("N2", -1, "M1"), # N2 inhibits M1
]
)
G.plot()
from corneto.methods import runVanillaCarnival
# These are the measurements (e.g. TF activity from Decoupler).
# Positive values correspond to up-regulation and negative values
# with down-regulation. The bigger the absolute value is,
# the bigger the importance is
measurements = {"M1": 1, "M2": 1}
# Perturbations are the upstream nodes were the signal originates on,
# for example, ligands or receptors.
perturbations = {"I1": 1, "I2": 1}
# We run the `standard` carnival problem. This interface is similar
# to the old R function https://saezlab.github.io/CARNIVAL/reference/runVanillaCarnival.html
P, Gf = runVanillaCarnival(perturbations, measurements, G, betaWeight=0.1)
(CORNETO) Feb 08 11:47:15 AM - INFO : 2/2 inputs mapped to the graph
(CORNETO) Feb 08 11:47:15 AM - INFO : 2/2 outputs mapped to the graph
(CORNETO) Feb 08 11:47:15 AM - INFO : Pruning the graph with size: V x E = (6, 6)...
(CORNETO) Feb 08 11:47:15 AM - INFO : Finished. Final size: V x E = (6, 6).
(CORNETO) Feb 08 11:47:15 AM - INFO : 2/2 inputs after pruning.
(CORNETO) Feb 08 11:47:15 AM - INFO : 2/2 outputs after pruning.
(CORNETO) Feb 08 11:47:15 AM - INFO : Converting into a flow graph...
(CORNETO) Feb 08 11:47:15 AM - INFO : Creating a network flow problem...
(CORNETO) Feb 08 11:47:15 AM - INFO : Preprocess completed.
===============================================================================
CVXPY
v1.6.0
===============================================================================
(CVXPY) Feb 08 11:47:15 AM: Your problem has 86 variables, 197 constraints, and 0 parameters.
(CVXPY) Feb 08 11:47:15 AM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Feb 08 11:47:15 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Feb 08 11:47:15 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Feb 08 11:47:15 AM: Your problem is compiled with the CPP canonicalization backend.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Feb 08 11:47:15 AM: Compiling problem (target solver=SCIP).
(CVXPY) Feb 08 11:47:15 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCIP
(CVXPY) Feb 08 11:47:15 AM: Applying reduction Dcp2Cone
(CVXPY) Feb 08 11:47:15 AM: Applying reduction CvxAttr2Constr
(CVXPY) Feb 08 11:47:15 AM: Applying reduction ConeMatrixStuffing
(CVXPY) Feb 08 11:47:15 AM: Applying reduction SCIP
(CVXPY) Feb 08 11:47:15 AM: Finished problem compilation (took 2.937e-02 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Feb 08 11:47:15 AM: Invoking solver SCIP to obtain a solution.
presolving:
(round 1, fast) 36 del vars, 95 del conss, 0 add conss, 48 chg bounds, 1 chg sides, 23 chg coeffs, 0 upgd conss, 0 impls, 28 clqs
(round 2, fast) 54 del vars, 129 del conss, 0 add conss, 52 chg bounds, 9 chg sides, 31 chg coeffs, 0 upgd conss, 0 impls, 20 clqs
(round 3, fast) 63 del vars, 145 del conss, 0 add conss, 52 chg bounds, 11 chg sides, 33 chg coeffs, 0 upgd conss, 0 impls, 18 clqs
(round 4, fast) 65 del vars, 147 del conss, 0 add conss, 52 chg bounds, 15 chg sides, 37 chg coeffs, 0 upgd conss, 0 impls, 18 clqs
(round 5, fast) 67 del vars, 151 del conss, 0 add conss, 52 chg bounds, 15 chg sides, 37 chg coeffs, 0 upgd conss, 0 impls, 18 clqs
(round 6, exhaustive) 67 del vars, 151 del conss, 0 add conss, 52 chg bounds, 15 chg sides, 37 chg coeffs, 40 upgd conss, 0 impls, 18 clqs
(0.0s) probing cycle finished: starting next cycle
(round 7, exhaustive) 69 del vars, 151 del conss, 0 add conss, 52 chg bounds, 15 chg sides, 37 chg coeffs, 40 upgd conss, 28 impls, 46 clqs
(round 8, medium) 69 del vars, 155 del conss, 0 add conss, 52 chg bounds, 15 chg sides, 37 chg coeffs, 40 upgd conss, 28 impls, 46 clqs
(0.0s) probing cycle finished: starting next cycle
(0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -)
(0.0s) symmetry computation finished: 2 generators found (max: 1500, log10 of symmetry group size: 0.0) (symcode time: 0.00)
dynamic symmetry handling statistics:
orbitopal reduction: no components
orbital reduction: no components
lexicographic reduction: no permutations
handled 1 out of 1 symmetry components
(round 9, exhaustive) 69 del vars, 155 del conss, 3 add conss, 52 chg bounds, 15 chg sides, 37 chg coeffs, 40 upgd conss, 28 impls, 46 clqs
presolving (10 rounds: 10 fast, 5 medium, 4 exhaustive):
69 deleted vars, 155 deleted constraints, 3 added constraints, 52 tightened bounds, 0 added holes, 15 changed sides, 37 changed coefficients
28 implications, 46 cliques
presolved problem has 27 variables (19 bin, 0 int, 0 impl, 8 cont) and 45 constraints
16 constraints of type <varbound>
18 constraints of type <setppc>
9 constraints of type <linear>
2 constraints of type <logicor>
transformed objective value is always integral (scale: 0.1)
Presolving Time: 0.00
time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl.
p 0.0s| 1 | 0 | 1 | - | clique| 0 | 27 | 45 | 42 | 0 | 0 | 0 | 0 |-3.500000e+00 | 1.000000e-01 | Inf | unknown
p 0.0s| 1 | 0 | 1 | - | locks| 0 | 27 | 45 | 42 | 0 | 0 | 0 | 0 |-3.500000e+00 |-6.000000e-01 | 483.33%| unknown
p 0.0s| 1 | 0 | 1 | - | vbounds| 0 | 27 | 46 | 42 | 0 | 0 | 1 | 0 |-3.500000e+00 |-1.600000e+00 | 118.75%| unknown
p 0.0s| 1 | 0 | 1 | - | vbounds| 0 | 27 | 48 | 42 | 0 | 0 | 3 | 0 |-3.500000e+00 |-3.200000e+00 | 9.37%| unknown
0.0s| 1 | 0 | 25 | - | 1784k | 0 | 27 | 49 | 42 | 0 | 0 | 8 | 0 |-3.300000e+00 |-3.200000e+00 | 3.13%| unknown
r 0.0s| 1 | 0 | 25 | - |rounding| 0 | 27 | 49 | 39 | 0 | 0 | 8 | 0 |-3.300000e+00 |-3.300000e+00 | 0.00%| unknown
0.0s| 1 | 0 | 25 | - | 1784k | 0 | 27 | 49 | 39 | 0 | 0 | 8 | 0 |-3.300000e+00 |-3.300000e+00 | 0.00%| unknown
0.0s| 1 | 0 | 25 | - | 1784k | 0 | 27 | 49 | 39 | 0 | 0 | 8 | 0 |-3.300000e+00 |-3.300000e+00 | 0.00%| unknown
SCIP Status : problem is solved [optimal solution found]
Solving Time (sec) : 0.01
Solving Nodes : 1
Primal Bound : -3.30000000000000e+00 (5 solutions)
Dual Bound : -3.30000000000000e+00
Gap : 0.00 %
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Feb 08 11:47:15 AM: Problem status: optimal
(CVXPY) Feb 08 11:47:15 AM: Optimal value: 7.000e-01
(CVXPY) Feb 08 11:47:15 AM: Compilation took 2.937e-02 seconds
(CVXPY) Feb 08 11:47:15 AM: Solver (including time spent in interface) took 1.530e-02 seconds
(CORNETO) Feb 08 11:47:15 AM - INFO : Finished in 0.07 s.
P.objectives[0].value
array([0.])
P.objectives[1].value
np.float64(7.0)
# The method for plotting will be simplified in next versions
G.plot(
custom_edge_attr=cn.pl.edge_style(P, edge_var="edge_values_c0"),
custom_vertex_attr=cn.pl.vertex_style(P, Gf, vertex_var="vertex_values_c0"),
)
import pandas as pd
from corneto.methods.carnival import get_result, get_selected_edges
V, E = get_result(P, Gf)
pd.DataFrame(V)
| V | value | |
|---|---|---|
| 0 | M2 | 1.0 |
| 1 | N2 | -1.0 |
| 2 | N1 | 0.0 |
| 3 | I2 | 1.0 |
| 4 | M1 | 1.0 |
| 5 | I1 | 1.0 |
| 6 | _s | 1.0 |
| 7 | _pert_c0 | 1.0 |
| 8 | _meas_c0 | 0.0 |
| 9 | _t | 0.0 |
pd.DataFrame(E)
| E | value | |
|---|---|---|
| 0 | ((I1), (N1)) | 0.0 |
| 1 | ((N1), (M1)) | 0.0 |
| 2 | ((N1), (M2)) | 0.0 |
| 3 | ((I2), (N2)) | -1.0 |
| 4 | ((N2), (M2)) | 1.0 |
| 5 | ((N2), (M1)) | 1.0 |
| 6 | ((_s), (_pert_c0)) | 1.0 |
| 7 | ((_pert_c0), (I1)) | 1.0 |
| 8 | ((_pert_c0), (I2)) | 1.0 |
| 9 | ((M1), (_meas_c0)) | 1.0 |
| 10 | ((M2), (_meas_c0)) | 1.0 |
| 11 | ((_meas_c0), (_t)) | 0.0 |
| 12 | ((), (_s)) | 1.0 |
| 13 | ((_t), ()) | 0.0 |
G_sol = Gf.edge_subgraph(get_selected_edges(P, Gf))
G_sol.plot()
