Defining a disease model
Suppose that there are \(n\) genetic mutations of interests. A fully expressive model would need to define \(2^n\) birth rates corresponding to each possible subsets of mutations that can be present in a cell’s genome. While this is possible for small \(n\), it quickly becomes intractable. For this reason, we choose a more restricting parameterization.
Our model assumes that there is a directed acyclic graph (DAG) \(G\) governing the order in which mutations can occur. The vertices of \(G\) are the mutations and the edges constrain the order in which the mutations occur. In particular, suppose that a cell’s genotype contains mutations \(V' \subset V\). During division, we assume that a cell can acquire any mutation \(m\) for which there is an edge \((v',m) \in E\) such that \(v' \in V'\). Furthermore, we assume that each edge \(e = (v,w)\) has a mutation rate \(u_e\) which gives the probability that a cell with mutation \(v\) acquires mutation \(w\) during division. When crossing edge \(e\), the birth rate of the cell multiplied by \(s_e\) (the selective advantage).
To summarize the above paragraph, we require a DAG \(G\), probabilities for cross each edge \(u_e\), and selective advantages for crossing each edge \(s_e\).
The simplest way to encode the above information is with a matrix with \(|E|\) rows and \(4\) columns. Function progressionChain() generates a simple example where \(G\) is a linear chain with \(n\) vertices:
library(SITH)
set.seed(278115205)
n <- 5
G <- progressionChain(5)
G
#> Head Tail Mut. rate Selective advantage
#> [1,] 0 1 0.001 1
#> [2,] 1 2 0.001 1
#> [3,] 2 3 0.001 1
#> [4,] 3 4 0.001 1
#> [5,] 4 5 0.001 1
The matrix returned by progressionChain() is the input that is required to simulateTumor(). Of course, the user can change head and tail to create more complicated DAGs.
Once we have obtained a matrix of the desired form, we can plug it right into the disease model parameter simulateTumor().
G[,3] <- rep(0.02, nrow(G))
out <- simulateTumor(max_pop = 1000, disease_model=G)
#> Initializing structures ... ...
#> Simulation complete. Releasing memory ... ...
#> Simulated time is 83.5229 days
#> Simulation completed in 0.047088 s.
#> Writing results ... ...
The out list is exactly the same as the returned by simulateTumor(). In particular all of the downstream functions described in the first vignette apply to the output of simulateTumorUDT().
We defined the disease model G to be a linear chain with \(5\) vertices. We should thus expect mutation \(0\) to be less frequent than mutation \(1\), since \(0\) is required for \(1\) to occur.
out$muts
#> id count MAF
#> 1 0 1000 1.000
#> 2 1 323 0.323
#> 3 2 12 0.012
Let’s also take a look at a 2D cross section
Converting from igraph
To define complicated networks, it is probably easier to obtain the matrix G from an igraph object. We provide functionality to do this. Let’s start by whipping up a simple DAG:
if (!requireNamespace("igraph", quietly = TRUE)) {
stop("igraph needed to compile the vignette.",
call. = FALSE)
}
el <- matrix(as.character(c(0,1,0,2,2,3,3,4,2,5,1,5)), ncol = 2, byrow = TRUE)
iG <- igraph::graph_from_edgelist(el)
plot(iG)
We’ll assume that \(0\) is the initial mutation (similar to the original model).
Now we can convert this object into a matrix that can be used by simulateTumor() using progressionDAG_from_igraph().
G <- progressionDAG_from_igraph(iG)
G
#> Head Tail Mut. rate Selective advantage
#> [1,] 0 1 0.001 1
#> [2,] 0 2 0.001 1
#> [3,] 2 3 0.001 1
#> [4,] 3 4 0.001 1
#> [5,] 2 5 0.001 1
#> [6,] 1 5 0.001 1
And now we can simulate a tumor according to these constraints:
out <- simulateTumor(max_pop=5000,disease_model=G, verbose = FALSE)
out$genotypes
#> X1 X2 X3 count
#> 1 0 -1 -1 4744
#> 2 0 1 -1 176
#> 3 0 2 -1 73
#> 4 0 1 2 7
#> 5 0 1 5 0
Modeling loss-of-function mutations
A loss-of-function mutation usually requires both copies of the gene to be mutated. A tumor cell may gain a selective advantage if a loss-of-function mutation occurs in the tumor suppressor gene TP53. However, it is possible that the cell is at a disadvantage if only one of the two copies are mutated. This might be referred to as “crossing the fitness valley.”
simulateTumor() is well-suited for studying a model of loss-of-function mutations.
We can model mutations in TP53 as follows: state \(0\) corresponds to the wild-type allele, state \(1\) corresponds to a cell with \(1\) copy of TP53, and cell \(2\) corresponds to a cell with both copies of TP53 mutated. Furthermore, we will assume that a cell with only \(1\) mutated copy of TP53 has a selective disadvantage. However, a cell with mutations in both copies of TP53 has a selective advantage. For simplicity, let’s assume that both copies can’t be lost in a single event, so that state \(2\) cells can only arise from state \(1\) cells. Using the functionality described above, we can simulate this model in just a few lines of code.
# Make a chain: 0 -> 1 -> 2
G <- progressionChain(2)
#Give a disadvantage to a cell in state 1 and an advantage to a cell in state 2
G[1,4] <- 0.9; G[2,4] <- 1.4
G[,3] <- c(0.005,0.005)
G
#> Head Tail Mut. rate Selective advantage
#> [1,] 0 1 0.005 0.9
#> [2,] 1 2 0.005 1.4
#Simulate the model with N = 250000 cells
out <- simulateTumor(disease_model=G)
#> Initializing structures ... ...
#> Simulated time: 301.088 days. Population is 86602 cells.
#> Simulated time: 335.081 days. Population is 157408 cells.
#> Simulated time: 356.48 days. Population is 223625 cells.
#> Simulation complete. Releasing memory ... ...
#> Simulated time is 363.48 days
#> Simulation completed in 2.40094 s.
#> Writing results ... ...
Let’s see what a cross section looks like:
Cells with one copy of the gene mutated appear to spread out somewhat evenly while cells with both copies mutated are clustered together. Remember, if you are ever confused about which color is which just check:
out$color_scheme
#> [1] "#808080FF" "#934A4EFF" "#E191D5FF"
Let’s also check the MAFs
out$muts
#> id count MAF
#> 1 0 250000 1.000000
#> 2 1 166007 0.664028
#> 3 2 161213 0.644852
out$genotypes
#> X1 X2 X3 count
#> 1 0 -1 -1 83993
#> 2 0 1 -1 4794
#> 3 0 1 2 161213
Of course, all of the other functions that we used in the introduction vignette can be applied to analyze this model as well.
