Key features of OncoSimulR

As mentioned above, OncoSimulR is now a very general package for forward genetic simulation, with applicability well beyond tumor progression. This is a summary of some of its key features:

• You can specify arbitrary interactions between genes, with arbitrary fitness effects, with explicit support for:

  • Restrictions in the accumulations of mutations, as specified by Oncogenetic Trees (OTs), Conjunctive Bayesian Networks (CBNs), semimonotone progression networks, and XOR relationships.

  • Epistatic interactions including, but not limited to, synthetic viability and synthetic lethality.

  • Order effects.

• You can add passenger mutations.

• You can add mutator/antimutator effects.

• Fitness and mutation rates can be gene-specific.

• You can add arbitrary numbers of non-interacting genes with arbitrary fitness effects.

• you can allow for deviations from the OT, CBN, semimonotone, and XOR models, specifying a penalty for such deviations (the s_h parameter).

• You can conduct multiple simulations, and sample from them with different temporal schemes and using both whole tumor or single cell sampling.

• You can stop the simulations using a flexible combination of conditions: final time, number of drivers, population size, fixation of certain genotypes, and a stochastic stopping mechanism that depends on population size.

• Right now, three different models are available, two that lead to exponential growth, one of them loosely based on Bozic et al. (2010), and another that leads to logistic-like growth, based on McFarland et al. (2013).

• You can use large numbers of genes (e.g., see an example of 50000 in section @ref(mcf50070)).

• Simulations are generally very fast: I use C++ to implement the BNB algorithm (see sections @ref(bnbmutation) and @ref(bnbdensdep) for more detailed comments on the usage of this algorithm).

• You can obtain the true sequence of events and the phylogenetic relationships between clones (see section @ref(meaningclone) for the details of what we mean by “clone”).

• You can generate random fitness landscapes (under the House of Cards, Rough Mount Fuji, or additive models, or combinations of the former and under the NK model) and use those landscapes as input to the simulation functions.

• You can plot fitness landscapes.

• You can obtain statistics of evolutionary predictability from the simulations.

• You can now also use simulations with frequency-dependent fitness: fitness (birth rate) is not fixed for a genotype, but can be a function of the frequecies of the clones (see section @ref(fdf)). We can therefore use OncoSimulR to examine, via simulations, results from game theory and adaptive dynamics and study complex scenarios that are not amenable to analytical solutions. More generally, we can model predation and parasitism, cooperation and mutualism, and commensalism.

• It is possible to start the simulation with arbitrary initial composition (section @ref(minitmut)) and to simulate multiple species (section @ref(multispecies)). You can thus run simulations that involve both ecological and evolutionary processes involving inter-species relationships plus genetic restrictions in evolution.

• It is possible to simulate many different therapeutic interventions. Section @ref(timefdf) shows examples of interventions where certain genotypes change fitness (because of chemotherapy) at specified times. More generally, since fitness (birth rates) can be made a function of total populations sizes and/or frequencies (see section @ref(fdf)), many different arbitrary intervention schemes can be simulated. Possible models are, of course, not limited to cancer chemotherapy, but could include antibiotic treatment of bacteria, antiviral therapy, etc.

The table below, modified from the table at the Genetics Simulation Resources (GSR) page, provides a summary of the key features of OncoSimulR. (An explanation of the meaning of terms specific to the GSR table is available from https://popmodels.cancercontrol.cancer.gov/gsr/search/ or from the Genetics Simulation Resources table itself, by moving the mouse over each term).

Further details about the original motivation for wanting to simulate data this way in the context of tumor progression can be found in Diaz-Uriarte (2015), where additional comments about model parameters and caveats are discussed.

Are there similar programs? The Java program by Reiter et al. (2013), TTP, offers somewhat similar functionality to the previous version of OncoSimulR, but it is restricted to at most four drivers (whereas v.1 of OncoSimulR allowed for up to 64), you cannot use arbitrary CBNs or OTs (or XORs or semimonotone graphs) to specify restrictions, there is no allowance for passengers, and a single type of model (a discrete time Galton-Watson process) is implemented. The current functionality of OncoSimulR goes well beyond the the previous version (and, thus, also the TPT of Reiter et al. (2013)). We now allow you to specify all types of fitness effects in other general forward genetic simulators such as FFPopSim (Zanini & Neher, 2012), and some that, to our knowledge (e.g., order effects) are not available from any genetics simulator. In addition, the “Lego system” to flexibly combine different fitness specifications is also unique; by “Lego system” I mean that we can combine different pieces and blocks, similarly to what we do with Lego bricks. (I find this an intuitive and very graphical analogy, which I have copied from Hothorn et al. (2006) and Hothorn et al. (2008)). In a nutshell, salient features of OncoSimulR compared to other simulators are the unparalleled flexibility to specify fitness and mutator effects, with modules and order effects as particularly unique, and the options for sampling and stopping the simulations, particularly convenient in cancer evolution models. Also unique in this type of software is the addition of functions for simulating fitness landscapes and assessing evolutionary predictability.

What kinds of questions is OncoSimulR suited for?

OncoSimulR can be used to address questions that span from the effect of mutator genes in cancer to the interplay between fitness landscapes and mutation rates. The main types of questions that OncoSimulR can help address involve combinations of:

• Simulating asexual evolution (the oncoSimul* functions) where:

  • Fitness is:
    • A function of specific epistatic effects between genes
    • A function of order effects
    • A function of epistatic effects specified using DAGs/posets where these DAGs/posets:
      • Are user-specified
      • Generated randomly (simOGraph)
    • Any mapping between genotypes and fitness where this mapping is:
      • User-specified
      • Generated randomly from families of random fitness landscapes (rfitness)
    • A function of the frequency of other genotypes (i.e., frequency-dependent fitness), such as in adaptive dynamics (see section @ref(fdf) for more details). This also allows you to model competition, cooperation and mutualism, parasitism and predation, and commensalism between clones.
  • Mutation rates can:
    • Vary between genes
    • Be affected by other genes

• Examining times to evolutionarily or biomedically relevant events (fixation of genotypes, reaching a minimal size, acquiring a minimal number of driver genes, etc —specified with the stopping conditions to the oncoSimul* functions).

• Using different sampling schemes (samplePop) that are related to:

  • Assessing genotypes from single-cell vs. whole tumor (or whole population) with the typeSample argument
  • Genotyping error (propError argument)
  • Timing of samples (timeSample argument)
  • … and assessing the consequences of those on the observed genotypes and their diversity (sampledGenotypes) and any other inferences that depend on the observational process.
  • (OncoSimulR returns the abundances of all genotypes at each of the sampling points, so you are not restricted by what the samplePop function provides.)

• Tracking the genealogical relationships of clones (plotClonePhylog) and assessing evolutionary predictability (LOD, POM).

Some specific questions that you can address with the help of OncoSimulR are discussed in section @ref(whatfor).

A quick overview of the main functions and their relationships is shown in Figure @ref(fig:frelats), where we use italics for the type/class of R object and courier font for the name of the functions.

Relationships between the main functions in OncoSimulR.

Examples of questions that can be addressed with OncoSimulR

Most of the examples in the rest of this vignette, starting with those in @ref(quickexample), focus on the mechanics. Here, we will illustrate some problems in cancer genomics and evolutionary genetics where OncoSimulR could be of help. This section does not try to provide an answer to any of these questions (those would be full papers by themselves). Instead, this section simply tries to illustrate some kinds of questions where you can use OncoSimulR; of course, the possible uses of OncoSimulR are only limited by your ingenuity. Here, I will only use short snippets of working code as we are limited by time of execution; for real work you would want to use many more scenarios and many more simulations, you would use appropriate statistical methods to compare the output of runs, etc, etc, etc.

## Load the package
library(OncoSimulR)
## This is package OncoSimulR. If you are running it on an aarch64 (arm64) platform with a MacOS note that the package fails some tests in that platform + OS that I have no way of debugging. Please read file README_tests_kjohnson3_aarch64-apple-darwin20.txt in the tests directory. As of 2024-10-09, this platform is unsupported until we can properly debug it.

## 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),
                   onlyCancer = TRUE,
                   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).

