R/e06-engineAct.R

In Umpire: Simulating Realistic Gene Expression and Clinical Data

###############################################################
# An ENGINE WITH ACTIVITY allows for the possibility that some
# components (or genes) in an expression engine (or tissue) might
# be transcriptionally inactive.  Thus, the true biological signal
# described previously should really be viewed as a mixture
#    S_gi = z_g * delta_0 + (1 - z_g) * T_gi
# where
#    delta_0 = a point mass at zero
#    T_gi = a random variable supported on the positive real line
#    z_g ~ Binom(pi) defines the activity state (1 = on, 0 = off)
# Note that we typically work not with T_gi but with its logarithm
# to some appropriate base.  That is, the multivariate normal or
# independent normal blocks used to construct engines should be
# applied on the logarithmic scale.

setClass("EngineWithActivity",
         contains = "Engine",
         slots = c(active="logical",
                   base="numeric"))

## Generates an EngineWithActivity object.
EngineWithActivity <- function(active, components, base=2) {
  e <- Engine(components)
  if (length(active) == 1) {
    active <- rbinom(nComponents(e), 1, active)==1
  }
  new("EngineWithActivity", e, active=active, base=base)
}

setValidity("EngineWithActivity", function(object) {
  msg <- NULL
  rightSize <- length(object@active) == nComponents(object)
  if (!rightSize) {
    msg <- c(msg, "number of object components not equal length of active")
  }
  if (any(object@base < 0)) {
    msg <- c(msg, "base is negative")
  }
  if (is.null(msg)) { # pass
    msg <- TRUE
  }
  msg
})

# The 'rand' method for an EngineWithActivity is a little bit
# tricky, since we do two things at once. First, we use the
# 'base' slot to exponentiate the random variables generated by
# the underlying Engine on the log scale.  We treat base = 0 as
# a special case, which means that we should continue to work on
# the scale of the Engine.  Second, we mask any inactive component
# by replacing the generated values with 0.
#
# Note that this is terribly inefficient if we only have a single
# homogeneous population, since we generate a certain amount of
# data only to throw it away.  The power comes when we allow
# cancer disregulation to turn a block on or off, when the
# underlying data reappears.

setMethod("rand", "EngineWithActivity", function(object, n, ...) {
  x <- callNextMethod()
  if (object@base > 0) { # exponentiate the log signal
    x <- object@base^x
  }
  blockSizes <- unlist(lapply(object@components, nrow))
  pi <- rep(object@active, times=blockSizes)
  x * pi # mask signals from inactive gene-blocks
})

setMethod("summary", "EngineWithActivity", function(object, ...) {
  callNextMethod()
  cat(paste("Fraction of active genes", sum(object@active)))
})