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

 Package: OncoSimulR

 Type: Package

 Title: Forward Genetic Simulation of Cancer Progresion with Epistasis 

-Version: 2.3.8

-Date: 2016-07-08

+Version: 2.3.9

+Date: 2016-07-09

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

 		     email = "rdiaz02@gmail.com"),

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

@@ -17,12 +17,13 @@ Description: Functions for forward population genetic simulation in

     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

+    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

+    and other tumor progression models; plotting and sampling from

     single or multiple realizations of the simulations, including

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

-    relationships of the clones.

+    single-cell sampling; plotting the parent-child relationships of the

+    clones; generating random fitness landscapes (Rough Mount Fuji, House

+    of Cards, and additive models) and plotting them.

 biocViews: BiologicalQuestion, SomaticMutation

 License: GPL (>= 3)

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

@@ -22,6 +22,8 @@ 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,

+                                 only_accessible = FALSE,

+                                 accessible_th = 0,

                                  ...) {

     ## FIXME future:

@@ -37,7 +39,6 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

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

@@ -52,11 +53,17 @@ plotFitnessLandscape <- function(x, show_labels = TRUE,

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

     ## with OncoSimul.

+  

     tfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)

-

     mutated <- rowSums(tfm$gfm[, -ncol(tfm$gfm)])

     gaj <- genot_to_adj_mat(tfm$gfm)

+    if(only_accessible) {

+        gaj <- filter_inaccessible(gaj, accessible_th)

+        remaining <- as.numeric(colnames(gaj))

+        mutated <- mutated[remaining]

+        tfm$afe <- tfm$afe[remaining, , drop = FALSE]

+    }

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

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

@@ -171,6 +178,44 @@ plot.evalAllGenotypes <- plot.evalAllGenotypesMut <-

 ######################################################################

+## wrap_filter_inaccessible <- function(x, max_num_genotypes, accessible_th) {

+##     ## wrap it, for my consumption

+##     tfm <- to_Fitness_Matrix(x, max_num_genotypes = max_num_genotypes)

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

+##     gaj <- genot_to_adj_mat(tfm$gfm)

+##     gaj <- filter_inaccessible(gaj, accessible_th)

+##     remaining <- as.numeric(colnames(gaj))

+##     mutated <- mutated[remaining]

+##     tfm$afe <- tfm$afe[remaining, , drop = FALSE]

+##     return(list(remaining = remaining,

+##                 mutated = mutated,

+##                 tfm = tfm))

+## }

+

+count_accessible_g <- function(gfm, accessible_th) {

+    gaj <- genot_to_adj_mat(gfm)

+    gaj <- filter_inaccessible(gaj, accessible_th)

+    return(ncol(gaj) - 1)

+}

+

+filter_inaccessible <- function(x, accessible_th) {

+    ## Return an adjacency matrix with only accessible paths. x is the gaj

+    ## matrix created in the plots. The difference between genotypes

+    ## separated by a hamming distance of 1

+    colnames(x) <- rownames(x) <- 1:ncol(x)

+    while(TRUE) {

+        ## remove first column

+        ## We use fact that all(na.omit(x) < u) is true if all are NA

+        ## so inaccessible rows are removed and thus destination columns

+        wrm <- which(apply(x[, -1, drop = FALSE], 2,

+                           function(y) {all(na.omit(y) < accessible_th)})) + 1

+        if(length(wrm) == 0) break;

+        x <- x[-wrm, -wrm, drop = FALSE]

+    }

+    x[x < 0] <- NA

+    return(x)

+}

+

 generate_matrix_genotypes <- function(g) {

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

@@ -6,10 +6,11 @@ rfitness <- function(g, c= 0.5,

                                      ## reference. If "max" this is rep(1, g)

                      scale = NULL, ## a two-element vector: min and max

                      wt_is_1 = TRUE, ## wt has fitness 1

-                     log = FALSE ## log only makes sense if all values >

+                     log = FALSE, ## log only makes sense if all values >

                                  ## 0. scale with min > 0, and/or set

                                  ## wt_is_1 = TRUE

-                     ) {

+                     min_accessible_genotypes = 0,

+                     accessible_th = 0) {

     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene

     ## and Crona, 2014. And this allows moving from HoC to purely additive

     ## changing c and sd.

@@ -17,42 +18,57 @@ rfitness <- function(g, c= 0.5,

     ## FIXME future: do this with order too?

     ##    - do not generate a matrix of genotypes but a matrix of ordered genot.

     m <- generate_matrix_genotypes(g)

-    f_r <- rnorm(nrow(m), mean = 0, sd = sd)

-    if(inherits(reference, "character") && length(reference) == 1) {

-        if(reference == "random") {

-            reference <- m[sample(nrow(m), 1), ]

-        } else if(reference == "max") {

-            reference <- rep(1, g)

+    done <- FALSE

+    while(!done) {

+        f_r <- rnorm(nrow(m), mean = 0, sd = sd)

+        if(inherits(reference, "character") && length(reference) == 1) {

+            if(reference == "random") {

+                referenceI <- m[sample(nrow(m), 1), ]

+            } else if(reference == "max") {

+                referenceI <- rep(1, g)

+            } 

+        } else {

+            referenceI <- reference

+            }

+        d_reference <- apply(m, 1, function(x) sum(abs(x - referenceI)))

+        f_det <- -c * d_reference

+        ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona

+        fi <- f_r + f_det

+        

+        if(!is.null(scale)) {

+            fi <- (fi - min(fi))/(max(fi) - min(fi))

+            fi <- scale[1] + fi * (scale[2] - scale[1])

         }

-    }

-    d_reference <- apply(m, 1, function(x) sum(abs(x - reference)))

-    f_det <- -c * d_reference

-    ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona

-    fi <- f_r + f_det

-    

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

-        fi <- (fi - min(fi))/(max(fi) - min(fi))

-        fi <- scale[1] + fi * (scale[2] - scale[1])

-    }

-    if(wt_is_1) {

-        ## We need to shift to avoid ratios involving negative numbers and

-        ## we need to avoid having any fitness at 0, specially the wt.  If

-        ## any negative value, add the min, and shift by the magnitude of

-        ## the min to avoid any at 0.

+        if(wt_is_1) {

+            ## We need to shift to avoid ratios involving negative numbers and

+            ## we need to avoid having any fitness at 0, specially the wt.  If

+            ## any negative value, add the min, and shift by the magnitude of

+            ## the min to avoid any at 0.

-        ## If you use scale and wt_is_1, this will move the scale. It is

-        ## not possible to obtain a linear transformation that will keep

-        ## the min and max of the scale, and wt at 1.

-        min_fi <- min(fi)

-        if(min_fi < 0)

-            fi <- fi + 2 * abs(min(fi))

-        fi <- fi/fi[1]

-    }

-    

-    if(log) {

-        fi <- log(fi/fi[1]) + 1

+            ## If you use scale and wt_is_1, this will move the scale. It is

+            ## not possible to obtain a linear transformation that will keep

+            ## the min and max of the scale, and wt at 1.

+            min_fi <- min(fi)

+            if(min_fi < 0)

+                fi <- fi + 2 * abs(min(fi))

+            fi <- fi/fi[1]

+        }

+        

+        if(log) {

+            fi <- log(fi/fi[1]) + 1

+        }

+        m <- cbind(m, Fitness = fi)

+        if(min_accessible_genotypes > 0) {

+            if(count_accessible_g(m, accessible_th) >= min_accessible_genotypes) {

+                done <- TRUE

+            } else {

+                ## Cannot start again with a fitness column

+                m <- m[, -ncol(m), drop = FALSE]

+            }

+        } else {

+            done <- TRUE

+        }

     }

-    m <- cbind(m, Fitness = fi)

     class(m) <- c(class(m), "genotype_fitness_matrix")

     return(m)

 }

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

+Changes in version 2.3.9 (2017-07-0):

+	- Accessible genotypes in rfitness and plotFitnessLandscape

+

 Changes in version 2.3.8 (2017-07-08):

 	- PDBasline default is now 1.2.

@@ -16,28 +16,40 @@

 \usage{

 plotFitnessLandscape(x, show_labels = TRUE,

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

-                                   lty = c(1, 2, 3),

-                                   use_ggrepel = FALSE,

-                                   log = FALSE, max_num_genotypes = 2000, ...)

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

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

+                     use_ggrepel = FALSE,

+                     log = FALSE, max_num_genotypes = 2000,

+                     only_accessible = FALSE,

+                     accessible_th = 0,

+                     ...)

 \method{plot}{genotype_fitness_matrix}(x, show_labels = TRUE,

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

                                    lty = c(1, 2, 3),

                                    use_ggrepel = FALSE,

-                                   log = FALSE, max_num_genotypes = 2000, ...)

+                                   log = FALSE, max_num_genotypes = 2000,

+                                   only_accessible = FALSE,

+                                   accessible_th = 0,

+                                   ...)

 \method{plot}{evalAllGenotypes}(x, show_labels = TRUE,

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

                                    lty = c(1, 2, 3),

                                    use_ggrepel = FALSE,

-                                   log = FALSE, max_num_genotypes = 2000, ...)

+                                   log = FALSE, max_num_genotypes = 2000,

+                                   only_accessible = FALSE,

+                                   accessible_th = 0,

+                                   ...)

 \method{plot}{evalAllGenotypesMut}(x, show_labels = TRUE,

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

                                    lty = c(1, 2, 3),

                                    use_ggrepel = FALSE,

-                                   log = FALSE, max_num_genotypes = 2000, ...)

+                                   log = FALSE, max_num_genotypes = 2000,

+                                   only_accessible = FALSE,

+                                   accessible_th = 0,

+                                   ...)

 }

 \arguments{

@@ -89,6 +101,18 @@ plotFitnessLandscape(x, show_labels = TRUE,

   types of input, we make a call to \code{\link{evalAllGenotypes}}, and

   use this as the maximum.}

+\item{only_accessible}{If TRUE, show only accessible paths. A path is

+  considered accesible if, at each mutational step (i.e., with the

+  addition of each mutation) fitness increases by at least

+  \code{accessible_th}. If you set \code{only_accessible = TRUE}, the

+  number of genotypes displayed can be much smaller than the number of

+  existing genotypes if many of those genotypes are not accessible via

+  any path.}

+

+\item{accessible_th}{The threshold for the minimal change in fitness at

+  each mutation step (i.e., between successive genotypes) to be used if

+  \code{only_accessible = TRUE}.}

+

   \item{\dots}{

       Other arguments passed to \code{plot}. Not used for now.

   }

@@ -4,16 +4,15 @@

 \title{Generate random fitness.}

-\description{

-  Generate random fitness under a House of Cards, Rough Mount Fuji, or

-  additive model.

-}

+\description{ Generate random fitness landscapes under a House of Cards,

+  Rough Mount Fuji, or additive model.  }

 \usage{

 rfitness(g, c = 0.5, sd = 1, reference = "random", scale = NULL,

-         wt_is_1 = TRUE, log = FALSE)

+         wt_is_1 = TRUE, log = FALSE, min_accessible_genotypes = 0,

+         accessible_th = 0)

 }