## 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)
## Using old version of fitnessEffects. Transforming fitnessEffects
##             to last version.
##   Genotype Birth 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

## Set a small initSize, as o.w. unlikely to pass the valley
initS <- 10
## The number of replicates is tiny, for the sake of speed
## of creation of the vignette. Even fewer in Windows, since we run on a single
## core

if(.Platform$OS.type == "windows") {
    nruns <- 4
} else {
    nruns <- 10
}

od_sim <- oncoSimulPop(nruns, 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

## 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 and
## call your own copy of MAGELLAN's binary

## to_Magellan(rl, file = "rl1.txt") ## (Commented out here to avoid writing files)

## or use the binary copy provided with OncoSimulR
## see also below.

Magellan_stats(rl) ## (Commented out here to avoid writing files)
##      ngeno     npeaks     nsinks      gamma     gamma.        r.s 
##     64.000      2.000      1.000      0.769      0.854      0.372 
##    nchains     nsteps       nori      depth       magn       sign 
##      1.000      5.000      4.000      2.000      0.863      0.129 
##      rsign       f.1.       X.2.      f.3..     mode_f     outD_m 
##      0.008      0.916      0.010      0.074      1.000      0.346 
##     outD_v    steps_m    reach_m     fitG_m      opt_i mProbOpt_0 
##      1.587      3.052     13.597     32.452     60.000      0.128 
##    opt_i.1 mProbOpt_1 
##     63.000      0.872

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

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

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)

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

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     Birth
## 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
## Warning: ggrepel: 12 unlabeled data points (too many overlaps).
## Consider increasing max.overlaps

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

Trade-offs and what is OncoSimulR not well suited for

OncoSimulR is designed for complex fitness specifications and selection scenarios and uses forward-time simulations; the types of questions where OncoSimulR can be of help are discussed in sections @ref(generalwhatfor) and @ref(whatfor) and running time and space consumption of OncoSimulR are addressed in section @ref(timings). You should be aware that coalescent simulations, sometimes also called backward-time simulations, are much more efficient for simulating neutral data as well as some special selection scenarios (Carvajal-Rodriguez, 2010; Hoban et al., 2011; Yuan et al., 2012).

In addition, since OncoSimulR allows you to specify fitness with arbitrary epistatic and order effects, as well as mutator effects, you need to learn the syntax of how to specify those effects and you might be paying a performance penalty if your scenario does not require this complexity. For instance, in the model of Beerenwinkel, Antal, et al. (2007), the fitness of a genotype depends only on the total number of drivers mutated, but not on which drivers are mutated (and, thus, not on the epistatic interactions nor the order of accumulation of the drivers). This means that the syntax for specifying that model could probably be a lot simpler (e.g., specify s per driver).

But it also means that code written for just that case could probably run much faster. First, because fitness evaluation is easier. Second, and possibly much more important, because what we need to keep track of leads to much simpler and economic structures: we do not need to keep track of clones (where two cells are regarded as different clones if they differ anywhere in their genotype), but only of clone types or clone classes as defined by the number of mutated drivers, and keeping track of clones can be expensive —see sections @ref(timings) and @ref(trackindivs).

So for those cases where you do not need the full flexibility of OncoSimulR, special purpose software might be easier to use and faster to run. Of course, for some types of problems this special purpose software might not be available, though.

Random fitness landscapes, clonal competition, predictability, and the strong selection weak mutation (SSWM) regime

Many studies about evolutionary predictability (among other topics) focus on the strong selection, weak mutation regime, SSWM (Gillespie, 1984; Orr, 2002) (see overview in Krug (2019)). In this regime, mutations are rare (much smaller than the mutation rate times the population size) and selection is strong (much larger than 1/population size), so that the population consists of a single clone most of the time, and evolution proceeds by complete, successive clonal expansions of advantageous mutations.

We can easily simulate variations around these scenarios with OncoSimulR, moving away from the SSWM by increasing the population size, or changing the size of the fitness differences.

The examples below, not run for the sake of speed, play with population size and fitness differences. To make sure we use a similar fitness landscape, we use the same simulated fitness landscape, scaled differently, so that the differences in fitness between mutants are increased or decreased while keeping their ranking identical (and, thus, having the same set of accessible and inaccessible genotypes and paths over the landscape).

If you run the code, you will see that as we increase population size we move further away from the SSWM: the population is no longer composed of a single clone most of the time.

Before running the examples, and to show the effects quantitatively, we define a simple wrapper to compute a few statistics.

