deleted file mode 100644

@@ -1 +0,0 @@

-UnsupportedPlatforms: perceval, petty

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

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progresion with Epistasis 

-Version: 2.1.2

-Date: 2016-03-27

+Version: 2.1.3

+Date: 2016-04-07

 Authors@R: c(person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),

 		     email = "rdiaz02@gmail.com"),

 	      person("Mark", "Taylor", role = "ctb", email = "ningkiling@gmail.com"))

@@ -500,10 +500,23 @@ allFitnessEffects <- function(rT = NULL,

     }

     if(!is.null(noIntGenes)) {

+        if(inherits(noIntGenes, "character")) {

+            wm <- paste("noIntGenes is a character vector.",

+                        "This is probably not what you want, and will",

+                        "likely result in an error downstream.",

+                        "You can get messages like",

+                        " 'not compatible with requested type', and others.",

+                        "We are stopping.")

+            stop(wm)

+        }

+            

         mg <- max(geneModule[, "GeneNumID"])

         gnum <- seq_along(noIntGenes) + mg

         if(!is.null(names(noIntGenes))) {

             ng <- names(noIntGenes)

+            if( grepl(",", ng, fixed = TRUE) || grepl(">", ng, fixed = TRUE)

+                || grepl(":", ng, fixed = TRUE))

+                stop("The name of some noIntGenes contain a ',' or a '>' or a ':'")

             if(any(ng %in% geneModule[, "Gene"] ))

                 stop("A gene in noIntGenes also present in the other terms")

         } else {

@@ -527,7 +540,18 @@ allFitnessEffects <- function(rT = NULL,

     }

     if(!is.null(drvNames)) {

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

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

+        ## drivers should never be in the geneNoInt; Why!!!???

+        ## Catch the problem. This is an overkill,

+        ## so since we catch the issue, we could leave the geneNoInt. But

+        ## that should not be there in this call.

+        ## if(any(drvNames %in% geneNoInt$Gene)) {

+        ##     stop(paste("At least one gene in drvNames is a geneNoInt gene.",

+        ##                "That is not allowed.",

+        ##                "If that gene is a driver, pass it as gene in the epistasis",

+        ##                "component."))

+        ## }

+        ## drv <- getGeneIDNum(geneModule, NULL, drvNames)

     } else {

         drv <- geneModule$GeneNumID[-1]

     }

@@ -1,7 +1,20 @@

+Changes in version 2.2.0 (for BioC 3.3):

+	- Plots of genotypes.

+	- Stacked area and stream plots (code from Marc Taylor).

+	- Example of modules and no epistasis.

+	- Removed requirement of Root in geneToModule.

+	- More tests (and reorganized them)

+	- Miscell. improvements and typos fixed in documentation and vignette.

+	- Added mutationPropGrowth as argument.

+	- Some minor bug fixes and additional checks for user errors.

+

+Changes in version 2.1.3 (2016-04-04):

+	- Fixed sporadic bug in countDrivers

+

 Changes in version 2.1.2 (2016-03-27):

 	- Arguments to BNB_Algo5 explicit.

 	- Example of modules and no epistasis.

-	- Remove requirement of Root in geneToModule.

+	- Removed requirement of Root in geneToModule.

 	- More tests (and reorganized them)

 	- Miscell. improvements in documentation and vignette.

@@ -71,11 +71,14 @@ allFitnessEffects(rT = NULL, epistasis = NULL, orderEffects = NULL,

 \item{noIntGenes}{

   A numeric vector (optionally named) with the fitness coefficients of genes

   (only genes, not modules) that show no interactions. These genes

-      cannot be part of modules. But you can specify modules that have

-      no epistatic interactions. See examples and vignette.

-    }

-    

-    \item{geneToModule}{

+  cannot be part of modules. But you can specify modules that have

+  no epistatic interactions. See examples and vignette.

+

+  Of course, avoid using potentially confusing characters in the

+  names. In particular, "," and ">" are not allowed as gene names.

+}

+

+\item{geneToModule}{

   A named character vector that allows to match genes and modules. The

   names are the modules, and each of the values is a character vector

@@ -159,11 +162,21 @@ allFitnessEffects(rT = NULL, epistasis = NULL, orderEffects = NULL,

 }

-\note{Please, note that the meaning of the fitness effects in the

+\note{

+  Please, note that the meaning of the fitness effects in the

   McFarland model is not the same as in the original paper; the fitness

   coefficients are transformed to allow for a simpler fitness function

   as a product of terms. This differs with respect to v.1. See the

-  vignette for details.}

+  vignette for details.

+

+  The names of the genes and modules can be fairly arbitrary. But if you

+  try hard you can confuse the parser. For instance, using gene or

+  module names that contain "," or ":", or ">" is likely to get you into

+  trouble. Of course, you know you should not try to use those

+  characters because you know those characters have special meanings to

+  separate names or indicate epistasis or order relationships.  Right

+  now, using those characters as names is caught (and result in

+  stopping) if passed as names for noIntGenes.  }

 \author{ Ramon Diaz-Uriarte

 }

@@ -1,5 +1,4 @@

-

-//     Copyright 2013, 2014, 2015 Ramon Diaz-Uriarte

+//     Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte

 //     This program is free software: you can redistribute it and/or modify

 //     it under the terms of the GNU General Public License as published by

@@ -149,6 +148,7 @@ void remove_zero_sp_nr(std::vector<int>& sp_to_remove,

   }

 }

+

 inline void driverCounts(int& maxNumDrivers,

 			 int& totalPresentDrivers,

 			 std::vector<int>& countByDriver,

@@ -163,24 +163,28 @@ inline void driverCounts(int& maxNumDrivers,

   // We used to do count_NumDrivers and then whichDrivers

   maxNumDrivers = 0;

   int tmpdr = 0;

-  int driver_indx;

-  

+  int driver_indx = 0; // the index in the driver table

   for(int j = 0; j < returnGenotypes.ncol(); ++j) {

     tmpdr = 0;

+    driver_indx = 0;

     for(int i : drv) {

-      driver_indx = i - 1;

-      tmpdr += returnGenotypes(driver_indx, j);

-      countByDriver[driver_indx] += returnGenotypes(driver_indx, j);

+      tmpdr += returnGenotypes(i - 1, j);

+      countByDriver[driver_indx] += returnGenotypes(i - 1, j);

+      ++driver_indx;

     }

     if(tmpdr > maxNumDrivers) maxNumDrivers = tmpdr;

   }

+  if(returnGenotypes.ncol() > 0) {

+    STOPASSERT(driver_indx == static_cast<int>( countByDriver.size()));

+  } else {

+    STOPASSERT(driver_indx <= static_cast<int>( countByDriver.size()));

+  }

   for(size_t i = 0; i < countByDriver.size(); ++i) {

     if(countByDriver[i] > 0) {

       presentDrivers.push_back(i + 1);

       ++totalPresentDrivers;

     }

   }

-  

 }

@@ -1596,10 +1600,11 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

   double totPopSize = 0;

   std::vector<double> sampleTotPopSize;

   std::vector<double> sampleLargestPopSize;

-  std::vector<int> sampleMaxNDr; //The number of drivers in the population

-  // with the largest number of drivers; and this for each time sample

-  std::vector<int> sampleNDrLargestPop; //Number of drivers in population

-  // with largest size (at each time sample)

+  std::vector<int> sampleMaxNDr; //The largest number of drivers in any

+				 //genotype or clone at each  time sample

+  std::vector<int> sampleNDrLargestPop; //Number of drivers in the clone

+					// with largest size (at each time

+					// sample)

   sampleTotPopSize.reserve(initIt);

   sampleLargestPopSize.reserve(initIt);

   sampleMaxNDr.reserve(initIt);

@@ -1827,12 +1832,6 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

 	       countByDriver, presentDrivers,

 	       returnGenotypes, fitnessEffects.drv);

-  // nr_count_NumDrivers(maxNumDrivers, countByDriver,

-  //   		      returnGenotypes, fitnessEffects.drv);

-

-  // whichDrivers(totalPresentDrivers, occurringDrivers, countByDriver);

-

-  // std::map<int, std::string> intName = mapGenesIntToNames(fitnessEffects);

   std::vector<std::string> genotypesAsStrings =

     genotypesToNameString(uniqueGenotypes_vector_nr, genesInFitness, intName);

   std::string driversAsString =

@@ -1850,31 +1849,6 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

   fill_SStats(perSampleStats, sampleTotPopSize, sampleLargestPopSize,

   	      sampleLargestPopProp, sampleMaxNDr, sampleNDrLargestPop);

-

-  // // // debuggin: precompute things

-  // DP2(simulsDone);

-  // DP2(maxWallTime);

-  // DP2(hittedWallTime);

-  // DP2(outNS_i);

-  // DP2( sampleMaxNDr[outNS_i]);

-  // DP2(sampleNDrLargestPop[outNS_i]);

-  // DP2(sampleLargestPopSize[outNS_i]);

-  // DP2(sampleLargestPopProp[outNS_i]);

-  // DP2((runningWallTime > maxWallTime));

-  // here("after precomp");

-  // here("*******************************************");

-

-  // Rcpp::List returnGenotypesO = Rcpp::wrap(uniqueGenotypesV);

-

-  // if(keepPhylog) {

-  //   Rcpp::DataFrame PhylogDF = Rcpp::DataFrame::create(Named("parent") = phylog.parent,

-  // 						       Named("child") = phylog.child,

-  // 						       Named("time") = phylog.time);

-  // } else {

-  //   Rcpp::DataFrame PhylogDF = Rcpp::DataFrame::create(Named("parent") = NA,

-  // 						       Named("child") = NA,

-  // 						       Named("time") = NA);

-  // }

   return

     List::create(Named("pops.by.time") = outNS,

@@ -1920,9 +1894,7 @@ Rcpp::List nr_BNB_Algo5(Rcpp::List rFE,

 					       Named("UnrecoverExcept") = false)

 		 );

-  //  END_RCPP

-    

-    }

+}

 // Creating return object:

@@ -5,11 +5,18 @@

 //#define DEBUGZ

 // #define DEBUGV

-//#define DEBUGW

-

+#define DEBUGW

 #define DP1(x) {Rcpp::Rcout << "\n DEBUG2: I am at " << x << std::endl;}

 #define DP2(x) {Rcpp::Rcout << "\n DEBUG2: Value of " << #x << " = " << x << std::endl;}

+#define DP3(x, t){                                                    \

+    Rcpp::Rcout <<"\n DEBUG2:" ;                                      \

+    for(int xut = 0; xut < t; ++xut) Rcpp::Rcout << "\t ";            \

+    Rcpp::Rcout << "  I am at " << x << std::endl;}