@@ -48,7 +47,34 @@ rfitness(g, c = 0.5, sd = 1, reference = "random", scale = NULL,

 \item{log}{If TRUE, log-transform fitness.}

+

+\item{min_accessible_genotypes}{If larger than 0, the minimum number of

+  accessible genotypes in the fitness landscape. A genotype is

+  considered accessible if you can reach if from the wildtype by going

+  through at least one path where all changes in fitness are larger or

+  equal to \code{accessible_th}. The changes in fitness are considered

+  at each mutational step, i.e., at each addition of one mutation we

+  compute the difference between the genotype with \code{k + 1}

+  mutations minus the ancestor genotype with \code{k} mutations. Thus, a

+  genotype is considered accessible if there is at least one path where

+  fitness increases at each mutational step by at least

+  \code{accessible_th}.

+

+  If the condition is not satisfied, we continue generating random

+  fitness landscapes with the specified parameters until the condition

+  is satisfied.

 }

+

+\item{accessible_th}{The threshold for the minimal change in fitness at

+  each mutation step (i.e., between successive genotypes) that allows a

+  genotype to be regarded as accessible. This only applies if

+  \code{min_accessible_genotypes} is larger than 0.  So if you want to

+  allow small decreases in fitness in successive steps, use a small

+  negative value for \code{accessible_th}.  }

+

+} 

