Browse code

version 2.99.1: frequency-dependent fitness functionality

ramon diaz-uriarte (at Phelsuma) authored on 10/12/2020 11:41:53
Showing61 changed files

... ...
@@ -1,8 +1,8 @@
1 1
 Package: OncoSimulR
2 2
 Type: Package
3 3
 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 
4
-Version: 2.21.0
5
-Date: 2020-06-03
4
+Version: 2.99.1
5
+Date: 2020-12-10
6 6
 Authors@R: c(
7 7
 	      person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),	
8 8
  	   		     email = "rdiaz02@gmail.com"),
... ...
@@ -141,29 +141,32 @@ Author: Ramon Diaz-Uriarte [aut, cre],
141 141
 	Rosalia Palomino Cabrera [ctb],
142 142
 	Rafael Barrero Rodriguez [ctb],
143 143
 	Silvia Talavera Marcos [ctb],
144
-	Niklas Endres [ctb] 
144
+	Niklas Endres [ctb]  
145 145
 Maintainer: Ramon Diaz-Uriarte <rdiaz02@gmail.com>
146 146
 Description: Functions for forward population genetic simulation in
147 147
     asexual populations, with special focus on cancer progression.
148 148
     Fitness can be an arbitrary function of genetic interactions between
149 149
     multiple genes or modules of genes, including epistasis, order
150
-    restrictions in mutation accumulation, and order effects.  Mutation
151
-    rates can differ between genes, and we can include mutator/antimutator
152
-    genes (to model mutator phenotypes). Simulations use continuous-time
153
-    models and can include driver and passenger genes and modules.  Also
154
-    included are functions for: simulating random DAGs of the type found
155
-    in Oncogenetic Trees, Conjunctive Bayesian Networks, and other cancer
156
-    progression models; plotting and sampling from single or multiple
157
-    realizations of the simulations, including single-cell sampling;
158
-    plotting the parent-child relationships of the clones; generating
159
-    random fitness landscapes (from Rough Mount Fuji, House of Cards,
160
-    additive, NK, Ising, and Eggbox models) and plotting them.
150
+    restrictions in mutation accumulation, and order effects.  Fitness can
151
+    also be a function of the relative and absolute frequencies of other
152
+    genotypes (i.e., frequency-dependent fitness). Mutation rates can
153
+    differ between genes, and we can include mutator/antimutator genes (to
154
+    model mutator phenotypes). Simulating multi-species scenarios and
155
+    therapeutic interventions is also possible. Simulations use
156
+    continuous-time models and can include driver and passenger genes and
157
+    modules.  Also included are functions for: simulating random DAGs of
158
+    the type found in Oncogenetic Trees, Conjunctive Bayesian Networks,
159
+    and other cancer progression models; plotting and sampling from single
160
+    or multiple realizations of the simulations, including single-cell
161
+    sampling; plotting the parent-child relationships of the clones;
162
+    generating random fitness landscapes (Rough Mount Fuji, House of
163
+    Cards, additive, NK, Ising, and Eggbox models) and plotting them.
161 164
 biocViews: BiologicalQuestion, SomaticMutation
162 165
 License: GPL (>= 3)
163 166
 URL: https://github.com/rdiaz02/OncoSimul, https://popmodels.cancercontrol.cancer.gov/gsr/packages/oncosimulr/
164 167
 BugReports: https://github.com/rdiaz02/OncoSimul/issues
165
-Depends: R (>= 3.3.0)
166
-Imports: Rcpp (>= 0.12.4), parallel, data.table, graph, Rgraphviz, gtools, igraph, methods, RColorBrewer, grDevices, car, dplyr, smatr, ggplot2, ggrepel
168
+Depends: R (>= 3.5.0)
169
+Imports: Rcpp (>= 0.12.4), parallel, data.table, graph, Rgraphviz, gtools, igraph, methods, RColorBrewer, grDevices, car, dplyr, smatr, ggplot2, ggrepel, stringr
167 170
 Suggests: BiocStyle, knitr, Oncotree, testthat (>= 1.0.0), rmarkdown, bookdown, pander
168 171
 LinkingTo: Rcpp
169 172
 VignetteBuilder: knitr
... ...
@@ -12,8 +12,9 @@ export("oncoSimulPop", "oncoSimulIndiv", "samplePop",
12 12
        "to_Magellan",
13 13
        "Magellan_stats",
14 14
        "sampledGenotypes"
15
-       , "POM", "LOD"
16
-       , "diversityPOM", "diversityLOD"
15
+     , "POM", "LOD"
16
+     , "diversityPOM", "diversityLOD"
17
+##     , "meanCompositionPop", "compositionPop2"
17 18
        )
18 19
 
19 20
 S3method(plot, oncosimul)
... ...
@@ -56,8 +57,8 @@ importFrom("gtools", combinations, permutations, mixedorder)
56 57
 ## importFrom("compare", compare)
57 58
 importFrom("graphics", "axis", "box", "legend", "matplot", "par", "polygon")
58 59
 importFrom("methods", "as")
59
-importFrom("stats", "na.omit", "runif", "smooth.spline")
60
-importFrom("utils", "type.convert")
60
+importFrom("stats", "na.omit", "runif", "smooth.spline", "setNames")
61
+importFrom("utils", "type.convert", "combn")
61 62
 importFrom("RColorBrewer", "brewer.pal")
62 63
 importFrom("grDevices", "colorRampPalette", "hsv", "rainbow")
