## ----echo=FALSE,results='hide',error=FALSE------------------------------- require(knitr, quietly = TRUE) opts_knit$set(concordance = TRUE) options(width = 68) ##opts_knit$set(stop_on_error = 2L) ## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"------------------ BiocStyle::latex() ## ----results="hide"------------------------------------------------------ library(OncoSimulR) library(graph) library(igraph) igraph_options(vertex.color = "SkyBlue2") ## ----echo=FALSE, results='hide'------------------------------------------ options(width = 68) ## ------------------------------------------------------------------------ packageVersion("OncoSimulR") ## ----fig.width=6.5, fig.height=10---------------------------------------- ## 1. Fitness effects: here we specify a ## epistatic model with modules. sa <- 0.1 sb <- -0.2 sab <- 0.25 sac <- -0.1 sbc <- 0.25 sv2 <- allFitnessEffects(epistasis = c("-A : B" = sb, "A : -B" = sa, "A : C" = sac, "A:B" = sab, "-A:B:C" = sbc), geneToModule = c( "A" = "a1, a2", "B" = "b", "C" = "c")) evalAllGenotypes(sv2, order = FALSE, addwt = TRUE) ## 2. Simulate the data. Here we use the "McFL" model and set explicitly ## parameters for mutation rate, final and initial sizes, etc. RNGkind("Mersenne-Twister") set.seed(983) ep1 <- oncoSimulIndiv(sv2, model = "McFL", mu = 5e-6, sampleEvery = 0.02, keepEvery = 0.5, initSize = 2000, finalTime = 3000, onlyCancer = FALSE) ## ----iep1x1,fig.width=6.5, fig.height=9.5-------------------------------- ## 3. We will not analyze those data any further. We will only plot them. ## For the sake of a small plot, we thin the data. par(mfrow = c(2, 1)) plot(ep1, show = "drivers", xlim = c(0, 1500), thinData = TRUE, thinData.keep = 0.5) ## Increase ylim and legend.ncols to avoid overlap of ## legend with rest of figure plot(ep1, show = "genotypes", ylim = c(0, 4500), legend.ncols = 4, xlim = c(0, 1500), thinData = TRUE, thinData.keep = 0.5) ## ------------------------------------------------------------------------ ## 1. Fitness effects: pancr <- allFitnessEffects( data.frame(parent = c("Root", rep("KRAS", 4), "SMAD4", "CDNK2A", "TP53", "TP53", "MLL3"), child = c("KRAS","SMAD4", "CDNK2A", "TP53", "MLL3", rep("PXDN", 3), rep("TGFBR2", 2)), s = 0.1, sh = -0.9, typeDep = "MN")) ## How does it look like? plot(pancr) ## ----fig.width=6.5, fig.height=10---------------------------------------- ## 2. Simulate from it. set.seed(1) ## Fix the seed, so we can repeat it ep2 <- oncoSimulIndiv(pancr, model = "McFL", mu = 1e-6, sampleEvery = 0.02, keepEvery = 1, initSize = 1000, finalTime = 10000, onlyCancer = FALSE) ## ----iep2x2,fig.width=6.5, fig.height=9---------------------------------- ## 3. What genotypes and drivers we get? And play with limits ## to show only parts of the data. We also thin them. par(mfrow = c(2, 1)) par(cex = 0.7) plot(ep2, show = "genotypes", xlim = c(2000, 4000), ylim = c(0, 2400), thinData = TRUE, thinData.keep = 0.5) plot(ep2, show = "drivers", addtot = TRUE, thinData = TRUE, thinData.keep = 0.5) ## ------------------------------------------------------------------------ ai1 <- evalAllGenotypes(allFitnessEffects( noIntGenes = c(0.05, -.2, .1)), order = FALSE) ## ------------------------------------------------------------------------ ai1 ## ------------------------------------------------------------------------ all(ai1[, "Fitness"] == c( (1 + .05), (1 - .2), (1 + .1), (1 + .05) * (1 - .2), (1 + .05) * (1 + .1), (1 - .2) * (1 + .1), (1 + .05) * (1 - .2) * (1 + .1))) ## ------------------------------------------------------------------------ (ai2 <- evalAllGenotypes(allFitnessEffects( noIntGenes = c(0.05, -.2, .1)), order = TRUE, addwt = TRUE)) ## ----fig.height=4-------------------------------------------------------- data(examplesFitnessEffects) plot(examplesFitnessEffects[["cbn1"]]) ## ------------------------------------------------------------------------ cs <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), s = 0.1, sh = -0.9, typeDep = "MN") cbn1 <- allFitnessEffects(cs) ## ----fig.height=3-------------------------------------------------------- plot(cbn1) ## ----fig.height=5-------------------------------------------------------- plot(cbn1, "igraph") ## ----fig.height=5-------------------------------------------------------- library(igraph) ## to make the reingold.tilford layout available plot(cbn1, "igraph", layout = layout.reingold.tilford) ## ------------------------------------------------------------------------ gfs <- evalAllGenotypes(cbn1, order = FALSE) gfs[1:15, ] ## ------------------------------------------------------------------------ c1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)), sh = c(rep(0, 4), c(-.1, -.2), c(-.05, -.06, -.07)), typeDep = "MN") try(fc1 <- allFitnessEffects(c1)) ## ------------------------------------------------------------------------ c1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)), sh = c(rep(0, 4), c(-.9, -.9), rep(-.95, 3)), typeDep = "MN") cbn2 <- allFitnessEffects(c1) ## ------------------------------------------------------------------------ gcbn2 <- evalAllGenotypes(cbn2, order = FALSE) gcbn2[1:10, ] ## ------------------------------------------------------------------------ gcbn2o <- evalAllGenotypes(cbn2, order = TRUE, max = 1956) gcbn2o[1:10, ] ## ------------------------------------------------------------------------ all.equal( gcbn2[c(1:21, 22, 28, 41, 44, 56, 63 ) , "Fitness"], c(1.01, 1.02, 0.1, 1.03, 1.04, 0.05, 1.01 * c(1.02, 0.1, 1.03, 1.04, 0.05), 1.02 * c(0.10, 1.03, 1.04, 0.05), 0.1 * c(1.03, 1.04, 0.05), 1.03 * c(1.04, 0.05), 1.04 * 0.05, 1.01 * 1.02 * 1.1, 1.01 * 0.1 * 0.05, 1.03 * 1.04 * 0.05, 1.01 * 1.02 * 1.1 * 0.05, 1.03 * 1.04 * 1.2 * 0.1, ## notice this 1.01 * 1.02 * 1.03 * 1.04 * 1.1 * 1.2 )) ## ------------------------------------------------------------------------ gcbn2[56, ] ## this is d, e, g, c all.equal(gcbn2[56, "Fitness"], 1.03 * 1.04 * 1.2 * 0.10) ## ------------------------------------------------------------------------ s1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)), sh = c(rep(0, 4), c(-.9, -.9), rep(-.95, 3)), typeDep = "SM") smn1 <- allFitnessEffects(s1) ## ----fig.height=3-------------------------------------------------------- plot(smn1) ## ------------------------------------------------------------------------ gsmn1 <- evalAllGenotypes(smn1, order = FALSE) ## ------------------------------------------------------------------------ gcbn2[c(8, 12, 22), ] gsmn1[c(8, 12, 22), ] gcbn2[c(20:21, 28), ] gsmn1[c(20:21, 28), ] ## ------------------------------------------------------------------------ x1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)), sh = c(rep(0, 4), c(-.9, -.9), rep(-.95, 3)), typeDep = "XMPN") xor1 <- allFitnessEffects(x1) ## ----fig.height=3-------------------------------------------------------- plot(xor1) ## ------------------------------------------------------------------------ gxor1 <- evalAllGenotypes(xor1, order = FALSE) ## ------------------------------------------------------------------------ gxor1[c(22, 41), ] c(1.01 * 1.02 * 0.1, 1.03 * 1.04 * 0.05) ## ------------------------------------------------------------------------ gxor1[28, ] 1.01 * 1.1 * 1.2 ## ------------------------------------------------------------------------ gxor1[44, ] 1.01 * 1.02 * 0.1 * 1.2 ## ------------------------------------------------------------------------ p3 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c", "f"), child = c("a", "b", "d", "e", "c", "c", "f", "f", "g", "g"), s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3), sh = c(rep(0, 4), c(-.9, -.9), c(-.95, -.95), c(-.99, -.99)), typeDep = c(rep("--", 4), "XMPN", "XMPN", "MN", "MN", "SM", "SM")) fp3 <- allFitnessEffects(p3) ## ----fig.height=3-------------------------------------------------------- plot(fp3) ## ----fig.height=6-------------------------------------------------------- plot(fp3, "igraph", layout.reingold.tilford) ## ------------------------------------------------------------------------ gfp3 <- evalAllGenotypes(fp3, order = FALSE) ## ------------------------------------------------------------------------ gfp3[c(9, 24, 29, 59, 60, 66, 119, 120, 126, 127), ] c(1.01 * 1.1, 1.03 * .05, 1.01 * 1.02 * 0.1, 0.1 * 0.05 * 1.3, 1.03 * 1.04 * 1.2, 1.01 * 1.02 * 0.1 * 0.05, 0.1 * 1.03 * 1.04 * 1.2 * 1.3, 1.01 * 1.02 * 0.1 * 1.03 * 1.04 * 1.2, 1.02 * 1.1 * 1.03 * 1.04 * 1.2 * 1.3, 1.01 * 1.02 * 1.03 * 1.04 * 0.1 * 1.2 * 1.3) ## ------------------------------------------------------------------------ s <- 0.2 sboth <- (1/(1 + s)) - 1 m0 <- allFitnessEffects(data.frame( parent = c("Root", "Root", "a1", "a2"), child = c("a1", "a2", "b", "b"), s = s, sh = -1, typeDep = "OR"), epistasis = c("a1:a2" = sboth)) evalAllGenotypes(m0, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ s <- 0.2 m1 <- allFitnessEffects(data.frame( parent = c("Root", "A"), child = c("A", "B"), s = s, sh = -1, typeDep = "OR"), geneToModule = c("Root" = "Root", "A" = "a1, a2", "B" = "b1")) evalAllGenotypes(m1, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ fnme <- allFitnessEffects(epistasis = c("A" = 0.1, "B" = 0.2), geneToModule = c("A" = "a1, a2", "B" = "b1")) evalAllGenotypes(fnme, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ fnme2 <- allFitnessEffects(epistasis = c("A" = 0.1, "B" = 0.2), geneToModule = c( "Root" = "Root", "A" = "a1, a2", "B" = "b1")) evalAllGenotypes(fnme, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ p4 <- data.frame(parent = c(rep("Root", 4), "A", "B", "D", "E", "C", "F"), child = c("A", "B", "D", "E", "C", "C", "F", "F", "G", "G"), s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3), sh = c(rep(0, 4), c(-.9, -.9), c(-.95, -.95), c(-.99, -.99)), typeDep = c(rep("--", 4), "XMPN", "XMPN", "MN", "MN", "SM", "SM")) fp4m <- allFitnessEffects(p4, geneToModule = c("Root" = "Root", "A" = "a1", "B" = "b1, b2", "C" = "c1", "D" = "d1, d2", "E" = "e1", "F" = "f1, f2", "G" = "g1")) ## ----fig.height=3-------------------------------------------------------- plot(fp4m) ## ----fig.height=3-------------------------------------------------------- plot(fp4m, expandModules = TRUE) ## ----fig.height=8-------------------------------------------------------- plot(fp4m, "igraph", layout = layout.reingold.tilford, expandModules = TRUE) ## ------------------------------------------------------------------------ gfp4 <- evalAllGenotypes(fp4m, order = FALSE, max = 1024) ## ------------------------------------------------------------------------ gfp4[c(12, 20, 21, 40, 41, 46, 50, 55, 64, 92, 155, 157, 163, 372, 632, 828), ] c(1.01 * 1.02, 1.02, 1.02 * 1.1, 0.1 * 1.3, 1.03, 1.03 * 1.04, 1.04 * 0.05, 0.05 * 1.3, 1.01 * 1.02 * 0.1, 1.02 * 1.1, 0.1 * 0.05 * 1.3, 1.03 * 0.05, 1.03 * 0.05, 1.03 * 1.04 * 1.2, 1.03 * 1.04 * 1.2, 1.02 * 1.1 * 1.03 * 1.04 * 1.2 * 1.3) ## ------------------------------------------------------------------------ o3 <- allFitnessEffects(orderEffects = c( "F > D > M" = -0.3, "D > F > M" = 0.4, "D > M > F" = 0.2, "D > M" = 0.1, "M > D" = 0.5), geneToModule = c("M" = "m", "F" = "f", "D" = "d") ) (ag <- evalAllGenotypes(o3, addwt = TRUE)) ## ------------------------------------------------------------------------ ofe1 <- allFitnessEffects(orderEffects = c("F > D" = -0.3, "D > F" = 0.4), geneToModule = c("F" = "f1, f2", "D" = "d1, d2") ) ag <- evalAllGenotypes(ofe1) ## ------------------------------------------------------------------------ ag[5:16,] ## ------------------------------------------------------------------------ ag[c(17, 39, 19, 29), ] ## ------------------------------------------------------------------------ ag[c(17, 39, 19, 29), "Fitness"] == c(1.4, 0.7, 1.4, 0.7) ## ------------------------------------------------------------------------ ag[c(43, 44),] ag[c(43, 44), "Fitness"] == c(1.4, 1.4) ## ------------------------------------------------------------------------ all(ag[41:52, "Fitness"] == 1.4) ## ------------------------------------------------------------------------ all(ag[53:64, "Fitness"] == 0.7) ## ------------------------------------------------------------------------ ofe2 <- allFitnessEffects(orderEffects = c("F > D" = -0.3, "D > F" = 0.4), geneToModule = c("F" = "f1, f2, f3", "D" = "d1, d2") ) ag2 <- evalAllGenotypes(ofe2, max = 325) ## ------------------------------------------------------------------------ all(ag2[grep("^d.*f.*", ag2[, 1]), "Fitness"] == 1.4) all(ag2[grep("^f.*d.*", ag2[, 1]), "Fitness"] == 0.7) oe <- c(grep("^f.*d.*", ag2[, 1]), grep("^d.*f.*", ag2[, 1])) all(ag2[-oe, "Fitness"] == 1) ## ------------------------------------------------------------------------ foi1 <- allFitnessEffects( orderEffects = c("D>B" = -0.2, "B > D" = 0.3), noIntGenes = c("A" = 0.05, "C" = -.2, "E" = .1)) ## ------------------------------------------------------------------------ foi1[c("geneModule", "long.geneNoInt")] ## ------------------------------------------------------------------------ agoi1 <- evalAllGenotypes(foi1, max = 325) head(agoi1) ## ------------------------------------------------------------------------ rn <- 1:nrow(agoi1) names(rn) <- agoi1[, 1] agoi1[rn[LETTERS[1:5]], "Fitness"] == c(1.05, 1, 0.8, 1, 1.1) ## ------------------------------------------------------------------------ agoi1[grep("^A > [BD]$", names(rn)), "Fitness"] == 1.05 agoi1[grep("^C > [BD]$", names(rn)), "Fitness"] == 0.8 agoi1[grep("^E > [BD]$", names(rn)), "Fitness"] == 1.1 agoi1[grep("^[BD] > A$", names(rn)), "Fitness"] == 1.05 agoi1[grep("^[BD] > C$", names(rn)), "Fitness"] == 0.8 agoi1[grep("^[BD] > E$", names(rn)), "Fitness"] == 1.1 ## ------------------------------------------------------------------------ all.equal(agoi1[230:253, "Fitness"] , rep((1 - 0.2) * 1.05 * 0.8 * 1.1, 24)) ## ------------------------------------------------------------------------ s <- 0.2 sv <- allFitnessEffects(epistasis = c("-A : B" = -1, "A : -B" = -1, "A:B" = s)) ## ------------------------------------------------------------------------ (asv <- evalAllGenotypes(sv, order = FALSE, addwt = TRUE)) ## ------------------------------------------------------------------------ evalAllGenotypes(sv, order = TRUE, addwt = TRUE) ## ------------------------------------------------------------------------ evalAllGenotypes(sv, order = FALSE, addwt = TRUE, model = "Bozic") ## ------------------------------------------------------------------------ s <- 0.2 svB <- allFitnessEffects(epistasis = c("-A : B" = -Inf, "A : -B" = -Inf, "A:B" = s)) evalAllGenotypes(svB, order = FALSE, addwt = TRUE, model = "Bozic") ## ------------------------------------------------------------------------ s <- 1 svB1 <- allFitnessEffects(epistasis = c("-A : B" = -Inf, "A : -B" = -Inf, "A:B" = s)) evalAllGenotypes(svB1, order = FALSE, addwt = TRUE, model = "Bozic") s <- 3 svB3 <- allFitnessEffects(epistasis = c("-A : B" = -Inf, "A : -B" = -Inf, "A:B" = s)) evalAllGenotypes(svB3, order = FALSE, addwt = TRUE, model = "Bozic") ## ------------------------------------------------------------------------ sa <- -0.1 sb <- -0.2 sab <- 0.25 sv2 <- allFitnessEffects(epistasis = c("-A : B" = sb, "A : -B" = sa, "A:B" = sab), geneToModule = c( "A" = "a1, a2", "B" = "b")) evalAllGenotypes(sv2, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ evalAllGenotypes(sv2, order = TRUE, addwt = TRUE) ## ------------------------------------------------------------------------ sa <- 0.1 sb <- 0.2 sab <- -0.8 sm1 <- allFitnessEffects(epistasis = c("-A : B" = sb, "A : -B" = sa, "A:B" = sab)) evalAllGenotypes(sm1, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ evalAllGenotypes(sm1, order = TRUE, addwt = TRUE) ## ------------------------------------------------------------------------ sa <- 0.2 sb <- 0.3 sab <- 0.7 e2 <- allFitnessEffects(epistasis = c("A: -B" = sa, "-A:B" = sb, "A : B" = sab)) evalAllGenotypes(e2, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ s2 <- ((1 + sab)/((1 + sa) * (1 + sb))) - 1 e3 <- allFitnessEffects(epistasis = c("A" = sa, "B" = sb, "A : B" = s2)) evalAllGenotypes(e3, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ sa <- 0.1 sb <- 0.15 sc <- 0.2 sab <- 0.3 sbc <- -0.25 sabc <- 0.4 sac <- (1 + sa) * (1 + sc) - 1 E3A <- allFitnessEffects(epistasis = c("A:-B:-C" = sa, "-A:B:-C" = sb, "-A:-B:C" = sc, "A:B:-C" = sab, "-A:B:C" = sbc, "A:-B:C" = sac, "A : B : C" = sabc) ) evalAllGenotypes(E3A, order = FALSE, addwt = FALSE) ## ------------------------------------------------------------------------ sa <- 0.1 sb <- 0.15 sc <- 0.2 sab <- 0.3 Sab <- ( (1 + sab)/((1 + sa) * (1 + sb))) - 1 Sbc <- ( (1 + sbc)/((1 + sb) * (1 + sc))) - 1 Sabc <- ( (1 + sabc)/( (1 + sa) * (1 + sb) * (1 + sc) * (1 + Sab) * (1 + Sbc) ) ) - 1 E3B <- allFitnessEffects(epistasis = c("A" = sa, "B" = sb, "C" = sc, "A:B" = Sab, "B:C" = Sbc, ## "A:C" = sac, ## not needed now "A : B : C" = Sabc) ) evalAllGenotypes(E3B, order = FALSE, addwt = FALSE) ## ------------------------------------------------------------------------ all(evalAllGenotypes(E3A, order = FALSE, addwt = FALSE) == evalAllGenotypes(E3B, order = FALSE, addwt = FALSE)) ## ------------------------------------------------------------------------ sa <- 0.2 sb <- 0.3 sab <- 0.7 em <- allFitnessEffects(epistasis = c("A: -B" = sa, "-A:B" = sb, "A : B" = sab), geneToModule = c("A" = "a1, a2", "B" = "b1, b2")) evalAllGenotypes(em, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ s2 <- ((1 + sab)/((1 + sa) * (1 + sb))) - 1 em2 <- allFitnessEffects(epistasis = c("A" = sa, "B" = sb, "A : B" = s2), geneToModule = c("A" = "a1, a2", "B" = "b1, b2") ) evalAllGenotypes(em2, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ fnme <- allFitnessEffects(epistasis = c("A" = 0.1, "B" = 0.2), geneToModule = c("A" = "a1, a2", "B" = "b1, b2, b3")) evalAllGenotypes(fnme, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ fnme <- allFitnessEffects(epistasis = c("A" = 0.1, "B" = 0.2, "A : B" = 0.0), geneToModule = c("A" = "a1, a2", "B" = "b1, b2, b3")) evalAllGenotypes(fnme, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ p4 <- data.frame(parent = c(rep("Root", 4), "A", "B", "D", "E", "C", "F"), child = c("A", "B", "D", "E", "C", "C", "F", "F", "G", "G"), s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3), sh = c(rep(0, 4), c(-.9, -.9), c(-.95, -.95), c(-.99, -.99)), typeDep = c(rep("--", 4), "XMPN", "XMPN", "MN", "MN", "SM", "SM")) oe <- c("C > F" = -0.1, "H > I" = 0.12) sm <- c("I:J" = -1) sv <- c("-K:M" = -.5, "K:-M" = -.5) epist <- c(sm, sv) modules <- c("Root" = "Root", "A" = "a1", "B" = "b1, b2", "C" = "c1", "D" = "d1, d2", "E" = "e1", "F" = "f1, f2", "G" = "g1", "H" = "h1, h2", "I" = "i1", "J" = "j1, j2", "K" = "k1, k2", "M" = "m1") set.seed(1) ## for repeatability noint <- rexp(5, 10) names(noint) <- paste0("n", 1:5) fea <- allFitnessEffects(rT = p4, epistasis = epist, orderEffects = oe, noIntGenes = noint, geneToModule = modules) ## ----fig.height=6.5------------------------------------------------------ plot(fea) ## ----fig.height=6.5------------------------------------------------------ plot(fea, "igraph") ## ----fig.height=6.5------------------------------------------------------ plot(fea, expandModules = TRUE) ## ----fig.height=7.------------------------------------------------------- plot(fea, "igraph", expandModules = TRUE) ## ------------------------------------------------------------------------ evalGenotype("k1 > i1 > h2", fea) ## 0.5 evalGenotype("k1 > h1 > i1", fea) ## 0.5 * 1.12 evalGenotype("k2 > m1 > h1 > i1", fea) ## 1.12 evalGenotype("k2 > m1 > h1 > i1 > c1 > n3 > f2", fea) ## 1.12 * 0.1 * (1 + noint[3]) * 0.05 * 0.9 ## ------------------------------------------------------------------------ randomGenotype <- function(fe, ns = NULL) { gn <- setdiff(c(fe$geneModule$Gene, fe$long.geneNoInt$Gene), "Root") if(is.null(ns)) ns <- sample(length(gn), 1) return(paste(sample(gn, ns), collapse = " > ")) } set.seed(2) ## for reproducibility evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: k2 > i1 > c1 > n1 > m1 ## Individual s terms are : 0.0755182 -0.9 ## Fitness: 0.107552 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: n2 > h1 > h2 ## Individual s terms are : 0.118164 ## Fitness: 1.11816 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: d2 > k2 > c1 > f2 > n4 > m1 > n3 > f1 > b1 > g1 > n5 > h1 > j2 ## Individual s terms are : 0.0145707 0.0139795 0.0436069 0.02 0.1 0.03 -0.95 0.3 -0.1 ## Fitness: 0.0725829 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: h2 > c1 > f1 > n2 > b2 > a1 > n1 > i1 ## Individual s terms are : 0.0755182 0.118164 0.01 0.02 -0.9 -0.95 -0.1 0.12 ## Fitness: 0.00624418 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: h2 > j1 > m1 > d2 > i1 > b2 > k2 > d1 > b1 > n3 > n1 > g1 > h1 > c1 > k1 > e1 > a1 > f1 > n5 > f2 ## Individual s terms are : 0.0755182 0.0145707 0.0436069 0.01 0.02 -0.9 0.03 0.04 0.2 0.3 -1 -0.1 0.12 ## Fitness: 0 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: n1 > m1 > n3 > i1 > j1 > n5 > k1 ## Individual s terms are : 0.0755182 0.0145707 0.0436069 -1 ## Fitness: 0 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: d2 > n1 > g1 > f1 > f2 > c1 > b1 > d1 > k1 > a1 > b2 > i1 > n4 > h2 > n2 ## Individual s terms are : 0.0755182 0.118164 0.0139795 0.01 0.02 -0.9 0.03 -0.95 0.3 -0.5 ## Fitness: 0.00420528 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: j1 > f1 > j2 > a1 > n4 > c1 > n3 > k1 > d1 > h1 ## Individual s terms are : 0.0145707 0.0139795 0.01 0.1 0.03 -0.95 -0.5 ## Fitness: 0.0294308 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: n5 > f2 > f1 > h2 > n4 > c1 > n3 > b1 ## Individual s terms are : 0.0145707 0.0139795 0.0436069 0.02 0.1 -0.95 ## Fitness: 0.0602298 evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE) ## Genotype: h1 > d1 > f2 ## Individual s terms are : 0.03 -0.95 ## Fitness: 0.0515 ## ------------------------------------------------------------------------ K <- 5 sd <- 0.1 sdp <- 0.15 sp <- 0.05 bauer <- data.frame(parent = c("Root", rep("p", K)), child = c("p", paste0("s", 1:K)), s = c(sd, rep(sdp, K)), sh = c(0, rep(sp, K)), typeDep = "MN") fbauer <- allFitnessEffects(bauer) ## ----fig.height=3-------------------------------------------------------- plot(fbauer) ## ------------------------------------------------------------------------ (b1 <- evalAllGenotypes(fbauer, order = FALSE))[1:10, ] ## ------------------------------------------------------------------------ (b2 <- evalAllGenotypes(fbauer, order = TRUE, max = 2000))[1:15, ] ## ------------------------------------------------------------------------ length(table(b1$Fitness)) length(table(b2$Fitness)) ## ----echo=FALSE, fig.height=4, fig.width=4------------------------------- df1 <- data.frame(x = c(1, 1.2, 1.4), f = c(1, 1.2, 1.2), names = c("wt", "A", "B")) plot(df1[, 2] ~ df1[, 1], axes = TRUE, xlab= "", ylab = "Fitness", xaxt = "n", yaxt = "n", ylim = c(1, 1.21)) segments(1, 1, 1.2, 1.2) segments(1, 1, 1.4, 1.2) text(1, 1, "wt", pos = 4) text(1.2, 1.2, "A", pos = 2) text(1.4, 1.2, "B", pos = 2) ## axis(1, tick = FALSE, labels = FALSE) ## axis(2, tick = FALSE, labels = FALSE) ## ------------------------------------------------------------------------ s <- 0.1 ## or whatever number m1a1 <- allFitnessEffects(data.frame(parent = c("Root", "Root"), child = c("A", "B"), s = s, sh = 0, typeDep = "MN")) evalAllGenotypes(m1a1, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ s <- 0.1 sab <- 0.3 m1a2 <- allFitnessEffects(epistasis = c("A:-B" = s, "-A:B" = s, "A:B" = sab)) evalAllGenotypes(m1a2, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ sab3 <- ((1 + sab)/((1 + s)^2)) - 1 m1a3 <- allFitnessEffects(data.frame(parent = c("Root", "Root"), child = c("A", "B"), s = s, sh = 0, typeDep = "MN"), epistasis = c("A:B" = sab3)) evalAllGenotypes(m1a3, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ all.equal(evalAllGenotypes(m1a2, order = FALSE, addwt = TRUE), evalAllGenotypes(m1a3, order = FALSE, addwt = TRUE)) ## ----echo=FALSE, fig.width=4, fig.height=4------------------------------- df1 <- data.frame(x = c(1, 1.2, 1.2, 1.4), f = c(1, 0.4, 0.3, 1.3), names = c("wt", "A", "B", "AB")) plot(df1[, 2] ~ df1[, 1], axes = TRUE, xlab= "", ylab = "Fitness", xaxt = "n", yaxt = "n", ylim = c(0.29, 1.32)) segments(1, 1, 1.2, 0.4) segments(1, 1, 1.2, 0.3) segments(1.2, 0.4, 1.4, 1.3) segments(1.2, 0.3, 1.4, 1.3) text(x = df1[, 1], y = df1[, 2], labels = df1[, "names"], pos = c(4, 2, 2, 2)) ## text(1, 1, "wt", pos = 4) ## text(1.2, 1.2, "A", pos = 2) ## text(1.4, 1.2, "B", pos = 2) ## ------------------------------------------------------------------------ sa <- -0.6 sb <- -0.7 sab <- 0.3 m1b1 <- allFitnessEffects(epistasis = c("A:-B" = sa, "-A:B" = sb, "A:B" = sab)) evalAllGenotypes(m1b1, order = FALSE, addwt = TRUE) ## ----echo=FALSE, fig.width=4, fig.height=4------------------------------- df1 <- data.frame(x = c(1, 1.2, 1.2, 1.4), f = c(1, 1.2, 0.7, 1.5), names = c("wt", "A", "B", "AB")) plot(df1[, 2] ~ df1[, 1], axes = TRUE, xlab = "", ylab = "Fitness", xaxt = "n", yaxt = "n", ylim = c(0.69, 1.53)) segments(1, 1, 1.2, 1.2) segments(1, 1, 1.2, 0.7) segments(1.2, 1.2, 1.4, 1.5) segments(1.2, 0.7, 1.4, 1.5) text(x = df1[, 1], y = df1[, 2], labels = df1[, "names"], pos = c(3, 3, 3, 2)) ## text(1, 1, "wt", pos = 4) ## text(1.2, 1.2, "A", pos = 2) ## text(1.4, 1.2, "B", pos = 2) ## ------------------------------------------------------------------------ sa <- 0.2 sb <- -0.3 sab <- 0.5 m1c1 <- allFitnessEffects(epistasis = c("A:-B" = sa, "-A:B" = sb, "A:B" = sab)) evalAllGenotypes(m1c1, order = FALSE, addwt = TRUE) ## ----echo=FALSE, fig.width=4.5, fig.height=3.5--------------------------- df1 <- data.frame(x = c(1, 2, 3, 4), f = c(1.1, 1, 0.95, 1.2), names = c("u", "wt", "i", "v")) plot(df1[, 2] ~ df1[, 1], axes = FALSE, xlab = "", ylab = "") par(las = 1) axis(2) axis(1, at = c(1, 2, 3, 4), labels = df1[, "names"], ylab = "") box() arrows(c(2, 2, 3), c(1, 1, 0.95), c(1, 3, 4), c(1.1, 0.95, 1.2)) ## text(1, 1, "wt", pos = 4) ## text(1.2, 1.2, "A", pos = 2) ## text(1.4, 1.2, "B", pos = 2) ## ------------------------------------------------------------------------ su <- 0.1 si <- -0.05 fvi <- 1.2 ## the fitnes of the vi mutant sv <- (fvi/(1 + si)) - 1 sui <- suv <- -1 od <- allFitnessEffects( data.frame(parent = c("Root", "Root", "i"), child = c("u", "i", "v"), s = c(su, si, sv), sh = -1, typeDep = "MN"), epistasis = c( "u:i" = sui, "u:v" = suv)) ## ----fig.width=3, fig.height=3------------------------------------------- plot(od) ## ------------------------------------------------------------------------ evalAllGenotypes(od, order = FALSE, addwt = TRUE) ## ----echo=FALSE, fig.width=4, fig.height=3------------------------------- df1 <- data.frame(x = c(1, 2, 3), f = c(1, 0.95, 1.2), names = c("wt", "1", "2")) plot(df1[, 2] ~ df1[, 1], axes = FALSE, xlab = "", ylab = "") par(las = 1) axis(2) axis(1, at = c(1, 2, 3), labels = df1[, "names"], ylab = "") box() segments(c(1, 2), c(1, 0.95), c(2, 3), c(0.95, 1.2)) ## text(1, 1, "wt", pos = 4) ## text(1.2, 1.2, "A", pos = 2) ## text(1.4, 1.2, "B", pos = 2) ## ----echo=FALSE, fig.width=4, fig.height=3------------------------------- df1 <- data.frame(x = c(1, 2, 3, 4), f = c(1, 0.95, 0.92, 1.2), names = c("wt", "1", "2", "3")) plot(df1[, 2] ~ df1[, 1], axes = FALSE, xlab = "", ylab = "") par(las = 1) axis(2) axis(1, at = c(1, 2, 3, 4), labels = df1[, "names"], ylab = "") box() segments(c(1, 2, 3), c(1, 0.95, 0.92), c(2, 3, 4), c(0.95, 0.92, 1.2)) ## text(1, 1, "wt", pos = 4) ## text(1.2, 1.2, "A", pos = 2) ## text(1.4, 1.2, "B", pos = 2) ## ------------------------------------------------------------------------ d1 <- -0.05 ## single mutant fitness 0.95 d2 <- -0.08 ## double mutant fitness 0.92 d3 <- 0.2 ## triple mutant fitness 1.2 s2 <- ((1 + d2)/(1 + d1)^2) - 1 s3 <- ( (1 + d3)/((1 + d1)^3 * (1 + s2)^3) ) - 1 w <- allFitnessEffects( data.frame(parent = c("Root", "Root", "Root"), child = c("A", "B", "C"), s = d1, sh = -1, typeDep = "MN"), epistasis = c( "A:B" = s2, "A:C" = s2, "B:C" = s2, "A:B:C" = s3)) ## ----fig.width=4, fig.height=4------------------------------------------- plot(w) ## ------------------------------------------------------------------------ evalAllGenotypes(w, order = FALSE, addwt = TRUE) ## ------------------------------------------------------------------------ wb <- allFitnessEffects( epistasis = c( "A" = d1, "B" = d1, "C" = d1, "A:B" = s2, "A:C" = s2, "B:C" = s2, "A:B:C" = s3)) evalAllGenotypes(wb, order = FALSE, addwt = TRUE) ## ---- fig.width=3, fig.height=3------------------------------------------ plot(wb) ## ------------------------------------------------------------------------ wc <- allFitnessEffects( epistasis = c( "A:-B:-C" = d1, "B:-C:-A" = d1, "C:-A:-B" = d1, "A:B:-C" = d2, "A:C:-B" = d2, "B:C:-A" = d2, "A:B:C" = d3)) evalAllGenotypes(wc, order = FALSE, addwt = TRUE) ## ----fig.width=4--------------------------------------------------------- pancr <- allFitnessEffects( data.frame(parent = c("Root", rep("KRAS", 4), "SMAD4", "CDNK2A", "TP53", "TP53", "MLL3"), child = c("KRAS","SMAD4", "CDNK2A", "TP53", "MLL3", rep("PXDN", 3), rep("TGFBR2", 2)), s = 0.1, sh = -0.9, typeDep = "MN")) plot(pancr) ## ----fig.height = 4------------------------------------------------------ rv1 <- allFitnessEffects(data.frame(parent = c("Root", "A", "KRAS"), child = c("A", "KRAS", "FBXW7"), s = 0.1, sh = -0.01, typeDep = "MN"), geneToModule = c("Root" = "Root", "A" = "EVC2, PIK3CA, TP53", "KRAS" = "KRAS", "FBXW7" = "FBXW7")) plot(rv1, expandModules = TRUE, autofit = TRUE) ## ----fig.height=6-------------------------------------------------------- rv2 <- allFitnessEffects(data.frame(parent = c("Root", "1", "2", "3", "4"), child = c("1", "2", "3", "4", "ELF3"), s = 0.1, sh = -0.01, typeDep = "MN"), geneToModule = c("Root" = "Root", "1" = "APC, FBXW7", "2" = "ATM, FAM123B, PIK3CA, TP53", "3" = "BRAF, KRAS, NRAS", "4" = "SMAD2, SMAD4, SOX9", "ELF3" = "ELF3")) plot(rv2, expandModules = TRUE, autofit = TRUE) ## ------------------------------------------------------------------------ K <- 5 sd <- 0.1 sdp <- 0.15 sp <- 0.05 bauer <- data.frame(parent = c("Root", rep("p", K)), child = c("p", paste0("s", 1:K)), s = c(sd, rep(sdp, K)), sh = c(0, rep(sp, K)), typeDep = "MN") fbauer <- allFitnessEffects(bauer) set.seed(1) ## Use fairly large mutation rate b1 <- oncoSimulIndiv(fbauer, mu = 5e-5, initSize = 1000) ## ----baux1,fig.width=6.5, fig.height=10---------------------------------- par(mfrow = c(3, 1)) ## First, drivers plot(b1, type = "line", addtot = TRUE) plot(b1, type = "stacked") plot(b1, type = "stream") ## ----baux2,fig.width=6.5, fig.height=10---------------------------------- par(mfrow = c(3, 1)) ## Next, genotypes plot(b1, show = "genotypes", type = "line") plot(b1, show = "genotypes", type = "stacked") plot(b1, show = "genotypes", type = "stream") ## ----fig.width=6--------------------------------------------------------- set.seed(456) nd <- 70 np <- 5000 s <- 0.1 sp <- 1e-3 spp <- -sp/(1 + sp) mcf1 <- allFitnessEffects(noIntGenes = c(rep(s, nd), rep(spp, np)), drv = seq.int(nd)) mcf1s <- oncoSimulIndiv(mcf1, model = "McFL", mu = 1e-7, detectionSize = 1e8, detectionDrivers = 100, sampleEvery = 0.02, keepEvery = 2, initSize = 2000, finalTime = 1000, onlyCancer = FALSE) summary(mcf1s) ## ----mcf1sx1,fig.width=6.5, fig.height=10-------------------------------- par(mfrow = c(2, 1)) ## I use thinData to make figures smaller and faster plot(mcf1s, addtot = TRUE, lwdClone = 0.9, log = "", thinData = TRUE, thinData.keep = 0.5) ## I also use here xlim, to focus only on a part of the ## data (and make it plot faster) plot(mcf1s, show = "drivers", type = "stacked", thinData = TRUE, thinData.keep = 0.5, xlim = c(600, 1000), legend.ncols = 2) ## ----eval=FALSE---------------------------------------------------------- ## ## set.seed(123) ## nd <- 70 ## np <- 50000 ## s <- 0.1 ## sp <- 1e-4 ## as we have many more passengers ## spp <- -sp/(1 + sp) ## mcfL <- allFitnessEffects(noIntGenes = c(rep(s, nd), rep(spp, np)), ## drv = seq.int(nd)) ## mcfLs <- oncoSimulIndiv(mcfL, ## model = "McFL", ## mu = 1e-7, ## detectionSize = 1e8, ## detectionDrivers = 100, ## sampleEvery = 0.02, ## keepEvery = 2, ## initSize = 1000, ## finalTime = 2000, ## onlyCancer = FALSE) ## ----mcflsx2,fig.width=6------------------------------------------------- data(mcfLs) plot(mcfLs, addtot = TRUE, lwdClone = 0.9, log = "", plotDiversity = TRUE) ## ------------------------------------------------------------------------ summary(mcfLs) ## number of passengers per clone summary(colSums(mcfLs$Genotypes[-(1:70), ])) ## ----mcflsx3------------------------------------------------------------- plot(mcfLs, type = "stacked", thinData = TRUE, thinData.keep = 0.5, plotDiversity = TRUE, xlim = c(0, 1000)) ## ------------------------------------------------------------------------ data(examplesFitnessEffects) names(examplesFitnessEffects) ## ------------------------------------------------------------------------ data(examplesFitnessEffects) evalAllGenotypes(examplesFitnessEffects$cbn1, order = FALSE)[1:10, ] sm <- oncoSimulIndiv(examplesFitnessEffects$cbn1, model = "McFL", mu = 5e-7, detectionSize = 1e8, detectionDrivers = 2, sampleEvery = 0.025, keepEvery = 5, initSize = 2000, onlyCancer = TRUE) summary(sm) ## ------------------------------------------------------------------------ set.seed(1234) evalAllGenotypes(examplesFitnessEffects$cbn1, order = FALSE, model = "Bozic")[1:10, ] sb <- oncoSimulIndiv(examplesFitnessEffects$cbn1, model = "Bozic", mu = 5e-6, detectionSize = 1e8, detectionDrivers = 4, sampleEvery = 2, initSize = 2000, onlyCancer = TRUE) summary(sb) ## ----sbx1,fig.width=6.5, fig.height=3.3---------------------------------- ## Show drivers, line plot par(cex = 0.75, las = 1) plot(sb,show = "drivers", type = "line", addtot = TRUE, plotDiversity = TRUE) ## ----sbx2,fig.width=6.5, fig.height=3.3---------------------------------- ## Drivers, stacked par(cex = 0.75, las = 1) plot(sb,show = "drivers", type = "stacked", plotDiversity = TRUE) ## ----sbx3,fig.width=6.5, fig.height=3.3---------------------------------- ## Drivers, stream par(cex = 0.75, las = 1) plot(sb,show = "drivers", type = "stream", plotDiversity = TRUE) ## ----sbx4,fig.width=6.5, fig.height=3.3---------------------------------- ## Genotypes, line plot par(cex = 0.75, las = 1) plot(sb,show = "genotypes", type = "line", plotDiversity = TRUE) ## ----sbx5,fig.width=6.5, fig.height=3.3---------------------------------- ## Genotypes, stacked par(cex = 0.75, las = 1) plot(sb,show = "genotypes", type = "stacked", plotDiversity = TRUE) ## ----sbx6,fig.width=6.5, fig.height=3.3---------------------------------- ## Genotypes, stream par(cex = 0.75, las = 1) plot(sb,show = "genotypes", type = "stream", plotDiversity = TRUE) ## ----fig.width=6--------------------------------------------------------- set.seed(4321) tmp <- oncoSimulIndiv(examplesFitnessEffects[["o3"]], model = "McFL", mu = 5e-5, detectionSize = 1e8, detectionDrivers = 3, sampleEvery = 0.025, max.num.tries = 10, keepEvery = 5, initSize = 2000, finalTime = 6000, onlyCancer = FALSE) ## ----tmpmx1,fig.width=6.5, fig.height=4.1-------------------------------- par(las = 1, cex = 0.85) plot(tmp, addtot = TRUE, log = "", plotDiversity = TRUE) ## ----tmpmx2,fig.width=6.5, fig.height=4.1-------------------------------- par(las = 1, cex = 0.85) plot(tmp, type = "stacked", plotDiversity = TRUE, ylim = c(0, 5500), legend.ncols = 4) ## ------------------------------------------------------------------------ evalAllGenotypes(examplesFitnessEffects[["o3"]], addwt = TRUE) ## ----tmpmx3,fig.width=6.5, fig.height=10--------------------------------- par(mfrow = c(2, 1)) plot(tmp, show = "genotypes", ylim = c(0, 5500), legend.ncols = 3) plot(tmp, show = "genotypes", type = "line", ylim = c(1, 6000)) ## ------------------------------------------------------------------------ set.seed(15) tmp <- oncoSimulIndiv(examplesFitnessEffects[["o3"]], model = "McFL", mu = 5e-5, detectionSize = 1e8, detectionDrivers = 3, sampleEvery = 0.015, max.num.tries = 10, keepEvery = 5, initSize = 2000, finalTime = 20000, onlyCancer = FALSE, extraTime = 1500) tmp ## ----tmpmdx5,fig.width=6.5, fig.height=4--------------------------------- par(las = 1, cex = 0.85) plot(tmp, addtot = TRUE, log = "", plotDiversity = TRUE) ## ----tmpmdx6,fig.width=6.5, fig.height=4--------------------------------- par(las = 1, cex = 0.85) plot(tmp, type = "stacked", plotDiversity = TRUE, legend.ncols = 4, ylim = c(0, 5200)) ## ----tmpmdx7,fig.width=6.5, fig.height=5.3------------------------------- ## Improve telling appart the most abundant ## genotypes by sorting colors ## differently via breakSortColors ## Modify ncols of legend, so it is legible by not overlapping ## with plot par(las = 1, cex = 0.85) plot(tmp, show = "genotypes", breakSortColors = "distave", plotDiversity = TRUE, legend.ncols = 4, ylim = c(0, 5300)) ## ------------------------------------------------------------------------ i1 <- allFitnessEffects(noIntGenes = c(1)) evalAllGenotypes(i1, order = FALSE, addwt = TRUE, model = "Bozic") i1_b <- oncoSimulIndiv(i1, model = "Bozic") ## ------------------------------------------------------------------------ evalAllGenotypes(i1, order = FALSE, addwt = TRUE, model = "Exp") i1_e <- oncoSimulIndiv(i1, model = "Exp") summary(i1_e) ## ----eval=FALSE---------------------------------------------------------- ## ## Convert the data ## lb1 <- OncoSimulWide2Long(b1) ## ## ## Install the streamgraph package from github and load ## library(devtools) ## devtools::install_github("hrbrmstr/streamgraph") ## library(streamgraph) ## ## ## Stream plot for Genotypes ## sg_legend(streamgraph(lb1, Genotype, Y, Time, scale = "continuous"), ## show=TRUE, label="Genotype: ") ## ## ## Staked area plot and we use the pipe ## streamgraph(lb1, Genotype, Y, Time, scale = "continuous", ## offset = "zero") %>% sg_legend(show=TRUE, label="Genotype: ") ## ------------------------------------------------------------------------ pancrPop <- oncoSimulPop(10, pancr, detectionSize = 1e7, keepEvery = 10, mc.cores = 2) summary(pancrPop) ## ------------------------------------------------------------------------ pancrSPop <- samplePop(pancrPop) pancrSPop ## ------------------------------------------------------------------------ pancrSamp <- oncoSimulSample(10, pancr) pancrSamp ## ------------------------------------------------------------------------ ## This code will only be evaluated under Windows if(.Platform$OS.type == "windows") try(pancrError <- oncoSimulPop(10, pancr, initSize = 1e-5, detectionSize = 1e7, keepEvery = 10, mc.cores = 2)) ## ------------------------------------------------------------------------ ## Do not run under Windows if(.Platform$OS.type != "windows") pancrError <- oncoSimulPop(10, pancr, initSize = 1e-5, detectionSize = 1e7, keepEvery = 10, mc.cores = 2) ## ----eval=FALSE---------------------------------------------------------- ## pancrError[[1]] ## ----eval=FALSE---------------------------------------------------------- ## pancrError[[1]][1] ## ----fig.height=6-------------------------------------------------------- o3init <- allFitnessEffects(orderEffects = c( "M > D > F" = 0.99, "D > M > F" = 0.2, "D > M" = 0.1, "M > D" = 0.9), noIntGenes = c("u" = 0.01, "v" = 0.01, "w" = 0.001, "x" = 0.0001, "y" = -0.0001, "z" = -0.001), geneToModule = c("M" = "m", "F" = "f", "D" = "d") ) oneI <- oncoSimulIndiv(o3init, model = "McFL", mu = 5e-5, finalTime = 500, detectionDrivers = 3, onlyCancer = FALSE, initSize = 1000, keepPhylog = TRUE, initMutant = c("m > u > d") ) plotClonePhylog(oneI, N = 0) ## ospI <- oncoSimulPop(4, o3init, model = "Exp", mu = 5e-5, finalTime = 500, detectionDrivers = 3, onlyCancer = TRUE, initSize = 10, keepPhylog = TRUE, initMutant = c("d > m > z"), mc.cores = 2 ) op <- par(mar = rep(0, 4), mfrow = c(2, 2)) plotClonePhylog(ospI[[1]]) plotClonePhylog(ospI[[2]]) plotClonePhylog(ospI[[3]]) plotClonePhylog(ospI[[4]]) par(op) ossI <- oncoSimulSample(4, o3init, model = "Exp", mu = 5e-5, finalTime = 500, detectionDrivers = 2, onlyCancer = TRUE, initSize = 10, initMutant = c("z > d"), thresholdWhole = 1 ## check presence of initMutant ) ## No phylogeny is kept with oncoSimulSample, but look at the ## OcurringDrivers and the sample ossI$popSample ossI$popSummary[, "OccurringDrivers", drop = FALSE] ## ------------------------------------------------------------------------ set.seed(15) tmp <- oncoSimulIndiv(examplesFitnessEffects[["o3"]], model = "McFL", mu = 5e-5, detectionSize = 1e8, detectionDrivers = 3, sampleEvery = 0.015, max.num.tries = 10, keepEvery = 5, initSize = 2000, finalTime = 20000, onlyCancer = FALSE, extraTime = 1500, keepPhylog = TRUE) tmp ## ------------------------------------------------------------------------ plotClonePhylog(tmp, N = 0) ## ------------------------------------------------------------------------ plotClonePhylog(tmp, N = 1) ## ----pcpkeepx1----------------------------------------------------------- plotClonePhylog(tmp, N = 1, keepEvents = TRUE) ## ------------------------------------------------------------------------ plotClonePhylog(tmp, N = 1, timeEvents = TRUE) ## ----fig.keep="none"----------------------------------------------------- get.adjacency(plotClonePhylog(tmp, N = 1, returnGraph = TRUE)) ## ------------------------------------------------------------------------ set.seed(456) mcf1s <- oncoSimulIndiv(mcf1, model = "McFL", mu = 1e-7, detectionSize = 1e8, detectionDrivers = 100, sampleEvery = 0.02, keepEvery = 2, initSize = 2000, finalTime = 1000, onlyCancer = FALSE, keepPhylog = TRUE) ## ------------------------------------------------------------------------ plotClonePhylog(mcf1s, N = 1) ## ------------------------------------------------------------------------ plotClonePhylog(mcf1s, N = 1, timeEvents = TRUE) ## ------------------------------------------------------------------------ plotClonePhylog(mcf1s, N = 1, t = c(800, 1000)) ## ------------------------------------------------------------------------ plotClonePhylog(mcf1s, N = 1, t = c(900, 1000), timeEvents = TRUE) ## ----fig.keep="none"----------------------------------------------------- g1 <- plotClonePhylog(mcf1s, N = 1, t = c(900, 1000), returnGraph = TRUE) ## ------------------------------------------------------------------------ plot(g1) ## ----eval=FALSE---------------------------------------------------------- ## op <- par(ask = TRUE) ## while(TRUE) { ## tmp <- oncoSimulIndiv(smn1, model = "McFL", ## mu = 5e-5, finalTime = 500, ## detectionDrivers = 3, ## onlyCancer = FALSE, ## initSize = 1000, keepPhylog = TRUE) ## plotClonePhylog(tmp, N = 0) ## } ## par(op) ## ------------------------------------------------------------------------ oi <- allFitnessEffects(orderEffects = c("F > D" = -0.3, "D > F" = 0.4), noIntGenes = rexp(5, 10), geneToModule = c("F" = "f1, f2, f3", "D" = "d1, d2") ) oiI1 <- oncoSimulIndiv(oi, model = "Exp") oiP1 <- oncoSimulPop(4, oi, keepEvery = 10, mc.cores = 2, keepPhylog = TRUE) ## ----fig.height=9-------------------------------------------------------- op <- par(mar = rep(0, 4), mfrow = c(2, 1)) plotClonePhylog(oiP1[[1]]) plotClonePhylog(oiP1[[2]]) par(op) ## ----fig.height=3-------------------------------------------------------- ## Node 2 and 3 depend on 1, and 4 depends on no one p1 <- cbind(c(1L, 1L, 0L), c(2L, 3L, 4L)) plotPoset(p1, addroot = TRUE) ## ----fig.height=3-------------------------------------------------------- ## A simple way to create a poset where no gene (in a set of 15) depends ## on any other. p4 <- cbind(0L, 15L) plotPoset(p4, addroot = TRUE) ## ----fig.height=3-------------------------------------------------------- pancreaticCancerPoset <- cbind(c(1, 1, 1, 1, 2, 3, 4, 4, 5), c(2, 3, 4, 5, 6, 6, 6, 7, 7)) storage.mode(pancreaticCancerPoset) <- "integer" plotPoset(pancreaticCancerPoset, names = c("KRAS", "SMAD4", "CDNK2A", "TP53", "MLL3","PXDN", "TGFBR2")) ## ----echo=FALSE,results='hide',error=FALSE------------------------------- options(width=60) ## ------------------------------------------------------------------------ ## use poset p1101 data(examplePosets) p1101 <- examplePosets[["p1101"]] ## Bozic Model b1 <- oncoSimulIndiv(p1101, keepEvery = 15) summary(b1) ## ----pb2bothx1,fig.height=5.5, fig.width=5.5----------------------------- b2 <- oncoSimulIndiv(p1101, keepEvery = 1) summary(b2) plot(b2) ## ----pbssttx1,eval=FALSE------------------------------------------------- ## plot(b2, type = "stacked") ## ----echo=FALSE,eval=TRUE------------------------------------------------ set.seed(1) ## for repeatability. Once I saw it not reach cancer. ## ------------------------------------------------------------------------ m2 <- oncoSimulIndiv(examplePosets[["p1101"]], model = "McFL", numPassengers = 0, detectionDrivers = 8, mu = 5e-7, initSize = 4000, sampleEvery = 0.025, finalTime = 25000, keepEvery = 5, detectionSize = 1e6) ## ----m2x1,fig.width=6.5, fig.height=10----------------------------------- par(mfrow = c(2, 1)) plot(m2, addtot = TRUE, log = "", thinData = TRUE, thinData.keep = 0.5) plot(m2, type = "stacked", thinData = TRUE, thinData.keep = 0.5) ## ------------------------------------------------------------------------ b3 <- oncoSimulIndiv(p1101, onlyCancer = FALSE) summary(b3) b4 <- oncoSimulIndiv(p1101, onlyCancer = FALSE) summary(b4) ## ----b3b4x1ch1, fig.width=8, fig.height=4-------------------------------- par(mfrow = c(1, 2)) par(cex = 0.8) ## smaller font plot(b3) plot(b4) ## ----ch2----------------------------------------------------------------- p1 <- oncoSimulPop(4, p1101, mc.cores = 2) par(mfrow = c(2, 2)) plot(p1, ask = FALSE) ## ----p1multx1,eval=FALSE------------------------------------------------- ## par(mfrow = c(2, 2)) ## plot(p1, type = "stream", ask = FALSE) ## ----p1multstx1,eval=FALSE----------------------------------------------- ## par(mfrow = c(2, 2)) ## plot(p1, type = "stacked", ask = FALSE) ## ------------------------------------------------------------------------ m1 <- oncoSimulPop(100, examplePosets[["p1101"]], numPassengers = 0, mc.cores = 2) ## ------------------------------------------------------------------------ genotypes <- samplePop(m1) ## ----fxot1,fig.width=4, fig.height=4------------------------------------- colSums(genotypes)/nrow(genotypes) require(Oncotree) ot1 <- oncotree.fit(genotypes) plot(ot1) ## ----fxot2,fig.width=4, fig.height=4------------------------------------- genotypesSC <- samplePop(m1, typeSample = "single") colSums(genotypesSC)/nrow(genotypesSC) ot2 <- oncotree.fit(genotypesSC) plot(ot2) ## ------------------------------------------------------------------------ ## No seed fixed, so reruns will give different DAGs. (a1 <- simOGraph(10)) library(graph) ## for simple plotting plot(as(a1, "graphNEL")) ## ------------------------------------------------------------------------ sessionInfo()