#### 4.1.1: improved and faster intervention tests

ramon diaz-uriarte (at Phelsuma) authored on 24/11/2022 13:24:15
Showing 1 changed files
 ... ... @@ -12772,6 +12772,39 @@ percentage_eliminated = (post_int_tot_pop/pre_int_tot_pop)*100 12772 12772  paste0("The percentage of population has decreased by ", percentage_eliminated, "%") 12773 12773   12774 12774   12775 + 12776 + 12777 + 12778 +### Differences between intervening on the total population or over specific genotypes: when do each occur? {#diff_timing_interv}  12779 + 12780 +Suppose an intervention that happens at time unit 10 (and, for the sake of simplicity, suppose we have set sampleEvery = 1). When, in the "What Happens" you specify something like  12781 + 12782 +{r idif1, eval = FALSE} 12783 + 12784 +N = 0.2 * N 12785 + 12786 + 12787 + 12788 +the total population size at time unit 10 is 0.2 times the population you had at the immediately previous sampling period; in this case, total population size at time 10 will be 0.2 the total population size at time 9. You can easily check this looking at the pops.by.time object (beware if you are not keeping all the sampling period; in case of doubt, and if you want to check this, make sure keepEvery is set to sampleEvery).  12789 + 12790 + 12791 +If you do, instead,  12792 + 12793 +{r idif2, eval = FALSE} 12794 + 12795 +n_A = 0.2 * n_A 12796 + 12797 + 12798 + 12799 +you will not see that n_A at time 10 is 0.2 n_A at time 9. The way the code works is: after we have done all the updates, etc, we change the n_A by the requested one. Thus, n_A at time 10 is not 0.2 the n_A at time 9, but 0.2 the n_A that you would have seen at time 10 had you not done an intervention. 12800 + 12801 + 12802 + 12803 + 12804 + 12805 + 12806 + 12807 + 12775 12808  ## Intervening in Rock-Paper-Scissors model in bacterial comunity {#rockscissorsints} 12776 12809   12777 12810  This example is taken from the example @\ref(rockscissors) inspired by @kerr2002a. Here it is described the relationship between 3 strains of _Escherichia coli_ bacteria, that turns out 

#### 3.99.12

ramon diaz-uriarte (at Phelsuma) authored on 19/10/2022 19:10:17
Showing 1 changed files
 ... ... @@ -23,7 +23,8 @@ classoption: a4paper 23 23  geometry: margin=3cm 24 24  fontsize: 12pt 25 25  bibliography: OncoSimulR.bib 26 -biblio-style: "apalike" 26 +biblio-style: apa 27 +csl: apa.csl 27 28  link-citations: true 28 29  vignette: > 29 30  %\VignetteIndexEntry{OncoSimulR: forward genetic simulation in asexual populations with arbitrary epistatic interactions and a focus on modeling tumor progression.} ... ... @@ -1365,7 +1366,11 @@ starting from BioConductor 3.2 (and, of course, available too from 1365 1366  development versions of BioC). So, if you are using the current stable or 1366 1367  development version of BioConductor, or you grab the sources from GitHub 1367 1368  () you are using what we call 1368 -version 2. **The functionality of version has been removed.** 1369 +version 2. **The functionality of version 1 has been removed.** 1370 + 1371 +Version 3 (for BioConductor 3.13) made frequency dependent fitness available in the stable version.  1372 + 1373 +Version 4 (BioConductor 3.16) introduces interventions and the possibility to specify, separately, birth and death (including frequency dependence). 1369 1374   1370 1375   1371 1376   ... ... @@ -3484,7 +3489,7 @@ Exp and Bozic models, or depends on the population size). This is 3484 3489  also shown in Table \@ref(tab:osrfeatures), in the rows for "Fitness 3485 3490  components", under "Evolutionary Features". 3486 3491   3487 -In the case of frequency-dependent fitness simulations (\@ref(fdf)), the 3492 +In the case of frequency-dependent fitness simulations (see section \@ref(fdf)), the 3488 3493  fitness effects must be reevaluated frequently so that birth rate, death 3489 3494  rate, or both, depending the model used, are updated. To do this it is 3490 3495  necessary to use a short step to reevaluate fitness; this is done using a ... ... @@ -3565,7 +3570,7 @@ as possible. Thus, there are two main ways of specifying fitness effects: 3565 3570   3566 3571  * Explicitly passing to OncoSimulR a mapping of genotypes to 3567 3572  fitness. Here you specify the fitness of each genotype. We will refer to 3568 - this as the **explicit mapping of genotypes to fitness**. 3573 + this as the **explicit mapping of genotypes to fitness**. This includes frequency-dependent fitness (section \@ref(fdf)). 3569 3574   3570 3575   3571 3576   ... ... @@ -9856,7 +9861,13 @@ s2fd <- oncoSimulIndiv(afe4, 9856 9861  seed = NULL,  9857 9862  errorHitMaxTries = FALSE,  9858 9863  errorHitWallTime = FALSE) 9859 -plot(s2fd, show = "genotypes") 9864 + 9865 +## In the Mac ARM64 architecture, the above 9866 +## run leads to an exception, which is really odd. 9867 +## While that is debugged, use try to prevent 9868 +## failure of the plot to abort vignette building. 9869 + 9870 +try(plot(s2fd, show = "genotypes")) 9860 9871   9861 9872   9862 9873  {r echo=FALSE}

#### 3.99.10; faster vignette, onlyCancer=FALSE, plot fails with sensible error if simulation had exception

