R/PathFlow.R

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

# Returns a list where each element is an invariant containing an array of the members of the invariants

#' Calculates minimimal invariants from a incidence matrix.
#'
#' To get place invariants, use incidence matrix (rows corresponding to places and columns to transitions)
#' To get transition invariants, use transpose of incidence matrix. 
#' 
#' @param S A matrix relating relationships between places and transitions.
#'
#' @return A named list of invariants
getMinimalInvariants <- function(S){
  nullMat <- round(MASS::Null(S), 4)
  if(length(nullMat) == 0) stop("Basis is empty.")
  rownames(nullMat) <- rownames(S)
  invariants <- list()
  for(j in 1:ncol(nullMat)){
    nullVec <- nullMat[, j]
    valCounts <- table(nullVec)
    valCounts <- valCounts[valCounts != 1]
    valCounts <- valCounts[names(valCounts) != 0] #clear vals that rounded to 0
    valCounts <- valCounts[names(valCounts) > 0] #clear negative elements
    if(length(valCounts) > 0){
      for(i in 1:length(valCounts)){
        val <- valCounts[i]
        invariant <- names(nullVec[nullVec == names(val)])
        invariants <- c(invariants, list(invariant))
      }  
    } 
  }
  lengths <- sapply(invariants, length)
  invariants <- invariants[order(lengths)]
  if(length(invariants) > 1){
    for(i in 1:(length(invariants) - 1)){
      for(j in (i+1):length(invariants)){
        invariants[[j]] <- setdiff(invariants[[j]], invariants[[i]])
      }
    }
  }
  invariants <- invariants[which(sapply(invariants, length) > 0)]
  invariants
}

#' Converts a bipartite biochemical network graph into a matrix reflecting pre transition edges,
#' a matrix reflecting post transition edges, and an incidence matrix, which equals post - pre.
#' 
#' @param g An igraph object. This must have a vertex attribute called "type" that labels each 
#' vertex as either a "place" or "transition"
#'
#' @return A named list containing the "pre", "post", and "S" matrices
getPetriNetMatrices <- function(g){
  places <- V(g)[V(g)$type == "place"]$name
  transitions <- V(g)[V(g)$type == "transition"]$name
  adjMat <- as.matrix(get.adjacency(g))
  pre <- adjMat[places, transitions]
  post <- t(adjMat[transitions, places])
  S <- post - pre
  list(pre = pre, post = post, S = S)
}

#' Find overlapping p-invariants
#' 
#' One of the sources of edges in the pathway map that is implied by a 
#' biochemical network is cases of p-invariants that share common vertices.
#' These 'overlapping' invariants will share a common edge in the pathway map.
#' @param pInvariants a list where each element contains a character array
#' of the names of vertices in a network.  These are returned by the
#'  \code{getMinimalInvariants} function.
#'
#' @return A list, where each element is a pair of indices of connected
#' pInvariants in the final matrix.
getOverlapEdges <- function(pInvariants){
  # Calculate all pairs of pInvariants. 
  # Convert combination matrix to data.frame 
  # so can treat as list in lapply ---------------------------
  pairwise.combinations <- as.data.frame(combn(1:length(pInvariants), 2))
  overlap.edges <- lapply(pairwise.combinations, function(pair){
    intersection <- intersect(pInvariants[[pair[1]]], pInvariants[[pair[2]]]) 
    check.intersection <- length(intersection) > 0
    if(check.intersection) return(pair)
  })
  overlap.edges <- overlap.edges[!unlist(lapply(overlap.edges, is.null))]
  names(overlap.edges) <- NULL
  overlap.edges
}

checkTransitionConnectsToPlace <- function(t, p, g){
  if(!(V(g)[p]$type == "place" && V(g)[t]$type == "transition")){
    stop("'p' wasn't a place or 't' wasnt a transition.")
  }
  connected.places <- V(g)[nei(V(g)[t], mode="out")]
  p %in% connected.places
}

checkFromP2P <- function(g, p1, p2, tInvariants){
  v1.out.transitions <- V(g)[nei(V(g)[v1], mode="out")]
  if(any(v1.out.transitions$type != "transition")) {
    stop("Expected transitions but got places.")
  }
  t.invariant.intersection <- lapply(tInvariants, function(tInvariant){
    intersect(v1.out.transitions$name, tInvariant)
  })
  is.empty <- !any(unlist(lapply(t.invariant.intersection, 
                                   function(item) length(item) > 0)))
  if(is.empty) return(FALSE)
  # Now look at each group of transitions ---------------------------
  check.connected.to.p2 <- lapply(t.invariant.intersection, function(intersection){
    # Check if any of those transitions are connected to p2---------------------------
    any(unlist(lapply(intersection, checkTransitionConnectsToPlace, p = p2, g = g)))
  })
  any(unlist(check.connected.to.p2))
}

check2PlacesConnected <- function(g, p1, p2, tInvariants){
  p1.to.p2 <- checkFromP2P(g, p1, p2, tInvariants)
  p2.to.p1 <- checkFromP2P(g, p2, p1, tInvariants)
  check.connected <- p1.to.p2  || p2.to.p1
  if(check.connected) message(p1, " and ", p2, " are connected!")
  check.connected
}

check2PInvariantsConnection <- function(pInvariant1, pInvariant2, tInvariants, g){
  place.pairs <- expand.grid(pInvariant1, pInvariant2, stringsAsFactors = F)
  place.pairs <- as.data.frame(t(place.pairs), stringsAsFactors = F)
  check.connected.pair <- any(unlist(lapply(place.pairs, function(place.pair){
    check2PlacesConnected(g, place.pair[1], place.pair[2], tInvariants)
  })))
  if(check.connected.pair) message("Found a connection between pInvariants")
  check.connected.pair
}

getConnectedEdges <- function(pInvariants, tInvariants, g){
  # Calculate all pairs of pInvariants. 
  # Convert combination matrix to data.frame 
  # so can treat as list in lapply ---------------------------
  pairwise.combinations <- as.data.frame(combn(1:length(pInvariants), 2))
  output <- lapply(pairwise.combinations, function(pair){
    pInvariant1 <- pInvariants[[pair[1]]]
    pInvariant2 <- pInvariants[[pair[2]]]
    is.connected <- check2PInvariantsConnection(pInvariant1, pInvariant2, tInvariants, g)
    if(is.connected) return(pair) else return(NULL)
  })
  output
}