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

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progresion with Epistasis 

-Version: 2.3.2

-Date: 2016-04-14

+Version: 2.3.3

+Date: 2016-06-23

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

 		     email = "rdiaz02@gmail.com"),

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

@@ -13,21 +13,23 @@ Description: Functions for forward population genetic simulation in

     asexual populations, with special focus on cancer progression.

     Fitness can be an arbitrary function of genetic interactions between

     multiple genes or modules of genes, including epistasis, order

-    restrictions in mutation accumulation, and order effects.  Simulations

+    restrictions in mutation accumulation, and order effects.  Mutation

+    rates can differ between genes, and we can include mutator/antimutator

+    genes (to model mutator phenotypes). Simulations

     use continuous-time models and can include driver and passenger genes

     and modules.  Also included are functions for simulating random DAGs

     of the type found in Oncogenetic Tress, Conjunctive Bayesian Networks,

     and other tumor progression models, and for plotting and sampling from

     single or multiple realizations of the simulations, including

-    single-cell sampling, as well as functions for plotting the true

-    phylogenetic relationships of the clones.

+    single-cell sampling, as well as functions for plotting the parent-child

+    relationships of the clones.

 biocViews: BiologicalQuestion, SomaticMutation

 License: GPL (>= 3)

 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, methods, RColorBrewer, grDevices

-Suggests: BiocStyle, knitr, Oncotree, testthat

+Depends: R (>= 3.3.0)

+Imports: Rcpp (>= 0.12.4), parallel, data.table, graph, Rgraphviz, gtools, igraph, methods, RColorBrewer, grDevices, car, dplyr, smatr, ggplot2, ggrepel

+Suggests: BiocStyle, knitr, Oncotree, testthat (>= 1.0.0)

 LinkingTo: Rcpp

 VignetteBuilder: knitr

@@ -2,8 +2,13 @@ useDynLib(OncoSimulR, .registration=TRUE)

 export("oncoSimulPop", "oncoSimulIndiv", "samplePop",

        "plotPoset", "oncoSimulSample", "allFitnessEffects",

-       "evalGenotype", "evalAllGenotypes", "simOGraph"

-     , "plotClonePhylog", "OncoSimulWide2Long"

+       "evalGenotype", "evalAllGenotypes", "simOGraph",

+       "plotClonePhylog", "OncoSimulWide2Long",

+       "allMutatorEffects", "evalAllGenotypesMut",

+       "evalGenotypeMut", "evalGenotypeFitAndMut",

+       "evalAllGenotypesFitAndMut",

+       "rfitness",

+       "plotFitnessLandscape"

        )

 S3method(plot, oncosimul)

@@ -13,8 +18,14 @@ S3method(plot, oncosimulpop)

 S3method(summary, oncosimulpop)

 S3method(print, oncosimulpop)

 S3method(plot, fitnessEffects)

+S3method(plot, genotype_fitness_matrix)

+S3method(plot, evalAllGenotypes)

+S3method(plot, evalAllGenotypesMut)

+import(ggplot2)

+importFrom("ggrepel", geom_text_repel, geom_label_repel)

