inittime <- Sys.time()
### Since some tests are slow and some tests are very fragile, for now I
### leave date() and seed()
## Uses same code as test.mutator, but with fitness landscape
## specification. So we use this convenience since we had
## always an allFitnessEffects object from a different specification
## but with many genes cannot be run as looooooots of genotypes
fe2fl <- function(x) {
allFitnessEffects(
genotFitness =
OncoSimulR:::allGenotypes_to_matrix(
evalAllGenotypes(x, addwt = TRUE)))
}
### This test takes about 10 seconds
cat(paste("\n Starting test.flfast-mutator.R test at", date()))
date()
## RNGkind("L'Ecuyer-CMRG") ## for the mclapplies
mutsPerClone <- function(x, per.pop.mean = TRUE) {
perCl <- function(z)
unlist(lapply(z$GenotypesWDistinctOrderEff, length))
perCl2 <- function(z)
mean(unlist(lapply(z$GenotypesWDistinctOrderEff, length)))
if(per.pop.mean)
unlist(lapply(x, function(u) perCl2(u)))
else
lapply(x, function(u) perCl(u))
}
## Minifunctions for testing
medianNClones <- function(x) {
median(summary(x)$NumClones)
}
NClones <- function(x) {
summary(x)$NumClones
}
enom <- function(name, mu, ni = no, pp = pops) {
## expected without a given name for init
ii <- which(names(mu) == name)
out <- ni * pp * mu[-ii]
out[order(names(out))]
}
pnom <- function(name, mu, ni = no, pp = pops) {
ee <- enom(name, mu, ni, pp)
ee/sum(ee)
}
snom <- function(name, out) {
## observed without the init
cs <- colSums(OncoSimulR:::geneCounts(out))
ii <- which(names(cs) == name)
cs <- cs[-ii]
cs[order(names(cs))]
}
sm <- function(name, out) {
## totals for a given gene
cs <- colSums(OncoSimulR:::geneCounts(out))
ii <- which(names(cs) == name)
cs[ii]
}
totalind <- function(out) {
## total num indivs
sum(unlist(lapply(out, function(x) x$TotalPopSize)))
}
smA <- function(out) {
## totals counts for all. So total mutated over all.
cs <- colSums(OncoSimulR:::geneCounts(out))
sum(cs)
}
smAPi <- function(out) {
## smAnom but divided by number of individuals.
smA(out)/totalind(out)
}
smAnom <- function(out, name) {
## totals counts for all. So total mutated over all.
## but remove one, the fixed starting one.
cs <- colSums(OncoSimulR:::geneCounts(out))
ii <- which(names(cs) == name)
cs <- cs[-ii]
sum(cs)
}
smAnomPi <- function(out, name) {
## smAnom but divided by number of individuals.
smAnom(out, name)/totalind(out)
}
## The threshold is rather arbitrary. Of course, you would not expect most
## of them to pass if they were false. But note that we cannot set too
## tiny a threshold unless we increase number of repetitions a lot, and
## that makes the test run for very long.
p.value.threshold <- 0.01
date()
test_that("This should not crash", {
## This used to crash because of not nulling the empty mutator effects
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
"b : c" = 0.5),
noIntGenes = c("e" = 0.1),
drvNames = c("a", "b", "c"))
moo <- rep(1e-5, 4)
names(moo) <- c("a", "b", "c", "e")
expect_output(print(oncoSimulIndiv(fe2fl(fe),
mu = moo,
finalTime = 1,
mutationPropGrowth = FALSE,
initSize = 1000,
onlyCancer = FALSE)),
"Individual OncoSimul",
fixed = TRUE)
muvar2 <- c("U" = 1e-5, "z" = 1e-5, "e" = 1e-5, "m" = 1e-5, "D" = 1e-5)
ni1 <- rep(0.02, 5)
names(ni1) <- names(muvar2)
fe1 <- allFitnessEffects(noIntGenes = ni1)
no <- 1e5
reps <- 10
bb <- oncoSimulIndiv(fe2fl(fe1),
mu = muvar2, onlyCancer = FALSE, detectionProb = NA,
initSize = no,
finalTime = 1,
seed = NULL
)
expect_output(print(bb),
"Individual OncoSimul",
fixed = TRUE)
expect_output(print(oncoSimulPop(4,
fe2fl(fe1),
mu = muvar2,
onlyCancer = FALSE, detectionProb = NA,
initSize = 1000,
finalTime = 1,
seed = NULL, mc.cores = 2)),
"Population of OncoSimul",
fixed = TRUE)
expect_output(print(oncoSimulPop(4,
fe2fl(fe),
mu = moo,
onlyCancer = FALSE, detectionProb = NA,
initSize = 1000,
finalTime = 1,
seed = NULL, mc.cores = 2)),
"Population of OncoSimul",
fixed = TRUE)
})
test_that("eval fitness and mut OK", {
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
"b : c" = 0.5),
noIntGenes = c("e" = 0.1))
fm <- allMutatorEffects(noIntGenes = c("a" = 10,
"c" = 5))
expect_output(ou <- evalGenotypeFitAndMut("a", fe2fl(fe), fm, verbose = TRUE),
#"10",
fixed = TRUE)
expect_identical(ou, c(1, 10))
expect_identical(evalGenotypeFitAndMut("b", fe2fl(fe), fm),
c(1, 1))
expect_identical(evalGenotypeFitAndMut("e", fe2fl(fe), fm),
c(1.1, 1))
expect_identical(evalGenotypeFitAndMut("b, e", fe2fl(fe), fm),
c(1.1, 1))
expect_equal(evalGenotypeFitAndMut("a, b, e", fe2fl(fe), fm),
c(1.3 * 1.1, 10))
expect_equal(evalGenotypeFitAndMut("a, b, c, e", fe2fl(fe), fm),
c(1.3 * 1.5 * 1.1, 10 * 5))
})
test_that("expect output oncoSimulIndiv", {
fe <- allFitnessEffects(noIntGenes = c("a" = 0.2,
"c" = 0.4,
"d" = 0.6,
"e" = 0.1))
fm <- allMutatorEffects(noIntGenes = c("a" = 10,
"c" = 5))
expect_output(print(oncoSimulIndiv(fe2fl(fe), muEF = fm, sampleEvery = 0.01,
keepEvery = 5)),
"Individual OncoSimul trajectory",
fixed = TRUE)
expect_output(print(oncoSimulIndiv(fe2fl(fe), muEF = fm, sampleEvery = 0.01,
keepEvery = 5)),
"Individual OncoSimul trajectory",
fixed = TRUE)
})
## test_that("eval mut genotypes", {
## fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
## "b : c" = 0.5),
## noIntGenes = c("e" = 0.1),
## drvNames = c(letters[1:3]))
## fm <- allMutatorEffects(noIntGenes = c("a" = 10,
## "c" = 5))
## expect_identical(evalAllGenotypesMut(fm)[, 2],
## c(10, 5, 50))
## expect_identical(evalGenotypeMut("a", fm),
## 10)
## expect_identical(evalGenotypeMut("c", fm),
## 5)
## expect_identical(evalGenotypeMut("c, a", fm),
## 50)
## expect_identical(evalGenotypeMut("a, c", fm),
## 50)
## expect_identical(evalGenotypeMut("a > c", fm),
## 50)
## expect_identical(evalGenotypeMut("c > a", fm),
## 50)
## expect_error(evalGenotypeMut("b", fm),
## "genotype contains NAs or genes not in fitnessEffects/mutatorEffects",
## fixed = TRUE)
## })
## test_that("eval mut genotypes, echo", {
## fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
## "b : c" = 0.5),
## noIntGenes = c("e" = 0.1),
## drvNames = c(letters[1:3]))
## fm <- allMutatorEffects(noIntGenes = c("a" = 10,
## "c" = 5))
## expect_output(evalGenotypeMut("a, c", fm, echo = TRUE),
## "Mutation rate product", fixed = TRUE)
## })
test_that("evaluating genotype and mutator, Bozic", {
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
"b : c" = 0.5),
drvNames = letters[1:3])
fm <- allMutatorEffects(noIntGenes = c("a" = 10,
"c" = 5))
expect_warning(
ou <- evalAllGenotypesFitAndMut(fe2fl(fe), fm, order = FALSE,
model = "Bozic"),
"Bozic model passing a fitness landscape",
fixed = TRUE)
## expect_equivalent(ou[, 2],
## c(1, 1, 1, 1 - .3, 1, 1 -.5, (1 -.3) * (1 - .5))
## )
## expect_equivalent(ou[, 3],
## c(10, 1, 5, 10, 10 * 5, 5, 10 * 5)
## )
})
test_that("mut and fitness both needed when needed", {
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
"b : c" = 0.5),
noIntGenes = c("e" = 0.1),
drvNames = c(letters[1:3]))
fm <- allMutatorEffects(noIntGenes = c("a" = 10,
"c" = 5))
expect_error(evalGenotypeFitAndMut("a, b", fe2fl(fe)),
'argument "mutatorEffects" is missing',
fixed = TRUE)
expect_error(evalGenotypeFitAndMut("a, b", fm),
"Genotype contains NAs or genes not in fitnessEffects/mutatorEffects",
fixed = TRUE)
expect_error(evalGenotypeFitAndMut("a, b", mutatorEffects = fm),
'argument "fitnessEffects" is missing',
fixed = TRUE)
expect_error(evalAllGenotypesFitAndMut(fe2fl(fe), order = TRUE),
'argument "mutatorEffects" is missing',
fixed = TRUE)
expect_error(evalAllGenotypesFitAndMut(fe2fl(fe), order = FALSE),
'argument "mutatorEffects" is missing',
fixed = TRUE)
expect_error(evalAllGenotypesFitAndMut(fm),
'argument "mutatorEffects" is missing',
fixed = TRUE)
expect_error(evalAllGenotypesFitAndMut(mutatorEffects = fm),
'argument "fitnessEffects" is missing',
fixed = TRUE)
})
test_that("we evaluate the WT", {
## Is fitness of wildtype always 0? Really? Evaluate it.
## It is: see evalGenotypeFitness
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
"b : c" = 0.5),
noIntGenes = c("e" = 0.1))
## expect_warning(
null <- capture.output(
ou <- OncoSimulR:::evalRGenotype(vector(mode = "integer",
length = 0),
fe2fl(fe),
0,
TRUE,
FALSE,
"evalGenotype",
0)## ,
## "WARNING: you have evaluated fitness/mutator status of a genotype of length zero",
## fixed = TRUE)
)
expect_identical(ou, 1)
})
test_that("we evaluate the WT, 2", {
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
"b : c" = 0.5),
noIntGenes = c("e" = 0.1))
fm <- OncoSimulR:::allMutatorEffects(noIntGenes = c("a" = 10,
"c" = 5))
## expect_warning(
null <- capture.output(
ou2 <- OncoSimulR:::evalRGenotypeAndMut(
vector(mode = "integer", length = 0),#rG
fe2fl(fe),#rFE
fm,#muEF
0,#spPop
OncoSimulR:::matchGeneIDs(fm, fe)$Reduced,#fullmutator_
TRUE,#verbose
FALSE,#prodneg
0)## , #currentTime
## "WARNING: you have evaluated fitness of a genotype of length zero.",
## fixed = TRUE)
)
expect_identical(ou2, c(1, 1))
})
test_that("fails if genes in mutator not in fitness", {
fe <- allFitnessEffects(noIntGenes = c("a" = 0.12,
"c" = 0.14,
"d" = 0.16,
"e" = 0.11))
fm4 <- allMutatorEffects(noIntGenes = c("a" = .010,
"b" = .03,
"d" = .08,
"c" = .05))
expect_error(oncoSimulIndiv(fe2fl(fe), muEF = fm4),
"Genes in mutatorEffects not present in fitnessEffects",
fixed = TRUE)
})
test_that("We cannot pass mutator/fitness objects to the wrong functions", {
fe2 <- allFitnessEffects(noIntGenes =
c(a1 = 0.1, a2 = 0.2,
b1 = 0.01, b2 = 0.3, b3 = 0.2,
c1 = 0.3, c2 = -0.2))
fm2 <- allMutatorEffects(epistasis = c("A" = 5,
"B" = 10,
"C" = 3),
geneToModule = c("A" = "a1, a2",
"B" = "b1, b2, b3",
"C" = "c1, c2"))
expect_error(evalAllGenotypesMut(fe2fl(fe2)),
"You are trying to get the mutator effects of a fitness specification.",
fixed = TRUE)
expect_error(evalAllGenotypesMut(fe2fl(fe2)),
"You are trying to get the mutator effects of a fitness specification.",
fixed = TRUE)
expect_error(evalAllGenotypes(fm2),
"You are trying to get the fitness of a mutator specification.",
fixed = TRUE)
expect_error(evalGenotypeMut("a1, b2", fe2fl(fe2)),
"You are trying to get the mutator effects of a fitness specification.",
fixed = TRUE)
expect_error(evalGenotype("a2, c2", fm2),
"You are trying to get the fitness of a mutator specification.",
fixed = TRUE)
})
test_that("evaluating genotype and mutator", {
fe <- allFitnessEffects(epistasis = c("a : b" = 0.3,
"b : c" = 0.5),
noIntGenes = c("e" = 0.1),
drvNames = letters[1:3])
fm <- allMutatorEffects(noIntGenes = c("a" = 10,
"c" = 5))
ou <- evalAllGenotypesFitAndMut(fe2fl(fe), fm, order = FALSE)
expect_equivalent(dplyr::filter(ou, Genotype == "a")[, c(2, 3)],
c(1, 10))
expect_equivalent(dplyr::filter(ou, Genotype == "b")[, c(2, 3)],
c(1, 1))
expect_equivalent(dplyr::filter(ou, Genotype == "c")[, c(2, 3)],
c(1, 5))
expect_equivalent(dplyr::filter(ou, Genotype == "e")[, c(2, 3)],
c(1.1, 1))
expect_equivalent(dplyr::filter(ou, Genotype == "a, b, c")[, c(2, 3)],
c(1.3 * 1.5, 10 * 5))
expect_equivalent(dplyr::filter(ou, Genotype == "b, c, e")[, c(2, 3)],
c(1.5 * 1.1, 5))
expect_equivalent(dplyr::filter(ou, Genotype == "a, b, c, e")[, c(2, 3)],
c(1.3 * 1.5 * 1.1, 10 * 5))
oo <- evalAllGenotypesFitAndMut(fe2fl(fe), fm, order = TRUE)
expect_equivalent(dplyr::filter(oo, Genotype == "a")[, c(2, 3)],
c(1, 10))
expect_equivalent(dplyr::filter(oo, Genotype == "b")[, c(2, 3)],
c(1, 1))
expect_equivalent(dplyr::filter(oo, Genotype == "c")[, c(2, 3)],
c(1, 5))
expect_equivalent(dplyr::filter(oo, Genotype == "e")[, c(2, 3)],
c(1.1, 1))
expect_equivalent(dplyr::filter(oo, Genotype == "a > b > c")[, c(2, 3)],
c(1.3 * 1.5, 10 * 5))
expect_equivalent(dplyr::filter(oo, Genotype == "b > c > e")[, c(2, 3)],
c(1.5 * 1.1, 5))
expect_equivalent(dplyr::filter(oo, Genotype == "e > b > c")[, c(2, 3)],
c(1.5 * 1.1, 5))
expect_equivalent(dplyr::filter(oo, Genotype == "a > b > c > e")[, c(2, 3)],
c(1.3 * 1.5 * 1.1, 10 * 5))
})
date()
test_that("Mutator, several modules differences, fitness eval", {
## the basis of what we do below, but fewer genes here
ln <- 2
m1 <- 5
ni <- rep(0, 3 * ln)
gna <- paste0("a", 1:ln)
gnb <- paste0("b", 1:ln)
gnc <- paste0("c", 1:ln)
names(ni) <- c(gna, gnb, gnc)
gn1 <- paste(c(gna, gnb, gnc), collapse = ", ")
gna <- paste(gna, collapse = ", ")
gnb <- paste(gnb, collapse = ", ")
gnc <- paste(gnc, collapse = ", ")
mut1 <- allMutatorEffects(epistasis = c("A" = m1),
geneToModule = c("A" = gn1))
mut2 <- allMutatorEffects(epistasis = c("A" = m1,
"B" = m1,
"C" = m1),
geneToModule = c("A" = gna,
"B" = gnb,
"C" = gnc))
f1 <- allFitnessEffects(noIntGenes = ni)
e1 <- evalAllGenotypesFitAndMut(fe2fl(f1), mut1, order = FALSE,
addwt = TRUE)
e2 <- evalAllGenotypesFitAndMut(fe2fl(f1), mut2, order = FALSE,
addwt = TRUE)
expect_identical(e1$Fitness, rep(1, 64))
expect_identical(e1$MutatorFactor, c(1, rep(5, 63)))
expect_identical(e2$Fitness, rep(1, 64))
expect_identical(
e2$MutatorFactor,
c(1, rep(5, 7),
rep(25, 8), 5,
rep(25, 4), 5,
rep(25, 5),
rep(125, 4), 25, 25,
rep(125, 4), rep(25, 6),
rep(125, 4), 25,
rep(125, 8), 25,
rep(125, 7))
)
})
date()
date()
## we can't do this, as over 10^60 genotypes
## test_that("McFL: Relative ordering of number of clones with mut prop growth and init and scrambled names", {
## max.tries <- 4
## for(tries in 1:max.tries) {
## ## Can occasionally blow up with pE.f: pE not finite.
## cat("\n x2gh: a runif is", runif(1), "\n")
## pops <- 50
## ft <- 1
## lni <- 200
## no <- 1e3
## ni <- c(5, 2, rep(0, lni))
## ## scramble around names
## names(ni) <- c("thisistheagene",
## "thisisthebgene",
## replicate(lni,
## paste(sample(letters, 12), collapse = "")))
## ni <- ni[order(names(ni))]
## fe <- allFitnessEffects(noIntGenes = ni)
## fm1 <- allMutatorEffects(noIntGenes = c("thisistheagene" = 8))
## mpg <- oncoSimulPop(pops, fe2fl(fe), muEF = fm1,
## finalTime = ft,
## mutationPropGrowth = TRUE,
## initSize = no, model = "McFL",
## initMutant = "thisistheagene",
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## mnpg <- oncoSimulPop(pops, fe2fl(fe), muEF = fm1,
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no, model = "McFL",
## initMutant = "thisistheagene",
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## pg <- oncoSimulPop(pops, fe2fl(fe),
## finalTime = ft,
## mutationPropGrowth = TRUE,
## initSize = no, model = "McFL",
## initMutant = "thisistheagene",
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## npg <- oncoSimulPop(pops, fe2fl(fe),
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no, model = "McFL",
## initMutant = "thisistheagene",
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## ## ## I once saw a weird thing
## ## expect_true(var(summary(mpg)$NumClones) > 1e-4)
## ## expect_true(var(summary(mnpg)$NumClones) > 1e-4)
## ## expect_true(var(summary(pg)$NumClones) > 1e-4)
## ## expect_true(var(summary(npg)$NumClones) > 1e-4)
## ## These are the real tests
## T1 <- ( wilcox.test(summary(mpg)$NumClones,
## summary(mnpg)$NumClones, alternative = "greater")$p.value < p.value.threshold)
## T2 <- (wilcox.test(summary(mpg)$NumClones,
## summary(pg)$NumClones, alternative = "greater")$p.value < p.value.threshold)
## T3 <- ( wilcox.test(summary(mnpg)$NumClones,
## summary(npg)$NumClones, alternative = "greater")$p.value < p.value.threshold)
## T4 <- ( wilcox.test(summary(pg)$NumClones,
## summary(npg)$NumClones, alternative = "greater")$p.value < p.value.threshold)
## T5 <- (t.test(mutsPerClone(mpg), mutsPerClone(mnpg), alternative = "greater")$p.value < p.value.threshold)
## T6 <- (t.test(mutsPerClone(mpg), mutsPerClone(pg), alternative = "greater")$p.value < p.value.threshold)
## T7 <- (t.test(mutsPerClone(mnpg), mutsPerClone(npg), alternative = "greater")$p.value < p.value.threshold)
## T8 <- (t.test(mutsPerClone(pg), mutsPerClone(npg), alternative = "greater")$p.value < p.value.threshold)
## if( T1 && T3 && T4 && T5 && T6 && T7 && T8 ) break;
## }
## cat(paste("\n done tries", tries, "\n"))
## expect_true(T1 && T3 && T4 && T5 && T6 && T7 && T8)
## })
## date()
##### Comparisons against expected freqs, using a chi-square
## If any mu is very large or any lni is very large, it can fail unless
## pops is very large. And having a large mutator effect is like having a
## very large mu. We want to use very small finalTime: It is birth and
## rate that compound processes and of course we have non-independent
## sampling (overdispersion) which can make chisq a bad idea.
## Thus, we use a tiny final time so we are basically getting just
## mutation events. We need to use a large number of pops to try to avoid
## empty cells with low mutation frequencies.
## We will play with the mutator effects. Note also that here mutator is
## specified passing a vector of same size as genome. It would be faster
## to use just the name of mutator gene.
## date()
## test_that("Expect freq genotypes, mutator and var mut rates", {
## ## We test that mutator does not affect expected frequencies of
## ## mutated genes: they are given by the mutation rate of each gene.
## max.tries <- 4
## for(tries in 1:max.tries) {
## cat("\n u6: a runif is", runif(1), "\n")
## pops <- 20
## ft <- .0001
## lni <- 80
## no <- 5e7
## ni <- c(0, 0, 0, rep(0, lni))
## ## scramble around names
## names(ni) <- c("hereisoneagene",
## "oreoisasabgene",
## "nnhsisthecgene",
## replicate(lni,
## paste(sample(letters, 12), collapse = "")))
## ni <- ni[order(names(ni))]
## fe <- allFitnessEffects(noIntGenes = ni)
## ## of course, passing a mutator of 1 makes everything slow.
## mutator1 <- rep(1, lni + 3)
## ## pg1 <- rep(1e-5, lni + 3)
## pg1 <- runif(lni + 3, min = 1e-7, max = 1e-4) ## max should not be
## ## huge here as mutator
## ## is 34. Can get beyond
## ## 1
## names(mutator1) <- sample(names(ni))
## names(pg1) <- sample(names(ni))
## mutator1["oreoisasabgene"] <- 34 ## 53
## m1 <- allMutatorEffects(noIntGenes = mutator1)
## ## have something with much larger mutation rate
## pg1["hereisoneagene"] <- 1e-3 ## 1e-3
## m1.pg1.b <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = pg1,
## muEF = m1,
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## initMutant ="oreoisasabgene",
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## expect_true(sm("oreoisasabgene", m1.pg1.b) == totalind(m1.pg1.b))
## enom("oreoisasabgene", pg1, no, pops)
## snom("oreoisasabgene", m1.pg1.b)
## p.fail <- 1e-3
## T1 <- (chisq.test(snom("oreoisasabgene", m1.pg1.b),
## p = pnom("oreoisasabgene", pg1, no, pops))$p.value > p.fail)
## if(T1) break;
## }
## cat(paste("\n done tries", tries, "\n"))
## expect_true(T1)
## })
## date()
## date()
## test_that("McFL, Expect freq genotypes, mutator and var mut rates", {
## max.tries <- 4
## for(tries in 1:max.tries) {
## ## We test that mutator does not affect expected frequencies of
## ## mutated genes: they are given by the mutation rate of each gene.
## cat("\n mcfu6: a runif is", runif(1), "\n")
## pops <- 20
## ft <- .0001
## lni <- 80
## no <- 5e7
## ni <- c(0, 0, 0, rep(0, lni))
## ## scramble around names
## names(ni) <- c("hereisoneagene",
## "oreoisasabgene",
## "nnhsisthecgene",
## replicate(lni,
## paste(sample(letters, 12), collapse = "")))
## ni <- ni[order(names(ni))]
## fe <- allFitnessEffects(noIntGenes = ni)
## ## of course, passing a mutator of 1 makes everything slow.
## mutator1 <- rep(1, lni + 3)
## ## pg1 <- rep(1e-5, lni + 3)
## pg1 <- runif(lni + 3, min = 1e-7, max = 1e-4) ## max should not be
## ## huge here as mutator
## ## is 34. Can get beyond
## ## 1
## names(mutator1) <- sample(names(ni))
## names(pg1) <- sample(names(ni))
## mutator1["oreoisasabgene"] <- 47
## m1 <- allMutatorEffects(noIntGenes = mutator1)
## ## have something with much larger mutation rate
## pg1["hereisoneagene"] <- 1e-3 ## 1e-3
## m1.pg1.b <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = pg1,
## muEF = m1,
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## model = "McFL",
## initMutant ="oreoisasabgene",
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## expect_true(sm("oreoisasabgene", m1.pg1.b) == totalind(m1.pg1.b))
## enom("oreoisasabgene", pg1, no, pops)
## snom("oreoisasabgene", m1.pg1.b)
## p.fail <- 1e-3
## T1 <- (chisq.test(snom("oreoisasabgene", m1.pg1.b),
## p = pnom("oreoisasabgene", pg1, no, pops))$p.value > p.fail)
## summary(m1.pg1.b)[, c(1:3, 8:9)]
## if(T1 ) break;
## }
## cat(paste("\n done tries", tries, "\n"))
## expect_true( T1 )
## })
## date()
## date()
## test_that("McFL: Same mu vector, different mutator; diffs in number muts, tiny t", {
## max.tries <- 4
## for(tries in 1:max.tries) {
## ## Here, there is no reproduction or death. Just mutation. And no double
## ## mutants either.
## ## We test:
## ## - mutator increases mutation rates as seen in:
## ## - number of clones created
## ## - number of total mutation events
## cat("\n nm2: a runif is", runif(1), "\n")
## pops <- 20
## ft <- .0001
## lni <- 100
## no <- 1e7
## fi <- rep(0, lni)
## muvector <- rep(5e-6, lni)
## ## scrambling names
## names(fi) <- replicate(lni,
## paste(sample(letters, 12), collapse = ""))
## names(muvector) <- sample(names(fi))
## ## choose something for mutator
## mutator10 <- mutator100 <- fi[5]
## mutator10[] <- 10
## mutator100[] <- 100
## fe <- allFitnessEffects(noIntGenes = fi)
## m10 <- allMutatorEffects(noIntGenes = mutator10)
## m100 <- allMutatorEffects(noIntGenes = mutator100)
## pop10 <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = muvector,
## muEF = m10,
## model = "McFL",
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## initMutant = names(mutator10),
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## pop100 <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = muvector,
## muEF = m100,
## model = "McFL",
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## initMutant = names(mutator10),
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## ## number of total mutations
## T1 <- (smAnomPi(pop10, names(mutator10)) < smAnomPi(pop100, names(mutator100)))
## ## number of clones
## T2 <- (wilcox.test(NClones(pop100), NClones(pop10), alternative = "greater")$p.value < p.value.threshold)
## T3 <- (t.test(mutsPerClone(pop100), mutsPerClone(pop10), alternative = "greater")$p.value < p.value.threshold)
## if( T1 && T2 && T3 ) break;
## }
## cat(paste("\n done tries", tries, "\n"))
## expect_true(T1 && T2 && T3 )
## })
## date()
## date()
## test_that("McFL: Same mu vector, different mutator; diffs in number muts, larger t", {
## ## reproduction, death, and double and possibly triple mutants. We
## ## decrease init pop size to make this fast.
## max.tries <- 4
## for(tries in 1:max.tries) {
## cat("\n nm3: a runif is", runif(1), "\n")
## pops <- 20
## ft <- 1
## lni <- 100
## no <- 1e5
## fi <- rep(0, lni)
## muvector <- rep(5e-6, lni)
## ## scrambling names
## names(fi) <- replicate(lni,
## paste(sample(letters, 12), collapse = ""))
## names(muvector) <- sample(names(fi))
## ## choose something for mutator
## mutator10 <- mutator100 <- fi[5]
## mutator10[] <- 10
## mutator100[] <- 100
## fe <- allFitnessEffects(noIntGenes = fi)
## m10 <- allMutatorEffects(noIntGenes = mutator10)
## m100 <- allMutatorEffects(noIntGenes = mutator100)
## pop10 <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = muvector,
## muEF = m10,
## model = "McFL",
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## initMutant = names(mutator10),
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## pop100 <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = muvector,
## muEF = m100,
## model = "McFL",
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## initMutant = names(mutator10),
## onlyCancer = FALSE, detectionProb = NA, seed = NULL, mc.cores = 2)
## ## number of total mutations
## expect_true(smAnomPi(pop10, names(mutator10)) < smAnomPi(pop100, names(mutator100)))
## ## number of clones
## T1 <- (wilcox.test(NClones(pop100), NClones(pop10), alternative = "greater")$p.value < p.value.threshold)
## T2 <- (t.test(mutsPerClone(pop100), mutsPerClone(pop10), alternative = "greater")$p.value < p.value.threshold)
## if(T1 && T2) break;
## }
## cat(paste("\n done tries", tries, "\n"))
## expect_true(T1 && T2 )
## })
## date()
## date()
## test_that(" MCFL Init with different mutators", {
## max.tries <- 4
## for(tries in 1:max.tries) {
## cat("\n mcfl_z2: a runif is", runif(1), "\n")
## pops <- 20
## ft <- .005
## lni <- 50
## no <- 1e7
## ni <- c(0, 0, 0, rep(0, lni))
## ## scramble around names
## names(ni) <- c("hereisoneagene",
## "oreoisasabgene",
## "nnhsisthecgene",
## replicate(lni,
## paste(sample(letters, 12), collapse = "")))
## ni <- ni[order(names(ni))]
## fe <- allFitnessEffects(noIntGenes = ni)
## mutator1 <- mutator2 <- rep(1, lni + 3)
## pg1 <- rep(5e-6, lni + 3)
## ## scramble names of mutator and per-gene too
## names(mutator1) <- sample(names(ni))
## names(pg1) <- sample(names(ni))
## mutator1["hereisoneagene"] <- 100
## mutator1["oreoisasabgene"] <- 1
## mutator1["nnhsisthecgene"] <- 0.01
## m1 <- allMutatorEffects(noIntGenes = mutator1)
## m1.pg1.a <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = pg1,
## muEF = m1,
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## model = "McFL",
## initMutant = "hereisoneagene",
## sampleEvery = 0.01, keepEvery = 5, seed = NULL,
## onlyCancer = FALSE, detectionProb = NA, mc.cores = 2)
## m1.pg1.b <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = pg1,
## muEF = m1,
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## model = "McFL",
## initMutant = "oreoisasabgene",
## sampleEvery = 0.01, keepEvery = 5, seed = NULL,
## onlyCancer = FALSE, detectionProb = NA, mc.cores = 2)
## m1.pg1.c <- oncoSimulPop(pops,
## fe2fl(fe),
## mu = pg1,
## muEF = m1,
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## model = "McFL",
## initMutant = "nnhsisthecgene",
## sampleEvery = 0.01, keepEvery = 5, seed = NULL,
## onlyCancer = FALSE, detectionProb = NA, mc.cores = 2)
## T1 <- (wilcox.test(NClones(m1.pg1.a), NClones(m1.pg1.b), alternative = "greater")$p.value < p.value.threshold)
## T2 <- (wilcox.test(NClones(m1.pg1.b), NClones(m1.pg1.c), alternative = "greater")$p.value < p.value.threshold)
## T3 <- (t.test(mutsPerClone(m1.pg1.a), mutsPerClone(m1.pg1.b), alternative = "greater")$p.value < p.value.threshold)
## T4 <- (t.test(mutsPerClone(m1.pg1.b), mutsPerClone(m1.pg1.c), alternative = "greater")$p.value < p.value.threshold)
## T5 <- (smAnomPi(m1.pg1.a, "hereisoneagene") >
## smAnomPi(m1.pg1.b, "oreoisasabgene"))
## T6 <- (smAnomPi(m1.pg1.b, "oreoisasabgene") >
## smAnomPi(m1.pg1.c, "nnhsisthecgene"))
## summary(m1.pg1.a)[, c(1:3, 8:9)]
## summary(m1.pg1.b)[, c(1:3, 8:9)]
## summary(m1.pg1.c)[, c(1:3, 8:9)]
## if(T1 && T2 && T3 && T4 && T5 && T6) break;
## }
## cat(paste("\n done tries", tries, "\n"))
## expect_true(T1 && T2 && T3 && T4 && T5 && T6)
## })
## date()
## date()
## test_that("Mutator, several modules differences", {
## max.tries <- 4
## for(tries in 1:max.tries) {
## cat("\n mmdSM1: a runif is", runif(1), "\n")
## reps <- 60
## no <- 5e3
## ft <- 50 ## you need it large enough to get enough hits
## mu <- 1e-5
## ln <- 100
## m1 <- 5 ## if this is too large, easy to get it to blow.
## ## Having three renders it too unpredictable in time.
## ## ni <- rep(0, 3 * ln)
## ## gna <- paste0("a", 1:ln)
## ## gnb <- paste0("b", 1:ln)
## ## gnc <- paste0("c", 1:ln)
## ## names(ni) <- c(gna, gnb, gnc)
## ## gn1 <- paste(c(gna, gnb, gnc), collapse = ", ")
## ## gna <- paste(gna, collapse = ", ")
## ## gnb <- paste(gnb, collapse = ", ")
## ## gnc <- paste(gnc, collapse = ", ")
## ## mut1 <- allMutatorEffects(epistasis = c("A" = m1),
## ## geneToModule = c("A" = gn1))
## ## mut2 <- allMutatorEffects(epistasis = c("A" = m1,
## ## "B" = m1,
## ## "C" = m1),
## ## geneToModule = c("A" = gna,
## ## "B" = gnb,
## ## "C" = gnc))
## ni <- rep(0, 2 * ln)
## gna <- paste0("a", 1:ln)
## gnb <- paste0("b", 1:ln)
## names(ni) <- c(gna, gnb )
## gn1 <- paste(c(gna, gnb), collapse = ", ")
## gna <- paste(gna, collapse = ", ")
## gnb <- paste(gnb, collapse = ", ")
## mut1 <- allMutatorEffects(epistasis = c("A" = m1),
## geneToModule = c("A" = gn1))
## mut2 <- allMutatorEffects(epistasis = c("A" = m1,
## "B" = m1),
## geneToModule = c("A" = gna,
## "B" = gnb))
## f1 <- allFitnessEffects(noIntGenes = ni)
## b1 <- oncoSimulPop(reps,
## fe2fl(f1),
## mu = mu,
## muEF = mut1,
## onlyCancer = FALSE, detectionProb = NA,
## initSize = no,
## finalTime = ft,
## seed = NULL, mc.cores = 2
## )
## gc()
## b2 <- oncoSimulPop(reps,
## fe2fl(f1),
## mu = mu,
## muEF = mut2,
## onlyCancer = FALSE, detectionProb = NA,
## initSize = no,
## finalTime = ft,
## seed = NULL, mc.cores = 2
## )
## gc()
## ## p.value.threshold <- 0.05 ## diffs here small, unless large reps
## ## summary(b2)[, c(1:3, 8:9)]
## ## summary(b1)[, c(1:3, 8:9)]
## ## mean(mutsPerClone(b2))
## ## mean(mutsPerClone(b1))
## ## This is, of course, affected by sampling only at end: we do not see
## ## the many intermediate events.
## T1 <- ( wilcox.test(summary(b2)$NumClones,
## summary(b1)$NumClones, alternative = "greater")$p.value < p.value.threshold)
## T2 <- (t.test(mutsPerClone(b2), mutsPerClone(b1), alternative = "greater")$p.value < p.value.threshold)
## ## it very rarely fails; what are the p-values?
## print(suppressWarnings(wilcox.test(summary(b2)$NumClones,
## summary(b1)$NumClones, alternative = "greater")$p.value))
## print(suppressWarnings(t.test(mutsPerClone(b2), mutsPerClone(b1), alternative = "greater")$p.value))
## if( T1 && T2 ) break;
## }
## cat(paste ("\n done tries ", tries, "\n"))
## expect_true( (T1 && T2) )
## })
## date()
## date()
## test_that("Mutator and mutPropGrowth, mcfl", {
## max.tries <- 4
## for(tries in 1:max.tries) {
## ## we stop on size
## ## Note that names of modules are different. Just for fun.
## cat("\n mmpg_mcfl: a runif is", runif(1), "\n")
## reps <- 14
## no <- 5e3
## ds <- 1.5e4
## ft <- 1000
## mu <- 2e-5
## ln <- 50
## m1 <- 1 ## if this is too large, easy to get it to blow.
## m50 <- 50
## gna <- paste0("a", 1:ln)
## gna <- paste(gna, collapse = ", ")
## f1 <- allFitnessEffects(epistasis = c("A" = 1),
## geneToModule = c("A" = gna))
## mut1 <- allMutatorEffects(epistasis = c("M" = m1),
## geneToModule = c("M" = gna))
## mut50 <- allMutatorEffects(epistasis = c("M" = m50),
## geneToModule = c("M" = gna))
## m1.pg <- oncoSimulPop(reps,
## fe2fl(f1),
## mu = mu,
## muEF = mut1,
## mutationPropGrowth = TRUE,
## onlyCancer = FALSE, detectionProb = NA,
## initSize = no,
## finalTime = ft,
## detectionSize = ds,
## sampleEvery = 0.05,
## keepEvery = 5,
## seed = NULL, mc.cores = 2, model = "McFL"
## )
## ## gc()
## m1.npg <- oncoSimulPop(reps,
## f1,
## mu = mu,
## muEF = mut1,
## mutationPropGrowth = FALSE,
## onlyCancer = FALSE, detectionProb = NA,
## initSize = no,
## finalTime = ft,
## detectionSize = ds,
## sampleEvery = 0.05,
## keepEvery = 5,
## seed = NULL, mc.cores = 2, model = "McFL"
## )
## ##gc()
## m50.pg <- oncoSimulPop(reps,
## fe2fl(f1),
## mu = mu,
## muEF = mut50,
## mutationPropGrowth = TRUE,
## onlyCancer = FALSE, detectionProb = NA,
## initSize = no,
## finalTime = ft,
## detectionSize = ds,
## sampleEvery = 0.05,
## keepEvery = 5,
## seed = NULL, mc.cores = 2, model = "McFL"
## )
## ## gc()
## m50.npg <- oncoSimulPop(reps,
## fe2fl(f1),
## mu = mu,
## muEF = mut50,
## mutationPropGrowth = FALSE,
## onlyCancer = FALSE, detectionProb = NA,
## initSize = no,
## finalTime = ft,
## detectionSize = ds,
## sampleEvery = 0.05,
## keepEvery = 5,
## seed = NULL, mc.cores = 2, model = "McFL"
## )
## ##gc()
## summary(m1.pg)[, c(1:3, 8:9)]
## summary(m50.pg)[, c(1:3, 8:9)]
## summary(m1.npg)[, c(1:3, 8:9)]
## summary(m50.npg)[, c(1:3, 8:9)]
## ## Over mutator, as we have mutPropGrowth, clones, etc, increase
## T1 <- ( wilcox.test(summary(m1.pg)$NumClones,
## summary(m1.npg)$NumClones,
## alternative = "greater")$p.value < p.value.threshold)
## T2 <- (t.test(mutsPerClone(m1.pg),
## mutsPerClone(m1.npg),
## alternative = "greater")$p.value < p.value.threshold)
## T3 <- ( wilcox.test(summary(m50.pg)$NumClones,
## summary(m50.npg)$NumClones,
## alternative = "greater")$p.value < p.value.threshold)
## T4 <- (t.test(mutsPerClone(m50.pg),
## mutsPerClone(m50.npg),
## alternative = "greater")$p.value < p.value.threshold)
## ## Over mutPropGrowth, as we increase mutator, clones, etc, increase
## T5 <- ( wilcox.test(summary(m50.pg)$NumClones,
## summary(m1.pg)$NumClones,
## alternative = "greater")$p.value < p.value.threshold)
## T6 <- (t.test(mutsPerClone(m50.pg),
## mutsPerClone(m1.pg),
## alternative = "greater")$p.value < p.value.threshold)
## T7 <- ( wilcox.test(summary(m50.npg)$NumClones,
## summary(m1.npg)$NumClones,
## alternative = "greater")$p.value < p.value.threshold)
## T8 <- (t.test(mutsPerClone(m50.npg),
## mutsPerClone(m1.npg),
## alternative = "greater")$p.value < p.value.threshold)
## if( T1 && T2 && T3 && T4 && T5 && T6 && T7 && T8) break;
## }
## cat(paste("\n done tries", tries, "\n"))
## expect_true(T1 && T2 && T3 && T4 && T5 && T6 && T7 && T8)
## })
## date()
######################################################################
######################################################################
## ## Some checks of C++ code with mutator and per-gene mutation rates
######################################################################
######################################################################
## set.seed(2)
## ft <- 600
## lni <- 10
## no <- 5e3
## ni <- c(0, 0, 0, rep(0, lni))
## ## scramble around names
## names(ni) <- c("hereisoneagene",
## "oreoisasabgene",
## "nnhsisthecgene",
## replicate(lni,
## paste(sample(letters, 12), collapse = "")))
## ni <- ni[order(names(ni))]
## ni <- sample(ni)
## fe <- allFitnessEffects(noIntGenes = ni)
## mutator1 <- rep(1, lni + 3)
## pg1 <- rep(1e-7, lni + 3) ## what breaks is using 1e-7 instead of 1e-9
## ## scramble names of mutator and per-gene too
## names(mutator1) <- sample(names(ni))
## names(pg1) <- sample(names(ni))
## mutator1["hereisoneagene"] <- 100
## mutator1["oreoisasabgene"] <- 5
## mutator1["nnhsisthecgene"] <- 0.01
## m1 <- allMutatorEffects(noIntGenes = mutator1)
## pg1["hereisoneagene"] <- 1e-5
## pg1["oreoisasabgene"] <- 1e-7
## pg1["nnhsisthecgene"] <- 1e-10
## m1.pg1.a <- oncoSimulIndiv(fe,
## mu = pg1,
## muEF = m1,
## finalTime = ft,
## mutationPropGrowth = FALSE,
## initSize = no,
## initMutant = "hereisoneagene",
## onlyCancer = FALSE,
## verbosity = 6)
## m1.pg1.a
## ## Next is handy to compare
## mu.nop <- function(p, mu) {
## nni <- names(ni)[p]
## pmu <- which(names(mu) %in% nni)
## mu[-pmu]
## }
## ## genes a and b are 1 and 10
## ni[c(1, 10)]
## ##### iteration 928,
## ## parent:
## 100 * sum(mu.nop(c(10), pg1))
## ## child
## 100 * sum(mu.nop(c(10, 3), pg1))
## ### iteration 255
## ## parent
## 100 * 5 * sum(mu.nop(c(10, 1), pg1))
## ## child
## 100 * 5 * sum(mu.nop(c(10, 1, 11), pg1))
## ## ## Some checks on C++ of mutationPropGrowth and mutator
## ni <- rep(0.4, 20)
## names(ni) <- c("a", "b", "c", "d", paste0("n", 1:16))
## fe <- allFitnessEffects(noIntGenes = ni)
## fm6 <- allMutatorEffects(noIntGenes = c("a" = 15,
## "b" = 15,
## "c" = 15,
## "d" = 15))
## set.seed(5) ## if you use seed of 2, it blows up with huge birth
## ## rate and many mutations > 1
## ## But then, 1.4^20 is > 800!
## oncoSimulIndiv(fe, muEF = fm6, finalTime =30,
## mutationPropGrowth = TRUE,
## initSize = 1e4,
## mu = 1e-06,
## verbosity = 6,
## onlyCancer = FALSE, seed = NULL)
## ## ###### Iteration 40.
## ## ## mutation
## ## ## parent
## ## 20 * 1e-06
## ## ## child
## ## 1.4 * 15 * 1e-06 * 19
## ## ###### Iteration 39
## ## ## mutation
## ## ## parent
## ## 1.4 * 15 * 1e-06 * 19
## ## ## chlid
## ## 1.4 * 1.4 * 15 * 1e-06 * 18
## ## If you want to check the internal running of C++, use verbosity >=
## ## 3. For instance, at iteration 3 a new species is created from the
## ## wildtype. The new clone has mutated gene 1. Thus the new mutation rate
## ## should be
## 10 * 3 * 1e-6 ## 10 for the effect of a, 3 for the number of remaining
## ## genes, and 1e-6 for the baseline mutation rate
## That is 3e-5, as shown.
## At iteration 9 we create a mutant at 3 from wildtype and its mutation rate is
## 30 * 3 * 1e-6
## At iteration 254 we create species 1,2,4 from 1,2. The parent mutation is
## 10 * 20 * 2 * 1e-6 = 4e-4
## and the child's is
## 10 * 20 * 40 * 1 * 1e-6 = 0.008
## set.seed(1)
## fe <- allFitnessEffects(noIntGenes = c("a" = 0.12,
## "b" = 0.14,
## "c" = 0.16,
## "d" = 0.11))
## fm6 <- allMutatorEffects(noIntGenes = c("a" = 10,
## "b" = 20,
## "c" = 30,
## "d" = 40))
## oncoSimulIndiv(fe, muEF = fm6, finalTime =100,
## mutationPropGrowth = FALSE,
## initSize = 1e5,
## verbosity = 6,
## onlyCancer = FALSE,
## seed = NULL)
cat(paste("\n Finished test.flfast-mutator.R test at", date(), "\n"))
cat(paste(" Took ", round(difftime(Sys.time(), inittime, units = "secs"), 2), "\n\n"))
rm(inittime)