63 64
 importFrom("dplyr", "full_join", "left_join", "right_join", "%>%", "mutate",
... ...
@@ -65,8 +66,7 @@ importFrom("dplyr", "full_join", "left_join", "right_join", "%>%", "mutate",
65 66
 importFrom("smatr", "ma") ## for major axis regression in some tests
66 67
 importFrom("car", "linearHypothesis")
67 68
 ## importFrom("nem", "transitive.reduction") ## now in file nem_transitive_reduction.R
69
+importFrom("stringr", "regex", "str_extract_all", "str_replace_all")
68 70
 ## importFrom("slam", "simple_triplet_zero_matrix", ## "colapply_simple_triplet_matrix",
69 71
 ##            "col_sums")
70
-
71
-
72
-
72
+## importFrom("graphics", "segments", "stripchart")
... ...
@@ -1,4 +1,4 @@
1
-## Copyright 2016, 2017 Ramon Diaz-Uriarte
1
+## Copyright 2016-2021 Ramon Diaz-Uriarte
2 2
 
3 3
 ## This program is free software: you can redistribute it and/or modify
4 4
 ## it under the terms of the GNU General Public License as published by
... ...
@@ -1,4 +1,4 @@
1
-## Copyright 2013, 2014, 2015, 2016, 2017 Ramon Diaz-Uriarte
1
+## Copyright 2013-2021 Ramon Diaz-Uriarte
2 2
 
3 3
 ## This program is free software: you can redistribute it and/or modify
4 4
 ## it under the terms of the GNU General Public License as published by
... ...
@@ -46,7 +46,7 @@ oncoSimulSample <- function(Nindiv,
46 46
                             initSize = 500,
47 47
                             s = 0.1,
48 48
                             sh = -1,
49
-                            K = initSize/(exp(1) - 1),
49
+                            K = sum(initSize)/(exp(1) - 1),
50 50
                             minDetectDrvCloneSz = "auto",
51 51
                             extraTime = 0,
52 52
                             finalTime = 0.25 * 25 * 365,
... ...
@@ -364,7 +364,7 @@ oncoSimulPop <- function(Nindiv,
364 364
                          initSize = 500,
365 365
                          s = 0.1,
366 366
                          sh = -1,
367
-                         K = initSize/(exp(1) - 1),
367
+                         K = sum(initSize)/(exp(1) - 1),
368 368
                          keepEvery = sampleEvery, 
369 369
                          minDetectDrvCloneSz = "auto",
370 370
                          extraTime = 0,
... ...
@@ -454,7 +454,7 @@ oncoSimulIndiv <- function(fp,
454 454
                            initSize = 500,
455 455
                            s = 0.1,
456 456
                            sh = -1,
457
-                           K = initSize/(exp(1) - 1),
457
+                           K = sum(initSize)/(exp(1) - 1),
458 458
                            keepEvery = sampleEvery,
459 459
                            minDetectDrvCloneSz = "auto",
460 460
                            extraTime = 0,
... ...
@@ -506,12 +506,15 @@ oncoSimulIndiv <- function(fp,
506 506
                           "Exp" = "exp",
507 507
                           "McFarlandLog" = "mcfarlandlog",
508 508
                           "McFL" = "mcfarlandlog",
509
+                          "McFarlandLogD" = "mcfarlandlogd",
510
+                          "McFLD" = "mcfarlandlogd",
509 511
                           stop("No valid value for model")
510 512
                           )
511
-    if(initSize < 1)
512
-        stop("initSize < 1")
513
+    if(max(initSize) < 1)
514
+        stop("all initSize < 1")
513 515
     
514
-    if( (K < 1) && (model %in% c("McFL", "McFarlandLog") )) {
516
+    if( (K < 1) && (model %in% c("McFL", "McFarlandLog",
517
+                                 "McFLD", "McFarlandLogD") )) {
515 518
         stop("Using McFarland's model: K cannot be < 1")
516 519
     }       ##  if ( !(model %in% c("McFL", "McFarlandLog") )) {
517 520
             ## K <- 1 ## K is ONLY used for McFarland; set it to 1, to avoid
... ...
@@ -526,7 +529,7 @@ oncoSimulIndiv <- function(fp,
526 529
     }
527 530
     birth <- -99
528 531
 
529
-    if( (typeFitness == "mcfarlandlog") &&
532
+    if( (typeFitness %in% c("mcfarlandlog", "mcfarlandlogd")) &&
530 533
        (sampleEvery > 0.05)) {
531 534
         warning("With the McFarland model you often want smaller sampleEvery")
532 535
     }
... ...
@@ -534,8 +537,13 @@ oncoSimulIndiv <- function(fp,
534 537
     if(minDetectDrvCloneSz == "auto") {
535 538
         if(model %in% c("Bozic", "Exp") )
536 539
             minDetectDrvCloneSz <- 0
537
-        else if (model %in% c("McFL", "McFarlandLog"))
538
-            minDetectDrvCloneSz <- initSize
540
+        else if (model %in% c("McFL", "McFarlandLog",
541
+                              "McFLD", "McFarlandLogD")) {
542
+            ## if(length(initSize) > 1) {
543
+            ##     message("length(initSize) > 1; using the sum of values for minDetectDrvCloneSz")
544
+            ## }
545
+            minDetectDrvCloneSz <- sum(initSize)
546
+        }
539 547
         ## minDetectDrvCloneSz <- eFinalMf(initSize, s, detectionDrivers)
540 548
         else
541 549
             stop("Unknown model")
... ...
@@ -1,20 +1,20 @@
1 1
 # This file was generated by Rcpp::compileAttributes
2 2
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
3 3
 
4
-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_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list) {
5
-    .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_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list)
4
+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_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list) {
5
+    .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_, n2, p2, PDBaseline, cPDetect_i, checkSizePEvery, AND_DrvProbExit, fixation_list)
6 6
 }
7 7
 
8 8
 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) {
9 9
     .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)
10 10
 }
11 11
 
12
-evalRGenotype <- function(rG, rFE, verbose, prodNeg, calledBy_) {
13
-    .Call('OncoSimulR_evalRGenotype', PACKAGE = 'OncoSimulR', rG, rFE, verbose, prodNeg, calledBy_)
12
+evalRGenotype <- function(rG, rFE, spPop, verbose, prodNeg, calledBy_, currentTime) {
13
+    .Call('OncoSimulR_evalRGenotype', PACKAGE = 'OncoSimulR', rG, rFE, spPop, verbose, prodNeg, calledBy_, currentTime)
14 14
 }
15 15
 
16
-evalRGenotypeAndMut <- function(rG, rFE, muEF, full2mutator_, verbose, prodNeg) {
17
-    .Call('OncoSimulR_evalRGenotypeAndMut', PACKAGE = 'OncoSimulR', rG, rFE, muEF, full2mutator_, verbose, prodNeg)
16
+evalRGenotypeAndMut <- function(rG, rFE, muEF, spPop, full2mutator_, verbose, prodNeg, currentTime) {
17
+    .Call('OncoSimulR_evalRGenotypeAndMut', PACKAGE = 'OncoSimulR', rG, rFE, muEF, spPop, full2mutator_, verbose, prodNeg, currentTime)
18 18
 }
19 19
 
20 20
 accessibleGenotypes <- function(y, f, numMut, th) {
... ...
@@ -1,4 +1,4 @@
1
-## Copyright 2013, 2014, 2015 Ramon Diaz-Uriarte
1
+## Copyright 2013-2021 Ramon Diaz-Uriarte
2 2
 
3 3
 ## This program is free software: you can redistribute it and/or modify
4 4
 ## it under the terms of the GNU General Public License as published by
... ...
@@ -1,4 +1,4 @@
1
-## Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte
1
+## Copyright 2013-2021 Ramon Diaz-Uriarte
2 2
 
3 3
 ## This program is free software: you can redistribute it and/or modify
4 4
 ## it under the terms of the GNU General Public License as published by
... ...
@@ -125,146 +125,234 @@ to_Fitness_Matrix <- function(x, max_num_genotypes) {
125 125
 ## rfitness, do nothing. Well, it actually does something.
126 126
 
127 127
 
128
-to_genotFitness_std <- function(x, simplify = TRUE,
128
+##Modified
129
+to_genotFitness_std <- function(x,
130
+                                frequencyDependentFitness = FALSE,
131
+                                frequencyType = frequencyType,
132
+                                simplify = TRUE,
129 133
                                 min_filter_fitness = 1e-9,
130 134
                                 sort_gene_names = TRUE) {
131
-    ## Would break with output from allFitnessEffects and
132
-    ## output from allGenotypeAndMut
135
+  ## Would break with output from allFitnessEffects and
136
+  ## output from allGenotypeAndMut
133 137
 
134
-    ## For the very special and weird case of
135
-    ## a matrix but only a single gene so with a 0 and 1
136
-    ## No, this is a silly and meaningless case.
137
-    ## if( ( ncol(x) == 2 ) && (nrow(x) == 1) && (x[1, 1] == 1) ) {
138
+  ## For the very special and weird case of
139
+  ## a matrix but only a single gene so with a 0 and 1
140
+  ## No, this is a silly and meaningless case.
141
+  ## if( ( ncol(x) == 2 ) && (nrow(x) == 1) && (x[1, 1] == 1) ) {
138 142
 
139
-    ## } else  blabla:
143
+  ## } else  blabla:
140 144
 
141
-    if(! (inherits(x, "matrix") || inherits(x, "data.frame")) )
142
-        stop("Input must inherit from matrix or data.frame.")
143 145
 
144
-    ## if((ncol(x) > 2) && !(inherits(x, "matrix"))
145
-    ##     stop(paste0("Genotype fitness input either two-column data frame",
146
-    ##          " or a numeric matrix with > 2 columns."))
147
-    ## if( (ncol(x) > 2) && (nrow(x) == 1) )
148
-    ##     stop(paste0("It looks like you have a matrix for a single genotype",
149
-    ##                 " of a single gene. For this degenerate cases use",
150
-    ##                 " a data frame specification."))
146
+
147
+  ## if((ncol(x) > 2) && !(inherits(x, "matrix"))
148
+  ##     stop(paste0("Genotype fitness input either two-column data frame",
149
+  ##          " or a numeric matrix with > 2 columns."))
150
+  ## if( (ncol(x) > 2) && (nrow(x) == 1) )
151
+  ##     stop(paste0("It looks like you have a matrix for a single genotype",
152
+  ##                 " of a single gene. For this degenerate cases use",
153
+  ##                 " a data frame specification."))
154
+
155
+
156
+
157
+    if(! (inherits(x, "matrix") || inherits(x, "data.frame")) )
158
+      stop("Input must inherit from matrix or data.frame.")
151 159
 
152 160
     if(ncol(x) > 2) {
161
+      if (!frequencyDependentFitness){
153 162
         if(inherits(x, "matrix")) {
154
-            if(!is.numeric(x))
155
-                stop("A genotype fitness matrix/data.frame must be numeric.")
163
+          if(!is.numeric(x))
164
+            stop("A genotype fitness matrix/data.frame must be numeric.")
156 165
         } else if(inherits(x, "data.frame")) {
157
-            if(!all(unlist(lapply(x, is.numeric))))
158
-                stop("A genotype fitness matrix/data.frame must be numeric.")
166
+          if(!all(unlist(lapply(x, is.numeric))))
167
+            stop("A genotype fitness matrix/data.frame must be numeric.")
159 168
         }
160
-
161
-        ## We are expecting here a matrix of 0/1 where columns are genes
162
-        ## except for the last column, that is Fitness
163
-        ## Of course, can ONLY work with epistastis, NOT order
164
-        ## return(genot_fitness_to_epistasis(x))
165
-        if(any(duplicated(colnames(x))))
166
-            stop("duplicated column names")
167
-
168
-        cnfl <- which(colnames(x)[-ncol(x)] == "")
169
-        if(length(cnfl)) {
170
-            freeletter <- setdiff(LETTERS, colnames(x))[1]
171
-            if(length(freeletter) == 0) stop("Renaiming failed")
172
-            warning("One column named ''. Renaming to ", freeletter)
173
-            colnames(x)[cnfl] <- freeletter
169
+      }else{
170
+        if(!inherits(x, "data.frame"))
171
+           stop("Input must inherit from data.frame.")
172
+        if(ncol(x) == 0){
173
+          stop("You have an empty data.frame")
174 174
         }
175
-        if(!is.null(colnames(x)) && sort_gene_names) {
176
-            ncx <- ncol(x)
177
-            cnx <- colnames(x)[-ncx]
178
-            ocnx <- gtools::mixedorder(cnx)
179
-            if(!(identical(cnx[ocnx], cnx))) {
180
-                message("Sorting gene column names alphabetically")
181
-                x <- cbind(x[, ocnx, drop = FALSE], Fitness = x[, (ncx)])
182
-            }
175
+        if(!all(unlist(lapply(x[,-ncol(x)], is.numeric)))){
176
+          stop("All columns except last one must be numeric.")
183 177
         }
184
-
185
-        if(is.null(colnames(x))) {
186
-            ncx <- (ncol(x) - 1)
187
-            message("No column names: assigning gene names from LETTERS")
188
-            if(ncx > length(LETTERS))
189
-                stop("More genes than LETTERS; please give gene names",
190
-                     " as you see fit.")
191
-            colnames(x) <- c(LETTERS[1:ncx], "Fitness")
178
+        if(is.factor(x[, ncol(x)])) {
179
+          warning("Last column of genotype fitness is a factor. ",
180
+                  "Converting to character.")
181
+          x[, ncol(x)] <- as.character(x[, ncol(x)])
192 182
         }
183
+        if(!all(unlist(lapply(x[, ncol(x)], is.character)))){
184
+          stop("All elements in last column must be character.")
185
+        }
186
+      }
187
+
188
+      ## We are expecting here a matrix of 0/1 where columns are genes
189
+      ## except for the last column, that is Fitness
190
+      ## Of course, can ONLY work with epistastis, NOT order
191
+      ## return(genot_fitness_to_epistasis(x))
192
+      if(any(duplicated(colnames(x))))
193
+        stop("duplicated column names")
194
+
195
+      cnfl <- which(colnames(x)[-ncol(x)] == "")
196
+      if(length(cnfl)) {
197
+        freeletter <- setdiff(LETTERS, colnames(x))[1]
198
+        if(length(freeletter) == 0) stop("Renaiming failed")
199
+        warning("One column named ''. Renaming to ", freeletter)
200
+        colnames(x)[cnfl] <- freeletter
201
+      }
202
+      if(!is.null(colnames(x)) && sort_gene_names) {
203
+        ncx <- ncol(x)
204
+        cnx <- colnames(x)[-ncx]
205
+        ocnx <- gtools::mixedorder(cnx)
206
+        if(!(identical(cnx[ocnx], cnx))) {
207
+          message("Sorting gene column names alphabetically")
208
+          x <- cbind(x[, ocnx, drop = FALSE], Fitness = x[, (ncx)])
209
+        }
210
+      }
211
+
212
+      if(is.null(colnames(x))) {
213
+        ncx <- (ncol(x) - 1)
214
+        message("No column names: assigning gene names from LETTERS")
215
+        if(ncx > length(LETTERS))
216
+          stop("More genes than LETTERS; please give gene names",
217
+               " as you see fit.")
218
+        colnames(x) <- c(LETTERS[1:ncx], "Fitness")
219
+      }
220
+      if(!all(as.matrix(x[, -ncol(x)]) %in% c(0, 1) ))
221
+        stop("First ncol - 1 entries not in {0, 1}.")
193 222
 
194
-        if(!all(as.matrix(x[, -ncol(x)]) %in% c(0, 1) ))
195
-            stop("First ncol - 1 entries not in {0, 1}.")
196 223
     } else {
197
-        if(!inherits(x, "data.frame"))
198
-            stop("genotFitness: if two-column must be data frame")
199
-        ## Make sure no factors
200
-        if(is.factor(x[, 1])) {
201
-            warning("First column of genotype fitness is a factor. ",
202
-                    "Converting to character.")
203
-            x[, 1] <- as.character(x[, 1])
204
-            }
205
-        ## Make sure no numbers
206
-        if(any(is.numeric(x[, 1])))
207
-            stop(paste0("genotFitness: first column of data frame is numeric.",
208
-                        " Ambiguous and suggests possible error. If sure,",
209
-                        " enter that column as character"))
210
-
211
-        omarker <- any(grepl(">", x[, 1], fixed = TRUE))
212
-        emarker <- any(grepl(",", x[, 1], fixed = TRUE))
213
-        nogoodepi <- any(grepl(":", x[, 1], fixed = TRUE))
214
-        ## if(omarker && emarker) stop("Specify only epistasis or order, not both.")
215
-        if(nogoodepi && emarker) stop("Specify the genotypes separated by a ',', not ':'.")
216
-        if(nogoodepi && !emarker) stop("Specify the genotypes separated by a ',', not ':'.")
217
-        ## if(nogoodepi && omarker) stop("If you want order, use '>' and if epistasis ','.")
218
-        ## if(!omarker && !emarker) stop("You specified neither epistasis nor order")
219
-        if(omarker) {
220
-            ## do something. To be completed
221
-            stop("This code not yet ready")
222
-            ## You can pass to allFitnessEffects genotype -> fitness mappings that
223
-            ## involve epistasis and order. But they must have different
224
-            ## genes. Otherwise, it is not manageable.
224
+
225
+      if(!inherits(x, "data.frame"))
226
+        stop("genotFitness: if two-column must be data frame")
227
+      if(ncol(x) == 0){
228
+        stop("You have an empty data.frame")
229
+      }
230
+      ## Make sure no factors
231
+      if(is.factor(x[, 1])) {
232
+        warning("First column of genotype fitness is a factor. ",
233
+                "Converting to character.")
234
+        x[, 1] <- as.character(x[, 1])
235
+      }
236
+      ## Make sure no numbers
237
+      if(any(is.numeric(x[, 1])))
238
+        stop(paste0("genotFitness: first column of data frame is numeric.",
239
+                    " Ambiguous and suggests possible error. If sure,",
240
+                    " enter that column as character"))
241
+
242
+      omarker <- any(grepl(">", x[, 1], fixed = TRUE))
243
+      emarker <- any(grepl(",", x[, 1], fixed = TRUE))
244
+      nogoodepi <- any(grepl(":", x[, 1], fixed = TRUE))
245
+      ## if(omarker && emarker) stop("Specify only epistasis or order, not both.")
246
+      if(nogoodepi && emarker) stop("Specify the genotypes separated by a ',', not ':'.")
247
+      if(nogoodepi && !emarker) stop("Specify the genotypes separated by a ',', not ':'.")
248
+      ## if(nogoodepi && omarker) stop("If you want order, use '>' and if epistasis ','.")
249
+      ## if(!omarker && !emarker) stop("You specified neither epistasis nor order")
250
+      if(omarker) {
251
+        ## do something. To be completed
252
+        stop("This code not yet ready")
253
+        ## You can pass to allFitnessEffects genotype -> fitness mappings that
254
+        ## involve epistasis and order. But they must have different
255
+        ## genes. Otherwise, it is not manageable.
256
+      }
257
+      if( emarker || ( (!omarker) && (!emarker) && (!nogoodepi)) ) {
258
+        ## the second case above corresponds to passing just single letter genotypes
259
+        ## as there is not a single marker
260
+        x <- x[, c(1, 2), drop = FALSE]
261
+        if(!all(colnames(x) == c("Genotype", "Fitness"))) {
262
+          message("Column names of object not Genotype and Fitness.",
263
+                  " Renaming them assuming that is what you wanted")
264
+          colnames(x) <- c("Genotype", "Fitness")
265
+        }
266
+        if((!omarker) && (!emarker) && (!nogoodepi)) {
267
+          message("All single-gene genotypes as input to to_genotFitness_std")
225 268
         }
226
-        if( emarker || ( (!omarker) && (!emarker) && (!nogoodepi)) ) {
227
-            ## the second case above corresponds to passing just single letter genotypes
228
-            ## as there is not a single marker
229
-            x <- x[, c(1, 2), drop = FALSE]
230
-            if(!all(colnames(x) == c("Genotype", "Fitness"))) {
231
-                message("Column names of object not Genotype and Fitness.",
232
-                        " Renaming them assuming that is what you wanted")
233
-                colnames(x) <- c("Genotype", "Fitness")
234
-            }
235
-            if((!omarker) && (!emarker) && (!nogoodepi)) {
236
-                message("All single-gene genotypes as input to to_genotFitness_std")
237
-            }
238
-            ## Yes, we need to do this to  scale the fitness and put the "-"
239
-            x <- allGenotypes_to_matrix(x)
269
+        ## Yes, we need to do this to  scale the fitness and put the "-"
270
+        if(frequencyDependentFitness){
271
+          anywt <- which(x[, 1] == "WT")
272
+          if (length(anywt) > 1){
273
+              stop("WT should not appear more than once in fitness specification")
274
+              ## stop("WT must appear once.")
275
+          }
276
+          ## if(length(anywt) == 0) {
277
+          ##     x <- rbind(data.frame(Genotype = "WT", Fitness = "0"),
278
+          ##                x)
279
+          ## }
280
+          if(is.factor(x[, ncol(x)])) {
281
+            warning("Second column of genotype fitness is a factor. ",
282
+                    "Converting to character.")
283
+            x[, ncol(x)] <- as.character(x[, ncol(x)])
284
+          }
240 285
         }
286
+
287
+        x <- allGenotypes_to_matrix(x, frequencyDependentFitness)
288
+      }
241 289
     }
242 290
     ## And, yes, scale all fitnesses by that of the WT
243
-    whichroot <- which(rowSums(x[, -ncol(x), drop = FALSE]) == 0)
244
-    if(length(whichroot) == 0) {
291
+
292
+    if (!frequencyDependentFitness){
293
+      whichroot <- which(rowSums(x[, -ncol(x), drop = FALSE]) == 0)
294
+      if(length(whichroot) == 0) {
245 295
         warning("No wildtype in the fitness landscape!!! Adding it with fitness 1.")
246 296
         x <- rbind(c(rep(0, ncol(x) - 1), 1), x)
247
-    } else if(x[whichroot, ncol(x)] != 1) {
297
+      } else if(x[whichroot, ncol(x)] != 1) {
248 298
         warning("Fitness of wildtype != 1.",
249 299
                 " Dividing all fitnesses by fitness of wildtype.")
250 300
         vwt <- x[whichroot, ncol(x)]
251 301
         x[, ncol(x)] <- x[, ncol(x)]/vwt
302
+      }
252 303
     }
304
+
253 305
     if(any(is.na(x)))
254 306
         stop("NAs in fitness matrix")
255
-    ## Make sure correct class
256
-    if(is.data.frame(x)) x <- as.matrix(x)
257
-    stopifnot(inherits(x, "matrix"))
258
-    ## if(simplify) {
259
-    ##     return(x[x[, ncol(x)] > min_filter_fitness, , drop = FALSE])
260
-    ## } else {
261
-    ##     return(x)
262
-    ## }
263
-    if(simplify) {
264
-        x <- x[x[, ncol(x)] > min_filter_fitness, , drop = FALSE]
307
+    if(!frequencyDependentFitness) {
308
+        if(is.data.frame(x)) 
309
+        x <- as.matrix(x)
310
+        stopifnot(inherits(x, "matrix"))
311
+
312
+        if(simplify) {
313
+           x <- x[x[, ncol(x)] > min_filter_fitness, , drop = FALSE]  
314
+        }
315
+        class(x) <- c("matrix", "genotype_fitness_matrix")
316
+    } else { ## frequency-dependent fitness
317
+        ## RDU: what is this for? It converts to numbers.
318
+        ## But it does not work reliably. It does not work when given as "n_"
319
+        ## We will do this at end. in full_FDF_spec
320
+        ## conversionTable_o <- conversionTable(ncol(x) - 1, frequencyType)
321
+        
322
+        ## x[, ncol(x)] <- sapply(x[, ncol(x)],
323
+        ##                        function(x){ findAndReplace(x, conversionTable_o)})
324
+        
325
+    if(frequencyType == "auto"){
326
+      ch <- paste(as.character(x[, ncol(x)]), collapse = "")
327
+      if( grepl("f_", ch, fixed = TRUE) ){
328
+        frequencyType = "rel"
329
+        pattern <- stringr::regex("f_(\\d*_*)*")
330
+        
331
+      } else if ( grepl("n_", ch, fixed = TRUE) ){
332
+        frequencyType = "abs"
333
+        pattern <- stringr::regex("n_(\\d*_*)*")
334
+        
335
+      } else { stop("No pattern found when frequencyType set to 'auto'") }
336
+        
337
+    } else if(frequencyType == "abs"){
338
+            pattern <- stringr::regex("n_(\\d*_*)*")
339
+        } else {
340
+            pattern <- stringr::regex("f_(\\d*_*)*")
341
+        }
342
+
343
+        regularExpressionVector <-
344
+            unique(unlist(lapply(x[, ncol(x)],
345
+                                 function(z) {stringr::str_extract_all(string = z,
346
+                                                                       pattern = pattern)})))
347
+        
348
+        if((!all(regularExpressionVector %in% fVariablesN(ncol(x) - 1, frequencyType))) |
349
+           !(length(intersect(regularExpressionVector,
350
+                              fVariablesN(ncol(x) - 1, frequencyType)) >= 1) )){
351
+            stop("There are some errors in fitness column")
352
+        }
353
+
265 354
     }
266
-    class(x) <- c("matrix", "genotype_fitness_matrix")
267
-    return(x)
355
+  return(x)
268 356
 }
269 357
 
270 358
 ## Deprecated after flfast
... ...
@@ -416,7 +504,8 @@ genot_fitness_to_epistasis <- function(x) {
416 504
 
417 505
 
418 506
 
419
-allGenotypes_to_matrix <- function(x) {
507
+allGenotypes_to_matrix <- function(x,
508
+                                   frequencyDependentFitness = FALSE) {
420 509
     ## Makes no sense to allow passing order: the matrix would have
421 510
     ## repeated rows. A > B and B > A both have exactly A and B
422 511
 
... ...
@@ -424,38 +513,69 @@ allGenotypes_to_matrix <- function(x) {
424 513
     ## a matrix with 0/1 in a column for each gene and a final column of
425 514
     ## Fitness
426 515
 
427
-    if(is.factor(x[, 1])) {
428
-        warning("First column of genotype fitness is a factor. ",
429
-                "Converting to character.")
516
+    if (is.factor(x[, 1])) {
517
+        warning(
518
+            "First column of genotype fitness is a factor. ",
519
+            "Converting to character."
520
+        )
430 521
         x[, 1] <- as.character(x[, 1])
431 522
     }
523
+    if (frequencyDependentFitness) {
524
+        if (is.factor(x[, ncol(x)])) {
525
+            warning(
526
+                "Second column of genotype fitness is a factor. ",
527
+                "Converting to character."
528
+            )
529
+            x[, ncol(x)] <- as.character(x[, ncol(x)])
530
+        }
531
+    }
532
+
432 533
     ## A WT can be specified with string "WT"
433 534
     anywt <- which(x[, 1] == "WT")
434
-    if(length(anywt) > 1) stop("More than 1 WT")
435
-    if(length(anywt) == 1) {
535
+    if (length(anywt) > 1) stop("More than 1 WT")
536
+    if (length(anywt) == 1) {
436 537
         fwt <- x[anywt, 2]
437 538
         x <- x[-anywt, ]
438 539
         ## Trivial case of passing just a WT?
439 540
     } else {
440
-        warning("No WT genotype. Setting its fitness to 1.")
441
-        fwt <- 1
541
+        if (!frequencyDependentFitness) {
542
+            warning("No WT genotype. Setting its fitness to 1.")
543
+            fwt <- 1
544
+        } else {
545
+            fwt <- NA
546
+            ##   message("No WT genotype in FDF: setting it to 0.")
547
+        }
442 548
     }
443
-    splitted_genots <- lapply(x$Genotype,
444
-                             function(z) nice.vector.eo(z, ","))
549
+    splitted_genots <- lapply(
550
+        x$Genotype,
551
+        function(z) nice.vector.eo(z, ",")
552
+    )
445 553
 
446 554
     all_genes <- sort(unique(unlist(splitted_genots)))
447 555
 
448 556
     m <- matrix(0, nrow = length(splitted_genots), ncol = length(all_genes))
449
-    the_match <- lapply(splitted_genots,
450
-                        function(z) match(z, all_genes))
557
+    the_match <- lapply(
558
+        splitted_genots,
559
+        function(z) match(z, all_genes)
560
+    )
451 561
     ## A lot simpler with a loop
452
-    for(i in 1:nrow(m)) {
562
+    for (i in 1:nrow(m)) {
453 563
         m[i, the_match[[i]]] <- 1
454 564
     }
455 565
     m <- cbind(m, x[, 2])
456 566
     colnames(m) <- c(all_genes, "Fitness")
457
-    m <- rbind(c(rep(0, length(all_genes)), fwt),
458
-               m)
567
+    if(!is.na(fwt))
568
+        m <- rbind(c(rep(0, length(all_genes)), fwt), m)
569
+
570
+    if (frequencyDependentFitness) {
571
+        m <- as.data.frame(m)
572
+        m[, 1:length(all_genes)] <- apply(
573
+            m[, 1:length(all_genes)],
574
+            2,
575
+            as.numeric
576
+        )
577
+        m[, ncol(m)] <- as.character(m[, ncol(m)])
578
+    }
459 579
     ## Ensure sorted
460 580
     ## m <- data.frame(m)
461 581
     rs <- rowSums(m[, -ncol(m), drop = FALSE])
... ...
@@ -466,6 +586,7 @@ allGenotypes_to_matrix <- function(x) {
466 586
 
467 587
 
468 588
 
589
+
469 590
 Magellan_stats <- function(x, max_num_genotypes = 2000,
470 591
                            verbose = FALSE,
471 592
                            use_log = FALSE,
... ...
@@ -1,4 +1,4 @@
1
-## Copyright 2013, 2014, 2015 Ramon Diaz-Uriarte
1
+## Copyright 2013-2021 Ramon Diaz-Uriarte
2 2
 
3 3
 ## This program is free software: you can redistribute it and/or modify
4 4
 ## it under the terms of the GNU General Public License as published by
... ...
@@ -1,4 +1,4 @@
1
-## Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte
1
+## Copyright 2013-2021 Ramon Diaz-Uriarte
2 2
 
3 3
 ## This program is free software: you can redistribute it and/or modify
4 4
 ## it under the terms of the GNU General Public License as published by
... ...
@@ -41,29 +41,41 @@ allMutatorEffects <- function(epistasis = NULL,
41 41
 
42 42
 evalAllGenotypesMut <- function(mutatorEffects,
43 43
                                 max = 256,
44
-                                addwt = FALSE) {
44
+                                addwt = FALSE,
45
+                                spPopSizes = NULL,
46
+                                currentTime = 0
47
+                                ) {
45 48
     evalAllGenotypesORMut(
46 49
         fmEffects = mutatorEffects,
47 50
         order = FALSE,
48 51
         max = max,
49 52
         model = "",
50
-        calledBy_= "evalGenotypeMut"
53
+        spPopSizes = spPopSizes,
54
+        calledBy_= "evalGenotypeMut",
55
+        currentTime = currentTime
51 56
     )
52 57
 }
53 58
 
54
-evalGenotypeMut <- function(genotype, mutatorEffects,
55
-                         verbose = FALSE,
56
-                         echo = FALSE) {
59
+evalGenotypeMut <- function(genotype, 
60
+                            mutatorEffects,
61
+                            spPopSizes = NULL,
62
+                            verbose = FALSE,
63
+                            echo = FALSE,
64
+                            currentTime = 0
65
+                            ) {
66
+    
57 67
     if(inherits(mutatorEffects, "fitnessEffects"))
58 68
         stop("You are trying to get the mutator effects of a fitness specification. ",
59 69
              "You did not pass an object of class mutatorEffects.")
60 70
     evalGenotypeORMut(genotype = genotype,
61
-                       fmEffects = mutatorEffects,
62
-                       verbose = verbose,
63
-                       echo = echo,
64
-                       model  = "" ,
65
-                       calledBy_= "evalGenotypeMut"
66
-                       )
71
+                      fmEffects = mutatorEffects,
72
+                      spPopSizes = spPopSizes,
73
+                      verbose = verbose,
74
+                      echo = echo,
75
+                      model  = "" ,
76
+                      calledBy_= "evalGenotypeMut",
77
+                      currentTime = currentTime
78
+                      )
67 79
 
68 80
 }
69 81
 
... ...
@@ -1,5 +1,5 @@
1 1
 
2
-## Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte
2
+## Copyright 2013-2021 Ramon Diaz-Uriarte
3 3
 
4 4
 ## This program is free software: you can redistribute it and/or modify
5 5
 ## it under the terms of the GNU General Public License as published by
... ...
@@ -414,6 +414,185 @@ getGeneIDNum <- function(geneModule, geneNoInt, fitnessLandscape_gene_id,
414 414
     )
415 415
 }
416 416
 
417
+## Next two in crate_flvars_fitvars.
418
+
419
+## ## generate fitnesLanscapeVariables for FDF
420
+## ## fitness landscape as data frame with cols as genes,
421
+## ## and last column is genotype
422
+## create_flvars <- function(x, frequencyType) {
423
+##     x <- x[, -ncol(x), drop = FALSE]
424
+##     pasted <- apply(x, 1, function(z) paste(which(z == 1), collapse = "_"))
425
+##     npasted <- apply(x, 1, function(z) paste(colnames(x)[which(z == 1)], collapse = "_"))
426
+##     if(frequencyType == "abs") {
427
+##         tmp <- paste0("n_", pasted)
428
+##     } else {
429
+##         tmp <- paste0("f_", pasted)
430
+##     }
431
+##     names(tmp) <- npasted
432
+##     return(tmp)
433
+## }
434
+
435
+## ## fitness specification with letters -> fitness specification with numbers
436
+## names_to_nums_fitness <- function(fitness, genotFitness) {
437
+##     from_subst_pattern <- colnames(genotFitness[, -ncol(genotFitness)])
438
+##     to_subst_pattern <- as.character(1:(ncol(genotFitness) - 1))
439
+##     names(to_subst_pattern) <- from_subst_pattern
440
+##     return(stringr::str_replace_all(fitness, to_subst_pattern))
441
+## }
442
+
443
+## genotFitnes and frequency type -> fitnesLanscapeVariables for FDF and
444
+##          fitness with numbers, not names
445
+##   Done in a single function since both operations make
446
+##   the same assumptions
447
+create_flvars_fitvars <- function(genotFitness, frequencyType) {
448
+    x <- genotFitness[, -ncol(genotFitness), drop = FALSE]
449
+    pasted <- apply(x, 1, function(z) paste(sort(which(z == 1)), collapse = "_"))
450
+    npasted <- apply(x, 1, function(z) paste(sort(colnames(x)[which(z == 1)]), collapse = "_"))
451
+    if(frequencyType == "abs") {
452
+        prefix <- "n_"
453
+        prefixre <- "^n_"
454
+    } else {
455
+        prefix <- "f_"
456
+        prefixre <- "^f_"
457
+    }
458
+    flvars <- paste0(prefix, pasted)
459
+    names(flvars) <- npasted
460
+
461
+    ## make sure we get f_1_2 and not f_2_1, etc
462
+    flvars2 <- flvars
463
+    names(flvars2) <- paste0(prefix, names(flvars))
464
+
465
+    rmwt <- which(flvars2 == prefix)
466
+    if(length(rmwt)) flvars2 <- flvars2[-rmwt] ## rm this.
467
+
468
+    ## Need to rev the vector, to ensure larger patterns come first
469
+    ## and to place "f_" as last.
470
+    rflvars2 <- rev(flvars2)
471
+    count_seps <- stringr::str_count(rflvars2, stringr::fixed("_"))
472
+    
473
+    if(any(diff(count_seps) > 0)) {
474
+        warning("flvars not ordered?",
475
+                "Check the conversion of gene names to numbers in fitness spec")
476
+        rflvars2 <- rflvars2[order(count_seps, decreasing = TRUE)]
477
+    }
478
+
479
+    ## Users can pass in many possible orderings. Get all.
480
+    full_rflvars <- all_orders_fv(rflvars2, prefix, prefixre)
481
+    Fitness_as_fvars <- stringr::str_replace_all(genotFitness$Fitness,
482
+                                                 stringr::fixed(full_rflvars))
483
+    return(list(flvars = flvars,
484
+                Fitness_as_fvars = Fitness_as_fvars))
485
+}
486
+
487
+
488
+## named vector, prefixes -> all possible combinations of the name
489
+all_orders_fv <- function(x, prefix, prefixre) {
490
+    f1 <- function(zz, prefix, prefixre) {
491
+        z <- gsub(prefixre, "", names(zz))
492
+        spl <- strsplit(z, "_")[[1]]
493
+        if(length(spl) <= 1) {
494
+            return(zz)
495
+        } else {
496
+            pp <- gtools::permutations(length(spl), length(spl), spl)
497
+            ppo <- apply(pp, 1, function(u) paste(u, collapse = "_"))
498
+            ppoo <- rep(zz, length(ppo))
499
+            names(ppoo) <- paste0(prefix, ppo)
500
+            return(ppoo)
501
+        }
502
+    }
503
+
504
+    ## I cannot use lapply, as it strips the names
505
+    out <- vector(mode = "character", length = 0)
506
+    for(i in 1:length(x)) {
507
+        out <- c(out, f1(x[i], prefix, prefixre))
508
+    }
509
+    return(out)
510
+}
511
+
512
+## ## character vector, named replacement -> replaced vector
513
+## ## named_replace: names are the (fixed) pattern, value the replacement
514
+## ## stringr::str_replace_all seems too smart and does not respect order
515
+## my_mgsub <- function(x, named_replace) {
516
+##     nn <- names(named_replace)
517
+##     xr <- x
518
+##     for(i in 1:length(named_replace)) {
519
+##         xr <- gsub(nn[i], named_replace[i], xr, fixed = TRUE)
520
+##     }
521
+##     return(xr)
522
+## }
523
+
524
+
525
+##New function
526
+fVariablesN <- function (g, frequencyType) {
527
+
528
+  if (is.null(g) | g == 0)
529
+    stop("Number of genes must be integer > 0")
530
+
531
+  combinationsList <- list()
532
+  for (i in 0:g) {
533
+    combinationsList <- append(combinationsList,
534
+                                 combn(1:g, i, list, simplify = TRUE))
535
+  }
536
+
537
+  if (frequencyType == "abs"){
538
+    fsVector <-sapply(sapply(combinationsList,
539
+                             function(x) paste0(x, collapse = "_")),
540
+                      function(x) paste0("n_", x))
541
+  }else{
542
+    fsVector <-sapply(sapply(combinationsList,
543
+                             function(x) paste0(x, collapse = "_")),
544
+                      function(x) paste0("f_", x))
545
+  }
546
+
547
+  return (fsVector)
548
+}
549
+
550
+fVariablesL <- function (g, frequencyType) {
551
+
552
+  if (is.null(g) | g == 0)
553
+    stop("Number of genes must be integer > 0")
554
+
555
+  if(g > length(LETTERS))
556
+    stop(paste0("Number of genes must be < length(LETTERS).",
557
+                " Please specify variables with numbers"))
558
+
559
+  combinationsList <- list()
560
+  for (i in 0:g) {
561
+    combinationsList <- append(combinationsList,
562
+                               combn(LETTERS[1:g], i, list, simplify = TRUE))
563
+  }
564
+
565
+  if (frequencyType == "abs"){
566
+    fsVector <-sapply(sapply(combinationsList,
567
+                             function(x) paste0(x, collapse = "_")),
568
+                      function(x) paste0("n_", x))
569
+  }else{
570
+    fsVector <-sapply(sapply(combinationsList,
571
+                             function(x) paste0(x, collapse = "_")),
572
+                      function(x) paste0("f_", x))
573
+  }
574
+
575
+  return (fsVector)
576
+}
577
+
578
+## ## Assuming we are using the full fitness landscapes (i.e., none of
579
+## ## setting genotypes with 0 fitness are absent from the table)
580
+## conversionTable <- function(g, frequencyType){
581
+##   df <- data.frame(let = fVariablesL(g, frequencyType)[-1],
582
+##                    num = fVariablesN(g, frequencyType)[-1],
583
+##                    stringsAsFactors = FALSE)
584
+##   return (df)
585
+## }
586
+
587
+## findAndReplace <- function(str, conversionTable_input){
588
+
589
+##   pattern <- rev(setNames(as.character(conversionTable_input$num),
590
+##                       conversionTable_input$let))
591
+
592
+##   str <- stringr::str_replace_all(string = str,
593
+##                                   pattern = pattern)
594
+##   return(str)
595
+## }
417 596
 
418 597
 allFitnessORMutatorEffects <- function(rT = NULL,
419 598
                                        epistasis = NULL,
... ...
@@ -424,109 +603,117 @@ allFitnessORMutatorEffects <- function(rT = NULL,
424 603
                                        keepInput = TRUE,
425 604
                                        genotFitness = NULL,
426 605
                                        ## refFE = NULL,
427
-                                       calledBy = NULL) {
428
-    ## From allFitnessEffects. Generalized so we deal with Fitness
429
-    ## and mutator.
430
-    
431
-    ## restrictions: the usual rt
606
+                                       calledBy = NULL,
607
+                                       frequencyDependentFitness = FALSE,
608
+                                       frequencyType = "freq_dep_not_used"){
609
+                                       #spPopSizes = NULL) {
610
+  ## From allFitnessEffects. Generalized so we deal with Fitness
611
+  ## and mutator.
432 612
 
433
-    ## epistasis: as it says, with the ":"
613
+  ## restrictions: the usual rt
434 614
 
435
-    ## orderEffects: the ">"
436
-    
437
-    ## All of the above can be genes or can be modules (if you pass a
438
-    ## geneToModule)
615
+  ## epistasis: as it says, with the ":"
439 616
 
440
-    ## rest: rest of genes, with fitness
617
+  ## orderEffects: the ">"
441 618
 
619
+  ## All of the above can be genes or can be modules (if you pass a
620
+  ## geneToModule)
442 621
 
443
-    ## For epistasis and order effects we create the output object but
444
-    ## missing the numeric ids of genes. With rT we do it in one go, as we
445
-    ## already know the mapping of genes to numeric ids. We could do the
446
-    ## same in epistasis and order, but we would be splitting twice
447
-    ## (whereas for rT extracting the names is very simple).
622
+  ## rest: rest of genes, with fitness
448 623
 
449
-    ## called appropriately?
450
-    if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
451
-        stop("How did you call this function?. Bug.")
452
-    
453
-    if(calledBy == "allMutatorEffects") {
454
-        ## very paranoid check
455
-        if( !is.null(rT) || !is.null(orderEffects) ||
456
-            !is.null(drvNames) || !is.null(genotFitness))
457
-            stop("allMutatorEffects called with forbidden arguments.",
458
-                 "Is this an attempt to subvert the function?")
459
-    }
460
-    
624
+
625
+  ## For epistasis and order effects we create the output object but
626
+  ## missing the numeric ids of genes. With rT we do it in one go, as we
627
+  ## already know the mapping of genes to numeric ids. We could do the
628
+  ## same in epistasis and order, but we would be splitting twice
629
+  ## (whereas for rT extracting the names is very simple).
630
+
631
+  ## called appropriately?
632
+  if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
633
+    stop("How did you call this function?. Bug.")
634
+
635
+  if(calledBy == "allMutatorEffects") {
636
+    ## very paranoid check
637
+    if( !is.null(rT) || !is.null(orderEffects) ||
638
+        !is.null(drvNames) || !is.null(genotFitness))
639
+      stop("allMutatorEffects called with forbidden arguments.",
640
+           "Is this an attempt to subvert the function?")
641
+  }
642
+
643
+  if(!frequencyDependentFitness) {
461 644
     rtNames <- NULL
462 645
     epiNames <- NULL
463 646
     orNames <- NULL
464 647
     if(!is.null(rT)) {
465
-        ## This is really ugly, but to prevent the stringsAsFactors I need it here:
466
-        rT$parent <- as.character(rT$parent)
467
-        rT$child <- as.character(rT$child)
468
-        rT$typeDep <- as.character(rT$typeDep)
469
-        rtNames <- unique(c(rT$parent, rT$child))
648
+      ## This is really ugly, but to prevent the stringsAsFactors I need it here:
649
+      rT$parent <- as.character(rT$parent)
650
+      rT$child <- as.character(rT$child)
651
+      rT$typeDep <- as.character(rT$typeDep)
652
+      rtNames <- unique(c(rT$parent, rT$child))
470 653
     }
654
+    
655
+    #if(!is.null(spPopSizes))
656
+      #warning("spPopSizes will be considered NULL if frequencyDependentFitness = FALSE")
657
+    
471 658
     if(!is.null(epistasis)) {
472
-        long.epistasis <- to.long.epist.order(epistasis, ":")
473
-        ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
474
-        ## deal with the possible negative signs
475
-        epiNames <- setdiff(unique(
476
-            unlist(lapply(long.epistasis,
477
-                          function(x) lapply(x$ids,
478
-                                             function(z) strsplit(z, "^-"))))),
479
-                            "")
659
+      long.epistasis <- to.long.epist.order(epistasis, ":")
660
+      ## epiNames <- unique(unlist(lapply(long.epistasis, function(x) x$ids)))
661
+      ## deal with the possible negative signs
662
+      epiNames <- setdiff(unique(
663
+        unlist(lapply(long.epistasis,
664
+                      function(x) lapply(x$ids,
665
+                                         function(z) strsplit(z, "^-"))))),
666
+        "")
480 667
     } else {
481
-        long.epistasis <- list()
668
+      long.epistasis <- list()
482 669
     }
483 670
     if(!is.null(orderEffects)) {
484
-        long.orderEffects <- to.long.epist.order(orderEffects, ">")
485
-        orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
671
+      long.orderEffects <- to.long.epist.order(orderEffects, ">")
672
+      orNames <- unique(unlist(lapply(long.orderEffects, function(x) x$ids)))
486 673
     } else {
487
-        long.orderEffects <- list()
674
+      long.orderEffects <- list()
488 675
     }
489 676
     allModuleNames <- unique(c(rtNames, epiNames, orNames))
490 677
     if(is.null(geneToModule)) {
491
-        gMOneToOne <- TRUE
492
-        geneToModule <- geneModuleNull(allModuleNames)
678
+      gMOneToOne <- TRUE
679
+      geneToModule <- geneModuleNull(allModuleNames)
493 680
     } else {
494
-        gMOneToOne <- FALSE
495
-        if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
496
-            stop(paste("Some values in geneToModule not present in any of",
497
-                       " rT, epistasis, or order effects"))
498
-        if(any(is.na(match(allModuleNames, names(geneToModule)))))
499
-            stop(paste("Some values in rT, epistasis, ",
500
-                       "or order effects not in geneToModule"))
681
+      gMOneToOne <- FALSE
682
+      if(any(is.na(match(setdiff(names(geneToModule), "Root"), allModuleNames))))
683
+        stop(paste("Some values in geneToModule not present in any of",
684
+                   " rT, epistasis, or order effects"))
685
+      if(any(is.na(match(allModuleNames, names(geneToModule)))))
686
+        stop(paste("Some values in rT, epistasis, ",
687
+                   "or order effects not in geneToModule"))
501 688
     }
502 689
     geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
503
-    
690
+
504 691
     idm <- unique(geneModule$ModuleNumID)
505 692
     names(idm) <- unique(geneModule$Module)
506 693
 
507 694
     if(!is.null(rT)) {
508
-        checkRT(rT)
509
-        long.rt <- to.long.rt(rT, idm)
695
+      checkRT(rT)
696
+      long.rt <- to.long.rt(rT, idm)
510 697
     } else {
511
-        long.rt <- list() ## yes, we want an object of length 0
698
+      long.rt <- list() ## yes, we want an object of length 0
512 699
     }
513 700
 
514 701
     ## Append the numeric ids to epistasis and order
515 702
     if(!is.null(epistasis)) {
516
-        long.epistasis <- lapply(long.epistasis,
517
-                                 function(x)
518
-                                     addIntID.epist.order(x, idm,
519
-                                                          sort = TRUE,
520
-                                                          sign = TRUE))
703
+      long.epistasis <- lapply(long.epistasis,
704
+                               function(x)
705
+                                 addIntID.epist.order(x, idm,
706
+                                                      sort = TRUE,
707
+                                                      sign = TRUE))
521 708
     }
522 709
     if(!is.null(orderEffects)) {
523
-        long.orderEffects <- lapply(long.orderEffects,
524
-                                    function(x)
525
-                                        addIntID.epist.order(x, idm,
526
-                                                             sort = FALSE,
527
-                                                             sign = FALSE))
710
+      long.orderEffects <- lapply(long.orderEffects,
711
+                                  function(x)
712
+                                    addIntID.epist.order(x, idm,
713
+                                                         sort = FALSE,
714
+                                                         sign = FALSE))
528 715
     }
529
-    
716
+
530 717
     if(!is.null(noIntGenes)) {
531 718
         if(inherits(noIntGenes, "character")) {
532 719
             wm <- paste("noIntGenes is a character vector.",
... ...
@@ -561,69 +748,69 @@ allFitnessORMutatorEffects <- function(rT = NULL,
561 748
                                 s = noIntGenes,
562 749
                                 stringsAsFactors = FALSE)
563 750
     } else {
564
-        geneNoInt <- data.frame()
751
+      geneNoInt <- data.frame()
565 752
     }
566
-    
753
+
567 754
     if(is.null(genotFitness)) {
568
-        genotFitness <- matrix(NA, nrow = 0, ncol = 1)
569
-        fitnessLandscape_df <- data.frame()
570
-        fitnessLandscape_gene_id <- data.frame()
755
+      genotFitness <- matrix(NA, nrow = 0, ncol = 1)
756
+      fitnessLandscape_df <- data.frame()
757
+      fitnessLandscape_gene_id <- data.frame()
571 758
     } else {
572
-        ## Yes, I am duplicating stuff for now.
573
-        ## This makes life simpler in C++:
574
-        ## In the map, the key is the genotype name, as
575
-        ## cnn <- colnames(genotFitness)[-ncol(genotFitness)]
576
-        cnn <- 1:(ncol(genotFitness) - 1)
577
-        gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
578
-                     function(x) paste(cnn[as.logical(x)],
579
-                                       collapse = ", "))
580
-        ## rownames(genotFitness) <- gfn
581
-        fitnessLandscape_df <-
582
-            data.frame(Genotype = gfn,
583
-                       Fitness = genotFitness[, ncol(genotFitness)],
584
-                       stringsAsFactors = FALSE)
585
-        fitnessLandscape_gene_id <- data.frame(
586
-            Gene = colnames(genotFitness)[-ncol(genotFitness)],
587
-            GeneNumID = cnn,
588
-            stringsAsFactors = FALSE)
589
-        
759
+      ## Yes, I am duplicating stuff for now.
760
+      ## This makes life simpler in C++:
761
+      ## In the map, the key is the genotype name, as
762
+      ## cnn <- colnames(genotFitness)[-ncol(genotFitness)]
763
+      cnn <- 1:(ncol(genotFitness) - 1)
764
+      gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
765
+                   function(x) paste(cnn[as.logical(x)],
766
+                                     collapse = ", "))
767
+      ## rownames(genotFitness) <- gfn
768
+      fitnessLandscape_df <-
769
+        data.frame(Genotype = gfn,
770
+                   Fitness = genotFitness[, ncol(genotFitness)],
771
+                   stringsAsFactors = FALSE)
772
+      fitnessLandscape_gene_id <- data.frame(
773
+        Gene = colnames(genotFitness)[-ncol(genotFitness)],
774
+        GeneNumID = cnn,
775
+        stringsAsFactors = FALSE)
776
+
590 777
     }
591
-    
778
+
592 779
     if( (length(long.rt) + length(long.epistasis) + length(long.orderEffects) +
593
-             nrow(geneNoInt) + nrow(genotFitness)) == 0)
594
-        stop("You have specified nothing!")
780
+         nrow(geneNoInt) + nrow(genotFitness)) == 0)
781
+      stop("You have specified nothing!")
595 782
 
596 783
     if(calledBy == "allFitnessEffects") {
597
-        if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
598
-            graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
599
-        } else {
600
-            graphE <- NULL
601
-        }
602
-    } else {
784
+      if((length(long.rt) + length(long.epistasis) + length(long.orderEffects)) > 1) {
785
+        graphE <- fitnessEffectsToIgraph(rT, epistasis, orderEffects)
786
+      } else {
603 787
         graphE <- NULL
788
+      }
789
+    } else {
790
+      graphE <- NULL
604 791
     }
605 792
     if(!is.null(drvNames)) {
606
-        drv <- unique(getGeneIDNum(geneModule, geneNoInt, fitnessLandscape_gene_id,
607
-                                   drvNames))
608
-        ## drivers should never be in the geneNoInt; Why!!!???
609
-        ## Catch the problem. This is an overkill,
610
-        ## so since we catch the issue, we could leave the geneNoInt. But
611
-        ## that should not be there in this call.
612
-        ## if(any(drvNames %in% geneNoInt$Gene)) {
613
-        ##     stop(paste("At least one gene in drvNames is a geneNoInt gene.",
614
-        ##                "That is not allowed.",
615
-        ##                "If that gene is a driver, pass it as gene in the epistasis",
616
-        ##                "component."))
617
-        ## }
618
-        ## drv <- getGeneIDNum(geneModule, NULL, drvNames)
793
+      drv <- unique(getGeneIDNum(geneModule, geneNoInt, fitnessLandscape_gene_id,
794
+                                 drvNames))
795
+      ## drivers should never be in the geneNoInt; Why!!!???
796
+      ## Catch the problem. This is an overkill,
797
+      ## so since we catch the issue, we could leave the geneNoInt. But
798
+      ## that should not be there in this call.
799
+      ## if(any(drvNames %in% geneNoInt$Gene)) {
800
+      ##     stop(paste("At least one gene in drvNames is a geneNoInt gene.",
801
+      ##                "That is not allowed.",
802
+      ##                "If that gene is a driver, pass it as gene in the epistasis",
803
+      ##                "component."))
804
+      ## }
805
+      ## drv <- getGeneIDNum(geneModule, NULL, drvNames)
619 806
     } else {
620
-        ## we used to have this default
621
-        ## drv <- geneModule$GeneNumID[-1]
622
-        drv <- vector(mode = "integer", length = 0L)
807
+      ## we used to have this default
808
+      ## drv <- geneModule$GeneNumID[-1]
809
+      drv <- vector(mode = "integer", length = 0L)
623 810
     }
624
-  
811
+
625 812
     if(!keepInput) {
626
-        rT <- epistasis <- orderEffects <- noIntGenes <- NULL
813
+      rT <- epistasis <- orderEffects <- noIntGenes <- NULL
627 814
     }
628 815
 
629 816
     out <- list(long.rt = long.rt,
... ...
@@ -641,14 +828,114 @@ allFitnessORMutatorEffects <- function(rT = NULL,
641 828
                 noIntGenes = noIntGenes,
642 829
                 fitnessLandscape = genotFitness,
643 830
                 fitnessLandscape_df = fitnessLandscape_df,
644
-                fitnessLandscape_gene_id = fitnessLandscape_gene_id
645
-                )
831
+                fitnessLandscape_gene_id = fitnessLandscape_gene_id,
832
+                fitnessLandscapeVariables = vector(mode = "character", length = 0L),
833
+                frequencyDependentFitness = frequencyDependentFitness,
834
+                frequencyType = frequencyType)
835
+                #spPopSizes = vector(mode = "integer", length = 0L)
836
+    
646 837
     if(calledBy == "allFitnessEffects") {
647
-        class(out) <- c("fitnessEffects")
838
+      class(out) <- c("fitnessEffects")
648 839
     } else if(calledBy == "allMutatorEffects") {
649
-        class(out) <- c("mutatorEffects")
840
+      class(out) <- c("mutatorEffects")
650 841
     }
651
-    return(out)
842
+  } else {
843
+
844
+    if(is.null(genotFitness)) {
845
+      #genotFitness <- matrix(NA, nrow = 0, ncol = 1)
846
+      #fitnessLandscape_df <- data.frame()
847
+      #fitnessLandscape_gene_id <- data.frame()
848
+      stop("You have a null genotFitness in a frequency dependent fitness situation.")
849
+    } else {
850
+      cnn <- 1:(ncol(genotFitness) - 1)
851
+      gfn <- apply(genotFitness[, -ncol(genotFitness), drop = FALSE], 1,
852
+                   function(x) paste(cnn[as.logical(x)],
853
+                                     collapse = ", "))
854
+      ## rownames(genotFitness) <- gfn
855
+      fitnessLandscape_df <-
856
+        data.frame(Genotype = gfn,
857
+                   Fitness = genotFitness[, ncol(genotFitness)],
858
+                   stringsAsFactors = FALSE)
859
+
860
+      attr(fitnessLandscape_df,'row.names') <-
861
+        as.integer(attr(fitnessLandscape_df,'row.names'))
862
+
863
+      fitnessLandscape_gene_id <- data.frame(
864
+        Gene = colnames(genotFitness)[-ncol(genotFitness)],
865
+        GeneNumID = cnn,
866
+        stringsAsFactors = FALSE)
867
+
868
+      if(frequencyType == "auto"){
869
+        ch <- paste(as.character(fitnessLandscape_df[, ncol(fitnessLandscape_df)]), collapse = "")
870
+        #print(ch)
871
+        if( grepl("f_", ch, fixed = TRUE) ){
872
+          frequencyType = "rel"
873
+        } else{
874
+          frequencyType = "abs"
875
+        }
876
+      } else { frequencyType = frequencyType }
877
+      ## Wrong: assumes all genotypes in fitness landscape
878
+      ## fitnessLandscapeVariables = fVariablesN(ncol(genotFitness) - 1, frequencyType)
879
+      stopifnot(identical(genotFitness$Fitness, fitnessLandscape_df$Fitness))
880
+      flvars_and_fitvars <- create_flvars_fitvars(genotFitness, frequencyType)
881
+      fitnessLandscapeVariables <- flvars_and_fitvars$flvars
882
+      Fitness_as_fvars <- flvars_and_fitvars$Fitness_as_fvars
883
+    }
884
+
885
+    if(!is.null(drvNames)) {
886
+      drv <- unique(getGeneIDNum(geneModule, geneNoInt, fitnessLandscape_gene_id,
887
+                                 drvNames))
888
+    } else {
889
+      drv <- vector(mode = "integer", length = 0L)
890
+    }
891
+    defaultGeneModuleDF <- data.frame(Gene = "Root",
892
+                                      Module = "Root",
893
+                                      GeneNumID = 0,
894
+                                      ModuleNumID = 0,
895
+                                      stringsAsFactors = FALSE)
896
+      ## Create, for the user, a single data frame with everything.
897
+      ## This is what C++ should consume
898
+
899
+   
900
+     
901
+      ## This ought to allow to pass fitness spec as letters. Preserve original
902
+      Fitness_original_as_letters <- fitnessLandscape_df$Fitness
903
+      fitnessLandscape_df$Fitness <- Fitness_as_fvars
904
+      
905
+      full_FDF_spec <- cbind(genotFitness[, -ncol(genotFitness)]
906
+                 , Genotype_letters = genotype_letterslabel(genotFitness[, -ncol(genotFitness)])
907
+                 , Genotype_fvars = fitnessLandscapeVariables ## used in C++
908
+                 , Fitness_as_fvars = Fitness_as_fvars
909
+                 , Fitness_as_letters = Fitness_original_as_letters
910
+                   )
911
+      rownames(full_FDF_spec) <- 1:nrow(full_FDF_spec)
912
+      
913
+    out <- list(long.rt = list(),
914
+                long.epistasis = list(),
915
+                long.orderEffects = list(),
916
+                long.geneNoInt = data.frame(),
917
+                geneModule = defaultGeneModuleDF, ##Trick to pass countGenesFe>2,
918
+                gMOneToOne = TRUE,
919
+                geneToModule = c(Root = "Root"),
920
+                graph = list(),
921
+                drv = drv,
922
+                rT = NULL,
923
+                epistasis = NULL,
924
+                orderEffects = NULL,
925
+                noIntGenes = NULL,
926
+                fitnessLandscape = genotFitness,
927
+                fitnessLandscape_df = fitnessLandscape_df,
928
+                fitnessLandscape_gene_id = fitnessLandscape_gene_id,
929
+                fitnessLandscapeVariables = fitnessLandscapeVariables,
930
+                frequencyDependentFitness = frequencyDependentFitness,
931
+                frequencyType = frequencyType,
932
+                full_FDF_spec = full_FDF_spec
933
+                #spPopSizes = spPopSizes
934
+              )
935
+
936
+    class(out) <- c("fitnessEffects")
937
+  }
938
+  return(out)
652 939
 }
653 940
 
654 941
 ## Former version, with fitness landscape
... ...
@@ -663,13 +950,13 @@ allFitnessORMutatorEffects <- function(rT = NULL,
663 950
 ##                                        calledBy = NULL) {
664 951
 ##     ## From allFitnessEffects. Generalized so we deal with Fitness
665 952
 ##     ## and mutator.
666
-    
953
+
667 954
 ##     ## restrictions: the usual rt
668 955
 
669 956
 ##     ## epistasis: as it says, with the ":"
670 957
 
671 958
 ##     ## orderEffects: the ">"
672
-    
959
+
673 960
 ##     ## All of the above can be genes or can be modules (if you pass a
674 961
 ##     ## geneToModule)
675 962
 
... ...
@@ -685,14 +972,14 @@ allFitnessORMutatorEffects <- function(rT = NULL,
685 972
 ##     ## called appropriately?
686 973
 ##     if( !(calledBy %in% c("allFitnessEffects", "allMutatorEffects") ))
687 974
 ##         stop("How did you call this function?. Bug.")
688
-    
975
+
689 976
 ##     if(calledBy == "allMutatorEffects") {
690 977
 ##         ## very paranoid check
691 978
 ##         if( !is.null(rT) || !is.null(orderEffects) || !is.null(drvNames))
692 979
 ##             stop("allMutatorEffects called with forbidden arguments.",
693 980
 ##                  "Is this an attempt to subvert the function?")
694 981
 ##     }
695
-    
982
+
696 983
 ##     rtNames <- NULL
697 984
 ##     epiNames <- NULL
698 985
 ##     orNames <- NULL
... ...
@@ -735,7 +1022,7 @@ allFitnessORMutatorEffects <- function(rT = NULL,
735 1022
 ##                        "or order effects not in geneToModule"))
736 1023
 ##     }
737 1024
 ##     geneModule <- gm.to.geneModuleL(geneToModule, one.to.one = gMOneToOne)
738
-    
1025
+
739 1026
 ##     idm <- unique(geneModule$ModuleNumID)
740 1027
 ##     names(idm) <- unique(geneModule$Module)
741 1028
 
... ...
@@ -761,7 +1048,7 @@ allFitnessORMutatorEffects <- function(rT = NULL,
761 1048
 ##                                                              sort = FALSE,
762 1049
 ##                                                              sign = FALSE))
763 1050
 ##     }
764
-    
1051
+
765 1052
 ##     if(!is.null(noIntGenes)) {
766 1053
 ##         if(inherits(noIntGenes, "character")) {
767 1054
 ##             wm <- paste("noIntGenes is a character vector.",
... ...
@@ -772,7 +1059,7 @@ allFitnessORMutatorEffects <- function(rT = NULL,
772 1059
 ##                         "We are stopping.")
773 1060
 ##             stop(wm)
774 1061
 ##         }
775
-            
1062
+
776 1063
 ##         mg <- max(geneModule[, "GeneNumID"])
777 1064
 ##         gnum <- seq_along(noIntGenes) + mg
778 1065
 ##         if(!is.null(names(noIntGenes))) {
... ...
@@ -829,7 +1116,7 @@ allFitnessORMutatorEffects <- function(rT = NULL,
829 1116
 ##         ## drv <- geneModule$GeneNumID[-1]
830 1117
 ##         drv <- vector(mode = "integer", length = 0L)
831 1118
 ##     }
832
-    
1119
+
833 1120
 ##     if(!keepInput) {
834 1121
 ##         rT <- epistasis <- orderEffects <- noIntGenes <- NULL
835 1122
 ##     }
... ...
@@ -845,7 +1132,7 @@ allFitnessORMutatorEffects <- function(rT = NULL,
845 1132
 ##                 rT = rT,
846 1133
 ##                 epistasis = epistasis,
847 1134
 ##                 orderEffects = orderEffects,
848
-##                 noIntGenes = noIntGenes                
1135
+##                 noIntGenes = noIntGenes
849 1136
 ##                 )
850 1137
 ##     if(calledBy == "allFitnessEffects") {
851 1138
 ##         class(out) <- c("fitnessEffects")
... ...
@@ -855,7 +1142,6 @@ allFitnessORMutatorEffects <- function(rT = NULL,
855 1142
 ##     return(out)
856 1143
 ## }
857 1144
 
858
-
859 1145
 allFitnessEffects <- function(rT = NULL,
860 1146
                               epistasis = NULL,
861 1147
                               orderEffects = NULL,
... ...
@@ -863,24 +1149,68 @@ allFitnessEffects <- function(rT = NULL,
863 1149
                               geneToModule = NULL,
864 1150
                               drvNames = NULL,
865 1151
                               genotFitness = NULL,
866
-                              keepInput = TRUE) {
1152
+                              keepInput = TRUE,
1153
+                              frequencyDependentFitness = FALSE,
1154
+                              frequencyType = NA) {
1155
+                              #spPopSizes = NULL) {
867 1156
 
868
-    if(!is.null(genotFitness)) {
869
-        if(!is.null(rT) || !is.null(epistasis) ||
870
-           !is.null(orderEffects) || !is.null(noIntGenes) ||
871
-           !is.null(geneToModule)) {
872
-            stop("You have a non-null genotFitness.",
873
-                 " If you pass the complete genotype to fitness mapping",
874
-                 " you cannot pass any of rT, epistasis, orderEffects",
875
-                 " noIntGenes or geneToModule.")
1157
+    if(!frequencyDependentFitness) {
1158
+        
1159
+        if(!is.na(frequencyType)){
1160
+            warning("frequencyType set to NA")
876 1161
         }
1162
+        ## this is a kludge, but we must pass something not NA and not NULL
1163
+        ## to the C++ code
1164
+        frequencyType = "freq_dep_not_used"
877 1165
 
878
-        genotFitness_std <- to_genotFitness_std(genotFitness, simplify = TRUE)
879
-        ## epistasis <- from_genotype_fitness(genotFitness)
1166
+    if(!is.null(genotFitness)) {
1167
+      if(!is.null(rT) || !is.null(epistasis) ||
1168
+         !is.null(orderEffects) || !is.null(noIntGenes) ||
1169
+         !is.null(geneToModule)) {
1170
+        stop("You have a non-null genotFitness.",
1171
+             " If you pass the complete genotype to fitness mapping",
1172
+             " you cannot pass any of rT, epistasis, orderEffects",
1173
+             " noIntGenes or geneToModule.")
1174
+      }
1175
+
1176
+      genotFitness_std <- to_genotFitness_std(genotFitness,
1177
+                                              frequencyDependentFitness = FALSE,
1178
+                                              frequencyType = frequencyType,
1179
+                                              simplify = TRUE)
1180
+      ## epistasis <- from_genotype_fitness(genotFitness)
880 1181
     } else {
881
-        genotFitness_std <- NULL
1182
+      genotFitness_std <- NULL
882 1183
     }
883 1184
     allFitnessORMutatorEffects(
1185
+      rT = rT,
1186
+      epistasis = epistasis,
1187
+      orderEffects = orderEffects,
1188
+      noIntGenes = noIntGenes,
1189
+      geneToModule = geneToModule,
1190
+      drvNames = drvNames,
1191
+      keepInput = keepInput,
1192
+      genotFitness = genotFitness_std,
1193
+      calledBy = "allFitnessEffects",
1194
+      frequencyDependentFitness = FALSE,
1195
+      frequencyType = frequencyType)
1196
+      #spPopSizes = spPopSizes)
1197
+
1198
+  } else {
1199
+
1200
+    if(!(frequencyType %in% c('abs', 'rel', 'auto'))){
1201
+      #set frequencyType = "auto" in case you did not specify 'rel' or 'abs'
1202
+      frequencyType = "auto"
1203
+      message("frequencyType set to 'auto'")
1204
+    }
1205
+
1206
+    if(is.null(genotFitness)) {
1207
+      stop("You have a null genotFitness in a frequency dependent fitness situation.")
1208
+    } else {
1209
+      genotFitness_std <- to_genotFitness_std(genotFitness,
1210
+                                              frequencyDependentFitness = TRUE,
1211
+                                              frequencyType = frequencyType,
1212
+                                              simplify = TRUE)
1213
+      allFitnessORMutatorEffects(
884 1214
         rT = rT,
885 1215
         epistasis = epistasis,
886 1216
         orderEffects = orderEffects,
... ...
@@ -889,7 +1219,12 @@ allFitnessEffects <- function(rT = NULL,
889 1219
         drvNames = drvNames,
890 1220
         keepInput = keepInput,
891 1221
         genotFitness = genotFitness_std,
892
-        calledBy = "allFitnessEffects")
1222
+        calledBy = "allFitnessEffects",
1223
+        frequencyDependentFitness = TRUE,
1224
+        frequencyType = frequencyType)
1225
+        #spPopSizes = spPopSizes)
1226
+    }
1227
+  }
893 1228
 }
894 1229
 
895 1230
 ## Former version
... ...
@@ -936,7 +1271,7 @@ allFitnessEffects <- function(rT = NULL,
936 1271
 ##     ## epistasis: as it says, with the ":"
937 1272
 
938 1273
 ##     ## orderEffects: the ">"
939
-    
1274
+
940 1275
 ##     ## All of the above can be genes or can be modules (if you pass a
941 1276
 ##     ## geneToModule)
942 1277
 
... ...
@@ -949,7 +1284,7 @@ allFitnessEffects <- function(rT = NULL,
949 1284
 ##     ## same in epistasis and order, but we would be splitting twice
950 1285
 ##     ## (whereas for rT extracting the names is very simple).
951 1286
 
952
-    
1287
+
953 1288
 ##     rtNames <- NULL
954 1289
 ##     epiNames <- NULL
955 1290
 ##     orNames <- NULL
... ...
@@ -969,7 +1304,7 @@ allFitnessEffects <- function(rT = NULL,
969 1304
 ##                           function(x) lapply(x$ids,
970 1305
 ##                                              function(z) strsplit(z, "^-"))))),
971 1306
 ##                             "")
972
-        
1307
+
973 1308
 ##     } else {
974 1309
 ##         long.epistasis <- list()
975 1310
 ##     }
... ...
@@ -1019,7 +1354,7 @@ allFitnessEffects <- function(rT = NULL,
1019 1354
 ##                                                              sort = FALSE,
1020 1355
 ##                                                              sign = FALSE))
1021 1356
 ##     }
1022
-    
1357
+
1023 1358
 ##     if(!is.null(noIntGenes)) {
1024 1359
 ##         mg <- max(geneModule[, "GeneNumID"])
1025 1360
 ##         gnum <- seq_along(noIntGenes) + mg
... ...
@@ -1067,7 +1402,7 @@ allFitnessEffects <- function(rT = NULL,
1067 1402
 ##                 rT = rT,
1068 1403
 ##                 epistasis = epistasis,
1069 1404
 ##                 orderEffects = orderEffects,
1070
-##                 noIntGenes = noIntGenes                
1405
+##                 noIntGenes = noIntGenes
1071 1406
 ##                 )
1072 1407
 ##     class(out) <- c("fitnessEffects")
1073 1408
 ##     return(out)
... ...
@@ -1078,7 +1413,7 @@ allFitnessEffects <- function(rT = NULL,
1078 1413
 ## rtAndGeneModule <- function(mdeps, gM = NULL) {
1079 1414
 ##     ## To show a table of restrictions when there are modules. Do not use
1080 1415
 ##     ## for anything else. Maybe as intermediate to plotting.
1081
-    
1416
+
1082 1417
 ##     ## Specify restriction table of modules and a mapping of modules to
1083 1418
 ##     ## genes. gM is a named vector; names are modules, values are elements
1084 1419
 ##     ## of each module.
... ...
@@ -1099,7 +1434,7 @@ allFitnessEffects <- function(rT = NULL,
1099 1434
 ##             stop("Some values in child not from a known module")
1100 1435
 ##         if(any(is.na(match(names(gM), c(mdeps[, 1], mdeps[, 2])))))
1101 1436
 ##             stop("Some values in module in neither parent or child")
1102
-        
1437
+
1103 1438
 ##         parent <- gM[mdeps[, 1]]
1104 1439
 ##         child <- gM[mdeps[, 2]]
1105 1440
 ##         df <- data.frame(parent = parent,
... ...
@@ -1139,10 +1474,12 @@ allFitnessEffects <- function(rT = NULL,
1139 1474
 
1140 1475
 evalGenotypeORMut <- function(genotype,
1141 1476
                               fmEffects,
1477
+                              spPopSizes = spPopSizes,
1142 1478
                               verbose = FALSE,
1143 1479
                               echo = FALSE,
1144 1480
                               model = "",
1145
-                              calledBy_= NULL) {
1481
+                              calledBy_= NULL,
1482
+                              currentTime = currentTime) {
1146 1483
     ## genotype can be a vector of integers, that are the exact same in
1147 1484
     ## the table of fmEffects or a vector of strings, or a vector (a
1148 1485
     ## string) with genes separated by "," or ">"
... ...
@@ -1159,12 +1496,20 @@ evalGenotypeORMut <- function(genotype,
1159 1496
     if( (model %in% c("Bozic", "bozic1", "bozic2")) &&
1160 1497
         (nrow(fmEffects$fitnessLandscape_df) > 0)) {
1161 1498
         warning("Bozic model passing a fitness landscape will not work",
1162
-                    " for now.")
1499
+                " for now.")
1163 1500
     }
1164
-    
1501
+
1502
+    ## This will avoid errors is evalRGenotype where spPopSizes = NULL  
1503
+    if (!fmEffects$frequencyDependentFitness) {
1504
+        spPopSizes = 0
1505
+    } else {
1506
+        spPopSizes <- match_spPopSizes(spPopSizes, fmEffects)
1507
+    }
1508
+
1165 1509
     if(echo)
1166 1510
         cat(paste("Genotype: ", genotype))
1167
-    if(!is.integer(genotype)) {
1511
+
1512
+    if(is.character(genotype)) {
1168 1513
         if(length(grep(">", genotype))) {
1169 1514
             genotype <- nice.vector.eo(genotype, ">")
1170 1515
         } else if(length(grep(",", genotype))) {
... ...
@@ -1176,27 +1521,76 @@ evalGenotypeORMut <- function(genotype,
1176 1521
         all.g.names <- c(fmEffects$geneModule$Gene,
1177 1522
                          fmEffects$long.geneNoInt$Gene,
1178 1523
                          fmEffects$fitnessLandscape_gene_id$Gene)
1179
-        genotype <- all.g.nums[match(genotype, all.g.names)]
1524
+        if((length(genotype) == 1) && (genotype %in% c("", "WT") )) {
1525
+            genotype <- 0
1526
+        } else {
1527
+            genotype <- all.g.nums[match(genotype, all.g.names)]
1528
+        }
1529
+    } else {
1530
+        all.g.nums <- c(fmEffects$geneModule$GeneNumID,
1531
+                        fmEffects$long.geneNoInt$GeneNumID,
1532
+                        fmEffects$fitnessLandscape_gene_id$GeneNumID)
1533
+        
1534
+        if(!all(sapply(genotype,  function(x) x %in% all.g.nums))){
1535
+            stop("Genotype as vector of numbers contains genes not in fitnessEffects/mutatorEffects.")
1536
+        }
1180 1537
     }
1181
-    if(any(is.na(genotype)))
1182
-        stop("genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1183
-    if(!length(genotype))
1184
-        stop("genotypes must have at least one mutated gene")
1185
-    if(any(genotype < 0)) {
1186
-        stop(paste("genotypes cannot contain negative values.",
1187
-                   "If you see this message, you found a bug."))
1538
+
1539
+    ## FIXME: About the WT: in fmEffects$geneModule$Gene we have "Root". We might
1540
+    ## want to refactor that and turn all those to WT. It would avoid the need for
1541
+    ## the if above to catch the WT or ""
1542
+    
1543
+    if(!fmEffects$frequencyDependentFitness){
1544
+        
1545
+        if( any(is.na(genotype)) ){
1546
+            stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1547
+        }
1548
+        
1549
+        if((!length(genotype))){
1550
+            stop("Genotypes must have at least one mutated gene")
1551
+        }
1552
+        if(any(genotype < 0)) {
1553
+            stop(paste("Genotypes cannot contain negative values.",
1554
+                       "If you see this message, you found a bug."))
1555
+        }
1556
+        if(length(genotype) == 1 && genotype == 0){
1557
+            ## We do not handle WT for now
1558
+            stop("Genotype cannot be 0. We do not handle WT on its own in non-freq-dep: its fitness is 1 by decree.")
1559
+        }
1560
+        
1561
+        if(any(genotype == 0)){
1562
+            stop("Genotype cannot contain any 0.")
1563
+        }
1564
+        
1565
+    } else {
1566
+        if(length(genotype) == 1 && is.na(genotype)){
1567
+            stop("Genotype contains NA or a gene not in fitnessEffects/mutatorEffects")
1568
+        } else if(length(genotype) == 1 && genotype == 0) {
1569
+            ## for WT
1570
+            genotype <- vector(mode = "integer", length = 0L)
1571
+        } else if(length(genotype) > 1){
1572
+            if( any(is.na(genotype)) ){
1573
+                stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1574
+            }
1575
+            if(any(genotype == 0)){
1576
+                stop("Genotype cannot contain any 0 if its length > 1")
1577
+            }
1578
+        }
1188 1579
     }
1580
+    
1189 1581
     if(model %in% c("Bozic", "bozic1", "bozic2") )
1190 1582
         prodNeg <- TRUE
1191 1583
     else
1192 1584
         prodNeg <- FALSE
1585
+
1193 1586
     ff <- evalRGenotype(rG = genotype,
1194 1587
                         rFE = fmEffects,
1588
+                        spPop = spPopSizes,
1195 1589
                         verbose = verbose,
1196 1590
                         prodNeg = prodNeg,
1197
-                        calledBy_ = calledBy_)
1591
+                        calledBy_ = calledBy_,
1592
+                        currentTime = currentTime)
1198 1593
 
1199
-    
1200 1594
     if(echo) {
1201 1595
         if(calledBy_ == "evalGenotype") {
1202 1596
             if(!prodNeg)
... ...
@@ -1206,35 +1600,53 @@ evalGenotypeORMut <- function(genotype,
1206 1600
         } else if(calledBy_ == "evalGenotypeMut") {
1207 1601
             cat(" Mutation rate product :", ff, "\n")
1208 1602
         }
1209
-        
1210
-    } 
1603
+
1604
+    }
1605
+
1211 1606
     return(ff)
1212 1607
 }
1213 1608
 
1214
-evalGenotype <- function(genotype, fitnessEffects,
1609
+evalGenotype <- function(genotype,
1610
+                         fitnessEffects,
1611
+                         spPopSizes = NULL,
1215 1612
                          verbose = FALSE,
1216 1613
                          echo = FALSE,
1217
-                         model = "") {
1614
+                         model = "",
1615
+                         currentTime = 0
1616
+                         ) {
1617
+  
1218 1618
     if(inherits(fitnessEffects, "mutatorEffects"))
1219 1619
          stop("You are trying to get the fitness of a mutator specification. ",
1220 1620
              "You did not pass an object of class fitnessEffects.")
1221 1621
 
1622
+   if (fitnessEffects$frequencyDependentFitness) {
1623
+      if (is.null(spPopSizes))
1624
+      stop("You have a NULL spPopSizes")
1625
+    if (!(length(spPopSizes) == nrow(fitnessEffects$fitnessLandscape)))
1626
+      stop("spPopSizes must be as long as number of genotypes")
1627
+   }
1628
+
1629
+
1222 1630
     evalGenotypeORMut(genotype = genotype,
1223
-                       fmEffects = fitnessEffects,
1224
-                       verbose = verbose,
1225
-                       echo = echo,
1226
-                       model  = model ,
1227
-                       calledBy_= "evalGenotype"
1228
-                       )
1631
+                      fmEffects = fitnessEffects,
1632
+                      spPopSizes = spPopSizes,
1633
+                      verbose = verbose,
1634
+                      echo = echo,
1635
+                      model  = model ,
1636
+                      calledBy_= "evalGenotype",
1637
+                      currentTime = currentTime)
1229 1638
 }
1230 1639
 
1231 1640
 
1232 1641
 evalGenotypeFitAndMut <- function(genotype,
1233 1642
                                   fitnessEffects,
1234 1643
                                   mutatorEffects,
1644
+                                  spPopSizes = NULL,
1235 1645
                                   verbose = FALSE,
1236 1646
                                   echo = FALSE,
1237
-                                  model = "") {
1647
+                                  model = "",
1648
+                                  currentTime = 0
1649
+                                  ) {
1238 1650
     
1239 1651
     ## Must deal with objects from previous, pre flfast, modifications
1240 1652
     if(!exists("fitnessLandscape_gene_id", where = fitnessEffects)) {
... ...
@@ -1246,11 +1658,26 @@ evalGenotypeFitAndMut <- function(genotype,
1246 1658
         warning("Bozic model passing a fitness landscape will not work",
1247 1659
                     " for now.")
1248 1660
     }
1661
+  
1662
+    if(fitnessEffects$frequencyDependentFitness) {
1663
+      if (is.null(spPopSizes))
1664
+        stop("You have a NULL spPopSizes")
1665
+      if (!(length(spPopSizes) == nrow(fitnessEffects$fitnessLandscape)))
1666
+          stop("spPopSizes must be as long as number of genotypes")
1667
+      spPopSizes <- match_spPopSizes(spPopSizes, fitnessEffects)
1668
+    }
1669
+  
1670
+    # This will avoid errors is evalRGenotype where spPopSizes = NULL  
1671
+    if (!fitnessEffects$frequencyDependentFitness) {
1672
+      spPopSizes = 0
1673
+    }
1674
+  
1249 1675
     prodNeg <- FALSE
1250 1676
     ## Next is from evalGenotypeAndMut
1251 1677
     if(echo)
1252 1678
         cat(paste("Genotype: ", genotype))
1253
-    if(!is.integer(genotype)) {
1679
+    
1680
+    if(is.character(genotype)) {
1254 1681
         if(length(grep(">", genotype))) {
1255 1682
             genotype <- nice.vector.eo(genotype, ">")
1256 1683
         } else if(length(grep(",", genotype))) {
... ...
@@ -1263,15 +1690,51 @@ evalGenotypeFitAndMut <- function(genotype,
1263 1690
                          fitnessEffects$long.geneNoInt$Gene,
1264 1691
                          fitnessEffects$fitnessLandscape_gene_id$Gene)
1265 1692
         genotype <- all.g.nums[match(genotype, all.g.names)]
1693
+        
1694
+    } else {
1695
+      all.g.nums <- c(fitnessEffects$geneModule$GeneNumID,
1696
+                      fitnessEffects$long.geneNoInt$GeneNumID,
1697
+                      fitnessEffects$fitnessLandscape_gene_id$GeneNumID)
1698
+      if(!all(sapply(genotype,  function(x)x %in% all.g.nums))){
1699
+        stop("Genotype as vector of numbers contains genes not in fitnessEffects/mutatorEffects.")
1700
+      }
1701
+    }
1702
+    
1703
+    if(!fitnessEffects$frequencyDependentFitness){
1704
+      
1705
+      if( any(is.na(genotype)) ){
1706
+        stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1707
+      }
1708
+      
1709
+      if((!length(genotype))){
1710
+        stop("Genotypes must have at least one mutated gene")
1266 1711
     }
1267
-    if(any(is.na(genotype)))
1268
-        stop("genotype contains NAs or genes not in fitnessEffects")
1269
-    if(!length(genotype))
1270
-        stop("genotypes must have at least one mutated gene")
1271 1712
     if(any(genotype < 0)) {
1272
-        stop(paste("genotypes cannot contain negative values.",
1713
+        stop(paste("Genotypes cannot contain negative values.",
1273 1714
                    "If you see this message, you found a bug."))
1274 1715
     }
1716
+      if(length(genotype) == 1 && genotype == 0){
1717
+        stop("Genotype cannot be 0.")
1718
+      }
1719
+      
1720
+      if(any(genotype == 0)){
1721
+        stop("Genotype cannot contain any 0.")
1722
+      }
1723
+      
1724
+    }else{
1725
+      if(length(genotype) == 1 && is.na(genotype)){
1726
+        stop("Genotype contains NA or a gene not in fitnessEffects/mutatorEffects")
1727
+      }else if(length(genotype) == 1 && genotype == 0){
1728
+        genotype <- vector(mode = "integer", length = 0L)
1729
+      }else if(length(genotype) > 1){
1730
+        if( any(is.na(genotype)) ){
1731
+          stop("Genotype contains NAs or genes not in fitnessEffects/mutatorEffects")
1732
+        }
1733
+        if(any(genotype == 0)){
1734
+          stop("Genotype cannot contain any 0 if its length > 1")
1735
+        }
1736
+      }
1737
+    }
1275 1738
 
1276 1739
     full2mutator_ <- matchGeneIDs(mutatorEffects,
1277 1740
                                   fitnessEffects)$Reduced
... ...
@@ -1279,10 +1742,14 @@ evalGenotypeFitAndMut <- function(genotype,
1279 1742
         prodNeg <- TRUE
1280 1743
     else
1281 1744
         prodNeg <- FALSE
1282
-    evalRGenotypeAndMut(genotype, fitnessEffects,
1283
-                        mutatorEffects, full2mutator_,
1745
+    evalRGenotypeAndMut(genotype,
1746
+                        fitnessEffects,
1747
+                        mutatorEffects,
1748
+                        spPopSizes,
1749
+                        full2mutator_,
1284 1750
                         verbose = verbose,
1285
-                        prodNeg = prodNeg)
1751
+                        prodNeg = prodNeg,
1752
+                        currentTime = currentTime)
1286 1753
 }
1287 1754
 
1288 1755
 ## evalGenotype <- function(genotype, fitnessEffects,
... ...
@@ -1341,15 +1808,53 @@ evalGenotypeFitAndMut <- function(genotype,
1341 1808
 
1342 1809
 ## I am here: simplify this
1343 1810
 
1811
+## fitnesEffects from FDF, spPopSizes -> reordered spPopSizes
1812
+##    Verify if named. If not, issue a warning.
1813
+##    If yes, check matches genotypes and reorder
1814
+match_spPopSizes <- function(sp, fe){
1815
+    if(!fe$frequencyDependentFitness) return(sp)
1816
+    if(is.null(names(sp))) {
1817
+        warning("spPopSizes unnamed: cannot check genotype names.")
1818
+        return(sp)
1819
+    }
1820
+    nns <- names(sp)
1821
+    nns <- sort_genes_genots(nns)
1822
+    names(sp) <- nns
1823
+    nnf <- fe$full_FDF_spec$Genotype_letters
1824
+    
1825
+    if( ( length(nns) != length(nnf) ) ||
1826
+        (!(identical(sort(nnf), sort(nns))) ) )
1827
+        stop("Genotype names in spPopSizes differ w.r.t. fitness landscape")
1828
+
1829
+    return(sp[nnf])
1830
+}
1831
+
1832
+## genotype names -> genotype names with names sorted
1833
+sort_genes_genots <- function(x) {
1834
+    return(
1835
+        unlist(
1836
+            lapply(
1837
+                lapply(strsplit(x, ", "), sort),
1838
+                function(x) paste(x, collapse = ", "))))
1839
+    ## Or this:
1840
+    ## lapply(names(nn3), function(x) paste(sort(strsplit(x, ", ")[[1]]),
1841
+    ## collapse = ", "))
1842
+}
1843
+
1844
+
1845
+        
1344 1846
 evalAllGenotypesORMut <- function(fmEffects,
1345
-                                  order = FALSE, max = 256,
1346
-                             addwt = FALSE,
1347
-                             model = "",
1348
-                             calledBy_ = "") {
1349
-##                             minimal = FALSE) {
1847
+                                  order = FALSE, 
1848
+                                  max = 256,
1849
+                                  addwt = FALSE,
1850
+                                  model = "",