+importFrom("stats", "rnorm")

 importFrom("data.table", rbindlist, .rbind.data.table) 

 importFrom(Rcpp, evalCpp)

 importFrom("igraph", igraph.to.graphNEL, graph.data.frame, V, E,

@@ -31,3 +42,11 @@ importFrom("stats", "na.omit", "runif", "smooth.spline")

 importFrom("utils", "type.convert")

 importFrom("RColorBrewer", "brewer.pal")

 importFrom("grDevices", "colorRampPalette", "hsv", "rainbow")

+importFrom("dplyr", "full_join", "left_join", "right_join", "%>%", "mutate",

+           "filter")

+importFrom("smatr", "ma") ## for major axis regression in some tests

+importFrom("car", "linearHypothesis")

+

+

+

+

@@ -23,6 +23,7 @@ oncoSimulSample <- function(Nindiv,

                             model = "Exp",

                             numPassengers = 0,

                             mu = 1e-6,

+                            muEF = NULL,

                             detectionSize = round(runif(Nindiv, 1e5, 1e8)),

                             detectionDrivers = {

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

@@ -62,6 +63,7 @@ oncoSimulSample <- function(Nindiv,

                             thresholdWhole = 0.5,

                             initMutant = NULL,

                             verbosity  = 1,

+                            showProgress = TRUE,

                             seed = "auto"){

     ## No longer using mclapply, because of the way we use the limit on

     ## the number of tries.

@@ -164,6 +166,7 @@ oncoSimulSample <- function(Nindiv,

                                model = model,

                                numPassengers = numPassengers,

                                mu = mu,

+                               muEF = muEF,

                                detectionSize = params[indiv, "detectionSize"],

                                detectionDrivers = params[indiv, "detectionDrivers"],

                                sampleEvery = sampleEvery,

@@ -195,6 +198,14 @@ oncoSimulSample <- function(Nindiv,

         numToRun <- (numToRun - 1)

         attemptsUsed <- attemptsUsed + tmp$other$attemptsUsed

         attemptsLeft <- (max.num.tries.total - attemptsUsed)

+        if(showProgress) {

+            cat("......  Done individual ", indiv,

+                ". Used ", attemptsUsed, "attempts.",

+                ". Running for ", as.double(difftime(Sys.time(),

+                                                    startTime, units = "secs")),

+                " secs.\n"

+                )

+        }

         indiv <- indiv + 1

         ## We need to check in exactly this order. Attempts left only

@@ -243,17 +254,30 @@ oncoSimulSample <- function(Nindiv,

 }

-samplePop <- function(x, timeSample = "last", typeSample = "whole",

+samplePop <- function(x, timeSample = "last",

+                      typeSample = "whole",

                       thresholdWhole = 0.5,

-                      geneNames = NULL) {

+                      geneNames = NULL,

+                      popSizeSample = NULL) {

+    ## timeSample <- match.arg(timeSample)

     gN <- geneNames

+    

+    if(!is.null(popSizeSample) && (length(popSizeSample) > 1) &&

+       (length(popSizeSample) != length(x))) {

+        message("length popSizeSample != number of subjects")

+        popSizeSample <- rep(popSizeSample, length.out = length(x))

+        }

+    ## A hack to prevent Map from crashing with a NULL

+    if(is.null(popSizeSample)) popSizeSample <- -99

     if(inherits(x, "oncosimulpop")) {

         z <- do.call(rbind,

-                     lapply(x,

-                            get.mut.vector,

-                            timeSample = timeSample,

-                            typeSample = typeSample,

-                            thresholdWhole = thresholdWhole))

+                     Map(get.mut.vector,

+                         x = x,

+                         timeSample = timeSample,

+                         typeSample = typeSample,

+                         thresholdWhole = thresholdWhole,

+                         popSizeSample = popSizeSample

+                         ))

         ## We need to check if the object is coming from v.2., to avoid

         ## having to force passing a vector of names

         if(is.null(gN) && (!is.null(x[[1]]$geneNames)))

@@ -262,7 +286,8 @@ samplePop <- function(x, timeSample = "last", typeSample = "whole",

         z <- get.mut.vector(x,

                             timeSample = timeSample,

                             typeSample = typeSample,

-                            thresholdWhole = thresholdWhole)

+                            thresholdWhole = thresholdWhole,

+                            popSizeSample = popSizeSample)

         dim(z) <- c(1, length(z))

         if(is.null(gN) && (!is.null(x$geneNames)))

             gN <- geneNames

@@ -285,6 +310,7 @@ oncoSimulPop <- function(Nindiv,

                          model = "Exp",

                          numPassengers = 0,

                          mu = 1e-6,

+                         muEF = NULL,

                          detectionSize = 1e8,

                          detectionDrivers = 4,

                          sampleEvery = ifelse(model %in% c("Bozic", "Exp"), 1,

@@ -329,6 +355,7 @@ oncoSimulPop <- function(Nindiv,

                         model = model,

                         numPassengers = numPassengers,

                         mu = mu,

+                        muEF = muEF,

                         detectionSize = detectionSize,

                         detectionDrivers = detectionDrivers,

                         sampleEvery = sampleEvery,

@@ -365,6 +392,7 @@ oncoSimulIndiv <- function(fp,

                            model = "Exp",

                            numPassengers = 0,

                            mu = 1e-6,

+                           muEF = NULL,

                            detectionSize = 1e8,

                            detectionDrivers = 4,

                            sampleEvery = ifelse(model %in% c("Bozic", "Exp"), 1,

@@ -434,8 +462,11 @@ oncoSimulIndiv <- function(fp,

             stop("Unknown model")

     }

+    if( (length(mu) > 1) && !inherits(fp, "fitnessEffects"))

+        stop("Per-gene mutation rates cannot be used with the old poset format")

+

     if(any(mu < 0)) {

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

+        stop("(at least one) mutation rate (mu) is negative")

     }

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

     ## legitimate case?

@@ -463,6 +494,8 @@ oncoSimulIndiv <- function(fp,

             stop(m)

         }

+        if(!is.null(muEF))

+            stop("Mutator effects cannot be especified with the old poset format")

         if(length(initMutant) > 1)

             stop("With the old poset, initMutant can only take a single value.")

         ## Seeding C++ is now much better in new version

@@ -507,7 +540,7 @@ oncoSimulIndiv <- function(fp,

                                      max.wall.time = max.wall.time,

                                      max.num.tries = max.num.tries,

                                      keepEvery = keepEvery,  

-                                     alpha = 0.0015,  

+                                     ## alpha = 0.0015,  

                                      sh = sh,

                                      K = K, 

                                      minDetectDrvCloneSz = minDetectDrvCloneSz,

@@ -558,7 +591,7 @@ oncoSimulIndiv <- function(fp,

                                         max.wall.time = max.wall.time,

                                         max.num.tries = max.num.tries,

                                         keepEvery = keepEvery,  

-                                        alpha = 0.0015,  

+                                        ## alpha = 0.0015,  

                                         K = K, 

                                         minDetectDrvCloneSz = minDetectDrvCloneSz,

                                         extraTime = extraTime,

@@ -566,7 +599,8 @@ oncoSimulIndiv <- function(fp,

                                         onlyCancer = onlyCancer,

                                         errorHitWallTime = errorHitWallTime,

                                         errorHitMaxTries = errorHitMaxTries,

-                                        keepPhylog = keepPhylog),

+                                        keepPhylog = keepPhylog,

+                                        MMUEF = muEF),

                   silent = !verbosity)

         objClass <- c("oncosimul", "oncosimul2")

     }

@@ -648,7 +682,7 @@ plot.oncosimulpop <- function(x, ask = TRUE,

                               show = "drivers", 

                               type = ifelse(show == "genotypes",

                                             "stacked", "line"),

-                              col = c(8, "orange", 6:1),

+                              col = "auto",

                               log = ifelse(type == "line", "y", ""),

                               ltyClone = 2:6,

                               lwdClone = 0.9,

@@ -1389,17 +1423,32 @@ plotClonePhylog <- function(x, N = 1, t = "last",

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

+closest_time <- function(x, size) {

+    ## Find the first time when a given size is reached

+    sizes <- rowSums(x$pops.by.time[, -1, drop = FALSE])

+    candidates <- which(sizes >= size)

+    if(length(candidates) == 0) {

+        warning(paste("Pop size never >= requested size.",

+                      "Thus, there is nothing to sample. You will get NAs"))

+        return(-99)

+    } else {

+        return(candidates[1])

+    }

+}

+

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

-    if(timeSample == "last") {

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

+    if( !is.null(popSizeSample) && (popSizeSample >= 0) )  {

+        the.time <- closest_time(tmp, popSizeSample)

+    } else 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"))

             the.time <- -99

         } else {

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

-          }

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

+        }

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

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

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

             } else {

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

               }

-      } else {

+    } else {

             stop("Unknown timeSample option")

         }

     return(the.time)

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

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

-                           thresholdWhole) {

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

+                           thresholdWhole, popSizeSample) {

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

     if(the.time < 0) { 

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

     } 

@@ -1441,10 +1490,10 @@ get.mut.vector <- function(x, timeSample, typeSample,

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

-        }

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

+    } else {

+        stop("Unknown typeSample option")

+    }

 }

@@ -1507,9 +1556,15 @@ oncoSimul.internal <- function(poset, ## restrict.table,

         restrict.table <- poset.to.restrictTable(poset)

     numDrivers <- nrow(restrict.table)

     numGenes <- (numDrivers + numPassengers)

-    

+    if(numGenes < 2)

+        stop("There must be at least two genes (loci) in the fitness effects.",

+             "If you only care about a degenerate case with just one,",

+             "you can enter a second gene",

+             "with fitness effect of zero.")

     if(numGenes > 64)

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

+        stop("Largest possible number of genes (loci) is 64 for version 1.",

+             "You are strongly encouraged to use the new specification",

+             "as in version 2.")

     ## These can never be set by the user

     ## if(initSize_species < 10) {

@@ -1567,8 +1622,8 @@ oncoSimul.internal <- function(poset, ## restrict.table,

         BNB_Algo5(restrictTable = rtC,

         numDrivers = numDrivers,

         numGenes = numGenes,

-        typeCBN_ = typeCBN,

-        birthRate = birth,

+        typeCBN_= typeCBN,

+        ## birthRate = birth, 

         s = s, 

         death = death,

         mu = mu,

@@ -1582,16 +1637,15 @@ oncoSimul.internal <- function(poset, ## restrict.table,

         verbosity = verbosity,

         speciesFS = speciesFS,

         ratioForce = ratioForce,

-        typeFitness_= typeFitness,

+        typeFitness_ = typeFitness,

         maxram = max.memory,

         mutationPropGrowth = as.integer(mutationPropGrowth),

         initMutant = initMutant,

         maxWallTime = max.wall.time,

         keepEvery = keepEvery,

-        alpha = alpha,

+        # alpha = alpha,

         sh = sh,

         K = K,

-        # end = TimeEvery# , 

         detectionDrivers = detectionDrivers,

         onlyCancer = onlyCancer,

         errorHitWallTime = errorHitWallTime,

@@ -1600,7 +1654,7 @@ oncoSimul.internal <- function(poset, ## restrict.table,

         minDetectDrvCloneSz = minDetectDrvCloneSz,

         extraTime = extraTime

     ),

-             NumDrivers = numDrivers

+    NumDrivers = numDrivers

              ))

 }

@@ -1,19 +1,24 @@

 # This file was generated by Rcpp::compileAttributes

 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

-nr_BNB_Algo5 <- function(rFE, mu, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, alpha, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog) {

-    .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, alpha, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog)

+nr_BNB_Algo5 <- function(rFE, mu, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_) {

+    .Call('OncoSimulR_nr_BNB_Algo5', PACKAGE = 'OncoSimulR', rFE, mu, death, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant_, maxWallTime, keepEvery, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime, keepPhylog, MMUEF, full2mutator_)

 }

-BNB_Algo5 <- function(restrictTable, numDrivers, numGenes, typeCBN_, birthRate, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, alpha, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime) {

-    .Call('OncoSimulR_BNB_Algo5', PACKAGE = 'OncoSimulR', restrictTable, numDrivers, numGenes, typeCBN_, birthRate, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, alpha, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime)

+BNB_Algo5 <- function(restrictTable, numDrivers, numGenes, typeCBN_, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime) {

+    .Call('OncoSimulR_BNB_Algo5', PACKAGE = 'OncoSimulR', restrictTable, numDrivers, numGenes, typeCBN_, s, death, mu, initSize, sampleEvery, detectionSize, finalTime, initSp, initIt, seed, verbosity, speciesFS, ratioForce, typeFitness_, maxram, mutationPropGrowth, initMutant, maxWallTime, keepEvery, sh, K, detectionDrivers, onlyCancer, errorHitWallTime, maxNumTries, errorHitMaxTries, minDetectDrvCloneSz, extraTime)

 }

-## readFitnessEffects <- function(rFE, echo) {

-##     invisible(.Call('OncoSimulR_readFitnessEffects', PACKAGE = 'OncoSimulR', rFE, echo))

-## }

+readFitnessEffects <- function(rFE, echo) {

+     invisible(.Call('OncoSimulR_readFitnessEffects', PACKAGE = 'OncoSimulR', rFE, echo))

+}

+

+evalRGenotype <- function(rG, rFE, verbose, prodNeg, calledBy_) {

+    .Call('OncoSimulR_evalRGenotype', PACKAGE = 'OncoSimulR', rG, rFE, verbose, prodNeg, calledBy_)

+}

-evalRGenotype <- function(rG, rFE, verbose, prodNeg) {

-    .Call('OncoSimulR_evalRGenotype', PACKAGE = 'OncoSimulR', rG, rFE, verbose, prodNeg)

+evalRGenotypeAndMut <- function(rG, rFE, muEF, full2mutator_, verbose, prodNeg) {

+    .Call('OncoSimulR_evalRGenotypeAndMut', PACKAGE = 'OncoSimulR', rG, rFE, muEF, full2mutator_, verbose, prodNeg)

 }

+

new file mode 100644

@@ -0,0 +1,398 @@

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

+## the Free Software Foundation, either version 3 of the License, or

+## (at your option) any later version.

+

+## This program is distributed in the hope that it will be useful,

+## but WITHOUT ANY WARRANTY; without even the implied warranty of

+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

+## GNU General Public License for more details.

+

+## You should have received a copy of the GNU General Public License

+## along with this program.  If not, see <http://www.gnu.org/licenses/>.

+

+## plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

+##     plot.genotype_fitness_matrix <- plotFitnessLandscape 

+

+## FIXME: show only accessible paths? 

+

+plotFitnessLandscape <- function(x, show_labels = TRUE,

+                                 col = c("green4", "red", "yellow"),

+                                 lty = c(1, 2, 3), use_ggrepel = FALSE,

+                                 log = FALSE, max_num_genotypes = 2000,

+                                 ...) {

+

+    ## FIXME future:

+    

+   

+

+    ## - Allow passing order effects. Change "allGenotypes_to_matrix"

+    ##   below? Probably not, as we cannot put order effects as a

+    ##   matrix. Do something else, like allow only order effects if from

+    ##   and allFitnessEffects object.

+

+    ## - Allow selection only some genotypes or alleles

+

+    ## Allow selecting only paths that involve a particular genotype in

+    ## adjacency matrix of genotype i go at row i and column i.  Follow back

+    ## all entries in row i, and their previous, and forward all column i.

+

+

+    ## gfm: genotype fitness matrix

+    ## afe: all fitness effects

+

+    ## We seem to go back and forth, but we need to ensure the afe and the

+    ## gfm are coherent. Since there are many ways to make them fail

+    ## (e.g., a user passing the wrong order, or incomplete matrices, etc)

+    ## we do it the following way.

+

+    ## Yes, we set WT to 1 as we call from_genotype_fitness on

+    ## matrices and genotype_fitness_matrix objects.

+    ## We do that because we call evalAllGenotypes to

+    ## get the string representation, etc. And this is for use

+    ## with OncoSimul.

+

+    ## FIXME: passing fitness as a 2 column data frame

+    

+    if( (inherits(x, "genotype_fitness_matrix")) ||

+        ( (is.matrix(x) || is.data.frame(x)) && (ncol(x) > 2) ) ) {

+        ## Why this? We go back and forth twice. We need both things. We

+        ## could construct the afe below by appropriately pasting the

+        ## columns names

+        afe <- evalAllGenotypes(allFitnessEffects(

+            epistasis = from_genotype_fitness(x)),

+            order = FALSE, addwt = TRUE, max = max_num_genotypes)

+        ## Might not be needed with the proper gfm object (so gmf <- x)

+        ## but is needed if arbitrary matrices.

+        gfm <- allGenotypes_to_matrix(afe) 

+    } else if(inherits(x, "fitnessEffects")) {

+        if(!is.null(x$orderEffects) )

+            stop("We cannot yet deal with order effects")

+        afe <- evalAllGenotypes(x,

+                                order = FALSE,

+                                addwt = TRUE, max = max_num_genotypes)

+        gfm <- allGenotypes_to_matrix(afe)

+    } else if( (inherits(x, "evalAllGenotypes")) ||

+               (inherits(x, "evalAllGenotypesMut"))) {

+        if(any(grepl(">", x[, 1], fixed = TRUE)))

+            stop("We cannot deal with order effects yet.")

+        x <- x[, c(1, 2)]

+        if(x[1, "Genotype"] != "WT") {

+            ## Yes, because we expect this present below

+            x <- rbind(data.frame(Genotype = "WT",

+                                  Fitness = 1,

+                                  stringsAsFactors = FALSE),

+                       x)

+        }

+        afe <- x

+        ## in case we pass an evalAllgenotypesfitandmut

+        gfm <- allGenotypes_to_matrix(afe)

+    } else if(is.data.frame(x)) {

+        ## Assume a two-column data frame of genotypes as character

+        ## vectors and fitness

+        if(colnames(x)[2] != "Fitness")

+            stop("We cannot guess what you are passing here") 

+        afe <- evalAllGenotypes(allFitnessEffects(genotFitness = x),

+                                order = FALSE, addwt = TRUE,

+                                max = max_num_genotypes)

+        gfm <- allGenotypes_to_matrix(afe)

+    } else {

+       stop("We cannot guess what you are passing here") 

+    }

+   

+

+    

+    ## if(inherits(x, "genotype_fitness_matrix")) {

+    ##     ## Why this? We go back and forth twice. We need both things. We

+    ##     ## could construct the afe below by appropriately pasting the

+    ##     ## columns names

+    ##     afe <- evalAllGenotypes(allFitnessEffects(

+    ##         epistasis = from_genotype_fitness(x)),

+    ##         order = FALSE, addwt = TRUE, max = 2000)

+    ##     gfm <- allGenotypes_to_matrix(afe) ## should not be needed? gfm <- x

+    ## } else if(inherits(x, "fitnessEffects")) {

+    ##     if(!is.null(x$orderEffects) )

+    ##         stop("We cannot yet deal with order effects")

+    ##     afe <- evalAllGenotypes(x,

+    ##                             order = FALSE,

+    ##                             addwt = TRUE, max = 2000)

+    ##     gfm <- allGenotypes_to_matrix(x)

+    ## } else {

+    ##     lc <- ncol(x)

+    ##     ## detect an appropriately formatted matrix

+    ##     if( (is.data.frame(x) || is.matrix(x)) ) {

+    ##         if(ncol(x) > 2) {

+    ##             ## We cannot be sure all genotypes are present

+    ##             ## but that is not our business?

+    ##             ## colnames(x)[lc] <- "Fitness"

+    ##             ## So this must be a matrix

+    ##             afe <- evalAllGenotypes(allFitnessEffects(

+    ##                 epistasis = from_genotype_fitness(x)),

+    ##                 order = FALSE, addwt = TRUE, max = 2000)

+    ##             gfm <- allGenotypes_to_matrix(afe)

+    ##         } else {

+    ##             ## This must be a data frame

+    ##             if(!is.data.frame(x))

+    ##                 stop("How are you passing a matrix here?")

+    ##             if(colnames(x)[lc] != "Fitness")

+    ##                 stop("We cannot guess what you are passing here")

+    ##             ## We are passed an allFitnessEffects output

+    ##             ## Will be simpler later as we will know immediately

+    ##             if(x[1, "Genotype"] != "WT") {

+    ##                 ## Yes, because we expect this present below

+    ##                 x <- rbind(data.frame(Genotype = "WT",

+    ##                                       Fitness = 1,

+    ##                                       stringsAsFactors = FALSE),

+    ##                            x)

+    ##             }

+    ##             x <- afe

+    ##             gfm <- allGenotypes_to_matrix(afe)

+    ##         }

+    ##     } else {

+    ##         stop("This is not allowed input")

+    ##     }

+    ## }

+

+    mutated <- rowSums(gfm[, -ncol(gfm)])

+    gaj <- genot_to_adj_mat(gfm)

+    vv <- which(!is.na(gaj), arr.ind = TRUE)

+    

+    ## plot(x = mutated, y = e1$Fitness, ylab = "Fitness",

+    ##      xlab = "", type = "n", axes = FALSE)

+    ## box()

+    ## axis(2)

+    ## text(x = mutated, y = x$Fitness, labels = x$Genotype)

+

+    ## The R CMD CHEKC notes about no visible binding for global variable

+

+    x_from <- y_from <- x_to <- y_to <- Change <- muts <-

+        label <- fitness <- Type <- NULL

+                

+                

+    dd <- data.frame(muts = mutated,

+                     fitness = afe$Fitness,

+                     label = afe$Genotype)

+    cl <- gaj[vv]

+    sg <- data.frame(x_from = mutated[vv[, 1]],

+                     y_from = afe$Fitness[vv[, 1]],

+                     x_to = mutated[vv[, 2]],

+                     y_to = afe$Fitness[vv[, 2]],

+                     Change = factor(ifelse(cl == 0, "Neutral",

+                                     ifelse(cl > 0, "Gain", "Loss"))))

+    ## From http://stackoverflow.com/a/17257422

+    number_ticks <- function(n) {function(limits) pretty(limits, n)}

+    

+    p <- ggplot() +

+        xlab("") + ylab("Fitness") + 

+        theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),

+              panel.grid.minor.x = element_blank()) +

+        geom_segment(data = sg,

+                       aes(x = x_from, y = y_from,

+                           xend = x_to, yend = y_to,

+                           colour = Change,

+                           lty = Change)) + scale_color_manual("Change",

+                                                               values = col)

+    if(log) {

+        p <- p + scale_y_log10(breaks = number_ticks(5))

+    } else {

+        p <- p + scale_y_continuous(breaks = number_ticks(5))

+    }

+    peaks_valleys <- peak_valley(gaj)

+    maxF <- peaks_valleys$peak

+    minF <- peaks_valleys$valley

+

+    ddMM <- dd

+    ddMM$Type <- NA

+    ddMM$Type[minF] <- "Sink"

+    ddMM$Type[maxF] <- "Peak"

+    ddMM <- ddMM[ddMM$Type %in% c("Sink", "Peak"), ]

+    

+    if(show_labels && use_ggrepel) {

+        p <- p + geom_text_repel(data = dd[-c(minF, maxF), ],

+                                 aes(x = muts, y = fitness, label = label)) +

+            geom_label_repel(data = ddMM,

+                             aes(x = muts, y = fitness, label = label, fill = Type),

+                             color = "black",

+                             fontface = "bold")

+            ## geom_label_repel(data = dd[maxF, ], aes(x = muts, y = fitness, label = label),

+            ##                  color = "black",

+            ##                  fontface = "bold", fill = col[1]) +

+            ## geom_label_repel(data = dd[minF, ], aes(x = muts, y = fitness, label = label),

+            ##                  color = "black",

+            ##                  fontface = "bold", fill = col[2]) 

+        

+        }

+    else {

+        if(!show_labels) {

+            ## Use same code, but make label empty

+            ddMM$label <- ""

+        }

+            

+        p <- p + geom_text(data = dd[-c(minF, maxF), ],

+                           aes(x = muts, y = fitness, label = label),

+                           vjust = -0.2, hjust = "inward") +

+            geom_label(data = ddMM,

+                       aes(x = muts, y = fitness, label = label, fill = Type),

+                       vjust = -0.2, hjust = "inward", color = "black",

+                       fontface = "bold")

+            ## geom_label(data = dd[maxF, ], aes(x = muts, y = fitness, label = label),

+            ##            vjust = -0.2, hjust = "inward", color = "black",

+            ##            fontface = "bold", fill = col[1]) +

+            ## geom_label(data = dd[minF, ], aes(x = muts, y = fitness, label = label),

+            ##            vjust = -0.2, hjust = "inward", color = "black",

+            ##            fontface = "bold", fill = col[2])

+    }

+    p <- p + scale_fill_manual("Local\nmax/min",  values = col)

+    p

+}

+

+

+plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

+    plot.genotype_fitness_matrix <-

+        plot_fitness_landscape <-

+            plotFitnessLandscape

+

+

+######################################################################

+######################################################################

+######################################################################

+####

+####            Internal functions

+####

+######################################################################

+######################################################################

+######################################################################

+

+

+

+generate_matrix_genotypes <- function(g) {

+    ## FIXME future: do this for order too? Only if rfitness for order.

+    ## Given a number of genes, generate all possible genotypes.

+    

+    if(g > 20) stop("This would generate more than one million genotypes")

+    ## g: number of genes

+    f1 <- function(n) {

+        lapply(seq.int(n), function(x) combinations(n = n, r = x))

+    }

+    genotNums <- f1(g)

+    list.of.vectors <- function(y) {

+        ## there's got to be a simpler way

+        lapply(unlist(lapply(y, function(x) {apply(x, 1, list)}),

+                      recursive = FALSE),

+               function(m) m[[1]])

+    }

+    genotNums <- list.of.vectors(genotNums)

+    

+    v <- rep(0, g)

+    mat <- matrix(unlist(lapply(genotNums,

+                                function(x) {

+                                    v[x] <- 1

+                                    return(v)

+                                }

+                                )), ncol = g, byrow = TRUE)

+    mat <- rbind(rep(0, g), mat)

+    colnames(mat) <- LETTERS[1:g]

+    return(mat)

+}

+

+

+

+genot_to_adj_mat <- function(x) {

+    ## FIXME this can take about 23% of the time of the ggplot call.

+    ## But them, we are quickly constructing a 2000*2000 matrix

+    ## Given a genotype matrix, as given by allGenotypes_to_matrix, produce a

+    ## directed adjacency matrix between genotypes connected (i.e., those

+    ## that differ by gaining one mutation) with value being the

+    ## difference in fitness between destination and origin

+

+    ## Make sure sorted, so ancestors always before descendants

+    rs0 <- rowSums(x[, -ncol(x)])

+    x <- x[order(rs0), ]

+    rm(rs0)

+    

+    y <- x[, -ncol(x)]

+    f <- x[, ncol(x)]

+    rs <- rowSums(y)

+

+    ## Move this to C++?

+    adm <- matrix(NA, nrow = length(rs), ncol = length(rs))

+    for(i in 1:length(rs)) { ## i is the current genotype

+        candidates <- which(rs == (rs[i] + 1))

+        for(j in candidates) {

+            sumdiff <- sum(abs(y[j, ] - y[i, ]))

+            ## if(sumdiff < 0) stop("eh?")

+            if(sumdiff == 1) adm[i, j] <- (f[j] - f[i])

+        }

+    }

+    return(adm)

+}

+

+

+peak_valley <- function(x) {

+    ## FIXME: when there are no identical entries, all this

+    ## can be simplified a lot. Actually, there could be a way

+    ## to only use the slow code on paths with 0.

+    ## But this does not seem to be the CPU hog.

+    ## x is an adjacency matrix where i,j is fitness_j - fitness_i.  Thus,

+    ## negative values means j has less fitness than the incoming.  And if

+    ## all not NA entries of a column are negative, then column j is is

+    ## sink candidate

+    ## Thus we locate local maxima and minima

+

+    ## Some of this complicated as we need to detect paths like -3, 0, 0,

+    ## 3.  The -3 and the two 0 are local minima with identical values.

+

+    ## We assume crucially that ancestors always have smaller indexes in

+    ## the matrix than descendants.

+

+    ## Valleys

+    bad_back <- vector("integer", nrow(x))

+    for(i in ncol(x):1) {

+        if(any(x[i, ] < 0, na.rm = TRUE) || bad_back[i]) {

+            ## this node is bad. Any ancestor with fitness >= is bad

+            bad_back[i] <- 1

+            reach_b <- which(x[, i] <= 0)

+            bad_back[reach_b] <- 1

+        }

+    }

+    bad_fwd <- vector("integer", nrow(x))

+    for(i in 1:nrow(x)) {

+        if( any(x[, i] > 0, na.rm = TRUE) || bad_fwd[i] ) {

+            ## this node is bad. Any descendant with fitness >= is bad

+            bad_fwd[i] <- 1

+            reach_f <- which(x[i, ] >= 0)

+            bad_fwd[reach_f] <- 1

+        }

+    }

+    bad <- union(which(bad_back == 1), which(bad_fwd == 1))

+    candidate <- which(apply(x, 2, function(z) all(z <= 0, na.rm = TRUE)))

+    valley <- setdiff(candidate, bad)

+

+    null <- suppressWarnings(try({rm(bad_back, bad_fwd, reach_b, reach_f, candidate)},

+        silent = TRUE))

+    ## Peaks

+    bad_back <- vector("integer", nrow(x))

+    for(i in ncol(x):1) {

+        if(any(x[i, ] > 0, na.rm = TRUE) || bad_back[i]) {

+            ## this node is bad. Any ancestor with fitness >= is bad

+            bad_back[i] <- 1

+            reach_b <- which(x[, i] >= 0)

+            bad_back[reach_b] <- 1

+        }

+    }

+    bad_fwd <- vector("integer", nrow(x))

+    for(i in 1:nrow(x)) {

+        if( any(x[, i] < 0, na.rm = TRUE) || bad_fwd[i] ) {

+            ## this node is bad. Any descendant with fitness >= is bad

+            bad_fwd[i] <- 1

+            reach_f <- which(x[i, ] <= 0)

+            bad_fwd[reach_f] <- 1

+        }

+    }

+    bad <- union(which(bad_back == 1), which(bad_fwd == 1))

+    candidate <- which(apply(x, 2, function(z) all(z >= 0, na.rm = TRUE)))

+    peak <- setdiff(candidate, bad)

+    return(list(peak = peak, valley = valley))

+}

new file mode 100644

@@ -0,0 +1,180 @@

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

+## the Free Software Foundation, either version 3 of the License, or

+## (at your option) any later version.

+

+## This program is distributed in the hope that it will be useful,

+## but WITHOUT ANY WARRANTY; without even the implied warranty of

+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

+## GNU General Public License for more details.

+

+## You should have received a copy of the GNU General Public License

+## along with this program.  If not, see <http://www.gnu.org/licenses/>.

+

+

+

+## Functions that allow passing a matrix or data frame of mappings

+## genotype -> fitness so this is taken as input in fitnessEffects.

+

+

+from_genotype_fitness <- function(x) {

+    ## Would break with output from allFitnessEffects and

+    ## output from allGenotypeAndMut

+    if(ncol(x) > 2) {

+        ## Of course, can ONLY work with epistastis, NOT order

+        return(genot_fitness_to_epistasis(x))

+    } else {

+        ## Make sure no factors

+        if(is.factor(x[, 1])) x[, 1] <- as.character(x[, 1])

+        omarker <- any(grepl(">", x[, 1], fixed = TRUE))

+        emarker <- any(grepl(",", x[, 1], fixed = TRUE))

+        nogoodepi <- any(grepl(":", x[, 1], fixed = TRUE))

+        ## if(omarker && emarker) stop("Specify only epistasis or order, not both.")

+        if(nogoodepi && emarker) stop("Specify the genotypes separated by a ',', not ':'.")

+        if(nogoodepi && !emarker) stop("Specify the genotypes separated by a ',', not ':'.")

+        ## if(nogoodepi && omarker) stop("If you want order, use '>' and if epistasis ','.")

+        ## if(!omarker && !emarker) stop("You specified neither epistasis nor order")

+        if(omarker) {

+            ## do something. To be completed