Browse code

2.99.92: rfitness with a three-element vector for scale and SSWM examples

ramon diaz-uriarte (at Phelsuma) authored on 27/04/2021 13:58:19
Showing8 changed files

... ...
@@ -1,8 +1,8 @@
1 1
 Package: OncoSimulR
2 2
 Type: Package
3 3
 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 
4
-Version: 2.99.9
5
-Date: 2021-04-22
4
+Version: 2.99.92
5
+Date: 2021-04-27
6 6
 Authors@R: c(
7 7
 	      person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),	
8 8
  	   		     email = "rdiaz02@gmail.com"),
... ...
@@ -202,48 +202,66 @@ rfitness <- function(g, c= 0.5,
202 202
             ## rm(m1)
203 203
             ## rm(fl1)
204 204
         }
205
-
206
-        if(!is.null(scale)) {
207
-            fi <- (fi - min(fi))/(max(fi) - min(fi))
208
-            fi <- scale[1] + fi * (scale[2] - scale[1])
209
-        }
210
-        if(wt_is_1 == "divide") {
211
-            ## We need to shift to avoid ratios involving negative numbers and
212
-            ## we need to avoid having any fitness at 0, specially the wt.  If
213
-            ## any negative value, add the min, and shift by the magnitude of
214
-            ## the min to avoid any at 0.
215
-
216
-            ## If you use scale and wt_is_1, this will move the scale. It is
217
-            ## not possible to obtain a linear transformation that will keep
218
-            ## the min and max of the scale, and wt at 1.
219
-            min_fi <- min(fi)
220
-            if(min_fi < 0)
221
-                fi <- fi + 2 * abs(min(fi))
222
-            fi <- fi/fi[1]
223
-        } else if (wt_is_1 == "subtract") {
224
-            fi <- fi - fi[1] + 1.0
225
-        } else if(wt_is_1 == "force") {
226
-            fi[1] <- 1.0
205
+        if(!(length(scale) == 3)) {
227 206
             if(!is.null(scale)) {
228
-                if( (1 > scale[2]) || (1 < scale[1]))
229
-                    warning("Using wt_is_1 = force and scale, but scale does ",
230
-                            "not include 1")
207
+                fi <- (fi - min(fi))/(max(fi) - min(fi))
208
+                fi <- scale[1] + fi * (scale[2] - scale[1])
209
+            }
210
+            if(wt_is_1 == "divide") {
211
+                ## We need to shift to avoid ratios involving negative numbers and
212
+                ## we need to avoid having any fitness at 0, specially the wt.  If
213
+                ## any negative value, add the min, and shift by the magnitude of
214
+                ## the min to avoid any at 0.
215
+
216
+                ## If you use scale and wt_is_1, this will move the scale. It is
217
+                ## not possible to obtain a linear transformation that will keep
218
+                ## the min and max of the scale, and wt at 1.
219
+                min_fi <- min(fi)
220
+                if(min_fi < 0)
221
+                    fi <- fi + 2 * abs(min(fi))
222
+                fi <- fi/fi[1]
223
+            } else if (wt_is_1 == "subtract") {
224
+                fi <- fi - fi[1] + 1.0
225
+            } else if(wt_is_1 == "force") {
226
+                fi[1] <- 1.0
227
+                if(!is.null(scale)) {
228
+                    if( (1 > scale[2]) || (1 < scale[1]))
229
+                        warning("Using wt_is_1 = force and scale, but scale does ",
230
+                                "not include 1")
231
+                }
231 232
             }
233
+        } else { ## a length-3 scale
234
+            if(scale[2] > scale[3]) warning("In scale, minimum fitness > wildtype")
235
+            fiwt <- fi[1]
236
+            prod_above <- (scale[1] - scale[3]) / (max(fi) - fiwt)
237
+            prod_below <- (scale[3] - scale[2]) / (fiwt - min(fi))
238
+            fi_above <- which(fi >= fiwt)
239
+            fi_below <- which(fi < fiwt)
240
+            fi[fi_above] <- ((fi[fi_above] - fiwt) * prod_above) + scale[3]
241
+            fi[fi_below] <- ((fi[fi_below] - fiwt) * prod_below) + scale[3]
242
+            fi[1] <- scale[3]
232 243
         }
244
+        
233 245
         if(log) {
246
+            if(length(scale == 3)) {
247
+                wt_is_1 <- "no"
248
+                message("You passed a three-element scale argument. ",
249
+                        "Setting wt_is_1 to no ",
250
+                        "to avoid modifying your requested value for WT.")
251
+            }
234 252
             ## If you want logs, you certainly do not want
235 253
             ## the log of a negative number
236 254
             fi[fi < 0] <- 0
237 255
             if(wt_is_1 == "no") {
238 256
                 fi <- log(fi)
239
-            } else {
240
-                ## by decree, fitness of wt is 1. So shift everything
241
-                fi <- log(fi) + 1
242
-            }
243
-            ## former expression, but it was more confusing
257
+                } else {
258
+                    ## by decree, fitness of wt is 1. So shift everything
259
+                    fi <- log(fi) + 1
260
+                }
261
+                ## former expression, but it was more confusing
244 262
             ## fi <- log(fi/fi[1]) + 1
245
-            
246 263
         }