## oncoSimul object  -> measures of clonal interference
##    they are not averaged over time. One value for sampled time
clonal_interf_per_time <- function(x) {
    x <- x$pops.by.time
    y <- x[, -1, drop = FALSE]
    shannon <- apply(y, 1, OncoSimulR:::shannonI)
    tot <- rowSums(y)
    half_tot <- tot * 0.5
    five_p_tot <- tot * 0.05
    freq_most_freq <- apply(y/tot, 1, max)
    single_more_half <- rowSums(y > half_tot)
    ## whether more than 1 clone with more than 5% pop.
    how_many_gt_5p <- rowSums(y > five_p_tot)
    several_gt_5p <- (how_many_gt_5p > 1)
    return(cbind(shannon, ## Diversity of clones
                 freq_most_freq, ## Frequency of the most freq. clone
                 single_more_half, ## Any clone with a frequency > 50%?
                 several_gt_5p, ## Are there more than 1 clones with
                                ## frequency > 5%?
                 how_many_gt_5p ## How many clones are there with
                                ## frequency > 5%
                 ))
}

set.seed(1)
r7b <- rfitness(7, scale = c(1.2, 0, 1))

## Large pop sizes: clonal interference
(sr7b <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b),
                       model = "McFL",
                       mu = 1e-6,
                       onlyCancer = FALSE,
                       finalTime = 400,
                       initSize = 1e7,
                       keepEvery = 4,
                       detectionSize = 1e10))

plot(sr7b, show = "genotypes")

colMeans(clonal_interf_per_time(sr7b))

## Small pop sizes: a single clone most of the time
(sr7c <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b),
                       model = "McFL",
                       mu = 1e-6,
                       onlyCancer = FALSE,
                       finalTime = 60000,
                       initSize = 1e3,
                       keepEvery = 4,
                       detectionSize = 1e10))

plot(sr7c, show = "genotypes")

colMeans(clonal_interf_per_time(sr7c))



## Even smaller fitness differences, but large pop. sizes
set.seed(1); r7b2 <- rfitness(7, scale = c(1.05, 0, 1))

(sr7b2 <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b2),
                       model = "McFL",
                       mu = 1e-6,
                       onlyCancer = FALSE,
                       finalTime = 3500,
                       initSize = 1e7,
                       keepEvery = 4,
                       detectionSize = 1e10))
sr7b2
plot(sr7b2, show = "genotypes")
colMeans(clonal_interf_per_time(sr7b2))


## Increase pop size further
(sr7b3 <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b2),
                       model = "McFL",
                       mu = 1e-6,
                       onlyCancer = FALSE,
                       finalTime = 1500,
                       initSize = 1e8,
                       keepEvery = 4,
                       detectionSize = 1e10))
sr7b3
plot(sr7b3, show = "genotypes")
colMeans(clonal_interf_per_time(sr7b3))

Steps for using OncoSimulR

Using this package will often involve the following steps:

1. Specify fitness effects: sections @ref(specfit) and @ref(litex).

2. Simulate cancer progression: section @ref(simul). You can simulate for a single individual or subject or for a set of subjects. You will need to:

   • Decide on a model. This basically amounts to choosing a model with exponential growth (“Exp” or “Bozic”) or a model with carrying capacity (“McFL”). If exponential growth, you can choose whether the the effects of mutations operate on the death rate (“Bozic”) or the birth rate (“Exp”)¹.

   • Specify other parameters of the simulation. In particular, decide when to stop the simulation, mutation rates, etc.

   Of course, at least for initial playing around, you can use the defaults.

3. Sample from the simulated data and do something with those simulated data (e.g., fit an OT model to them, examine diversity or time until cancer, etc). Most of what you do with the data, however, is outside the scope of this package and this vignette.

Before anything else, let us load the package in case it was not yet loaded. We also explicitly load graph and igraph for the vignette to work (you do not need that for your usual interactive work). And I set the default color for vertices in igraph.

library(OncoSimulR)
library(graph)
library(igraph)
igraph_options(vertex.color = "SkyBlue2")

To be explicit, what version are we running?

packageVersion("OncoSimulR")
## [1] '4.9.0'

Two quick examples of fitness specifications

Following @ref(steps) we will run two very minimal examples. First a model with a few genes and epistasis:

## 1. Fitness effects: here we specify an
##    epistatic model with modules.
sa <- 0.1
sb <- -0.2
sab <- 0.25
sac <- -0.1
sbc <- 0.25
sv2 <- allFitnessEffects(epistasis = c("-A : B" = sb,
                                       "A : -B" = sa,
                                       "A : C" = sac,
                                       "A:B" = sab,
                                       "-A:B:C" = sbc),
                         geneToModule = c(
                             "A" = "a1, a2",
                             "B" = "b",
                             "C" = "c"),
                         drvNames = c("a1", "a2", "b", "c"))
evalAllGenotypes(sv2, addwt = TRUE)
##        Genotype Birth
## 1            WT 1.000
## 2            a1 1.100
## 3            a2 1.100
## 4             b 0.800
## 5             c 1.000
## 6        a1, a2 1.100
## 7         a1, b 1.250
## 8         a1, c 0.990
## 9         a2, b 1.250
## 10        a2, c 0.990
## 11         b, c 1.000
## 12    a1, a2, b 1.250
## 13    a1, a2, c 0.990
## 14     a1, b, c 1.125
## 15     a2, b, c 1.125
## 16 a1, a2, b, c 1.125

## 2. Simulate the data. Here we use the "McFL" model and set
##    explicitly parameters for mutation rate, initial size, size
##    of the population that will end the simulations, etc

RNGkind("Mersenne-Twister")
set.seed(983)
ep1 <- oncoSimulIndiv(sv2, model = "McFL",
                      mu = 5e-6,
                      sampleEvery = 0.025,
                      keepEvery = 0.5,
                      initSize = 2000,
                      finalTime = 3000,
                      onlyCancer = FALSE)

## 3. We will not analyze those data any further. We will only plot
## them.  For the sake of a small plot, we thin the data.
plot(ep1, show = "drivers", xlim = c(0, 1500),
     thinData = TRUE, thinData.keep = 0.5)

Plot of drivers of an epistasis simulation.

As a second example, we will use a model where we specify restrictions in the order of accumulation of mutations using a DAG with the pancreatic cancer poset in Gerstung, Eriksson, et al. (2011) (see more details in section @ref(pancreas)):

## 1. Fitness effects:
pancr <- allFitnessEffects(
    data.frame(parent = c("Root", rep("KRAS", 4),
                   "SMAD4", "CDNK2A",
                   "TP53", "TP53", "MLL3"),
               child = c("KRAS","SMAD4", "CDNK2A",
                   "TP53", "MLL3",
                   rep("PXDN", 3), rep("TGFBR2", 2)),
               s = 0.1,
               sh = -0.9,
               typeDep = "MN"),
    drvNames = c("KRAS", "SMAD4", "CDNK2A", "TP53",
                 "MLL3", "TGFBR2", "PXDN"))

