https://github.com/NathanWeinstein/Angiogenesis-Model
Tip revision: 0f67e1a0ffdaaf4991b588b1e970566c4b01b1f4 authored by NathanWeinstein on 08 November 2017, 17:38:49 UTC
Update ECexplorer.R
Update ECexplorer.R
Tip revision: 0f67e1a
RobustnessAnalisis.R
# R CMD BATCH MutantECexplorer.R
# Execute from R by typing "source("MutantECexplorer.R")"
# We will use BoolNet
if(!require(BoolNet)){
install.packages("BoolNet")
}
library(BoolNet)
source("matrixPlotter.R")
source("ECexplorerTools.R")
# A function to analize the robustness of the update rules of a network
getUpdateRuleRobustness <- function(net, nSamples = 1000){
genes <- net$genes
sensitivities <- c()
for(g in genes){
r <- perturbTrajectories(net, measure = "sensitivity", numSamples = nSamples, flipBits = 1, gene = g)
sensitivities <- append(sensitivities,r$value)
#print(genes[i])
#print(r$value)
}
names(sensitivities) <- genes
return(sensitivities)
}
# A function to calculate the robustness of Phalanx, Stalk and Tip ECs.
getECrobustness <- function(net, numSamples = 10000, flipBits = 1){
tip <- 0
stalk <- 0
phalanx <- 0
ctip <- 0
cstalk <- 0
cphalanx <- 0
ctotal <- 0
for(i in 1:numSamples){
# Generate a random state
si <- floor(runif(length(net$genes), 0, 1.999))
names(si) <- net$genes
# Generate a perturbed version of the random state
sii <- si
for(j in 1:flipBits){
pj <- floor(runif(1, 1, length(net$genes) + 0.999))
if(sii[pj] == 0){
sii[pj] <- 1
}else{
sii[pj] <- 0
}
}
# Get and sort the atractor of the random state and the perturbed state
attr <- getAttractor(net, si)
kind <- sortAtractor(analyzeAttractor(net, attr))
attrp <- getAttractor(net, sii)
kindp <- sortAtractor(analyzeAttractor(net, attrp))
# Store the result
if(kind == kindp){
ctotal <- ctotal + 1
}
if(kind == "Tip"){
tip <- tip + 1
if(kindp == "Tip"){
ctip <- ctip + 1
}
}else if(kind == "Stalk"){
stalk <- stalk + 1
if(kindp == "Stalk"){
cstalk <- cstalk + 1
}
}else if(kind == "Phalanx"){
phalanx <- phalanx + 1
if(kindp == "Phalanx"){
cphalanx <- cphalanx + 1
}
}
}
rtip <- ctip/tip
rstalk <- cstalk/stalk
rphalanx <- cphalanx/phalanx
rtotal <- ctotal/numSamples
return(c(rtip = rtip, rstalk = rstalk, rphalanx = rphalanx, rtotal = rtotal))
}
## First verify that the models are in the working directory
WD <- getwd()
simple_model <- paste(WD, "angiogenesis.txt", sep = "/")
if(file.exists(simple_model)){
}else{
print("The file containing the simplified model is not in the working directory")
}
# Load the simplified network from a file
net <- loadNetwork(simple_model)
# 1) Robustnes of attractor detemination to molecular activation noise
attractor_robustness <- perturbTrajectories(net, measure = "attractor",numSamples = 1000, flipBits = 1)
# 2) Robustness of the Tip, Stalk and Phalanx EC behavior
ECbehavior_robustness <- getECrobustness(net, numSamples = 1000, flipBits = 1)
# 3) Sensitivity of the components of the update function
update_rule_sensitivity <- getUpdateRuleRobustness(net, nSamples = 1000)
# Save the results
robustness <- list("attractor" = attractor_robustness, "ECbehavior" = ECbehavior_robustness, "sensitivities" = update_rule_sensitivity)
save(robustness, file = "robustness.Rdata")