@@ -1,8 +1,8 @@

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progresion with Epistasis 

-Version: 1.99.4

-Date: 2015-06-22

+Version: 1.99.5

+Date: 2015-06-25

 Author: Ramon Diaz-Uriarte. Melissa O'Neill for randutils.

 Maintainer: Ramon Diaz-Uriarte <rdiaz02@gmail.com>

 Description: Functions for forward population genetic simulation in

@@ -356,9 +356,25 @@ getNamesID <- function(fp) {

 }

-getDrv <- function(geneModule, geneNoInt, drv) {

-    indicesM <- sort(match( drv, geneModule$Gene))

-    indicesI <- sort(match( drv, geneNoInt$Gene))

+getGeneIDNum <- function(geneModule, geneNoInt, drv, sort = TRUE) {

+    ## It returns the genes, as NumID, in the given vector with names drv

+    ## initMutant uses this, for simplicity, without sorting, but noInt

+    ## are always sorted

+

+    ## Also used for the drivers with sort = TRUE

+

+    ## Yes, we must do it twice because we do not know before hand which

+    ## is which. This makes sure no NA. Period.

+    if(any(is.na( match(drv, c(geneModule$Gene, geneNoInt$Gene))))) {

+        stop(paste("For driver or initMutant you have passed genes",

+                   "not in the fitness table."))

+    }

+    

+    indicesM <- as.vector(na.omit(match( drv, geneModule$Gene)))

+    indicesI <- as.vector(na.omit(sort(match( drv, geneNoInt$Gene))))

+    if(sort) {

+        indicesM <- sort(indicesM)

+    }

     return(c(

         geneModule$GeneNumID[indicesM],

         geneNoInt$GeneNumID[indicesI])

@@ -489,7 +505,7 @@ allFitnessEffects <- function(rT = NULL,

     }

     if(!is.null(drvNames)) {

-        drv <- getDrv(geneModule, geneNoInt, drvNames)

+        drv <- getGeneIDNum(geneModule, geneNoInt, drvNames)

     } else {

         drv <- geneModule$GeneNumID[-1]

     }

@@ -831,7 +847,15 @@ nr_oncoSimul.internal <- function(rFE,

                    "as created, for instance, with function",

                    "allFitnessEffects"))

     if(!is.null(initMutant)) {

-        initMutant <- getDrv(rFE$geneModule, initMutant)

+       if(length(grep(">", initMutant))) {

+            initMutant <- nice.vector.eo(initMutant, ">")

+        } else if(length(grep(",", initMutant))) {

+            initMutant <- nice.vector.eo(initMutant, ",")

+        }

+        initMutant <- getGeneIDNum(rFE$geneModule,

+                             rFE$long.geneNoInt,

+                                   initMutant,

+                                   FALSE)

     } else {

         initMutant <- vector(mode = "integer")

     }

@@ -1,3 +1,6 @@

+Changes in version 1.99.5 (2015-06-25):

+	- Fixed bug in initMutant, added tests, and vignette section.

+

 Changes in version 1.99.4 (2015-06-22):

 	- Tried randutils, from O'Neill. Won't work with gcc-4.6.

 	- bool issue in Windows/gcc-4.6.

@@ -267,10 +267,12 @@ to be detected.  For \code{oncoSimulSample} this can be a vector.

 %% \item{seed}{The seed to use for the C++ random number generator. If not

 %%   passed explicitly (the default) then chosen from R.}

-\item{initMutant}{For v.2, a vector of the mutations of the initial

-  clone for the simulations. For v.1, the single mutation of the initial

-  clone for the simulations. The default (if you pass nothing) is to

-  start the simulation from the wildtype genotype with nothing mutated.}

+\item{initMutant}{For v.2, a string with the mutations of the initial

+  mutant, if any. This is the same format as for

+  \code{\link{evalGenotype}}. For v.1, the single mutation of the

+  initial clone for the simulations. The default (if you pass nothing)

+  is to start the simulation from the wildtype genotype with nothing

+  mutated.}

 \item{max.num.tries}{Only applies when \code{onlyCancer = TRUE}. What is

@@ -857,9 +857,8 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

     Genotypes[0] = createNewGenotype(wtGenotype(),

 				     initMutant,

 				     fitnessEffects,

-				     ran_gen); // FIXME: nr, here. What is a "wt

-					// genotype"? Does it have "0"

-					// mutated, or nothing. Nothing.

+				     ran_gen,

+				     false); 

     if(typeModel == TypeModel::beerenwinkel) {

       popParams[0].death = 1.0; //note same is in McFarland.

@@ -1251,7 +1250,8 @@ static void nr_innerBNB(const fitnessEffectsAll& fitnessEffects,

 	newGenotype = createNewGenotype(Genotypes[nextMutant],

 					newMutations,

 					fitnessEffects,

-					ran_gen);

+					ran_gen,

+					true);

 	// nr_change

 	// newGenotype = Genotypes[nextMutant];

@@ -517,9 +517,13 @@ std::map<int, std::string> mapGenesIntToNames(const fitnessEffectsAll& fe) {

 Genotype createNewGenotype(const Genotype& parent,

 			   const std::vector<int>& mutations,

 			   const fitnessEffectsAll& fe,

-			   std::mt19937& ran_gen

+			   std::mt19937& ran_gen,

 			   //randutils::mt19937_rng& ran_gen

-			   ) {

+			   bool random = true) {

+  // random: if multiple mutations, randomly shuffle the ordered ones?

+  // This is the way to go if chromothripsis, but not if we give an

+  // initial mutant

+  

   Genotype newGenot = parent;

   std::vector<int> tempOrder; // holder for multiple muts if order.

   bool sort_rest = false;

@@ -553,7 +557,9 @@ Genotype createNewGenotype(const Genotype& parent,

   // If there is order but multiple simultaneous mutations

   // (chromothripsis), we randomly insert them

-  if(tempOrder.size() > 1)

+

+  // FIXME: initMutant cannot use this!! we give the order!!!

+  if( (tempOrder.size() > 1) && random)

     shuffle(tempOrder.begin(), tempOrder.end(), ran_gen);

   // The new randutils engine:

   // if(tempOrder.size() > 1)

@@ -174,9 +174,9 @@ void obtainMutations(const Genotype& parent,

 Genotype createNewGenotype(const Genotype& parent,

 			   const std::vector<int>& mutations,

 			   const fitnessEffectsAll& fe,

-			   std::mt19937& ran_gen

+			   std::mt19937& ran_gen,

 			   //randutils::mt19937_rng& ran_gen

-			   );

+			   bool random);

 std::vector<double> evalGenotypeFitness(const Genotype& ge,

 					const fitnessEffectsAll& F);

new file mode 100644

@@ -0,0 +1,162 @@

+test_that("initMutant crashes",

+          {

+              o3 <- allFitnessEffects(orderEffects = c(

+                                          "M > D > F" = 0.99,

+                                          "D > M > F" = 0.2,

+                                          "D > M"     = 0.1,

+                                          "M > D"     = 0.9),

+                                      noIntGenes = c(0.01, 0.01),

+                                      geneToModule =

+                                          c("Root" = "Root",

+                                            "M" = "m",

+                                            "F" = "f",

+                                            "D" = "d") )

+              expect_error(oncoSimulIndiv(o3, model = "McFL",

+                      mu = 5e-5, finalTime = 500,

+                      detectionDrivers = 3,

+                      onlyCancer = FALSE,

+                      initSize = 1000,

+                      keepPhylog = TRUE

+                     , initMutant = c("b, a") ## also "a" and "b" separate

+                                          ))

+              expect_error(oncoSimulIndiv(o3, model = "McFL",

+                      mu = 5e-5, finalTime = 500,

+                      detectionDrivers = 3,

+                      onlyCancer = FALSE,

+                      initSize = 1000,

+                      keepPhylog = TRUE

+                     , initMutant = c("b", "a") ## also "a" and "b" separate

+                                          ))

+              expect_error(oncoSimulIndiv(o3, model = "McFL",

+                      mu = 5e-5, finalTime = 500,

+                      detectionDrivers = 3,

+                      onlyCancer = FALSE,

+                      initSize = 1000,

+                      keepPhylog = TRUE

+                     , initMutant = c("1", "2") ## also "a" and "b" separate

+                                          ))

+          })

+

+

+test_that("initMutant lexicog order",

+          {

+              o3 <- allFitnessEffects(orderEffects = c(

+                                          "M > D > F" = 0.99,

+                                          "D > M > F" = 0.2,

+                                          "D > M"     = 0.1,

+                                          "M > D"     = 0.9),

+                                      noIntGenes = c(0.01, 0.01),

+                                      geneToModule =

+                                          c("Root" = "Root",

+                                            "M" = "m",

+                                            "F" = "f",

+                                            "D" = "d") )

+              tmp <- oncoSimulIndiv(o3, model = "McFL",

+                      mu = 5e-5, finalTime = 500,

+                      detectionDrivers = 3,

+                      onlyCancer = FALSE,

+                      initSize = 1000,

+                      keepPhylog = TRUE

+                     , initMutant = c("d > m")

+                      )

+              cn <- colnames(igraph::get.adjacency(plotClonePhylog(tmp, N = 0,

+                                                                   returnGraph = TRUE),

+                                                   sparse = FALSE))

+              expect_true( "d, m_" %in% cn )

+              expect_true( "d, m, f_" %in% cn )

+              expect_false( "m, d_" %in% cn )

+              expect_false( "m, d, f_" %in% cn )

+          })

+

+

+test_that("initMutant lexicog order with noint",

+          {

+              o3 <- 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, "z" = 0.01),

+                                      geneToModule =

+                                          c("Root" = "Root",

+                                            "M" = "m",

+                                            "F" = "f",

+                                            "D" = "d") )

+              tmp <- oncoSimulIndiv(o3, model = "McFL",

+                      mu = 5e-5, finalTime = 500,

+                      detectionDrivers = 3,

+                      onlyCancer = FALSE,

+                      initSize = 1000,

+                      keepPhylog = TRUE

+                     , initMutant = c("d > m > z")

+                      )

+              cn <- colnames(igraph::get.adjacency(plotClonePhylog(tmp, N = 0,

+                                                                   returnGraph = TRUE),

+                                                   sparse = FALSE))

+              expect_true( "d, m_z" %in% cn )

+              expect_true( "d, m, f_z" %in% cn )

+              expect_false( "m, d_" %in% cn )

+              expect_false( "m, d, f_" %in% cn )

+          })

+

+

+test_that("initMutant non lexicog order",

+          {

+              o3 <- 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, "z" = 0.01),

+                                      geneToModule =

+                                          c("Root" = "Root",

+                                            "M" = "m",

+                                            "F" = "f",

+                                            "D" = "d") )

+              tmp <- oncoSimulIndiv(o3, model = "McFL",

+                      mu = 5e-5, finalTime = 500,

+                      detectionDrivers = 3,

+                      onlyCancer = FALSE,

+                      initSize = 1000,

+                      keepPhylog = TRUE

+                     , initMutant = c("m > d")

+                      )

+              cn <- colnames(igraph::get.adjacency(plotClonePhylog(tmp, N = 0,

+                                                                   returnGraph = TRUE),

+                                                   sparse = FALSE))

+              expect_false( "d, m_" %in% cn )

+              expect_false( "d, m, f_" %in% cn )

+              expect_true( "m, d_" %in% cn )

+              expect_true( "m, d, f_" %in% cn )

+          })

+

+

+test_that("initMutant non lexicog order",

+          {

+              o3 <- 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, "z" = 0.01),

+                                      geneToModule =

+                                          c("Root" = "Root",

+                                            "M" = "m",

+                                            "F" = "f",

+                                            "D" = "d") )

+              tmp <- oncoSimulIndiv(o3, model = "McFL",

+                      mu = 5e-5, finalTime = 500,

+                      detectionDrivers = 3,

+                      onlyCancer = FALSE,

+                      initSize = 1000,

+                      keepPhylog = TRUE

+                     , initMutant = c("m > u > d")

+                      )

+              cn <- colnames(igraph::get.adjacency(plotClonePhylog(tmp, N = 0,

+                                                                   returnGraph = TRUE),

+                                                   sparse = FALSE))

+              expect_false( "d, m_" %in% cn )

+              expect_false( "d, m, f_" %in% cn )

+              expect_true( "m, d_u" %in% cn )

+              expect_true( "m, d, f_u" %in% cn )

+          })