## Plot the DAG of the fitnessEffects object
plot(pancr)

Plot of DAG corresponding to fitnessEffects object.

## 2. Simulate from it. We change several possible options.

set.seed(1) ## Fix the seed, so we can repeat it
## We set a small finalTime to speed up the vignette


ep2 <- oncoSimulIndiv(pancr, model = "McFL",
                     mu = 1e-6,
                     sampleEvery = 0.02,
                     keepEvery = 1,
                     initSize = 1000,
                     finalTime = 20000,
                     detectionDrivers = 3,
                     onlyCancer = FALSE)

## 3. What genotypes and drivers we get? And play with limits
##    to show only parts of the data. We also aggressively thin
##    the data.
par(cex = 0.7)
plot(ep2, show = "genotypes", xlim = c(500, 1800),
     ylim = c(0, 2400),
     thinData = TRUE, thinData.keep = 0.3)

Plot of genotypes of a simulation from a DAG.

The rest of this vignette explores all of those functions and arguments in much more detail.

Citing OncoSimulR and other documentation

In R, you can do

citation("OncoSimulR")
## If you use OncoSimulR, please cite the OncoSimulR
## Bioinformatics paper.  OncoSimulR has been used in three
## large comparative studies of methods to infer restrictions
## in the order of accumulation of mutations (cancer
## progression models) published in PLoS Computational Biology,
## Bioinformatics and BMC Bioinformatics; you might want to
## cite those too, if appropriate, such as when referring to
## using evolutionary simulations to assess oncogenetic
## tree/cancer progression methods performance.
## 
##   R Diaz-Uriarte. OncoSimulR: genetic simulation with
##   arbitrary epistasis and mutator genes in asexual
##   populations.  2017. Bioinformatics, 33, 1898--1899.
##   https://doi.org/10.1093/bioinformatics/btx077.
## 
##   R Diaz-Uriarte and C. Vasallo. Every which way? On
##   predicting tumor evolution using cancer progression models
##   2019 PLoS Computational Biology
##   https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007246
## 
##   R Diaz-Uriarte. Cancer progression models and fitness
##   landscapes: a many-to-many relationship 2017
##   Bioinformatics.
##   https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btx663/
## 
##   R Diaz-Uriarte. Identifying restrictions in the order of
##   accumulation of mutations during tumor progression:
##   effects of passengers, evolutionary models, and sampling
##   2015. BMC Bioinformatics, 16(41).
## 
## To see these entries in BibTeX format, use
## 'print(<citation>, bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.

which will tell you how to cite the package. Please, do cite the Bionformatics paper if you use the package in publications.

This is the URL for the Bioinformatics paper: https://doi.org/10.1093/bioinformatics/btx077 (there is also an early preprint at bioRxiv, but it should now point to the Bioinformatics paper).

Testing, code coverage, and other examples

OncoSimulR includes more than 2000 tests that are run at every check cycle. These tests provide a code coverage of more than 90% including both the C++ and R code. Another set of over 500 long-running (several hours) tests can be run on demand (see directory ‘/tests/manual’). In addition to serving as test cases, some of that code also provides further examples of usage.

Versions

