@@ -1,9 +1,9 @@

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progresion with Epistasis 

-Version: 1.99.5

-Date: 2015-06-25

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

+Version: 1.99.6

+Date: 2015-09-26

+Author: Ramon Diaz-Uriarte. 

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

 Description: Functions for forward population genetic simulation in

     asexual populations, with special focus on cancer progression.

@@ -19,7 +19,8 @@ Description: Functions for forward population genetic simulation in

     phylogenetic relationships of the clones.

 biocViews: BiologicalQuestion, SomaticMutation

 License: GPL (>= 3)

-URL: https://github.com/rdiaz02/OncoSimul, https://popmodels.cancercontrol.cancer.gov/gsr/packages/oncosimulr/, http://ligarto.org/rdiaz

+URL: https://github.com/rdiaz02/OncoSimul, https://popmodels.cancercontrol.cancer.gov/gsr/packages/oncosimulr/

+BugReports: https://github.com/rdiaz02/OncoSimul/issues

 Depends: R (>= 3.1.0)

 Imports: Rcpp (>= 0.11.1), parallel, data.table, graph, Rgraphviz, gtools, igraph

 Suggests: BiocStyle, knitr, Oncotree, testthat

@@ -34,6 +34,8 @@ oncoSimulSample <- function(Nindiv,

                                 } else {

                                     nd <- (2 : round(0.75 * max(fp)))

                                 }

+                                if(length(nd) == 1) ## for sample

+                                       nd <- c(nd, nd)

                                 sample(nd, Nindiv,

                                        replace = TRUE)

                             },

@@ -67,7 +69,7 @@ oncoSimulSample <- function(Nindiv,

     if(keepPhylog)

         warning(paste("oncoSimulSample does not return the phylogeny",

                       "for now, so there is little point in storing it."))

-    

+

     if(max.num.tries.total < Nindiv)

         stop(paste("You have requested something impossible: ",

                    "max.num.tries.total < Nindiv"))

@@ -84,13 +86,16 @@ oncoSimulSample <- function(Nindiv,

                          detectionSize = detectionSize,

                          detectionDrivers = detectionDrivers)[, -1, drop = FALSE]

+    ## FIXME: we are not triggering an error, just a message. This is on

+    ## purpose, since some of these conditions DO provide useful

+    ## output. Do we want these to be errors?

     f.out.attempts <- function() {

         message("Run out of attempts")

         return(list(

             popSummary = NA,

             popSample = NA,

             HittedMaxTries = TRUE,

-            hittedWallTime = FALSE,

+            HittedWallTime = FALSE,

             UnrecoverExcept = FALSE

         ))    

     }    

@@ -112,7 +117,7 @@ oncoSimulSample <- function(Nindiv,

             popSummary = NA,

             popSample = NA,

             HittedMaxTries = TRUE,

-            hittedWallTime = FALSE,

+            HittedWallTime = FALSE,

             UnrecoverExcept = FALSE

         ))    

     }    

