get.sim.trajectories {Karen}

Simulate a clonal tracking dataset from a given cell differentiation network.

Description

This function simulates clone-specific trajectories for a cell differentiation network associated to a set of (constrained) biochemical reactions, cell types, and missing/latent cell types.

Usage

get.sim.trajectories(
  rct.lst,
  constr.lst = NULL,
  latSts.lst,
  ct.lst,
  th,
  S,
  nCL,
  X0,
  s2 = 1e-08,
  r0 = 0,
  r1 = 0,
  f = 0,
  ntps,
  trunc = FALSE
)

Arguments

rct.lst

A list of biochemical reactions defining the cell differentiation network. A differentiation move from cell type "A" to cell type "B" must be coded as "A->B" Duplication of cell "A" must be coded as "A->1" Death of cell "A" must be coded as "A->0".

constr.lst

(defaults to NULL, when no constraints are needed) List of linear constraints that must be applied to the biochemical reactions. For example, if we need the constraint "A->B = B->C + B->D", this must be coded using the following syntax c("theta\[\'A->B\'\]=(theta\[\'B->C\'\] + theta\[\'B->D\'\])").

latSts.lst

List of the latent cell types. If for example counts are not available for cell types "A" and "B", then latSts.lst = c("A", "B").

ct.lst

List of all the cell types involved in the network formulation. For example, if the network is defined by the biochemical reactions are A->B" and "A->C", then ct.lst = c("A", "B", "C").

th

The vector parameter that must be used for simulation. The length of th equals the number of unconstrained reactions plus 2 (for the noise parameters (\rho_0, \rho_1)). Only positive parameters can be provided.

S

The length of each trajectory.

nCL

An integer defining the number of distinct clones.

X0

A p-dimensional vector for the initial condition of the cell types, where p is the number of distinct cell types provided in ct.lst.

s2

(defaults to 1e-8) A positive value for the overall noise variance.

r0

(defaults to 0) A positive value for the intercept defining the noise covariance matrix R_k = \rho_0 + \rho_1G_kX_k).

r1

(defaults to 0) A positive value for the slope defining the noise covariance matrix R_k = \rho_0 + \rho_1G_kX_k).

f

(defaults to 0) The fraction of measurements that must be considered as missing/latent.

ntps

Number of time points to consider from the whole simulated clonal tracking dataset.

trunc

(defaults to FALSE) Logical, indicating whether sampling from a truncated multivariate normal must be performed.

Value

A list containing the following:

• "X"The simulated process.

• "Y"The simulated noisy-corrupted measurements.

Examples

rcts <- c("HSC->T", ## reactions
          "HSC->M",
          "T->0",
          "M->0")

cnstr <- NULL
latsts <- "HSC" ## latent cell types

ctps <- unique(setdiff(c(sapply(rcts, function(r){ ## all cell types
  as.vector(unlist(strsplit(r, split = "->", fixed = TRUE)))
}, simplify = "array")), c("0", "1")))

## simulation parameters:
S <- 100 ## trajectories length
nCL <- 2 ## number of clones
X0 <- rep(0, length(ctps)) ## initial condition
names(X0) <- ctps
X0["HSC"] <- 1
ntps <- 5 ## number of time-points
f_NA <- 0 ## fraction of observed data

th.true <- c(1.9538674, 1.0559815, 0.7232172, 0.7324133) ## dynamic parameters
names(th.true) <- rcts
s2.true <- 1e-8 ## additonal noise
r0.true <- .1 ## intercept noise parameter
r1.true <- .01 ## slope noise parameter

## simulate trajectories:
XY <- get.sim.trajectories(rct.lst = rcts,
                           constr.lst = cnstr,
                           latSts.lst = latsts,
                           ct.lst = ctps,
                           th = th.true,
                           S = S,
                           nCL = nCL,
                           X0 = X0,
                           s2 = s2.true,
                           r0 = r0.true,
                           r1 = r1.true,
                           f = f_NA,
                           ntps = ntps,
                           trunc = FALSE)

XY$X ## process
XY$Y ## measurements