Cacher is the code snippet organizer for pro developers

# The following simulates data from a multilayer perceptron model of a protein kinase cascade, 
# a key signal processing mechanism in cells.
# Use the generatePerceptronGraph function to visualize the structure of the cascade.
# Each node represents a protein that can be "activated" or "unactivated".
# The value for each node is the proportion of activated over total amount of the protein.
# The rate of change of that proportion is regulated by an ODE, 
# where the amount increases depending on the proportion of activated nodes in the previous layer,
# and decreases depending on a rate of inactivation.
# This models phosphorylation (activation) and dephosphorylation (inactivation) via a phosphotase 
# with simple first order kinetics.
# The simulation can produce and visualize the dynamics of a system given a set of parameters and starting values,
# and can produce replicates of "snapshots" -- single point observations from the steady state trajectory of 
# each node, which can be used with multivariate statistical models designed to analyze real life 
# proteomics data.
# There is an intervention option, which prevents activation of a given node.  This is intended to 
# mimic an experimental perturbation, which would be used to infer causality in the system.
# Technical notes: degredation (decrease in total amount of protein over time) is not modeled as it is assumed
# degredation does not occur on the time scale of cell signaling.  The amount of phosphotase for a given protein
# is assumed constant and not explicitly modeled.

generatePerceptronGraph <- function(numLayers, layerWidth){
  # Generates a graphNel structure representing a multi-layer perceptron
  # Starts and ends with single nodes "R" and "T", with multiple "layers" in between.
  # numLayers is the number of Layers, layerWidth is the number of nodes in the layer
  require(graph)
  layerNames <- LETTERS[1:numLayers]
  layer <- layerNames[1]
  layerList <- lapply(layerNames, function(layer) {
    paste(rep(layer, each = 2), 1:layerWidth, sep = "")
  })
  edgeMat <- NULL
  for(i in 2:length(layerList)){
    for(j in 1:length(layerList[[i]])){
      for(k in 1:length(layerList[[i - 1]])){
        edgeMat <- rbind(edgeMat, c(layerList[[i -1 ]][j], layerList[[i]][k]))
      }
    }
  }
  top <- c("R1", "R2")
  edgeMat <- rbind(do.call("rbind", lapply(top, 
                                           function(nd) t(rbind(nd, layerList[[1]])))), 
                   edgeMat)
  edgeMat <- rbind(edgeMat, 
                   t(sapply(layerList[[numLayers]], function(nd) c(nd, "T"))))
  dimnames(edgeMat) <- NULL
  ftM2graphNEL(edgeMat)
}

generatePerceptronEdgeMatrix <- function(numLayers, layerWidth){
  # Generates an edge matrix structure representing a multi-layer perceptron
  # Starts and ends with single nodes "R" and "T", with multiple "layers" in between.
  # numLayers is the number of Layers, layerWidth is the number of nodes in the layer
  require(graph)
  layerNames <- LETTERS[1:numLayers]
  layer <- layerNames[1]
  layerList <- lapply(layerNames, function(layer) {
    paste(rep(layer, each = 2), 1:layerWidth, sep = "")
  })
  top <- c("R1", "R2")
  edgeMat <- NULL
  for(i in 2:length(layerList)){
    for(j in 1:length(layerList[[i]])){
      for(k in 1:length(layerList[[i - 1]])){
        edgeMat <- rbind(edgeMat, c(layerList[[i -1 ]][j], layerList[[i]][k]))
      }
    }
  }
  edgeMat <- rbind(do.call("rbind", lapply(top, 
                                           function(nd) t(rbind(nd, layerList[[1]])))), 
                   edgeMat)
  edgeMat <- rbind(edgeMat, 
                   t(sapply(layerList[[numLayers]], function(nd) c(nd, "T"))))
  dimnames(edgeMat) <- NULL
  edgeMat
}

iParents <- function(g, v){
  #A simple "parents" function that recovers the parents of a vertex 
  #in a igraph class graph
  #g is a graph of class igraph
  #v is  a igraph.vs vertex
  require(igraph)
  V(g)[nei(v, mode = "in")]
}

getRate <- function(v.name, g){
  v <- V(g)[v.name]
  v.parents <- iParents(g, v)
  rateLaw <- sum(E(g)[v.parents %->% v]$phosParams * v.parents$on * v$off, 
                 -1 * v$taseParams * v$on)
  rateLaw
}

getRateLaw <- function(t, y, parms){
  V(g)$on <- y[1:vcount(g)] #Populate initial values values
  V(g)$off <- y[(vcount(g) + 1):(2 * vcount(g))]
  E(g)$phosParams <- parms[phosParams] #Add the parameters directly to graph object  
  V(g)$taseParams <- parms[taseParams]
  rl <- sapply(V(g)$name, getRate, g = g) #Pass g directly to rate function along with params.
  c(rl, -1 *rl)
}