@@ -197,7 +202,7 @@ oncoSimulSample <- function(Nindiv,

             return(f.out.time.cpp())

         } else if( indiv > Nindiv ) {

             if(verbosity > 0)

-                message(paste("Successfully sampled ", Nindiv, " individuals"))

+                message(paste0("Successfully sampled ", Nindiv, " individuals"))

             class(pop) <- "oncosimulpop"

             if(inherits(fp, "fitnessEffects")) {

                 geneNames <- names(getNamesID(fp))

@@ -215,6 +220,7 @@ oncoSimulSample <- function(Nindiv,

                 UnrecoverExcept = FALSE

             ))

         } else if( attemptsLeft <= 0 ) {

+              ## it is very unlikely this will ever happen. 

             return(f.out.attempts())

         } else  if( as.double(difftime(Sys.time(), startTime, units = "secs"))

                    > max.wall.time.total ) {

@@ -400,6 +406,8 @@ oncoSimulIndiv <- function(fp,

     if(mu < 0) {

         stop("mutation rate (mu) is negative")

     }

+    ## We do not test for equality to 0. That might be a weird but

+    ## legitimate case?

     if( (keepEvery > 0) & (keepEvery < sampleEvery)) {

         keepEvery <- sampleEvery

@@ -542,9 +550,12 @@ oncoSimulIndiv <- function(fp,

 }

 summary.oncosimul <- function(object, ...) {

-

-    if(object$HittedWallTime || object$HittedMaxTries ||

-       object$other$UnrecoverExcept)

+    

+    if(object$other$UnrecoverExcept) ## yes, when bailing out from

+                                     ## except. can have just minimal

+                                     ## content

+        return(NA)

+    else if (object$HittedWallTime || object$HittedMaxTries)

         return(NA)

     else {

         tmp <- object[c("NumClones", "TotalPopSize", "LargestClone",

@@ -588,8 +599,8 @@ summary.oncosimulpop <- function(object, ...) {

 }

 print.oncosimulpop <- function(x, ...) {

-    cat("\nPopulation of OncoSimul trajectories of ",

-        length(x), " individuals. Call :\n")

+    cat("\nPopulation of OncoSimul trajectories of",

+        length(x), "individuals. Call :\n")

     print(attributes(x)$call)

     cat("\n")

     print(summary(x))

@@ -737,9 +748,10 @@ plotPoset <- function(x, names = NULL, addroot = FALSE,

         box()

 }

-plotAdjMat <- function(adjmat) {

-    plot(as(adjmat, "graphNEL"))

-}

+## this function seems to never be used

+## plotAdjMat <- function(adjmat) {

+##     plot(as(adjmat, "graphNEL"))

+## }

@@ -827,38 +839,6 @@ plotClonePhylog <- function(x, N = 1, t = "last",

-## plotClonePhylog <- function(x, timeEvent = FALSE,

-##                             showEvents = TRUE,

-##                             fixOverlap = TRUE) {

-##     if(!inherits(x, "oncosimul2"))

-##         stop("Phylogenetic information is only stored with v >=2")

-##     if(nrow(x$other$PhylogDF) == 0)

-##         stop("It seems you run the simulation with keepPhylog= FALSE")

-##     ## requireNamespace("igraph")

-##     df <- x$other$PhylogDF

-##     if(!showEvents) {

-##         df <- df[!duplicated(df[, c(1, 2)]), ]

-##     }

-##     g <- igraph::graph.data.frame(df)

-##     l0 <- igraph::layout.reingold.tilford(g)

-##     if(!timeEvent) {

-##         plot(g, layout = l0)

-##     } else {

-##         l1 <- l0

-##         indexAppear <- match(V(g)$name, as.character(df[, 2]))

-##         firstAppear <- df$time[indexAppear]

-##         firstAppear[1] <- 0

-##         l1[, 2] <- (max(firstAppear) - firstAppear)

-##         if(fixOverlap) {

-##             dx <- which(duplicated(l1[, 1]))

-##             if(length(dx)) {

-##                 ra <- range(l1[, 1])

-##                 l1[dx, 1] <- runif(length(dx), ra[1], ra[2])

-##             }

-##         }

-##         plot(g, layout = l1)         

-##     }

-## }

@@ -866,58 +846,73 @@ plotClonePhylog <- function(x, N = 1, t = "last",

 ############# The rest are internal functions

-get.mut.vector.whole <- function(tmp, timeSample = "last", threshold = 0.5) {

-    ## Obtain, from  results from a simulation run, the vector

-    ## of 0/1 corresponding to each gene.

-    

-    ## threshold is the min. proportion for a mutation to be detected

-    ## We are doing whole tumor sampling here, as in Sprouffske

-    ## timeSample: do we sample at end, or at a time point, chosen

-    ## randomly, from all those with at least one driver?

-    

-    

+get.the.time.for.sample <- function(tmp, timeSample) {

     if(timeSample == "last") {

-        return(as.numeric(

-            (tcrossprod(tmp$pops.by.time[nrow(tmp$pops.by.time), -1],

-                        tmp$Genotypes)/tmp$TotalPopSize) > threshold))

+        if(tmp$TotalPopSize == 0) {

+            warning(paste("Final population size is 0.",

+                          "Thus, there is nothing to sample with",

+                          "sampling last. You will get NAs"))

+            the.time <- -99

+        } else {

+              the.time <- nrow(tmp$pops.by.time)

+          }

     } else if (timeSample %in% c("uniform", "unif")) {

-        the.time <- sample(which(tmp$PerSampleStats[, 4] > 0), 1)

-        pop <- tmp$pops.by.time[the.time, -1]

-        popSize <- tmp$PerSampleStats[the.time, 1]

-        return( as.numeric((tcrossprod(pop,

-                                       tmp$Genotypes)/popSize) > threshold) )

-    }

+          candidate.time <- which(tmp$PerSampleStats[, 4] > 0)

+          

+          if (length(candidate.time) == 0) {

+              warning(paste("There is not a single sampled time",

+                            "at which there are any mutants.",

+                            "Thus, no uniform sampling possible.",

+                            "You will get NAs"))

+              the.time <- -99

+              ## return(rep(NA, nrow(tmp$Genotypes)))

+          } else if (length(candidate.time) == 1) {

+                message("Only one sampled time period with mutants.")

+                the.time <- candidate.time

+            } else {

+                  the.time <- sample(candidate.time, 1)

+              }

+      } else {

+            stop("Unknown timeSample option")

+        }

+    return(the.time)

 }

-get.mut.vector.singlecell <- function(tmp, timeSample = "last") {

-    ## No threshold, as single cell.

-

-    ## timeSample: do we sample at end, or at a time point, chosen

-    ## randomly, from all those with at least one driver?

+get.mut.vector <- function(x, timeSample, typeSample,

+                           thresholdWhole) {

+    the.time <- get.the.time.for.sample(x, timeSample)

+    if(the.time < 0) { 

+        return(rep(NA, nrow(x$Genotypes)))

+    } 

+    pop <- x$pops.by.time[the.time, -1]

-    if(timeSample == "last") {

-        the.time <- nrow(tmp$pops.by.time)

-    } else if (timeSample %in% c("uniform", "unif")) {

-        the.time <- sample(which(tmp$PerSampleStats[, 4] > 0), 1)

+    if(all(pop == 0)) {

+        stop("You found a bug: this should never happen")

     }

-    pop <- tmp$pops.by.time[the.time, -1]

-    ##       popSize <- tmp$PerSampleStats[the.time, 1]

-    ## genot <- sample(seq_along(pop), 1, prob = pop)

-    return(tmp$Genotypes[, sample(seq_along(pop), 1, prob = pop)])

+    

+    if(typeSample %in% c("wholeTumor", "whole")) {

+        popSize <- x$PerSampleStats[the.time, 1]

+        return( as.numeric((tcrossprod(pop,

+                                       x$Genotypes)/popSize) > thresholdWhole) )

+    } else if (typeSample %in%  c("singleCell", "single")) {

+

+          return(x$Genotypes[, sample(seq_along(pop), 1, prob = pop)])

+      } else {

+            stop("Unknown typeSample option")

+        }

 }

-get.mut.vector <- function(x, timeSample = "whole", typeSample = "last",

-                           thresholdWhole = 0.5) {

-    if(typeSample %in% c("wholeTumor", "whole")) {

-        get.mut.vector.whole(x, timeSample = timeSample,

-                             threshold = thresholdWhole)

-    } else if(typeSample %in%  c("singleCell", "single")) {

-        get.mut.vector.singlecell(x, timeSample = timeSample)

-    }

-}

+

+

+

+

+

+

+

+

 oncoSimul.internal <- function(poset, ## restrict.table,

@@ -973,13 +968,13 @@ oncoSimul.internal <- function(poset, ## restrict.table,

     if(numGenes > 64)

         stop("Largest possible number of genes is 64")

-    

-    if(initSize_species < 10) {

-        warning("initSize_species too small?")

-    }

-    if(initSize_iter < 100) {

-        warning("initSize_iter too small?")

-    }

+    ## These can never be set by the user

+    ## if(initSize_species < 10) {

+    ##     warning("initSize_species too small?")

+    ## }

+    ## if(initSize_iter < 100) {

+    ##     warning("initSize_iter too small?")

+    ## }

     ## numDrivers <- nrow(restrict.table)

     if(length(unique(restrict.table[, 1])) != numDrivers)

@@ -1067,25 +1062,32 @@ oncoSimul.internal <- function(poset, ## restrict.table,

 }

+eFinalMf <- function(initSize, s, j) {

+    ## Expected final sizes for McF, when K is set to the default.

+    # j is number of drivers

+    ## as it says, with no passengers

+    ## Set B(d) = D(N)

+    K <- initSize/(exp(1) - 1)

+    return(K * (exp( (1 + s)^j) - 1))

+}

-

-

-create.muts.by.time <- function(tmp) { ## tmp is the output from Algorithm5

-    if(tmp$NumClones > 1) {

-        NumMutations <- apply(tmp$Genotypes, 2, sum)

-        muts.by.time <- cbind(tmp$pops.by.time[, c(1), drop = FALSE],

-                              t(apply(tmp$pops.by.time[, -c(1),

-                                                       drop = FALSE], 1,

-                                      function(x) tapply(x,

-                                                         NumMutations, sum))))

-        colnames(muts.by.time)[c(1)] <- "Time"

-    } else {

-        muts.by.time <- tmp$pops.by.time

-    }

-    return(muts.by.time)

-} 

+## We are not using this anymore

+## create.muts.by.time <- function(tmp) { ## tmp is the output from Algorithm5

+##     if(tmp$NumClones > 1) {

+##         NumMutations <- apply(tmp$Genotypes, 2, sum)

+##         muts.by.time <- cbind(tmp$pops.by.time[, c(1), drop = FALSE],

+##                               t(apply(tmp$pops.by.time[, -c(1),

+##                                                        drop = FALSE], 1,

+##                                       function(x) tapply(x,

+##                                                          NumMutations, sum))))

+##         colnames(muts.by.time)[c(1)] <- "Time"

+##     } else {

+##         muts.by.time <- tmp$pops.by.time

+##     }

+##     return(muts.by.time)

+## } 

 create.drivers.by.time <- function(tmp, ndr) {

@@ -1203,8 +1205,8 @@ plotClones <- function(z, ndr = NULL, na.subs = TRUE,

 plotDrivers0 <- function(x,

                          ndr,

-                         timescale = 4,

-                         trim.no.drivers = TRUE,

+                         timescale = 1,

+                         trim.no.drivers = FALSE,

                          addtot = TRUE,

                          addtotlwd = 2,

                          na.subs = TRUE, log = "y", type = "l",

@@ -1213,6 +1215,8 @@ plotDrivers0 <- function(x,

                          ...) {

     ## z <- create.drivers.by.time(x, numDrivers)

     z <- create.drivers.by.time(x, ndr)

+    ## we can now never enter here because trim.no.drivers is always FALSE

+    ## in call.

     if(trim.no.drivers && x$MaxDriversLast) {

         fi <- which(apply(z[, -c(1, 2), drop = FALSE], 1,

                           function(x) sum(x) > 0))[1]

@@ -1222,11 +1226,16 @@ plotDrivers0 <- function(x,

     if(na.subs){

         y[y == 0] <- NA

     }

-    if(timescale != 1) {

-        time <- timescale * z[, 1]

-    } else {

-        time <- z[, 1]

-    }

+

+    ## Likewise, we can never enter here now as timescale fixed at 1. And

+    ## this is silly.

+    ## if(timescale != 1) {

+    ##     time <- timescale * z[, 1]

+    ## } else {

+    ##     time <- z[, 1]

+    ## }

+    time <- timescale * z[, 1]

+    

     if(nrow(y) <= 2) type <- "b"

     matplot(x = time,

             y = y,

@@ -1243,61 +1252,193 @@ plotDrivers0 <- function(x,

 }

-rtNoDep <- function(numdrivers) {

-    ## create a restriction table with no dependencies

-    x <- matrix(nrow = numdrivers, ncol = 3)

-    x[, 1] <- 1:numdrivers

-    x[, 2] <- 0

-    x[, 3] <- -9

-    return(x)

-}

-## simulate from generative model. This might not be fully correct!!!

+## No longer used

+## rtNoDep <- function(numdrivers) {

+##     ## create a restriction table with no dependencies

+##     x <- matrix(nrow = numdrivers, ncol = 3)

+##     x[, 1] <- 1:numdrivers

+##     x[, 2] <- 0

+##     x[, 3] <- -9

+##     return(x)

+## }

-simposet <- function(poset, p) {

-    ## if (length(parent.nodes) != length (child.nodes)){

-    ##     print("An Error Occurred")

-    ## }

-    ##    else {

-    num.genes <- max(poset) - 1 ## as root is not a gene

-    genotype <-t(c(1, rep(NA, num.genes)))

-    colnames(genotype) <- as.character(0:num.genes)

+

+## Simulate from generative model. This is a quick function, and is most

+## likely wrong! Never used for anything.

+

+## simposet <- function(poset, p) {

+##     ## if (length(parent.nodes) != length (child.nodes)){

+##     ##     print("An Error Occurred")

+##     ## }

+##     ##    else {

+##     num.genes <- max(poset) - 1 ## as root is not a gene

+##     genotype <-t(c(1, rep(NA, num.genes)))

+##     colnames(genotype) <- as.character(0:num.genes)

-    poset$runif <- runif(nrow(poset))

-    ## this.relation.prob.OK could be done outside, but having it inside

-    ## the loop would allow to use different thresholds for different

-    ## relationships

-    for (i in (1:nrow(poset))) {

-        child <- poset[i, 2]

-        this.relation.prob.OK <- as.numeric(poset[i, "runif"] > p)

-        the.parent <- genotype[ poset[i, 1] ] ## it's the value of parent in genotype. 

-        if (is.na(genotype[child])){

-            genotype[child] <- this.relation.prob.OK * the.parent  

-        }

-        else

-            genotype[child] <- genotype[child]*(this.relation.prob.OK * the.parent)

-    }

-    ##    }

+##     poset$runif <- runif(nrow(poset))

+##     ## this.relation.prob.OK could be done outside, but having it inside

+##     ## the loop would allow to use different thresholds for different

+##     ## relationships

+##     for (i in (1:nrow(poset))) {

+##         child <- poset[i, 2]

+##         this.relation.prob.OK <- as.numeric(poset[i, "runif"] > p)

+##         the.parent <- genotype[ poset[i, 1] ] ## it's the value of parent in genotype. 

+##         if (is.na(genotype[child])){

+##             genotype[child] <- this.relation.prob.OK * the.parent  

+##         }

+##         else

+##             genotype[child] <- genotype[child]*(this.relation.prob.OK * the.parent)

+##     }

+##     ##    }

-    return(genotype)

-}

-

-

-

-

-eFinalMf <- function(initSize, s, j) {

-    ## Expected final sizes for McF, when K is set to the default.

-    # j is number of drivers

-    ## as it says, with no passengers

-    ## Set B(d) = D(N)

-    K <- initSize/(exp(1) - 1)

-    return(K * (exp( (1 + s)^j) - 1))

-}

+##     return(genotype)

+## }

 ## to plot and adjacency matrix in this context can do

 ## plotPoset(intAdjMatToPoset(adjMat))

 ## where intAdjMatToPoset is from best oncotree code: generate-random-trees.

 ## No! the above is simpler

+

+

+

+

+## get.mut.vector.whole <- function(tmp, timeSample = "last", threshold = 0.5) {

+##     ## Obtain, from  results from a simulation run, the vector

+##     ## of 0/1 corresponding to each gene.

+    

+##     ## threshold is the min. proportion for a mutation to be detected

+##     ## We are doing whole tumor sampling here, as in Sprouffske

+

+##     ## timeSample: do we sample at end, or at a time point, chosen

+##     ## randomly, from all those with at least one driver?

+    

+##     if(timeSample == "last") {

+##         if(tmp$TotalPopSize == 0)

+##             warning(paste("Final population size is 0.",

+##                           "Thus, there is nothing to sample with ",

+##                           "sampling last. You will get NAs"))

+##         return(as.numeric(

+##             (tcrossprod(tmp$pops.by.time[nrow(tmp$pops.by.time), -1],

+##                         tmp$Genotypes)/tmp$TotalPopSize) > threshold))

+##     } else if (timeSample %in% c("uniform", "unif")) {

+##           candidate.time <- which(tmp$PerSampleStats[, 4] > 0)

+          

+##           if (length(candidate.time) == 0) {

+##               warning(paste("There is not a single sampled time",

+##                             "at which there are any mutants.",

+##                             "Thus, no uniform sampling possible.",

+##                             "You will get NAs"))

+##               return(rep(NA, nrow(tmp$Genotypes)))

+##           } else if (length(candidate.time) == 1) {

+##                 the.time <- candidate.time

+##             } else {

+##                   the.time <- sample(candidate.time, 1)

+##               }

+##           pop <- tmp$pops.by.time[the.time, -1]

+##           popSize <- tmp$PerSampleStats[the.time, 1]

+##           ## if(popSize == 0)

+##           ##     warning(paste("Population size at this time is 0.",

+##           ##                   "Thus, there is nothing to sample at this time point.",

+##           ##                   "You will get NAs"))

+##           return( as.numeric((tcrossprod(pop,

+##                                        tmp$Genotypes)/popSize) > threshold) )

+##       }

+## }

+

+

+

+##           the.time <- sample(which(tmp$PerSampleStats[, 4] > 0), 1)

+##           if(length(the.time) == 0) {

+##               warning(paste("There are no clones with drivers at any time point.",

+##                             "No uniform sampling possible.",

+##                             "You will get a vector of NAs."))

+##             return(rep(NA, nrow(tmp$Genotypes)))  

+##           }

+## get.mut.vector.singlecell <- function(tmp, timeSample = "last") {

+##     ## No threshold, as single cell.

+

+##     ## timeSample: do we sample at end, or at a time point, chosen

+##     ## randomly, from all those with at least one driver?

+    

+##     if(timeSample == "last") {

+##         the.time <- nrow(tmp$pops.by.time)

+##     } else if (timeSample %in% c("uniform", "unif")) {

+##          candidate.time <- which(tmp$PerSampleStats[, 4] > 0)

+         

+##          if (length(candidate.time) == 0) {

+##              warning(paste("There is not a single sampled time",

+##                            "at which there are any mutants.",

+##                            "Thus, no uniform sampling possible.",

+##                            "You will get NAs"))

+##              return(rep(NA, nrow(tmp$Genotypes)))

+##          } else if (length(candidate.time) == 1) {

+##                the.time <- candidate.time

+##            } else {

+##                  the.time <- sample(candidate.time, 1)

+##              }

+

+##      }

+##     pop <- tmp$pops.by.time[the.time, -1]

+##     ##       popSize <- tmp$PerSampleStats[the.time, 1]

+##     ## genot <- sample(seq_along(pop), 1, prob = pop)

+##     if(all(pop == 0)) {

+##         warning(paste("All clones have a population size of 0",

+##                       "at the chosen time. Nothing to sample.",

+##                       "You will get NAs"))

+##         return(rep(NA, nrow(tmp$Genotypes)))

+##     } else {

+##           return(tmp$Genotypes[, sample(seq_along(pop), 1, prob = pop)])

+##       }

+## }

+

+

+## get.mut.vector <- function(x, timeSample = "whole", typeSample = "last",

+##                            thresholdWhole = 0.5) {

+##     if(typeSample %in% c("wholeTumor", "whole")) {

+##         get.mut.vector.whole(x, timeSample = timeSample,

+##                              threshold = thresholdWhole)

+##     } else if(typeSample %in%  c("singleCell", "single")) {

+##         get.mut.vector.singlecell(x, timeSample = timeSample)

+##     }

+## }

+

+

+

+

+

+## plotClonePhylog <- function(x, timeEvent = FALSE,

+##                             showEvents = TRUE,

+##                             fixOverlap = TRUE) {

+##     if(!inherits(x, "oncosimul2"))

+##         stop("Phylogenetic information is only stored with v >=2")

+##     if(nrow(x$other$PhylogDF) == 0)

+##         stop("It seems you run the simulation with keepPhylog= FALSE")

+##     ## requireNamespace("igraph")

+##     df <- x$other$PhylogDF

+##     if(!showEvents) {

+##         df <- df[!duplicated(df[, c(1, 2)]), ]

+##     }

+##     g <- igraph::graph.data.frame(df)

+##     l0 <- igraph::layout.reingold.tilford(g)

+##     if(!timeEvent) {

+##         plot(g, layout = l0)

+##     } else {

+##         l1 <- l0

+##         indexAppear <- match(V(g)$name, as.character(df[, 2]))

+##         firstAppear <- df$time[indexAppear]

+##         firstAppear[1] <- 0

+##         l1[, 2] <- (max(firstAppear) - firstAppear)

+##         if(fixOverlap) {

+##             dx <- which(duplicated(l1[, 1]))

+##             if(length(dx)) {

+##                 ra <- range(l1[, 1])

+##                 l1[dx, 1] <- runif(length(dx), ra[1], ra[2])

+##             }

+##         }

+##         plot(g, layout = l1)         

+##     }

+## }

@@ -252,22 +252,24 @@ OTtoPoset <- function(x) {

 ##     return(am[cn, cn])

 ## }

-sortAdjMat <- function(am) {

-    ## If column names, except Root, are integers, sort as integers. O.w.,

-    ## general lexicog. sort.

-    cn <- colnames(am)

-    rootpos <- grep("^Root$", cn) 

-    if(length(rootpos) != 1)

-        stop("No root in adj mat, or multiple Roots")

-    cn0 <- colnames(am)[-rootpos]

-    namesInts <- type.convert(cn0, as.is = TRUE)

-    if(is.integer(namesInts)) {

-        cn <- c("Root", sort(namesInts))

-    } else {

-        cn <- c("Root", sort(cn0))

-    }

-    return(am[cn, cn])

-}

+

+## No longer used

+## sortAdjMat <- function(am) {

+##     ## If column names, except Root, are integers, sort as integers. O.w.,

+##     ## general lexicog. sort.

+##     cn <- colnames(am)

+##     rootpos <- grep("^Root$", cn) 

+##     if(length(rootpos) != 1)

+##         stop("No root in adj mat, or multiple Roots")

+##     cn0 <- colnames(am)[-rootpos]

+##     namesInts <- type.convert(cn0, as.is = TRUE)

+##     if(is.integer(namesInts)) {

+##         cn <- c("Root", sort(namesInts))

+##     } else {

+##         cn <- c("Root", sort(cn0))

+##     }

+##     return(am[cn, cn])

+## }

@@ -98,6 +98,7 @@ gm.to.geneModuleL <- function(gm, one.to.one) {

     if(length(unique(geneMod$Gene)) != nrow(geneMod)) {

         stop("Are there identical gene names in different modules?")

     }

+    ## I think this is unreacheable now. Caught earlier.

     if(length(unique(geneMod$GeneNumID)) != nrow(geneMod)) {

         stop("Are there identical gene names in different modules?")

     }

@@ -193,7 +194,8 @@ to.long.rt <- function(rt, idm) {

         z$childNumID <- idm[z$child]

         z$parentsNumID <- idm[z$parents]

         if( any(is.na(z$parentsNumID)) ||

-           any(is.na(z$childNumID)) ) {

+            any(is.na(z$childNumID)) ) {

+            ## I think this can no longer be reached ever. Caught before.

             stop(paste("An ID is NA:",

                        "Is a gene part of two different modules?",

                        "(That includes being by itself and part",

@@ -531,46 +533,46 @@ allFitnessEffects <- function(rT = NULL,

 }

-

-rtAndGeneModule <- function(mdeps, gM = NULL) {

-    ## To show a table of restrictions when there are modules. Do not use

-    ## for anything else. Maybe as intermediate to plotting.

+## No longer used

+## rtAndGeneModule <- function(mdeps, gM = NULL) {

+##     ## To show a table of restrictions when there are modules. Do not use

+##     ## for anything else. Maybe as intermediate to plotting.

-    ## Specify restriction table of modules and a mapping of modules to

-    ## genes. gM is a named vector; names are modules, values are elements

-    ## of each module.

-

-    ## We do nothing important if gM is NULL except checks

-

-    ## If there are modules, the table shows the individual genes.

-    checkRT(mdeps)

-    ## if(ncol(mdeps) != 5)

-    ##     stop("mdeps must be of exactly 5 columns")

-    ## if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))

-    ##     stop(paste("Column names of mdeps not of appropriate format. ",

-    ##                "Should be parent, child, s, sh, typeDep"))

-    if(!is.null(gM)) {

-        if(any(is.na(match(mdeps[ , 1], names(gM)))))

-            stop("Some values in parent not from a known module")

-        if(any(is.na(match(mdeps[ , 2], names(gM)))))

-            stop("Some values in child not from a known module")

-        if(any(is.na(match(names(gM), c(mdeps[, 1], mdeps[, 2])))))

-            stop("Some values in module in neither parent or child")

+##     ## Specify restriction table of modules and a mapping of modules to

+##     ## genes. gM is a named vector; names are modules, values are elements

+##     ## of each module.

+

+##     ## We do nothing important if gM is NULL except checks

+

+##     ## If there are modules, the table shows the individual genes.

+##     checkRT(mdeps)

+##     ## if(ncol(mdeps) != 5)

+##     ##     stop("mdeps must be of exactly 5 columns")

+##     ## if(!identical(colnames(mdeps), c("parent", "child", "s", "sh", "typeDep")))

+##     ##     stop(paste("Column names of mdeps not of appropriate format. ",

+##     ##                "Should be parent, child, s, sh, typeDep"))

+##     if(!is.null(gM)) {

+##         if(any(is.na(match(mdeps[ , 1], names(gM)))))

+##             stop("Some values in parent not from a known module")

+##         if(any(is.na(match(mdeps[ , 2], names(gM)))))

+##             stop("Some values in child not from a known module")

+##         if(any(is.na(match(names(gM), c(mdeps[, 1], mdeps[, 2])))))

+##             stop("Some values in module in neither parent or child")

-        parent <- gM[mdeps[, 1]]

-        child <- gM[mdeps[, 2]]

-        df <- data.frame(parent = parent,

-                         child = child,

-                         s = mdeps$s,

-                         sh = mdeps$sh,

-                         typeDep = mdeps$typeDep,

-                         stringsAsFactors = FALSE)

-    } else {

-        df <- mdeps

-    }

-    rownames(df) <- seq.int(nrow(df))

-    return(df)

-}

+##         parent <- gM[mdeps[, 1]]

+##         child <- gM[mdeps[, 2]]

+##         df <- data.frame(parent = parent,

+##                          child = child,

+##                          s = mdeps$s,

+##                          sh = mdeps$sh,

+##                          typeDep = mdeps$typeDep,

+##                          stringsAsFactors = FALSE)

+##     } else {

+##         df <- mdeps

+##     }

+##     rownames(df) <- seq.int(nrow(df))

+##     return(df)

+## }

 ## wrap.test.rt <- function(rt, gM = NULL) {

 ##     ## FIXME add epistasis and orderEffects

@@ -579,18 +581,18 @@ rtAndGeneModule <- function(mdeps, gM = NULL) {

 ##     wrap_test_rt(lrt$long.rt, lrt$geneModule)

 ## }

-

-wrap.readFitnessEffects <- function(rt, epi, oe, ni, gm, echo = TRUE) {

-    tt <- allFitnessEffects(rt, epi, oe, ni, gm)

-    readFitnessEffects(tt, echo = echo)

-    ## readFitnessEffects(tt$long.rt,

-    ##                    tt$long.epistasis,

-    ##                    tt$long.orderEffects,

-    ##                    tt$long.geneNoInt,

-    ##                    tt$geneModule,

-    ##                    tt$gMOneToOne,

-    ##                    echo = TRUE)

-}

+## No longer used

+## wrap.readFitnessEffects <- function(rt, epi, oe, ni, gm, echo = TRUE) {

+##     tt <- allFitnessEffects(rt, epi, oe, ni, gm)

+##     readFitnessEffects(tt, echo = echo)

+##     ## readFitnessEffects(tt$long.rt,

+##     ##                    tt$long.epistasis,

+##     ##                    tt$long.orderEffects,

+##     ##                    tt$long.geneNoInt,

+##     ##                    tt$geneModule,

+##     ##                    tt$gMOneToOne,

+##     ##                    echo = TRUE)

+## }

 evalGenotype <- function(genotype, fitnessEffects,

                          verbose = FALSE,

@@ -659,8 +661,8 @@ evalAllGenotypes <- function(fitnessEffects, order = TRUE, max = 256,

     nn <- nrow(fitnessEffects$geneModule) -1  + nrow(fitnessEffects$long.geneNoInt)

     tnn <- tot(nn)

     if(tnn > max) {

-        m <- paste("There are ", tnn, "genotypes.")

-        m <- paste(m, " This is larger than max.")

+        m <- paste("There are", tnn, "genotypes.")

+        m <- paste(m, "This is larger than max.")

         m <- paste(m, "Adjust max and rerun if you want")

         stop(m)

     }

@@ -859,16 +861,22 @@ nr_oncoSimul.internal <- function(rFE,

     } else {

         initMutant <- vector(mode = "integer")

     }

-    if(initSize_species < 10) {

-        warning("initSize_species too small?")

-    }

-    if(initSize_iter < 100) {

-        warning("initSize_iter too small?")

-    }

+    ## these are never user options

+    ## if(initSize_species < 10) {

+    ##     warning("initSize_species too small?")

+    ## }

+    ## if(initSize_iter < 100) {

+    ##     warning("initSize_iter too small?")

+    ## }

     if(typeFitness %in% c("bozic1", "bozic2")) {

         thesh <- unlist(lapply(rFE$long.rt, function(x) x$sh))

-        thes <- unlist(lapply(rFE$long.rt, function(x) x$s))

+        ## thes <- unlist(lapply(rFE$long.rt, function(x) x$s))

+        thes <- unlist(c(lapply(rFE$long.rt, function(x) x$s),

+                         lapply(rFE$long.epistasis, function(x) x$s),

+                         lapply(rFE$long.orderEffects, function(x) x$s),

+                         rFE$long.geneNoInt$s))

+        

         if(any(thes > 1 )) {

             m <- paste("You are using a Bozic model with",

                        "the new restriction specification, and you have",

@@ -876,13 +884,16 @@ nr_oncoSimul.internal <- function(rFE,

                        "But that is the same as setting that to 1:",

                        "we obviously cannot allow negative death rates,",

                        "nor problems derived from multiplying odd or even",

-                       "numbers of negative numbers.")

+                       "numbers of negative numbers.",

+                       "You can set those > 1 to 1, but then you will",

+                       "eventually get the simulations to abort when the",

+                       "death rate becomes 0.")

             stop(m)

         }

-        if(any(thesh == -1)) {

+        if(any( (thesh <= -1) & (thesh > -Inf))) {

             m <- paste("You are using a Bozic model with",

                        "the new restriction specification, and you have",

-                       "at least one sh of -1. Maybe you mean -Inf?")

+                       "at least one sh <= -1. Maybe you mean -Inf?")

             warning(m)

         }

         if(any(thes == 1 )) {

Binary files a/data/examplesFitnessEffects.RData and b/data/examplesFitnessEffects.RData differ

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

+Changes in version 1.99.6 (2015-09-26):

+	- Improved test coverage and removed stringsAsfactors from tests