Browse code

1.99.7.: initMutant available in oncoSimulSample/Pop

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/OncoSimulR@108872 bc3139a8-67e5-0310-9ffc-ced21a209358

Ramon Diaz-Uriarte authored on 27/09/2015 12:39:23
Showing 8 changed files

... ...
@@ -1,8 +1,8 @@
1 1
 Package: OncoSimulR
2 2
 Type: Package
3 3
 Title: Forward Genetic Simulation of Cancer Progresion with Epistasis 
4
-Version: 1.99.6
5
-Date: 2015-09-26
4
+Version: 1.99.7
5
+Date: 2015-09-27
6 6
 Author: Ramon Diaz-Uriarte. 
7 7
 Maintainer: Ramon Diaz-Uriarte <rdiaz02@gmail.com>
8 8
 Description: Functions for forward population genetic simulation in
... ...
@@ -56,9 +56,10 @@ oncoSimulSample <- function(Nindiv,
56 56
                             ## well, obviously they are errors
57 57
                             ## errorHitWallTime = TRUE,
58 58
                             ## errorHitMaxTries = TRUE,
59
-                            verbosity  = 1,
60 59
                             typeSample = "whole",
61 60
                             thresholdWhole = 0.5,
61
+                            initMutant = NULL,
62
+                            verbosity  = 1,
62 63
                             seed = "auto"){
63 64
     ## No longer using mclapply, because of the way we use the limit on
64 65
     ## the number of tries.
... ...
@@ -174,6 +175,7 @@ oncoSimulSample <- function(Nindiv,
174 175
                                errorHitWallTime = TRUE,
175 176
                                errorHitMaxTries = TRUE,
176 177
                                seed = seed,
178
+                               initMutant = initMutant,
177 179
                                keepPhylog = keepPhylog)
178 180
         
179 181
         if(tmp$other$UnrecoverExcept) {
... ...
@@ -211,8 +213,10 @@ oncoSimulSample <- function(Nindiv,
211 213
             }
212 214
             return(list(
213 215
                 popSummary = summary(pop),
214
-                popSample = samplePop(pop, typeSample = typeSample,
215
-                    thresholdWhole = thresholdWhole, geneNames = geneNames),
216
+                popSample = samplePop(pop, timeSample = "last",
217
+                                      typeSample = typeSample,
218
+                                      thresholdWhole = thresholdWhole,
219
+                                      geneNames = geneNames),
216 220
                 attemptsUsed = attemptsUsed,
217 221
                 probCancer = Nindiv/attemptsUsed,
218 222
                 HittedMaxTries = FALSE,
... ...
@@ -294,6 +298,7 @@ oncoSimulPop <- function(Nindiv,
294 298
                          max.num.tries = 500,
295 299
                          errorHitWallTime = TRUE,
296 300
                          errorHitMaxTries = TRUE,
301
+                         initMutant = NULL,
297 302
                          verbosity  = 0,
298 303
                          mc.cores = detectCores(),
299 304
                          seed = "auto") {
... ...
@@ -328,7 +333,8 @@ oncoSimulPop <- function(Nindiv,
328 333
                         errorHitWallTime = errorHitWallTime,
329 334
                         errorHitMaxTries = errorHitMaxTries,
330 335
                         verbosity = verbosity,
331
-                        seed = seed, keepPhylog = keepPhylog
336
+                        seed = seed, keepPhylog = keepPhylog,
337
+                        initMutant = initMutant
332 338
                     ),
333 339
                     mc.cores = mc.cores
334 340
                     )
... ...
@@ -895,7 +901,7 @@ get.mut.vector <- function(x, timeSample, typeSample,
895 901
     if(typeSample %in% c("wholeTumor", "whole")) {
896 902
         popSize <- x$PerSampleStats[the.time, 1]
897 903
         return( as.numeric((tcrossprod(pop,
898
-                                       x$Genotypes)/popSize) > thresholdWhole) )
904
+                                       x$Genotypes)/popSize) >= thresholdWhole) )
899 905
     } else if (typeSample %in%  c("singleCell", "single")) {
900 906
 
901 907
           return(x$Genotypes[, sample(seq_along(pop), 1, prob = pop)])
... ...
@@ -1,3 +1,6 @@
1
+Changes in version 1.99.7 (2015-09-27):
2
+	- initMutant available in oncoSimulPop and oncoSimulSample
3
+
1 4
 Changes in version 1.99.6 (2015-09-26):
2 5
 	- Improved test coverage and removed stringsAsfactors from tests
3 6
 	- Consistent handling of corner cases in Bozic
... ...
@@ -51,7 +51,9 @@ oncoSimulPop(Nindiv, fp, model = "Exp", numPassengers = 30, mu = 1e-6,
51 51
                 max.num.tries = 500,
52 52
                 errorHitWallTime = TRUE,
53 53
                 errorHitMaxTries = TRUE,
54
-                verbosity = 0, mc.cores = detectCores(),
54
+                initMutant = NULL,
55
+                verbosity = 0,
56
+                mc.cores = detectCores(),
55 57
                 seed = "auto")
56 58
 
57 59
 
... ...
@@ -90,9 +92,10 @@ oncoSimulSample(Nindiv,
90 92
                 max.memory = 2000,
91 93
                 max.wall.time.total = 600,
92 94
                 max.num.tries.total = 500 * Nindiv,
93
-                verbosity  = 1,
94 95
                 typeSample = "whole",
95 96
                 thresholdWhole = 0.5,
97
+                initMutant = NULL,
98
+                verbosity  = 1,
96 99
                 seed = "auto")
97 100
 
98 101
 
... ...
@@ -102,6 +105,7 @@ oncoSimulSample(Nindiv,
102 105
   
103 106
   \item{Nindiv}{Number of individuals or number of different
104 107
     trajectories to simulate.}
108
+
105 109
   \item{fp}{
106 110
     
107 111
     Either a poset that specifies the order restrictions (see
... ...
@@ -372,15 +376,14 @@ to be detected.  For \code{oncoSimulSample} this can be a vector.
372 376
   you set it to "auto", then if you are using v.1, the behavior is the
373 377
   same as if you set it to NULL (a seed will be generated in R and
374 378
   passed to C++) but if you are using v.2, a random seed will be
375
-  produced in C++. %% using the randutils code from M.E.O'Neill.
376
-  If you need
377
-  reproducibility, either pass a value or set it to NULL (setting it to
378
-  NULL will make the C++ seed reproducible if you use the same seed in R
379
-  via \code{set.seed}). However, even using the same value of
379
+  produced in C++. %% using the randutils code from M.E.O'Neill.  If you
380
+  need reproducibility, either pass a value or set it to NULL (setting
381
+  it to NULL will make the C++ seed reproducible if you use the same
382
+  seed in R via \code{set.seed}). However, even using the same value of
380 383
   \code{seed} is unlikely to give the exact same results between
381
-  platforms and compilers.  Note that the defaults are not
382
-  the same in \code{oncoSimulIndiv}, \code{oncoSimulPop} and
383
-  \code{oncoSimulSample}. }
384
+  platforms and compilers.  Moreover, note that the defaults for
385
+  \code{seed} are not the same in \code{oncoSimulIndiv},
386
+  \code{oncoSimulPop} and \code{oncoSimulSample}. }
384 387
 
385 388
 
386 389
 
... ...
@@ -32,7 +32,7 @@ samplePop(x, timeSample = "last", typeSample = "whole",
32 32
     in the simulation between the time when the first driver appeared
33 33
     and the final time period. "unif" means that it is almost sure that
34 34
     different individuals will be sampled at different times. "last"
35
-    does not guarantee that different individuals will be sample at the
35
+    does not guarantee that different individuals will be sampled at the
36 36
     same time unit, only that all will be sampled in the last time unit
37 37
     of their simulation.
38 38
   }
... ...
@@ -18,7 +18,8 @@ test_that("initMutant crashes",
18 18
                       initSize = 1000,
19 19
                       keepPhylog = TRUE
20 20
                      , initMutant = c("b, a") ## also "a" and "b" separate
21
-                                          ))
21
+                      ),
22
+                      "genes not in the fitness table")
22 23
               expect_error(oncoSimulIndiv(o3, model = "McFL",
23 24
                       mu = 5e-5, finalTime = 500,
24 25
                       detectionDrivers = 3,
... ...
@@ -26,7 +27,8 @@ test_that("initMutant crashes",
26 27
                       initSize = 1000,
27 28
                       keepPhylog = TRUE
28 29
                      , initMutant = c("b", "a") ## also "a" and "b" separate
29
-                                          ))
30
+                      ),
31
+                      "genes not in the fitness table")
30 32
               expect_error(oncoSimulIndiv(o3, model = "McFL",
31 33
                       mu = 5e-5, finalTime = 500,
32 34
                       detectionDrivers = 3,
... ...
@@ -34,7 +36,8 @@ test_that("initMutant crashes",
34 36
                       initSize = 1000,
35 37
                       keepPhylog = TRUE
36 38
                      , initMutant = c("1", "2") ## also "a" and "b" separate
37
-                                          ))
39
+                      ),
40
+                      "genes not in the fitness table")
38 41
           })