+#define DP4(x, t){                                                    \

+    Rcpp::Rcout <<"\n DEBUG2:" ;                                      \

+    for(int xut = 0; xut < t; ++xut) Rcpp::Rcout << "\t ";            \

+    Rcpp::Rcout << "  Value of " << #x << " = " << x << std::endl; }

 /* void here(std::string x) { */

 /*   Rcpp::Rcout << "\n DEBUG: HERE at " << x << std::endl; */

new file mode 100644

@@ -0,0 +1,59 @@

+cat(paste("\n Starting driverCounts long at", date()))

+

+date()

+test_that("Runs without crashes", {

+    

+    iteration <- 1

+    ni <- runif(10, min = -0.01, max = 0.1)

+    names(ni) <- paste0("g", 2:11)

+    ## each single one

+    for(i in 2:11) {

+        cat("\n doing iteration", iteration, "\n")

+        ni[] <- runif(10, min = -0.01, max = 0.1)

+        ni <- sample(ni)

+        drvN <- paste0("g", i)

+        fe31 <- allFitnessEffects(noIntGenes = ni,

+                                  drvNames = drvN)

+        expect_silent(uu <- oncoSimulPop(20,

+                                         fe31, 

+                                         mu = 1e-6,

+                                         initSize = 1e5,

+                                         model = "McFL",

+                                         detectionSize = 5e6,

+                                         finalTime = 5,

+                                         keepEvery = 1,

+                                         onlyCancer = FALSE,

+                                         mc.cores = 2,

+                                         sampleEvery = 0.03

+                                         ))

+        iteration <- iteration + 1

+    }

+    ## all pairs, trios, etc, shuffled order

+    for(n in 2:10) {

+        m <- combinations(10, n, 2:11)

+        for(j in 1:nrow(m)) {

+            cat("\n doing iteration", iteration, "\n")

+            ni[] <- runif(10, min = -0.01, max = 0.1)

+            ni <- sample(ni)

+            drvN <- paste0("g", sample(m[j, ]))

+            fe31 <- allFitnessEffects(noIntGenes = ni,

+                                  drvNames = drvN)

+            expect_silent(uu <- oncoSimulPop(10,

+                                             fe31, 

+                                             mu = 1e-6,

+                                             initSize = 1e5,

+                                             model = "McFL",

+                                             detectionSize = 5e6,

+                                             finalTime = 5,

+                                             keepEvery = 1,

+                                             onlyCancer = FALSE,

+                                             mc.cores = 2,

+                                             sampleEvery = 0.03

+                                             ))

+            iteration <- iteration + 1

+        }

+    }

+

+})

