1.3.1 Recovering restrictions in the order of accumulation of mutations

This is a question that was addressed, for instance, in Diaz-Uriarte (2015): do methods that try to infer restrictions in the order of accumulation of mutations (e.g., Szabo and Boucher 2008; Gerstung et al. 2009; Ramazzotti et al. 2015) work well under different evolutionary models and with different sampling schemes? (This issue is also touched upon in section 13.2).

A possible way to examine that question would involve:

• generating random DAGs that encode restrictions;
• simulating cancer evolution using those DAGs;
• sampling the data and adding different levels of noise to the sampled data;
• running the inferential method;
• comparing the inferred DAG with the original, true, one.

## For reproducibility
set.seed(2)
RNGkind("L'Ecuyer-CMRG")

## Simulate a DAG
g1 <- simOGraph(4, out = "rT")

## Simulate 10 evolutionary trajectories
s1 <- oncoSimulPop(10, allFitnessEffects(g1, drvNames = 1:4),
                   mc.cores = 2, ## adapt to your hardware
                   seed = NULL) ## for reproducibility of vignette

## Sample those data uniformly, and add noise
d1 <- samplePop(s1, timeSample = "unif", propError = 0.1)
## 
##  Subjects by Genes matrix of 10 subjects and 4 genes.

## You would now run the appropriate inferential method and
## compare observed and true. For example

## require(Oncotree)
## fit1 <- oncotree.fit(d1)

## Now, you'd compare fitted and original. This is well beyond 
## the scope of this document (and OncoSimulR itself).

1.3.2 Sign epistasis and probability of crossing fitness valleys

This question, and the question in the next section (1.3.3), encompass a wide range of issues that have been addressed in evolutionary genetics studies and which include from detailed analysis of simple models with a few uphill paths and valleys as in Weissman et al. (2009) or Ochs and Desai (2015), to questions that refer to larger, more complex fitness landscapes as in Szendro, Franke, et al. (2013) or Franke et al. (2011) (see below).

Using as an example Ochs and Desai (2015) (we will see this example again in section 5.3, where we cover different ways of specifying fitness), we could specify the fitness landscape and run simulations until fixation (with argument fixation to oncoSimulPop —see more details in section 6.3.3 and 6.3.4, again with this example). We would then examine the proportion of genotypes fixed under different scenarios. And we can extend this example by adding mutator genes:

## For reproducibility
set.seed(2)
RNGkind("L'Ecuyer-CMRG")

## Specify fitness effects. 

## Numeric values arbitrary, but set the intermediate genotype en
## route to ui as mildly deleterious so there is a valley.

## As in Ochs and Desai, the ui and uv genotypes
## can never appear. 

u <- 0.2; i <- -0.02; vi <- 0.6; ui <- uv <- -Inf

od <- allFitnessEffects(
    epistasis = c("u" = u,  "u:i" = ui,
                  "u:v" = uv, "i" = i,
                  "v:-i" = -Inf, "v:i" = vi))

## For the sake of extending this example, also turn i into a
## mutator gene

odm <- allMutatorEffects(noIntGenes = c("i" = 50))

## How do mutation and fitness look like for each genotype?
evalAllGenotypesFitAndMut(od, odm, addwt = TRUE)
##   Genotype Fitness MutatorFactor
## 1       WT   1.000             1
## 2        i   0.980            50
## 3        u   1.200             1
## 4        v   0.000             1
## 5     i, u   0.000            50
## 6     i, v   1.568            50
## 7     u, v   0.000             1
## 8  i, u, v   0.000            50

Ochs and Desai explicitly say “Each simulated population was evolved until either the uphill genotype or valley-crossing genotype fixed.” So we will use fixation.

## Set a small initSize, as o.w. unlikely to pass the valley
initS <- 10
## The number of replicates is tiny, 10, for the sake of speed
## of creation of the vignette
od_sim <- oncoSimulPop(10, od, muEF = odm,
                       fixation = c("u", "i, v"), initSize = initS,
                       model = "McFL",
                       mu = 1e-4, detectionDrivers = NA, 
                       finalTime = NA,
                       detectionSize = NA, detectionProb = NA,
                       onlyCancer = TRUE, 
                       mc.cores = 2, ## adapt to your hardware
                       seed = NULL) ## for reproducibility
## What is the frequency of each final genotype?
sampledGenotypes(samplePop(od_sim))
## 
##  Subjects by Genes matrix of 10 subjects and 3 genes.
##   Genotype Freq
## 1     i, v    4
## 2        u    6
## 
##  Shannon's diversity (entropy) of sampled genotypes: 0.6730117

1.3.3 Predictability of evolution in complex fitness landscapes

Focusing now on predictability in more general fitness landscapes, we would run simulations under random fitness landscapes with varied ruggedness, and would then examine the evolutionary predictability of the trajectories with measures such as “Lines of Descent” and “Path of the Maximum” (Szendro, Franke, et al. 2013) and the diversity of the sampled genotypes under different sampling regimes (see details in section 10).

## For reproducibility
set.seed(7)
RNGkind("L'Ecuyer-CMRG")


## Repeat the following loop for different combinations of whatever
## interests you, such as number of genes, or distribution of the
## c and sd (which affect how rugged the landscape is), or 
## reference genotype, or evolutionary model, or stopping criterion, 
## or sampling procedure, or ...

##  Generate a random fitness landscape, from the Rough Mount
##  Fuji model, with g genes, and c ("slope" constant) and
##  reference chosen randomly (reference is random by default and 
##  thus not specified below). Require a minimal number of  
##  accessible genotypes

g <- 6
c <- runif(1, 1/5, 5)
rl <- rfitness(g, c = c, min_accessible_genotypes = g)

## Plot it if you want; commented here as it takes long for a
## vignette

## plot(rl)

## Obtain landscape measures from Magellan. Export to Magellan
to_Magellan(rl, file = "rl1.txt")

## (Getting the statistics from Magellan requires either 
## calling the web app (http://wwwabi.snv.jussieu.fr/public/Magellan/) 
## or asking Magellan's authors for the software. This is of course
## beyond the scope of the example and the package.)

## Simulate evolution in that landscape many times (here just 10)
simulrl <- oncoSimulPop(10, allFitnessEffects(genotFitness = rl),
                        keepPhylog = TRUE, keepEvery = 1,
                        initSize = 4000,
                        seed = NULL, ## for reproducibility
                        mc.cores = 2) ## adapt to your hardware

## Obtain measures of evolutionary predictability
diversityLOD(LOD(simulrl))
## [1] 1.05492
diversityPOM(POM(simulrl))
## [1] 1.05492
sampledGenotypes(samplePop(simulrl, typeSample = "whole"))
## 
##  Subjects by Genes matrix of 10 subjects and 6 genes.
##   Genotype Freq
## 1        A    4
## 2     A, F    1
## 3     E, F    4
## 4        F    1
## 
##  Shannon's diversity (entropy) of sampled genotypes: 1.19355

1.3.4 Mutator and antimutator genes

The effects of mutator and antimutator genes have been examined both in cancer genetics (Nowak 2006; Tomlinson, Novelli, and Bodmer 1996) and in evolutionary genetics (Gerrish et al. 2007), and are related to wider issues such as Muller’s ratchet and the evolution of sex. There are, thus, a large range of questions related to mutator and antimutator genes.

One question addressed in Tomlinson, Novelli, and Bodmer (1996) concerns under what circumstances mutator genes are likely to play a role in cancer progression. For instance, Tomlinson, Novelli, and Bodmer (1996) find that an increased mutation rate is more likely to matter if the number of required mutations in driver genes needed to reach cancer is large and if the mutator effect is large.

We might want to ask, then, how long it takes before to reach cancer under different scenarios. Time to reach cancer is stored in the component FinalTime of the output. We would specify different numbers and effects of mutator genes (argument muEF). We would also change the criteria for reaching cancer and in our case we can easily do that by specifying different numbers in detectionDrivers. Of course, we would also want to examine the effects of varying numbers of mutators, drivers, and possibly fitness consequences of mutators. Below we assume mutators are neutral and we assume there are no additional genes with deleterious mutations, but this need not be so, of course (see also Tomlinson, Novelli, and Bodmer 1996; Gerrish et al. 2007; Christopher McFarland, Mirny, and Korolev 2014).

Let us run an example. For the sake of simplicity, we assume no epistatic interactions.

sd <- 0.1 ## fitness effect of drivers
sm <- 0 ## fitness effect of mutator
nd <- 20 ## number of drivers
nm <- 5  ## number of mutators
mut <- 10 ## mutator effect

fitnessGenesVector <- c(rep(sd, nd), rep(sm, nm))
names(fitnessGenesVector) <- 1:(nd + nm)
mutatorGenesVector <- rep(mut, nm)
names(mutatorGenesVector) <- (nd + 1):(nd + nm)

ft <- allFitnessEffects(noIntGenes = fitnessGenesVector,
                        drvNames = 1:nd)
mt <- allMutatorEffects(noIntGenes = mutatorGenesVector)

Now, simulate using the fitness and mutator specification. We fix the number of drivers to cancer, and we stop when those numbers of drivers are reached. Since we only care about the time it takes to reach cancer, not the actual trajectories, we set keepEvery = NA:

## For reproducibility
set.seed(2)
RNGkind("L'Ecuyer-CMRG")


ddr <- 4
st <- oncoSimulPop(4, ft, muEF = mt,
                   detectionDrivers = ddr,
                   finalTime = NA,
                   detectionSize = NA,
                   detectionProb = NA,
                   onlyCancer = TRUE,
                   keepEvery = NA, 
                   mc.cores = 2, ## adapt to your hardware
                   seed = NULL) ## for reproducibility

## How long did it take to reach cancer?
unlist(lapply(st, function(x) x$FinalTime))
## [1]  370  141 1793  282

(Incidentally, notice that it is easy to get OncoSimulR to throw an exception if you accidentally specify a huge mutation rate when all mutator genes are mutated: see section 12.7.)

1.3.5 Epistatic interactions between drivers and passengers in cancer and the consequences of order effects

K <- 4
sp <- 1e-5
sdp <- 0.015
sdplus <- 0.05
sdminus <- 0.1

cnt <- (1 + sdplus)/(1 + sdminus)
prod_cnt <- cnt - 1
bauer <- data.frame(parent = c("Root", rep("D", K)),
                    child = c("D", paste0("s", 1:K)),
                    s = c(prod_cnt, rep(sdp, K)),
                    sh = c(0, rep(sp, K)),
                    typeDep = "MN")
fbauer <- allFitnessEffects(bauer)
(b1 <- evalAllGenotypes(fbauer, order = FALSE, addwt = TRUE))
##             Genotype   Fitness
## 1                 WT 1.0000000
## 2                  D 0.9545455
## 3                 s1 1.0000100
## 4                 s2 1.0000100
## 5                 s3 1.0000100
## 6                 s4 1.0000100
## 7              D, s1 0.9688636
## 8              D, s2 0.9688636
## 9              D, s3 0.9688636
## 10             D, s4 0.9688636
## 11            s1, s2 1.0000200
## 12            s1, s3 1.0000200
## 13            s1, s4 1.0000200
## 14            s2, s3 1.0000200
## 15            s2, s4 1.0000200
## 16            s3, s4 1.0000200
## 17         D, s1, s2 0.9833966
## 18         D, s1, s3 0.9833966
## 19         D, s1, s4 0.9833966
## 20         D, s2, s3 0.9833966
## 21         D, s2, s4 0.9833966
## 22         D, s3, s4 0.9833966
## 23        s1, s2, s3 1.0000300
## 24        s1, s2, s4 1.0000300
## 25        s1, s3, s4 1.0000300
## 26        s2, s3, s4 1.0000300
## 27     D, s1, s2, s3 0.9981475
## 28     D, s1, s2, s4 0.9981475
## 29     D, s1, s3, s4 0.9981475
## 30     D, s2, s3, s4 0.9981475
## 31    s1, s2, s3, s4 1.0000400
## 32 D, s1, s2, s3, s4 1.0131198

## How does the fitness landscape look like?
plot(b1, use_ggrepel = TRUE) ## avoid overlapping labels

## For reproducibility
set.seed(2)
RNGkind("L'Ecuyer-CMRG")

totalpops <- 5
initSize <- 100
sb1 <- oncoSimulPop(totalpops, fbauer, model = "Exp",
                    initSize = initSize,
                    onlyCancer = FALSE, 
                    mc.cores = 2, ## adapt to your hardware
                    seed = NULL) ## for reproducibility
                    
## What proportion of the simulations reach 4x initSize?
sum(summary(sb1)[, "TotalPopSize"] > (4 * initSize))/totalpops
## [1] 0.2

totalpops <- 5
initSize <- 100
sb2 <- oncoSimulPop(totalpops, fbauer, model = "Exp",
                    initSize = initSize,
                    onlyCancer = TRUE,
                    detectionSize = 10 * initSize,
                    mc.cores = 2, ## adapt to your hardware
                    seed = NULL) ## for reproducibility
                
## How long did it take to reach cancer?
unlist(lapply(sb2, function(x) x$FinalTime))
## [1] 416 354 339 445 215

## Order effects involving three genes.

## Genotype "D, M" has different fitness effects
## depending on whether M or D mutated first.
## Ditto for genotype "F, D, M".

## Meaning of specification: X > Y means
## that X is mutated before Y.


o3 <- allFitnessEffects(orderEffects = c(
                            "F > D > M" = -0.3,
                            "D > F > M" = 0.4,
                            "D > M > F" = 0.2,
                            "D > M"     = 0.1,
                            "M > D"     = 0.5))

## With the above specification, let's double check
## the fitness of the possible genotypes

(oeag <- evalAllGenotypes(o3, addwt = TRUE, order = TRUE))
##     Genotype Fitness
## 1         WT    1.00
## 2          D    1.00
## 3          F    1.00
## 4          M    1.00
## 5      D > F    1.00
## 6      D > M    1.10
## 7      F > D    1.00
## 8      F > M    1.00
## 9      M > D    1.50
## 10     M > F    1.00
## 11 D > F > M    1.54
## 12 D > M > F    1.32
## 13 F > D > M    0.77
## 14 F > M > D    1.50
## 15 M > D > F    1.50
## 16 M > F > D    1.50

## For reproducibility
set.seed(2)
RNGkind("L'Ecuyer-CMRG")

totalpops <- 5
soe1 <- oncoSimulPop(totalpops, o3, model = "Exp",
                    initSize = 500,
                    onlyCancer = FALSE,
                    mc.cores = 2, ## adapt to your hardware
                    seed = NULL) ## for reproducibility
                    
## What proportion of the simulations do not end up extinct?
sum(summary(soe1)[, "TotalPopSize"] > 0)/totalpops
## [1] 0.4

1.7.1 HTML and PDF versions of the vignette

A PDF version of this vignette is available from https://rdiaz02.github.io/OncoSimul/pdfs/OncoSimulR.pdf. And an HTML version from https://rdiaz02.github.io/OncoSimul/OncoSimulR.html. These files should correspond to the most recent, GitHub version, of the package (i.e., they might include changes not yet available from the BioConductor package).

2.1.1 Changing fitness: \(s=0.1\) and \(s=0.05\)

In the above fitness specification the fitness effect of each gene (when its restrictions are satisfied) is \(s = 0.1\) (see section 3.2 for details). Here we rerun all the above benchmarks using \(s= 0.05\) (the results from these benchmarks are available as data(benchmark_1_0.05)) and results are shown below in Table 2.6.

As expected, having a smaller \(s\) leads to slower processes in most cases, since it takes longer to reach the exiting conditions sooner. Particularly noticeable are the runs for the McFL models (notice the increases in population size and number of iterations —see also below).

That is not the case, however, for exp5 and exp6 (and exp5_noc and exp6_noc). When running with \(s=0.05\) the simulations exit at a later time (see column “Final Time”) but they exit with smaller population sizes. Here we have an interaction between sampling frequency, speed of growth of the population, mutation events and number of clones. In populations that grow much faster mutation events will happen more often (which will trigger further iterations of the algorithm); in addition, more new clones will be created, even if they only exist for short times and become extinct by the following sampling period (so they are not reflected in the pops.by.time matrix). These differences are proportionally larger the larger the rate of growth of the population. Thus, they are larger between, say, the exp5 at \(s=0.1\) and \(s=0.05\) than between the exp4 at the two different \(s\): the exp5 exit conditions can only be satisfied at much larger population sizes so at populations sizes when growth is much faster (recall we are dealing with exponential growth).

Recall also that with the default settings in detectionProb, we assess the exiting condition every 20 time periods (argument checkSizePEvery); this means that for fast growing populations, the increase in population size between successive checks of the exit conditions will be much larger (this phenomenon is also discussed in section 2.2.1).

Thus, what is happening in the exp5 and exp6 with \(s=0.1\) is that close to the time the exit conditions could be satisfied, they are growing very fast, accumulating mutants, and incurring in additional iterations. They exit sooner in terms of time periods, but they do much more work before arriving there.

The setting of checkSizePEvery is also having a huge effect on the McFL model simulations (the number of iterations is \(>10^6\)). Even more than in the previous section, checking the exiting condition more often (smaller checkSizePEvery) would probably be justified here (notice also the very large final times) and would lead to a sharp decrease in number of iterations and, thus, running time.

The moral here is that in complex simulations like this (and most simulations are complex), the effects of some parameters (\(s\) in this case) might look counter-intuitive at first. Thus the need to “experiment before launching a large number of simulations”.

2.2.1 Common use cases, set 1.

We will use the Exp and the McFL models, run with different parameters. The script is provided as ‘benchmark_2.R’, under ‘/inst/miscell’, with output in the ‘miscell-files/vignette_bench_Rout’ directory of the main OncoSimul repository at https://github.com/rdiaz02/OncoSimul. The data are available as data(benchmark_2).

For the Exp model the call will be

oncoSimulPop(Nindiv,
            fitness,
            detectionProb = NA, 
            detectionSize = 1e6,
            initSize = 500,
            detectionDrivers = NA,
            keepPhylog = TRUE,
            model = "Exp",
            errorHitWallTime = FALSE,
            errorHitMaxTries = FALSE,
            finalTime = 5000,
            onlyCancer = FALSE,
            mc.cores = 1,
            sampleEvery = 0.5,
            keepEvery = 1)

