CellNOpt Documentation Center (CNODocs) (0.1.6)

3.1. Example combining SBMLqual, CNORdt, CNORode, CNORfeeder packages

This tutorial is for advanced users. It combines the CellNOptR, CNORode and CNORfeeder packages. The goal of this tutorial is to provide the material to reproduce some results shown or discussed about in the SBMLqual paper, which is available in the arxiv.

The first part reproduces one of figure shown in the paper. This figure was reproduced by the different software used in the paper using the same SBMLqual model as input. Here, we show how one can use CellNOptR oftware to simulate the same dynamical states as the one shown in the paper.

The second part is more about CellNOptR software itself: we deal with the creation of some data sets (with CNORode), optimisation (with CNORdt) and usage of CNORfeeder to find missing links.

Warning

In order to reproduce the examples and results provided here below, you must use some R packages available on cellnopt website rather than the current version available on BioConductor website, which are currently older than latest version required to reproduce the results. The examples were generated with these package/version:

  • CellNOptR “1.8.2”
  • CNORdt “1.5.1”
  • CNORfeeder “1.2.1”
  • CNORode “1.3.3”

3.1.1. Reading a SBMLqual and simulate using CNORdt

In CellNOptR, it is possible to load a SBMLqual file since version 1.7.8. Given the model, one can then decide to use different formalisms to simulate the dynamical states.

In this section, we provide a script called 2013/sbml/simulate.R that loads a SBMLqual and performs a simulation using another package called CNORdt. This package provide a discrete time formalism to simulate the dynamical states.

You can download the script from the link above and execute it as follows:

R script --no-save --no-restore < simulate.R

The SBMLqual file can be downloaded here: 2013/sbml/data/ModelV5.xml. The file is expected to be found in the directory: ./data. Otherwise you can change the script.

Each combination of stimuli lead to a different dynamical states. There are 4 combinations and therefore 4 pictures are generated. Here is one of them for which TNFa and EGF are stimulated:

_images/simulation_egf1_tnfa1.png

3.1.2. Creating data from SBMLQual and optimising a PKN on the data set.

In this section, we will do the following:

  1. Using the same SBMLQual model as above, which we will call True model, we first generate the dynamical states using the ODE formalism (CNORode). The states are saved in MIDAS format.
  2. Perform an optimisation (using CNORdt formalism) of the True model against the data for sanity check (the optimised model should be identical to the True model !!) and they fit almost perfect.
  3. Then, we will create a Prior Knowledge Network (PKN) that is slightly different from the True model.
  4. Using the PKN model and the data, we optimise the PKN model using a CNORdt formalism. We will see that there are discrepancies due to the differences between the True model and the PKN, in particular that there is a missing link in the PKN.
  5. Finally, we will use CNORfeeder to infer the missing link based on the data. An optimisation will show that the CNORfeeder found the missing link.

This example is based on the PKN and True model available in http://iopscience.iop.org/1478-3975/9/4/045003

Here is the code to generate the image corresponding to the True model:

library(CellNOptR)
model = readSBMLQual("ModelV5.xml")

library(CNORdt)
data(CNOlistPB, package="CNORdt")

plotModel(model, CNOlistPB)

Note

the CNOlistPB variable is optional and used solely for setting colors on nodes

Note

In the network plot, you can see that (1) there are 2 AND gates; they do not have any corresponding OR gates (2) there is a link from traf2 to mkk7.

_images/TrueModelSBML.png

3.1.2.1. Creating the data

To create the data using an ODE formalism, you will need a set of parameters available in 2013/sbml/data/paramsODE.RData.

Warning

if the SBML is changed, the order of the species may be different an you would need to map the parameters stored in the RData file above with your new SBML file.

library(CNORode) # for the simulation
library(CNORdt)  # for the data

data(CNOlistPB, package="CNORdt")

# Read the model with a dummy link on ph self loop
model = readSBMLQual("ModelV5.xml")

# load a variable called params that is required for the ODE run
load("paramsODE.RData")

timeSignals = seq(0,30,1)
CNOlistPB$timeSignals = seq(0,30,2)   # a trick to keep the original dt

# somehow stimuli set to exactly zero leads to NA so, replace them with
# tiny value but must set back to zero later on
CNOlistPB$valueStimuli[1,1] = 0.01
CNOlistPB$valueStimuli[1,2] = 0.01

