STASNet
````r diagram="@startuml skinparam usecase { BackGroundColor<< S3 Class >> #eeeeee
BackGroundColor<< Helper Functions >> #88ee88 } (ModelSet) << S3 Class >> (MRAmodel) << S3 Class >> (profileLikelihoodResults) << List >> (getCombinationMatrix) << Helper Functions >> (importModel) --> MRAmodel (rebuildModel) --> MRAmodel (createModel)---->MRAmodel (MRAmodel)-->(plotModelPrediction) (MRAmodel)-->(plotModelAccuracy) (MRAmodel)-->(simulateModel) (MRAmodel)-->(plotModelScores) (createModelSet) --> ModelSet (ModelSet)-->(extractSubmodels) (ModelSet)-->(plotParameters) (ModelSet)-->(plotResiduals) (extractSubmodels) --> (MRAmodel) (simulateModel)-->(plotModelSimulation) (getCombinationMatrix)-->(simulateModel) (MRAmodel)-->(profileLikelihood) (profileLikelihood)-->(profileLikelihoodResults) (profileLikelihoodResults)-->(addPLinfos) (MRAmodel)-->(addPLinfos) (addPLinfos)-->(MRAmodel) (MRAmodel)-->(suggestExtension) (MRAmodel)-->(selectMinimalModel) @enduml"
write(diagram, "use_diagram.uml") system2("java", " -jar ~/bin/Java/plantuml.jar use_diagram.uml")
![Use diagram](use_diagram.png) **NOTE** Various scripts are available in the inst/ folder of `STASNet` to run the fitting procedure from the command line Loading the library: ```r library(STASNet)
# Structure as adjacency list structure = rbind(c("N1", "N2"), c("N2", "N3"), c("N2", "N5"), c("N5", "N4"), c("N3", "N4")) basal = c("N1", "N2", "N3", "N4", "N5") midas_data = as.data.frame(matrix(0, nrow=4, ncol=7, dimnames=list(NULL, c("ID:type", "TR:N1", "TR:N2i", "DA:ALL", "DV:N2", "DV:N3", "DV:N4")))) midas_data[,1] = c("c", "t", "t", "t") # Type of experiment midas_data[,c(2,3)] = cbind(c(0,1,1,0), c(0,0,1,1)) # Perturbations midas_data[,c(5:7)] = cbind(c(111, 250, 187, 98), c(100, 150, 147, 76), c(103, 600, 217, 54)) # Experimental values print(midas_data)
Create the model using the data and the network layout as prior knowledge.
NOTE Depending on the number of initialisations ('inits') and the size of the network it is possible that the global optimum is not reached. Results can thus vary between runs of createModel()
for big networks with an insufficient number of initialisations.
model = createModel(model_links = structure, basal_file = basal, data.stimulation = midas_data, inits = 100, model_name = "example")
#setwd(system.file("extdata", package="STASNet")) #model = createModel("network.tab", "basal.dat", "HCT116_MIDAS.csv", inits=100, perform_plots=T, # method="randomlhs"); #TODO: replace the HCT116 dataset with the trimed dataset used in the package
We can then visualise the resulting network. Be aware that the fitted values are paths and not directly link,
thus the values displayed on a link should not be taken as is (use printParameters()
to get the actual paths).
printParameters(model) plotModelGraph(model)
We can the evaluate how the simulation compares to the data.
plotModelAccuracy(model) plotModelScores(model) simulation = plotModelSimulation(model, with_data= TRUE,log_axis = T) print(simulation)
Profile likelihood aims at providing confidence intervals for the parameters. niplotPL()
provides a convenient way to plot the likelihood and profiles of parameters.
profiles = profileLikelihood(model, nb_points = 1000, nb_cores = 8) model = addPLinfos(model, profiles) niplotPL(profiles, file_plots = FALSE)
Note the correlation between N2->N3 and N2->N4->N5 which indicates a structural non identifiability, i.e. the experimental design does not allow to differentiate between them. Note also the symmetry around 0 for the inhibitor iN2 due to the fact that it is always considered negative for the fit.
The parameter confidence intervals can be visualised.
plotModelParameters(model)
We can then simulate with confidence intervals.
simulation = plotModelSimulation(model, with_data= TRUE,log_axis = T)
Next to the parameterisation also the structural sanity of the network can be assessed by searching for superfluous or missing connections.
First we try to reduce the network in oder to test for links that do not significantly contribute to the overall fit
red_model = selectMinimalModel(model)
In the model two links are found to be neglectable and the thusly reduced model is returned as output of selectMinimalModel()
After haveing reduce th model it is also of interest to search for missing links that would improve the overall fit. In contrast to the reduction this is done by returning a list which contains teh statistics of each link so that the extension can be then manually curated and compared to available biological insights from previous experiments or literature.
#TODO make ModelStructure to accept a list of all names to always preserve unconnected nodes in the network ext_mat = suggestExtension(red_model)
A model set is a joint fitting of several dataset with the same model and some links varying independently between models.
mset = createModelSet("network.tab", "basal.tab", c("Widr_"))
mset = addVariableParameters(mset)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.