inst/doc/CancerInSilico.R

In CancerInSilico: An R interface for computational modeling of tumor progression

## ----include=FALSE, cache=FALSE-----------------------------------------------
library(CancerInSilico)
library(gplots)
library(Rtsne)
library(viridis)
library(rgl)

## -----------------------------------------------------------------------------
simple_mod <- suppressMessages(inSilicoCellModel(initialNum=30, runTime=72,
    density=0.1, outputIncrement=24, randSeed=123))

## -----------------------------------------------------------------------------
plotCells(simple_mod, time=0)
plotCells(simple_mod, time=36)
plotCells(simple_mod, time=72)

## -----------------------------------------------------------------------------
# hours in simulation
times <- 0:simple_mod@runTime

# plot number of cells over time
nCells <- sapply(times, getNumberOfCells, model=simple_mod)
plot(times, nCells, type="l", xlab="hour", ylab="number of cells")

# plot population density over time
den <- sapply(times, getDensity, model=simple_mod)
plot(times, den, type="l", xlab="hour", ylab="population density")

## -----------------------------------------------------------------------------
drug <- new("Drug", name="Drug_A", timeAdded=24,
    cycleLengthEffect=function(type, length) length * 2)
drug_mod <- suppressMessages(inSilicoCellModel(initialNum=30, runTime=72,
    density=0.1, drugs=c(drug), outputIncrement=24, randSeed=123))

# hours in simulation
times <- 0:simple_mod@runTime

# plot number of cells over time
nCells <- sapply(times, getNumberOfCells, model=simple_mod)
nCells_drug <- sapply(times, getNumberOfCells, model=drug_mod)
plot(times, nCells, type="l", xlab="hour", ylab="number of cells")
lines(times, nCells_drug, type="l", xlab="hour", ylab="number of cells",
    col="red")

## -----------------------------------------------------------------------------
type_A <- new("CellType", name="A", minCycle=16, cycleLength=function() 16)
fast_cells_mod <- suppressMessages(inSilicoCellModel(initialNum=30, runTime=72,
    density=0.1, cellTypes=c(type_A), outputIncrement=24, randSeed=123))

# hours in simulation
times <- 0:fast_cells_mod@runTime

# plot number of cells over time
nCells <- sapply(times, getNumberOfCells, model=simple_mod)
nCells_fast <- sapply(times, getNumberOfCells, model=fast_cells_mod)
plot(times, nCells, type="l", xlab="hour", ylab="number of cells")
lines(times, nCells_fast, type="l", xlab="hour", ylab="number of cells",
    col="red")

## -----------------------------------------------------------------------------
type_B <- new("CellType", name="B", size=1, minCycle=16,
    cycleLength=function() 16 + rexp(1,1/4))
type_C <- new("CellType", name="C", size=1, minCycle=32,
    cycleLength=function() 32 + rexp(1,1/4))
two_types_mod <- suppressMessages(inSilicoCellModel(initialNum=30, runTime=72,
    density=0.1, cellTypes=c(type_B, type_C), cellTypeInitFreq=c(0.4,0.6),
    outputIncrement=24, randSeed=123))

## -----------------------------------------------------------------------------
getTypeBProportion <- function(time)
{
    N <- getNumberOfCells(two_types_mod, time)
    sum(sapply(1:N, function(i) getCellType(two_types_mod, time, i) == 1)) / N
}
times <- 0:two_types_mod@runTime
Bprop <- sapply(times, getTypeBProportion)
plot(times, Bprop, type="l", xlab="hour", ylab="type B proportion")

## -----------------------------------------------------------------------------
mitosisGeneNames <- paste("m_", letters[1:20], sep="")
mitosisExpression <- function(model, cell, time)
{
    ifelse(getCellPhase(model, time, cell) == "M", 1, 0)
}

pwyMitosis <- new("Pathway", genes=mitosisGeneNames,
    expressionScale=mitosisExpression)

## -----------------------------------------------------------------------------
contactInhibitionGeneNames <- paste("ci_", letters[1:15], sep="")
contactInhibitionExpression <- function(model, cell, time)
{
    getLocalDensity(model, time, cell, 3.3)
}
pwyContactInhibition <- new("Pathway", genes=contactInhibitionGeneNames,
    expressionScale=contactInhibitionExpression)

## -----------------------------------------------------------------------------
# create simulated data set
allGenes <- c(mitosisGeneNames, contactInhibitionGeneNames)
geneMeans <- 2 + rexp(length(allGenes), 1/20)
data <- t(pmax(sapply(geneMeans, rnorm, n=25, sd=2), 0))
rownames(data) <- allGenes

# calibrate pathways
pwyMitosis <- calibratePathway(pwyMitosis, data)
pwyContactInhibition <- calibratePathway(pwyContactInhibition, data)