@@ -57,6 +57,9 @@ BiocStyle::latex()

 %%}

+%% FIXME: the homozigos, etc.

+

+

 \bioctitle[\textit{OncoSimulR: genetic simulation with arbitrary

   epistasis}]{OncoSimulR: forward genetic simulation in asexual

   populations with arbitrary epistatic interactions and a focus on modeling

@@ -79,8 +82,12 @@ BiocStyle::latex()

 %%     \gitAbbrevHash)}}

+ %% I use x.y.z and gitinfo does not deal with this well in gitVtags et al.

+%% \date{\gitAuthorDate\ {\footnotesize (Version:\gitVtags, Rev: \gitAbbrevHash)}}

+

 \date{\gitAuthorDate\ {\footnotesize (Rev: \gitAbbrevHash)}}

- 

+

+

 \begin{document}

 \maketitle

@@ -138,10 +145,11 @@ core of the code is implemented in C++, providing for fast execution.

 Finally, to help with simulation studies, code to simulate random graphs

 of the kind often seen in CBN, OTs, etc, is also available.

+\subsection{Key features of OncoSimulR}\label{key}

-Therefore, OncoSimulR is now a very general package for forward genetic

-simulation, with applicability well beyond tumor progression. This is a

-summary of some of the key features:

+As mentioned above, OncoSimulR is now a very general package for forward