In this vignette and the documentation I often refer to version 1 (v.1) and version 2 of OncoSimulR. Version 1 is the version available up to, and including, BioConductor v. 3.1. Version 2 of OncoSimulR is available starting from BioConductor 3.2 (and, of course, available too from development versions of BioC). So, if you are using the current stable or development version of BioConductor, or you grab the sources from GitHub (https://github.com/rdiaz02/OncoSimul) you are using what we call version 2. The functionality of version 1 has been removed.

Version 3 (for BioConductor 3.13) made frequency dependent fitness available in the stable version.

Version 4 (BioConductor 3.16) introduces interventions and the possibility to specify, separately, birth and death (including frequency dependence).

Exp and McFL with “detectionProb” and pancreas example

To get familiar with some of they factors that affect time and size, we will use the fitness specification from section @ref(quickexample), with the detectionProb stopping mechanism (see @ref(detectprob)). We will use the two main growth models (exponential and McFarland). Each model will be run with two settings of keepEvery. With keepEvery = 1 (runs exp1 and mc1), population samples are stored at time intervals of 1 (even if most of the clones in those samples later become extinct). With keepEvery = NA (runs exp2 and mc2) no intermediate population samples are stored, so clones that become extinct at any sampling period are pruned and only the existing clones at the end of the simulation are returned (see details in @ref(prune)).

Will run 100 simulations. The results I show are for a laptop with an 8-core Intel Xeon E3-1505M CPU, running Debian GNU/Linux (the results from these benchmarks are available as data(benchmark_1)).

## Specify fitness
pancr <- allFitnessEffects(
    data.frame(parent = c("Root", rep("KRAS", 4),
                   "SMAD4", "CDNK2A",
                   "TP53", "TP53", "MLL3"),
               child = c("KRAS","SMAD4", "CDNK2A",
                   "TP53", "MLL3",
                   rep("PXDN", 3), rep("TGFBR2", 2)),
               s = 0.1,
               sh = -0.9,
               typeDep = "MN"),
    drvNames = c("KRAS", "SMAD4", "CDNK2A", "TP53",
                 "MLL3", "TGFBR2", "PXDN"))

Nindiv <- 100 ## Number of simulations run.
              ## Increase this number to decrease sampling variation

## keepEvery = 1
t_exp1 <- system.time(
    exp1 <- oncoSimulPop(Nindiv, pancr,
                         onlyCancer = TRUE,
                         detectionProb = "default",
                         detectionSize = NA,
                         detectionDrivers = NA,
                            finalTime = NA,
                            keepEvery = 1,
                            model = "Exp",
                            mc.cores = 1))["elapsed"]/Nindiv


t_mc1 <- system.time(
    mc1 <- oncoSimulPop(Nindiv, pancr,
                        onlyCancer = TRUE,
                        detectionProb = "default",
                        detectionSize = NA,
                        detectionDrivers = NA,
                           finalTime = NA,
                           keepEvery = 1,
                           model = "McFL",
                           mc.cores = 1))["elapsed"]/Nindiv

## keepEvery = NA
t_exp2 <- system.time(
    exp2 <- oncoSimulPop(Nindiv, pancr,
                         onlyCancer = TRUE,
                         detectionProb = "default",
                         detectionSize = NA,
                            detectionDrivers = NA,
                            finalTime = NA,
                            keepEvery = NA,
                            model = "Exp",
                            mc.cores = 1))["elapsed"]/Nindiv


t_mc2 <- system.time(
    mc2 <- oncoSimulPop(Nindiv, pancr,
                        onlyCancer = TRUE,
                        detectionProb = "default",
                        detectionSize = NA,
                           detectionDrivers = NA,
                           finalTime = NA,
                           keepEvery = NA,
                           model = "McFL",
                           mc.cores = 1))["elapsed"]/Nindiv

We can obtain times, sizes of objects, and summaries of numbers of clones, iterations, and final times doing, for instance:

cat("\n\n\n t_exp1 = ", t_exp1, "\n")
object.size(exp1)/(Nindiv * 1024^2)
cat("\n\n")
summary(unlist(lapply(exp1, "[[", "NumClones")))
summary(unlist(lapply(exp1, "[[", "NumIter")))
summary(unlist(lapply(exp1, "[[", "FinalTime")))
summary(unlist(lapply(exp1, "[[", "TotalPopSize")))

The above runs yield the following:

The above table shows that a naive comparison (looking simply at execution time) might conclude that the McFL model is much, much slower than the Exp model. But that is not the complete story: using the detectionProb stopping mechanism (see @ref(detectprob)) will lead to stopping the simulations very quickly in the exponential model because as soon as a clone with fitness  > 1 appears it starts growing exponentially. In fact, we can see that the number of iterations and the final time are much smaller in the Exp than in the McFL model. We will elaborate on this point below (section @ref(common1)), when we discuss the setting for checkSizePEvery (here left at its default value of 20): 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.

This table also shows that the keepEvery = NA setting, which was in effect in simulations exp2 and mc2, can make a difference especially for the McFL models, as seen by the median number of clones and the size of the returned object. Models exp2 and mc2 do not store any intermediate population samples so clones that become extinct at any sampling period are pruned and only the existing clones at the end of the simulation are returned. In contrast, models exp1 and mc1 store population samples at time intervals of 1 (keepEvery = 1), even if many of those clones eventually become extinct. We will return to this issue below as execution time and object size depend strongly on the number of clones tracked.

We can run the exponential model again modifying the arguments of the detectionProb mechanism; in two of the models below (exp3 and exp4) no detection can take place unless populations are at least 100 times larger than the initial population size, and probability of detection is 0.1 with a population size 1,000 times larger than the initial one (PDBaseline = 5e4, n2 = 5e5). In the other two models (exp5 and exp6), no detection can take place unless populations are at least 1,000 times larger than the initial population size, and probability of detection is 0.1 with a population size 100,000 times larger than the initial one (PDBaseline = 5e5, n2 = 5e7)². In runs exp3 and exp5 we set keepEvery = 1 and in runs exp4 and exp6 we set keepEvery = NA.

t_exp3 <- system.time(
    exp3 <- oncoSimulPop(Nindiv, pancr,
                         onlyCancer = TRUE,
                         detectionProb = c(PDBaseline = 5e4,
                                           p2 = 0.1, n2 = 5e5,
                                           checkSizePEvery = 20),
                         detectionSize = NA,
                         detectionDrivers = NA,
                            finalTime = NA,
                            keepEvery = 1,
                            model = "Exp",
                            mc.cores = 1))["elapsed"]/Nindiv

t_exp4 <- system.time(
    exp4 <- oncoSimulPop(Nindiv, pancr,
                         onlyCancer = TRUE,
                         detectionProb = c(PDBaseline = 5e4,
                                           p2 = 0.1, n2 = 5e5,
                                           checkSizePEvery = 20),
                         detectionSize = NA,
                         detectionDrivers = NA,
                            finalTime = NA,
                            keepEvery = NA,
                            model = "Exp",
                            mc.cores = 1))["elapsed"]/Nindiv



t_exp5 <- system.time(
    exp5 <- oncoSimulPop(Nindiv, pancr,
                         onlyCancer = TRUE,
                         detectionProb = c(PDBaseline = 5e5,
                                           p2 = 0.1, n2 = 5e7),
                         detectionSize = NA,
                         detectionDrivers = NA,
                            finalTime = NA,
                            keepEvery = 1,
                            model = "Exp",
                            mc.cores = 1))["elapsed"]/Nindiv

t_exp6 <- system.time(
    exp6 <- oncoSimulPop(Nindiv, pancr,
                         onlyCancer = TRUE,
                         detectionProb = c(PDBaseline = 5e5,
                                           p2 = 0.1, n2 = 5e7),
                         detectionSize = NA,
                         detectionDrivers = NA,
                            finalTime = NA,
                            keepEvery = NA,
                            model = "Exp",
                            mc.cores = 1))["elapsed"]/Nindiv

As above, keepEvery = NA (in exp4 and exp6) leads to much smaller object sizes and slightly smaller numbers of clones and execution times. Changing the exiting conditions (by changing detectionProb arguments) leads to large increases in number of iterations (in this case by factors of about 15x to 25x) and a corresponding increase in execution time as well as much larger population sizes (in some cases  > 10¹⁰).

In some of the runs of exp5 and exp6 we get the (recoverable) exception message from the C++ code: Recoverable exception ti set to DBL_MIN. Rerunning, which is related to those simulations reaching total population sizes  > 10¹⁰; we return to this below (section @ref(popgtzx)). You might also wonder why total and median population sizes are so large in these two runs, given the exiting conditions. One of the reasons is that we are using the default checkSizePEvery = 20, so the interval between successive checks of the exiting condition is large; this is discussed at greater length in section @ref(common1).

All the runs above used the default value onlyCancer = TRUE. This means that simulations will be repeated until the exiting conditions are reached (see details in section @ref(endsimul)) and, therefore, any simulation that ends up in extinction will be repeated. This setting can thus have a large effect on the exponential models, because when the initial population size is not very large and we start from the wildtype, it is not uncommon for simulations to become extinct (when birth and death rates are equal and the population size is small, it is easy to reach extinction before a mutation in a gene that increases fitness occurs). But this is rarely the case in the McFarland model (unless we use really tiny initial population sizes) because of the dependency of death rate on total population size (see section @ref(mcfl)).

The number of attempts until cancer was reached in the above models is shown in Table @ref(tab:bench1c) (the values can be obtained from any of the above runs doing, for instance, median(unlist(lapply(exp1, function(x) x$other$attemptsUsed))) ):

The McFL models finish in a single attempt. The exponential model simulations where we can exit with small population sizes (exp1, exp2) need many fewer attempts to reach cancer than those where large population sizes are required (exp3 to exp6). There is no relevant different among those last four, which is what we would expect: a population that has already reached a size of 50,000 cells from an initial population size of 500 is obviously a growing population where there is at least one mutant with positive fitness; thus, it unlikely to go extinct and therefore having to grow up to at least 500,000 will not significantly increase the risk of extinction.

We will now rerun all of the above models with argument onlyCancer = FALSE. The results are shown in Table @ref(tab:timing3) (note that the differences between this table and Table @ref(tab:bench1) for the McFL models are due only to sampling variation).

Now most simulations under the exponential model end up in extinction, as seen by the median population size of 0 (but not all, as the mean and max. population size are clearly away from zero). Consequently, simulations under the exponential model are now faster (and the size of the average returned object is smaller). Of course, whether one should run simulations with onlyCancer = TRUE or onlyCancer = FALSE will depend on the question being asked (see, for example, section @ref(exbauer) for a question where we will naturally want to use onlyCancer = FALSE).

To make it easier to compare results with those of the next section, Table @ref(tab:allr1bck) shows all the runs so far.

Several “common use cases” runs

Let us now execute some simulations under more usual conditions. We will use seven different fitness specifications: the pancreas example, two random fitness landscapes, and four sets of independent genes (200 to 4000 genes) with fitness effects randomly drawn from exponential distributions:

pancr <- allFitnessEffects(
    data.frame(parent = c("Root", rep("KRAS", 4),
                   "SMAD4", "CDNK2A",
                   "TP53", "TP53", "MLL3"),
               child = c("KRAS","SMAD4", "CDNK2A",
                   "TP53", "MLL3",
                   rep("PXDN", 3), rep("TGFBR2", 2)),
               s = 0.1,
               sh = -0.9,
               typeDep = "MN"),
    drvNames = c("KRAS", "SMAD4", "CDNK2A", "TP53",
                 "MLL3", "TGFBR2", "PXDN"))


## Random fitness landscape with 6 genes
## At least 50 accessible genotypes
rfl6 <- rfitness(6, min_accessible_genotypes = 50)
attributes(rfl6)$accessible_genotypes ## How many accessible
rf6 <- allFitnessEffects(genotFitness = rfl6)


## Random fitness landscape with 12 genes
## At least 200 accessible genotypes
rfl12 <- rfitness(12, min_accessible_genotypes = 200)
attributes(rfl12)$accessible_genotypes ## How many accessible
rf12 <- allFitnessEffects(genotFitness = rfl12)




## Independent genes; positive fitness from exponential distribution
## with mean around 0.1, and negative from exponential with mean
## around -0.02. Half of genes positive fitness effects, half
## negative.

ng <- 200 re_200 <- allFitnessEffects(noIntGenes = c(rexp(ng/2, 10),
-rexp(ng/2, 50)))

ng <- 500
re_500 <- allFitnessEffects(noIntGenes = c(rexp(ng/2, 10),
                                           -rexp(ng/2, 50)))

ng <- 2000
re_2000 <- allFitnessEffects(noIntGenes = c(rexp(ng/2, 10),
                                            -rexp(ng/2, 50)))

ng <- 4000
re_4000 <- allFitnessEffects(noIntGenes = c(rexp(ng/2, 10),
                                            -rexp(ng/2, 50)))

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)

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)