And for McFL:

initSize <- 1000
oncoSimulPop(Nindiv,
              fitness,
               detectionProb = c(
                   PDBaseline = 1.4 * initSize,
                   n2 = 2 * initSize,
                   p2 = 0.1,
                   checkSizePEvery = 4),
               initSize = initSize,
               detectionSize = NA,
               detectionDrivers = NA,
               keepPhylog = TRUE,
               model = "McFL",
               errorHitWallTime = FALSE,
               errorHitMaxTries = FALSE,
               finalTime = 5000,
               max.wall.time = 10,
               onlyCancer = FALSE,
               mc.cores = 1,
               keepEvery = 1)

For the exponential model we will stop simulations when populations have \(>10^6\) cells (simulations start from 500 cells). For the McFarland model we will use the detectionProb mechanism (see section 6.3.2 for details); we could have used as stopping mechanism detectionSize = 2 * initSize (which would be basically equivalent to reaching cancer, as argued in (McFarland et al. 2013)) but we want to provide further examples under the detectionProb mechanism. We will start from 1000 cells, not 500 (starting from 1000 we almost always reach cancer in a single run).

Why not use the detectionProb mechanism with the Exp models? Because it can be hard to intuitively understand what are reasonable settings for the parameters of the detectionProb mechanism when used in a population that is growing exponentially, especially if different genes have very different effects on fitness. Moreover, we are using fitness specifications that are very different (compare the fitness landscape of six genes, the pancreas specification, and the fitness specification with 4000 genes with fitness effects drawn from an exponential distribution —re_4000). In contrast, the detectionProb mechanism might be simpler to reason about in a population that is growing under a model of carrying capacity with possibly large periods of stasis. Let us emphasize that it is not that the detectionProb mechanism does not make sense with the Exp model; it is simply that the parameters might need finer adjustment for them to make sense, and in these benchmarks we are dealing with widely different fitness specifications.

Note also that we specify checkSizePEvery = 4 (instead of the default, which is 20). Why? Because the fitness specifications where fitness effects are drawn from exponential distributions (re_200 to re_4000 above) include many genes (well, up to 4000) some of them with possibly very large effects. In these conditions, simulations can run very fast in the sense of “units of time”. If we check exiting conditions every 20 units the population could have increased its size several orders of magnitude in between checks (this is also discussed in sections 2.1.1 and 6.3.2). You can verify this by running the script with other settings for checkSizePEvery (and being aware that large settings might require you to wait for a long time). To ensure that populations have really grown, we have increased the setting of PDBaseline so that no simulation can be considered for stopping unless its size is 1.4 times larger than initSize.

In all cases we use keepEvery = 1 and keepPhylog = TRUE (so we store the population sizes of all clones every 1 time unit and we keep the complete genealogy of clones). Finally, we run all models with errorHitWallTime = FALSE and errorHitMaxTries = FALSE so that we can see results even if stopping conditions are not met.

The results of the benchmarks, using 100 individual simulations, are shown in Table 2.7.

In most cases, simulations run reasonably fast (under 0.1 seconds per individual simulation) and the returned objects are small. I will only focus on a few cases.

The McFL model with random fitness landscape rf12 and with pancr does not satisfy the conditions of detectionProb in most cases: its median final time is 5000, which was the maximum final time specified. This suggests that the fitness landscape is such that it is unlikely that we will reach population sizes \(> 1400\) (remember we the setting for PDBaseline) before 5000 time units. There is nothing particular about using a fitness landscape of 12 genes and other runs in other 12-gene random fitness landscapes do not show this pattern. However, complex fitness landscapes might be such that genotypes of high fitness (those that allow reaching a large population size quickly) are not easily accessible³ so reaching them might take a long time. This does not affect the exponential model in the same way because, well, because there is exponential growth in that model: any genotype with fitness \(>1\) will grow exponentially (of course, at possibly very different rates). You might want to play with the script and modify the call to rfitness (using different values of reference and c, for instance) to have simpler paths to a maximum or modify the call to oncoSimulPop (with, say, finalTime to much larger values). Some of these issues are related to more general questions about fitness landscapes and accessibility (see section 1.3.2 and references therein).

You could also set onlyCancer = TRUE. This might make sense if you are interested in only seeing simulations that “reach cancer” (where “reach cancer” means reaching a state you define as a function of population size or drivers). However, if you are exploring fitness landscapes, onlyCancer = TRUE might not always be reasonable as reaching a particular population size, for instance, might just not be possible under some fitness landscapes (this phenomenon is of course not restricted to random fitness landscapes —see also section 2.3.3).

As we anticipated above, the detectionProb mechanism has to be used with care: some of the simulations run in very short “time units”, such as those for the fitness specifications with 2000 and 4000 genes. Having used a checkSizePEvery = 20 probably would not have made sense.

Finally, it is interesting that in the cases examined here, the two slowest running simulations are from “Exp”, with fitnesses re_2000 and re_4000 (and the third slowest is also Exp, under re_500). These are also the cases with the largest number of clones. Why? In the “Exp” model there is no competition, and fitness specifications re_2000 and re_4000 have genomes with many genes with positive fitness contributions. It is thus very easy to obtain, from the wildtype ancestor, a large number of clones all of which have birth rates \(>1\) and, thus, clones that are unlikely to become extinct.

2.2.2 Common use cases, set 2.

We will now rerun the simulations above changing the following:

• finalTime set to 25000.
• onlyCancer set to TRUE.
• The “Exp” models will stop when population size \(> 10^5\).

This is in script ‘benchmark_3.R’, under ‘/inst/miscell’, with output in the ‘miscell-files/vignette_bench_Rout’ directory of the main OncoSimul repository at https://github.com/rdiaz02/OncoSimul. The data are available as data(benchmark_3).

Since we increased the maximum final time and forced runs to “reach cancer” the McFL run with the pancreas fitness specification takes a bit longer because it also has to do a larger number of iterations. Interestingly, notice that the median final time is close to 10000, so the runs in 2.2.1 with maximum final time of 5000 would have had a hard time finishing with onlyCancer = TRUE.

Forcing simulations to “reach cancer” and just random differences between the random fitness landscape also affect the McFL run under rf12: final time is below 5000 and the median number of iterations is about half of what was above.

Finally, by stopping the Exp simulations at \(10^5\), simulations with re_2000 and re_4000 finish now in much shorter times (but they still take longer than their McFL counterparts) and the number of clones created is much smaller.

2.3.1 Exponential model with 10,000 and 50,000 genes

ng <- 10000
u <- allFitnessEffects(noIntGenes = c(rep(0.1, ng/2), 
                                      rep(-0.1, ng/2)))

t_e_10000 <- system.time(
    e_10000 <- oncoSimulPop(5, u, model = "Exp", mu = 1e-7,
                            detectionSize = 1e6,
                            detectionDrivers = NA,
                            detectionProb = NA,
                            keepPhylog = TRUE,
                            onlyCancer = FALSE,
                            mutationPropGrowth = TRUE,
                            mc.cores = 1))

t_e_10000
##    user  system elapsed 
##   4.368   0.196   4.566 
  
summary(e_10000)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      5017      1180528       415116       143    7547
## 2      3726      1052061       603612       131    5746
## 3      4532      1100721       259510       132    6674
## 4      4150      1283115       829728        99    6646
## 5      4430      1139185       545958       146    6748

print(object.size(e_10000), units = "MB")
## 863.9 Mb

t_e_10000b <- system.time(
    e_10000b <- oncoSimulPop(5,
                             u,
                             model = "Exp",
                             mu = 1e-7,
                             detectionSize = 1e6,
                             detectionDrivers = NA,
                             detectionProb = NA,
                             keepPhylog = TRUE,
                             onlyCancer = FALSE,
                             keepEvery = NA,
                             mutationPropGrowth = TRUE,
                             mc.cores = 1
                             ))

t_e_10000b
##    user  system elapsed 
##   5.484   0.100   5.585
  
summary(e_10000b)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      2465      1305094       727989        91    6447
## 2      2362      1070225       400329       204    8345
## 3      2530      1121164       436721       135    8697
## 4      2593      1206293       664494       125    8149
## 5      2655      1186994       327835       191    8572

print(object.size(e_10000b), units = "MB")
## 488.3 Mb

ng <- 50000
u <- allFitnessEffects(noIntGenes = c(rep(0.1, ng/2), 
                                      rep(-0.1, ng/2)))
t_e_50000 <- system.time(
    e_50000 <- oncoSimulPop(5,
                            u,
                            model = "Exp",
                            mu = 1e-7,
                            detectionSize = 1e6,
                            detectionDrivers = NA,
                            detectionProb = NA,
                            keepPhylog = TRUE,
                            onlyCancer = FALSE,
                            keepEvery = NA,
                            mutationPropGrowth = FALSE,
                            mc.cores = 1
                            ))


t_e_50000
##    user  system elapsed
##  44.192   1.684  45.891 

summary(e_50000)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      7367      1009949       335455     75.00   18214
## 2      8123      1302324       488469     63.65   17379
## 3      8408      1127261       270690     72.57   21144
## 4      8274      1138513       318152     80.59   20994
## 5      7520      1073131       690814     70.00   18569

print(object.size(e_50000), units = "MB")
## 7598.6 Mb

ng <- 50000
u <- allFitnessEffects(noIntGenes = c(rep(0.1, ng/2), 
                                      rep(-0.1, ng/2)))
t_e_50000np <- system.time(
    e_50000np <- oncoSimulPop(5,
                              u,
                              model = "Exp",
                              mu = 1e-7,
                              detectionSize = 1e6,
                              detectionDrivers = NA,
                              detectionProb = NA,
                              keepPhylog = TRUE,
                              onlyCancer = FALSE,
                              keepEvery = 1,
                              mutationPropGrowth = FALSE,
                              mc.cores = 1
                              ))

t_e_50000np
##   user  system elapsed
## 42.316   2.764  45.079 

summary(e_50000np)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1     13406      1027949       410074     71.97   19469
## 2     12469      1071325       291852     66.00   17834
## 3     11821      1089834       245720     90.00   16711
## 4     14008      1165168       505607     77.61   19675
## 5     14759      1074621       205954     87.68   20597

print(object.size(e_50000np), units = "MB")
## 12748.4 Mb

ng <- 50000
u <- allFitnessEffects(noIntGenes = c(rep(0.1, ng/2), 
                                      rep(-0.1, ng/2)))

t_e_50000c <- system.time(
    e_50000c <- oncoSimulPop(5,
                             u,
                             model = "Exp",
                             mu = 1e-7,
                             detectionSize = 1e6,
                             detectionDrivers = NA,
                             detectionProb = NA,
                             keepPhylog = TRUE,
                             onlyCancer = FALSE,
                             keepEvery = NA,
                             mutationPropGrowth = TRUE,
                             mc.cores = 1
                             ))

t_e_50000c
##    user  system elapsed 
## 84.228   2.416  86.665

summary(e_50000c)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1     11178      1241970       344479     84.74   27137
## 2     12820      1307086       203544     91.94   33448
## 3     10592      1126091       161057     83.81   26064
## 4     11883      1351114       148986     65.68   25396
## 5     10518      1101392       253523     99.79   26082

print(object.size(e_50000c), units = "MB")
## 10904.9 Mb

r1 <- oncoSimulIndiv(u,
                     model = "Exp",
                     mu = 1e-7,
                     detectionSize = 1e6,
                     detectionDrivers = NA,
                     detectionProb = NA,
                     keepPhylog = TRUE,
                     onlyCancer = FALSE,
                     mutationPropGrowth = TRUE
                     )
summary(r1)[c(1, 8)]
##   NumClones  FinalTime
## 1      3887        345

print(object.size(r1), units = "MB")
## 160 Mb

## Size of the two largest objects inside:
sizes <- lapply(r1, function(x) object.size(x)/(1024^2))
sort(unlist(sizes), decreasing = TRUE)[1:2]
## Genotypes pops.by.time 
##       148.28        10.26 

dim(r1$Genotypes)
## [1] 10000  3887

2.3.2 McFarland model with 50,000 genes; the effect of keepEvery

We show an example of McFarland’s model with 50,000 genes in section 6.5.3. We will show here a few more examples with those many genes but with a different fitness specification and changing several other settings.

ng <- 50000
u <- allFitnessEffects(noIntGenes = c(rep(0.1, ng/2), 
                                      rep(-0.1, ng/2)))

t_mc_50000_nmpg <- system.time(
    mc_50000_nmpg <- oncoSimulPop(5,
                                  u,
                                  model = "McFL",
                                  mu = 1e-7,
                                  detectionSize = 1e6,
                                  detectionDrivers = NA,
                                  detectionProb = NA,
                                  keepPhylog = TRUE,
                                  onlyCancer = FALSE,
                                  keepEvery = NA,
                                  mutationPropGrowth = FALSE,
                                  mc.cores = 1
                                  ))
t_mc_50000_nmpg
##   user  system elapsed 
##  30.46    0.54   31.01 


summary(mc_50000_nmpg)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      1902      1002528       582752     284.2   31137
## 2      2159      1002679       404858     274.8   36905
## 3      2247      1002722       185678     334.5   42429
## 4      2038      1009606       493574     218.4   32519
## 5      2222      1004661       162628     291.0   38470

print(object.size(mc_50000_nmpg), units = "MB")
## 2057.6 Mb

t_mc_50000_nmpg_k <- system.time(
    mc_50000_nmpg_k <- oncoSimulPop(5,
                                    u,
                                    model = "McFL",
                                    mu = 1e-7,
                                    detectionSize = 1e6,
                                    detectionDrivers = NA,
                                    detectionProb = NA,
                                    keepPhylog = TRUE,
                                    onlyCancer = FALSE,
                                    keepEvery = 1,
                                    mutationPropGrowth = FALSE,
                                    mc.cores = 1
                                    ))

t_mc_50000_nmpg_k
##    user  system elapsed 
##  30.000   1.712  31.714 

summary(mc_50000_nmpg_k)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      8779      1000223       136453     306.7   38102
## 2      7442      1006563       428150     345.3   35139
## 3      8710      1003509       224543     252.3   35659
## 4      8554      1002537       103889     273.7   36783
## 5      8233      1003171       263005     301.8   35236

print(object.size(mc_50000_nmpg_k), units = "MB")
## 8101.4 Mb

ng <- 50000
u <- allFitnessEffects(noIntGenes = c(rep(0.1, ng/2), 
                                      rep(-0.1, ng/2)))

t_mc_50000_nmpg_3e6 <- system.time(
    mc_50000_nmpg_3e6 <- oncoSimulPop(5,
                                      u,
                                      model = "McFL",
                                      mu = 1e-7,
                                      detectionSize = 3e6,
                                      detectionDrivers = NA,
                                      detectionProb = NA,
                                      keepPhylog = TRUE,
                                      onlyCancer = FALSE,
                                      keepEvery = NA,
                                      mutationPropGrowth = FALSE,
                                      mc.cores = 1
                                      ))
t_mc_50000_nmpg_3e6
##    user  system elapsed 
##  77.240   1.064  78.308 

summary(mc_50000_nmpg_3e6)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      5487      3019083       836793     304.5   65121
## 2      4812      3011816       789146     286.3   53087
## 3      4463      3016896      1970957     236.6   45918
## 4      5045      3028142       956026     360.3   63464
## 5      4791      3029720       916692     358.1   55012

print(object.size(mc_50000_nmpg_3e6), units = "MB")
## 4759.3 Mb
                                                   

t_mc_50000_nmpg_5mu <- system.time(
    mc_50000_nmpg_5mu <- oncoSimulPop(5,
                                      u,
                                      model = "McFL",
                                      mu = 5e-7,
                                      detectionSize = 1e6,
                                      detectionDrivers = NA,
                                      detectionProb = NA,
                                      keepPhylog = TRUE,
                                      onlyCancer = FALSE,
                                      keepEvery = NA,
                                      mutationPropGrowth = FALSE,
                                      mc.cores = 1
                                      ))

t_mc_50000_nmpg_5mu
##    user  system elapsed 
## 167.332   1.796 169.167 

summary(mc_50000_nmpg_5mu)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      7963      1004415       408352     99.03   57548
## 2      8905      1010751       120155    130.30   74738
## 3      8194      1005465       274661     96.98   58546
## 4      9053      1014049       119943    112.23   75379
## 5      8982      1011817        95047     99.95   76757

print(object.size(mc_50000_nmpg_5mu), units = "MB")
## 8314.4 Mb

t_mc_50000_nmpg_5mu_k <- system.time(
    mc_50000_nmpg_5mu_k <- oncoSimulPop(5,
                                        u,
                                        model = "McFL",
                                        mu = 5e-7,
                                        detectionSize = 1e6,
                                        detectionDrivers = NA,
                                        detectionProb = NA,
                                        keepPhylog = TRUE,
                                        onlyCancer = FALSE,
                                        keepEvery = 1,
                                        mutationPropGrowth = FALSE,
                                        mc.cores = 1
                                        ))
                                                   
t_mc_50000_nmpg_5mu_k
##    user  system elapsed 
## 174.404   5.068 179.481 

summary(mc_50000_nmpg_5mu_k)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1     25294      1001597       102766     123.4   74524
## 2     23766      1006679       223010     124.3   71808
## 3     21755      1001379       203638     114.8   62609
## 4     24889      1012103       161003     119.3   75031
## 5     21844      1002927       255388     108.8   64556

print(object.size(mc_50000_nmpg_5mu_k), units = "MB")
## 22645.8 Mb
 

t_mc_50000 <- system.time(
    mc_50000 <- oncoSimulPop(5,
                             u,
                             model = "McFL",
                             mu = 1e-7,
                             detectionSize = 1e6,
                             detectionDrivers = NA,
                             detectionProb = NA,
                             keepPhylog = TRUE,
                             onlyCancer = FALSE,
                             keepEvery = NA,
                             mutationPropGrowth = TRUE,
                             mc.cores = 1
                             ))

t_mc_50000
##    user  system elapsed 
## 303.352   2.808 306.223 