39 42
 
40 43
 
... ...
@@ -160,3 +163,186 @@ test_that("initMutant non lexicog order",
160 163
               expect_true( "m, d_u" %in% cn )
161 164
               expect_true( "m, d, f_u" %in% cn )
162 165
           })
166
+
167
+
168
+test_that("initMutant with oncoSimulSample", {
169
+    o3init <- allFitnessEffects(orderEffects = c(
170
+                            "M > D > F" = 0.99,
171
+                            "D > M > F" = 0.2,
172
+                            "D > M"     = 0.1,
173
+                            "M > D"     = 0.9),
174
+                        noIntGenes = c("u" = 0.01, 
175
+                                       "v" = 0.01,
176
+                                       "w" = 0.001,
177
+                                       "x" = 0.0001,
178
+                                       "y" = -0.0001,
179
+                                       "z" = -0.001),
180
+                        geneToModule =
181
+                            c("Root" = "Root",
182
+                              "M" = "m",
183
+                              "F" = "f",
184
+                              "D" = "d") )
185
+    ossI <- oncoSimulSample(4, 
186
+                        o3init, model = "Exp",
187
+                        mu = 5e-5, finalTime = 500,
188
+                        detectionDrivers = 2,
189
+                        onlyCancer = TRUE,
190
+                        initSize = 10,
191
+                        initMutant = c("z > d"),
192
+                        thresholdWhole = 1 ## check presence of initMutant
193
+                        )
194
+    ssp <- ossI$popSample
195
+    expect_equal(ssp[, c("d", "z")],
196
+                 matrix(1, nrow = 4, ncol = 2,
197
+                        dimnames = list(NULL, c("d", "z"))))
198
+    expect_false(sum(ssp) == prod(dim(ssp))) ## we don't just have all of
199
+                                             ## them, which would turn the
200
+                                             ## previous into irrelevant
201
+    expect_equal(grep("d",
202
+                      as.character(ossI$popSummary$OccurringDrivers)), 1:4)
203
+})
204
+
205
+
206
+test_that("initMutant with oncoSimulSample, 2", {
207
+    o3init <- allFitnessEffects(orderEffects = c(
208
+                            "M > D > F" = 0.99,
209
+                            "D > M > F" = 0.2,
210
+                            "D > M"     = 0.1,
211
+                            "M > D"     = 0.9,
212
+                            "M > A"     = 0.25),
213
+                        noIntGenes = c("u" = 0.01, 
214
+                                       "v" = 0.01,
215
+                                       "w" = 0.001,
216
+                                       "x" = 0.0001,
217
+                                       "y" = -0.0001,
218
+                                       "z" = -0.001),
219
+                        geneToModule =
220
+                            c("Root" = "Root",
221
+                              "A" = "a",
222
+                              "M" = "m",
223
+                              "F" = "f",
224
+                              "D" = "d") )
225
+    ossI <- oncoSimulSample(4, 
226
+                        o3init, model = "Exp",
227
+                        mu = 5e-5, finalTime = 500,
228
+                        detectionDrivers = 3,
229
+                        onlyCancer = TRUE,
230
+                        initSize = 10,
231
+                        initMutant = c("z > a"),
232
+                        thresholdWhole = 1 ## check presence of initMutant
233
+                        )
234
+    ssp <- ossI$popSample
235
+    expect_equal(ssp[, c("a", "z")],
236
+                 matrix(1, nrow = 4, ncol = 2,
237
+                        dimnames = list(NULL, c("a", "z"))))
238
+    expect_false(sum(ssp) == prod(dim(ssp))) ## we don't just have all of
239
+                                             ## them, which would turn the
240
+                                             ## previous into irrelevant
241
+    expect_equal(grep("a",
242
+                      as.character(ossI$popSummary$OccurringDrivers)), 1:4)
243
+})
244
+
245
+
246
+test_that("initMutant with oncoSimulPop", {
247
+    o3init <- allFitnessEffects(orderEffects = c(
248
+                            "M > D > F" = 0.99,
249
+                            "D > M > F" = 0.2,
250
+                            "D > M"     = 0.1,
251
+                            "M > D"     = 0.9),
252
+                        noIntGenes = c("u" = 0.01, 
253
+                                       "v" = 0.01,
254
+                                       "w" = 0.001,
255
+                                       "x" = 0.0001,
256
+                                       "y" = -0.0001,
257
+                                       "z" = -0.001),
258
+                        geneToModule =
259
+                            c("Root" = "Root",
260
+                              "M" = "m",
261
+                              "F" = "f",
262
+                              "D" = "d") )
263
+    ospI <- oncoSimulPop(4, 
264
+                        o3init, model = "Exp",
265
+                        mu = 5e-5, finalTime = 1000,
266
+                        detectionDrivers = 3,
267
+                        onlyCancer = TRUE,
268
+                        keepPhylog = TRUE,
269
+                        initSize = 10,
270
+                        initMutant = c("d > m > y"),
271
+                        mc.cores = 2
272
+                        )
273
+    genepos <- match(c("d", "m", "y"), ospI[[1]]$geneNames)
274
+    expect_true( all(
275
+        unlist(lapply(ospI,
276
+                      function(x) x$Genotypes[genepos, , drop = FALSE])) == 1))
277
+    ## make sure not all 1, by error, which would render previous useless
278
+    expect_false( all( 
279
+        unlist(lapply(ospI,
280
+                      function(x) x$Genotypes)) == 1))
281
+    expect_true(all(unlist(lapply(ospI,
282
+                                  function(x) (
283
+                                      lapply(x$GenotypesWDistinctOrderEff,
284
+                                             function(z) genepos %in% z))))))
285
+    ## Like before, check no silly things
286
+    allgenepos <- seq_along(ospI[[1]]$geneNames)
287
+    expect_false(all(unlist(lapply(ospI,
288
+                                  function(x) (
289
+                                      lapply(x$GenotypesWDistinctOrderEff,
290
+                                             function(z) allgenepos %in% z))))))
291
+    expect_true(all(
292
+        lapply(ospI,
293
+               function(x)
294
+                   as.character(x$other$PhylogDF[1, 1])) == "d, m_y"))
295
+})
296
+
297
+
298
+
299
+test_that("initMutant with oncoSimulPop, 2", {
300
+    o3init <- allFitnessEffects(orderEffects = c(
301
+                            "M > D > F" = 0.99,
302
+                            "D > M > F" = 0.2,
303
+                            "D > M"     = 0.1,
304
+                            "M > D"     = 0.9),
305
+                        noIntGenes = c("u" = 0.01, 
306
+                                       "v" = 0.01,
307
+                                       "w" = 0.001,
308
+                                       "x" = 0.0001,
309
+                                       "y" = -0.0001,
310
+                                       "z" = -0.001),
311
+                        geneToModule =
312
+                            c("Root" = "Root",
313
+                              "M" = "m",
314
+                              "F" = "f",
315
+                              "D" = "d") )
316
+    ospI <- oncoSimulPop(4, 
317
+                        o3init, model = "Exp",
318
+                        mu = 5e-5, finalTime = 70,
319
+                        detectionDrivers = 4, ## yes, reach end
320
+                        onlyCancer = FALSE,
321
+                        keepPhylog = TRUE,
322
+                        initSize = 10,
323
+                        initMutant = c("m > v > d"),
324
+                        mc.cores = 2
325
+                        )
326
+    genepos <- match(c("d", "m", "v"), ospI[[1]]$geneNames)
327
+    expect_true( all(
328
+        unlist(lapply(ospI,
329
+                      function(x) x$Genotypes[genepos, , drop = FALSE])) == 1))
330
+    ## make sure not all 1, by error, which would render previous useless
331
+    expect_false( all( 
332
+        unlist(lapply(ospI,
333
+                      function(x) x$Genotypes)) == 1))
334
+    expect_true(all(unlist(lapply(ospI,
335
+                                  function(x) (
336
+                                      lapply(x$GenotypesWDistinctOrderEff,
337
+                                             function(z) genepos %in% z))))))
338
+    ## Like before, check no silly things
339
+    allgenepos <- seq_along(ospI[[1]]$geneNames)
340
+    expect_false(all(unlist(lapply(ospI,
341
+                                  function(x) (
342
+                                      lapply(x$GenotypesWDistinctOrderEff,
343
+                                             function(z) allgenepos %in% z))))))
344
+    expect_true(all(
345
+        lapply(ospI,
346
+               function(x)
347
+                   as.character(x$other$PhylogDF[1, 1])) == "m, d_v"))
348
+})
... ...
@@ -3323,21 +3323,65 @@ o3init <- allFitnessEffects(orderEffects = c(
3323 3323
                             "D > M > F" = 0.2,
3324 3324
                             "D > M"     = 0.1,
3325 3325
                             "M > D"     = 0.9),
3326
-                        noIntGenes = c("u" = 0.01, "z" = 0.01),
3326
+                        noIntGenes = c("u" = 0.01, 
3327
+                                       "v" = 0.01,
3328
+                                       "w" = 0.001,
3329
+                                       "x" = 0.0001,
3330
+                                       "y" = -0.0001,
3331
+                                       "z" = -0.001),
3327 3332
                         geneToModule =
3328 3333
                             c("Root" = "Root",
3329 3334
                               "M" = "m",
3330 3335
                               "F" = "f",
3331 3336
                               "D" = "d") )