sim_data = plotLBodeModelSim(CNOlistPB, model, ode_parameters=params,
    timeSignals=timeSignals, maxStepSize=0.035, show=F)
CNOlistPB$timeSignals = timeSignals

CNOlistPB$valueStimuli[1,1] = 0.
CNOlistPB$valueStimuli[1,2] = 0.

# finally converts the simdata to a cnolist to be saved in a file.
cnolist = simdata2cnolist(sim_data, CNOlistPB, model)

# keep only species of interest
species2select = c('raf1','erk', 'ap1', 'gsk3','p38', 'nfkb')
conditions2select = c(1,2,3,4,5,6,7,8,9,10)
for (i in seq_along(cnolist@signals)){
    cnolist@signals[i][[1]] =
cnolist@signals[i][[1]][conditions2select,species2select]
}
plot(cnolist)
_images/plotcnolist.png

Finally, let us save the new data set into a file:

writeMIDAS(cnolist, "data.csv", overwrite=T)

3.1.2.2. Checking

Let us first check that the model optimised against the data gives back the SBMLqual model.

library(CNORdt)
model = readSBMLQual("ModelV5.xml")
cnolist = CNOlist("data.csv")

# test that we can retrieve the exact data/model with itself
upperB = 10
lowerB= 0.8
boolUpdates = 30
bString=rep(1,30)

# optimised with itself to check and get the best RMSE
# Here we can set the maxtime to only 20 seconds
opt1 = gaBinaryDT(CNOlist=cnolist, model=model,
initBstring=rep(1,length(model$reacID)),
    boolUpdates=boolUpdates,maxTime=20, lowerB=lowerB, upperB=upperB,
    stallGenMax=20)

We can check that the best model is the model itself (all bits are on) and that the RMSE is small:

> opt1$bString
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> opt1$bScore
0.00142731

Finally, the following plot shows that errors across conditions and readouts are all small (white) by using this code:

cutAndPlotResultsDT(CNOlist=cnolist, model=model, upperB=upperB,
    lowerB=lowerB, bString=opt1$bString,    boolUpdates=boolUpdates)
_images/fit_check.png

3.1.2.3. Defining a PKN model

We use the original True model stored in the SBMLqual (image on top) but alter it in the following ways:

  1. We remove the 2 AND gates since we don’t know a priori the logical gates

    and add back the OR gates.

  2. We add the link tnfr–>pi3k

  3. We add the link pi3k–>rac–>map3k1

  4. We remove the link traf2->ask1->mkk7

The PKN is in SIF format available here 2013/sbml/PKN_noask1.sif

3.1.2.4. Optimising the PKN against data using CNORdt

First, we can check indeed that the True model can be used as a PKN to be optimised against the data generated. In theory, the final model should be identical to the True model.

pknmodel2 = readSIF("PKN_noask1.sif")
model = preprocessing(cnolist, pknmodel2, compression=F, expansion=T)
plotModel(model, cnolist)

# let us try this one where the additional links
# TNFR->PI3K and PI3K-RAC-Mapk31 are off as well as the unneeded AND and OR
gates.

# If you run the optimisation long enough, you'll get the optimised model to
# have the following bitstring:
optbs = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,0,0,0,1,0,0,0,1,0)

# It gives an RMSE = 0.031

opt1 = gaBinaryDT(CNOlist=cnolist, model=model,
    initBstring=optbs,boolUpdates=boolUpdates,maxTime=1000, lowerB=lowerB,
    upperB=upperB, stallGenMax=200, elitism=2, popSize=100, sizeFac=1e-4)

Finally, the following plot shows that errors across conditions and readouts are all small (white) by using this code:

cutAndPlotResultsDT(CNOlist=cnolist, model=model, upperB=upperB,
    lowerB=lowerB, bString=optbs,    boolUpdates=boolUpdates)

The fit is not good for the AP1 species. It means that

  1. The additional links that were added have been correctly removed (e.g., tnfr->pi3k)
  2. The OR gates added have been removed while the correct required AND gates (only) were kept
  3. The missing link on ASK1 prevents the AP1 to be fitted correctly when TNFa is on.
_images/fit_bad.png