summary(mc_50000)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1     13928      1010815       219814     210.9   91255
## 2     12243      1003267       214189     178.1   67673
## 3     13880      1014131       124354     161.4   88322
## 4     14104      1012941        75521     205.7   98583
## 5     12428      1005594       232603     167.4   70359

print(object.size(mc_50000), units = "MB")
## 12816.6 Mb

2.3.3 Examples with \(s = 0.05\)

A script with the above runs but using \(s=0.05\) instead of \(s=0.1\) is available from the repository (‘miscell-files/vignette_bench_Rout/large_num_genes_0.05.Rout’). I will single out a couple of cases here.

First, we repeat the run shown in section 2.3.2.5:

t_mc_50000_nmpg_5mu_k <- system.time(
    mc_50000_nmpg_5mu_k <- oncoSimulPop(2,
                                        u,
                                        model = "McFL",
                                        mu = 5e-7,
                                        detectionSize = 1e6,
                                        detectionDrivers = NA,
                                        detectionProb = NA,
                                        keepPhylog = TRUE,
                                        onlyCancer = FALSE,
                                        keepEvery = 1,
                                        mutationPropGrowth = FALSE,
                                        mc.cores = 1
                                        ))
t_mc_50000_nmpg_5mu_k
##    user  system elapsed 
## 305.512   5.164 310.711 
 
summary(mc_50000_nmpg_5mu_k)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1     61737      1003273       104460  295.8731  204214
## 2     65072      1000540       133068  296.6243  210231
 
print(object.size(mc_50000_nmpg_5mu_k), units = "MB")
## 24663.6 Mb

Note we use only two replicates, since those two already lead to a 24 GB returned object as we are tracking more than 60,000 clones, more than twice those with \(s=0.1\). The reason for the difference in number of clones and iterations is of course the change from \(s=0.1\) to \(s=0.05\): under the McFarland model to reach population sizes of \(10^6\) starting from an equilibrium population of 500 we need about 43 mutations (whereas only about 22 are needed if \(s=0.1\)⁵).

Next, let us rerun 2.3.2.1:

t_mc_50000_nmpg <- system.time(
    mc_50000_nmpg <- oncoSimulPop(5,
                                  u,
                                  model = "McFL",
                                  mu = 1e-7,
                                  detectionSize = 1e6,
                                  detectionDrivers = NA,
                                  detectionProb = NA,
                                  keepPhylog = TRUE,
                                  onlyCancer = FALSE,
                                  keepEvery = NA,
                                  mutationPropGrowth = FALSE,
                                  mc.cores = 1
                                  ))
t_mc_50000_nmpg
##    user  system elapsed 
## 111.236   0.596 111.834 
 
summary(mc_50000_nmpg)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      2646      1000700       217188   734.475  108566
## 2      2581      1001626       209873   806.500  107296
## 3      2903      1001409       125148   841.700  120859
## 4      2310      1000146       473948   906.300   91519
## 5      2704      1001290       448409   838.800  103556
 
print(object.size(mc_50000_nmpg), units = "MB")
## 2638.3 Mb

Using \(s=0.05\) leads to a large increase in final time and number of iterations. However, as we are using the keepEvery = NA setting, the increase in number of clones tracked and in size of returned object is relatively small.

2.3.4 The different consequences of keepEvery = NA in the Exp and McFL models

We have seen that keepEvery = NA often leads to much smaller returned objects when using the McFarland model than when using the Exp model. Why? Because in the McFarland model there is strong competition and there can be complete clonal sweeps so that in extreme cases a single clone might be all that is left after some time. This is not the case in the exponential models.

Of course, the details depend on the difference in fitness effects between different genotypes (or clones). In particular, we have seen several examples where even with keepEvery=NA there are a lot of clones in the McFL models. In those examples many clones had identical fitness (the fitness effects of all genes with positive fitness was the same, and ditto for the genes with negative fitness effects), so no clone ends up displacing all the others.

2.3.5 Are we keeping the complete history (genealogy) of the clones?

Yes we are if we run with keepPhylog = TRUE, regardless of the setting for keepEvery. As explained in section 12.2, OncoSimulR prunes clones that never had a population size larger than zero at any sampling period (so they are not reflected in the pops.by.time matrix in the output). And when we set keepEvery = NA we are telling OncoSimulR to discard all sampling periods except the very last one (i.e., the pops.by.time matrix contains only the clones with 1 or more cells at the end of the simulation).

keepPhylog operates differently: it records the exact time at which a clone appeared and the clone that gave rise to it. This information is kept regardless of whether or not those clones appear in the pops.by.time matrix.

Keeping the complete genealogy might be of limited use if the pops.by.time matrix only contains the very last period. However, you can use plotClonePhylog and ask to be shown only clones that exist in the very last period (while of course showing all of their ancestors, even if those are now extinct —i.e., regardless of their abundance).

For instance, in run 2.3.1.3 we could have looked at the information stored about the genealogy of clones by doing (we look at the first “individual” of the simulation, of the five “individuals” we simulated):

head(e_50000[[1]]$other$PhylogDF)
##   parent child   time
## 1         3679 0.8402
## 2         4754 1.1815
## 3        20617 1.4543
## 4        15482 2.3064
## 5         4431 3.7130
## 6        41915 4.0628

tail(e_50000[[1]]$other$PhylogDF)
##                          parent                            child time
## 20672               3679, 20282               3679, 20282, 22359 75.0
## 20673        3679, 17922, 22346        3679, 17922, 22346, 35811 75.0
## 20674                2142, 3679                2142, 3679, 25838 75.0
## 20675        3679, 17922, 19561        3679, 17922, 19561, 43777 75.0
## 20676 3679, 15928, 19190, 20282 3679, 15928, 19190, 20282, 49686 75.0
## 20677         2142, 3679, 16275         2142, 3679, 16275, 24201 75.0

where each row corresponds to one event of appearance of a new clone, the column labeled “parent” are the mutated genes in the parent, and the column labeled “child” are the mutated genes in the child.

And we could plot the genealogical relationships of clones that have a population size of at least one in the last period (again, while of course showing all of their ancestors, even if those are now extinct —i.e., regardless of their current numbers) doing:

plotClonePhylog(e_50000[[1]]) ## plot not shown

What is the cost of keep the clone genealogies? In terms of time it is minor. In terms of space, and as shown in the example above, we can end up storing a data frame with tends of thousands of rows and three columns (two factors, one float). In the example above the size of that data frame is approximately 2 MB for a single simulation. This is much smaller than the pops.by.time or Genotypes matrices, but it can quickly build up if you routinely launch, say, 1000 simulations via oncoSimulPop. That is why the default is keepPhylog = FALSE as this information is not needed as often as that in the other two matrices (pops.by.time and Genotypes).

3.1.1 Explicit mapping of genotypes to fitness

Conceptually, the simplest way to specify fitness is to specify the mapping of all genotypes to fitness explicitly. An example will make this clear. Let’s suppose you have a simple two-gene scenario, so a total of four genotypes, and you have a data frame with genotypes and fitness, where genoytpes are specified as character vectors, with mutated genes separated by commas:

m4 <- data.frame(G = c("WT", "A", "B", "A, B"), F = c(1, 2, 3, 4))

Now, let’s give that to the allFitnessEffects function:

fem4 <- allFitnessEffects(genotFitness = m4)
## Warning in to_genotFitness_std(genotFitness, simplify = TRUE): First
## column of genotype fitness is a factor. Converting to character.
## Column names of object not Genotype and Fitness. Renaming them assuming that is what you wanted

(The message is just telling you what the program guessed you wanted.)

That’s it. You can try to plot that fitnessEffects object

try(plot(fem4))

In this case, you probably want to plot the fitness landscape.

plotFitnessLandscape(evalAllGenotypes(fem4))

You can also check what OncoSimulR thinks the fitnesses are, with the evalAllGenotypes function that we will use repeatedly below (of course, here we should see the same fitnesses we entered):

evalAllGenotypes(fem4, addwt = TRUE)
##   Genotype Fitness
## 1       WT       1
## 2        A       2
## 3        B       3
## 4     A, B       4

And you can plot the fitness landscape:

plotFitnessLandscape(evalAllGenotypes(fem4))

To specify the mapping you can also use a matrix (or data frame) with \(g + 1\) columns; each of the first \(g\) columns contains a 1 or a 0 indicating that the gene of that column is mutated or not. Column \(g+ 1\) contains the fitness values. And you do not even need to specify all the genotypes (we will assume that the missing genotypes have fitness 1):

m6 <- cbind(c(1, 1), c(1, 0), c(2, 3))
fem6 <- allFitnessEffects(genotFitness = m6)
## No column names: assigning gene names from LETTERS
## Warning in to_genotFitness_std(genotFitness, simplify = TRUE): No
## wildtype in the fitness landscape!!! Adding it with fitness 1.
evalAllGenotypes(fem6, addwt = TRUE)
##   Genotype Fitness
## 1       WT       1
## 2        A       3
## 3        B       0
## 4     A, B       2
## plot(fem6)

plotFitnessLandscape(evalAllGenotypes(fem6))

This way of giving a fitness specification to OncoSimulR might be ideal if you directly generate random mappings of genotypes to fitness (or random fitness landscapes), as we will do in section 9. Specially when the fitness landscape contains many non-viable genotypes (which are considered those with fitness —birth rate— \(<1e-9\)) this can result in considerable savings as we only need to check the fitness of the viable genotypes in a table (a C++ map). Note, however, that using the Bozic model with the fitness landscape specification is not tested. In addition, for speed, missing genotypes from the fitness landscape specification are taken to be non-viable genotypes (beware!! this is a breaking change relative to versions < 2.9.1)⁶.

3.1.2 How to specify fitness effects with the lego system

An alternative general approach followed in many genetic simulators is to specify how particular combinations of alleles modify the wildtype genotype or the genotype that contains the individual effects of the interacting genes (e.g., see equation 1 in the supplementary material for FFPopSim (Zanini and Neher 2012)). For example, if we specify that a mutation in “A” contributes 0.04, a mutation in “B” contributes 0.03, and the double mutation “A:B” contributes 0.1, that means that the fitness of the “A, B” genotype (the genotype with A and B mutated) is that of the wildtype (1, by default), plus (actually, times —see section 3.2— but plus on the log scale) the effects of having A mutated, plus (times) the effects of having B mutated, plus (times) the effects of “A:B” both being mutated.

We will see below that with the “lego system” it is possible to do something very similar to the explicit mapping of section 3.1.1. But this will sometimes require a more cumbersome notation (and sometimes also will require your doing some math). We will see examples in sections 3.7.1, 3.7.2 and 3.7.3 or the example in 5.4.2. But then, if we can be explicit about (at least some of) the mappings \(genotype \rightarrow fitness\), how are these procedures different? When you use the “lego system” you can combine both a partial explicit mapping of genotypes to fitness with arbitrary fitness effects of other genes/modules. In other words, with the “lego system” OncoSimulR makes it simple to be explicit about the mapping of specific genotypes, while also using the “how this specific effects modifies previous effects” logic, leading to a flexible specification. This also means that in many cases the same fitness effects can be specified in several different ways.

Most of the rest of this section is devoted to explaining how to combine those pieces. Before that, however, we need to discuss the fitness model we use.

3.2.1 McFarland parameterization

In the original model of McFarland et al. (2013), the effects of drivers contribute to the numerator of the birth rate, and those of the (deleterious) passengers to the denominator as: \(\frac{(1 + s)^d}{(1 + s_p)^p}\), where \(d\) and \(p\) are, respectively, the total number of drivers and passengers in a genotype, and here the fitness effects of all drivers is the same (\(s\)) and that of all passengers the same too (\(s_p\)). Note that, as written above, and as explicitly said in McFarland et al. (2013) (see p. 2911) and CD McFarland (2014) (see p. 9), “(…) \(s_p\) is the fitness disadvantage conferred by a passenger”. In other words, the larger the \(s_p\) the more deleterious the passenger.

This is obvious, but I make it explicit because in our parameterization a positive \(s\) means fitness advantage, whereas fitness disadvantages are associated with negative \(s\). Of course, if you rewrite the above expression as \(\frac{(1 + s)^d}{(1 - s_p)^p}\) then we are back to the “positive means fitness advantage and negative means fitness disadvantage”.

As CD McFarland (2014) explains (see p. 9, bottom), we can rewrite the above expression so that there are no terms in the denominator. McFarland writes it as (I copy verbatim from the fourth and fifth lines from the bottom on his p. 9) \((1 + s_d)^{n_d} (1 - s_p^{'})^{n_p}\) where \(s_p^{'} = s_p/(1 + s_p)\).

However, if we want to express everything as products (no ratios) and use the “positive s means advantage and negative s means disadvantage” rule, we want to write the above expression as \((1 + s_d)^{n_d} (1 + s_{pp})^{n_p}\) where \(s_{pp} = -s_p/(1 + s_p)\). And this is actually what we do in v.2. There is an example, for instance, in section 6.5.2 where you will see:

sp <- 1e-3 
spp <- -sp/(1 + sp)

so we are going from the “(…) \(s_p\) is the fitness disadvantage conferred by a passenger” in McFarland et al. (2013) (p. 2911) and CD McFarland (2014) (p. 9) to the expression where we have a product \(\prod (1 + s_i)\), with the “positive s means advantage and negative s means disadvantage” rule. This reparameterization applies to v.2. In v.1 we use the same parameterization as in the original one in McFarland et al. (2013), but with the “positive s means advantage and negative s means disadvantage” rule (so we are using expression \(\frac{(1 + s)^d}{(1 - s_p)^p}\)).

For death rate, we use the expression that McFarland et al. (2013) (see their p. 2911) use “(…) for large cancers (grown to \(10^6\) cells)”: \(D(N) = \log(1 + N/K)\) where \(K\) is the initial equilibrium population size. As the authors explain, for large N/K the above expression “(…) recapitulates Gompertzian dynamics observed experimentally for large tumors”.

By default, OncoSimulR uses a value of \(K=initSize/(e^{1} - 1)\) so that the starting population is at equilibrium.

3.2.2 No viability of clones and types of models

For all models where fitness affects directly the birth rate (all except Bozic), if you specify that some event (say, mutating gene A) has \(s_A \le -1\), if that event happens then birth rate becomes zero. This is taken to indicate that the clone is not even viable and thus disappears immediately without any chance for mutation⁸.

Models based on Bozic, however, have a birth rate of 1 and mutations affect the death rate. In this case, a death rate larger than birth rate, per se, does not signal immediate extinction and, moreover, even for death rates that are a few times larger than birth rates, the clone could mutate before becoming extinct⁹.

In general, if you want to identify some mutations or some combinations of mutations as leading to immediate extinction (i.e., no viability), of the affected clone, set it to \(-\infty\) as this would work even if how birth rates of 0 are handled changes. Most examples below evaluate fitness by its effects on the birth rate. You can see one where we do it both ways in Section 3.11.1.

3.4.1 AND, OR, XOR relationships

The literature on Oncogenetic trees, CBNs, etc, has used graphs as a way of showing the restrictions in the order in which mutations can accumulate. The meaning of “convergent arrows” in these graphs, however, differs. In Figure 1 of Korsunsky et al. (2014) we are shown a simple diagram that illustrates the three basic different meanings of convergent arrows using two parental nodes. We will illustrate it here with three. Suppose we focus on node “g” in the following figure (we will create it shortly)

data(examplesFitnessEffects)
plot(examplesFitnessEffects[["cbn1"]])

• In relationships of the type used in Conjunctive Bayesian Networks (CBN) (e.g., Gerstung et al. 2009), we are modeling an AND relationship, also called CMPN by Korsunsky et al. (2014) or monotone relationship by Farahani and Lagergren (2013). If the relationship in the graph is fully respected, then “g” will only appear if all of “c”, “d”, and “e” are already mutated.

• Semimonotone relationships sensu Farahani and Lagergren (2013) or DMPN sensu Korsunsky et al. (2014) are OR relationships: “g” will appear if one or more of “c”, “d”, or “e” are already mutated.

• XMPN relationships (Korsunsky et al. 2014) are XOR relationships: “g” will be present only if exactly one of “c”, “d”, or “e” is present.

Note that Oncogenetic trees (Desper et al. 1999; Szabo and Boucher 2008) need not deal with the above distinctions, since the DAGs are trees: no node has more than one incoming connection or more than one parent¹⁰.

To have a flexible way of specifying all of these restrictions, we will want to be able to say what kind of dependency each child node has on its parents.

3.4.2 Fitness effects

Those DAGs specify dependencies and, as explained in Diaz-Uriarte (2015), it is simple to map them to a simple evolutionary model: any set of mutations that does not conform to the restrictions encoded in the graph will have a fitness of 0. However, we might not want to require absolute compliance with the DAG. This means we might want to allow deviations from the DAG with a corresponding penalization that is, however, not identical to setting fitness to 0 (again, see Diaz-Uriarte 2015). This we can do by being explicit about the fitness effects of the deviations from the restrictions encoded in the DAG. We will use below a column of s for the fitness effect when the restrictions are satisfied and a column of sh when they are not. (See also 3.2 for the details about the meaning of the fitness effects).

That way of specifying fitness effects makes it also trivial to use the model in Hjelm, Höglund, and Lagergren (2006) where all mutations might be allowed to occur, but the presence of some mutations increases the probability of occurrence of other mutations. For example, the values of sh could be all small positive ones (or for mildly deleterious effects, small negative numbers), while the values of s are much larger positive numbers.

3.4.3 Extended posets

In version 1 of this package we used posets in the sense of Beerenwinkel, Eriksson, and Sturmfels (2007) and Gerstung et al. (2009), as explained in section 13.1 and in the help for poset. Here, we continue using two columns, that specify parents and children, but we add columns for the specific values of fitness effects (both s and sh —i.e., fitness effects for what happens when restrictions are and are not satisfied) and for the type of dependency as explained in section 3.4.1.

We can now illustrate the specification of different fitness effects using DAGs.

3.4.4 DAGs: A first conjunction (AND) example

cs <-  data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"),
                 child = c("a", "b", "d", "e", "c", "c", rep("g", 3)),
                 s = 0.1,
                 sh = -0.9,
                 typeDep = "MN")