+

+

 \details{

   The model used here follows the Rough Mount Fuji model in Szendro et

new file mode 100644

@@ -0,0 +1,33 @@

+test_that("exercise accessible genotypes code", {

+    ## Because of compiler differences, this might only test the relevant

+    ## code in Linux. But conditions will be true in all systems.

+    

+    RNGkind("Mersenne-Twister")

+

+    set.seed(3)

+    r1 <- rfitness(4, reference = rep(0, 4), c = 2.1,

+                   min_accessible_genotypes = 1)

+    expect_true(max(r1[-1, "Fitness"]) >= 1.0)

+    

+    set.seed(3)

+    r1 <- rfitness(4, reference = rep(0, 4), c = 2.1,

+                   min_accessible_genotypes = 1,

+                   accessible_th = 0.06)

+    expect_true(max(r1[-1, "Fitness"]) >= 1.06)

+

+    set.seed(3)

+    r2 <- rfitness(4, reference = rep(0, 4), c = 2.1,

+                   min_accessible_genotypes = 2,

+                   accessible_th = 0.06)

+    expect_true(max(r2[-1, "Fitness"]) >= 1.06)

+

+    ## Exercising plotting code

+    plot(r2)

+    plot(r2, only_accessible = TRUE)

+    plot(r2, only_accessible = TRUE, accessible_th = 0.06)

+    ## this only show WT

+    plot(r2, only_accessible = TRUE, accessible_th = 0.2)

+   

+}) 

+

+set.seed(NULL)

@@ -4317,10 +4317,25 @@ samplePop(pancrPop)

 Now a simple multiple call to \texttt{oncoSimulIndiv} wrapped inside

 \Rfunction{mclapply}; this is basically the same we just did above. We set

-the class of the object to allow direct usage of \Rfunction{samplePop}.

+the class of the object to allow direct usage of

+\Rfunction{samplePop}. (Note: in Windows \texttt{mc.cores > 1} is not

+supported, so for the vignette to run in Windows, Linux, and Mac we

+explicitly set it here in the call to \Rfunction{mclapply}. For regular

+usage, you will not need to do this; just use whatever is appropriate for

+your operating system and number of cores. As well, we do not need any of

+this with \Rfunction{oncoSimulPop} because the code inside

+\Rfunction{oncoSimulPop} already takes care of setting \texttt{mc.cores}

+to 1 in Windows).

 <<>>=

 library(parallel)

+

+if(.Platform$OS.type == "windows") {

+    mc.cores <- 1

+} else {

+    mc.cores <- 2

+}

+

 p2 <- mclapply(1:4, function(x) oncoSimulIndiv(pancr,

                                                detectionSize = 1e7,

                                                keepEvery = 10),