## -----------------------------------------------------------------------------
params <- new("GeneExpressionParams")
params@randSeed <- 123 # control this for reporducibility
params@nCells <- 30 # sample 30 cells at each time point to measure activity
params@sampleFreq <- 6 # measure activity every 6 hours

pwys <- c(pwyMitosis, pwyContactInhibition)
pwyActivity <- inSilicoGeneExpression(simple_mod, pwys, params)$pathways

## -----------------------------------------------------------------------------
# mitosis
plot(seq(0,72,6), pwyActivity[[1]], type="l", col="orange", ylim=c(0,1))
# contact inhibition
lines(seq(0,72,6), pwyActivity[[2]], col="blue")

## -----------------------------------------------------------------------------
pwyMitosis@expressionScale = function(model, cell, time)
{
    window <- c(max(time - 2, 0), min(time + 2, model@runTime))
    a1 <- getAxisLength(model, window[1], cell)
    a2 <- getAxisLength(model, window[2], cell)
    if (is.na(a1)) a1 <- 0 # in case cell was just born
    return(ifelse(a2 < a1, 1, 0))
}
pwys <- c(pwyMitosis, pwyContactInhibition)
pwyActivity <- inSilicoGeneExpression(simple_mod, pwys, params)$pathways
# mitosis
plot(seq(0,72,6), pwyActivity[[1]], type="l", col="orange", ylim=c(0,1))
# contact inhibition
lines(seq(0,72,6), pwyActivity[[2]], col="blue")

## -----------------------------------------------------------------------------
pwyMitosis@transformMidpoint = 0.1  
pwyMitosis@transformSlope = 5 / 0.1
pwys <- c(pwyMitosis, pwyContactInhibition)
pwyActivity <- inSilicoGeneExpression(simple_mod, pwys, params)$pathways
# mitosis
plot(seq(0,72,6), pwyActivity[[1]], type="l", col="orange", ylim=c(0,1))
# contact inhibition
lines(seq(0,72,6), pwyActivity[[2]], col="blue")

## -----------------------------------------------------------------------------
params@RNAseq <- FALSE # generate microarray data
params@singleCell <- FALSE # generate bulk data
params@perError <- 0.1 # parameter for simulated noise

pwys <- c(pwyMitosis, pwyContactInhibition)
ge <- inSilicoGeneExpression(simple_mod, pwys, params)$expression

## -----------------------------------------------------------------------------
ndx <- apply(ge, 1, var) == 0 # remove zero variance rows
heatmap.2(ge[!ndx,], 
    col=greenred, scale="row",
    trace="none", hclust=function(x) hclust(x,method="complete"),
    distfun=function(x) as.dist((1-cor(t(x)))/2), 
    Colv=FALSE, dendrogram="row",
    RowSideColors = ifelse(rownames(ge[!ndx,]) %in%
        mitosisGeneNames, "orange", "blue"),
    labRow = FALSE, labCol = seq(0,72,6),
    main="Bulk Gene Expression from Simple Cell Simulation")

## -----------------------------------------------------------------------------
# gene names
B_genes <- paste("b.", letters[1:20], sep="")
C_genes <- paste("c.", letters[1:20], sep="")

# pathway behavior
pwy_B <- new("Pathway", genes=B_genes, expressionScale=
    function(model, cell, time) ifelse(getCellType(model, time, cell)==1, 1, 0))
pwy_C <- new("Pathway", genes=C_genes, expressionScale=
    function(model, cell, time) ifelse(getCellType(model, time, cell)==2, 1, 0))

# calibrate pathways
geneMeans <- 2 + rexp(length(c(B_genes, C_genes)), 1/20)
data <- t(pmax(sapply(geneMeans, rnorm, n=25, sd=2), 0))
rownames(data) <- c(B_genes, C_genes)
pwy_B <- calibratePathway(pwy_B, data)
pwy_C <- calibratePathway(pwy_C, data)

## -----------------------------------------------------------------------------
params@RNAseq <- TRUE
params@singleCell <- TRUE
params@dropoutPresent <- TRUE
ge <- inSilicoGeneExpression(two_types_mod, c(pwy_B, pwy_C), params)$expression

## -----------------------------------------------------------------------------
cells <- unname(sapply(colnames(ge), function(x) strsplit(x,"_")[[1]][1]))
cells <- as.numeric(gsub("c", "", cells))
type <- sapply(cells, getCellType, model=two_types_mod,
    time=two_types_mod@runTime)
type[type==1] <- "red"
type[type==2] <- "blue"

pca <- prcomp(ge, center=FALSE, scale.=FALSE)
plot(pca$rotation[,c(1,2)], col=type)