cbn1 <- allFitnessEffects(cs)

(We skip one letter, just to show that names need not be consecutive or have any particular order.)

We can get a graphical representation using the default “graphNEL”

plot(cbn1)

or one using “igraph”:

plot(cbn1, "igraph")

Since we have a parent and children, the reingold.tilford layout is probably the best here, so you might want to use that:

library(igraph) ## to make the reingold.tilford layout available
plot(cbn1, "igraph", layout = layout.reingold.tilford)

And what is the fitness of all genotypes?

gfs <- evalAllGenotypes(cbn1, order = FALSE, addwt = TRUE)

gfs[1:15, ]
##    Genotype Fitness
## 1        WT    1.00
## 2         a    1.10
## 3         b    1.10
## 4         c    0.10
## 5         d    1.10
## 6         e    1.10
## 7         g    0.10
## 8      a, b    1.21
## 9      a, c    0.11
## 10     a, d    1.21
## 11     a, e    1.21
## 12     a, g    0.11
## 13     b, c    0.11
## 14     b, d    1.21
## 15     b, e    1.21

You can verify that for each genotype, if a mutation is present without all of its dependencies present, you get a \((1 - 0.9)\) multiplier, and you get a \((1 + 0.1)\) multiplier for all the rest with its direct parents satisfied. For example, genotypes “a”, or “b”, or “d”, or “e” have fitness \((1 + 0.1)\), genotype “a, b, c” has fitness \((1 + 0.1)^3\), but genotype “a, c” has fitness \((1 + 0.1) (1 - 0.9) = 0.11\).

3.4.5 DAGs: A second conjunction example

Let’s try a first attempt at a somewhat more complex example, where the fitness consequences of different genes differ.

c1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"),
                 child = c("a", "b", "d", "e", "c", "c", rep("g", 3)),
                 s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)),
                 sh = c(rep(0, 4), c(-.1, -.2), c(-.05, -.06, -.07)),
                 typeDep = "MN")

try(fc1 <- allFitnessEffects(c1))

If you try this, you’ll get an error. There is an error because the “sh” varies within a child, and we do not allow that for a poset-type specification, as it is ambiguous. If you need arbitrary fitness values for arbitrary combinations of genotypes, you can specify them using epistatic effects as in section 3.7 and order effects as in section 3.6.

Why do we need to specify as many “s” and “sh” as there are rows (or a single one, that gets expanded to those many) when the “s” and “sh” are properties of the child node, not of the edges? Because, for ease, we use a data.frame.

We fix the error in our specification. Notice that the “sh” is not set to \(-1\) in these examples. If you want strict compliance with the poset restrictions, you should set \(sh = -1\) or, better yet, \(sh = -\infty\) (see section 3.2.2), but having an \(sh > -1\) will lead to fitnesses that are \(> 0\) and, thus, is a way of modeling small deviations from the poset (see discussion in Diaz-Uriarte 2015).

Note that for those nodes that depend only on “Root” the type of dependency is irrelevant.

c1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"),
                 child = c("a", "b", "d", "e", "c", "c", rep("g", 3)),
                 s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)),
                 sh = c(rep(0, 4), c(-.9, -.9), rep(-.95, 3)),
                 typeDep = "MN")

cbn2 <- allFitnessEffects(c1)

We could get graphical representations but the figures would be the same as in the example in section 3.4.4, since the structure has not changed, only the numeric values.

What is the fitness of all possible genotypes? Here, order of events per se does not matter, beyond that considered in the poset. In other words, the fitness of genotype “a, b, c” is the same no matter how we got to “a, b, c”. What matters is whether or not the genes on which each of “a”, “b”, and “c” depend are present or not (I only show the first 10 genotypes)

gcbn2 <- evalAllGenotypes(cbn2, order = FALSE)
gcbn2[1:10, ]
##    Genotype Fitness
## 1         a  1.0100
## 2         b  1.0200
## 3         c  0.1000
## 4         d  1.0300
## 5         e  1.0400
## 6         g  0.0500
## 7      a, b  1.0302
## 8      a, c  0.1010
## 9      a, d  1.0403
## 10     a, e  1.0504

Of course, if we were to look at genotypes but taking into account order of occurrence of mutations, we would see no differences

gcbn2o <- evalAllGenotypes(cbn2, order = TRUE, max = 1956)
gcbn2o[1:10, ]
##    Genotype Fitness
## 1         a  1.0100
## 2         b  1.0200
## 3         c  0.1000
## 4         d  1.0300
## 5         e  1.0400
## 6         g  0.0500
## 7     a > b  1.0302
## 8     a > c  0.1010
## 9     a > d  1.0403
## 10    a > e  1.0504

(The \(max = 1956\) is there so that we show all the genotypes, even if they are more than 256, the default.)

You can check the output and verify things are as they should. For instance:

all.equal(
        gcbn2[c(1:21, 22, 28, 41, 44, 56, 63 ) , "Fitness"],
        c(1.01, 1.02, 0.1, 1.03, 1.04, 0.05,
          1.01 * c(1.02, 0.1, 1.03, 1.04, 0.05),
          1.02 * c(0.10, 1.03, 1.04, 0.05),
          0.1 * c(1.03, 1.04, 0.05),
          1.03 * c(1.04, 0.05),
          1.04 * 0.05,
          1.01 * 1.02 * 1.1,
          1.01 * 0.1 * 0.05,
          1.03 * 1.04 * 0.05,
          1.01 * 1.02 * 1.1 * 0.05,
          1.03 * 1.04 * 1.2 * 0.1, ## notice this
          1.01 * 1.02 * 1.03 * 1.04 * 1.1 * 1.2
          ))
## [1] TRUE

A particular one that is important to understand is genotype with mutated genes “c, d, e, g”:

gcbn2[56, ] 
##      Genotype  Fitness
## 56 c, d, e, g 0.128544
all.equal(gcbn2[56, "Fitness"], 1.03 * 1.04 * 1.2 * 0.10)
## [1] TRUE

where “g” is taken as if its dependencies are satisfied (as “c”, ``d“, and”e" are present) even when the dependencies of “c” are not satisfied (and that is why the term for “c” is 0.9).

3.4.6 DAGs: A semimonotone or “OR” example

We will reuse the above example, changing the type of relationship:

s1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"),
                 child = c("a", "b", "d", "e", "c", "c", rep("g", 3)),
                 s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)),
                 sh = c(rep(0, 4), c(-.9, -.9), rep(-.95, 3)),
                 typeDep = "SM")

smn1 <- allFitnessEffects(s1)

It looks like this (where edges are shown in blue to denote the semimonotone relationship):

plot(smn1)

gsmn1 <- evalAllGenotypes(smn1, order = FALSE)

Having just one parental dependency satisfied is now enough, in contrast to what happened before. For instance:

gcbn2[c(8, 12, 22), ]
##    Genotype Fitness
## 8      a, c 0.10100
## 12     b, c 0.10200
## 22  a, b, c 1.13322
gsmn1[c(8, 12, 22), ]
##    Genotype Fitness
## 8      a, c 1.11100
## 12     b, c 1.12200
## 22  a, b, c 1.13322

gcbn2[c(20:21, 28), ]
##    Genotype Fitness
## 20     d, g 0.05150
## 21     e, g 0.05200
## 28  a, c, g 0.00505
gsmn1[c(20:21, 28), ]
##    Genotype Fitness
## 20     d, g  1.2360
## 21     e, g  1.2480
## 28  a, c, g  1.3332

3.4.7 An “XMPN” or “XOR” example

Again, we reuse the example above, changing the type of relationship:

x1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"),
                 child = c("a", "b", "d", "e", "c", "c", rep("g", 3)),
                 s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)),
                 sh = c(rep(0, 4), c(-.9, -.9), rep(-.95, 3)),
                 typeDep = "XMPN")

xor1 <- allFitnessEffects(x1)

It looks like this (edges in red to denote the “XOR” relationship):

plot(xor1)

gxor1 <- evalAllGenotypes(xor1, order = FALSE)

Whenever “c” is present with both “a” and “b”, the fitness component for “c” will be \((1 - 0.1)\). Similarly for “g” (if more than one of ``d“,”e“, or”c" is present, it will show as \((1 - 0.05)\)). For example:

gxor1[c(22, 41), ] 
##    Genotype Fitness
## 22  a, b, c 0.10302
## 41  d, e, g 0.05356
c(1.01 * 1.02 * 0.1, 1.03 * 1.04 * 0.05)
## [1] 0.10302 0.05356

However, having just both “a” and “b” is identical to the case with CBN and the monotone relationship (see sections 3.4.5 and 3.4.6). If you want the joint presence of “a” and “b” to result in different fitness than the product of the individual terms, without considering the presence of “c”, you can specify that using general epistatic effects (section 3.7).

We also see a very different pattern compared to CBN (section 3.4.5) here:

gxor1[28, ] 
##    Genotype Fitness
## 28  a, c, g  1.3332
1.01 * 1.1 * 1.2
## [1] 1.3332

as exactly one of the dependencies for both “c” and “g” are satisfied.

But

gxor1[44, ] 
##      Genotype  Fitness
## 44 a, b, c, g 0.123624
1.01 * 1.02 * 0.1 * 1.2
## [1] 0.123624

is the result of a \(0.1\) for “c” (and a \(1.2\) for “g” that has exactly one of its dependencies satisfied).

3.4.8 Posets: the three types of relationships

p3 <- data.frame(
    parent = c(rep("Root", 4), "a", "b", "d", "e", "c", "f"),
    child = c("a", "b", "d", "e", "c", "c", "f", "f", "g", "g"),
    s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3),
    sh = c(rep(0, 4), c(-.9, -.9), c(-.95, -.95), c(-.99, -.99)),
    typeDep = c(rep("--", 4), 
                "XMPN", "XMPN", "MN", "MN", "SM", "SM"))
fp3 <- allFitnessEffects(p3)

This is how it looks like:

plot(fp3)

We can also use “igraph”:

plot(fp3, "igraph", layout.reingold.tilford)

gfp3 <- evalAllGenotypes(fp3, order = FALSE)

Let’s look at a few:

gfp3[c(9, 24, 29, 59, 60, 66, 119, 120, 126, 127), ]
##                Genotype   Fitness
## 9                  a, c 1.1110000
## 24                 d, f 0.0515000
## 29              a, b, c 0.1030200
## 59              c, f, g 0.0065000
## 60              d, e, f 1.2854400
## 66           a, b, c, f 0.0051510
## 119       c, d, e, f, g 0.1671072
## 120    a, b, c, d, e, f 0.1324260
## 126    b, c, d, e, f, g 1.8749428
## 127 a, b, c, d, e, f, g 0.1721538

c(1.01 * 1.1, 1.03 * .05, 1.01 * 1.02 * 0.1, 0.1 * 0.05 * 1.3,
  1.03 * 1.04 * 1.2, 1.01 * 1.02 * 0.1 * 0.05,
  0.1 * 1.03 * 1.04 * 1.2 * 1.3,
  1.01 * 1.02 * 0.1 * 1.03 * 1.04 * 1.2,
  1.02 * 1.1 * 1.03 * 1.04 * 1.2 * 1.3,
  1.01 * 1.02 * 1.03 * 1.04 * 0.1 * 1.2 * 1.3)
##  [1] 1.1110000 0.0515000 0.1030200 0.0065000 1.2854400 0.0051510
##  [7] 0.1671072 0.1324260 1.8749428 0.1721538

As before, looking at the order of mutations makes no difference (look at the test directory to see a test that verifies this assertion).

3.5.1 What does a module provide

Modules can provide very compact ways of specifying relationships when you want to, well, model the existence of modules. For simplicity suppose there is a module, “A”, made of genes “a1” and “a2”, and a module “B”, made of a single gene “b1”. Module “B” can mutate if module “A” is mutated, but mutating both “a1” and “a2” provides no additional fitness advantage compared to mutating only a single one of them. We can specify this as:

s <- 0.2
sboth <- (1/(1 + s)) - 1
m0 <- allFitnessEffects(data.frame(
    parent = c("Root", "Root", "a1", "a2"),
    child = c("a1", "a2", "b", "b"),
    s = s,
    sh = -1,
    typeDep = "OR"),
                        epistasis = c("a1:a2" = sboth))
evalAllGenotypes(m0, order = FALSE, addwt = TRUE)
##    Genotype Fitness
## 1        WT    1.00
## 2        a1    1.20
## 3        a2    1.20
## 4         b    0.00
## 5    a1, a2    1.20
## 6     a1, b    1.44
## 7     a2, b    1.44
## 8 a1, a2, b    1.44

Note that we need to add an epistasis term, with value “sboth” to capture the idea of “mutating both”a1" and “a2” provides no additional fitness advantage compared to mutating only a single one of them“; see details in section 3.7.

Now, specify it using modules:

s <- 0.2
m1 <- allFitnessEffects(data.frame(
    parent = c("Root", "A"),
    child = c("A", "B"),
    s = s,
    sh = -1,
    typeDep = "OR"),
                        geneToModule = c("Root" = "Root",
                                         "A" = "a1, a2",
                                         "B" = "b1"))
evalAllGenotypes(m1, order = FALSE, addwt = TRUE)
##     Genotype Fitness
## 1         WT    1.00
## 2         a1    1.20
## 3         a2    1.20
## 4         b1    0.00
## 5     a1, a2    1.20
## 6     a1, b1    1.44
## 7     a2, b1    1.44
## 8 a1, a2, b1    1.44

This captures the ideas directly. The typing savings here are small, but they can be large with modules with many genes.

3.5.2 Specifying modules

How do you specify modules? The general procedure is simple: you pass a vector that makes explicit the mapping from modules to sets of genes. We just saw an example. There are several additional examples such as 3.5.3, 3.6.2, 3.7.4.

It is important to note that, once you specify modules, we expect all of the relationships (except those that involve the non interacting genes) to be specified as modules. Thus, all elements of the epistasis, posets (the DAGs) and order effects components should be specified in terms of modules. But you can, of course, specify a module as containing a single gene (and a single gene with the same name as the module).

What about the “Root” node? If you use a “restriction table”, that restriction table (that DAG) must have a node named “Root” and in the mapping of genes to module there must be a first entry that has a module and gene named “Root”, as we saw above with geneToModule = c("Root" = "Root", .... We force you to do this to be explicit about the “Root” node. This is not needed (thought it does not hurt) with other fitness specifications. For instance, if we have a model with two modules, one of them with two genes (see details in section 3.8) we do not need to pass a “Root” as in

fnme <- allFitnessEffects(epistasis = c("A" = 0.1,
                                        "B" = 0.2),
                          geneToModule = c("A" = "a1, a2",
                                           "B" = "b1"))
evalAllGenotypes(fnme, order = FALSE, addwt = TRUE)
##     Genotype Fitness
## 1         WT    1.00
## 2         a1    1.10
## 3         a2    1.10
## 4         b1    1.20
## 5     a1, a2    1.10
## 6     a1, b1    1.32
## 7     a2, b1    1.32
## 8 a1, a2, b1    1.32

but it is also OK to have a “Root” in the geneToModule:

fnme2 <- allFitnessEffects(epistasis = c("A" = 0.1,
                                        "B" = 0.2),
                          geneToModule = c(
                              "Root" = "Root",
                              "A" = "a1, a2",
                              "B" = "b1"))
evalAllGenotypes(fnme, order = FALSE, addwt = TRUE)
##     Genotype Fitness
## 1         WT    1.00
## 2         a1    1.10
## 3         a2    1.10
## 4         b1    1.20
## 5     a1, a2    1.10
## 6     a1, b1    1.32
## 7     a2, b1    1.32
## 8 a1, a2, b1    1.32

3.5.3 Modules and posets again: the three types of relationships and modules

We use the same specification of poset, but add modules. To keep it manageable, we only add a few genes for some modules, and have some modules with a single gene. Beware that the number of genotypes is starting to grow quite fast, though. We capitalize to differentiate modules (capital letters) from genes (lowercase with a number), but this is not needed.

p4 <- data.frame(
    parent = c(rep("Root", 4), "A", "B", "D", "E", "C", "F"),
    child = c("A", "B", "D", "E", "C", "C", "F", "F", "G", "G"),
    s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3),
    sh = c(rep(0, 4), c(-.9, -.9), c(-.95, -.95), c(-.99, -.99)),
    typeDep = c(rep("--", 4), 
                "XMPN", "XMPN", "MN", "MN", "SM", "SM"))

fp4m <- allFitnessEffects(
    p4,
    geneToModule = c("Root" = "Root", "A" = "a1",
                     "B" = "b1, b2", "C" = "c1",
                     "D" = "d1, d2", "E" = "e1",
                     "F" = "f1, f2", "G" = "g1"))

By default, plotting shows the modules:

plot(fp4m)

but we can show the gene names instead of the module names:

plot(fp4m, expandModules = TRUE)

or

plot(fp4m, "igraph", layout = layout.reingold.tilford, 
     expandModules = TRUE)

We obtain the fitness of all genotypes in the usual way:

gfp4 <- evalAllGenotypes(fp4m, order = FALSE, max = 1024)

Let’s look at a few of those:

gfp4[c(12, 20, 21, 40, 41, 46, 50, 55, 64, 92,
       155, 157, 163, 372, 632, 828), ]
##                   Genotype  Fitness
## 12                  a1, b2 1.030200
## 20                  b1, b2 1.020000
## 21                  b1, c1 1.122000
## 40                  c1, g1 0.130000
## 41                  d1, d2 1.030000
## 46                  d2, e1 1.071200
## 50                  e1, f1 0.052000
## 55                  f2, g1 0.065000
## 64              a1, b2, c1 0.103020
## 92              b1, b2, c1 1.122000
## 155             c1, f2, g1 0.006500
## 157             d1, d2, f1 0.051500
## 163             d1, f1, f2 0.051500
## 372         d1, d2, e1, f2 1.285440
## 632     d1, d2, e1, f1, f2 1.285440
## 828 b2, c1, d1, e1, f2, g1 1.874943