+date()

+cat(paste("\n Ending driverCounts long at", date()))

@@ -1,5 +1,10 @@

 cat(paste("\n Starting at mutPropGrowth long", date(), "\n"))

+## Note that many of the tests below where we do a test and have a

+## comparison like "whatever$p.value > p.fail " are, of course, expected

+## to fail with prob. ~ p.fail even if things are perfectly OK.

+

+

 ## Why this does not really reflect what we want, and why number of clones

 ## is better that capture the idea of "more mutations". NumClones reflects

 ## the creation of a new clone, something that happens whenever there is a

@@ -50,7 +55,7 @@ mutsPerClone <- function(x, per.pop.mean = TRUE) {

 RNGkind("L'Ecuyer-CMRG") ## for the mclapplies

 ## If crashes I want to see where: thus output seed.

-

+p.value.threshold <- 0.001

 ## this tests takes 10 seconds. Moved to long. So this was in the standard

 ## ones.

@@ -60,7 +65,7 @@ test_that("mutPropGrowth diffs with s> 0", {

     set.seed(pseed)

     cat("\n mgp1: the seed is", pseed, "\n")

     ft <- 4 ## 2.7

-    pops <- 100

+    pops <- 150

     lni <- 150 ## 100

     no <- 1e3 ## 5e1 

     ni <- c(2, rep(0, lni)) ## 2 ## 4 ## 5

@@ -96,10 +101,12 @@ test_that("mutPropGrowth diffs with s> 0", {

     ## summary(mutsPerClone(nca))

     ## summary(mutsPerClone(nca2))

     ## The real comparison

-    expect_true( median(summary(nca)$NumClones) >

-                 median(summary(nca2)$NumClones))

-    expect_true( mean(mutsPerClone(nca)) >

-                 mean(mutsPerClone(nca2)))

+    expect_true( wilcox.test(summary(nca)$NumClones,

+                             summary(nca2)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(mutsPerClone(nca),

+                             mutsPerClone(nca2),

+                             alternative = "greater")$p.value < p.value.threshold)

 })

 cat("\n", date(), "\n")

@@ -127,7 +134,7 @@ test_that("mutPropGrowth no diffs with s = 0", {

     set.seed(pseed)

     cat("\n mgp1ND: the seed is", pseed, "\n")

     ft <- 4 ## 2.7

-    pops <- 200

+    pops <- 250

     lni <- 150 ## 100

     no <- 1e3 ## 5e1 

     ni <- c(0, rep(0, lni)) ## 2 ## 4 ## 5

@@ -224,7 +231,7 @@ test_that("mutPropGrowth no diffs with s = 0, oncoSimulSample", {

     set.seed(pseed)

     cat("\n oss1: the seed is", pseed, "\n")

     ft <- 3.5 ## 4

-    pops <- 200

+    pops <- 250

     lni <- 200 ## 150

     no <- 1e3 ## 5e1 

     ni <- c(0, rep(0, lni)) ## 2 ## 4 ## 5

@@ -290,7 +297,7 @@ test_that("mutPropGrowth no diffs with s = 0, oncoSimulSample, McFL", {

     set.seed(pseed)

     cat("\n ossmcfND1: the seed is", pseed, "\n")

     ft <- 40 ## 4

-    pops <- 200

+    pops <- 250

     lni <- 200 ## 150

     no <- 1e3 ## 5e1 

     ni <- c(0, rep(0, lni)) ## 2 ## 4 ## 5

@@ -351,14 +358,6 @@ test_that("mutPropGrowth no diffs with s = 0, oncoSimulSample, McFL", {

 cat("\n", date(), "\n")

-

-

-

-

-

-

-

-

 ## The tests below can occasionally fail (but that probability decreases

 ## as we increase number of pops), as they should.

 ## Some of these tests take some time. For now, show times.

@@ -388,7 +387,7 @@ test_that("oncoSimulSample Without initmutant and modules", {

     pseed <- sample(9999999, 1)

     set.seed(pseed)

     cat("\n osS: the seed is", pseed, "\n")

-    pops <- 60

+    pops <- 150

     lni <- 1 ## no fitness effects genes

     fni <- 50 ## fitness effects genes

     no <- 1e4 

@@ -450,10 +449,12 @@ test_that("oncoSimulSample Without initmutant and modules", {

     (mutsPerClone2 <- rowSums(b2$popSample))

     summary(mutsPerClone1)

     summary(mutsPerClone2)

-    expect_true( mean(mutsPerClone2) >

-                 mean(mutsPerClone1))

-    expect_true( median(b2$popSummary[, "NumClones"]) >

-                 median(b1$popSummary[, "NumClones"]))

+    expect_true( wilcox.test(mutsPerClone2,

+                             mutsPerClone1,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(b2$popSummary[, "NumClones"],

+                             b1$popSummary[, "NumClones"],

+                             alternative = "greater")$p.value < p.value.threshold)

 })

 cat("\n", date(), "\n")

@@ -463,7 +464,7 @@ test_that("oncoSimulSample Without initmutant and modules, McFL", {

     pseed <- sample(9999999, 1)

     set.seed(pseed)

     cat("\n osSMcFL: the seed is", pseed, "\n")

-    pops <- 60

+    pops <- 150

     lni <- 1 ## no fitness effects genes

     fni <- 50 ## fitness effects genes

     no <- 1e4 ## note we use only 10 in the other example below

@@ -525,14 +526,179 @@ test_that("oncoSimulSample Without initmutant and modules, McFL", {

     (mutsPerClone2 <- rowSums(b2$popSample))

     summary(mutsPerClone1)

     summary(mutsPerClone2)

-    expect_true( mean(mutsPerClone2) >

-                 mean(mutsPerClone1))

-    expect_true( median(b2$popSummary[, "NumClones"]) >

-                 median(b1$popSummary[, "NumClones"]))

+    expect_true( wilcox.test(mutsPerClone2,

+                             mutsPerClone1,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(b2$popSummary[, "NumClones"],

+                             b1$popSummary[, "NumClones"],

+                             alternative = "greater")$p.value < p.value.threshold)

 })

 cat("\n", date(), "\n")

+

+## Same as above, but we stop in detectionSize, so no differences are

+## attributable just to differences in population size.

+

+

+cat("\n", date(), "\n")

+test_that("oncoSimulSample Without initmutant and modules, fixed size", {

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n osSFPS: the seed is", pseed, "\n")

+    pops <- 150

+    lni <- 1 ## no fitness effects genes

+    fni <- 50 ## fitness effects genes

+    no <- 1e4 

+    ft <- 9  #4 

+    s3 <- 2.5 

+    mu <- 1e-5 

+    ## noInt have no fitness effects, but can accumulate mutations

+    ni <- rep(0, lni)

+    ## Those with fitness effects in one module, so

+    ## neither fitness nor mut. rate blow up

+    gn <- paste(paste0("a", 1:fni), collapse = ", ")

+    f3 <- allFitnessEffects(epistasis = c("A" = s3),

+                            geneToModule = c("A" = gn),

+                            noIntGenes = ni)

+    x <- 1e-9 ## so basically anything that appears once

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n osSFPSa: the seed is", pseed, "\n")

+    b1 <- oncoSimulSample(pops,

+                          f3,

+                          mu = mu,

+                          mutationPropGrowth = FALSE,

+                          finalTime =ft,

+                          initSize = no,

+                          onlyCancer = FALSE,

+                          sampleEvery = 0.01,

+                          detectionSize = 6e4,

+                          detectionDrivers = 99,

+                          seed =NULL,

+                          thresholdWhole = x)

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n osSFPSb: the seed is", pseed, "\n")

+    b2 <- oncoSimulSample(pops,

+                         f3,

+                         mu = mu,

+                         mutationPropGrowth = TRUE,

+                         finalTime =ft,

+                         initSize = no,

+                         onlyCancer = FALSE,

+                         sampleEvery = 0.01,

+                          detectionSize = 6e4,

+                          detectionDrivers = 99,

+                          seed =NULL,

+                         thresholdWhole = x)

+    b1$popSummary[1:5, c(1:3, 8:9)]

+    summary(b1$popSummary[, "NumClones"])

+    summary(b1$popSummary[, "TotalPopSize"])

+    b2$popSummary[1:5, c(1:3, 8:9)]

+    summary(b2$popSummary[, "NumClones"])

+    summary(b2$popSummary[, "TotalPopSize"])

+    ## cc1 and cc2 should all be smaller than pops, or you are maxing

+    ## things and not seeing patterns

+    (cc1 <- colSums(b1$popSample))

+    (cc2 <- colSums(b2$popSample))

+    ## Of course, these are NOT really mutationsPerClone: we collapse over

+    ## whole population.

+    (mutsPerClone1 <- rowSums(b1$popSample))

+    (mutsPerClone2 <- rowSums(b2$popSample))

+    summary(mutsPerClone1)

+    summary(mutsPerClone2)

+    expect_true( wilcox.test(mutsPerClone2,

+                             mutsPerClone1,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(b2$popSummary[, "NumClones"],

+                             b1$popSummary[, "NumClones"],

+                             alternative = "greater")$p.value < p.value.threshold)

+})

+cat("\n", date(), "\n")

+

+

+cat("\n", date(), "\n")

+test_that("oncoSimulSample Without initmutant and modules, McFL, fixed size", {

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n osSFPSMcFL: the seed is", pseed, "\n")

+    pops <- 150

+    lni <- 1 ## no fitness effects genes

+    fni <- 50 ## fitness effects genes

+    no <- 1e4 ## note we use only 10 in the other example below

+    ft <- 10  ##4 

+    s3 <- 2.5 

+    mu <- 1e-5 

+    ## noInt have no fitness effects, but can accumulate mutations

+    ni <- rep(0, lni)

+    ## Those with fitness effects in one module, so

+    ## neither fitness nor mut. rate blow up

+    gn <- paste(paste0("a", 1:fni), collapse = ", ")

+    f3 <- allFitnessEffects(epistasis = c("A" = s3),

+                            geneToModule = c("A" = gn),

+                            noIntGenes = ni)

+    x <- 1e-9 ## 1/no

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n osSFPSMcFLa: the seed is", pseed, "\n")

+    b1 <- oncoSimulSample(pops,

+                          f3,

+                          mu = mu,

+                          mutationPropGrowth = FALSE,

+                          finalTime =ft,

+                          initSize = no,

+                          onlyCancer = FALSE,

+                          sampleEvery = 0.01,

+                          detectionSize = 2.5e4,

+                          detectionDrivers = 99,

+                          seed =NULL,

+                          model = "McFL",

+                          thresholdWhole = x)

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n osSFPSMcFLb: the seed is", pseed, "\n")

+    b2 <- oncoSimulSample(pops,

+                         f3,

+                         mu = mu,

+                         mutationPropGrowth = TRUE,

+                         finalTime =ft,

+                         initSize = no,

+                         onlyCancer = FALSE,

+                         sampleEvery = 0.01,

+                          detectionSize = 2.5e4,

+                          detectionDrivers = 99,

+                         seed =NULL,

+                         model = "McFL",

+                         thresholdWhole = x)

+    b1$popSummary[1:5, c(1:3, 8:9)]

+    summary(b1$popSummary[, "NumClones"])

+    summary(b1$popSummary[, "TotalPopSize"])

+    b2$popSummary[1:5, c(1:3, 8:9)]

+    summary(b2$popSummary[, "NumClones"])

+    summary(b2$popSummary[, "TotalPopSize"])

+    ## cc1 and cc2 should all be smaller than pops, or you are maxing

+    ## things and not seeing patterns

+    (cc1 <- colSums(b1$popSample))

+    (cc2 <- colSums(b2$popSample))

+    (mutsPerClone1 <- rowSums(b1$popSample))

+    (mutsPerClone2 <- rowSums(b2$popSample))

+    summary(mutsPerClone1)

+    summary(mutsPerClone2)

+    expect_true( wilcox.test(mutsPerClone2,

+                             mutsPerClone1,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(b2$popSummary[, "NumClones"],

+                             b1$popSummary[, "NumClones"],

+                             alternative = "greater")$p.value < p.value.threshold)

+})

+cat("\n", date(), "\n")

+

+

+

+

+

+

 ## This generally works, but not always. Because: a) starting with a you

 ## can get a mutation in b. Or starting with b you can get a mutation in

 ## a. When either happens is stochastic. And we are also mixing the

@@ -546,7 +712,7 @@ test_that("Ordering of number of clones with mutpropgrowth", {

     pseed <- sample(9999999, 1)

     set.seed(pseed)

     cat("\n omp1: the seed is", pseed, "\n")

-    pops <- 100

+    pops <- 150

     lni <- 200

     no <- 5e3

     ni <- c(5, 2, rep(0, lni))

@@ -590,12 +756,15 @@ test_that("Ordering of number of clones with mutpropgrowth", {

     expect_true(var(summary(nca2)$NumClones) > 1e-4)

     expect_true(var(summary(ncb2)$NumClones) > 1e-4)

     ## The real comparison

-    expect_true( median(summary(nca)$NumClones) >

-                 median(summary(nca2)$NumClones))

-    expect_true( median(summary(nca)$NumClones) >

-                 median(summary(ncb)$NumClones))

-    expect_true( median(summary(ncb)$NumClones) >

-                 median(summary(ncb2)$NumClones))

+    expect_true( wilcox.test(summary(nca)$NumClones,

+                             summary(nca2)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(nca)$NumClones,

+                             summary(ncb)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(ncb)$NumClones,

+                             summary(ncb2)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

 })

 cat("\n", date(), "\n")

@@ -604,7 +773,7 @@ test_that("Ordering of number of clones with mutpropgrowth, McFL", {

     pseed <- sample(9999999, 1)

     set.seed(pseed)

     cat("\n omp2: the seed is", pseed, "\n")

-    pops <- 100

+    pops <- 150

     lni <- 200

     no <- 5e3

     ni <- c(5, 2, rep(0, lni))

@@ -648,12 +817,15 @@ test_that("Ordering of number of clones with mutpropgrowth, McFL", {

     expect_true(var(summary(nca2)$NumClones) > 1e-4)

     expect_true(var(summary(ncb2)$NumClones) > 1e-4)

     ## The real comparison

-    expect_true( median(summary(nca)$NumClones) >

-                 median(summary(nca2)$NumClones))

-    expect_true( median(summary(nca)$NumClones) >

-                 median(summary(ncb)$NumClones))

-    expect_true( median(summary(ncb)$NumClones) >

-                 median(summary(ncb2)$NumClones))

+    expect_true( wilcox.test(summary(nca)$NumClones,

+                             summary(nca2)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(nca)$NumClones,

+                             summary(ncb)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(ncb)$NumClones,

+                             summary(ncb2)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

 })

 ## Psss ideas.  The idea here is that you hit the module with fitness

@@ -665,7 +837,7 @@ test_that("Without initmutant", {

     pseed <- sample(9999999, 1)

     set.seed(pseed)

     cat("\n s3: the seed is", pseed, "\n")

-    pops <- 200

+    pops <- 250

     lni <- 1 ## no fitness effects genes

     fni <- 50 ## fitness effects genes

     no <- 1e3

@@ -704,16 +876,19 @@ test_that("Without initmutant", {

                          seed = NULL, mc.cores = 2)

     ## summary(s3.g)[, c(1, 2, 3, 8, 9)]

     ## summary(s3.ng)[, c(1, 2, 3, 8, 9)]

-    expect_true( mean(mutsPerClone(s3.g)) >

-                 mean(mutsPerClone(s3.ng)))

-    expect_true( median(summary(s3.g)$NumClones) >

-                 median(summary(s3.ng)$NumClones))

+    expect_true( wilcox.test(mutsPerClone(s3.g),

+                             mutsPerClone(s3.ng),

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(s3.g)$NumClones,

+                             summary(s3.ng)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

 gc() 

 })

 cat("\n", date(), "\n")

 cat("\n", date(), "\n")

 test_that("Without initmutant, 2", {

+    

     ## More of the above. Use smaller s2 and smaller mutation, but then to

     ## see it reliably you need large ft and we also increase

     ## init. pop. size.

@@ -722,9 +897,9 @@ test_that("Without initmutant, 2", {

     cat("\n s2: the seed is", pseed, "\n")

     s2 <- 1.0

     ft <- 15

-    pops <- 200

+    pops <- 150

     lni <- 1 ## no fitness effects genes

-    fni <- 50 ## fitness effects genes

+    fni <- 200 ## fitness effects genes

     no <- 1e4

     mu <- 5e-6 ## easier to see

     ## noInt have no fitness effects, but can accumulate mutations

@@ -757,13 +932,16 @@ test_that("Without initmutant, 2", {

                          initSize = no, keepEvery = 1,

                          onlyCancer = FALSE,

                          seed = NULL, mc.cores = 2)

-    ## summary(s2.g)[, c(1, 2, 3, 8, 9)]

-    ## summary(s2.ng)[, c(1, 2, 3, 8, 9)]

-    expect_true( mean(mutsPerClone(s2.g)) >

-                 mean(mutsPerClone(s2.ng)))

-    expect_true( median(summary(s2.g)$NumClones) >

-                 median(summary(s2.ng)$NumClones))

+    summary(s2.g)[, c(1, 2, 3, 8, 9)]

+    summary(s2.ng)[, c(1, 2, 3, 8, 9)]

+    expect_true( wilcox.test(mutsPerClone(s2.g),

+                             mutsPerClone(s2.ng),

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(s2.g)$NumClones,

+                             summary(s2.ng)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

 gc() 

+

 })

 cat("\n", date(), "\n")

@@ -774,7 +952,7 @@ test_that("McFL: Without initmutant", {

     cat("\n mcfls2: the seed is", pseed, "\n")

     s2 <- 2.0

     ft <- 250

-    pops <- 200 ## 200

+    pops <- 200

     lni <- 1 ## no fitness effects genes

     fni <- 50 ## fitness effects genes

     no <- 1e3

@@ -811,10 +989,12 @@ test_that("McFL: Without initmutant", {

                          seed = NULL, mc.cores = 2)

     ## summary(s2.g)[, c(1, 2, 3, 8, 9)]

     ## summary(s2.ng)[, c(1, 2, 3, 8, 9)]

-    expect_true( mean(mutsPerClone(s2.g)) >

-                 mean(mutsPerClone(s2.ng)))

-    expect_true( median(summary(s2.g)$NumClones) >

-                 median(summary(s2.ng)$NumClones))

+    expect_true( wilcox.test(mutsPerClone(s2.g),

+                             mutsPerClone(s2.ng),

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(s2.g)$NumClones,

+                             summary(s2.ng)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

     ## summary(mutsPerClone(s2.g))

     ## summary(mutsPerClone(s2.ng))

     ## summary(summary(s2.g)$NumClones)

@@ -824,6 +1004,216 @@ gc()

 cat("\n", date(), "\n")

+

+

+

+### As before, but fixing final population size.

+cat("\n", date(), "\n")

+test_that("detectionSize. Without initmutant", {

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n s3FPS: the seed is", pseed, "\n")

+    pops <- 200

+    lni <- 1 ## no fitness effects genes

+    fni <- 50 ## fitness effects genes

+    no <- 1e3

+    ft <- 10 ## 5

+    s3 <- 3.0

+    mu <- 5e-5 ## easier to see

+    ## noInt have no fitness effects, but can accumulate mutations

+    ni <- rep(0, lni)

+    ## Those with fitness effects in one module, so

+    ## neither fitness nor mut. rate blow up

+    gn <- paste(paste0("a", 1:fni), collapse = ", ")

+    f3 <- allFitnessEffects(epistasis = c("A" = s3),

+                            geneToModule = c("A" = gn),

+                            noIntGenes = ni)

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n s3FPSa: the seed is", pseed, "\n")

+    s3.ng <- oncoSimulPop(pops,

+                          f3,

+                          mu = mu,

+                          mutationPropGrowth = FALSE,

+                          detectionDrivers = 9999,

+                          detectionSize = 5e5,

+                          sampleEvery = 0.01,

+                          keepEvery = 1,

+                          finalTime =ft,

+                          initSize = no,

+                          onlyCancer = FALSE,

+                          seed = NULL, mc.cores = 2)

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n s3FPSb: the seed is", pseed, "\n")

+    s3.g <- oncoSimulPop(pops,

+                         f3,

+                         mu = mu,

+                         mutationPropGrowth = TRUE,

+                         detectionDrivers = 9999,

+                          detectionSize = 5e5,

+                          sampleEvery = 0.01,

+                          keepEvery = 1,

+                         finalTime =ft,

+                         initSize = no,

+                         onlyCancer = FALSE,

+                         seed = NULL, mc.cores = 2)

+    ## summary(s3.g)[, c(1, 2, 3, 8, 9)]

+    ## summary(s3.ng)[, c(1, 2, 3, 8, 9)]

+    expect_true( wilcox.test(mutsPerClone(s3.g),

+                             mutsPerClone(s3.ng),

+                             alternative = "greater")$p.value < p.value.threshold)

+    expect_true( wilcox.test(summary(s3.g)$NumClones,

+                             summary(s3.ng)$NumClones,

+                             alternative = "greater")$p.value < p.value.threshold)

+    summary(summary(s3.g)[, 2])

+    summary(summary(s3.ng)[, 2])

+    

+gc() 

+})

+cat("\n", date(), "\n")

+

+cat("\n", date(), "\n")

+test_that("detectionSize. Without initmutant, 2", {

+    

+    ## More of the above. Use smaller s2 and smaller mutation, but then to

+    ## see it reliably you need large final popsize

+    ## Can fail sometimes because differences are small

+    pseed <- sample(9999999, 1)

+    set.seed(pseed)

+    cat("\n s2FPS: the seed is", pseed, "\n")

+    s2 <- 1.0

+    ft <- 500  ## very large, so we stop on size

+    pops <- 200 

+    fni <- 300 ## fitness effects genes

+    no <- 1e3

+    mu <- 5e-6 

+    ## noInt have no fitness effects, but can accumulate mutations

+    ## Those with fitness effects in one module, so

+    ## neither fitness nor mut. rate blow up

+    gn <- paste(paste0("a", 1:fni), collapse = ", ")