_CancerInSilico-1.Rmd

---
title: "CancerInSilico"
output: html_notebook
---

An R interface for computational modeling of tumor progression

https://bioconductor.org/packages/release/bioc/html/CancerInSilico.html



```r
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("CancerInSilico", version = "3.8")
library(CancerInSilico)
```



```r
browseVignettes("CancerInSilico")
```


https://bioconductor.org/packages/release/bioc/vignettes/CancerInSilico/inst/doc/CancerInSilico.html

# Running a Cell Simulation



## Run Simple Simulation

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

## Plot CellModel Object



```r
plotCells(simple_mod, time=0)
```

```r
plotCells(simple_mod, time=36)
```

```r
plotCells(simple_mod, time=72)
```


## Query Cell Information




```r
# 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")
```

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

# Drugs

```r
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))
```



```r
# 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")
```


# Cell Types

## Adding a Single Cell Type



```r
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))
```


```r
# 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")
```


## Adding Multiple Cell Types



```r
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))
```


## Getting Cell Type



```r
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")
```


# Pathways


```r
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)
```



```r
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)
```


## Calibrate Gene Expression Range

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


## Generate Pathway Activity

```r
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
```


## Visualize Pathway Activity

```r
# 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")
```


## Accounting for Model Effects




```r
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")
```



## Normalize Pathway Activity


```r
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")
```


# Simulating Bulk Gene Expression Data

## Simulating Microarray Data

```r
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
```

## Visualize Bulk Gene Expression Data

```r
ndx <- apply(ge, 1, var) == 0 # remove zero variance rows
gplots::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")
```

# Simulating Single Cell Gene Expression Data

## Cell Type Pathways


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


## Simulating Single Cell RNA-seq

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


## Visualize Single Cell Data

```r
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)
```