## ----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)