Nothing
R/readSBML.r
In SBMLR: SBML-R Interface and Analysis Tools
Defines functions
getRuleLeaves
ML2R
"readSBML" <-
function(filename)
{ # takes SBML in filename.xml and maps it to a SBML class model
# using both Sax and DOM (for mathml) based parsing.
sbmlHandler <- function ()
{ # first block here sets up the parent environment used by all handler functions
sbml<-"x" # "x" is just a starting string value
modelid<-"x" #storing model id
lnotes<-NULL
compartments <- list()
reactLaws <- list()
species <- list()
rules<- list()
reactions<- list()
globalParameters=list()
reactants=NULL
products=NULL
modifiers=NULL
currRxnID=NULL
parameters=NULL # local to rate law
parameterIDs=NULL # local to rate law
globalParameterIDs=NULL
notes=FALSE; reactant=FALSE; product=FALSE
law=FALSE; parameter=FALSE; math=FALSE
# # VV's additions
# modelname<-"x" #storing the model name if separately given
# paramException <- 1
# globalParamException <- 1
# ParametersList<- list() #list of parameter objects
.startElement <- function(name, atts, ...) {
# cat("Start: Name =",name," ",paste(names(atts),atts,sep=" = "),"\n")
if(name=="sbml") sbml<<-atts
if(name=="annotation") print("skipping annotation")
if(name=="model") {modelid<<-atts[["id"]]} # VV replaces this with ...
# if(name=="model")
# { numitems <- length(atts)
# if(numitems < 1) #if model does not contain a name/id, we give it an arbitrary one.
# modelid[[1]]<<-"BioModel"
# else if(numitems == 1) # if only one attribute supplied
# {
# if(is.character(atts[[1]])) #if the attribute is a string, it must be model name. Copy to also id.
# {
# modelname<<-atts[[1]] # store model name
# modelid<<-atts[[1]] # store as model id (copy Essentially)
# }
# }
# else if(numitems > 1) #both Id and name of model are supplied, read both
# {
# modelname<<-atts[["name"]] # first element is model name
# modelid <<-atts[["id"]] # second element has to be model id
# }
# }
# if(name=="compartment") {compartments[[atts[1]]]<<-atts}
if(name=="compartment")
if( "id" %in% names(atts)) compartments[[atts["id"]]]<<-atts
# if(name=="species") {species[[atts[1]]]<<-atts }
if(name=="species")
if( "id" %in% names(atts)) species[[atts["id"]]]<<-atts
# if(name=="assignmentRule") rules[[atts[1]]]$idOutput<<-atts[[1]]
if(name=="assignmentRule") rules[[atts["variable"]]]$idOutput<<-atts[["variable"]]
# if(name=="reaction") {reactions[[atts[1]]]$id<<-atts[[1]]
# reactions[[atts[1]]]$reversible<<-as.logical(atts[[2]])
# currRxnID<<-atts[1]}
if(name=="reaction")
{
lstnames <- names(atts)
numitems <- length(lstnames)
nameslist <- list()
id <- "x"
reverse <- FALSE
name <- "x"
count <- 1
while( count <= numitems )
{
switch(lstnames[[count]],
"id" = { id = atts[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"reversible" = { reverse = as.logical(atts[[count]]) ;nameslist[[length(nameslist)+1]] <- "reversible" },
"name" = { name = as.character(atts[[count]]); nameslist[[length(nameslist)+1]] <- "name"}
)
count <- count + 1
}
reactions[[atts["id"]]]$id<<-id
reactions[[atts["id"]]]$reversible<<-reverse
# if(reverse) reactions[[atts["id"]]]$reversible<<-reverse #carry in R only if true
currRxnID<<-atts["id"]
#reactions[[atts["id"]]]$id<<-atts[["id"]]
#reactions[[atts[1]]]$reversible<<-as.logical(atts[[2]])
#currRxnID<<-atts[1]
}
if(name=="listOfReactants") reactant<<-TRUE
if(name=="listOfProducts") product<<-TRUE
if(name=="kineticLaw") law<<-TRUE
if(name=="math") math<<-TRUE
if((name=="speciesReference")&reactant)
reactants<<-c(reactants,species=atts[["species"]])
# reactants<<-c(reactants,species=atts[[1]])
if((name=="speciesReference")&product)
products<<-c(products,species=atts[["species"]])
# products<<-c(products,species=atts[[1]])
if(name=="modifierSpeciesReference")
modifiers<<-c(modifiers,species=atts[["species"]])
if((name=="parameter")&law){
parameterIDs<<-c(parameterIDs,atts[["id"]])
parameters<<-c(parameters,atts[["value"]])}
# if((name=="parameter")&law) #parameter encountered within a kinetic law definition
# {
# values <- names(atts)
# if( "id" %in% values) parameterIDs<<-c(parameterIDs,atts[["id"]])
# else {
# cat('Parameter parsed without id. Setting default id', '\n')
# parameterIDs<<-c(parameterIDs, paste("default", paramException))
# paramException <- paramException + 1
# }
# if( "value" %in% values) parameters<<-c(parameters,atts[["value"]])
# else {
# cat('Warning..Parsing parameter without value. Setting it as 0.' ,'\n')
# parameters<<-c(parameters, as.numeric(0))
# }
# }
if((name=="parameter")&(!law)){
globalParameterIDs<<-c(globalParameterIDs,atts[["id"]])
globalParameters<<-c(globalParameters,as.numeric(atts[["value"]]))}
# if((name=="parameter")&!law) #parameter encountered outside a kinetic law definition - So in globalparamslist
# {
# #cat("within parameters:", atts[["id"]], atts[["value"]], "\n")
# values <- names(atts)
# if( "id" %in% values) {
# globalParameterIDs<<-c(globalParameterIDs,atts[["id"]])
# ParametersList[[atts["id"]]] <<- atts #our new list of Parameter Objects
# }
# else {
# cat('Global Parameter parsed without id. Setting default id', '\n')
# tempParamId <- paste("Globaldefault", globalParamException)
# globalParameterIDs<<-c(globalParameterIDs, tempParamId)
# ParametersList[[tempParamId]] <<- atts
# globalParamException <- globalParamException + 1
# }
# if( "value" %in% values) globalParameters<<-c(globalParameters,as.numeric(atts[["value"]]))
# else {
# cat('Warning..Parsing Global parameter without value. Setting it to 0.', '\n')
# globalParameters<<-c(globalParameters, as.numeric(0))
# }
# } # end if param in law
} # end .startElement()
.endElement <- function(name) {
if(name=="listOfReactants") reactant<<-FALSE
if(name=="listOfProducts") product<<-FALSE
if(name=="kineticLaw") law<<-FALSE
if(name=="math") math<<-FALSE
if((name=="listOfParameters")&(!law)) names(globalParameters)<<-globalParameterIDs
if(name=="reaction") {
names(reactants)<<-NULL
names(modifiers)<<-NULL
names(products)<<-NULL
reactions[[currRxnID]]$reactants<<-reactants
reactions[[currRxnID]]$modifiers<<-modifiers
reactions[[currRxnID]]$products<<-products
parameters<<-as.numeric(parameters)
names(parameters)<<-parameterIDs
reactions[[currRxnID]]$parameters<<-parameters
reactants<<-NULL;products<<-NULL
modifiers<<-NULL;parameters<<-NULL;parameterIDs<<-NULL }
}
.text <- function(x, ...) {
if (!math) lnotes<<-c(lnotes,x)
# cat("Txt:", x,"\n")
}
getModel<- function()
{
# fixComps=function(x) {
# lst=list(x[[1]],as.numeric(x[[2]]) );
# names(lst)<-names(x);
# lst
# }
# VV replaces fixComps with the following:
fixComps=function(x)
{
lstnames <- names(x)
count <- 1
numit <- length(lstnames)
id <- "x"
size <- 0
name <- "x"
nameslist <- list()
while( count <= numit )
{
switch(lstnames[[count]],
"id" = { id = x[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"size" = { size = as.numeric(x[[count]]) ;nameslist[[length(nameslist)+1]] <- "size" },
"name" = { name = as.character(x[[count]]); nameslist[[length(nameslist)+1]] <- "name"}
)
count = count + 1
}
#cat("Compartment id:", id,"\n", "Compartment size:" , size, "\n","Compartment name:", name, "\n")
if(numit == 2) # only 2 attributes present. We need to find them.
{ if(id == "x") #id not set but name and size are.
id <- "default"
else if(name == "x") #name not set, we copy the id.
name <- id
else if(size== "0") #size not set
size <- 1 #arbitrary setting as 1
lst = list(id,size,name)
names(lst)<-c("id","size","name")
lst
} else if(numit == 3) # 3 attributes/items present.
{ lst = list(id,size,name)
names(lst)<-c("id","size","name")
lst
}
}
# fixSpecies=function(x) {
# lst=list(x[[1]],as.numeric(x[[2]]),x[[3]],as.logical(x[[4]]));
# names(lst)<-c("id","ic","compartment","bc");
# lst
# }
# VV replaces fixSpecies with the following
fixSpecies=function(x)
{
#cat (names(x), "\n")
#cat(toString(x) , "\n")
numitems <- length(x)
lstnames <- names(x)
count <-1
id <- "x" #species Id
ic <- 0 #species initial concentration
compart <- "def" #species compartment
bc <- FALSE #species boundary condition
name <- "def"
nameslist <- list()
while( count <= numitems)
{
switch(lstnames[[count]],
"id" = { id <- x[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"name" = { name <- x[[count]]; nameslist[[length(nameslist)+1]] <- "name"},
"initialConcentration" = { ic <- as.numeric(x[[count]]) ;nameslist[[length(nameslist)+1]] <- "ic" },
"compartment" = { compart <- as.character(x[[count]]); nameslist[[length(nameslist)+1]] <- "compartment"},
"boundaryCondition" = { bc <- as.logical(x[[count]]); nameslist[[length(nameslist)+1]] <- "bc"}
)
count = count + 1
}
#lst = list(id,ic,compart,bc, name)
lst = list(id,as.numeric(ic), compart, as.logical(bc))
names(lst) <- c("id","ic","compartment","bc")
#names(lst)<-c("id","ic","compartment","bc", "name");
lst
}
# and VV adds in fixParams
fixParams=function(x)
{
numitems <- length(x)
lstnames <- names(x)
count <-1
id <- "x" #Parameter Id
value <- 0 #Parameter value
name <- "def"
constant <- FALSE
nameslist <- list()
while( count <= numitems)
{
switch(lstnames[[count]],
"id" = { id <- x[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"name" = { name <- x[[count]]; nameslist[[length(nameslist)+1]] <- "name"},
"value" = { value <- as.numeric(x[[count]]) ;nameslist[[length(nameslist)+1]] <- "value" },
"constant" = { constant <- as.logical(x[[count]]) ; nameslist[[length(nameslist)+1]] <- "constant"}
)
count = count + 1
}
lst = list(id,as.numeric(value))
names(lst) <- c("id","value")
lst
}
compartments=sapply(compartments,fixComps, simplify = FALSE)
#species=t(sapply(species,fixSpecies, simplify = TRUE)[2:4,]) # this changes the species model structure for better looks in R dumps
species=sapply(species,fixSpecies, simplify = FALSE) # this keeps the better looks in the SBMLR model definition file
# ParametersList = sapply(ParametersList, fixParams, simplify = FALSE) #VV: building params list
# list(sbml=sbml,id=modelid[[1]], notes=lnotes,compartments=compartments, # VV, not clear how ParametersList differs from globalParameters
# species=species,globalParameters=globalParameters, ParametersList=ParametersList, rules=rules,reactions=reactions)
list(sbml=sbml,id=modelid[[1]], notes=lnotes,compartments=compartments, # TR may revert to this??
species=species,globalParameters=globalParameters, rules=rules,reactions=reactions) # returns values accrued in parent env
}
list(.startElement = .startElement, .endElement = .endElement,
.text = .text, # , dom = function() {con}
getModel = getModel
) # function returns a list of functions, each with a common parent environment = stuff before function definitions
}
# END handler definition
# NOTE: though handlers are neat, one must question if the added baggage is worth it, i.e. compare to read.SBML of older versions
# *********************************************************************************
# The next block of three functions converts mathML XMLnode objects into R expression objects
# This approach is better than the old read.SBML approach in that the parenthesis overkill is avoided!
mathml2R <-function(node) {UseMethod("mathml2R", node)}
mathml2R.XMLDocument <-function(doc) {return(mathml2R(doc$doc$children))}
# mathml2R.default<-function(children)
# { expr <- expression() # this gets used when a "list" of children nodes are sent in
# for(i in children) expr=c(expr, mathml2R(i))
# return(expr)
# }
mathml2R.default<-function(children)
{ expr <- expression() # this gets used when a "list" of children nodes are sent in
n=length(children)
# cat("into default length n is ",n,"\n")
# for(i in children) expr=c(expr, mathml2R(i))
for(i in 1:n) expr=c(expr, mathml2R(children[[i]]))
if (n>3) {#print("n>3") # this fixes libsbml problem that times is not binary
if (expr[[1]]=="*") expr[[1]]=as.name("prod") # in R, prod takes arb # of args
if (expr[[1]]=="+") expr[[1]]=as.name("sum") # similary for sum
}
# print(children)
# print(expr)
# print("leaving default")
return(expr)
}
mathml2R.XMLNode <-function(node){
nm <- xmlName(node)
# cat("XMLNode: node name is ",nm," and the node class is",class(node),"\n")
if(nm=="power"||nm == "divide"||nm =="times"||nm=="plus"||nm=="minus") {
op <- switch(nm, power="^", divide="/",times="*",plus="+",minus="-")
val <- as.name(op)
} else if((nm == "ci")|(nm == "cn")) {
if(nm == "ci") val <- as.name(node$children[[1]]$value)
if(nm == "cn") val <- as.numeric(node$children[[1]]$value)
} else if(nm == "apply") {
val <- mathml2R(node$children)
mode(val) <- "call"
} else {cat("error: nm =",nm," not in set!\n")}
return(as.expression(val))
}
# ********** END the mathML2R block of method based on node type codes *************************
if(!require(XML)) print("Error in Read.SBML(): First Install the XML package http://www.omegahat.org/RSXML")
edoc <- xmlEventParse(filename,handlers=sbmlHandler(),ignoreBlanks = TRUE)
model=edoc$getModel() # SAX approach using the handler. Output of getModel() in edoc list is what we want.
doc <- xmlTreeParse(filename,ignoreBlanks = TRUE) # use DOM just for rules and reactions
model$htmlNotes=doc$doc$children$sbml[["model"]][["notes"]]
rules=doc$doc$children$sbml[["model"]][["listOfRules"]]
reactions=doc$doc$children$sbml[["model"]][["listOfReactions"]]
globalParameters=names(model$globalParameters)
nRules=length(rules)
#print(nRules)
#nRules=0
if (nRules>0){
# ruleIDs=NULL
for (i in 1:nRules) # for( i in 1:(nRules-1)) # VV stops 1 shy of end????
{ # assume they are assignment rules
mathml<-rules[[i]][["math"]][[1]]
model$rules[[i]]$mathmlLaw=mathml
e<-mathml2R(mathml)
model$rules[[i]]$exprLaw<-e[[1]]
model$rules[[i]]$strLaw<-gsub(" ","",toString(e[1]))
leaves<-getRuleLeaves(mathml)
r<-model$rules[[i]]$inputs<-setdiff(leaves,globalParameters) # must deduce inputs by substracting global params
# r=model$rules[[i]]$inputs
#model$rules[[i]]$idOutput=xmlAttrs(rules[[i]])[["variable"]][[1]] # moved to handler to preset the rules list length!
model$rules[[i]]$law=makeLaw(r,NULL,model$rules[[i]]$exprLaw)
# ruleIDs[i]<-model$rules[[i]]$idOutput
}
# names(model$rules)<-sapply(model$rules,function(x) x$idOutput)
}
nReactions=length(reactions)
if (nReactions>0){
# rIDs=NULL;
for (i in 1:nReactions)
{
model$reactions[[i]]$mathmlLaw=reactions[[i]][["kineticLaw"]][["math"]][[1]]
e=mathml2R(reactions[[i]][["kineticLaw"]][["math"]][[1]])
model$reactions[[i]]$exprLaw=e[[1]]
#print(e[[1]])
#print(toString(e[1]))
model$reactions[[i]]$strLaw=gsub(" ","",toString(e[1]))
#paste(as.character(model$reactions[[i]]$expr)[c(2,1,3)],collapse=""))
r=model$reactions[[i]]$reactants
# VV wants to add products in here, perhaps for reversible reactions
# r=c(model$reactions[[i]]$reactants, model$reactions[[i]]$products) #build using both reactants and products objects. TODO - add compartments
p=names(model$reactions[[i]]$parameters)
m=model$reactions[[i]]$modifiers
e=model$reactions[[i]]$exprLaw
model$reactions[[i]]$law=makeLaw(c(r,m),p,e)
# rIDs[i]<-model$reactions[[i]]$id
}
# This is for indexing by names/IDs of reactions
# names(model$reactions)<-rIDs
# names(model$reactions)<-sapply(model$reactions,function(x) x$id)
}
class(model)<-"SBMLR"
model
}
# the following is called by both readSBML and readSBMLR so it outside where both can reach it.
# Note that keeing it here instead of in a separate file => no need to document it
"makeLaw"<-function(r,p,e){
# takes reactant list r, parameter list p and rate law R expression e
# and makes a reaction rate law function out of them.
lawTempl=function(r,p=NULL){ }
i=2
for (j in seq(along=p)){
# if(!is.null(p))
# for (j in 1:length(p)){
body(lawTempl)[[i]]<-call("=",as.name(p[j]),call("[",as.name("p"),p[j]))
i=i+1}
# for (j in 1:length(r)){
for (j in seq(along=r)){
body(lawTempl)[[i]]<-call("=",as.name(r[j]),call("[",as.name("r"),r[j]))
i=i+1}
body(lawTempl)[[i]]<-e
lawTempl
}
#
#"makeLaw"<-function(inputs,pars,lawCall)
# function(inputs,lawCall,pars=NULL)
# with(as.list(c(inputs,pars)),lawCall)
#
# The next two functions are used by rules and were taken straight from read.SBML
# The idea is that SBML doesn't provide a list of atoms/leaves with rules, so we have to create them
# to place them in their model slots, and to use them to create the R function definition for the rule
# using makeLaw with a null for parameters, since they are passed global for rules.
ML2R<- function(type) # map MathML operator symbols into R symbols
switch(type,
"times" = "*",
"divide" = "/",
"plus" = "+",
"minus" = "-",
"power" = "^",
"exp" = "exp",
"ln" = "log",
"not found") # end definition of ML2R
getRuleLeaves<-function(math)
{ n=length(math)
S=c(NULL)
op=ML2R(xmlName(math[[1]]))
for (j in 2:n )
if ((xmlName(math[[j]])=="ci")|(xmlName(math[[j]])=="cn")) S=c(S,as.character(xmlValue(math[[j]]))) else
S=c(S,Recall(math[[j]]) )
S
}
Try the SBMLR package in your browser
Any scripts or data that you put into this service are public.
SBMLR documentation built on Nov. 8, 2020, 6:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions getRuleLeaves ML2R
"readSBML" <-
function(filename)
{ # takes SBML in filename.xml and maps it to a SBML class model
# using both Sax and DOM (for mathml) based parsing.
sbmlHandler <- function ()
{ # first block here sets up the parent environment used by all handler functions
sbml<-"x" # "x" is just a starting string value
modelid<-"x" #storing model id
lnotes<-NULL
compartments <- list()
reactLaws <- list()
species <- list()
rules<- list()
reactions<- list()
globalParameters=list()
reactants=NULL
products=NULL
modifiers=NULL
currRxnID=NULL
parameters=NULL # local to rate law
parameterIDs=NULL # local to rate law
globalParameterIDs=NULL
notes=FALSE; reactant=FALSE; product=FALSE
law=FALSE; parameter=FALSE; math=FALSE
# # VV's additions
# modelname<-"x" #storing the model name if separately given
# paramException <- 1
# globalParamException <- 1
# ParametersList<- list() #list of parameter objects
.startElement <- function(name, atts, ...) {
# cat("Start: Name =",name," ",paste(names(atts),atts,sep=" = "),"\n")
if(name=="sbml") sbml<<-atts
if(name=="annotation") print("skipping annotation")
if(name=="model") {modelid<<-atts[["id"]]} # VV replaces this with ...
# if(name=="model")
# { numitems <- length(atts)
# if(numitems < 1) #if model does not contain a name/id, we give it an arbitrary one.
# modelid[[1]]<<-"BioModel"
# else if(numitems == 1) # if only one attribute supplied
# {
# if(is.character(atts[[1]])) #if the attribute is a string, it must be model name. Copy to also id.
# {
# modelname<<-atts[[1]] # store model name
# modelid<<-atts[[1]] # store as model id (copy Essentially)
# }
# }
# else if(numitems > 1) #both Id and name of model are supplied, read both
# {
# modelname<<-atts[["name"]] # first element is model name
# modelid <<-atts[["id"]] # second element has to be model id
# }
# }
# if(name=="compartment") {compartments[[atts[1]]]<<-atts}
if(name=="compartment")
if( "id" %in% names(atts)) compartments[[atts["id"]]]<<-atts
# if(name=="species") {species[[atts[1]]]<<-atts }
if(name=="species")
if( "id" %in% names(atts)) species[[atts["id"]]]<<-atts
# if(name=="assignmentRule") rules[[atts[1]]]$idOutput<<-atts[[1]]
if(name=="assignmentRule") rules[[atts["variable"]]]$idOutput<<-atts[["variable"]]
# if(name=="reaction") {reactions[[atts[1]]]$id<<-atts[[1]]
# reactions[[atts[1]]]$reversible<<-as.logical(atts[[2]])
# currRxnID<<-atts[1]}
if(name=="reaction")
{
lstnames <- names(atts)
numitems <- length(lstnames)
nameslist <- list()
id <- "x"
reverse <- FALSE
name <- "x"
count <- 1
while( count <= numitems )
{
switch(lstnames[[count]],
"id" = { id = atts[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"reversible" = { reverse = as.logical(atts[[count]]) ;nameslist[[length(nameslist)+1]] <- "reversible" },
"name" = { name = as.character(atts[[count]]); nameslist[[length(nameslist)+1]] <- "name"}
)
count <- count + 1
}
reactions[[atts["id"]]]$id<<-id
reactions[[atts["id"]]]$reversible<<-reverse
# if(reverse) reactions[[atts["id"]]]$reversible<<-reverse #carry in R only if true
currRxnID<<-atts["id"]
#reactions[[atts["id"]]]$id<<-atts[["id"]]
#reactions[[atts[1]]]$reversible<<-as.logical(atts[[2]])
#currRxnID<<-atts[1]
}
if(name=="listOfReactants") reactant<<-TRUE
if(name=="listOfProducts") product<<-TRUE
if(name=="kineticLaw") law<<-TRUE
if(name=="math") math<<-TRUE
if((name=="speciesReference")&reactant)
reactants<<-c(reactants,species=atts[["species"]])
# reactants<<-c(reactants,species=atts[[1]])
if((name=="speciesReference")&product)
products<<-c(products,species=atts[["species"]])
# products<<-c(products,species=atts[[1]])
if(name=="modifierSpeciesReference")
modifiers<<-c(modifiers,species=atts[["species"]])
if((name=="parameter")&law){
parameterIDs<<-c(parameterIDs,atts[["id"]])
parameters<<-c(parameters,atts[["value"]])}
# if((name=="parameter")&law) #parameter encountered within a kinetic law definition
# {
# values <- names(atts)
# if( "id" %in% values) parameterIDs<<-c(parameterIDs,atts[["id"]])
# else {
# cat('Parameter parsed without id. Setting default id', '\n')
# parameterIDs<<-c(parameterIDs, paste("default", paramException))
# paramException <- paramException + 1
# }
# if( "value" %in% values) parameters<<-c(parameters,atts[["value"]])
# else {
# cat('Warning..Parsing parameter without value. Setting it as 0.' ,'\n')
# parameters<<-c(parameters, as.numeric(0))
# }
# }
if((name=="parameter")&(!law)){
globalParameterIDs<<-c(globalParameterIDs,atts[["id"]])
globalParameters<<-c(globalParameters,as.numeric(atts[["value"]]))}
# if((name=="parameter")&!law) #parameter encountered outside a kinetic law definition - So in globalparamslist
# {
# #cat("within parameters:", atts[["id"]], atts[["value"]], "\n")
# values <- names(atts)
# if( "id" %in% values) {
# globalParameterIDs<<-c(globalParameterIDs,atts[["id"]])
# ParametersList[[atts["id"]]] <<- atts #our new list of Parameter Objects
# }
# else {
# cat('Global Parameter parsed without id. Setting default id', '\n')
# tempParamId <- paste("Globaldefault", globalParamException)
# globalParameterIDs<<-c(globalParameterIDs, tempParamId)
# ParametersList[[tempParamId]] <<- atts
# globalParamException <- globalParamException + 1
# }
# if( "value" %in% values) globalParameters<<-c(globalParameters,as.numeric(atts[["value"]]))
# else {
# cat('Warning..Parsing Global parameter without value. Setting it to 0.', '\n')
# globalParameters<<-c(globalParameters, as.numeric(0))
# }
# } # end if param in law
} # end .startElement()
.endElement <- function(name) {
if(name=="listOfReactants") reactant<<-FALSE
if(name=="listOfProducts") product<<-FALSE
if(name=="kineticLaw") law<<-FALSE
if(name=="math") math<<-FALSE
if((name=="listOfParameters")&(!law)) names(globalParameters)<<-globalParameterIDs
if(name=="reaction") {
names(reactants)<<-NULL
names(modifiers)<<-NULL
names(products)<<-NULL
reactions[[currRxnID]]$reactants<<-reactants
reactions[[currRxnID]]$modifiers<<-modifiers
reactions[[currRxnID]]$products<<-products
parameters<<-as.numeric(parameters)
names(parameters)<<-parameterIDs
reactions[[currRxnID]]$parameters<<-parameters
reactants<<-NULL;products<<-NULL
modifiers<<-NULL;parameters<<-NULL;parameterIDs<<-NULL }
}
.text <- function(x, ...) {
if (!math) lnotes<<-c(lnotes,x)
# cat("Txt:", x,"\n")
}
getModel<- function()
{
# fixComps=function(x) {
# lst=list(x[[1]],as.numeric(x[[2]]) );
# names(lst)<-names(x);
# lst
# }
# VV replaces fixComps with the following:
fixComps=function(x)
{
lstnames <- names(x)
count <- 1
numit <- length(lstnames)
id <- "x"
size <- 0
name <- "x"
nameslist <- list()
while( count <= numit )
{
switch(lstnames[[count]],
"id" = { id = x[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"size" = { size = as.numeric(x[[count]]) ;nameslist[[length(nameslist)+1]] <- "size" },
"name" = { name = as.character(x[[count]]); nameslist[[length(nameslist)+1]] <- "name"}
)
count = count + 1
}
#cat("Compartment id:", id,"\n", "Compartment size:" , size, "\n","Compartment name:", name, "\n")
if(numit == 2) # only 2 attributes present. We need to find them.
{ if(id == "x") #id not set but name and size are.
id <- "default"
else if(name == "x") #name not set, we copy the id.
name <- id
else if(size== "0") #size not set
size <- 1 #arbitrary setting as 1
lst = list(id,size,name)
names(lst)<-c("id","size","name")
lst
} else if(numit == 3) # 3 attributes/items present.
{ lst = list(id,size,name)
names(lst)<-c("id","size","name")
lst
}
}
# fixSpecies=function(x) {
# lst=list(x[[1]],as.numeric(x[[2]]),x[[3]],as.logical(x[[4]]));
# names(lst)<-c("id","ic","compartment","bc");
# lst
# }
# VV replaces fixSpecies with the following
fixSpecies=function(x)
{
#cat (names(x), "\n")
#cat(toString(x) , "\n")
numitems <- length(x)
lstnames <- names(x)
count <-1
id <- "x" #species Id
ic <- 0 #species initial concentration
compart <- "def" #species compartment
bc <- FALSE #species boundary condition
name <- "def"
nameslist <- list()
while( count <= numitems)
{
switch(lstnames[[count]],
"id" = { id <- x[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"name" = { name <- x[[count]]; nameslist[[length(nameslist)+1]] <- "name"},
"initialConcentration" = { ic <- as.numeric(x[[count]]) ;nameslist[[length(nameslist)+1]] <- "ic" },
"compartment" = { compart <- as.character(x[[count]]); nameslist[[length(nameslist)+1]] <- "compartment"},
"boundaryCondition" = { bc <- as.logical(x[[count]]); nameslist[[length(nameslist)+1]] <- "bc"}
)
count = count + 1
}
#lst = list(id,ic,compart,bc, name)
lst = list(id,as.numeric(ic), compart, as.logical(bc))
names(lst) <- c("id","ic","compartment","bc")
#names(lst)<-c("id","ic","compartment","bc", "name");
lst
}
# and VV adds in fixParams
fixParams=function(x)
{
numitems <- length(x)
lstnames <- names(x)
count <-1
id <- "x" #Parameter Id
value <- 0 #Parameter value
name <- "def"
constant <- FALSE
nameslist <- list()
while( count <= numitems)
{
switch(lstnames[[count]],
"id" = { id <- x[[count]]; nameslist[[length(nameslist)+1]] <- "id"},
"name" = { name <- x[[count]]; nameslist[[length(nameslist)+1]] <- "name"},
"value" = { value <- as.numeric(x[[count]]) ;nameslist[[length(nameslist)+1]] <- "value" },
"constant" = { constant <- as.logical(x[[count]]) ; nameslist[[length(nameslist)+1]] <- "constant"}
)
count = count + 1
}
lst = list(id,as.numeric(value))
names(lst) <- c("id","value")
lst
}
compartments=sapply(compartments,fixComps, simplify = FALSE)
#species=t(sapply(species,fixSpecies, simplify = TRUE)[2:4,]) # this changes the species model structure for better looks in R dumps
species=sapply(species,fixSpecies, simplify = FALSE) # this keeps the better looks in the SBMLR model definition file
# ParametersList = sapply(ParametersList, fixParams, simplify = FALSE) #VV: building params list
# list(sbml=sbml,id=modelid[[1]], notes=lnotes,compartments=compartments, # VV, not clear how ParametersList differs from globalParameters
# species=species,globalParameters=globalParameters, ParametersList=ParametersList, rules=rules,reactions=reactions)
list(sbml=sbml,id=modelid[[1]], notes=lnotes,compartments=compartments, # TR may revert to this??
species=species,globalParameters=globalParameters, rules=rules,reactions=reactions) # returns values accrued in parent env
}
list(.startElement = .startElement, .endElement = .endElement,
.text = .text, # , dom = function() {con}
getModel = getModel
) # function returns a list of functions, each with a common parent environment = stuff before function definitions
}
# END handler definition
# NOTE: though handlers are neat, one must question if the added baggage is worth it, i.e. compare to read.SBML of older versions
# *********************************************************************************
# The next block of three functions converts mathML XMLnode objects into R expression objects
# This approach is better than the old read.SBML approach in that the parenthesis overkill is avoided!
mathml2R <-function(node) {UseMethod("mathml2R", node)}
mathml2R.XMLDocument <-function(doc) {return(mathml2R(doc$doc$children))}
# mathml2R.default<-function(children)
# { expr <- expression() # this gets used when a "list" of children nodes are sent in
# for(i in children) expr=c(expr, mathml2R(i))
# return(expr)
# }
mathml2R.default<-function(children)
{ expr <- expression() # this gets used when a "list" of children nodes are sent in
n=length(children)
# cat("into default length n is ",n,"\n")
# for(i in children) expr=c(expr, mathml2R(i))
for(i in 1:n) expr=c(expr, mathml2R(children[[i]]))
if (n>3) {#print("n>3") # this fixes libsbml problem that times is not binary
if (expr[[1]]=="*") expr[[1]]=as.name("prod") # in R, prod takes arb # of args
if (expr[[1]]=="+") expr[[1]]=as.name("sum") # similary for sum
}
# print(children)
# print(expr)
# print("leaving default")
return(expr)
}
mathml2R.XMLNode <-function(node){
nm <- xmlName(node)
# cat("XMLNode: node name is ",nm," and the node class is",class(node),"\n")
if(nm=="power"||nm == "divide"||nm =="times"||nm=="plus"||nm=="minus") {
op <- switch(nm, power="^", divide="/",times="*",plus="+",minus="-")
val <- as.name(op)
} else if((nm == "ci")|(nm == "cn")) {
if(nm == "ci") val <- as.name(node$children[[1]]$value)
if(nm == "cn") val <- as.numeric(node$children[[1]]$value)
} else if(nm == "apply") {
val <- mathml2R(node$children)
mode(val) <- "call"
} else {cat("error: nm =",nm," not in set!\n")}
return(as.expression(val))
}
# ********** END the mathML2R block of method based on node type codes *************************
if(!require(XML)) print("Error in Read.SBML(): First Install the XML package http://www.omegahat.org/RSXML")
edoc <- xmlEventParse(filename,handlers=sbmlHandler(),ignoreBlanks = TRUE)
model=edoc$getModel() # SAX approach using the handler. Output of getModel() in edoc list is what we want.
doc <- xmlTreeParse(filename,ignoreBlanks = TRUE) # use DOM just for rules and reactions
model$htmlNotes=doc$doc$children$sbml[["model"]][["notes"]]
rules=doc$doc$children$sbml[["model"]][["listOfRules"]]
reactions=doc$doc$children$sbml[["model"]][["listOfReactions"]]
globalParameters=names(model$globalParameters)
nRules=length(rules)
#print(nRules)
#nRules=0
if (nRules>0){
# ruleIDs=NULL
for (i in 1:nRules) # for( i in 1:(nRules-1)) # VV stops 1 shy of end????
{ # assume they are assignment rules
mathml<-rules[[i]][["math"]][[1]]
model$rules[[i]]$mathmlLaw=mathml
e<-mathml2R(mathml)
model$rules[[i]]$exprLaw<-e[[1]]
model$rules[[i]]$strLaw<-gsub(" ","",toString(e[1]))
leaves<-getRuleLeaves(mathml)
r<-model$rules[[i]]$inputs<-setdiff(leaves,globalParameters) # must deduce inputs by substracting global params
# r=model$rules[[i]]$inputs
#model$rules[[i]]$idOutput=xmlAttrs(rules[[i]])[["variable"]][[1]] # moved to handler to preset the rules list length!
model$rules[[i]]$law=makeLaw(r,NULL,model$rules[[i]]$exprLaw)
# ruleIDs[i]<-model$rules[[i]]$idOutput
}
# names(model$rules)<-sapply(model$rules,function(x) x$idOutput)
}
nReactions=length(reactions)
if (nReactions>0){
# rIDs=NULL;
for (i in 1:nReactions)
{
model$reactions[[i]]$mathmlLaw=reactions[[i]][["kineticLaw"]][["math"]][[1]]
e=mathml2R(reactions[[i]][["kineticLaw"]][["math"]][[1]])
model$reactions[[i]]$exprLaw=e[[1]]
#print(e[[1]])
#print(toString(e[1]))
model$reactions[[i]]$strLaw=gsub(" ","",toString(e[1]))
#paste(as.character(model$reactions[[i]]$expr)[c(2,1,3)],collapse=""))
r=model$reactions[[i]]$reactants
# VV wants to add products in here, perhaps for reversible reactions
# r=c(model$reactions[[i]]$reactants, model$reactions[[i]]$products) #build using both reactants and products objects. TODO - add compartments
p=names(model$reactions[[i]]$parameters)
m=model$reactions[[i]]$modifiers
e=model$reactions[[i]]$exprLaw
model$reactions[[i]]$law=makeLaw(c(r,m),p,e)
# rIDs[i]<-model$reactions[[i]]$id
}
# This is for indexing by names/IDs of reactions
# names(model$reactions)<-rIDs
# names(model$reactions)<-sapply(model$reactions,function(x) x$id)
}
class(model)<-"SBMLR"
model
}
# the following is called by both readSBML and readSBMLR so it outside where both can reach it.
# Note that keeing it here instead of in a separate file => no need to document it
"makeLaw"<-function(r,p,e){
# takes reactant list r, parameter list p and rate law R expression e
# and makes a reaction rate law function out of them.
lawTempl=function(r,p=NULL){ }
i=2
for (j in seq(along=p)){
# if(!is.null(p))
# for (j in 1:length(p)){
body(lawTempl)[[i]]<-call("=",as.name(p[j]),call("[",as.name("p"),p[j]))
i=i+1}
# for (j in 1:length(r)){
for (j in seq(along=r)){
body(lawTempl)[[i]]<-call("=",as.name(r[j]),call("[",as.name("r"),r[j]))
i=i+1}
body(lawTempl)[[i]]<-e
lawTempl
}
#
#"makeLaw"<-function(inputs,pars,lawCall)
# function(inputs,lawCall,pars=NULL)
# with(as.list(c(inputs,pars)),lawCall)
#
# The next two functions are used by rules and were taken straight from read.SBML
# The idea is that SBML doesn't provide a list of atoms/leaves with rules, so we have to create them
# to place them in their model slots, and to use them to create the R function definition for the rule
# using makeLaw with a null for parameters, since they are passed global for rules.
ML2R<- function(type) # map MathML operator symbols into R symbols
switch(type,
"times" = "*",
"divide" = "/",
"plus" = "+",
"minus" = "-",
"power" = "^",
"exp" = "exp",
"ln" = "log",
"not found") # end definition of ML2R
getRuleLeaves<-function(math)
{ n=length(math)
S=c(NULL)
op=ML2R(xmlName(math[[1]]))
for (j in 2:n )
if ((xmlName(math[[j]])=="ci")|(xmlName(math[[j]])=="cn")) S=c(S,as.character(xmlValue(math[[j]]))) else
S=c(S,Recall(math[[j]]) )
S
}
Try the SBMLR package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.