+genetic simulation, with applicability well beyond tumor progression. This

+is a summary of some of the key features:

 \begin{itemize}

@@ -177,12 +185,16 @@ summary of some of the key features:

     \item Code in C++ is available (though not yet callable from R) for

       using several other models, including the one from Beerenwinkel and

       collaborators \cite{Beerenwinkel2007b}.

+      

     \item You can use very large numbers of genes (e.g., see an example of

       50000 in section \ref{mcf50070} ).

     \item Simulations are generally very fast as I use C++ to implement

       the BNB algorithm.

+    \item You can obtain the true sequence of events and the phylogenetic

+      relationships between clones.

+      

 \end{itemize}

@@ -244,11 +256,32 @@ need that for your usual interactive work).

 <<results="hide">>=

 library(OncoSimulR)

-packageVersion("OncoSimulR")

 library(graph)

 library(igraph)

 @ 

+To be explicit, what version are we running?

+<<>>=

+packageVersion("OncoSimulR")

+@ 

+

+

+\subsection{A temporary note about versions}\label{versions}

+

+In this vignette and the documentation I often refer to version 1 (v.1)

+and version 2 of OncoSimulR. Version 1 is the version available up to, and

+including, BioConductor v.\ 3.1. Version 2 of OncoSimulR is available for