getODEFunction <- function(g, y, parms){
  function(t, y, parms) list(getRateLaw(t, y, parms))
}

simTrajectory <- function(fList, y, parms, times){
  require(deSolve)
  out <- ode(y, times, func=fList, parms=parms)
  out2 <- out[, 1:(vcount(g) + 1)]
  class(out2) <- class(out)
  return(out2)
}

testParams <- function(g, parms, stopTime, lineColor, int=NA){
  y.on <- runif(vcount(g))
  if(!is.na(int)){#Simulate the intervention by adjusting starting values and rates
    intIdx <- V(g)$name == int
    y.on[intIdx] <- 0
    nullIdx <- grepl(paste("->", int, sep=""), names(parms))
    parms[grepl(paste("->", int, sep=""), names(parms))] <- 0
  }
  y.off <- 1 - y.on
  names(y.on) <- paste(V(g)$name, ".on", sep="")
  names(y.off) <- paste(V(g)$name, ".off", sep="")
  y <- c(y.on, y.off)
  odeFunctionList <- getODEFunction(g, y, parms)
  trajectory <- simTrajectory(odeFunctionList, y, parms, seq(0, stopTime, .1))
  colnames(trajectory)[ 2:ncol(trajectory)] <- V(g)$name
  plot(trajectory, ylim=c(0,1), col = lineColor)
}

simCells <- function(g, parms, snapshot, m, int=NA){
  # The int argument can take the name of a node.
  # That node will be targeted for intervention.
  # The starting value for that node is set to 0,
  # as is its rate of activation.
  dataMat <- matrix(0, nrow = m, ncol = vcount(g) , m)
  for(i in 1:m){
    msg = paste(" replicate:", i)
    prefix = "Observational data"
    y.on <- runif(vcount(g))
    if(!is.na(int)){#Simulate the intervention by adjusting starting values and rates
      intIdx <- V(g)$name == int
      y.on[intIdx] <- 0
      nullIdx <- grepl(paste("->", int, sep=""), names(parms))
      parms[grepl(paste("->", int, sep=""), names(parms))] <- 0
      prefix <- paste("Intervention in", int, sep="")
    }
    message(paste(prefix, msg))
    y.off <- 1 - y.on
    names(y.on) <- paste(V(g)$name, ".on", sep="")
    names(y.off) <- paste(V(g)$name, ".off", sep="")
    y <- c(y.on, y.off)
    odeFunctionList <- getODEFunction(g, y, parms)
    trajectory <- simTrajectory(odeFunctionList, y, parms, seq(snapshot - .1, snapshot, .01))
    cellSim <- trajectory[nrow(trajectory), -1]
    dataMat[i, ] <- cellSim
  }
  colnames(dataMat) <- names(cellSim)
  return(dataMat)
}

getRateLaw <- function(t, y, parms){
  V(g)$on <- y[1:vcount(g)] #Populate initial values values
  V(g)$off <- y[(vcount(g) + 1):(2 * vcount(g))]
  E(g)$phosParams <- parms[phosParams] #Add the parameters directly to graph object  
  V(g)$taseParams <- parms[taseParams]
  rl <- sapply(V(g)$name, getRate, g = g) #Pass g directly to rate function along with params.
  c(rl, -1 *rl)
}

# # Simulate Data from a graph with 3 middle layers (plus an input layer) 
# ## Visualize graph structure 
# g <- generatePerceptronGraph(3, 2)
# plot(g)
# ## Recreate the graph using the igraph class
# require(igraph)
# em <- generatePerceptronEdgeMatrix(numLayers = 3, 2)
# g <- graph.data.frame(as.data.frame(em), directed = T)
# ## Set the parameters
# phosParams <- apply(em, 1, function(row) paste(row, collapse="->"))
# taseParams <- paste(V(g)$name, ".alpha", sep="")
# paramNames <- c(phosParams, taseParams)
# parms <- c(rep(.21, ecount(g)), 0, 0, rep(.3, vcount(g) - 2))
# names(parms) <- paramNames
# ## Visualize system
# testParams(g, parms)
# ## Select a "snapshot" -- a timepoint to grab from the system
# snapshot <- 30 
# ## Select number of replicates per sample
# m <- 5000
# ## Simulate observational data
# obsdata <- simCells(g, parms, 30, m)
# ## Simulate interventional data
# intDataList <- lapply(V(g)$name, function(int) simCells(g, parms, 15, m, int = int))