264
+        
247 265
         if(truncate_at_0) {
248 266
             ## yes, truncate but add noise to prevent identical
249 267
             fi[fi <= 0] <- runif(sum(fi <= 0),
... ...
@@ -1,3 +1,7 @@
1
+Changes in version 2.99.92 (2021-04-27):
2
+	- rfitness: scale can take a three-element vector.
3
+	- Vignette: examples (not run) for deviations from SSWM.
4
+
1 5
 Changes in version 2.99.9 (2021-04-22):
2 6
 	- Fixed date typo in one citation.
3 7
 	- We were inconsistent, allowing some examples of one gene.
... ...
@@ -55,9 +55,36 @@ additive models.}
55 55
   genotypes with that number of mutations have equal probability of
56 56
   being the reference). }
57 57
 
58
-\item{scale}{Either NULL (nothing is done) or a two-element vector. If a
59
-  two-element vector, fitness is re-scaled between \code{scale[1]} (the
60
-  minimum) and \code{scale[2]} (the maximum).}
58
+\item{scale}{Either NULL (nothing is done) or a two- or three-element
59
+  vector.
60
+
61
+  If a two-element vector, fitness is re-scaled between
62
+  \code{scale[1]} (the minimum) and \code{scale[2]} (the maximum) and,
63
+  later, if you have selected it, \code{wt_is_1} will be enforced.
64
+
65
+  If you pass a three element vector, fitness is re-scaled so that the
66
+  new maximum fitness is \code{scale[1]}, the new minimum is
67
+  \code{scale[2]} and the new wildtype is \code{scale[3]}. If you pass a
68
+  three element vector, none of the \code{wt_is_1} options apply in this
69
+  case, to ensure you obtain the range you want. If you want the
70
+  wildtype to be one, pass it as the third element of the vector.
71
+
72
+  As a consequence of using a three element vector, the amount of
73
+  stretching/compressing (i.e., scaling) of fitness values larger than
74
+  that of the wildtype will likely be different from the scaling of
75
+  fitness values smaller than that of the wildtype.  In other words,
76
+  this argument allows you to change the spread of the positive and
77
+  negative fitness values (and you can make this difference extreme and
78
+  make most fitness values less than wildtype be 0 by using a huge
79
+  negative number --huge in absolute value-- for \code{scale[2]} if you
80
+  then truncate at 0 --see \code{truncate_at_9}).
81
+
82
+  Using a three element vector is probably the most natural way of
83
+  changing the scale and range of fitness.
84
+
85
+  See also \code{log} if you want the log-transformed values to respect
86
+  the scale.
87
+}
61 88
 
62 89
 \item{wt_is_1}{If "divide" the fitness of all genotypes is
63 90
   divided by the fitness of the wildtype (after possibly adding a value
... ...
@@ -83,14 +110,28 @@ additive models.}
83 110
   option can easily lead to landscapes with no accessible genotypes
84 111
   (even if you also use \code{scale}).
85 112
 
86
-  If "no", the fitness of the wildtype is not modified.  }
113
+  If "no", the fitness of the wildtype is not modified.
114
+
115
+  This option has no effect if you pass a three-element vector for
116
+  \code{scale}. Using a three-element vector for \code{scale} is
117
+  probably the most natural way of changing the scale and range of
118
+  fitness while setting the wildtype to value of your choice.
119
+  
120
+}
87 121
 