+the current BioConductor development branch, 3.2. If you look at the

+package version, however, it currently shows as 1.99.x (where x should be

+a number $\ge 4$). That is because of the versioning scheme of

+BioConductor. This will become 2.0 in the next release of BioConductor.

+

+Summary: if you are using the current development version of BioConductor,

+or you grab the sources from github

+(\Burl{https://github.com/rdiaz02/OncoSimul}) you are using what we call

+version 2. If you see a version in the package that says ``1.99.x'', where

+``x'' is any number, you are too.

+

 \section{Specifying fitness effects}\label{specfit}

@@ -1875,6 +1908,50 @@ evalGenotype(randomGenotype(fea), fea, echo = TRUE, verbose = TRUE)

 @ 

+\subsection{Homozygosity, heterozygosity, oncogenes, tumor suppressors}\label{oncog}

+

+We are using what is conceptually a single linear chromosome. However, you

+can use it to model scenarios where the numbers of copies affected matter,

+by properly duplicating the genes. 

+

+Suppose we have a tumor suppressor gene, G, with two copies, one from Mom

+and one from Dad. We can have a table like:

+

+

+\begin{tabular} {c c c}

+  $G_M$ & $G_D$ & Fitness \\

+  \hline

+  wt&wt& 1 \\

+  wt&M& 1 \\

+  M&wt& 1\\

+  M&M& $(1 + s)$\\

+  \hline

+\end{tabular}

+

+where $s > 0$, meaning that you need two hits, one in each copy, to

+trigger the clonal expansion.

+

+

+What about oncogenes? A simple model is that one single hit leads to

+clonal expansion and additional hits lead to no additional changes, as in

+this table for gene O, where again the M or D subscript denotes the copy

+from Mom or from Dad:

+

+\begin{tabular} {c c c}

+  $O_M$ & $O_D$ & Fitness \\

+  \hline

+  wt&wt& 1 \\

+  wt&M&  $(1 + s)$\\

+  M&wt& $(1 + s)$\\

+  M&M& $(1 + s)$\\

+  \hline

+\end{tabular}

+

+If you have multiple copies you can proceed similarly. As you can see,

+these are nothing but special cases of synthetic mortality (\ref{sl}),

+synthetic viability (\ref{sv}) and epistasis (\ref{epi}).

+

+

 \section{Specifying fitness effects: some examples from the literature}\label{litex}

 \subsection{Bauer et al}\label{bauer}

@@ -3227,6 +3304,41 @@ thus in the output you can have two columns with both genes mutated.

 %% more meaningful names.

+\subsection{Can I start the simulation from a specific mutant?}\label{initmut}

+

+You bet. In v.1 you can only give the initial mutant as one with a single

+mutated gene. In version 2, however, you can specify the genotype for the

+initial mutant with the same flexibility as in

+\Rfunction{evalGenotype}. Here we show a couple of examples (we use the

+representation of the phylogeny, discussed in section \ref{phylog} of the

+clones so that you can see which clones appear, and from which).

+

+

+<<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, "z" = 0.01),

+                        geneToModule =

+                            c("Root" = "Root",

+                              "M" = "m",

+                              "F" = "f",

+                              "D" = "d") )