3332
-tmpI <- oncoSimulIndiv(o3init, model = "McFL",
3337
+
3338
+oneI <- oncoSimulIndiv(o3init, model = "McFL",
3333 3339
                        mu = 5e-5, finalTime = 500,
3334 3340
                        detectionDrivers = 3,
3335 3341
                        onlyCancer = FALSE,
3336 3342
                        initSize = 1000,
3337 3343
                        keepPhylog = TRUE,
3338 3344
                        initMutant = c("m > u > d")
3339
-                      )
3340
-plotClonePhylog(tmpI, N = 0)
3345
+                       )
3346
+plotClonePhylog(oneI, N = 0)
3347
+
3348
+
3349
+## 
3350
+ospI <- oncoSimulPop(4, 
3351
+                    o3init, model = "Exp",
3352
+                    mu = 5e-5, finalTime = 500,
3353
+                    detectionDrivers = 3,
3354
+                    onlyCancer = TRUE,
3355
+                    initSize = 10,
3356
+                    keepPhylog = TRUE,
3357
+                    initMutant = c("d > m > z"),
3358
+                    mc.cores = 2
3359
+                    )
3360
+
3361
+op <- par(mar = rep(0, 4), mfrow = c(2, 2))
3362
+plotClonePhylog(ospI[[1]])
3363
+plotClonePhylog(ospI[[2]])
3364
+plotClonePhylog(ospI[[3]])
3365
+plotClonePhylog(ospI[[4]])
3366
+par(op)
3367
+
3368
+
3369
+ossI <- oncoSimulSample(4, 
3370
+                        o3init, model = "Exp",
3371
+                        mu = 5e-5, finalTime = 500,
3372
+                        detectionDrivers = 2,
3373
+                        onlyCancer = TRUE,
3374
+                        initSize = 10,
3375
+                        initMutant = c("z > d"),
3376
+                        thresholdWhole = 1 ## check presence of initMutant
3377
+                    )
3378
+
3379
+## No phylogeny is kept with oncoSimulSample, but look at the 
3380
+## OcurringDrivers and the sample
3381
+
3382
+ossI$popSample
3383
+ossI$popSummary[, "OccurringDrivers", drop = FALSE]
3384
+
3341 3385
 