88 122
 
89 123
 \item{log}{If TRUE, log-transform fitness. Actually, there are two
90 124
   cases: if \code{wt_is_1 = "no"} we simply log the fitness values;
91 125
   otherwise, we log the fitness values and add a 1, thus shifting all
92 126
   fitness values, because by decree the fitness (birth rate) of the
93
-  wildtype must be 1.}
127
+  wildtype must be 1.
128
+
129
+  If you pass a three-element vector for scale, you will want to pass
130
+  \code{exp(desired_max)}, \code{exp(desired_min)}, and
131
+  \code{exp(desired_wildtype)} to the \code{scale} argument. (We first
132
+  scale values in the original scale and then log them). In this case,
133
+  we ignore whatever you passed as \code{wt_is_1}, setting \code{wt_is_1
134
+  = "no"} to avoid modifying your requested value for the wildtype.}
94 135
 
95 136
 \item{min_accessible_genotypes}{If not NULL, the minimum number of
96 137
   accessible genotypes in the fitness landscape. A genotype is
... ...
@@ -1,3 +1,4 @@
1
+## This does test some options, but mainly exercises others
1 2
 inittime <- Sys.time()
2 3
 cat(paste("\n Starting exercise-rfitness at", date(), "\n"))
3 4
 test_that("Expect output", {
... ...
@@ -60,45 +61,137 @@ test_that("Error if wrong arguments", {
60 61
 test_that("Additive, Ising, Eggbox, Full exercising", {
61 62
 
62 63
     expect_output(print(rfitness(4, model = "Ising")), "Fitness", fixed = TRUE)
63
-expect_output(print(rfitness(4, model = "Eggbox")), "Fitness", fixed = TRUE)
64
-expect_output(print(rfitness(4, model = "Full")), "Fitness", fixed = TRUE)
65
-expect_output(print(rfitness(4, model = "Additive")), "Fitness", fixed = TRUE)
66
-expect_output(print(rfitness(4, model = "Ising", i = 0.0002, I = 0.5,
67
-                             circular = TRUE)), "Fitness", fixed = TRUE)
68
-expect_output(print(rfitness(4, model = "Ising", i = 2, I = 0,
69
-                             circular = TRUE)), "Fitness", fixed = TRUE)
70
-expect_output(print(rfitness(4, model = "Ising", i = 2, I = -2,
71
-                             circular = TRUE)), "Fitness", fixed = TRUE)
72
-expect_output(print(rfitness(4, model = "Ising", i = 2, I = 0.5,
73
-                             circular = FALSE)), "Fitness", fixed = TRUE)
74
-expect_output(print(rfitness(4, model = "Eggbox", e = 2, E = 0)),
75
-              "Fitness", fixed = TRUE)
76
-expect_output(print(rfitness(4, model = "Eggbox", e = 2, E = 2.4)),
77
-              "Fitness", fixed = TRUE)
78
-expect_output(print(rfitness(4, model = "Eggbox", e = 0, E = 2.4)),
79
-              "Fitness", fixed = TRUE)
80
-expect_output(print(rfitness(4, model = "Full", i = 0.002, I = 2,
81
-                K = 2, r = TRUE,
82
-                p = 0.2, P = 0.3, o = 0.3, O = 1)), "Fitness", fixed = TRUE)
83
-expect_error(rfitness(4, model = "Full", K = 5), "It makes no sense",
84
-              fixed = TRUE)
85
-expect_output(print(rfitness(4, model = "Full", i = 0.002, I = 2,
86
-                K = 2, r = TRUE,
87
-                p = 0.2, P = 0.3, o = 0.3, O = 1,
88
-                s = 0.5, S=0, d = 0.1,
89
-                e = 1, E = 1.5,
90
-                H = 0.5)), "Fitness", fixed = TRUE)
91
-expect_output(print(rfitness(4, model = "Additive", mu = 0, sd = 1)),
92
-              "Fitness", fixed = TRUE)
93
-expect_output(print(rfitness(4, model = "Additive", mu = 1, sd = 0)),
94
-              "Fitness", fixed = TRUE)
95
-expect_output(print(rfitness(4, model = "Additive", mu = 3, sd = 1)),
96
-              "Fitness", fixed = TRUE)
97
-expect_output(print(rfitness(4, model = "Additive", mu = -4, sd = 10)),
98
-              "Fitness", fixed = TRUE)
64
+    expect_output(print(rfitness(4, model = "Eggbox")), "Fitness", fixed = TRUE)
65
+    expect_output(print(rfitness(4, model = "Full")), "Fitness", fixed = TRUE)
66
+    expect_output(print(rfitness(4, model = "Additive")), "Fitness", fixed = TRUE)
67
+    expect_output(print(rfitness(4, model = "Ising", i = 0.0002, I = 0.5,
68
+                                 circular = TRUE)), "Fitness", fixed = TRUE)
69
+    expect_output(print(rfitness(4, model = "Ising", i = 2, I = 0,
70
+                                 circular = TRUE)), "Fitness", fixed = TRUE)
71
+    expect_output(print(rfitness(4, model = "Ising", i = 2, I = -2,
72
+                                 circular = TRUE)), "Fitness", fixed = TRUE)
73
+    expect_output(print(rfitness(4, model = "Ising", i = 2, I = 0.5,
74
+                                 circular = FALSE)), "Fitness", fixed = TRUE)
75
+    expect_output(print(rfitness(4, model = "Eggbox", e = 2, E = 0)),
76
+                  "Fitness", fixed = TRUE)
77
+    expect_output(print(rfitness(4, model = "Eggbox", e = 2, E = 2.4)),
78
+                  "Fitness", fixed = TRUE)
79
+    expect_output(print(rfitness(4, model = "Eggbox", e = 0, E = 2.4)),
80
+                  "Fitness", fixed = TRUE)
81
+    expect_output(print(rfitness(4, model = "Full", i = 0.002, I = 2,
82
+                                 K = 2, r = TRUE,
83
+                                 p = 0.2, P = 0.3, o = 0.3, O = 1)), "Fitness", fixed = TRUE)
84
+    expect_error(rfitness(4, model = "Full", K = 5), "It makes no sense",
85
+                 fixed = TRUE)
86
+    expect_output(print(rfitness(4, model = "Full", i = 0.002, I = 2,
87
+                                 K = 2, r = TRUE,
88
+                                 p = 0.2, P = 0.3, o = 0.3, O = 1,
89
+                                 s = 0.5, S=0, d = 0.1,
90
+                                 e = 1, E = 1.5,
91
+                                 H = 0.5)), "Fitness", fixed = TRUE)
92
+    expect_output(print(rfitness(4, model = "Additive", mu = 0, sd = 1)),
93
+                  "Fitness", fixed = TRUE)
94
+    expect_output(print(rfitness(4, model = "Additive", mu = 1, sd = 0)),
95
+                  "Fitness", fixed = TRUE)
96
+    expect_output(print(rfitness(4, model = "Additive", mu = 3, sd = 1)),
97
+                  "Fitness", fixed = TRUE)
98
+    expect_output(print(rfitness(4, model = "Additive", mu = -4, sd = 10)),
99
+                  "Fitness", fixed = TRUE)
99 100
 
100 101
 })
102
+ 
103
+test_that("Testing the three-element scale argument", {
104
+    for(j in 1:5) {
105
+        ## First, WT not equal to 1
106
+        ## Don't use too stupid arguments
107
+        scalev <- sort(runif(3, -2, 5), decreasing = TRUE)
108
+        scalev <- scalev[c(1, 3, 2)]
109
+        i <- round(100 * runif(1))
110
+        set.seed(i)
111
+        ra <- rfitness(7, truncate_at_0 = FALSE, wt_is_1 = "no")
112
+        set.seed(i)
113
+        rb <- rfitness(7, scale = scalev, truncate_at_0 = FALSE)
114
+
115
+        rap <- ra[, "Fitness"][ra[, "Fitness"] >= ra[1, "Fitness"]]
116
+        ran <- ra[, "Fitness"][ra[, "Fitness"] < ra[1, "Fitness"]]
117
+        rbp <- rb[, "Fitness"][rb[, "Fitness"] >= scalev[3]]
118
+        rbn <- rb[, "Fitness"][rb[, "Fitness"] < scalev[3]]
119
+
120
+        expect_equivalent(cor(rap, rbp), 1)
121
+        expect_equivalent(cor(ran, rbn), 1)
122
+
123
+        expect_equivalent(max(rb[, "Fitness"]), scalev[1])
124
+        expect_equivalent(min(rb[, "Fitness"]), scalev[2])
125
+        expect_equivalent(rb[1, "Fitness"], scalev[3])
126
+
127
+        ## X, Y, Z: X: original, Y: after scaling
128
+        ## Just dealing with the positive
129
+        ## Do not assume WT to be 1
130
+
131
+        ## Write as original * slope + intercept
132
+        ## so you can compare against
133
+        ## lm(rbp ~ rap)
134
+        xM <- max(ra[, "Fitness"])
135
+        xW <- ra[1, "Fitness"]
136
+        y <- rap * (scalev[1] - scalev[3])/(xM - xW) +
137
+            (scalev[3] - xW * ((scalev[1] - scalev[3])/(xM - xW)))
138
+
139
+        expect_equivalent(rbp, y)
140
+        
141
+        set.seed(i)
142
+        rc <- rfitness(7, scale = exp(scalev), truncate_at_0 = FALSE, log = TRUE)
143
+        
144
+        rcp <- rc[, "Fitness"][rc[, "Fitness"] >= scalev[3]]
145
+        rcn <- rc[, "Fitness"][rc[, "Fitness"] < scalev[3]]
146
+
147
+        expect_equivalent(cor(rap, exp(rcp)), 1)
148
+        expect_equivalent(cor(ran, exp(rcn)), 1)
149
+
150
+        ## X, U, Z: X: original, U: after scaling; Z: after log
151
+        ## Just dealing with the positive
152
+        ## Do not assume WT to be 1
153
+        ## can compare against lm(exp(rcp) ~ rap)
101 154
 
155
+        u <- rap * (exp(scalev[1]) - exp(scalev[3]))/(xM - xW) +
156
+            (exp(scalev[3]) - xW * ((exp(scalev[1]) - exp(scalev[3]))/(xM - xW)))
157
+        expect_equivalent(u, exp(rcp))
158
+
159
+        ## This shows, by the way, that the log-transformed ain't linear function of
160
+        ## log of original
161
+        z <- log( (rap - xW) * ((exp(scalev[1]) - exp(scalev[3]))/(xM - xW)) + exp(scalev[3]))
162
+        expect_equivalent(z, rcp)
163
+
164
+        ############################################
165
+        ## Repeat setting WT = 1
166
+        set.seed(i)
167
+        rax <- rfitness(7, truncate_at_0 = FALSE, wt_is_1 = "subtract")
168
+        rapx <- rax[, "Fitness"][rax[, "Fitness"] >= 1]
169
+        ranx <- rax[, "Fitness"][rax[, "Fitness"] < 1]
170
+
171
+        expect_equivalent(cor(rapx, rbp), 1)
172
+        expect_equivalent(cor(ranx, rbn), 1)
173
+        
174
+        xMx <- max(rax[, "Fitness"])
175
+        xWx <- rax[1, "Fitness"]
176
+        yx <- rapx * (scalev[1] - scalev[3])/(xMx - xWx) +
177
+            (scalev[3] - xWx * ((scalev[1] - scalev[3])/(xMx - xWx)))
178
+
179
+        expect_equivalent(rbp, yx)
180
+
181
+        ## log
182
+        expect_equivalent(cor(rapx, exp(rcp)), 1)
183
+        expect_equivalent(cor(ranx, exp(rcn)), 1)
184
+
185
+        ux <- rapx * (exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx) +
186
+            (exp(scalev[3]) - xWx * ((exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx)))
187
+        expect_equivalent(ux, exp(rcp))
188
+
189
+        ## This shows, by the way, that the log-transformed ain't linear function of
190
+        ## log of original
191
+        zx <- log( (rapx - xWx) * ((exp(scalev[1]) - exp(scalev[3]))/(xMx - xWx)) + exp(scalev[3]))
192
+        expect_equivalent(zx, rcp)
193
+    }
194
+})
102 195
 
103 196
 cat(paste("\n Ending exercise-rfitness at", date(), "\n"))
104 197
 cat(paste("  Took ", round(difftime(Sys.time(), inittime, units = "secs"), 2), "\n\n"))
... ...
@@ -515,7 +515,8 @@ encompass a wide range of issues that have been addressed in evolutionary
515 515
 genetics studies and which include from detailed analysis of simple models
516 516
 with a few uphill paths and valleys as in @Weissman2009 or @Ochs2015, to
517 517
 questions that refer to larger, more complex fitness landscapes as in
518
-@szendro_predictability_2013 or @franke_evolutionary_2011 (see below).
518
+@szendro_predictability_2013 or @franke_evolutionary_2011 or
519
+@krug2019 (see below).
519 520
 
520 521
 
521 522
 Using as an example @Ochs2015 (we will see this example again in section
... ...
@@ -604,7 +605,8 @@ ruggedness, and would then examine the evolutionary predictability
604 605
 of the trajectories with measures such as "Lines of Descent" and
605 606
 "Path of the Maximum" [@szendro_predictability_2013] and the
606 607
 diversity of the sampled genotypes under different sampling regimes (see
607
-details in section \@ref(evolpredszend)).
608
+details in section \@ref(evolpredszend)). (See also related comments
609
+in section \@ref(sswm-rfitness)).
608 610
 
609 611
 
610 612
 
... ...
@@ -999,6 +1001,137 @@ course, for some types of problems this special purpose software might not
999 1001
 be available, though.
1000 1002
 
1001 1003
 
1004
+## Random fitness landscapes, clonal competition, predictability, and the strong selection weak mutation (SSWM) regime {#sswm-rfitness}
1005
+
1006
+Many studies about evolutionary predictability (among other topics)
1007
+focus on the strong selection, weak mutation regime, SSWM
1008
+[@gillespie1984; @orr_population_2002] (see overview in
1009
+@krug2019). In this regime, mutations are rare (much smaller than
1010
+the mutation rate times the population size) and selection is strong
1011
+(much larger than 1/population size), so that the population
1012
+consists of a single clone most of the time, and evolution proceeds
1013
+by complete, successive clonal expansions of advantageous mutations.
1014
+
1015
+We can easily simulate variations around these scenarios with
1016
+OncoSimulR, moving away from the SSWM by increasing the population
1017
+size, or changing the size of the fitness differences. 
1018
+
1019
+The examples below, not run for the sake of speed, play with
1020
+population size and fitness differences. To make sure we use a
1021
+similar fitness landscape, we use the same simulated fitness
1022
+landscape, scaled differently, so that the differences in fitness
1023
+between mutants are increased or decreased while keeping their
1024
+ranking identical (and, thus, having the same set of accessible and
1025
+inaccessible genotypes and paths over the landscape).
1026
+
1027
+If you run the code, you will see that as we increase population
1028
+size we move further away from the SSWM: the population is no longer
1029
+composed of a single clone most of the time.
1030
+
1031
+
1032
+Before running the examples, and to show the effects quantitatively,
1033
+we define a simple wrapper to compute a few statistics.
1034
+
1035
+```{r sswm-rfitness-funct}
1036
+## oncoSimul object  -> measures of clonal interference
1037
+##    they are not averaged over time. One value for sampled time
1038
+clonal_interf_per_time <- function(x) {
1039
+    x <- x$pops.by.time
1040
+    y <- x[, -1, drop = FALSE]
1041
+    shannon <- apply(y, 1, OncoSimulR:::shannonI)
1042
+    tot <- rowSums(y)
1043
+    half_tot <- tot * 0.5
1044
+    five_p_tot <- tot * 0.05
1045
+    freq_most_freq <- apply(y/tot, 1, max)
1046
+    single_more_half <- rowSums(y > half_tot)
1047
+    ## whether more than 1 clone with more than 5% pop.
1048
+    how_many_gt_5p <- rowSums(y > five_p_tot)
1049
+    several_gt_5p <- (how_many_gt_5p > 1)
1050
+    return(cbind(shannon, ## Diversity of clones
1051
+                 freq_most_freq, ## Frequency of the most freq. clone
1052
+                 single_more_half, ## Any clone with a frequency > 50%?
1053
+                 several_gt_5p, ## Are there more than 1 clones with
1054
+                                ## frequency > 5%?
1055
+                 how_many_gt_5p ## How many clones are there with
1056
+                                ## frequency > 5%
1057
+                 ))
1058
+}
1059
+```
1060
+
1061
+
1062
+```{r sswm-rfitness0, eval=FALSE}
1063
+set.seed(1)
1064
+r7b <- rfitness(7, scale = c(1.2, 0, 1))
1065
+
1066
+## Large pop sizes: clonal interference
1067
+(sr7b <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b),
1068
+                       model = "McFL",
1069
+                       mu = 1e-6,
1070
+                       onlyCancer = FALSE,
1071
+                       finalTime = 400,
1072
+                       initSize = 1e7,
1073
+                       keepEvery = 4,
1074
+                       detectionSize = 1e10))
1075
+
1076
+plot(sr7b, show = "genotypes")
1077
+
1078
+colMeans(clonal_interf_per_time(sr7b))
1079
+```
1080
+
1081
+```{r sswm-rfitness, eval=FALSE}
1082
+## Small pop sizes: a single clone most of the time
1083
+(sr7c <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b),
1084
+                       model = "McFL",
1085
+                       mu = 1e-6,
1086
+                       onlyCancer = FALSE,
1087
+                       finalTime = 60000,
1088
+                       initSize = 1e3,
1089
+                       keepEvery = 4,
1090
+                       detectionSize = 1e10))
1091
+
1092
+plot(sr7c, show = "genotypes")
1093
+
1094
+colMeans(clonal_interf_per_time(sr7c))
1095
+
1096
+
1097
+
1098
+## Even smaller fitness differences, but large pop. sizes
1099
+set.seed(1); r7b2 <- rfitness(7, scale = c(1.05, 0, 1))
1100
+
1101
+(sr7b2 <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b2),
1102
+                       model = "McFL",
1103
+                       mu = 1e-6,
1104
+                       onlyCancer = FALSE,
1105
+                       finalTime = 3500,
1106
+                       initSize = 1e7,
1107
+                       keepEvery = 4,
1108
+                       detectionSize = 1e10))
1109
+sr7b2
1110
+plot(sr7b2, show = "genotypes")
1111
+colMeans(clonal_interf_per_time(sr7b2))
1112
+
1113
+
1114
+## Increase pop size further
1115
+(sr7b3 <- oncoSimulIndiv(allFitnessEffects(genotFitness = r7b2),
1116
+                       model = "McFL",
1117
+                       mu = 1e-6,
1118
+                       onlyCancer = FALSE,
1119
+                       finalTime = 1500,
1120
+                       initSize = 1e8,
1121
+                       keepEvery = 4,
1122
+                       detectionSize = 1e10))
1123
+sr7b3
1124
+plot(sr7b3, show = "genotypes")
1125
+colMeans(clonal_interf_per_time(sr7b3))
1126
+
1127
+```
1128
+
1129
+
1130
+
1131
+
1132
+\clearpage
1133
+
1134
+
1002 1135
 
1003 1136
 
1004 1137
 ## Steps for using OncoSimulR {#steps}
... ...
@@ -8357,6 +8490,9 @@ Magellan_stats(rnk2)
8357 8490
 fitness landscapes; with frequency-dependent fitness, as in section
8358 8491
 \@ref{fdf} fitness landscapes as such are not defined.)
8359 8492
 
8493
+
8494
+
8495
+
8360 8496
 \clearpage
8361 8497
 
8362 8498
 
... ...
@@ -13637,8 +13773,12 @@ paste("Sum times of chunks = ", sum(unlist(all_times))/60, " minutes")
13637 13773
 
13638 13774
 # Funding
13639 13775
 
13640
-Supported by BFU2015-67302-R (MINECO/FEDER, EU) to
13641
-RDU. S. Sanchez-Carrillo partially supported by a "Beca de Colaboración" from Universidad Autónoma de Madrid (2017-2018); J. Antonio Miguel supported by PEJ-2019-AI/BMD-13961 from Comunidad de Madrid.
13776
+Supported by BFU2015-67302-R (MINECO/FEDER, EU) and
13777
+PID2019-111256RB-I00/AEI/10.13039/501100011033 to
13778
+RDU. S. Sanchez-Carrillo partially supported by a "Beca de
13779
+Colaboración" from Universidad Autónoma de Madrid (2017-2018);
13780
+J. Antonio Miguel supported by PEJ-2019-AI/BMD-13961 from Comunidad
13781
+de Madrid.
13642 13782
 
13643 13783
 
13644 13784
 ```{r, echo=FALSE, eval=TRUE}
... ...
@@ -929,3 +929,50 @@ volume = {274},
929 929
   year = {2013},
930 930
   pages = {1056–1062}
931 931
 }
932
+
933
+
934
+@article{krug2019,
935
+  title = {Accessibility Percolation in Random Fitness Landscapes},
936
+  author = {Krug, Joachim},
937
+  year = {2019},
938
+  month = mar,
939
+  url = {http://arxiv.org/abs/1903.11913},
940
+  urldate = {2021-04-23},
941
+  archiveprefix = {arXiv},
942
+  eprint = {1903.11913},
943
+  eprinttype = {arxiv},
944
+  journal = {arXiv:1903.11913 [math, q-bio]},
945
+}
946
+
947
+
948
+
949
+@article{orr_population_2002,
950
+  title = {The Population Genetics of Adaptation: The Adaptation of Dna Sequences},
951
+  shorttitle = {The Population Genetics of Adaptation},
952
+  author = {Orr, H. Allen},
953
+  year = {2002},
954
+  month = jul,
955
+  volume = {56},
956
+  pages = {1317--1330},
957
+  issn = {0014-3820},
958
+  doi = {10.1554/0014-3820(2002)056[1317:TPGOAT]2.0.CO;2},
959
+  url = {http://www.bioone.org/doi/abs/10.1554/0014-3820%282002%29056%5B1317%3ATPGOAT%5D2.0.CO%3B2},
960
+  journal = {Evolution},
961
+  number = {7}
962
+}
963
+
964
+
965
+
966
+@article{gillespie1984,
967
+  title = {Molecular {{Evolution Over}} the {{Mutational Landscape}}},
968
+  author = {Gillespie, John H.},
969
+  year = {1984},
970
+  volume = {38},
971
+  pages = {1116--1129},
972
+  publisher = {{[Society for the Study of Evolution, Wiley]}},
973
+  issn = {0014-3820},
974
+  doi = {10.2307/2408444},
975
+  url = {https://www.jstor.org/stable/2408444},
976
+  journal = {Evolution},
977
+  number = {5}
978
+}
... ...
@@ -1,15 +1,15 @@
1 1
 \usepackage[%
2
-		shash={e936e3f},
3
-		lhash={e936e3f8cc52b926f7b16691b9fd58aae19a465e},
2
+		shash={bb4ad3b},
3
+		lhash={bb4ad3b194b97eba03a197a4e6d432f1a1e9167f},
4 4
 		authname={ramon diaz-uriarte (at Phelsuma)},
5 5
 		authemail={rdiaz02@gmail.com},
6
-		authsdate={2021-03-20},
7
-		authidate={2021-03-20 15:18:00 +0100},
8
-		authudate={1616249880},
6
+		authsdate={2021-04-27},
7
+		authidate={2021-04-27 15:43:25 +0200},
8
+		authudate={1619531005},
9 9
 		commname={ramon diaz-uriarte (at Phelsuma)},
10 10
 		commemail={rdiaz02@gmail.com},
11
-		commsdate={2021-03-20},
12
-		commidate={2021-03-20 15:18:00 +0100},
13
-		commudate={1616249880},
14
-		refnames={ (HEAD -> freq-dep-fitness, origin/freq-dep-fitness)}
11
+		commsdate={2021-04-27},
12
+		commidate={2021-04-27 15:43:25 +0200},
13
+		commudate={1619531005},
14
+		refnames={ (HEAD -> freq-dep-fitness, origin/rfitness_scale_3, rfitness_scale_3)}
15 15
 	]{gitsetinfo}
16 16
\ No newline at end of file