c(1.01 * 1.02, 1.02, 1.02 * 1.1, 0.1 * 1.3, 1.03, 
  1.03 * 1.04, 1.04 * 0.05, 0.05 * 1.3,  
  1.01 * 1.02 * 0.1, 1.02 * 1.1, 0.1 * 0.05 * 1.3,
  1.03 * 0.05, 1.03 * 0.05, 1.03 * 1.04 * 1.2, 1.03 * 1.04 * 1.2, 
  1.02 * 1.1 * 1.03 * 1.04 * 1.2 * 1.3)
##  [1] 1.030200 1.020000 1.122000 0.130000 1.030000 1.071200 0.052000
##  [8] 0.065000 0.103020 1.122000 0.006500 0.051500 0.051500 1.285440
## [15] 1.285440 1.874943

3.6.1 Order effects: three-gene orders

Consider this case, where three specific three-gene orders and two two-gene orders (one of them a subset of one of the three) lead to different fitness compared to the wild-type. We add also modules, to show its usage (but just limit ourselves to using one gene per module here).

Order effects are specified using a \(x > y\), which means that that order effect is satisfied when module \(x\) is mutated before module \(y\).

o3 <- allFitnessEffects(orderEffects = c(
                            "F > D > M" = -0.3,
                            "D > F > M" = 0.4,
                            "D > M > F" = 0.2,
                            "D > M"     = 0.1,
                            "M > D"     = 0.5),
                        geneToModule =
                            c("M" = "m",
                              "F" = "f",
                              "D" = "d") )


(ag <- evalAllGenotypes(o3, addwt = TRUE, order = TRUE))
##     Genotype Fitness
## 1         WT    1.00
## 2          d    1.00
## 3          f    1.00
## 4          m    1.00
## 5      d > f    1.00
## 6      d > m    1.10
## 7      f > d    1.00
## 8      f > m    1.00
## 9      m > d    1.50
## 10     m > f    1.00
## 11 d > f > m    1.54
## 12 d > m > f    1.32
## 13 f > d > m    0.77
## 14 f > m > d    1.50
## 15 m > d > f    1.50
## 16 m > f > d    1.50

(The meaning of the notation in the output table is as follows: “WT” denotes the wild-type, or non-mutated clone. The notation \(x > y\) means that a mutation in “x” happened before a mutation in “y”. A genotype \(x > y\ \_\ z\) means that a mutation in “x” happened before a mutation in “y”; there is also a mutation in “z”, but that is a gene for which order does not matter).

The values for the first nine genotypes come directly from the fitness specifications. The 10th genotype matches \(D > F > M\) (\(= (1 + 0.4)\)) but also \(D > M\) (\((1 + 0.1)\)). The 11th matches \(D > M > F\) and \(D > M\). The 12th matches \(F > D > M\) but also \(D > M\). Etc.

3.6.2 Order effects and modules with multiple genes

Consider the following case:

ofe1 <- allFitnessEffects(
    orderEffects = c("F > D" = -0.3, "D > F" = 0.4),
    geneToModule =
        c("F" = "f1, f2",
          "D" = "d1, d2") )

ag <- evalAllGenotypes(ofe1, order = TRUE)

There are four genes, \(d1, d2, f1, f2\), where each \(d\) belongs to module \(D\) and each \(f\) belongs to module \(F\).

What to expect for cases such as \(d1 > f1\) or \(f1 > d1\) is clear, as shown in

ag[5:16,]
##    Genotype Fitness
## 5   d1 > d2     1.0
## 6   d1 > f1     1.4
## 7   d1 > f2     1.4
## 8   d2 > d1     1.0
## 9   d2 > f1     1.4
## 10  d2 > f2     1.4
## 11  f1 > d1     0.7
## 12  f1 > d2     0.7
## 13  f1 > f2     1.0
## 14  f2 > d1     0.7
## 15  f2 > d2     0.7
## 16  f2 > f1     1.0

Likewise, cases such as \(d1 > d2 > f1\) or \(f2 > f1 > d1\) are clear, because in terms of modules they map to $ D > F$ or \(F > D\): the observed order of mutation \(d1 > d2 > f1\) means that module \(D\) was mutated first and module \(F\) was mutated second. Similar for \(d1 > f1 > f2\) or \(f1 > d1 > d2\): those map to \(D > F\) and \(F > D\). We can see the fitness of those four case in:

ag[c(17, 39, 19, 29), ]
##        Genotype Fitness
## 17 d1 > d2 > f1     1.4
## 39 f2 > f1 > d1     0.7
## 19 d1 > f1 > d2     1.4
## 29 f1 > d1 > d2     0.7

and they correspond to the values of those order effects, where \(F > D = (1 - 0.3)\) and \(D > F = (1 + 0.4)\):

ag[c(17, 39, 19, 29), "Fitness"] == c(1.4, 0.7, 1.4, 0.7)
## [1] TRUE TRUE TRUE TRUE

What if we match several patterns? For example, \(d1 > f1 > d2 > f2\) and \(d1 > f1 > f2 > d2\)? The first maps to \(D > F > D > F\) and the second to \(D > F > D\). But since we are concerned with which one happened first and which happened second we should expect those two to correspond to the same fitness, that of pattern \(D > F\), as is the case:

ag[c(43, 44),]
##             Genotype Fitness
## 43 d1 > f1 > d2 > f2     1.4
## 44 d1 > f1 > f2 > d2     1.4
ag[c(43, 44), "Fitness"] == c(1.4, 1.4)
## [1] TRUE TRUE

More generally, that applies to all the patterns that start with one of the “d” genes:

all(ag[41:52, "Fitness"] == 1.4)
## [1] TRUE

Similar arguments apply to the opposite pattern, \(F > D\), which apply to all the possible gene mutation orders that start with one of the “f” genes. For example:

all(ag[53:64, "Fitness"] == 0.7)
## [1] TRUE

3.6.3 Order and modules with 325 genotypes

We can of course have more than two genes per module. This just repeats the above, with five genes (there are 325 genotypes, and that is why we pass the “max” argument to evalAllGenotypes, to allow for more than the default 256).

ofe2 <- allFitnessEffects(
    orderEffects = c("F > D" = -0.3, "D > F" = 0.4),
    geneToModule =
        c("F" = "f1, f2, f3",
          "D" = "d1, d2") )
ag2 <- evalAllGenotypes(ofe2, max = 325, order = TRUE)

We can verify that any combination that starts with a “d” gene and then contains at least one “f” gene will have a fitness of \(1+0.4\). And any combination that starts with an “f” gene and contains at least one “d” genes will have a fitness of \(1 - 0.3\). All other genotypes have a fitness of 1:

all(ag2[grep("^d.*f.*", ag2[, 1]), "Fitness"] == 1.4)
## [1] TRUE
all(ag2[grep("^f.*d.*", ag2[, 1]), "Fitness"] == 0.7)
## [1] TRUE
oe <- c(grep("^f.*d.*", ag2[, 1]), grep("^d.*f.*", ag2[, 1]))
all(ag2[-oe, "Fitness"] == 1)
## [1] TRUE

3.6.4 Order effects and genes without interactions

We will now look at both order effects and interactions. To make things more interesting, we name genes so that the ordered names do split nicely between those with and those without order effects (this, thus, also serves as a test of messy orders of names).

foi1 <- allFitnessEffects(
    orderEffects = c("D>B" = -0.2, "B > D" = 0.3),
    noIntGenes = c("A" = 0.05, "C" = -.2, "E" = .1))

You can get a verbose view of what the gene names and modules are (and their automatically created numeric codes) by:

foi1[c("geneModule", "long.geneNoInt")]
## $geneModule
##   Gene Module GeneNumID ModuleNumID
## 1 Root   Root         0           0
## 2    B      B         1           1
## 3    D      D         2           2
## 
## $long.geneNoInt
##   Gene GeneNumID     s
## A    A         3  0.05
## C    C         4 -0.20
## E    E         5  0.10

We can get the fitness of all genotypes (we set \(max = 325\) because that is the number of possible genotypes):

agoi1 <- evalAllGenotypes(foi1,  max = 325, order = TRUE)
head(agoi1)
##   Genotype Fitness
## 1        B    1.00
## 2        D    1.00
## 3        A    1.05
## 4        C    0.80
## 5        E    1.10
## 6    B > D    1.30

Now:

rn <- 1:nrow(agoi1)
names(rn) <- agoi1[, 1]

agoi1[rn[LETTERS[1:5]], "Fitness"] == c(1.05, 1, 0.8, 1, 1.1)
## [1] TRUE TRUE TRUE TRUE TRUE

According to the fitness effects we have specified, we also know that any genotype with only two mutations, one of which is either “A”, “C” “E” and the other is “B” or “D” will have the fitness corresponding to “A”, “C” or “E”, respectively:

agoi1[grep("^A > [BD]$", names(rn)), "Fitness"] == 1.05
## [1] TRUE TRUE
agoi1[grep("^C > [BD]$", names(rn)), "Fitness"] == 0.8
## [1] TRUE TRUE
agoi1[grep("^E > [BD]$", names(rn)), "Fitness"] == 1.1
## [1] TRUE TRUE
agoi1[grep("^[BD] > A$", names(rn)), "Fitness"] == 1.05
## [1] TRUE TRUE
agoi1[grep("^[BD] > C$", names(rn)), "Fitness"] == 0.8
## [1] TRUE TRUE
agoi1[grep("^[BD] > E$", names(rn)), "Fitness"] == 1.1
## [1] TRUE TRUE

We will not be playing many additional games with regular expressions, but let us check those that start with “D” and have all the other mutations, which occupy rows 230 to 253; fitness should be equal (within numerical error, because of floating point arithmetic) to the order effect of “D” before “B” times the other effects \((1 - 0.3) * 1.05 * 0.8 * 1.1 = 0.7392\)

all.equal(agoi1[230:253, "Fitness"] ,
          rep((1 - 0.2) * 1.05 * 0.8 * 1.1, 24))
## [1] TRUE

and that will also be the value of any genotype with the five mutations where “D” comes before “B” such as those in rows 260 to 265, 277, or 322 and 323, but it will be equal to \((1 + 0.3) * 1.05 * 0.8 * 1.1 = 1.2012\) in those where “B” comes before “D”. Analogous arguments apply to four, three, and two mutation genotypes.

3.7.1 Epistasis: two alternative specifications

We want the following mapping of genotypes to fitness:

Suppose that the actual numerical values are \(s_a = 0.2, s_b = 0.3, s_{ab} = 0.7\).

We specify the above as follows:

sa <- 0.2
sb <- 0.3
sab <- 0.7

e2 <- allFitnessEffects(epistasis =
                            c("A: -B" = sa,
                              "-A:B" = sb,
                              "A : B" = sab))