Can we use a large number of genes?

Yes. In fact, in OncoSimulR there is no pre-set limit on genome size. However, large numbers of genes can lead to unacceptably large returned object sizes and/or running time. We discuss several examples next that illustrate some of the major issues to consider. Another example with 50,000 genes is shown in section @ref(mcf50070).

We have seen in @ref(bench1) and @ref(common1) that for the Exp model, benchmark results using detectionProb require a lot of care and can be misleading. Here, we will fix initial population sizes (to 500) and all final population sizes will be set to  ≥ 10⁶. In addition, to avoid the confounding factor of the onlyCancer = TRUE argument, we will set it to FALSE, so we measure directly the time of individual runs.

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

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

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

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

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

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

Population sizes  ≥ 10¹⁰

We have already seen examples where population sizes reach 10⁸ to 10¹⁰, as in Tables @ref(tab:bench1b), @ref(tab:timing3), @ref(tab:timing3xf). What about even larger population sizes?

The C++ code will unconditionally alert if population sizes exceed 4 * 10¹⁵ as in those cases loosing precision (as we are using doubles) would be unavoidable, and we would also run into problems with the generation of binomial random variates (code that illustrates and discusses this problem is available in file “example-binom-problems.cpp”, in directory “/inst/miscell”). However, well before we reach 4 * 10¹⁵ we loose precision from other sources. One of the most noticeable ones is that when we reach population sizes around 10¹¹ the C++ code will often alert us by throwing exceptions with the message Recoverable exception ti set to DBL_MIN. Rerunning. I throw this exception because t_i, the random variable for time to next mutation, is less than DBL_MIN, the minimum representable floating-point number. This happens because, unless we use really tiny mutation rates, the time to a mutation starts getting closer to zero as population sizes grow very large. It might be possible to ameliorate these problems somewhat by using long doubles (instead of doubles) or special purpose libraries that provide more precision. However, this would make it harder to run the same code in different operating systems and would likely decrease execution speed on the rest of the common scenarios for which OncoSimulR has been designed.

The following code shows some examples where we use population sizes of 10¹⁰ or larger. Since we do not want simulations in the exponential model to end because of extinction, I use a fitness specification where all genes have a positive fitness effect and we start all simulations from a large population (to make it unlikely that the population will become extinct before cells mutate and start increasing in numbers). We set the maximum running time to 10 minutes. We keep the genealogy of the clones and use keepEvery = 1.

ng <- 50
u <- allFitnessEffects(noIntGenes = c(rep(0.1, ng)))

t_mc_k_50_1e11 <- system.time(
    mc_k_50_1e11 <- oncoSimulPop(5,
                                 u,
                                 model = "McFL",
                                 mu = 1e-7,
                                 detectionSize = 1e11,
                                 initSize = 1e5,
                                 detectionDrivers = NA,
                                 detectionProb = NA,
                                 keepPhylog = TRUE,
                                 onlyCancer = FALSE,
                                 mutationPropGrowth = FALSE,
                                 keepEvery = 1,
                                 finalTime = 5000,
                                 mc.cores = 1,
                                 max.wall.time = 600
                                 ))

## Recoverable exception ti set to DBL_MIN. Rerunning.
## Recoverable exception ti set to DBL_MIN. Rerunning.