3342 3386
 @ 
3343 3387
 
... ...
@@ -1,15 +1,15 @@
1 1
 \usepackage[%
2
-		shash={5ea2efc},
3
-		lhash={5ea2efc59f87434290924e8cd341fac5f5ee67ae},
2
+		shash={ea5bee4},
3
+		lhash={ea5bee4855766a04a04786475d269790a28c8373},
4 4
 		authname={ramon diaz-uriarte (at Bufo)},
5 5
 		authemail={rdiaz02@gmail.com},
6
-		authsdate={2015-09-26},
7
-		authidate={2015-09-26 14:52:30 +0200},
8
-		authudate={1443271950},
6
+		authsdate={2015-09-27},
7
+		authidate={2015-09-27 14:25:22 +0200},
8
+		authudate={1443356722},
9 9
 		commname={ramon diaz-uriarte (at Bufo)},
10 10
 		commemail={rdiaz02@gmail.com},
11
-		commsdate={2015-09-26},
12
-		commidate={2015-09-26 14:52:30 +0200},
13
-		commudate={1443271950},
11
+		commsdate={2015-09-27},
12
+		commidate={2015-09-27 14:25:22 +0200},
13
+		commudate={1443356722},
14 14
 		refnames={ (HEAD -> master)}
15 15
 	]{gitsetinfo}
16 16
\ No newline at end of file