evalAllGenotypes(e2, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT     1.0
## 2        A     1.2
## 3        B     1.3
## 4     A, B     1.7

That uses the “-” specification, so we explicitly exclude some patterns: with “A:-B” we say “A when there is no B”.

But we can also use a specification where we do not use the “-”. That requires a different numerical value of the interaction, because now, as we are rewriting the interaction term as genotype “A is mutant, B is mutant” the double mutant will incorporate the effects of “A mutant”, ``B mutant" and “both A and B mutants”. We can define a new \(s_2\) that satisfies \((1 + s_{ab}) = (1 + s_a) (1 + s_b) (1 + s_2)\) so \((1 + s_2) = (1 + s_{ab})/((1 + s_a) (1 + s_b))\) and therefore specify as:

s2 <- ((1 + sab)/((1 + sa) * (1 + sb))) - 1

e3 <- allFitnessEffects(epistasis =
                            c("A" = sa,
                              "B" = sb,
                              "A : B" = s2))
evalAllGenotypes(e3, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT     1.0
## 2        A     1.2
## 3        B     1.3
## 4     A, B     1.7

Note that this is the way you would specify effects with FFPopsim (Zanini and Neher 2012). Whether this specification or the previous one with “-” is simpler will depend on the model. For synthetic mortality and viability, I think the one using “-” is simpler to map genotype tables to fitness effects. See also section 3.7.2 and 3.7.3 and the example in section 5.4.2.

Finally, note that we can also specify some of these effects by combining the graph and the epistasis, as shown in section 5.2.1 or 5.4.2.

3.7.2 Epistasis with three genes and two alternative specifications

Suppose we have

where missing rows have a fitness of 1 (they have been deleted for conciseness). Note that the mutant for exactly A and C has a fitness that is the product of the individual terms (so there is no epistasis in that case).

sa <- 0.1
sb <- 0.15
sc <- 0.2
sab <- 0.3
sbc <- -0.25
sabc <- 0.4

sac <- (1 + sa) * (1 + sc) - 1

E3A <- allFitnessEffects(epistasis =
                            c("A:-B:-C" = sa,
                              "-A:B:-C" = sb,
                              "-A:-B:C" = sc,
                              "A:B:-C" = sab,
                              "-A:B:C" = sbc,
                              "A:-B:C" = sac,
                              "A : B : C" = sabc)
                                                )

evalAllGenotypes(E3A, order = FALSE, addwt = FALSE)
##   Genotype Fitness
## 1        A    1.10
## 2        B    1.15
## 3        C    1.20
## 4     A, B    1.30
## 5     A, C    1.32
## 6     B, C    0.75
## 7  A, B, C    1.40

We needed to pass the \(s_{ac}\) coefficient explicitly, even if it that term was just the product. We can try to avoid using the “-”, however (but we will need to do other calculations). For simplicity, I use capital ``S" in what follows where the letters differ from the previous specification:

sa <- 0.1
sb <- 0.15
sc <- 0.2
sab <- 0.3
Sab <- ( (1 + sab)/((1 + sa) * (1 + sb))) - 1
Sbc <- ( (1 + sbc)/((1 + sb) * (1 + sc))) - 1
Sabc <- ( (1 + sabc)/
          ( (1 + sa) * (1 + sb) * (1 + sc) *
            (1 + Sab) * (1 + Sbc) ) ) - 1

E3B <- allFitnessEffects(epistasis =
                             c("A" = sa,
                               "B" = sb,
                               "C" = sc,
                               "A:B" = Sab,
                               "B:C" = Sbc,
                               ## "A:C" = sac, ## not needed now
                               "A : B : C" = Sabc)
                                                )
evalAllGenotypes(E3B, order = FALSE, addwt = FALSE)
##   Genotype Fitness
## 1        A    1.10
## 2        B    1.15
## 3        C    1.20
## 4     A, B    1.30
## 5     A, C    1.32
## 6     B, C    0.75
## 7  A, B, C    1.40

The above two are, of course, identical:

all(evalAllGenotypes(E3A, order = FALSE, addwt = FALSE) == 
    evalAllGenotypes(E3B, order = FALSE, addwt = FALSE))
## [1] TRUE

We avoid specifying the “A:C”, as it just follows from the individual “A” and “C” terms, but given a specified genotype table, we need to do a little bit of addition and multiplication to get the coefficients.

3.7.3 Why can we specify some effects with a “-”?

Let’s suppose we want to specify the synthetic viability example seen before:

where “wt” denotes wild type and “M” denotes mutant.

If you want to directly map the above table to the fitness table for the program, to specify the genotype “A is wt, B is a mutant” you can specify it as “-A,B”, not just as “B”. Why? Because just the presence of a “B” is also compatible with genotype “A is mutant and B is mutant”. If you use “-” you are explicitly saying what should not be there so that “-A,B” is NOT compatible with “A, B”. Otherwise, you need to carefully add coefficients. Depending on what you are trying to model, different specifications might be simpler. See the examples in section 3.7.1 and 3.7.2. You have both options.

3.7.4 Epistasis: modules

There is nothing conceptually new, but we will show an example here:

sa <- 0.2
sb <- 0.3
sab <- 0.7

em <- allFitnessEffects(epistasis =
                            c("A: -B" = sa,
                              "-A:B" = sb,
                              "A : B" = sab),
                        geneToModule = c("A" = "a1, a2",
                                         "B" = "b1, b2"))
evalAllGenotypes(em, order = FALSE, addwt = TRUE)
##          Genotype Fitness
## 1              WT     1.0
## 2              a1     1.2
## 3              a2     1.2
## 4              b1     1.3
## 5              b2     1.3
## 6          a1, a2     1.2
## 7          a1, b1     1.7
## 8          a1, b2     1.7
## 9          a2, b1     1.7
## 10         a2, b2     1.7
## 11         b1, b2     1.3
## 12     a1, a2, b1     1.7
## 13     a1, a2, b2     1.7
## 14     a1, b1, b2     1.7
## 15     a2, b1, b2     1.7
## 16 a1, a2, b1, b2     1.7

Of course, we can do the same thing without using the “-”, as in section 3.7.1:

s2 <- ((1 + sab)/((1 + sa) * (1 + sb))) - 1

em2 <- allFitnessEffects(epistasis =
                            c("A" = sa,
                              "B" = sb,
                              "A : B" = s2),
                         geneToModule = c("A" = "a1, a2",
                                         "B" = "b1, b2")
                         )
evalAllGenotypes(em2, order = FALSE, addwt = TRUE)
##          Genotype Fitness
## 1              WT     1.0
## 2              a1     1.2
## 3              a2     1.2
## 4              b1     1.3
## 5              b2     1.3
## 6          a1, a2     1.2
## 7          a1, b1     1.7
## 8          a1, b2     1.7
## 9          a2, b1     1.7
## 10         a2, b2     1.7
## 11         b1, b2     1.3
## 12     a1, a2, b1     1.7
## 13     a1, a2, b2     1.7
## 14     a1, b1, b2     1.7
## 15     a2, b1, b2     1.7
## 16 a1, a2, b1, b2     1.7

3.9.1 A simple synthetic viability example

A simple and extreme example of synthetic viability is shown in the following table, where the joint mutant has fitness larger than the wild type, but each single mutant is lethal.

where “wt” denotes wild type and “M” denotes mutant.

We can specify this (setting \(s = 0.2\)) as (I play around with spaces, to show there is a certain flexibility with them):

s <- 0.2
sv <- allFitnessEffects(epistasis = c("-A : B" = -1,
                                      "A : -B" = -1,
                                      "A:B" = s))

Now, let’s look at all the genotypes (we use “addwt” to also get the wt, which by decree has fitness of 1), and disregard order:

(asv <- evalAllGenotypes(sv, order = FALSE, addwt = TRUE))
##   Genotype Fitness
## 1       WT     1.0
## 2        A     0.0
## 3        B     0.0
## 4     A, B     1.2

Asking the program to consider the order of mutations of course makes no difference:

evalAllGenotypes(sv, order = TRUE, addwt = TRUE)
##   Genotype Fitness
## 1       WT     1.0
## 2        A     0.0
## 3        B     0.0
## 4    A > B     1.2
## 5    B > A     1.2

Another example of synthetic viability is shown in section 5.2.2.

Of course, if multiple simultaneous mutations are not possible in the simulations, it is not possible to go from the wildtype to the double mutant in this model where the single mutants are not viable.

3.9.2 Synthetic viability, non-zero fitness, and modules

This is a slightly more elaborate case, where there is one module and the single mutants have different fitness between themselves, which is non-zero. Without the modules, this is the same as in Misra, Szczurek, and Vingron (2014), Figure 1b, which we go over in section 5.2.

where \(s_a, s_b < 0\) but \(s_{ab} > 0\).

sa <- -0.1
sb <- -0.2
sab <- 0.25
sv2 <- allFitnessEffects(epistasis = c("-A : B" = sb,
                             "A : -B" = sa,
                             "A:B" = sab),
                         geneToModule = c(
                             "A" = "a1, a2",
                             "B" = "b"))
evalAllGenotypes(sv2, order = FALSE, addwt = TRUE)
##    Genotype Fitness
## 1        WT    1.00
## 2        a1    0.90
## 3        a2    0.90
## 4         b    0.80
## 5    a1, a2    0.90
## 6     a1, b    1.25
## 7     a2, b    1.25
## 8 a1, a2, b    1.25

And if we look at order, of course it makes no difference:

evalAllGenotypes(sv2, order = TRUE, addwt = TRUE)
##       Genotype Fitness
## 1           WT    1.00
## 2           a1    0.90
## 3           a2    0.90
## 4            b    0.80
## 5      a1 > a2    0.90
## 6       a1 > b    1.25
## 7      a2 > a1    0.90
## 8       a2 > b    1.25
## 9       b > a1    1.25
## 10      b > a2    1.25
## 11 a1 > a2 > b    1.25
## 12 a1 > b > a2    1.25
## 13 a2 > a1 > b    1.25
## 14 a2 > b > a1    1.25
## 15 b > a1 > a2    1.25
## 16 b > a2 > a1    1.25

3.11.1 Synthetic viability using Bozic model

If we were to use the above specification with Bozic’s models, we might not get what we think we should get:

evalAllGenotypes(sv, order = FALSE, addwt = TRUE, model = "Bozic")
##   Genotype Death_rate
## 1       WT        1.0
## 2        A        2.0
## 3        B        2.0
## 4     A, B        0.8

What gives here? The simulation code would alert you of this (see section 3.11.2) in this particular case because there are ``-1“, which might indicate that this is not what you want. The problem is that you probably want the Death rate to be infinity (the birth rate was 0, so no clone viability, when we used birth rates —section 3.2.2).

Let us say so explicitly:

s <- 0.2
svB <- allFitnessEffects(epistasis = c("-A : B" = -Inf,
                                      "A : -B" = -Inf,
                                      "A:B" = s))
evalAllGenotypes(svB, order = FALSE, addwt = TRUE, model = "Bozic")
##   Genotype Death_rate
## 1       WT        1.0
## 2        A        Inf
## 3        B        Inf
## 4     A, B        0.8

Likewise, values of \(s\) larger than one have no effect beyond setting \(s = 1\) (a single term of \((1 - 1)\) will drive the product to 0, and as we cannot allow negative death rates negative values are set to 0):

s <- 1
svB1 <- allFitnessEffects(epistasis = c("-A : B" = -Inf,
                                       "A : -B" = -Inf,
                                       "A:B" = s))

evalAllGenotypes(svB1, order = FALSE, addwt = TRUE, model = "Bozic")
##   Genotype Death_rate
## 1       WT          1
## 2        A        Inf
## 3        B        Inf
## 4     A, B          0


s <- 3
svB3 <- allFitnessEffects(epistasis = c("-A : B" = -Inf,
                                       "A : -B" = -Inf,
                                       "A:B" = s))

evalAllGenotypes(svB3, order = FALSE, addwt = TRUE, model = "Bozic")
##   Genotype Death_rate
## 1       WT          1
## 2        A        Inf
## 3        B        Inf
## 4     A, B          0

Of course, death rates of 0.0 are likely to lead to trouble down the road, when we actually conduct simulations (see section 3.11.2).

3.11.2 Numerical issues with death rates of 0 in Bozic model

As we mentioned above (section 3.11.1) death rates of 0 can lead to trouble when using Bozic’s model:

i1 <- allFitnessEffects(noIntGenes = c(1, 0.5))
evalAllGenotypes(i1, order = FALSE, addwt = TRUE, 
                 model = "Bozic")
##   Genotype Death_rate
## 1       WT        1.0
## 2        1        0.0
## 3        2        0.5
## 4     1, 2        0.0
                 
i1_b <- oncoSimulIndiv(i1, model = "Bozic")
## Warning in nr_oncoSimul.internal(rFE = fp, birth = birth, death
## = death, : You are using a Bozic model with the new restriction
## specification, and you have at least one s of 1. If that gene is
## mutated, this will lead to a death rate of 0 and the simulations
## will abort when you get a non finite value.
## 
##  DEBUG2: Value of rnb = nan
## 
##  DEBUG2: Value of m = 1
## 
##  DEBUG2: Value of pe = 0
## 
##  DEBUG2: Value of pm = 0.999999
## 
##  this is spP
## 
##  popSize = 1
##  birth = 1
##  death = 0
##  W = 1
##  R = 1
##  mutation = 1e-06
##  timeLastUpdate = 1096.35
##  absfitness = -inf
##  numMutablePos =1
## 
##  Unrecoverable exception: Algo 2: retval not finite. Aborting.

Of course, there is no problem in using the above with other models:

evalAllGenotypes(i1, order = FALSE, addwt = TRUE, 
                 model = "Exp")
##   Genotype Fitness
## 1       WT     1.0
## 2        1     2.0
## 3        2     1.5
## 4     1, 2     3.0
i1_e <- oncoSimulIndiv(i1, model = "Exp")
summary(i1_e)
##   NumClones TotalPopSize LargestClone MaxNumDrivers MaxDriversLast
## 1         3    181659523    181588631             0              0
##   NumDriversLargestPop TotalPresentDrivers FinalTime NumIter
## 1                    0                   0       204     551
##   HittedWallTime HittedMaxTries errorMF minDMratio minBMratio
## 1          FALSE          FALSE      NA      5e+05      5e+05
##   OccurringDrivers
## 1

5.1.1 Using a DAG

We cannot find a simple mapping from their expressions to our fitness parameterization, but we can get fairly close by using a DAG; in this one, note the unusual feature of having one of the “s” terms (that for the driver dependency on root) be negative. Using the parameters given in the legend of their Figure 3 for \(s_p, S_D^+, S_D^-, S_{DP}\) and obtaining that negative value for the dependency of the driver on root we can do:

K <- 4
sp <- 1e-5
sdp <- 0.015
sdplus <- 0.05
sdminus <- 0.1
cnt <- (1 + sdplus)/(1 + sdminus)
prod_cnt <- cnt - 1
bauer <- data.frame(parent = c("Root", rep("D", K)),
                    child = c("D", paste0("s", 1:K)),
                    s = c(prod_cnt, rep(sdp, K)),
                    sh = c(0, rep(sp, K)),
                    typeDep = "MN")
fbauer <- allFitnessEffects(bauer)
(b1 <- evalAllGenotypes(fbauer, order = FALSE, addwt = TRUE))
##             Genotype   Fitness
## 1                 WT 1.0000000
## 2                  D 0.9545455
## 3                 s1 1.0000100
## 4                 s2 1.0000100
## 5                 s3 1.0000100
## 6                 s4 1.0000100
## 7              D, s1 0.9688636
## 8              D, s2 0.9688636
## 9              D, s3 0.9688636
## 10             D, s4 0.9688636
## 11            s1, s2 1.0000200
## 12            s1, s3 1.0000200
## 13            s1, s4 1.0000200
## 14            s2, s3 1.0000200
## 15            s2, s4 1.0000200
## 16            s3, s4 1.0000200
## 17         D, s1, s2 0.9833966
## 18         D, s1, s3 0.9833966
## 19         D, s1, s4 0.9833966
## 20         D, s2, s3 0.9833966
## 21         D, s2, s4 0.9833966
## 22         D, s3, s4 0.9833966
## 23        s1, s2, s3 1.0000300
## 24        s1, s2, s4 1.0000300
## 25        s1, s3, s4 1.0000300
## 26        s2, s3, s4 1.0000300
## 27     D, s1, s2, s3 0.9981475
## 28     D, s1, s2, s4 0.9981475
## 29     D, s1, s3, s4 0.9981475
## 30     D, s2, s3, s4 0.9981475
## 31    s1, s2, s3, s4 1.0000400
## 32 D, s1, s2, s3, s4 1.0131198

(We use “D” for “driver” or “primary driver”, as is it is called in the original paper, and “s” for secondary drivers, somewhat similar to passengers).

Note that what we specify as “typeDep” is irrelevant (MN, SMN, or XMPN make no difference).

This is the DAG:

plot(fbauer)

And if you compare the tabular output of evalAllGenotypes you can see that the values of fitness reproduces the fitness landscape that they show in their Figure 1. We can also use our plot for fitness landscapes:

plot(b1, use_ggrepel = TRUE)

5.1.2 Specifying fitness of genotypes directly

An alternative approach to specify the fitness, if the number of genotypes is reasonably small, is to directly evaluate fitness as given by their expressions. Then, use the genotFitness argument to allFitnessEffects.

We will create all possible genotypes; then we will write a function that gives the fitness of each genotype according to their expression; finally, we will call this function on the data frame of genotypes, and pass this data frame to allFitnessEffects.

m1 <- expand.grid(D = c(1, 0), s1 = c(1, 0), s2 = c(1, 0),
                  s3 = c(1, 0), s4 = c(1, 0))

fitness_bauer <- function(D, s1, s2, s3, s4, 
                          sp = 1e-5, sdp = 0.015, sdplus = 0.05,
                          sdminus = 0.1) {
    if(!D) {
        b <- 0.5 * ( (1 + sp)^(sum(c(s1, s2, s3, s4))))
    } else {
        b <- 0.5 * 
            (((1 + sdplus)/(1 + sdminus)  *
              (1 + sdp)^(sum(c(s1, s2, s3, s4)))))
    }
    fitness <- b - (1 - b)
    our_fitness <- 1 + fitness ## prevent negative fitness and
    ## make wt fitness = 1
    return(our_fitness)
}

m1$Fitness <- 
    apply(m1, 1, function(x) do.call(fitness_bauer, as.list(x)))

bauer2 <- allFitnessEffects(genotFitness = m1)

Now, show the fitness of all genotypes:

evalAllGenotypes(bauer2, order = FALSE, addwt = TRUE)
##             Genotype   Fitness
## 1                 WT 1.0000000
## 2                  D 0.9545455
## 3                 s1 1.0000100
## 4                 s2 1.0000100
## 5                 s3 1.0000100
## 6                 s4 1.0000100
## 7              D, s1 0.9688636
## 8              D, s2 0.9688636
## 9              D, s3 0.9688636
## 10             D, s4 0.9688636
## 11            s1, s2 1.0000200
## 12            s1, s3 1.0000200
## 13            s1, s4 1.0000200
## 14            s2, s3 1.0000200
## 15            s2, s4 1.0000200
## 16            s3, s4 1.0000200
## 17         D, s1, s2 0.9833966
## 18         D, s1, s3 0.9833966
## 19         D, s1, s4 0.9833966
## 20         D, s2, s3 0.9833966
## 21         D, s2, s4 0.9833966
## 22         D, s3, s4 0.9833966
## 23        s1, s2, s3 1.0000300
## 24        s1, s2, s4 1.0000300
## 25        s1, s3, s4 1.0000300
## 26        s2, s3, s4 1.0000300
## 27     D, s1, s2, s3 0.9981475
## 28     D, s1, s2, s4 0.9981475
## 29     D, s1, s3, s4 0.9981475
## 30     D, s2, s3, s4 0.9981475
## 31    s1, s2, s3, s4 1.0000400
## 32 D, s1, s2, s3, s4 1.0131198

Can we use modules in this example, if we use the “lego system”? Sure, as in any other case.

5.2.1 Example 1.a

In that figure it is evident that the fitness effect of “A” and “B” are the same. There are two different models depending on whether “AB” is just the product of both, or there is epistasis. In the first case probably the simplest is:

s <- 0.1 ## or whatever number
m1a1 <- allFitnessEffects(data.frame(parent = c("Root", "Root"),
                                     child = c("A", "B"),
                                     s = s,
                                     sh = 0,
                                     typeDep = "MN"))
evalAllGenotypes(m1a1, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT    1.00
## 2        A    1.10
## 3        B    1.10
## 4     A, B    1.21

If the double mutant shows epistasis, as we saw before (section 3.7.1) we have a range of options. For example:

s <- 0.1
sab <- 0.3
m1a2 <- allFitnessEffects(epistasis = c("A:-B" = s,
                                        "-A:B" = s,
                                        "A:B" = sab))
evalAllGenotypes(m1a2, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT     1.0
## 2        A     1.1
## 3        B     1.1
## 4     A, B     1.3

But we could also modify the graph dependency structure, and we have to change the value of the coefficient, since that is what multiplies each of the terms for “A” and “B”: \((1 + s_{AB}) = (1 + s)^2(1 + s_{AB3})\)

sab3 <- ((1 + sab)/((1 + s)^2)) - 1
m1a3 <- allFitnessEffects(data.frame(parent = c("Root", "Root"),
                                     child = c("A", "B"),
                                     s = s,
                                     sh = 0,
                                     typeDep = "MN"),
                          epistasis = c("A:B" = sab3))
evalAllGenotypes(m1a3, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT     1.0
## 2        A     1.1
## 3        B     1.1
## 4     A, B     1.3

And, obviously

all.equal(evalAllGenotypes(m1a2, order = FALSE, addwt = TRUE),
          evalAllGenotypes(m1a3, order = FALSE, addwt = TRUE))
## [1] TRUE

5.2.2 Example 1.b

This is a specific case of synthetic viability (see also section 3.9):

Here, \(S_A, S_B < 0\), \(S_B < 0\), \(S_{AB} > 0\) and \((1 + S_{AB}) (1 + S_A) (1 + S_B) > 1\).

As before, we can specify this in several different ways. The simplest is to specify all genotypes:

sa <- -0.6
sb <- -0.7
sab <- 0.3
m1b1 <- allFitnessEffects(epistasis = c("A:-B" = sa,
                                        "-A:B" = sb,
                                        "A:B" = sab))
evalAllGenotypes(m1b1, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT     1.0
## 2        A     0.4
## 3        B     0.3
## 4     A, B     1.3

We could also use a tree and modify the “sab” for the epistasis, as before (5.2.1).

5.2.3 Example 1.c

The final case, in figure 1.c of Misra et al., is just epistasis, where a mutation in one of the genes is deleterious (possibly only mildly), in the other is beneficial, and the double mutation has fitness larger than any of the other two.

Here we have that \(s_A > 0\), \(s_B < 0\), \((1 + s_{AB}) (1 + s_A) (1 + s_B) > (1 + s_{AB})\) so \(s_{AB} > \frac{-s_B}{1 + s_B}\)

As before, we can specify this in several different ways. The simplest is to specify all genotypes:

sa <- 0.2
sb <- -0.3
sab <- 0.5
m1c1 <- allFitnessEffects(epistasis = c("A:-B" = sa,
                                        "-A:B" = sb,
                                        "A:B" = sab))
evalAllGenotypes(m1c1, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT     1.0
## 2        A     1.2
## 3        B     0.7
## 4     A, B     1.5

We could also use a tree and modify the “sab” for the epistasis, as before (5.2.1).

5.4.1 Figure 1.a

where the “1” and “2” in the figure refer to the total number of mutations in two different loci. This is, therefore, very similar to the example in section 5.2.2. Here we have, in their notation, \(\delta_1 < 0\), fitness of single “A” or single “B” = \(1 + \delta_1\), \(S_{AB} > 0\), \((1 + S_{AB})(1 + \delta_1)^2 > 1\).

5.4.2 Figure 1.b

In their figure 1b they show

Where, as before, 1, 2, 3, denote the total number of mutations over three different loci and \(\delta_1 < 0\), \(\delta_2 < 0\), fitness of single mutant is \((1 + \delta_1)\), of double mutant is \((1 + \delta_2)\) so that \((1 + \delta_2) = (1 + \delta_1)^2 (1 + s_2)\) and of triple mutant is \((1 + \delta_3)\), so that \((1 + \delta_3) = (1 + \delta_1)^3 (1 + s_2)^3 (1 + s_3)\).

We can specify this combining a graph with epistasis:

d1 <- -0.05 ## single mutant fitness 0.95
d2 <- -0.08 ## double mutant fitness 0.92
d3 <- 0.2   ## triple mutant fitness 1.2

s2 <- ((1 + d2)/(1 + d1)^2) - 1
s3 <- ( (1 + d3)/((1 + d1)^3 * (1 + s2)^3) ) - 1

w <- allFitnessEffects(
    data.frame(parent = c("Root", "Root", "Root"),
               child = c("A", "B", "C"),
               s = d1,
               sh = -1,
               typeDep = "MN"),
    epistasis = c(
        "A:B" = s2,
        "A:C" = s2,
        "B:C" = s2,
        "A:B:C" = s3))

The model can be shown graphically as:

plot(w)

And fitness of all genotypes is:

evalAllGenotypes(w, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT    1.00
## 2        A    0.95
## 3        B    0.95
## 4        C    0.95
## 5     A, B    0.92
## 6     A, C    0.92
## 7     B, C    0.92
## 8  A, B, C    1.20

Alternatively, we can directly specify what each genotype adds to the fitness, given the included genotype. This is basically replacing the graph by giving each of “A”, “B”, and “C” directly:

wb <- allFitnessEffects(
    epistasis = c(
        "A" = d1,
        "B" = d1,
        "C" = d1,
        "A:B" = s2,
        "A:C" = s2,
        "B:C" = s2,
        "A:B:C" = s3))

evalAllGenotypes(wb, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT    1.00
## 2        A    0.95
## 3        B    0.95
## 4        C    0.95
## 5     A, B    0.92
## 6     A, C    0.92
## 7     B, C    0.92
## 8  A, B, C    1.20

The plot, of course, is not very revealing and we cannot show that there is a three-way interaction (only all three two-way interactions):

plot(wb)

As we have seen several times already (sections 3.7.1, 3.7.2, 3.7.3) we can also give the genotypes directly and, consequently, the fitness of each genotype (not the added contribution):

wc <- allFitnessEffects(
    epistasis = c(
        "A:-B:-C" = d1,
        "B:-C:-A" = d1,
        "C:-A:-B" = d1,
        "A:B:-C" = d2,
        "A:C:-B" = d2,
        "B:C:-A" = d2,
        "A:B:C" = d3))
evalAllGenotypes(wc, order = FALSE, addwt = TRUE)
##   Genotype Fitness
## 1       WT    1.00
## 2        A    0.95
## 3        B    0.95
## 4        C    0.95
## 5     A, B    0.92
## 6     A, C    0.92
## 7     B, C    0.92
## 8  A, B, C    1.20

6.3.1 Ending the simulations: conditions

• onlyCancer = TRUE. A simulation will be repeated until any one of the “reach cancer” conditions is met, if this happens before the simulation reaches finalTime¹¹. These conditions are:

  • Total population size becomes larger than detectionSize.
  • The number of drivers in any one genotype or clone becomes equal to, or larger than, detectionDrivers; note that this allows you to stop the simulation as soon as a specific genotype is found, by using exactly and only the genes that make that genotype as the drivers.
  • A gene or gene combination among those listed in fixation becomes fixed in the population (i.e., has a frequency is 1) (see details in (6.3.3 and 6.3.4).
  • The tumor is detected according to a stochastic detection mechanism, where the probability of “detecting the tumor” increases with population size; this is explained below (6.3.2) and is controlled by argument detectionProb.

  As we exit as soon as any of the exiting conditions is reached, if you only care about one condition, set the other to NA (see also section 6.3.2.1).

• onlyCancer = FALSE. A simulation will run only once, and will exit as soon as any of the above conditions are met or as soon as the total population size becomes zero or we reach finalTime.

As an example of onlyCancer = TRUE, focusing on the first two mechanisms, suppose you give detectionSize = 1e4 and detectionDrivers =3 (and you have detectionProb = NA). A simulation will exit as soon as it reaches a total population size of \(10^4\) or any clone has four drivers, whichever comes first (if any of these happen before finalTime).

In the onlyCancer = TRUE case, what happens if we reach finalTime (or the population size becomes zero) before any of the ``reach cancer" conditions have been fulfilled? The simulation will be repeated again, within the following limits:

• max.wall.time: the total wall time we allow an individual simulation to run;
• max.num.tries: the maximum number of times we allow a simulation to be repeated to reach cancer;
• max.wall.time.total and max.num.tries.total, similar to the above but over a set of simulations in function oncoSimulSample.

Incidentally, we keep track of the number of attempts used (the component other$attemptsUsed$) before we reach cancer, so you can estimate (as from a negative binomial sampling) the probability of reaching your desired end point under different scenarios.

The onlyCancer = FALSE case might be what you want to do when you examine general population genetics scenarios without focusing on possible sampling issues. To do this, set finalTime to the value you want and set onlyCancer = FALSE; in addition, set detectionProb to “NA” and detectionDrivers and detectionSize to “NA” or to huge numbers¹². In this scenario you simply collect the simulation output at the end of the run, regardless of what happened with the population (it became extinct, it did not reach a large size, it did not accumulate drivers, etc).

6.3.2 Stochastic detection mechanism: “detectionProb”

This is the process that is controlled by the argument detectionProb. Here the probability of tumor detection increases with the total population size. This is biologically a reasonable assumption: the larger the tumor, the more likely it is it will be detected.

At regularly spaced times during the simulation, we compute the probability of detection as a function of size and determine (by comparing against a random uniform number) if the simulation should finish. For simplicity, and to make sure the probability is bounded between 0 and 1, we use the function

where \(P(N)\) is the probability that a tumor with a population size \(N\) will be detected, and \(c\) (argument \(cPDetect\) in the oncoSimul* functions) controls how fast \(P(N)\) increases with increasing population size relative to a baseline, \(B\) (\(PDBaseline\) in the oncoSimul* functions); with \(B\) we both control the minimal population size at which this mechanism stats operating (because we will rarely want detection unless there is some meaningful increase of population size over initSize) and we model the increase in \(P(N)\) as a function of relative differences with respect to \(B\). (Note that this is a major change in version 2.9.9. Before version 2.9.9, the expression used was \(P(N) = 1 - e^{ -c ( N - B)}\), so we did not make the increase relative to \(B\); of course, you can choose an appropriate \(c\) to make different models comparable, but the expression used before 2.9.9 made it much harder to compare simulations with very different initial population sizes, as baselines are often naturall a function of initial population sizes.)

The \(P(N)\) refers to the probability of detection at each one of the occasions when we assess the probability of exiting. When, or how often, do we do that? When we assess probability of exiting is controlled by checkSizePEvery, which will often be much larger than sampleEvery¹³. Biologically, a way to think of checkSizePEvery is “time between doctor appointments”.

An important warning, though: for populations that are growing very, very fast or where some genes might have very large effects on fitness even a moderate checkSizePEvery of, say, 10, might be inappropriate, since populations could have increased by several orders of magnitude between successive checks. This issue is also discussed in section 2.1.1 and 2.2.

Finally, you can specify \(c\) (\(cPDetect\)) directly (you will need to set n2 and p2 to NA). However, it might be more intuitive to specify the pair n2, p2, such that \(P(n2) = p2\) (and from that pair we solve for the value of \(cPDetect\)).

You can get a feeling for the effects of these arguments by playing with the following code, that we do not execute here for the sake of speed. Here no mutation has any effect, but there is a non-zero probability of exiting as soon as the total population size becomes larger than the initial population size. So, eventually, all simulations will exit and, as we are using the McFarland model, population size will vary slightly around the initial population size.

gi2 <- rep(0, 5)
names(gi2) <- letters[1:5]
oi2 <- allFitnessEffects(noIntGenes = gi2)
s5 <- oncoSimulPop(200,
                   oi2,
                   model = "McFL",
                   initSize = 1000,
                   detectionProb = c(p2 = 0.1, 
                                     n2 = 2000, 
                                     PDBaseline = 1000,
                                     checkSizePEvery = 2), 
                   detectionSize = NA,
                   finalTime = NA, 
                   keepEvery = NA,
                   detectionDrivers = NA)
s5
hist(unlist(lapply(s5, function(x) x$FinalTime)))

As you decrease checkSizePEvery the distribution of “FinalTime” will resemble more and more an exponential distribution.

In this vignette, there are some further examples of using this mechanism in 6.5.4 and 6.5.2, with the default arguments.

6.3.3 Fixation of genes/gene combinations

In some cases we might be interested in running simulations until a particular set of genes, or gene combinations, reaches fixation. This exit condition might be more relevant than some of the above in many non-cancer-related evolutionary genetics scenarios.

Simulations will stop as soon as any of the genes or gene combinations in the vector (or list) fixation reaches a frequency of 1. These gene combinations might have non-zero intersection (i.e., they might share genes), and those genes need not be drivers. If we want simulations to only stop when fixation of those genes/gene combinations is reached, we will set all other stopping conditions to NA. It is, of course, up to you to ensure that those stopping conditions are reasonable (that they can be reached) and to use, or not, finalTime; otherwise, simulations will eventually abort (e.g., when max.wall.time or max.num.tries are reached). Since we are asking for fixation, the Exp or Bozic models will often not be appropriate here; instead, models with competition such as McFL are more appropriate.

We return here to the example from section 5.3.

u <- 0.2; i <- -0.02; vi <- 0.6; ui <- uv <- -Inf
od2 <- allFitnessEffects(
    epistasis = c("u" = u,  "u:i" = ui,
                  "u:v" = uv, "i" = i,
                  "v:-i" = -Inf, "v:i" = vi))

Ochs and Desai explain that “Each simulated population was evolved until either the uphill genotype or valley-crossing genotype fixed.” (see Ochs and Desai (2015), p.2, section “Simulations”). We will do the same here. We specify that we want to end the simulation when either the “u” or the “v, i” genotypes have reached fixation, by passing those genotype combinations as the fixation argument (in this example using fixation = c("u", "v") would have been equivalent, since the “v” genotype by itself has fitness of 0).

We want to be explicit that fixation will be the one and only condition for ending the simulations, and thus we set arguments detectionDrivers, finalTime, detectionSize and detectionProb explicitly to NA. (We set the number of repetitions only to 10 for the sake of speed when creating the vignette).

initS <- 20
## We use only a small number of repetitions for the sake
## of speed.
od100 <- oncoSimulPop(10, od2,
                      fixation = c("u", "v, i"),
                      model = "McFL",
                      mu = 1e-4,
                      detectionDrivers = NA,
                      finalTime = NA,
                      detectionSize = NA,
                      detectionProb = NA,
                      onlyCancer = TRUE,
                      initSize = initS, 
                      mc.cores = 2)

What is the frequency of each genotype among the simulations? (or, what is the frequency of fixation of each genotype?)

sampledGenotypes(samplePop(od100))                      
## 
##  Subjects by Genes matrix of 10 subjects and 3 genes.
##   Genotype Freq
## 1     i, v    3
## 2        u    7
## 
##  Shannon's diversity (entropy) of sampled genotypes: 0.6108643

Note the very large variability in reaching fixation

head(summary(od100)[, c(1:3, 8:9)])
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1         3           29           29  1502.825   60123
## 2         3           18           18  8321.300  332896
## 3         3           24           24  1728.200   69137
## 4         4           47           47  4289.600  171601
## 5         3           30           30  2003.950   80171
## 6         3           38           38   995.600   39831

6.3.4 Fixation of genotypes

Section 6.3.3 deals with the fixation of gene/gene combinations. What if you want fixation on specific genotypes? To give an example, suppose we have a five loci genotype and suppose that you want to stop the simulations only if genotypes “A”, “B, E”, or “A, B, C, D, E” reach fixation. You do not want to stop it if, say, genotype “A, B, E” reaches fixation. To specify genotypes, you prepend the genotype combinations with a ``_,’’, and that tells OncoSimulR that you want fixation of genotypes, not just gene combinations.

An example of the differences between the mechanisms can be seen from this code:

## Create a simple fitness landscape
rl1 <- matrix(0, ncol = 6, nrow = 9)
colnames(rl1) <- c(LETTERS[1:5], "Fitness")
rl1[1, 6] <- 1
rl1[cbind((2:4), c(1:3))] <- 1
rl1[2, 6] <- 1.4
rl1[3, 6] <- 1.32
rl1[4, 6] <- 1.32
rl1[5, ] <- c(0, 1, 0, 0, 1, 1.5)
rl1[6, ] <- c(0, 0, 1, 1, 0, 1.54)
rl1[7, ] <- c(1, 0, 1, 1, 0, 1.65)
rl1[8, ] <- c(1, 1, 1, 1, 0, 1.75)
rl1[9, ] <- c(1, 1, 1, 1, 1, 1.85)
class(rl1) <- c("matrix", "genotype_fitness_matrix")
## plot(rl1) ## to see the fitness landscape

## Gene combinations
local_max_g <- c("A", "B, E", "A, B, C, D, E")
## Specify the genotypes
local_max <- paste0("_,", local_max_g)

fr1 <- allFitnessEffects(genotFitness = rl1, drvNames = LETTERS[1:5])
initS <- 2000


######## Stop on gene combinations   #####
r1 <- oncoSimulPop(10,
                       fp = fr1,
                       model = "McFL",
                       initSize = initS,
                       mu = 1e-4,
                       detectionSize = NA,
                       sampleEvery = .03,
                       keepEvery = 1, 
                       finalTime = 50000,
                       fixation = local_max_g, 
                       detectionDrivers = NA,
                       detectionProb = NA,
                       onlyCancer = TRUE,
                       max.num.tries = 500,
                       max.wall.time = 20, 
                       errorHitMaxTries = TRUE,
                       keepPhylog = FALSE,
                       mc.cores = 2)
sp1 <- samplePop(r1, "last", "singleCell")
## 
##  Subjects by Genes matrix of 10 subjects and 5 genes.
sgsp1 <- sampledGenotypes(sp1)
## often you will stop on gene combinations that
## are not local maxima in the fitness landscape
sgsp1
##     Genotype Freq
## 1          A    4
## 2 A, B, C, D    3
## 3    A, C, D    2
## 4       B, E    1
## 
##  Shannon's diversity (entropy) of sampled genotypes: 1.279854
sgsp1$Genotype %in% local_max_g
## [1]  TRUE FALSE FALSE  TRUE

####### Stop on genotypes   ####

r2 <- oncoSimulPop(10,
                       fp = fr1,
                       model = "McFL",
                       initSize = initS,
                       mu = 1e-4,
                       detectionSize = NA,
                       sampleEvery = .03,
                       keepEvery = 1, 
                       finalTime = 50000,
                       fixation = local_max, 
                       detectionDrivers = NA,
                       detectionProb = NA,
                       onlyCancer = TRUE,
                       max.num.tries = 500,
                       max.wall.time = 20, 
                       errorHitMaxTries = TRUE,
                       keepPhylog = FALSE,
                       mc.cores = 2)
## All final genotypes should be local maxima                       
sp2 <- samplePop(r2, "last", "singleCell")
## 
##  Subjects by Genes matrix of 10 subjects and 5 genes.
sgsp2 <- sampledGenotypes(sp2)
sgsp2$Genotype %in% local_max_g
## [1] TRUE TRUE TRUE

6.3.5 Fixation: tolerance, number of periods, minimal size

In particular if you specify stopping on genotypes, you might want to think about three additional parameters: fixation_tolerance, min_successive_fixation, and fixation_min_size.

fixation_tolerance: fixation is considered to have happened if the genotype/gene combinations specified as genotypes/gene combinations for fixation have reached a frequency \(> 1 - fixation\_tolerance\). (The default is 0, so we ask for genotypes/gene combinations with a frequency of 1, which might not be what you want with large mutation rates and complex fitness landscape with genotypes of similar fitness.)

min_successive_fixation: during how many successive sampling periods the conditions of fixation need to be fulfilled before declaring fixation. These must be successive sampling periods without interruptions (i.e., a single period when the condition is not fulfilled will set the counter to 0). This can help to exclude short, transitional, local maxima that are quickly replaced by other genotypes. (The default is 50, but this is probably too small for ``real life’’ usage).

fixation_min_size: you might only want to consider fixation to have happened if a minimal size has been reached (this can help weed out local maxima that have fitness that is barely above that of the wild-type genotype). (The default is 0).

An example of using those options:

## Create a simple fitness landscape
rl1 <- matrix(0, ncol = 6, nrow = 9)
colnames(rl1) <- c(LETTERS[1:5], "Fitness")
rl1[1, 6] <- 1
rl1[cbind((2:4), c(1:3))] <- 1
rl1[2, 6] <- 1.4
rl1[3, 6] <- 1.32
rl1[4, 6] <- 1.32
rl1[5, ] <- c(0, 1, 0, 0, 1, 1.5)
rl1[6, ] <- c(0, 0, 1, 1, 0, 1.54)
rl1[7, ] <- c(1, 0, 1, 1, 0, 1.65)
rl1[8, ] <- c(1, 1, 1, 1, 0, 1.75)
rl1[9, ] <- c(1, 1, 1, 1, 1, 1.85)
class(rl1) <- c("matrix", "genotype_fitness_matrix")
## plot(rl1) ## to see the fitness landscape

## The local fitness maxima are
## c("A", "B, E", "A, B, C, D, E")

fr1 <- allFitnessEffects(genotFitness = rl1, drvNames = LETTERS[1:5])
initS <- 2000

## Stop on genotypes

r3 <- oncoSimulPop(10,
                  fp = fr1,
                  model = "McFL",
                  initSize = initS,
                  mu = 1e-4,
                  detectionSize = NA,
                  sampleEvery = .03,
                  keepEvery = 1, 
                  finalTime = 50000,
                  fixation = c(paste0("_,",
                                   c("A", "B, E", "A, B, C, D, E")),
                               fixation_tolerance = 0.1,
                               min_successive_fixation = 200,
                               fixation_min_size = 3000),
                  detectionDrivers = NA,
                  detectionProb = NA,
                  onlyCancer = TRUE,
                  max.num.tries = 500,
                  max.wall.time = 20, 
                  errorHitMaxTries = TRUE,
                  keepPhylog = FALSE,
                  mc.cores = 2)

6.3.6 Mixing stopping on gene combinations and genotypes

This would probably be awfully confusing and is not tested formally (though it should work). Let me know if you think this is an important feature. (Pull requests with tests welcome.)

6.5.1 Bauer’s example again

We will use the model of Bauer, Siebert, and Traulsen (2014) that we saw in section 5.1.

K <- 5
sd <- 0.1
sdp <- 0.15
sp <- 0.05
bauer <- data.frame(parent = c("Root", rep("p", K)),
                    child = c("p", paste0("s", 1:K)),
                    s = c(sd, rep(sdp, K)),
                    sh = c(0, rep(sp, K)),
                    typeDep = "MN")
fbauer <- allFitnessEffects(bauer, drvNames = "p")
set.seed(1)
## Use fairly large mutation rate
b1 <- oncoSimulIndiv(fbauer, mu = 5e-5, initSize = 1000,
                     finalTime = NA,
                     onlyCancer = TRUE,
                     detectionProb = "default")

We will now use a variety of plots

par(mfrow = c(3, 1))
## First, drivers
plot(b1, type = "line", addtot = TRUE)
plot(b1, type = "stacked")
plot(b1, type = "stream")

Figure 6.1: Three drivers’ plots of a simulation of Bauer’s model

par(mfrow = c(3, 1))
## Next, genotypes
plot(b1, show = "genotypes", type = "line")
plot(b1, show = "genotypes", type = "stacked")
plot(b1, show = "genotypes", type = "stream")

Figure 6.2: Three genotypes’ plots of a simulation of Bauer’s model

In this case, probably the stream plots are most helpful. Note, however, that (in contrast to some figures in the literature showing models of clonal expansion) the stream plot (or the stacked plot) does not try to explicitly show parent-descendant relationships, which would hardly be realistically possible in these plots (although the plots of phylogenies in section 8 could be of help).

6.5.2 McFarland model with 5000 passengers and 70 drivers

set.seed(678)
nd <- 70  
np <- 5000 
s <- 0.1  
sp <- 1e-3 
spp <- -sp/(1 + sp)
mcf1 <- allFitnessEffects(noIntGenes = c(rep(s, nd), rep(spp, np)),
                          drvNames = seq.int(nd))
mcf1s <-  oncoSimulIndiv(mcf1,
                         model = "McFL", 
                         mu = 1e-7,
                         detectionProb = "default",
                         detectionSize = NA, 
                         detectionDrivers = NA,
                         sampleEvery = 0.025,
                         keepEvery = 8,
                         initSize = 2000,
                         finalTime = 4000,
                         onlyCancer = FALSE)
summary(mcf1s)
##   NumClones TotalPopSize LargestClone MaxNumDrivers MaxDriversLast
## 1       618         3207         2841             3              3
##   NumDriversLargestPop TotalPresentDrivers FinalTime NumIter
## 1                    3                   5   1961.05   81003
##   HittedWallTime HittedMaxTries    errorMF minDMratio minBMratio
## 1          FALSE          FALSE 0.01289961   1584.625   1972.387
##    OccurringDrivers
## 1 5, 24, 30, 57, 67

par(mfrow  = c(2, 1))
## I use thinData to make figures smaller and faster
plot(mcf1s, addtot = TRUE, lwdClone = 0.9, log = "", 
     thinData = TRUE, thinData.keep = 0.5)
plot(mcf1s, show = "drivers", type = "stacked",
     thinData = TRUE, thinData.keep = 0.3,
     legend.ncols = 2)

With the above output (where we see there are over 500 different genotypes) trying to represent the genotypes makes no sense.

6.5.3 McFarland model with 50,000 passengers and 70 drivers: clonal competition

The next is too slow (takes a couple of minutes in an i5 laptop) and too big to run in a vignette, because we keep track of over 4000 different clones (which leads to a result object of over 800 MB):

set.seed(123)
nd <- 70  
np <- 50000 
s <- 0.1  
sp <- 1e-4 ## as we have many more passengers
spp <- -sp/(1 + sp)
mcfL <- allFitnessEffects(noIntGenes = c(rep(s, nd), rep(spp, np)),
                          drvNames = seq.int(nd))
mcfLs <-  oncoSimulIndiv(mcfL,
                         model = "McFL", 
                         mu = 1e-7,
                         detectionSize = 1e8, 
                         detectionDrivers = 100,
                         sampleEvery = 0.02,
                         keepEvery = 2,
                         initSize = 1000,
                         finalTime = 2000,
                         onlyCancer = FALSE)

But you can access the pre-stored results and plot them (beware: this object has been trimmed by removing empty passenger rows in the Genotype matrix)

data(mcfLs)
plot(mcfLs, addtot = TRUE, lwdClone = 0.9, log = "",
     thinData = TRUE, thinData.keep = 0.3,
     plotDiversity = TRUE)

The argument plotDiversity = TRUE asks to show a small plot on top with Shannon’s diversity index.

summary(mcfLs)
##   NumClones TotalPopSize LargestClone MaxNumDrivers MaxDriversLast
## 1      4458         1718          253             3              3
##   NumDriversLargestPop TotalPresentDrivers FinalTime NumIter
## 1                    3                  70      2000  113759
##   HittedWallTime HittedMaxTries    errorMF minDMratio minBMratio
## 1          FALSE          FALSE 0.01921737   184.1019   199.6085
##   OccurringDrivers
## 1   13, 38, 40, 69
## number of passengers per clone
summary(colSums(mcfLs$Genotypes[-(1:70), ]))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   4.000   6.000   5.673   7.750  13.000

Note that we see clonal competition between clones with the same number of drivers (and with different drivers, of course). We will return to this (section 6.5.5).

A stacked plot might be better to show the extent of clonal competition (plotting takes some time —a stream plot reveals similar patterns and is also slower than the line plot). I will aggressively thin the data for this plot so it is faster and smaller (but we miss some of the fine grain, of course):

plot(mcfLs, type = "stacked", thinData = TRUE, 
     thinData.keep = 0.2,
     plotDiversity = TRUE,
     xlim = c(0, 1000))

6.5.4 Simulation with a conjunction example

We will use several of the previous examples. Most of them are in file examplesFitnessEffects, where they are stored inside a list, with named components (names the same as in the examples above):

data(examplesFitnessEffects)
names(examplesFitnessEffects)
##  [1] "cbn1"   "cbn2"   "smn1"   "xor1"   "fp3"    "fp4m"   "o3"    
##  [8] "ofe1"   "ofe2"   "foi1"   "sv"     "svB"    "svB1"   "sv2"   
## [15] "sm1"    "e2"     "E3A"    "em"     "fea"    "fbauer" "w"     
## [22] "pancr"

We will simulate using the simple CBN-like restrictions of section 3.4.4 with two different models.

data(examplesFitnessEffects)
evalAllGenotypes(examplesFitnessEffects$cbn1, order = FALSE)[1:10, ]
##    Genotype Fitness
## 1         a    1.10
## 2         b    1.10
## 3         c    0.10
## 4         d    1.10
## 5         e    1.10
## 6         g    0.10
## 7      a, b    1.21
## 8      a, c    0.11
## 9      a, d    1.21
## 10     a, e    1.21
sm <-  oncoSimulIndiv(examplesFitnessEffects$cbn1,
                      model = "McFL", 
                      mu = 5e-7,
                      detectionSize = 1e8, 
                      detectionDrivers = 2,
                      detectionProb = "default",
                      sampleEvery = 0.025,
                      keepEvery = 5,
                      initSize = 2000,
                      onlyCancer = TRUE)
summary(sm)
##   NumClones TotalPopSize LargestClone MaxNumDrivers MaxDriversLast
## 1         6         2721         2005             2              2
##   NumDriversLargestPop TotalPresentDrivers FinalTime NumIter
## 1                    2                   4  2256.725   90282
##   HittedWallTime HittedMaxTries    errorMF minDMratio minBMratio
## 1          FALSE          FALSE 0.01501927   309423.6   333333.3
##   OccurringDrivers
## 1       a, b, c, e

set.seed(1234)
evalAllGenotypes(examplesFitnessEffects$cbn1, order = FALSE, 
                 model = "Bozic")[1:10, ]
##    Genotype Death_rate
## 1         a       0.90
## 2         b       0.90
## 3         c       1.90
## 4         d       0.90
## 5         e       0.90
## 6         g       1.90
## 7      a, b       0.81
## 8      a, c       1.71
## 9      a, d       0.81
## 10     a, e       0.81
sb <-  oncoSimulIndiv(examplesFitnessEffects$cbn1,
                       model = "Bozic", 
                       mu = 5e-6,
                      detectionProb = "default",
                       detectionSize = 1e8, 
                       detectionDrivers = 4,
                       sampleEvery = 2,
                       initSize = 2000,
                       onlyCancer = TRUE)
summary(sb)
##   NumClones TotalPopSize LargestClone MaxNumDrivers MaxDriversLast
## 1         7        25636        23975             2              2
##   NumDriversLargestPop TotalPresentDrivers FinalTime NumIter
## 1                    1                   5        88      53
##   HittedWallTime HittedMaxTries errorMF minDMratio minBMratio
## 1          FALSE          FALSE      NA   33333.33   33333.33
##   OccurringDrivers
## 1    a, b, c, e, g

As usual, we will use several plots here.

## Show drivers, line plot
par(cex = 0.75, las = 1)
plot(sb,show = "drivers", type = "line", addtot = TRUE,
     plotDiversity = TRUE)

## Drivers, stacked
par(cex = 0.75, las = 1)
plot(sb,show = "drivers", type = "stacked", plotDiversity = TRUE)

## Drivers, stream
par(cex = 0.75, las = 1)
plot(sb,show = "drivers", type = "stream", plotDiversity = TRUE)

## Genotypes, line plot
par(cex = 0.75, las = 1)
plot(sb,show = "genotypes", type = "line", plotDiversity = TRUE)

## Genotypes, stacked
par(cex = 0.75, las = 1)
plot(sb,show = "genotypes", type = "stacked", plotDiversity = TRUE)

## Genotypes, stream
par(cex = 0.75, las = 1)
plot(sb,show = "genotypes", type = "stream", plotDiversity = TRUE)

The above illustrates again that different types of plots can be useful to reveal different patterns in the data. For instance, here, because of the huge relative frequency of one of the clones/genotypes, the stacked and stream plots do not reveal the other clones/genotypes as we cannot use a log-transformed y-axis, even if there are other clones/genotypes present.

6.5.5 Simulation with order effects and McFL model

(We use a somewhat large mutation rate than usual, so that the simulation runs quickly.)

set.seed(4321)
tmp <-  oncoSimulIndiv(examplesFitnessEffects[["o3"]],
                       model = "McFL", 
                       mu = 5e-5,
                       detectionSize = 1e8, 
                       detectionDrivers = 3,
                       sampleEvery = 0.025,
                       max.num.tries = 10,
                       keepEvery = 5,
                       initSize = 2000,
                       finalTime = 6000,
                       onlyCancer = FALSE) 

We show a stacked and a line plot of the drivers:

par(las = 1, cex = 0.85)
plot(tmp, addtot = TRUE, log = "", plotDiversity = TRUE,
     thinData = TRUE, thinData.keep = 0.2)

par(las = 1, cex = 0.85)
plot(tmp, type = "stacked", plotDiversity = TRUE, 
     ylim = c(0, 5500), legend.ncols = 4,
     thinData = TRUE, thinData.keep = 0.2)

In this example (and at least under Linux, with both GCC and clang —random number streams in C++, and thus simulations, can differ between combinations of operating system and compiler), we can see that the mutants with three drivers do not get established when we stop the simulation at time 6000. This is one case where the summary statistics about number of drivers says little of value, as fitness is very different for genotypes with the same number of mutations, and does not increase in a simple way with drivers:

evalAllGenotypes(examplesFitnessEffects[["o3"]], addwt = TRUE,
                 order = TRUE)
##     Genotype Fitness
## 1         WT    1.00
## 2          d    1.00
## 3          f    1.00
## 4          m    1.00
## 5      d > f    1.00
## 6      d > m    1.10
## 7      f > d    1.00
## 8      f > m    1.00
## 9      m > d    1.50
## 10     m > f    1.00
## 11 d > f > m    1.54
## 12 d > m > f    1.32
## 13 f > d > m    0.77
## 14 f > m > d    1.50
## 15 m > d > f    1.50
## 16 m > f > d    1.50

A few figures could help:

plot(tmp, show = "genotypes", ylim = c(0, 5500), legend.ncols = 3,
     thinData = TRUE, thinData.keep = 0.5)

(When reading the figure legends, recall that genotype \(x > y\ \_\ z\) is one where a mutation in “x” happened before a mutation in “y”, and there is also a mutation in “z” for which order does not matter. Here, there are no genes for which order does not matter and thus there is nothing after the “_“).

In this case, the clones with three drivers end up displacing those with two by the time we stop; moreover, notice how those with one driver never really grow to a large population size, so we basically go from a population with clones with zero drivers to a population made of clones with two or three drivers:

set.seed(15)
tmp <-  oncoSimulIndiv(examplesFitnessEffects[["o3"]],
                       model = "McFL", 
                       mu = 5e-5,
                       detectionSize = 1e8, 
                       detectionDrivers = 3,
                       sampleEvery = 0.025,
                       max.num.tries = 10,
                       keepEvery = 5,
                       initSize = 2000,
                       finalTime = 20000,
                       onlyCancer = FALSE,
                       extraTime = 1500)
tmp
## 
## Individual OncoSimul trajectory with call:
##  oncoSimulIndiv(fp = examplesFitnessEffects[["o3"]], model = "McFL", 
##     mu = 5e-05, detectionSize = 1e+08, detectionDrivers = 3, 
##     sampleEvery = 0.025, initSize = 2000, keepEvery = 5, extraTime = 1500, 
##     finalTime = 20000, onlyCancer = FALSE, max.num.tries = 10)
## 
##   NumClones TotalPopSize LargestClone MaxNumDrivers MaxDriversLast
## 1         8         4222         2758             3              3
##   NumDriversLargestPop TotalPresentDrivers FinalTime NumIter
## 1                    3                   3   6259.25  251948
##   HittedWallTime HittedMaxTries    errorMF minDMratio minBMratio
## 1          FALSE          FALSE 0.01441608   6161.251   6666.667
##   OccurringDrivers
## 1          d, f, m
## 
## Final population composition:
##       Genotype    N
## 1           _     0
## 2         d _     0
## 3     d > f _     0
## 4         f _     0
## 5         m _     0
## 6     m > d _  1464
## 7 m > d > f _  2758
## 8     m > f _     0

use a drivers plot:

par(las = 1, cex = 0.85)
plot(tmp, addtot = TRUE, log = "", plotDiversity = TRUE,
     thinData = TRUE, thinData.keep = 0.5)

par(las = 1, cex = 0.85)
plot(tmp, type = "stacked", plotDiversity = TRUE,
     legend.ncols = 4, ylim = c(0, 5200), xlim = c(3400, 5000),
     thinData = TRUE, thinData.keep = 0.5)

Now show the genotypes explicitly:

## Improve telling apart the most abundant 
## genotypes by sorting colors
## differently via breakSortColors
## Modify ncols of legend, so it is legible by not overlapping
## with plot
par(las = 1, cex = 0.85)
plot(tmp, show = "genotypes", breakSortColors = "distave",
     plotDiversity = TRUE, legend.ncols = 4,
     ylim = c(0, 5300), xlim = c(3400, 5000),
     thinData = TRUE, thinData.keep = 0.5)

As before, the argument plotDiversity = TRUE asks to show a small plot on top with Shannon’s diversity index. Here, as before, the quick clonal expansion of the clone with two drivers leads to a sudden drop in diversity (for a while, the population is made virtually of a single clone). Note, however, that compared to section 6.5.3, we are modeling here a scenario with very few genes, and correspondingly very few possible genotypes, and thus it is not strange that we observe very little diversity.

(We have used extraTime to continue the simulation well past the point of detection, here specified as three drivers. Instead of specifying extraTime we can set the detectionDrivers value to a number larger than the number of existing possible drivers, and the simulation will run until finalTime if onlyCancer = FALSE.)

12.2.1 sampleEvery, keepPhylog, and pruning

At each sampling time (where sampleEvery determines the time units between sampling times) the abundance of all the clones with number of cells \(>0\) is recorded. This is the structure that at the end of the run is converted into the pops.by.time matrix.

Now, some clones might arise from mutation between successive population samples but these clones might be extinct by the time we take a population sample. These clones do not appear in the pops.by.time matrix because, as we just said, they have 0 cells at the time of sampling. Of course, some of these clones might appear again later and reach a size larger than 0 at some posterior sampling time; it is at this time when this/these clone(s) will appear in the pops.by.time matrix. This pruning of clones with 0 cells can allow considerable savings in computing time (OncoSimulR needs to track the genotype of clones, their population sizes, their birth, death, and mutation rates, their next mutation time and the last time they were updated and thus it is important that we only loop over structures with information that is really needed).

However, we still need to track clones as clones, not simply as classes such as “number of mutated genes”. Therefore, very large genomes can represent a problem if they lead to the creation and tracking of many different clones (even if they have the same number of mutated genes), as we have seen, for instance, in section 2.3. In this case, programs that only keep track of numbers of mutated genes or of drivers, not individual clones, can of course achieve better speed.

What about the genealogy? If you ask OncoSimulR to keep track of the complete parent-child relationships (keepPhylog = TRUE), you might see in the genealogy clones that are not present in pops.by.time if these are clones that never had a population size larger than 0 at any sampling time. To give an example, suppose that we will take population samples at times 0, 1, and 2. Clone A, with a population size larger than 0 at time 1, gives rise at time 1.5 to clone B; clone B then gives rise to clone C at time 1.8. Finally, suppose that at time 2 only clone C is alive. In other words, when we carry out the update of the population with Algorithm 5 from Mather, Hasty, and Tsimring (2012), clones A and B have size 0. Now, at time 1 clones B and C did not yet exist, and clone B is never alive at times 1 or 2. Thus, clone B is not present in pops.by.time. But we cannot remove clone B from our genealogy if we want to reflect the complete genealogy of C. Thus, pops.by.time will show only clones A and C (not B) but the complete genealogy will show clones A, B, C (and will show that B appeared from A at time 1.5 and C appeared from B at time 1.8). Since function plotClonePhylog offers a lot of flexibility with respect to what clones to show depending on their population sizes at different times, you can prevent being shown B, but its existence is there should you need it (see also 2.3.5).

13.1.1 Simulating progression in several subjects

To simulate the progression in a bunch of subjects (we will use only four, so as not to fill the vignette with plots) we can do, with the same settings as above:

p1 <- oncoSimulPop(4, p1101, mc.cores = 2)
par(mfrow = c(2, 2))
plot(p1, ask = FALSE)

We can also use stream and stacked plots, though they might not be as useful in this case. For the sake of keeping the vignette small, these are commented out.

par(mfrow = c(2, 2))
plot(p1, type = "stream", ask = FALSE)

par(mfrow = c(2, 2))
plot(p1, type = "stacked", ask = FALSE)