git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/OncoSimulR@118947 bc3139a8-67e5-0310-9ffc-ced21a209358
... | ... |
@@ -1,7 +1,7 @@ |
1 | 1 |
Package: OncoSimulR |
2 | 2 |
Type: Package |
3 | 3 |
Title: Forward Genetic Simulation of Cancer Progresion with Epistasis |
4 |
-Version: 2.3.5 |
|
4 |
+Version: 2.3.6 |
|
5 | 5 |
Date: 2016-06-25 |
6 | 6 |
Authors@R: c(person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"), |
7 | 7 |
email = "rdiaz02@gmail.com"), |
... | ... |
@@ -224,10 +224,14 @@ allGenotypes_to_matrix <- function(x) { |
224 | 224 |
} |
225 | 225 |
|
226 | 226 |
|
227 |
- |
|
228 |
-magellan_stats <- function(x, max_num_genotypes = 2000, |
|
227 |
+## Magellan_stats and Magellan_draw cannot be tested |
|
228 |
+## routinely, as they depend on external software |
|
229 |
+Magellan_stats <- function(x, max_num_genotypes = 2000, |
|
229 | 230 |
verbose = FALSE, |
230 |
- fl_statistics = "fl_statistics") { |
|
231 |
+ use_log = TRUE, |
|
232 |
+ short = TRUE, |
|
233 |
+ fl_statistics = "fl_statistics") { # nocov start |
|
234 |
+ ## I always use |
|
231 | 235 |
## if(!is.null(x) && is.null(file)) |
232 | 236 |
## stop("one of object or file name") |
233 | 237 |
## if(is.null(file)) |
... | ... |
@@ -235,30 +239,77 @@ magellan_stats <- function(x, max_num_genotypes = 2000, |
235 | 239 |
fnret <- tempfile() |
236 | 240 |
if(verbose) |
237 | 241 |
cat("\n Using input file", fn, " and output file ", fnret, "\n") |
242 |
+ |
|
243 |
+ if(use_log) { |
|
244 |
+ logarg <- "-l" |
|
245 |
+ } else { |
|
246 |
+ logarg <- NULL |
|
247 |
+ } |
|
248 |
+ if(short) { |
|
249 |
+ shortarg <- "-s" |
|
250 |
+ } else { |
|
251 |
+ shortarg <- NULL |
|
252 |
+ } |
|
253 |
+ |
|
238 | 254 |
to_Magellan(x, fn, max_num_genotypes = max_num_genotypes) |
239 |
- call_M <- system(paste(fl_statistics, fn, "-s", "-o", fnret)) |
|
240 |
- return(read.table(fnret, skip = 1, header = TRUE)[-1]) |
|
241 |
-} |
|
255 |
+ call_M <- system2(fl_statistics, args = paste(fn, shortarg, logarg, "-o", fnret)) |
|
256 |
+ if(short) { |
|
257 |
+ tmp <- as.vector(read.table(fnret, skip = 1, header = TRUE)[-1]) |
|
258 |
+ ## Make names more explicit, but check we have what we think we have |
|
259 |
+ stopifnot(length(tmp) == 23) |
|
260 |
+ stopifnot(identical(names(tmp), |
|
261 |
+ c("ngeno", "npeaks", "nsinks", "gamma", "gamma.", "r.s", |
|
262 |
+ "nchains", "nsteps", "nori", "depth", "magn", "sign", |
|
263 |
+ "rsign", "w.1.", "w.2.", "w.3..", "mode_w", "s.1.", |
|
264 |
+ "s.2.", "s.3..", "mode_s", "pairs_s", "outD_v"))) |
|
265 |
+ names(tmp) <- c("n_genotypes", "n_peaks", "n_sinks", "gamma", "gamma_star", |
|
266 |
+ "r/s","n_chains", "n_steps", "n_origins", "max_depth", |
|
267 |
+ "epist_magn", "epist_sign", "epist_rsign", |
|
268 |
+ "walsh_coef_1", "walsh_coef_2", "walsh_coef_3", "mode_walsh", |
|
269 |
+ "synerg_coef_1", "synerg_coef_2", "synerg_coef_3", "mode_synerg", |
|
270 |
+ "std_dev_pairs", "var_outdegree") |
|
271 |
+ } else { |
|
272 |
+ tmp <- readLines(fnret) |
|
273 |
+ } |
|
274 |
+ return(tmp) |
|
275 |
+} # nocov end |
|
276 |
+ |
|
277 |
+Magellan_draw <- function(x, max_num_genotypes = 2000, |
|
278 |
+ verbose = FALSE, |
|
279 |
+ args = "-f", |
|
280 |
+ fl_draw = "fl_draw", |
|
281 |
+ svg_open = "xdg-open", |
|
282 |
+ file_name = NULL) { # nocov start |
|
283 |
+ ## It always works by appending the name so file_name is without the .svg |
|
284 |
+ if(is.null(file_name)) { |
|
285 |
+ fn <- tempfile() |
|
286 |
+ } else { |
|
287 |
+ fn <- file_name |
|
288 |
+ } |
|
289 |
+ fn_out <- paste0(fn, ".svg") |
|
290 |
+ if(verbose) |
|
291 |
+ cat("\n Using input file", fn, " and output file ", fn_out, "\n") |
|
292 |
+ |
|
293 |
+ to_Magellan(x, fn, max_num_genotypes = max_num_genotypes) |
|
294 |
+ call_M <- system2(fl_draw, args = paste(fn, args), wait = FALSE) |
|
295 |
+ call_view <- system2(svg_open, args = fn_out, wait = FALSE, |
|
296 |
+ stdout = ifelse(verbose, "", FALSE), |
|
297 |
+ stderr = ifelse(verbose, "", FALSE)) |
|
298 |
+ |
|
299 |
+ invisible() |
|
300 |
+} # nocov end |
|
242 | 301 |
|
243 | 302 |
|
244 |
-## |
|
245 | 303 |
|
246 |
-## ## Example of Bozic issues |
|
304 |
+## ## Example of Bozic issues in conversions of fitnes |
|
247 | 305 |
## m1 <- cbind(c(0, 1), c(1, 0), c(2, 3)) |
248 |
- |
|
249 | 306 |
## m2 <- cbind(c(1, 1), c(1, 0), c(2, 3)) |
250 |
- |
|
251 | 307 |
## m3 <- data.frame(G = c("A, B", "A"), F = c(1, 2)) |
252 |
- |
|
253 | 308 |
## m4 <- data.frame(G = c("A, B", "A", "WT", "B"), F = c(3, 2, 1, 4)) |
254 |
- |
|
255 | 309 |
## m5 <- data.frame(G = c("A, B", "A", "WT", "B"), F = c(3, 2, 1, 0)) |
256 |
- |
|
257 | 310 |
## m6 <- data.frame(G = c("A, B", "A", "WT", "B"), F = c(3, 2.5, 2, 0)) |
258 | 311 |
|
259 |
- |
|
260 |
- |
|
261 | 312 |
## And no, it makes no sense to use any of this for mutator: in mutator I |
262 | 313 |
## directly have the multiplication factor of each gene. Which is likely |
263 |
-## what people want anyway. Add it later if needed. by using a ratio |
|
314 |
+## what people want anyway. Add it later if needed by using a ratio |
|
264 | 315 |
## instead of a "-" |
... | ... |
@@ -122,6 +122,11 @@ Ramon Diaz-Uriarte |
122 | 122 |
epistasis stats for the landscape, and several visual manipulation |
123 | 123 |
options. |
124 | 124 |
|
125 |
+ One feature of this function that is not available in MAGELLAN is |
|
126 |
+ showing genotype labels (i.e., annotated by gene names), which can be |
|
127 |
+ helpful if the different genotypes mean something to you. |
|
128 |
+ |
|
129 |
+ |
|
125 | 130 |
In addition to the above differences, another difference between this |
126 | 131 |
plot and those of MAGELLAN is \bold{how sinks/peaks of more than one |
127 | 132 |
genotype are dealt with}. This plot will show as sinks or peaks sets |
... | ... |
@@ -143,7 +148,7 @@ Ramon Diaz-Uriarte |
143 | 148 |
possibility of peaks/sinks made of more than one genotype actually |
144 | 149 |
makes code much simpler. |
145 | 150 |
|
146 |
- |
|
151 |
+ |
|
147 | 152 |
Sometimes not showing the any links that involve a decrease in fitness |
148 | 153 |
can help see non-accessible pathways (in strong selection, no multiple |
149 | 154 |
mutations, etc); do this by passing, for instance, an NA for the |
150 | 155 |
new file mode 100644 |
... | ... |
@@ -0,0 +1,10 @@ |
1 |
+test_that("Expect output", { |
|
2 |
+ r1 <- rfitness(4) |
|
3 |
+ expect_output(to_Magellan(r1, NULL)) |
|
4 |
+ cs <- data.frame(parent = c(rep("Root", 3), "a", "d", "c"), |
|
5 |
+ child = c("a", "b", "d", "e", "c", "f"), |
|
6 |
+ s = 0.1, |
|
7 |
+ sh = -0.9, |
|
8 |
+ typeDep = "MN") |
|
9 |
+ expect_output(to_Magellan(allFitnessEffects(cs), NULL)) |
|
10 |
+}) |
... | ... |
@@ -36,3 +36,63 @@ test_that("We fail with order", { |
36 | 36 |
test_that("We fail with wrong separator, :", { |
37 | 37 |
expect_error(OncoSimulR:::from_genotype_fitness(data.frame(G = "A : B", F = 1))) |
38 | 38 |
}) |
39 |
+ |
|
40 |
+ |
|
41 |
+test_that("Dividing fitness by wt and missing genot is 1", { |
|
42 |
+ |
|
43 |
+ m7 <- cbind(c(1, 1, 0), c(1, 0, 0), c(2, 3, 5)) |
|
44 |
+ expect_message(fem7 <- allFitnessEffects(genotFitness = m7), |
|
45 |
+ "Fitness of wildtype != 1. Dividing all fitnesses by fitness of wildtype.", |
|
46 |
+ fixed = TRUE) |
|
47 |
+ expect_equivalent(evalAllGenotypes(fem7, order = FALSE, addwt = TRUE)[, 2], |
|
48 |
+ c(1, 3/5, 1, 2/5)) |
|
49 |
+}) |
|
50 |
+ |
|
51 |
+test_that("The WT is added if absent", { |
|
52 |
+ m7 <- cbind(c(1, 1, 0), c(1, 0, 0), c(2, 3, 5)) |
|
53 |
+ fem7 <- allFitnessEffects(genotFitness = m7) |
|
54 |
+ ag <- evalAllGenotypes(fem7, order = FALSE) |
|
55 |
+ ## internal call |
|
56 |
+ expect_equivalent(OncoSimulR:::allGenotypes_to_matrix(ag)[, 3], |
|
57 |
+ c(1, 3/5, 1, 2/5)) |
|
58 |
+ ## the user visible, which is via plotFitnessLandscape -> to_Fitness_Matrix |
|
59 |
+ plot(ag) |
|
60 |
+}) |
|
61 |
+ |
|
62 |
+ |
|
63 |
+test_that("genotFitness not combined with others", { |
|
64 |
+ |
|
65 |
+ m7 <- cbind(c(1, 1, 0), c(1, 0, 0), c(2, 3, 5)) |
|
66 |
+ |
|
67 |
+ expect_error(allFitnessEffects(genotFitness = m7, |
|
68 |
+ noIntGenes = c(.1, .2)), |
|
69 |
+ "You have a non-null genotFitness.", |
|
70 |
+ fixed = TRUE) |
|
71 |
+ |
|
72 |
+ expect_error(allFitnessEffects(genotFitness = m7, |
|
73 |
+ epistasis = c("A" = .3, |
|
74 |
+ "B" = .1, |
|
75 |
+ "C" = .2)), |
|
76 |
+ "You have a non-null genotFitness.", |
|
77 |
+ fixed = TRUE) |
|
78 |
+ |
|
79 |
+ expect_error(allFitnessEffects(genotFitness = m7, |
|
80 |
+ orderEffects = c("G > H" = 1.1)), |
|
81 |
+ "You have a non-null genotFitness.", |
|
82 |
+ fixed = TRUE) |
|
83 |
+ |
|
84 |
+ |
|
85 |
+ p3 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c", "f"), |
|
86 |
+ child = c("a", "b", "d", "e", "c", "c", "f", "f", "g", "g"), |
|
87 |
+ s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3), |
|
88 |
+ sh = c(rep(0, 4), c(-.9, -.9), c(-.95, -.95), c(-.99, -.99)), |
|
89 |
+ typeDep = c(rep("--", 4), |
|
90 |
+ "XMPN", "XMPN", "MN", "MN", "SM", "SM")) |
|
91 |
+ |
|
92 |
+ expect_error(allFitnessEffects(rT = p3, genotFitness = m7), |
|
93 |
+ "You have a non-null genotFitness.", |
|
94 |
+ fixed = TRUE) |
|
95 |
+ |
|
96 |
+}) |
|
97 |
+ |
|
98 |
+ |
... | ... |
@@ -210,6 +210,35 @@ test_that("eval mut genotypes", { |
210 | 210 |
|
211 | 211 |
|
212 | 212 |
|
213 |
+test_that("eval mut genotypes, echo", { |
|
214 |
+ fe <- allFitnessEffects(epistasis = c("a : b" = 0.3, |
|
215 |
+ "b : c" = 0.5), |
|
216 |
+ noIntGenes = c("e" = 0.1), |
|
217 |
+ drvNames = c(letters[1:3])) |
|
218 |
+ fm <- allMutatorEffects(noIntGenes = c("a" = 10, |
|
219 |
+ "c" = 5)) |
|
220 |
+ expect_output(evalGenotypeMut("a, c", fm, echo = TRUE), |
|
221 |
+ "Mutation rate product", fixed = TRUE) |
|
222 |
+}) |
|
223 |
+ |
|
224 |
+ |
|
225 |
+test_that("evaluating genotype and mutator, Bozic", { |
|
226 |
+ fe <- allFitnessEffects(epistasis = c("a : b" = 0.3, |
|
227 |
+ "b : c" = 0.5), |
|
228 |
+ drvNames = letters[1:3]) |
|
229 |
+ fm <- allMutatorEffects(noIntGenes = c("a" = 10, |
|
230 |
+ "c" = 5)) |
|
231 |
+ ou <- evalAllGenotypesFitAndMut(fe, fm, order = FALSE, |
|
232 |
+ model = "Bozic") |
|
233 |
+ expect_equivalent(ou[, 2], |
|
234 |
+ c(1, 1, 1, 1 - .3, 1, 1 -.5, (1 -.3) * (1 - .5)) |
|
235 |
+ ) |
|
236 |
+ expect_equivalent(ou[, 3], |
|
237 |
+ c(10, 1, 5, 10, 10 * 5, 5, 10 * 5) |
|
238 |
+ ) |
|
239 |
+}) |
|
240 |
+ |
|
241 |
+ |
|
213 | 242 |
test_that("mut and fitness both needed when needed", { |
214 | 243 |
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3, |
215 | 244 |
"b : c" = 0.5), |
... | ... |
@@ -1,6 +1,10 @@ |
1 | 1 |
test_that("Exercise plotting", { |
2 | 2 |
r1 <- rfitness(4) |
3 | 3 |
expect_silent(plot(r1)) |
4 |
+ expect_silent(plot(r1, log = TRUE)) |
|
5 |
+ expect_silent(plot(r1, log = TRUE, use_ggrepel = TRUE)) |
|
6 |
+ expect_silent(plot(r1, log = TRUE, show_labels = FALSE)) |
|
7 |
+ |
|
4 | 8 |
|
5 | 9 |
## Specify fitness in a matrix, and plot it |
6 | 10 |
m5 <- cbind(c(0, 1, 0, 1), c(0, 0, 1, 1), c(1, 2, 3, 5.5)) |
... | ... |
@@ -12,14 +16,54 @@ test_that("Exercise plotting", { |
12 | 16 |
noIntGenes = c("e" = 0.1)) |
13 | 17 |
|
14 | 18 |
expect_silent(plot(evalAllGenotypes(fe, order = FALSE))) |
15 |
- |
|
19 |
+ |
|
16 | 20 |
## same as |
17 | 21 |
expect_silent(plotFitnessLandscape(evalAllGenotypes(fe, order = FALSE))) |
18 |
- |
|
22 |
+ ## more ggrepel |
|
23 |
+ expect_silent(plot(evalAllGenotypes(fe, order = FALSE), use_ggrepel = TRUE)) |
|
24 |
+}) |
|
25 |
+ |
|
26 |
+ |
|
27 |
+test_that("to_FitnessMatrix stops as it should", { |
|
28 |
+ x1 <- data.frame(a = 1:2, b = 1:2) |
|
29 |
+ expect_error(OncoSimulR:::to_Fitness_Matrix(x1, 2000), |
|
30 |
+ "We cannot guess what you are passing", |
|
31 |
+ fixed = TRUE) |
|
32 |
+ x2 <- list(a = 12, b = 13) |
|
33 |
+ expect_error(OncoSimulR:::to_Fitness_Matrix(x2, 2000), |
|
34 |
+ "We cannot guess what you are passing", |
|
35 |
+ fixed = TRUE) |
|
19 | 36 |
}) |
20 | 37 |
|
21 | 38 |
|
22 | 39 |
|
40 |
+test_that("to_FitnessMatrix can deal with df", { |
|
41 |
+ m4 <- data.frame(G = c("A, B", "A", "WT", "B"), |
|
42 |
+ Fitness = c(3, 2, 1, 4)) |
|
43 |
+ expect_message(OncoSimulR:::to_Fitness_Matrix(m4, 2000), |
|
44 |
+ "Column names of object", fixed = TRUE) |
|
45 |
+ m5 <- data.frame(G = c("A, B", "B"), |
|
46 |
+ Fitness = c(3, 2)) |
|
47 |
+ expect_message(OncoSimulR:::to_Fitness_Matrix(m5, 2000), |
|
48 |
+ "Column names of object", fixed = TRUE) |
|
49 |
+ x1 <- data.frame(a = c("A, B"), Fitness = 2) |
|
50 |
+ expect_message(OncoSimulR:::to_Fitness_Matrix(x1, 2000), |
|
51 |
+ "Column names of object", fixed = TRUE) |
|
52 |
+ x2 <- data.frame(a = c("A, B", "B"), Fitness = c(2, 3)) |
|
53 |
+ expect_message(OncoSimulR:::to_Fitness_Matrix(x2, 2000), |
|
54 |
+ "Column names of object", fixed = TRUE) |
|
55 |
+ x3 <- data.frame(a = c("A, B", "C"), Fitness = c(2, 3)) |
|
56 |
+ expect_message(OncoSimulR:::to_Fitness_Matrix(x3, 2000), |
|
57 |
+ "Column names of object", fixed = TRUE) |
|
58 |
+ ## Now, the user code |
|
59 |
+ expect_message(plotFitnessLandscape(x1)) |
|
60 |
+ expect_message(plotFitnessLandscape(x2)) |
|
61 |
+ expect_message(plotFitnessLandscape(x3)) |
|
62 |
+ expect_message(plotFitnessLandscape(m5)) |
|
63 |
+ expect_message(plotFitnessLandscape(m4)) |
|
64 |
+}) |
|
65 |
+ |
|
66 |
+ |
|
23 | 67 |
test_that("internal peak valley functions", { |
24 | 68 |
|
25 | 69 |
x <- matrix(NA, 14, 14) |
... | ... |
@@ -387,5 +387,22 @@ test_that("exercising sampling code, customSize", { |
387 | 387 |
|
388 | 388 |
|
389 | 389 |
|
390 |
+test_that("exercising sampling code, customSize", { |
|
391 |
+ cs <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), |
|
392 |
+ child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), |
|
393 |
+ s = 0.1, |
|
394 |
+ sh = -0.9, |
|
395 |
+ typeDep = "MN") |
|
396 |
+ cbn1 <- allFitnessEffects(cs, drvNames = c("a", "b", "c", "d", "e", "g")) |
|
397 |
+ o1 <- oncoSimulIndiv(cbn1, detectionSize = 1e8, |
|
398 |
+ onlyCancer = TRUE, |
|
399 |
+ detectionDrivers = 1, |
|
400 |
+ max.num.tries = 5000, |
|
401 |
+ sampleEvery = 0.03, |
|
402 |
+ keepEvery = 1000000) |
|
403 |
+ expect_message(samplePop(o1, timeSample = "unif"), |
|
404 |
+ "Only one sampled time period with mutants", fixed = TRUE) |
|
405 |
+}) |
|
406 |
+ |
|
390 | 407 |
|
391 | 408 |
cat(paste("\n Ending samplePop tests", date(), "\n")) |
... | ... |
@@ -1,15 +1,15 @@ |
1 | 1 |
\usepackage[% |
2 |
- shash={a8daef6}, |
|
3 |
- lhash={a8daef6cae6f3c218599ea59264b950e95040f54}, |
|
4 |
- authname={Ramon Diaz-Uriarte}, |
|
5 |
- authemail={rdiaz02@users.noreply.github.com}, |
|
6 |
- authsdate={2016-06-24}, |
|
7 |
- authidate={2016-06-24 10:14:31 +0200}, |
|
8 |
- authudate={1466756071}, |
|
9 |
- commname={Ramon Diaz-Uriarte}, |
|
10 |
- commemail={rdiaz02@users.noreply.github.com}, |
|
11 |
- commsdate={2016-06-24}, |
|
12 |
- commidate={2016-06-24 10:14:31 +0200}, |
|
13 |
- commudate={1466756071}, |
|
14 |
- refnames={ (HEAD, origin/master, origin/HEAD)} |
|
2 |
+ shash={3f0ca61}, |
|
3 |
+ lhash={3f0ca612b9bf91cbc9e99fcaa84a03de2fb1446d}, |
|
4 |
+ authname={ramon diaz-uriarte (at Bufo)}, |
|
5 |
+ authemail={rdiaz02@gmail.com}, |
|
6 |
+ authsdate={2016-06-25}, |
|
7 |
+ authidate={2016-06-25 14:02:01 +0200}, |
|
8 |
+ authudate={1466856121}, |
|
9 |
+ commname={ramon diaz-uriarte (at Bufo)}, |
|
10 |
+ commemail={rdiaz02@gmail.com}, |
|
11 |
+ commsdate={2016-06-25}, |
|
12 |
+ commidate={2016-06-25 14:02:01 +0200}, |
|
13 |
+ commudate={1466856121}, |
|
14 |
+ refnames={ (HEAD -> master, origin/master, origin/HEAD)} |
|
15 | 15 |
]{gitsetinfo} |
16 | 16 |
\ No newline at end of file |