+tmpI <- oncoSimulIndiv(o3init, model = "McFL",

+                       mu = 5e-5, finalTime = 500,

+                       detectionDrivers = 3,

+                       onlyCancer = FALSE,

+                       initSize = 1000,

+                       keepPhylog = TRUE,

+                       initMutant = c("m > u > d")

+                      )

+plotClonePhylog(tmpI, N = 0)

+

+@ 

+

 \section{Showing the true phylogenetic relationships of clones}\label{phylog}

@@ -3258,7 +3370,12 @@ tmp <-  oncoSimulIndiv(examplesFitnessEffects[["o3"]],

 tmp

 @ 

-We can plot the phylogenetic relationships of every clone ever created

+We can plot the phylogenetic relationships\footnote{There are several

+  packages in R devoted to phylogenetic inference and related issues. For

+  instance, \CRANpkg{ape}. I have not used that infrastructure because of

+  our very specific needs and circumstances; for instance, internal nodes

+  are observed, we can have networks instead of trees, and we have no

+  uncertainty about when events occurred.} of every clone ever created

 (with fitness larger than 0 ---clones without viability are never shown):

 <<>>=

@@ -3375,12 +3492,33 @@ key clones.

 %% plotClonePhylog(tmp, FALSE, FALSE, FALSE)

 %% @ 

+

+It is of course quite possible that, especially if we consider few genes,

+our phylogeny will be a network, not a tree, as the same child node can

+have multiple parents. You can play with this example, modified from one

+we saw before (section \ref{mn1}):

+

+<<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)

+@ 

+

+

 \subsection{Phylogenies from multiple runs}\label{phylogmult}

 If you use \Rfunction{oncoSimulPop} you can store and plot the phylogenies

 of the different runs:

 <<>>=

+

 oi <- allFitnessEffects(orderEffects =

                c("F > D" = -0.3, "D > F" = 0.4),

                noIntGenes = rexp(5, 10),

@@ -3397,10 +3535,13 @@ oiP1 <- oncoSimulPop(4, oi,

 @ 

 We will plot the first two:

-<<>>=

-par(mfrow = c(2, 1))

+<<fig.height=9>>=

+

+op <- par(mar = rep(0, 4), mfrow = c(2, 1))

 plotClonePhylog(oiP1[[1]])

 plotClonePhylog(oiP1[[2]])

+par(op)

+

 @ 

@@ -1,15 +1,15 @@

 \usepackage[%

-		shash={7bd55ca},

-		lhash={7bd55caf3447abdf6298a7264394e745f71cbc98},

+		shash={4be528e},

+		lhash={4be528e9d8681d84f611c60f30950975ab76155c},

 		authname={ramon diaz-uriarte (at Bufo)},

 		authemail={rdiaz02@gmail.com},

-		authsdate={2015-06-22},

-		authidate={2015-06-22 19:20:10 +0200},

-		authudate={1434993610},

+		authsdate={2015-06-24},

+		authidate={2015-06-24 09:36:48 +0200},

+		authudate={1435131408},

 		commname={ramon diaz-uriarte (at Bufo)},

 		commemail={rdiaz02@gmail.com},

-		commsdate={2015-06-22},

-		commidate={2015-06-22 19:20:10 +0200},

-		commudate={1434993610},

+		commsdate={2015-06-24},

+		commidate={2015-06-24 09:36:48 +0200},

+		commudate={1435131408},

 		refnames={ (HEAD, master)}

 	]{gitsetinfo}

\ No newline at end of file