inst/extdata/schoeberl-processing.R

In robertness/pathflow: Converts biochemical reaction networks to pathway maps.

#' A data processing function for the raw Schoeberl data.
#' Used to convert a reaction into a set of edges in a bipartite reaction network.
#' 
#' @param reaction a character in the form "A + B  -> A-B" 
#' or "A + B  <- A-B" (in the reversable case).
#' @return An reaction net (petri net)-style bipartite edge list 
parseBipartiteEdges <- function(reaction){
  direction <- ifelse(grepl("->", reaction), "->", 
                      ifelse(grepl("<-", reaction), "<-", "error"))
  if(direction == "error") stop("There was an error in parsing reaction direction.")
  product.reactant.split <- as.list(unlist(strsplit(reaction, direction)))
  product.reactant.split <- as.list(gsub(" ", "", product.reactant.split))
  product.reactant.split <- lapply(product.reactant.split, function(item){
    unlist(strsplit(item, "\\+"))}) 
  if(direction == "->"){
    direct.edges <- expand.grid(product.reactant.split, stringsAsFactors = F)
  } else {
    direct.edges <- expand.grid(product.reactant.split[c(2,1)], stringsAsFactors = F)
  }
  names(direct.edges) <- NULL
  from.edges <- cbind(direct.edges[, 1], reaction)
  to.edges <- cbind(reaction, direct.edges[, 2])
  colnames(to.edges) <- colnames(from.edges) <- NULL
  bipartite.edges <- unique(as.data.frame(rbind(from.edges, to.edges))) 
  bipartite.edges <- tryCatch('names<-'(bipartite.edges, c("from", "to")),
                              error = function(e) NULL)
  if(is.null(bipartite.edges)) stop("Reaction ", reaction, " was not processed quickly." )
  bipartite.edges$to <- as.character(bipartite.edges$to)
  bipartite.edges$from <- as.character(bipartite.edges$from)
  bipartite.edges
}

#' Given a choice of two vertices to remove, remove the vertex with the shortest
#' path to the output vertex.
#' @param g A graph object representing a petri-net style reaction network  
#' @param v1 a vertex name.
#' @param v2 a vertex name.
#' @param target a vertex name.
#' @return g modified such that either v1, v2 is removed, or neither is removed. 
deleteOneOfTwoVertices <- function(g, v1, v2){
  g.v1 <- delete.vertices(g, v2)
  pl.v1 <- as.numeric(shortest.paths(g.mod1, v1, target, mode = "out"))
  g.v2 <- delete.vertices(g, v1)
  pl.v2 <- as.numeric(shortest.paths(g.mod2, v2, target, mode = "out"))
  if(is.infinite(pl.v1) && is.infinite(pl.v2)){
    warning("Neither ", v1, " or ", v2, " have any path to the target.")
    return(g)
  }
  if(pl.v1 >= pl.v2) return(g.v1) else return(g.v2)
}

#' Eliminating reversable reaction
#' 
#' Reversable transactions cancel one another out in the incidence matrix.  
#' This compares both reactions (forward and backward) associated with the reversable
#' transaction, and removes the one that has the longerst path to the output vertex.
#' If both paths are equal in length, the forward transaction is selected.
#' 
#' @param g A graph object representing a petri-net style reaction network, with vertex
#' attributes: "name", "type" ('reaction' or 'product'), and "reaction_number" (a character, not a numeric).  
#' @param v a vertex name.
#' @return g modified such that either v1, v2 is removed, or neither is removed. 
makeOneWay <- function(g, v){
  if(V(g)[v]$type == "species") return(g)
  # if(!isBDownstreamOfA(g, v, V(g)[is.input])) return(g) # Only fix reactions on path from input 
  same_reaction <- V(g)[type == "reaction"][reaction_number == V(g)[v]$reaction_number]
  if(length(same_reaction) == 1) return(g)
  is_forward <- grepl("->", same_reaction$name)
  forward <- same_reaction[is_forward]$name
  backward <- same_reaction[!is_forward]$name
  g <- deleteOneOfTwoVertices(g, forward, backward, V(g)[is.output])
  g
}

#' When traversing a graph to remove reversable transactions, provide the next set of 
#' reactions vertices given the current vertex.
#' @param g A graph object representing a petri-net style reaction network, with vertex
#' attributes: "name", "type" ('reaction' or 'product'), and "reaction_number" (a character, not a numeric).  
#' @param v a vertex name.
#' @return A character array of vertex names.
getReactionAncestors <- function(g, v){
  ancestor.reactions <- iparents(g, iparents(g, v))
  if(any(V(g)[ancestor.reactions]$type != "reaction")){
    stop("Got a species when expecting a reaction.")
  }
  check.downstream <- unlist(lapply(ancestor.reactions, function(reaction){
    isBDownstreamOfA(g, reaction, input)
  }))
  next.reactions <- ancestor.reactions[check.downstream]
  next.reactions
}

getDownstreamReactions <- function(g, v){
 V(g)[type == "reaction"][getDownstreamNodes(g, v)] 
}