ramon diaz-uriarte (at Phelsuma) authored on 13/10/2022 13:32:19
Showing 1 changed files
 ... ... @@ -470,7 +470,7 @@ A possible way to examine that question would involve: 470 470  - comparing the inferred DAG with the original, true, one. 471 471   472 472   473 -{r, echo=FALSE} 473 +{r rzx0, echo=FALSE} 474 474  set.seed(2) 475 475  RNGkind("L'Ecuyer-CMRG") 476 476   ... ... @@ -485,6 +485,7 @@ g1 <- simOGraph(4, out = "rT") 485 485   486 486  ## Simulate 10 evolutionary trajectories 487 487  s1 <- oncoSimulPop(10, allFitnessEffects(g1, drvNames = 1:4), 488 + onlyCancer = TRUE, 488 489  mc.cores = 2, ## adapt to your hardware 489 490  seed = NULL) ## for reproducibility of vignette 490 491   ... ... @@ -658,6 +659,7 @@ Magellan_stats(rl) ## (Commented out here to avoid writing files) 658 659  ## Simulate evolution in that landscape many times (here just 10) 659 660  simulrl <- oncoSimulPop(10, allFitnessEffects(genotFitness = rl), 660 661  keepPhylog = TRUE, keepEvery = 1, 662 + onlyCancer = TRUE, 661 663  initSize = 4000, 662 664  seed = NULL, ## for reproducibility 663 665  mc.cores = 2) ## adapt to your hardware ... ... @@ -1180,14 +1182,14 @@ loaded. We also explicitly load r Biocpkg("graph") and 1180 1182  for your usual interactive work). And I set the default color for 1181 1183  vertices in igraph. 1182 1184   1183 -{r, results="hide",message=FALSE, echo=TRUE, include=TRUE} 1185 +{r rzx011, results="hide",message=FALSE, echo=TRUE, include=TRUE} 1184 1186  library(OncoSimulR) 1185 1187  library(graph) 1186 1188  library(igraph) 1187 1189  igraph_options(vertex.color = "SkyBlue2") 1188 1190    1189 1191   1190 -{r, echo=FALSE, results='hide'} 1192 +{r rzx02, echo=FALSE, results='hide'} 1191 1193  options(width = 68) 1192 1194    1193 1195   ... ... @@ -1203,7 +1205,7 @@ packageVersion("OncoSimulR") 1203 1205  Following \@ref(steps) we will run two very minimal examples. First a 1204 1206  model with a few genes and **epistasis**: 1205 1207   1206 -{r, fig.width=6.5, fig.height=10} 1208 +{r f7b3v, fig.width=6.5, fig.height=10} 1207 1209  ## 1. Fitness effects: here we specify an  1208 1210  ## epistatic model with modules. 1209 1211  sa <- 0.1 ... ... @@ -1283,6 +1285,8 @@ plot(pancr) 1283 1285   1284 1286  set.seed(1) ## Fix the seed, so we can repeat it 1285 1287  ## We set a small finalTime to speed up the vignette 1288 + 1289 + 1286 1290  ep2 <- oncoSimulIndiv(pancr, model = "McFL", 1287 1291  mu = 1e-6, 1288 1292  sampleEvery = 0.02, ... ... @@ -1291,6 +1295,7 @@ ep2 <- oncoSimulIndiv(pancr, model = "McFL", 1291 1295  finalTime = 20000, 1292 1296  detectionDrivers = 3, 1293 1297  onlyCancer = FALSE) 1298 + 1294 1299    1295 1300   1296 1301   ... ... @@ -1558,10 +1563,11 @@ Nindiv <- 100 ## Number of simulations run. 1558 1563   1559 1564  ## keepEvery = 1 1560 1565  t_exp1 <- system.time( 1561 - exp1 <- oncoSimulPop(Nindiv, pancr,  1562 - detectionProb = "default",  1563 - detectionSize = NA, 1564 - detectionDrivers = NA, 1566 + exp1 <- oncoSimulPop(Nindiv, pancr, 1567 + onlyCancer = TRUE, 1568 + detectionProb = "default",  1569 + detectionSize = NA, 1570 + detectionDrivers = NA, 1565 1571  finalTime = NA, 1566 1572  keepEvery = 1, 1567 1573  model = "Exp",  ... ... @@ -1569,10 +1575,11 @@ t_exp1 <- system.time( 1569 1575   1570 1576   1571 1577  t_mc1 <- system.time( 1572 - mc1 <- oncoSimulPop(Nindiv, pancr,  1573 - detectionProb = "default",  1574 - detectionSize = NA, 1575 - detectionDrivers = NA, 1578 + mc1 <- oncoSimulPop(Nindiv, pancr, 1579 + onlyCancer = TRUE, 1580 + detectionProb = "default",  1581 + detectionSize = NA, 1582 + detectionDrivers = NA, 1576 1583  finalTime = NA, 1577 1584  keepEvery = 1,  1578 1585  model = "McFL",  ... ... @@ -1580,9 +1587,10 @@ t_mc1 <- system.time( 1580 1587   1581 1588  ## keepEvery = NA 1582 1589  t_exp2 <- system.time( 1583 - exp2 <- oncoSimulPop(Nindiv, pancr,  1584 - detectionProb = "default",  1585 - detectionSize = NA, 1590 + exp2 <- oncoSimulPop(Nindiv, pancr, 1591 + onlyCancer = TRUE, 1592 + detectionProb = "default",  1593 + detectionSize = NA, 1586 1594  detectionDrivers = NA, 1587 1595  finalTime = NA, 1588 1596  keepEvery = NA,  ... ... @@ -1591,9 +1599,10 @@ t_exp2 <- system.time( 1591 1599   1592 1600   1593 1601  t_mc2 <- system.time( 1594 - mc2 <- oncoSimulPop(Nindiv, pancr,  1595 - detectionProb = "default",  1596 - detectionSize = NA, 1602 + mc2 <- oncoSimulPop(Nindiv, pancr, 1603 + onlyCancer = TRUE, 1604 + detectionProb = "default",  1605 + detectionSize = NA, 1597 1606  detectionDrivers = NA, 1598 1607  finalTime = NA, 1599 1608  keepEvery = NA, ... ... @@ -1705,11 +1714,12 @@ running time and space consumption. 1705 1714  {r timing2, eval = FALSE} 1706 1715  t_exp3 <- system.time( 1707 1716  exp3 <- oncoSimulPop(Nindiv, pancr,  1708 - detectionProb = c(PDBaseline = 5e4, 1709 - p2 = 0.1, n2 = 5e5, 1710 - checkSizePEvery = 20),  1711 - detectionSize = NA, 1712 - detectionDrivers = NA, 1717 + onlyCancer = TRUE, 1718 + detectionProb = c(PDBaseline = 5e4, 1719 + p2 = 0.1, n2 = 5e5, 1720 + checkSizePEvery = 20), 1721 + detectionSize = NA, 1722 + detectionDrivers = NA, 1713 1723  finalTime = NA, 1714 1724  keepEvery = 1,  1715 1725  model = "Exp",  ... ... @@ -1717,11 +1727,12 @@ t_exp3 <- system.time( 1717 1727   1718 1728  t_exp4 <- system.time( 1719 1729  exp4 <- oncoSimulPop(Nindiv, pancr,  1720 - detectionProb = c(PDBaseline = 5e4, 1721 - p2 = 0.1, n2 = 5e5, 1722 - checkSizePEvery = 20),  1723 - detectionSize = NA, 1724 - detectionDrivers = NA, 1730 + onlyCancer = TRUE, 1731 + detectionProb = c(PDBaseline = 5e4, 1732 + p2 = 0.1, n2 = 5e5, 1733 + checkSizePEvery = 20), 1734 + detectionSize = NA, 1735 + detectionDrivers = NA, 1725 1736  finalTime = NA, 1726 1737  keepEvery = NA,  1727 1738  model = "Exp",  ... ... @@ -1731,10 +1742,11 @@ t_exp4 <- system.time( 1731 1742   1732 1743  t_exp5 <- system.time( 1733 1744  exp5 <- oncoSimulPop(Nindiv, pancr,  1734 - detectionProb = c(PDBaseline = 5e5, 1735 - p2 = 0.1, n2 = 5e7),  1736 - detectionSize = NA, 1737 - detectionDrivers = NA, 1745 + onlyCancer = TRUE, 1746 + detectionProb = c(PDBaseline = 5e5, 1747 + p2 = 0.1, n2 = 5e7), 1748 + detectionSize = NA, 1749 + detectionDrivers = NA, 1738 1750  finalTime = NA, 1739 1751  keepEvery = 1,  1740 1752  model = "Exp",  ... ... @@ -1742,10 +1754,11 @@ t_exp5 <- system.time( 1742 1754   1743 1755  t_exp6 <- system.time( 1744 1756  exp6 <- oncoSimulPop(Nindiv, pancr,  1745 - detectionProb = c(PDBaseline = 5e5, 1746 - p2 = 0.1, n2 = 5e7),  1747 - detectionSize = NA, 1748 - detectionDrivers = NA, 1757 + onlyCancer = TRUE, 1758 + detectionProb = c(PDBaseline = 5e5, 1759 + p2 = 0.1, n2 = 5e7), 1760 + detectionSize = NA, 1761 + detectionDrivers = NA, 1749 1762  finalTime = NA, 1750 1763  keepEvery = NA,  1751 1764  model = "Exp",  ... ... @@ -3614,13 +3627,13 @@ total of four genotypes, and you have a data frame with genotypes 3614 3627  and fitness, where genoytpes are specified as character vectors, 3615 3628  with mutated genes separated by commas: 3616 3629   3617 -{r} 3630 +{r n73} 3618 3631  m4 <- data.frame(G = c("WT", "A", "B", "A, B"), F = c(1, 2, 3, 4)) 3619 3632    3620 3633   3621 3634  Now, let's give that to the allFitnessEffects function: 3622 3635   3623 -{r} 3636 +{r n74} 3624 3637  fem4 <- allFitnessEffects(genotFitness = m4) 3625 3638    3626 3639  (The message is just telling you what the program guessed you ... ... @@ -3628,14 +3641,14 @@ wanted.) 3628 3641   3629 3642   3630 3643  That's it. You can try to plot that fitnessEffects object 3631 -{r} 3644 +{r n75} 3632 3645  try(plot(fem4)) 3633 3646    3634 3647   3635 3648  In this case, you probably want to plot the fitness landscape.  3636 3649   3637 3650   3638 -{r, fig.width=6.5, fig.height = 6.5} 3651 +{r n76, fig.width=6.5, fig.height = 6.5} 3639 3652  plotFitnessLandscape(evalAllGenotypes(fem4)) 3640 3653    3641 3654   ... ... @@ -3643,14 +3656,14 @@ You can also check what OncoSimulR thinks the fitnesses are, with the 3643 3656  evalAllGenotypes function that we will use repeatedly below 3644 3657  (of course, here we should see the same fitnesses we entered): 3645 3658   3646 -{r} 3659 +{r n78} 3647 3660  evalAllGenotypes(fem4, addwt = TRUE) 3648 3661    3649 3662   3650 3663   3651 3664  And you can plot the fitness landscape: 3652 3665   3653 -{r} 3666 +{r n79} 3654 3667  plotFitnessLandscape(evalAllGenotypes(fem4)) 3655 3668    3656 3669   ... ... @@ -3662,14 +3675,14 @@ contains the fitness values. And you do not even need to specify all the 3662 3675  genotypes: the missing genotypes are assigned a fitness 0 ---except 3663 3676  for the WT genotype which, if missing, is assigned a fitness of 1: 3664 3677   3665 -{r} 3678 +{r n80} 3666 3679  m6 <- cbind(c(1, 1), c(1, 0), c(2, 3)) 3667 3680  fem6 <- allFitnessEffects(genotFitness = m6) 3668 3681  evalAllGenotypes(fem6, addwt = TRUE) 3669 3682  ## plot(fem6) 3670 3683   3671 3684   3672 -{r, fig.width=6.5, fig.height = 6.5} 3685 +{r n81, fig.width=6.5, fig.height = 6.5} 3673 3686  plotFitnessLandscape(evalAllGenotypes(fem6)) 3674 3687    3675 3688   ... ... @@ -3960,7 +3973,7 @@ mutated. The $s_i$ can come from any distribution you want. As an example 3960 3973  let's use three genes. We know there are no order effects, but we will 3961 3974  also see what happens if we examine genotypes as ordered. 3962 3975   3963 -{r} 3976 +{r nx82} 3964 3977   3965 3978  ai1 <- evalAllGenotypes(allFitnessEffects( 3966 3979  noIntGenes = c(0.05, -.2, .1), frequencyDependentFitness = FALSE), order = FALSE) ... ... @@ -3969,11 +3982,11 @@ ai1 <- evalAllGenotypes(allFitnessEffects( 3969 3982   3970 3983  We can easily verify the first results: 3971 3984   3972 -{r} 3985 +{r nx83} 3973 3986  ai1 3974 3987    3975 3988   3976 -{r} 3989 +{r nx84} 3977 3990  all(ai1[, "Fitness"] == c( (1 + .05), (1 - .2), (1 + .1), 3978 3991  (1 + .05) * (1 - .2), 3979 3992  (1 + .05) * (1 + .1), ... ... @@ -3985,7 +3998,7 @@ all(ai1[, "Fitness"] == c( (1 + .05), (1 - .2), (1 + .1), 3985 3998  And we can see that considering the order of mutations (see section 3986 3999  \@ref(oe)) makes no difference: 3987 4000   3988 -{r} 4001 +{r nx85} 3989 4002  (ai2 <- evalAllGenotypes(allFitnessEffects( 3990 4003  noIntGenes = c(0.05, -.2, .1)), order = TRUE, 3991 4004  addwt = TRUE)) ... ... @@ -4019,7 +4032,7 @@ that illustrates the three basic different meanings of convergent arrows 4019 4032  using two parental nodes. We will illustrate it here with three. Suppose 4020 4033  we focus on node "g" in the following figure (we will create it shortly) 4021 4034   4022 -{r, fig.height=4} 4035 +{r nx091, fig.height=4} 4023 4036  data(examplesFitnessEffects) 4024 4037  plot(examplesFitnessEffects[["cbn1"]]) 4025 4038    ... ... @@ -4103,7 +4116,7 @@ using DAGs. 4103 4116   4104 4117  ### DAGs: A first conjunction (AND) example {#cbn1} 4105 4118   4106 -{r} 4119 +{r nx092} 4107 4120   4108 4121  cs <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), 4109 4122  child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), ... ... @@ -4120,12 +4133,12 @@ have any particular order.) 4120 4133   4121 4134   4122 4135  We can get a graphical representation using the default "graphNEL" 4123 -{r, fig.height=3} 4136 +{r nx093, fig.height=3} 4124 4137  plot(cbn1) 4125 4138    4126 4139   4127 4140  or one using "igraph": 4128 -{r, fig.height=5} 4141 +{r nx094, fig.height=5} 4129 4142  plot(cbn1, "igraph") 4130 4143    4131 4144   ... ... @@ -4134,7 +4147,7 @@ plot(cbn1, "igraph") 4134 4147  Since we have a parent and children, the reingold.tilford layout is 4135 4148  probably the best here, so you might want to use that: 4136 4149   4137 -{r, fig.height=5} 4150 +{r nx095, fig.height=5} 4138 4151  library(igraph) ## to make the reingold.tilford layout available 4139 4152  plot(cbn1, "igraph", layout = layout.reingold.tilford) 4140 4153    ... ... @@ -4143,7 +4156,7 @@ plot(cbn1, "igraph", layout = layout.reingold.tilford) 4143 4156   4144 4157  And what is the fitness of all genotypes? 4145 4158   4146 -{r} 4159 +{r nx096} 4147 4160  gfs <- evalAllGenotypes(cbn1, order = FALSE, addwt = TRUE) 4148 4161   4149 4162  gfs[1:15, ] ... ... @@ -4163,7 +4176,7 @@ genotype "a, c" has fitness $(1 + 0.1) (1 - 0.9) = 0.11$. 4163 4176   4164 4177  Let's try a first attempt at a somewhat more complex example, where the 4165 4178  fitness consequences of different genes differ. 4166 -{r} 4179 +{r nx097} 4167 4180   4168 4181  c1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), 4169 4182  child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), ... ... @@ -4207,7 +4220,7 @@ are $> 0$ and, thus, is a way of modeling small deviations from the poset 4207 4220  Note that for those nodes that depend only on "Root" the type of 4208 4221  dependency is irrelevant. 4209 4222   4210 -{r} 4223 +{r nx098} 4211 4224  c1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), 4212 4225  child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), 4213 4226  s = c(0.01, 0.02, 0.03, 0.04, 0.1, 0.1, rep(0.2, 3)), ... ... @@ -4230,7 +4243,7 @@ we got to "a, b, c". What matters is whether or not the genes on which 4230 4243  each of "a", "b", and "c" depend are present or not (I only show the first 4231 4244  10 genotypes) 4232 4245   4233 -{r} 4246 +{r nx099} 4234 4247  gcbn2 <- evalAllGenotypes(cbn2, order = FALSE) 4235 4248  gcbn2[1:10, ] 4236 4249    ... ... @@ -4239,7 +4252,7 @@ gcbn2[1:10, ] 4239 4252  Of course, if we were to look at genotypes but taking into account order 4240 4253  of occurrence of mutations, we would see no differences 4241 4254   4242 -{r} 4255 +{r nx0100} 4243 4256  gcbn2o <- evalAllGenotypes(cbn2, order = TRUE, max = 1956) 4244 4257  gcbn2o[1:10, ] 4245 4258    ... ... @@ -4249,7 +4262,7 @@ if they are more than 256, the default.) 4249 4262   4250 4263  You can check the output and verify things are as they should. For instance: 4251 4264   4252 -{r} 4265 +{r nx0101} 4253 4266  all.equal( 4254 4267  gcbn2[c(1:21, 22, 28, 41, 44, 56, 63 ) , "Fitness"], 4255 4268  c(1.01, 1.02, 0.1, 1.03, 1.04, 0.05, ... ... @@ -4270,7 +4283,7 @@ all.equal( 4270 4283  A particular one that is important to understand is genotype with 4271 4284  mutated genes "c, d, e, g": 4272 4285   4273 -{r} 4286 +{r nx0102} 4274 4287  gcbn2[56, ]  4275 4288  all.equal(gcbn2[56, "Fitness"], 1.03 * 1.04 * 1.2 * 0.10) 4276 4289    ... ... @@ -4283,7 +4296,7 @@ satisfied (and that is why the term for "c" is 0.9). 4283 4296  ### DAGs: A semimonotone or "OR" example {#mn1} 4284 4297   4285 4298  We will reuse the above example, changing the type of relationship: 4286 -{r} 4299 +{r nx0103} 4287 4300   4288 4301  s1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), 4289 4302  child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), ... ... @@ -4297,12 +4310,12 @@ smn1 <- allFitnessEffects(s1) 4297 4310   4298 4311  It looks like this (where edges are shown in blue to denote the 4299 4312  semimonotone relationship): 4300 -{r, fig.height=3} 4313 +{r nx0104, fig.height=3} 4301 4314  plot(smn1) 4302 4315    4303 4316   4304 4317   4305 -{r} 4318 +{r nx0105} 4306 4319  gsmn1 <- evalAllGenotypes(smn1, order = FALSE) 4307 4320   4308 4321    ... ... @@ -4310,7 +4323,7 @@ gsmn1 <- evalAllGenotypes(smn1, order = FALSE) 4310 4323  Having just one parental dependency satisfied is now enough, in contrast 4311 4324  to what happened before. For instance: 4312 4325   4313 -{r} 4326 +{r nx0106} 4314 4327  gcbn2[c(8, 12, 22), ] 4315 4328  gsmn1[c(8, 12, 22), ] 4316 4329   ... ... @@ -4323,7 +4336,7 @@ gsmn1[c(20:21, 28), ] 4323 4336   4324 4337  Again, we reuse the example above, changing the type of relationship: 4325 4338   4326 -{r} 4339 +{r nx0107} 4327 4340   4328 4341  x1 <- data.frame(parent = c(rep("Root", 4), "a", "b", "d", "e", "c"), 4329 4342  child = c("a", "b", "d", "e", "c", "c", rep("g", 3)), ... ... @@ -4337,11 +4350,11 @@ xor1 <- allFitnessEffects(x1) 4337 4350   4338 4351   4339 4352  It looks like this (edges in red to denote the "XOR" relationship): 4340 -{r, fig.height=3} 4353 +{r nx0108, fig.height=3} 4341 4354  plot(xor1) 4342 4355    4343 4356   4344 -{r} 4357 +{r nx0109} 4345 4358   4346 4359  gxor1 <- evalAllGenotypes(xor1, order = FALSE) 4347 4360   ... ... @@ -4352,7 +4365,7 @@ Whenever "c" is present with both "a" and "b", the fitness component 4352 4365  for "c" will be $(1 - 0.1)$. Similarly for "g" (if more than one of 4353 4366  "d", "e", or "c" is present, it will show as $(1 - 0.05)$). For example: 4354 4367   4355 -{r} 4368 +{r nx0110} 4356 4369  gxor1[c(22, 41), ]  4357 4370  c(1.01 * 1.02 * 0.1, 1.03 * 1.04 * 0.05) 4358 4371    ... ... @@ -4368,7 +4381,7 @@ epistatic effects (section 4368 4381   4369 4382  We also see a very different pattern compared to CBN (section \@ref(cbn2)) 4370 4383  here: 4371 -{r} 4384 +{r nx0111} 4372 4385  gxor1[28, ]  4373 4386  1.01 * 1.1 * 1.2 4374 4387    ... ... @@ -4376,7 +4389,7 @@ gxor1[28, ] 4376 4389  as exactly one of the dependencies for both "c" and "g" are satisfied. 4377 4390   4378 4391  But  4379 -{r} 4392 +{r nx0112} 4380 4393  gxor1[44, ]  4381 4394  1.01 * 1.02 * 0.1 * 1.2 4382 4395    ... ... @@ -4386,7 +4399,7 @@ one of its dependencies satisfied). 4386 4399   4387 4400  ### Posets: the three types of relationships {#p3} 4388 4401   4389 -{r} 4402 +{r nx0113} 4390 4403   4391 4404  p3 <- data.frame( 4392 4405  parent = c(rep("Root", 4), "a", "b", "d", "e", "c", "f"), ... ... @@ -4399,18 +4412,18 @@ fp3 <- allFitnessEffects(p3) 4399 4412    4400 4413   4401 4414  This is how it looks like: 4402 -{r, fig.height=3} 4415 +{r nx0114, fig.height=3} 4403 4416  plot(fp3) 4404 4417    4405 4418   4406 4419  We can also use "igraph": 4407 4420   4408 -{r, fig.height=6} 4421 +{r nx0115, fig.height=6} 4409 4422  plot(fp3, "igraph", layout.reingold.tilford) 4410 4423    4411 4424   4412 4425   4413 -{r} 4426 +{r nx0116} 4414 4427   4415 4428  gfp3 <- evalAllGenotypes(fp3, order = FALSE) 4416 4429   ... ... @@ -4420,7 +4433,7 @@ gfp3 <- evalAllGenotypes(fp3, order = FALSE) 4420 4433   4421 4434  Let's look at a few: 4422 4435   4423 -{r} 4436 +{r nx0117} 4424 4437  gfp3[c(9, 24, 29, 59, 60, 66, 119, 120, 126, 127), ] 4425 4438   4426 4439  c(1.01 * 1.1, 1.03 * .05, 1.01 * 1.02 * 0.1, 0.1 * 0.05 * 1.3, ... ... @@ -4464,7 +4477,7 @@ there is a module, "A", made of genes "a1" and "a2", and a module 4464 4477  additional fitness advantage compared to mutating only a single one of 4465 4478  them. We can specify this as: 4466 4479   4467 -{r} 4480 +{r ny0111} 4468 4481  s <- 0.2 4469 4482  sboth <- (1/(1 + s)) - 1 4470 4483  m0 <- allFitnessEffects(data.frame( ... ... @@ -4486,7 +4499,7 @@ single one of them"; see details in section \@ref(epi). 4486 4499   4487 4500  Now, specify it using modules: 4488 4501   4489 -{r} 4502 +{r ny01112} 4490 4503  s <- 0.2 4491 4504  m1 <- allFitnessEffects(data.frame( 4492 4505  parent = c("Root", "A"), ... ... @@ -4538,7 +4551,7 @@ modules, one of them with two genes (see details in section 4538 4551  \@ref(mod-no-epi)) we do not need to pass a "Root" as in 4539 4552   4540 4553   4541 -{r} 4554 +{r ny0113} 4542 4555  fnme <- allFitnessEffects(epistasis = c("A" = 0.1, 4543 4556  "B" = 0.2), 4544 4557  geneToModule = c("A" = "a1, a2", ... ... @@ -4548,7 +4561,7 @@ evalAllGenotypes(fnme, order = FALSE, addwt = TRUE) 4548 4561   4549 4562  but it is also OK to have a "Root" in the geneToModule: 4550 4563   4551 -{r} 4564 +{r ny0114} 4552 4565  fnme2 <- allFitnessEffects(epistasis = c("A" = 0.1, 4553 4566  "B" = 0.2), 4554 4567  geneToModule = c( ... ... @@ -4570,7 +4583,7 @@ modules (capital letters) from genes (lowercase with a number), but this 4570 4583  is not needed. 4571 4584   4572 4585   4573 -{r} 4586 +{r ny0115} 4574 4587  p4 <- data.frame( 4575 4588  parent = c(rep("Root", 4), "A", "B", "D", "E", "C", "F"), 4576 4589  child = c("A", "B", "D", "E", "C", "C", "F", "F", "G", "G"), ... ... @@ -4589,19 +4602,19 @@ fp4m <- allFitnessEffects( 4589 4602   4590 4603  By default, plotting shows the modules: 4591 4604   4592 -{r, fig.height=3} 4605 +{r ny0116, fig.height=3} 4593 4606  plot(fp4m) 4594 4607    4595 4608   4596 4609  but we can show the gene names instead of the module names: 4597 4610   4598 -{r, fig.height=3} 4611 +{r ny0117, fig.height=3} 4599 4612  plot(fp4m, expandModules = TRUE) 4600 4613    4601 4614   4602 4615  or 4603 4616   4604 -{r, fig.height=6} 4617 +{r ny0118, fig.height=6} 4605 4618  plot(fp4m, "igraph", layout = layout.reingold.tilford,  4606 4619  expandModules = TRUE) 4607 4620   ... ... @@ -4609,13 +4622,13 @@ plot(fp4m, "igraph", layout = layout.reingold.tilford, 4609 4622   4610 4623  We obtain the fitness of all genotypes in the usual way: 4611 4624   4612 -{r} 4625 +{r ny0119} 4613 4626  gfp4 <- evalAllGenotypes(fp4m, order = FALSE, max = 1024) 4614 4627    4615 4628   4616 4629  Let's look at a few of those: 4617 4630   4618 -{r} 4631 +{r ny0120} 4619 4632  gfp4[c(12, 20, 21, 40, 41, 46, 50, 55, 64, 92, 4620 4633  155, 157, 163, 372, 632, 828), ] 4621 4634   ... ... @@ -4658,7 +4671,7 @@ its usage (but just limit ourselves to using one gene per module here). 4658 4671  Order effects are specified using a $x > y$, which means that that order 4659 4672  effect is satisfied when module $x$ is mutated before module $y$. 4660 4673   4661 -{r} 4674 +{r ny0121} 4662 4675  o3 <- allFitnessEffects(orderEffects = c( 4663 4676  "F > D > M" = -0.3, 4664 4677  "D > F > M" = 0.4, ... ... @@ -4695,7 +4708,7 @@ M$. The 12th matches$F > D > M$but also$D > M$. Etc. 4695 4708   4696 4709  Consider the following case: 4697 4710   4698 -{r} 4711 +{r ny0123} 4699 4712   4700 4713  ofe1 <- allFitnessEffects( 4701 4714  orderEffects = c("F > D" = -0.3, "D > F" = 0.4), ... ... @@ -4712,7 +4725,7 @@$D$and each$f$belongs to module$F$. 4712 4725   4713 4726  What to expect for cases such as$d1 > f1$or$f1 > d1$is clear, as shown in 4714 4727   4715 -{r} 4728 +{r ny0124} 4716 4729  ag[5:16,] 4717 4730    4718 4731   ... ... @@ -4723,14 +4736,14 @@ and module$F$was mutated second. Similar for$d1 > f1 > f2$or 4723 4736 $f1 > d1 > d2$: those map to$D > F$and$F > D$. We can see the fitness 4724 4737  of those four case in: 4725 4738   4726 -{r} 4739 +{r ny0125} 4727 4740  ag[c(17, 39, 19, 29), ] 4728 4741    4729 4742   4730 4743  and they correspond to the values of those order effects, where$F > D = 4731 4744  (1 - 0.3)$and$D > F = (1 + 0.4)$: 4732 4745   4733 -{r} 4746 +{r ny0126} 4734 4747  ag[c(17, 39, 19, 29), "Fitness"] == c(1.4, 0.7, 1.4, 0.7) 4735 4748    4736 4749   ... ... @@ -4740,14 +4753,14 @@$D > F > D$. But since we are concerned with which one happened first and 4740 4753  which happened second we should expect those two to correspond to the same 4741 4754  fitness, that of pattern$D > F$, as is the case: 4742 4755   4743 -{r} 4756 +{r ny0127} 4744 4757  ag[c(43, 44),] 4745 4758  ag[c(43, 44), "Fitness"] == c(1.4, 1.4) 4746 4759    4747 4760  More generally, that applies to all the patterns that start with one of 4748 4761  the "d" genes: 4749 4762   4750 -{r} 4763 +{r ny0128} 4751 4764  all(ag[41:52, "Fitness"] == 1.4) 4752 4765    4753 4766   ... ... @@ -4755,7 +4768,7 @@ Similar arguments apply to the opposite pattern,$F > D$, which apply to 4755 4768  all the possible gene mutation orders that start with one of the "f" 4756 4769  genes. For example: 4757 4770   4758 -{r} 4771 +{r ny0129} 4759 4772  all(ag[53:64, "Fitness"] == 0.7) 4760 4773    4761 4774   ... ... @@ -4768,7 +4781,7 @@ the above, with five genes (there are 325 genotypes, and that is why we 4768 4781  pass the "max" argument to evalAllGenotypes, to allow for 4769 4782  more than the default 256). 4770 4783   4771 -{r} 4784 +{r ny0130} 4772 4785   4773 4786  ofe2 <- allFitnessEffects( 4774 4787  orderEffects = c("F > D" = -0.3, "D > F" = 0.4), ... ... @@ -4785,7 +4798,7 @@ combination that starts with an "f" gene and contains at least one "d" 4785 4798  genes will have a fitness of$1 - 0.3$. All other genotypes have a 4786 4799  fitness of 1: 4787 4800   4788 -{r} 4801 +{r ny0131} 4789 4802  all(ag2[grep("^d.*f.*", ag2[, 1]), "Fitness"] == 1.4) 4790 4803  all(ag2[grep("^f.*d.*", ag2[, 1]), "Fitness"] == 0.7) 4791 4804  oe <- c(grep("^f.*d.*", ag2[, 1]), grep("^d.*f.*", ag2[, 1])) ... ... @@ -4801,7 +4814,7 @@ more interesting, we name genes so that the ordered names do split nicely 4801 4814  between those with and those without order effects (this, thus, also 4802 4815  serves as a test of messy orders of names). 4803 4816   4804 -{r} 4817 +{r ny0132} 4805 4818   4806 4819  foi1 <- allFitnessEffects( 4807 4820  orderEffects = c("D>B" = -0.2, "B > D" = 0.3), ... ... @@ -4812,14 +4825,14 @@ foi1 <- allFitnessEffects( 4812 4825  You can get a verbose view of what the gene names and modules are (and 4813 4826  their automatically created numeric codes) by: 4814 4827   4815 -{r} 4828 +{r ny0133} 4816 4829  foi1[c("geneModule", "long.geneNoInt")] 4817 4830    4818 4831   4819 4832  We can get the fitness of all genotypes (we set$max = 325$because that 4820 4833  is the number of possible genotypes): 4821 4834   4822 -{r} 4835 +{r ny0134} 4823 4836  agoi1 <- evalAllGenotypes(foi1, max = 325, order = TRUE) 4824 4837  head(agoi1) 4825 4838    ... ... @@ -4827,7 +4840,7 @@ head(agoi1) 4827 4840   4828 4841   4829 4842  Now: 4830 -{r} 4843 +{r ny0135} 4831 4844  rn <- 1:nrow(agoi1) 4832 4845  names(rn) <- agoi1[, 1] 4833 4846   ... ... @@ -4840,7 +4853,7 @@ genotype with only two mutations, one of which is either "A", "C" "E" and 4840 4853  the other is "B" or "D" will have the fitness corresponding to "A", "C" or 4841 4854  "E", respectively: 4842 4855   4843 -{r} 4856 +{r ny0136} 4844 4857  agoi1[grep("^A > [BD]$", names(rn)), "Fitness"] == 1.05 4845 4858  agoi1[grep("^C > [BD]$", names(rn)), "Fitness"] == 0.8 4846 4859  agoi1[grep("^E > [BD]$", names(rn)), "Fitness"] == 1.1 ... ... @@ -4856,7 +4869,7 @@ which occupy rows 230 to 253; fitness should be equal (within numerical 4856 4869  error, because of floating point arithmetic) to the order effect of "D" 4857 4870  before "B" times the other effects $(1 - 0.3) * 1.05 * 0.8 * 1.1 = 0.7392$ 4858 4871   4859 -{r} 4872 +{r ny0137} 4860 4873  all.equal(agoi1[230:253, "Fitness"] , 4861 4874  rep((1 - 0.2) * 1.05 * 0.8 * 1.1, 24)) 4862 4875    ... ... @@ -4890,7 +4903,7 @@ Suppose that the actual numerical values are $s_a = 0.2, s_b = 0.3, s_{ab} 4890 4903  = 0.7$. 4891 4904   4892 4905  We specify the above as follows:  4893 -{r} 4906 +{r ny0138} 4894 4907  sa <- 0.2 4895 4908  sb <- 0.3 4896 4909  sab <- 0.7 ... ... @@ -4915,7 +4928,7 @@ satisfies $(1 + s_{ab}) = (1 + s_a) (1 + s_b) (1 + s_2)$ so 4915 4928  $(1 + s_2) = (1 + s_{ab})/((1 + s_a) (1 + s_b))$ and therefore specify as: 4916 4929   4917 4930   4918 -{r} 4931 +{r ny0139} 4919 4932  s2 <- ((1 + sab)/((1 + sa) * (1 + sb))) - 1 4920 4933   4921 4934  e3 <- allFitnessEffects(epistasis = ... ... @@ -4965,7 +4978,7 @@ conciseness). Note that the mutant for exactly A and C has a fitness that 4965 4978  is the product of the individual terms (so there is no epistasis in that case). 4966 4979   4967 4980   4968 -{r} 4981 +{r ny0140} 4969 4982  sa <- 0.1 4970 4983  sb <- 0.15 4971 4984  sc <- 0.2 ... ... @@ -4996,7 +5009,7 @@ term was just the product. We can try to avoid using the "-", however 4996 5009  specification: 4997 5010   4998 5011   4999 -{r} 5012 +{r ny0141} 5000 5013   5001 5014  sa <- 0.1 5002 5015  sb <- 0.15 ... ... @@ -5023,7 +5036,7 @@ evalAllGenotypes(E3B, order = FALSE, addwt = FALSE) 5023 5036   5024 5037  The above two are, of course, identical: 5025 5038   5026 -{r} 5039 +{r ny0142} 5027 5040  all(evalAllGenotypes(E3A, order = FALSE, addwt = FALSE) ==  5028 5041  evalAllGenotypes(E3B, order = FALSE, addwt = FALSE)) 5029 5042    ... ... @@ -5069,7 +5082,7 @@ both options. 5069 5082  ### Epistasis: modules {#epimod} 5070 5083  There is nothing conceptually new, but we will show an example here: 5071 5084   5072 -{r} 5085 +{r ny0143} 5073 5086   5074 5087  sa <- 0.2 5075 5088  sb <- 0.3 ... ... @@ -5087,7 +5100,7 @@ evalAllGenotypes(em, order = FALSE, addwt = TRUE) 5087 5100   5088 5101  Of course, we can do the same thing without using the "-", as in section \@ref(e2): 5089 5102   5090 -{r} 5103 +{r ny0144} 5091 5104  s2 <- ((1 + sab)/((1 + sa) * (1 + sb))) - 1 5092 5105   5093 5106  em2 <- allFitnessEffects(epistasis = ... ... @@ -5122,7 +5135,7 @@ the modules in the epistasis component but have no term for ":" 5122 5135  (the colon). Let's see an example: 5123 5136   5124 5137   5125 -{r} 5138 +{r ny0145} 5126 5139   5127 5140  fnme <- allFitnessEffects(epistasis = c("A" = 0.1, 5128 5141  "B" = 0.2), ... ... @@ -5137,7 +5150,7 @@ In previous versions these was possible using the longer, still accepted 5137 5150  way of specifying a : with a value of 0, but this is no longer 5138 5151  needed: 5139 5152   5140 -{r} 5153 +{r ny0146} 5141 5154  fnme <- allFitnessEffects(epistasis = c("A" = 0.1, 5142 5155  "B" = 0.2, 5143 5156