t_mc_k_50_1e11
## user  system elapsed
## 613.612   0.040 613.664

summary(mc_k_50_1e11)[, c(1:3, 8, 9)]
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      5491 100328847809  44397848771  1019.950  942764
## 2      3194 100048090441  34834178374   789.675  888819
## 3      5745 100054219162  24412502660   927.950  929231
## 4      4017 101641197799  60932177160   750.725  480938
## 5      5393 100168156804  41659212367   846.250  898245

## print(object.size(mc_k_50_1e11), units = "MB")
## 177.8 Mb

We get to 10¹¹. But notice the exception with the warning about t_i. Notice also that this takes a long time and we run a very large number of iterations (getting close to one million in some cases).

Now the exponential model with detectionSize = 1e11:

t_exp_k_50_1e11 <- system.time(
    exp_k_50_1e11 <- oncoSimulPop(5,
                                  u,
                                  model = "Exp",
                                  mu = 1e-7,
                                  detectionSize = 1e11,
                                  initSize = 1e5,
                                  detectionDrivers = NA,
                                  detectionProb = NA,
                                  keepPhylog = TRUE,
                                  onlyCancer = FALSE,
                                  mutationPropGrowth = FALSE,
                                  keepEvery = 1,
                                  finalTime = 5000,
                                  mc.cores = 1,
                                  max.wall.time = 600,
                                  errorHitWallTime = FALSE,
                                  errorHitMaxTries = FALSE
                                  ))

## Recoverable exception ti set to DBL_MIN. Rerunning.
## Hitted wall time. Exiting.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Hitted wall time. Exiting.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Recoverable exception ti set to DBL_MIN. Rerunning.
## Hitted wall time. Exiting.
## Hitted wall time. Exiting.

t_exp_k_50_1e11
##     user   system  elapsed
## 2959.068    0.128 2959.556
try(summary(exp_k_50_1e11)[, c(1:3, 8, 9)])
##   NumClones TotalPopSize LargestClone FinalTime NumIter
## 1      6078  65172752616  16529682757  235.7590 1883438
## 2      5370 106476643712  24662446729  232.0000 2516675
## 3      2711  21911284363  17945303353  224.8608  543698
## 4      2838  13241462284   2944300245  216.8091  372298
## 5      7289  76166784312  10941729810  240.0217 1999489

print(object.size(exp_k_50_1e11), units = "MB")
## 53.5 Mb

Note that we almost reached max.wall.time (600 * 5 = 3000). What if we wanted to go up to 10¹²? We would not be able to do it in 10 minutes. We could set max.wall.time to a value larger than 600 to allow us to reach larger sizes but then we would be waiting for a possibly unacceptable time for simulations to finish. Moreover, this would eventually fail as simulations would keep hitting the t_i exception without ever being able to complete. Finally, even if we were very patient, hitting that t_i exception should make us worry about possible biases in the samples.

A summary of some determinants of running time and space consumption

To summarize this section, we have seen:

• Both McFL and Exp can be run in short times over a range of sizes for the detectionProb and detectionSize mechanisms using a complex fitness specification with moderate numbers of genes. These are the typical or common use cases of OncoSimulR.

• The keepEvery argument can have a large effect on time in the McFL models and specially on object sizes. If only the end result of the simulation is to be used, you should set keepEvery = NA.

• The distribution of fitness effects and the fitness landscape can have large effects on running times. Sometimes these are intuitive and simple to reason about, sometimes they are not as they interact with other factors (e.g., stopping mechanism, numbers of clones, etc). In general, there can be complex interactions between different settings, from mutation rate to fitness effects to initial size. As usual, test before launching a massive simulation.

• Simulations start to slow down and lead to a very large object size when we keep track of around 6000 to 10000 clones. Anything that leads to these patterns will slow down the simulations.

• OncoSimulR needs to keep track of genotypes (or clones), not just numbers of drivers and passengers, because it allows you to use complex fitness and mutation specifications that depend on specific genotypes. The keepEvery = NA is an approach to store only the minimal information needed, but it is unavoidable that during the simulations we might be forced to deal with many thousands of different clones.

Introduction to the specification of fitness effects

With OncoSimulR you can specify different types of effects on fitness:

• A special type of epistatic effect that is particularly amenable to be represented as a graph (a DAG). In this graph having, say, “B” be a child of “A” means that a mutation in B can only accumulate if a mutation in A is already present. This is what OT (Desper et al., 1999; Szabo & Boucher, 2008), CBN (Beerenwinkel, Eriksson, et al., 2007; Gerstung et al., 2009; Gerstung, Eriksson, et al., 2011), progression networks (Farahani & Lagergren, 2013), and other similar models (Korsunsky et al., 2014) generally mean. Details are provided in section @ref(posetslong). Note that this is not an order effect (discussed below): the fitness of a genotype from this DAGs is a function of whether or not the restrictions in the graph are satisfied, not the historical sequence of how they were satisfied.

• Effects where the order in which mutations are acquired matters, as illustrated in section @ref(oe). There is, in fact, empirical evidence of these effects (Ortmann et al., 2015). For instance, the fitness of genotype “A, B” would differ depending on whether A or B was acquired first (or, as in the actual example in (Ortmann et al., 2015), the fitness of the mutant with JAK2 and TET2 mutated will depend on which of the genes was mutated first).

• General epistatic effects (e.g., section @ref(epi)), including synthetic viability (e.g., section @ref(sv)) and synthetic lethality/mortality (e.g., section @ref(sl)).

• Genes that have independent effects on fitness (section @ref(noint)).

• Modules (see section @ref(modules0)) allow you to specify any of the above effects (except those for genes without interactions, as it would not make sense there) in terms of modules (sets of genes), not individual genes. We will introduce them right after @ref(posetslong), and we will continue using them thereafter.

A guiding design principle of OncoSimulR is to try to make the specification of those effects as simple as possible but also as flexible as possible. Thus, there are two main ways of specifying fitness effects:

• Combining different types of effects in a single specification. For instance, you can combine epistasis with order effects with no interaction genes with modules. What you would do here is specify the effects that different mutations (or their combinations) have on fitness (the fitness effects) and then have OncoSimulR take care of combining them as if each of these were lego pieces. We will refer to this as the lego system of fitness effects. (As explained above, I find this an intuitive and very graphical analogy, which I have copied from Hothorn et al. (2006) and Hothorn et al. (2008)).

• Explicitly passing to OncoSimulR a mapping of genotypes to fitness. Here you specify the fitness of each genotype. We will refer to this as the explicit mapping of genotypes to fitness. This includes frequency-dependent fitness (section @ref(fdf)).

