R/create-network.R

In davidkwca/inferTPBN: Infer Threshold PBNs from Timeseries Data.

##'  Create random threshold PBN.
##'
##'
##' This first uses its parameters to let BoolNet create a network with the
##' desired nodes, number of edges, topology, etc.
##' It then creates adds all the relevant parameters to make it usable as a
##' TPBN.
##'
##' @param inputProbabilities List of probabilities determining how many
##'   regulatory sets there are per gene.
##' @param ... Other parameters, used by BoolNet.
##' @return A network.
##' @seealso \code{\link{transformNetwork}} \code{\link{setupPBN}}
##' @examples
##' net <- createNetwork(inputProbabilities=1, 5, 2, topology="homogeneous")
createNetwork <- function(inputProbabilities=1, ...){
  net <- generateRandomNKNetwork(...)
  n <- length(net$genes)
  for (g in 1:n){
    ## To create a PBN, we need a list.
    net$interactions[[g]] <- list(net$interactions[[g]])

    ## We'll take the number of regulating sets according to inputProbabilities.
    numberInputs <-
      sample(1:length(inputProbabilities), size=1, prob=inputProbabilities)
    c <- runif(numberInputs)
    c <- c / sum(c)
    ## Instantiate each possible regulatory set randomly.
    for (i in 1:numberInputs){
      k <- length(net$interactions[[g]][[1]])
      net$interactions[[g]][[i]] <- net$interactions[[g]][[1]]
      net$interactions[[g]][[i]]$input <- sample(1:n, size=k)
      net$interactions[[g]][[i]]$probability <- c[i]
    }
  }
  return(net)
}
##' Add TPBN parameters to a BoolNet network.
##'
##'
##' Add perturbation probability, threshold parameters, and Boolean functions
##' to the network, to turn it into a TPBN.
##'
##' @param net A network generated by \code{\link{createNetwork}}.
##' @param p Noise strength.
##' @return A threshold PBN.
##' @seealso \code{\link{createNetwork}} \code{\link{A}} \code{\link{aToFunction}}
##' @examples
##' # Create a BoolNet network and transform it.
##' net <- generateRandomNKNetwork(5, 2, "homogeneous")
##' net <- transformNetwork(net)
transformNetwork <- function(net, p=0.01){
  ## Translates a BoolNet Network into a thresholded PBN
  net$flipProbability <- p
  net <- A(net)
  net <- aToFunction(net)
  return(net)
}

##' Add threshold parameters to network.
##'
##'
##' Threshold parameters determine whether an input will have a positive or
##' negative relationship with the target genes.
##'
##' @param net The network to be transformed.
##' @return The network with associated parameters.
A <- function(net){
  ## This function assigns every regulatory set of every gene a vector a,
  ## which will later determine the linear influence of the regulator.
  for (gene in net$genes){
    for (i in 1:length(net$interactions[gene][[1]])){
      input <- net$interactions[gene][[1]][[i]]$input
      net$interactions[gene][[1]][[i]]$a <-
        sample(c(-1, 1), length(input), replace=TRUE)
    }
  }
  return(net)
}
##' Turn threshold parameters into Boolean functions.
##'
##'
##' Take a network with a set of threshold parameters and compute the correct
##' Boolean functions belonging to these parameters.
##'
##' @param net A network which has threshold parameters.
##' @return A Network with logic functions associated.
##' @examples
##' net <- generateRandomNKNetwork(5, 2, "homogeneous")
##' net <- A(net)
##' net <- aToFunction(net)
aToFunction <- function(net){
  ## This translates all $a vectors into regulatory functions.
  for (gene in net$genes){
    for (i in 1:length(net$interactions[gene][[1]])){
      f <- aToFunctionGene(net, gene, i)
      net$interactions[gene][[1]][[i]]$func <- f
    }
  }
  return(net)
}

aToFunctionGene <- function(net, gene, i){
  ## This translates a specific gene regulatory function's $a into the function.
  a <- net$interactions[gene][[1]][[i]]$a
  return(a2function(a))
}

a2function <- function(a){
  ## This translates an arbitrary vector a into the output list for the
  ## associated function.
  k <- length(a)
  X <- list()

  ## Create the truth table on k inputs.
  for (j in 1:k){
    X[[j]] <- 0:1
  }
  X <- expand.grid(X)

  ## Compute the function, and encode inputs which lead to no change by -1.
  f.prelim <- apply(X, 1, function(x) sum(x * a))
  unchanged <- (f.prelim == 0)
  f <- as.numeric(f.prelim > 0)
  f[unchanged] <- -1
  return(f)
}

##' Simulate a threshold PBN.
##'
##'
##' From a given starting configurations, simulate the network Emax times for
##' Tmax steps in every single timeseries.
##'
##' @param net The network from which the timeseries is generated.
##' @param Tmax The number of time points per timeseries.
##' @param Emax The number of timeseries.
##' @param start If "random", gene values are initialized randomly.
##' @param start.values If start is not "random" these start values are used.
##' @param fix.genes Vector of genes to be kept fixed at a certain level.
##' @param fix.values Vector of values at which genes should be fixed.
##' @return A list of timeseries.
##' @examples
##' net <- createNetwork(inputProbabilities=1, 5, 2, "homogeneous")
##' ts.multi <- simulateNetwork(net, 10, 50) # Creates 50 timeseries à 10 points
simulateNetwork <- function(net, Tmax, Emax, start="random", start.values,
                            fix.genes=-1, fix.values=-1){
  n <- length(net$genes)
  ts.multi <-
    lapply(1:Emax,
           function(e){
             timeseries <- matrix(nrow=n, ncol=Tmax)
             if (start == "random"){
               timeseries[, 1] <- sample(c(0, 1), size=n, replace=TRUE)
             } else {
               timeseries[, 1] <- start.values
             }
             for (t in 2:Tmax){
               timeseries[, t] <-
                 sapply(1:n, function(g) predict(net, g, timeseries[, (t-1)]))
             }
             return(timeseries)
           }
           )
  return(ts.multi)
}

##' Generate gene's next state given current state.
##'
##'
##' Generate gene's next state given current state.
##'
##' @param net The network to be simulated.
##' @param g The gene whose value is to be predicted.
##' @param state The current state of the netwrok.
##' @return The simulated value of the gene.
predict <- function(net, g, state){
  output <- -1
  interactions <- net$interactions[[g]]

  ## Sample one of the controlling function according to its probability
  c <- sapply(interactions, function(x) x$probability)
  i <- sample(1:length(interactions), size=1, prob=c)

  ## Compute the output value
  input <- interactions[[i]]$input
  inputValues <- state[input]
  index <- unbinary(paste(inputValues, collapse="")) + 1
  outputOp <- interactions[[i]]$func[index]
  if (outputOp == -1){
    output <- state[g]
  } else {
    output <- outputOp
  }
  return(output)
}