Both approaches have advantages and disadvantages. Here I emphasize some relevant differences.

• With the lego system you can specify huge genomes with an enormous variety of interactions, since the possible genotypes are not constructed in advance. You would not be able to do this with the explicit mapping of genotypes to fitness if you wanted to, say, construct that mapping for a modest genotype of 500 genes (you’d have more genotypes than particles in the observable Universe).

• For many models/data you often intuitively start with the fitness of the genotypes, not the fitness consequences of the different mutations. In these cases, you’d need to do the math to specify the terms you want if you used the lego system so you’ll probably use the specification with the direct mapping genotype → fitness.

• Likewise, sometimes you already have a moderate size genotype → fitness mapping and you certainly do not want to do the math by hand: here the lego system would be painful to use.

• But sometimes we do think in terms of “the effects on fitness of such and such mutations are” and that immediately calls for the lego system, where you focus on the effects, and let OncoSimulR take care of doing the math of combining.

• If you want to use order effects, you must use the lego system (at least for now).

• If you want to specify modules, you must use the lego system (the explicit mapping of genotypes is, by its very nature, ill-suited for this).

• The lego system might help you see what your model really means: in many cases, you can obtain fairly succinct specifications of complex fitness models with just a few terms. Similarly, depending on what your emphasis is, you can often specify the same fitness landscape in several different ways.

Regardless of the route, you need to get that information into OncoSimulR’s functions. The main function we will use is allFitnessEffects: this is the function in charge of reading the fitness specifications. We also need to discuss how, what, and where you have to pass to allFitnessEffects.

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

fem4 <- allFitnessEffects(genotFitness = m4)
## Column names of object not Genotype and Birth Renaming them assuming that is what you wanted

try(plot(fem4))
## Error in plot.fitnessEffects(fem4) : 
##   This fitnessEffects object can not be ploted this way. It is probably one with fitness landscape specification,  so you might want to plot the fitness landscape instead.

plotFitnessLandscape(evalAllGenotypes(fem4))

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

plotFitnessLandscape(evalAllGenotypes(fem4))

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,
## frequencyDependentBirth = FALSE, : No wildtype in the fitness
## landscape!!! Adding it with birth 1.
evalAllGenotypes(fem6, addwt = TRUE)
##   Genotype Birth
## 1       WT     1
## 2        A     3
## 3        B     0
## 4     A, B     2
## plot(fem6)

plotFitnessLandscape(evalAllGenotypes(fem6))

r <- data.frame(Genotype = c("WT", "A", "B", "A, B"),
                 Fitness = c("10 * f_",
                             "10 * f_1",
                             "50 * f_2",
                             "200 * (f_1 + f_2) + 50 * f_1_2"))

afe <- allFitnessEffects(genotFitness = r,
                         frequencyDependentFitness = TRUE,
                         frequencyType = "rel")
## Warning in allFitnessEffects(genotFitness = r,
## frequencyDependentFitness = TRUE, : v2 functionality detected.
## Adapting to v3 functionality.

plotFitnessLandscape(evalAllGenotypes(afe,
                                      spPopSizes = c(WT = 2500, A = 2000,
                                                     B = 5500, "A, B" = 700)))
## Using old version of fitnessEffects. Transforming fitnessEffects
##             to last version.
## Using old version of fitnessEffects. Transforming fitnessEffects
##             to last version.

r <- data.frame(Genotype = c("WT", "A", "B", "A, B"),
                 Fitness = c("10 * f_",
                             "10 * f_A",
                             "50 * f_B",
                             "200 * (f_A + f_B) + 50 * f_A_B"))

Numeric values of fitness effects

We evaluate fitness using the usual (Beerenwinkel, Eriksson, et al., 2007; Datta et al., 2013; Gillespie, 1993; Zanini & Neher, 2012) multiplicative model: fitness is ∏(1 + s_i) where s_i is the fitness effect of gene (or gene interaction) i. In all models except Bozic, this fitness refers to the growth rate (the death rate being fixed to 1⁷). The original model of McFarland et al. (2013) has a slightly different parameterization, but you can go easily from one to the other (see section @ref(mcfl)).

For the Bozic model (Bozic et al., 2010), however, the birth rate is set to 1, and the death rate then becomes ∏(1 − s_i).

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

Genes without interactions

This is a simple scenario. Each gene i has a fitness effect s_i if mutated. The s_i can come from any distribution you want. As an example let’s use three genes. We know there are no order effects, but we will also see what happens if we examine genotypes as ordered.

ai1 <- evalAllGenotypes(allFitnessEffects(
    noIntGenes = c(0.05, -.2, .1), frequencyDependentFitness = FALSE), order = FALSE)
## Warning in allFitnessEffects(noIntGenes = c(0.05, -0.2, 0.1),
## frequencyDependentFitness = FALSE): v2 functionality detected.
## Adapting to v3 functionality.
## Using old version of fitnessEffects. Transforming fitnessEffects
##             to last version.
## Using old version of fitnessEffects. Transforming fitnessEffects
##             to last version.

We can easily verify the first results:

ai1
##   Genotype Fitness
## 1        1   1.050
## 2        2   0.800
## 3        3   1.100
## 4     1, 2   0.840
## 5     1, 3   1.155
## 6     2, 3   0.880
## 7  1, 2, 3   0.924

all(ai1[, "Fitness"]  == c( (1 + .05), (1 - .2), (1 + .1),
       (1 + .05) * (1 - .2),
       (1 + .05) * (1 + .1),
       (1 - .2) * (1 + .1),
       (1 + .05) * (1 - .2) * (1 + .1)))
## [1] TRUE

And we can see that considering the order of mutations (see section @ref(oe)) makes no difference:

(ai2 <- evalAllGenotypes(allFitnessEffects(
    noIntGenes = c(0.05, -.2, .1)), order = TRUE,
    addwt = TRUE))
##     Genotype Birth
## 1         WT 1.000
## 2          1 1.050
## 3          2 0.800
## 4          3 1.100
## 5      1 > 2 0.840
## 6      1 > 3 1.155
## 7      2 > 1 0.840
## 8      2 > 3 0.880
## 9      3 > 1 1.155
## 10     3 > 2 0.880
## 11 1 > 2 > 3 0.924
## 12 1 > 3 > 2 0.924
## 13 2 > 1 > 3 0.924
## 14 2 > 3 > 1 0.924
## 15 3 > 1 > 2 0.924
## 16 3 > 2 > 1 0.924