Browse code

2.17.7: no longer depends on nem; more models from MAGELLAN

ramon diaz-uriarte (at Phelsuma) authored on 30/01/2020 13:39:14
Showing11 changed files

... ...
@@ -1,13 +1,45 @@
1 1
 Package: OncoSimulR
2 2
 Type: Package
3 3
 Title: Forward Genetic Simulation of Cancer Progression with Epistasis 
4
-Version: 2.17.5
5
-Date: 2019-12-19
4
+Version: 2.17.7
5
+Date: 2020-01-30
6 6
 Authors@R: c(person("Ramon", "Diaz-Uriarte", role = c("aut", "cre"),
7 7
 		     email = "rdiaz02@gmail.com"),
8
-	      person("Mark", "Taylor", role = "ctb", email = "ningkiling@gmail.com"))
8
+	      person("Mark", "Taylor", role = "ctb", email = "ningkiling@gmail.com"),
9
+	      person("Guillermo", "Gorines Cordero", role = "ctb"),
10
+	      person("Ivan", "Lorca Alonso", role = "ctb"),
11
+	      person("Francisco", "Mu~noz Lopez", role = "ctb"),
12
+	      person("David", "Roncero Moro~no", role = "ctb"),
13
+	      person("Alvaro" , "Quevedo" , role = "ctb"),
14
+	      person("Pablo" ,"Perez" , role = "ctb"),
15
+	      person("Cristina" ,"Devesa" , role = "ctb"),
16
+	      person("Alejandro", "Herrador" , role = "ctb"),
17
+	      person("Holger", "Froehlich", role = "ctb"),
18
+	      person("Florian", "Markowetz", role = "ctb"),
19
+	      person("Achim", "Tresch", role = "ctb"),
20
+	      person("Theresa", "Niederberger", role = "ctb"),
21
+	      person("Christian", "Bender", role = "ctb"),
22
+	      person("Matthias", "Maneck", role = "ctb"),
23
+	      person("Claudio", "Lottaz", role = "ctb"),
24
+	      person("Tim", "Beissbarth", role = "ctb"))
9 25
 Author: Ramon Diaz-Uriarte [aut, cre],
10
-	Mark Taylor [ctb]
26
+	Mark Taylor [ctb],
27
+	Guillermo Gorines Cordero [ctb],
28
+	Ivan Lorca Alonso [ctb],
29
+	Francisco Mu\~noz Lopez [ctb],
30
+	David Roncero Moro\~no [ctb],
31
+	Alvaro Quevedo [ctb],
32
+	Pablo Perez [ctb],
33
+	Cristina Devesa [ctb],
34
+	Alejandro Herrador [ctb],
35
+	Holger Froehlich [ctb],
36
+	Florian Markowetz [ctb],
37
+	Achim Tresch [ctb],
38
+	Theresa Niederberger [ctb],
39
+	Christian Bender [ctb],
40
+	Matthias Maneck [ctb],
41
+	Claudio Lottaz [ctb],
42
+	Tim Beissbarth [ctb]
11 43
 Maintainer: Ramon Diaz-Uriarte <rdiaz02@gmail.com>
12 44
 Description: Functions for forward population genetic simulation in
13 45
     asexual populations, with special focus on cancer progression.
... ...
@@ -29,7 +61,7 @@ License: GPL (>= 3)
29 61
 URL: https://github.com/rdiaz02/OncoSimul, https://popmodels.cancercontrol.cancer.gov/gsr/packages/oncosimulr/
30 62
 BugReports: https://github.com/rdiaz02/OncoSimul/issues
31 63
 Depends: R (>= 3.3.0)
32
-Imports: Rcpp (>= 0.12.4), parallel, data.table, graph, Rgraphviz, gtools, igraph, methods, RColorBrewer, grDevices, car, dplyr, smatr, ggplot2, ggrepel, nem
64
+Imports: Rcpp (>= 0.12.4), parallel, data.table, graph, Rgraphviz, gtools, igraph, methods, RColorBrewer, grDevices, car, dplyr, smatr, ggplot2, ggrepel
33 65
 Suggests: BiocStyle, knitr, Oncotree, testthat (>= 1.0.0), rmarkdown, bookdown, pander
34 66
 LinkingTo: Rcpp
35 67
 VignetteBuilder: knitr
... ...
@@ -64,7 +64,7 @@ importFrom("dplyr", "full_join", "left_join", "right_join", "%>%", "mutate",
64 64
            "filter")
65 65
 importFrom("smatr", "ma") ## for major axis regression in some tests
66 66
 importFrom("car", "linearHypothesis")
67
-importFrom("nem", "transitive.reduction")
67
+## importFrom("nem", "transitive.reduction") ## now in file nem_transitive_reduction.R
68 68
 ## importFrom("slam", "simple_triplet_zero_matrix", ## "colapply_simple_triplet_matrix",
69 69
 ##            "col_sums")
70 70
 
... ...
@@ -83,7 +83,19 @@ simOGraph <- function(n, h = ifelse(n >= 4, 4, n),
83 83
     ## Prune to remove indirect connections
84 84
     if(multilevelParent & removeDirectIndirect) {
85 85
         ## adjMat <- transitiveReduction(adjMat)
86
-        trm <- nem::transitive.reduction(nem::transitive.closure(adjMat))
86
+        ## calling transitive.closure not necessary, as that
87
+        ## is called inside transitive.reduction
88
+        ## trm <- nem::transitive.reduction(nem::transitive.closure(adjMat))
89
+
90
+        ## could use relations package as
91
+        ## r1 <- relations::transitive_reduction(
92
+        ##       relations::transitive_closure(relations::as.relation(adjMat2)))
93
+        ## trm <- relation_incidence(r1)
94
+        ## but storage mode is double, etc, etc.
95
+        ## And would need to double check it is working as I think it is.
96
+        ## For now, trying to use nem's code directly
97
+        
98
+        trm <- nem_transitive.reduction(adjMat)
87 99
         stopifnot(all(trm %in% c(0L, 1L) ))
88 100
         storage.mode(trm) <- "integer"
89 101
         adjMat <- trm
... ...
@@ -100,6 +112,7 @@ simOGraph <- function(n, h = ifelse(n >= 4, 4, n),
100 112
     }
101 113
 }
102 114
 
115
+
103 116
 ## simOG <- simOGraph
104 117
 
105 118
 mysplitb <- function(n, gr) {
106 119
new file mode 100644
... ...
@@ -0,0 +1,139 @@
1
+###############################################
2
+## <code from nem package>
3
+
4
+## Title: (Dynamic) Nested Effects Models and Deterministic Effects
5
+##         Propagation Networks to reconstruct phenotypic hierarchies
6
+## Version: 2.60.0
7
+## Author: Holger Froehlich, Florian Markowetz, Achim Tresch, Theresa
8
+## Niederberger, Christian Bender, Matthias Maneck, Claudio Lottaz, Tim
9
+## Beissbarth
10
+## Maintainer: Holger Froehlich <frohlich at bit.uni-bonn.de>
11
+## https://bioconductor.org/packages/release/bioc/html/nem.html
12
+## License: GPL (>= 2)
13
+
14
+nem_transitive.reduction <- function(g){
15
+	if (!(class(g)%in%c("matrix","graphNEL"))) stop("Input must be an adjacency matrix or graphNEL object")
16
+	if(class(g) == "graphNEL"){
17
+		g = as(g, "matrix")		
18
+	}
19
+# 	if("Rglpk" %in% loadedNamespaces()){ # constraints müssen einzeln hinzugefügt und jedesmal das ILP gelöst werden. Danach müssen jedesmal die constraints überprüft werden und nicht mehr gebrauchte rausgeschmissen werden
20
+# 		g = abs(g) - diag(diag(g))
21
+# 		mat = matrix(0, ncol=sum(g), nrow=0)
22
+# 		idx = cbind(which(g == 1, arr.ind=T), 1:sum(g))		
23
+# 		for(y in 1:nrow(g)){
24
+# 			for(x in 1:nrow(g)){
25
+# 				if(g[x,y] != 0){
26
+# 					for(j in 1:nrow(g)){
27
+# 						if((g[y,j] != 0) && (g[x,j] != 0)){
28
+# 							mat.tmp = double(sum(g))
29
+# 							mat.tmp[idx[idx[,1] == x & idx[,2] == y, 3]] = 1
30
+# 							mat.tmp[idx[idx[,1] == y & idx[,2] == j, 3]] = 1
31
+# 							mat.tmp[idx[idx[,1] == x & idx[,2] == j, 3]] = -1
32
+# 							mat = rbind(mat, mat.tmp)
33
+# 						}						
34
+# 					}
35
+# 				}
36
+# 			}
37
+# 		}
38
+# 		
39
+# 		solve.problem = function(mat.tmp){
40
+# 			obj = rep(1, NCOL(mat.tmp))
41
+# 			rhs = 2
42
+# 			dir = ">="
43
+# 			sol = Rglpk_solve_LP(obj, mat.tmp, dir, rhs, max=TRUE, types=rep("B", NCOL(mat.tmp)))
44
+# 			print(sol)			
45
+# 			del = idx[which(sol$solution == 0), c(1, 2),drop=F]
46
+# 			for(i in 1:NROW(del)){
47
+# 				g[del[i,1], del[i,2]] = 0
48
+# 			}
49
+# 			g
50
+# 		}
51
+# 		
52
+# 		while(NROW(mat) > 0){
53
+# 			g = solve.problem(mat[1,,drop=F])
54
+# 			i = 1
55
+# 			while(i <= NROW(mat)){
56
+# 				idx.pos = which(mat[i,,drop=F] == 1)
57
+# 				idx.neg = which(mat[i,,drop=F] == -1)
58
+# 				xy = idx[idx[,3] %in% idx.pos, c(1,2)]
59
+# 				z = idx[idx[,3] == idx.neg, c(1,2)]
60
+# 				if(!(g[xy[1,1], xy[1,2]] == 1 & g[xy[2,1], xy[2,2]] == 1 & g[z[1], z[2]] == 1)) # remove resolved constraints
61
+# 					mat = mat[-i,,drop=F]
62
+# 				else
63
+# 					i = i + 1
64
+# 			}
65
+# 		}
66
+# 	}
67
+# 	else{
68
+# 	if(class(g) == "matrix"){		
69
+		# modified algorithm from Sedgewick book: just remove transitive edges instead of inserting them		
70
+    g = nem_transitive.closure(g, mat=TRUE) # bug fix: algorithm only works for transitively closed graphs!
71
+		g = g - diag(diag(g))
72
+		type = (g > 1)*1 - (g < 0)*1	
73
+		for(y in 1:nrow(g)){
74
+			for(x in 1:nrow(g)){
75
+				if(g[x,y] != 0){
76
+					for(j in 1:nrow(g)){
77
+						if((g[y,j] != 0) & sign(type[x,j])*sign(type[x,y])*sign(type[y,j]) != -1){ 
78
+						    g[x,j] = 0
79
+						}
80
+					}
81
+				}
82
+			}
83
+		}
84
+# 	}
85
+# 	else{		
86
+# 		nodenames=nodes(g)		
87
+# 		for(y in 1:length(nodes(g))){
88
+# 			edges = edgeL(g)
89
+# 			x = which(sapply(edges,function(l) y %in% unlist(l)))
90
+# 			j = unlist(edges[[y]])
91
+# 			cands = sapply(edges[x], function(e) list(intersect(unlist(e),j)))			
92
+# 			cands = cands[sapply(cands,length) > 0]
93
+# 			if(length(cands) > 0)
94
+# 				for(c in 1:length(cands)){ 				
95
+# 					jj = unlist(cands[c])					
96
+# 					g = removeEdge(rep(names(cands)[c],length(jj)),nodenames[jj],g)
97
+# 				}
98
+# 		}
99
+# 	}	
100
+	g		
101
+}
102
+
103
+
104
+
105
+nem_transitive.closure <- function(g,mat=FALSE,loops=TRUE){
106
+
107
+    if (!(class(g)%in%c("graphNEL","matrix"))) stop("Input must be either graphNEL object or adjacency matrix")
108
+    g <- as(g, "matrix")
109
+    
110
+    #-- adjacency matrix
111
+#     if (class(g)=="matrix"){
112
+		n <- ncol(g)
113
+		matExpIterativ <- function(x,pow,y=x,z=x,i=1) {
114
+		while(i < pow) {
115
+			z <- z %*% x
116
+			y <- y+z
117
+			i <- i+1
118
+		}
119
+		return(y)
120
+		}
121
+	
122
+		h <- matExpIterativ(g,n)
123
+		h <- (h>0)*1   
124
+		dimnames(h) <- dimnames(g)
125
+		if (!loops) diag(h) <- rep(0,n) else diag(h) <- rep(1,n)
126
+		if (!mat) h <- as(h,"graphNEL")	
127
+#     }
128
+
129
+# #     -- graphNEL object
130
+#     if (class(g)=="graphNEL"){
131
+#         tc <- RBGL::transitive.closure(g)    
132
+#         if (loops) tc$edges <- unique(cbind(tc$edges,rbind(tc$nodes,tc$nodes)),MARGIN=2)
133
+# 
134
+#         h <- ftM2graphNEL(ft=t(tc$edges),V=tc$nodes)
135
+#         if (mat) h <- as(h, "matrix")
136
+#     } 
137
+       
138
+    return(h)
139
+}
... ...
@@ -75,7 +75,22 @@ rfitness <- function(g, c= 0.5,
75 75
                      truncate_at_0 = TRUE,
76 76
                      K = 1,
77 77
                      r = TRUE,
78
-                     model = c("RMF", "NK")) {
78
+                     i = 0, # Ising, cost incompatibility
79
+                     I = -1, # Ising, sd for "i"
80
+                     circular = FALSE, # Ising, circular arrangement
81
+                     e = 0, # Eggbox, +/- e
82
+                     E = -1, # Eggbox, noise on "e"
83
+                     H = -1, # HoC stdev
84
+                     s = 0.1, # mean multiplivative
85
+                     S = -1, # SD multiplicative
86
+                     d = 0, # disminishing/increasing for multiplicative
87
+                     o = 0, # mean optimum
88
+                     O = -1, # sd optimum
89
+                     p = 0, # mean production for non 0 allele (optimum)
90
+                     P = -1, # sd for p
91
+                    
92
+                     model = c("RMF", "Additive", 
93
+                               "NK", "Ising", "Eggbox", "Full")) {
79 94
     ## Like Franke et al., 2011 and others of Krug. Very similar to Greene
80 95
     ## and Crona, 2014. And this allows moving from HoC to purely additive
81 96
     ## changing c and sd.
... ...
@@ -107,12 +122,46 @@ rfitness <- function(g, c= 0.5,
107 122
             f_det <- -c * d_reference
108 123
             ## f_det <- rowSums(m) * slope/nrow(m) ## this is Greene and Krona
109 124
             fi <- f_r + f_det
125
+        } else if (model == "Additive") {
126
+      ## get fitness effect for mutations in each gene
127
+      mutants <-rep(1,g)
128
+            ## FIXME: Why not just?
129
+            ## f_single_mut <- rnorm(g, mean = mu, sd = sd)
130
+            f_single_mut <- sapply(mutants, FUN = function(x) 
131
+                                            rnorm(x, mean = mu, sd = sd))
132
+      ## find which gene is mutated 
133
+      m2 <- m == 1
134
+      ## Sum the fitness effect of that mutation to generate a vector fi with
135
+      ## the fitness for each mutant condition
136
+      fi <- apply(m2, MARGIN = 1, FUN = function (x) sum(x*f_single_mut))
137
+      ## remove unnecessary variables
138
+      rm (f_single_mut, m2)
110 139
         } else if(model == "NK") {
111 140
             if(K >= g) stop("It makes no sense to have K >= g")
112 141
             argsnk <- paste0("-K ", K,
113 142
                              ifelse(r, " -r ", " "),
114 143
                              g, " 2")
115 144
             fl1 <- system2(fl_generate_binary(), args = argsnk, stdout = TRUE)[-1]
145
+        } else if (model == "Ising") {
146
+      argsIsing <- paste0("-i ", i, " -I ", I ,
147
+                          ifelse(circular, " -c ", " "),
148
+                          g, " 2")
149
+      fl1 <- system2(fl_generate_binary(), args = argsIsing, stdout = TRUE)[-1]
150
+    } else if (model == "Eggbox") {
151
+      argsEgg <- paste0("-e ", e, " -E ", E," ", g, " 2")
152
+      fl1 <- system2(fl_generate_binary(), args = argsEgg, stdout = TRUE)[-1]
153
+    } else if (model == "Full") {
154
+      if(K >= g) stop("It makes no sense to have K >= g")
155
+      argsFull <- paste0("-K ", K, ifelse(r, " -r ", " "),
156
+                         "-i ", i, " -I ", I , ifelse(circular, " -c ", " "),
157
+                         "-e ", e, " -E ", E, " ",
158
+                         "-H ", H, " ",
159
+                         "-s ", s, " -S ", S, " -d ", d, " ",
160
+                         "-o ", o, " -O ", O, " -p ", p, " -P ", P, " ",
161
+                         g, " 2")
162
+      fl1 <- system2(fl_generate_binary(), args = argsFull, stdout = TRUE)[-1]
163
+    }
164
+    if (model == "Eggbox" || model == "Ising" || model == "Full" || model == "NK") {
116 165
             fl1 <- matrix(
117 166
                 as.numeric(unlist(strsplit(paste(fl1, collapse = " "), " "))),
118 167
                 ncol = g + 1, byrow = TRUE)
... ...
@@ -1,3 +1,11 @@
1
+Changes in version 2.17.7 (2020-01-30):
2
+	- Removed dependency on nem, to be deprecated. Adding nem's
3
+	code (file nem_transitive_reduction.R).
4
+
5
+Changes in version 2.17.6 (2020-01-22):
6
+	- Further additions of MAGELLAN's functionality and Additive
7
+	model to rfitness
8
+
1 9
 Changes in version 2.17.5 (2019-12-19):
2 10
 	- Oooops: forgot Makevars.win
3 11
 	
... ...
@@ -1,11 +1,12 @@
1 1
 \name{rfitness}
2 2
 \alias{rfitness}
3
-
3
+\encoding{UTF-8}
4 4
 
5 5
 \title{Generate random fitness.}
6 6
 
7 7
 \description{ Generate random fitness landscapes under a House of Cards,
8
-  Rough Mount Fuji, additive model, and Kauffman's NK model.  }
8
+  Rough Mount Fuji (RMF), additive (multiplicative) model, Kauffman's NK
9
+  model, Ising model, Eggbox model and Full model}
9 10
 
10 11
 
11 12
 \usage{
... ...
@@ -14,7 +15,9 @@ rfitness(g, c = 0.5, sd = 1, mu = 1, reference = "random", scale = NULL,
14 15
          wt_is_1 = c("subtract", "divide", "force", "no"),
15 16
          log = FALSE, min_accessible_genotypes = NULL,
16 17
          accessible_th = 0, truncate_at_0 = TRUE,
17
-         K = 1, r = TRUE, model = c("RMF", "NK"))
18
+         K = 1, r = TRUE, i = 0, I = -1, circular = FALSE, e = 0, E = -1,
19
+         H = -1, s = 0.1, S = -1, d = 0, o = 0, O = -1, p = 0, P = -1, 
20
+         model = c("RMF", "Additive", "NK", "Ising", "Eggbox", "Full"))
18 21
 }
19 22
 
20 23
 
... ...
@@ -25,28 +28,32 @@ rfitness(g, c = 0.5, sd = 1, mu = 1, reference = "random", scale = NULL,
25 28
   \item{g}{Number of genes.}
26 29
 
27 30
   \item{c}{The decrease in fitness of a genotype per each unit increase
28
-    in Hamming distance from the reference genotype (see \code{reference}).}
31
+    in Hamming distance from the reference genotype for the RMF model
32
+    (see \code{reference}).}
29 33
 
30 34
   \item{sd}{The standard deviation of the random component (a normal
31
-  distribution of mean \code{mu} and standard deviation \code{sd}).}
35
+  distribution of mean \code{mu} and standard deviation \code{sd}) for
36
+  the RMF and additive models .}
32 37
 
33 38
 \item{mu}{The mean of the random component (a normal distribution of
34
-mean \code{mu} and standard deviation \code{sd}).}
35
-
36
-
37
-\item{reference}{The reference genotype: for the deterministic, additive
38
-  part, this is the genotype with maximal fitness, and all other
39
-  genotypes decrease their fitness by \code{c} for every unit of Hamming
40
-  distance from this reference. If "random" a genotype will be randomly
41
-  chosen as the reference. If "max" the genotype with all positions
42
-  mutated will be chosen as the reference. If you pass a vector (e.g.,
43
-  \code{reference = c(1, 0, 1, 0)}) that will be the reference genotype.
44
-  If "random2" a genotype will be randomly chosen as the reference. In
45
-  contrast to "random", however, not all genotypes have the same
46
-  probability of being chosen; here, what is equal is the probability
47
-  that the reference genotype has 1, 2, ..., g, mutations (and, once a
48
-  number mutations is chosen, all genotypes with that number of
49
-  mutations have equal probability of being the reference). }
39
+mean \code{mu} and standard deviation \code{sd}) for the RMF and
40
+additive models.}
41
+
42
+
43
+\item{reference}{The reference genotype: in the RMF model, for the
44
+  deterministic, additive part, this is the genotype with maximal
45
+  fitness, and all other genotypes decrease their fitness by \code{c}
46
+  for every unit of Hamming distance from this reference. If "random" a
47
+  genotype will be randomly chosen as the reference. If "max" the
48
+  genotype with all positions mutated will be chosen as the
49
+  reference. If you pass a vector (e.g., \code{reference = c(1, 0, 1,
50
+  0)}) that will be the reference genotype.  If "random2" a genotype
51
+  will be randomly chosen as the reference. In contrast to "random",
52
+  however, not all genotypes have the same probability of being chosen;
53
+  here, what is equal is the probability that the reference genotype has
54
+  1, 2, ..., g, mutations (and, once a number mutations is chosen, all
55
+  genotypes with that number of mutations have equal probability of
56
+  being the reference). }
50 57
 
51 58
 \item{scale}{Either NULL (nothing is done) or a two-element vector. If a
52 59
   two-element vector, fitness is re-scaled between \code{scale[1]} (the
... ...
@@ -127,9 +134,45 @@ mean \code{mu} and standard deviation \code{sd}).}
127 134
 
128 135
 \item{r}{For the NK model, whether interacting loci are chosen at random
129 136
   (\code{r = TRUE}) or are neighbors (\code{r = FALSE}).}
137
+\item{i}{For de Ising model, i is the mean cost for incompatibility with which
138
+  the genotype's fitness is penalized when in two adjacent genes, only one of 
139
+  them is mutated.}
140
+
141
+\item{I}{For the Ising model, I is the standard deviation for the cost 
142
+  incompatibility (i).}
143
+  
144
+\item{circular}{For the Ising model, whether there is a circular arrangement, 
145
+  where the last and the first genes are adjacent to each other.}
146
+
147
+\item{e}{For the Eggbox model, mean effect in fitness for the neighbor
148
+  locus +/- e.}
149
+  
150
+\item{E}{For the Eggbox model, noise added to the mean effect in fitness (e).}
151
+
152
+\item{H}{For Full models, standard deviation for the House of Cards model.}
153
+
154
+\item{s}{For Full models, mean of the fitness for the Multiplicative model.}
155
+
156
+\item{S}{For Full models, standard deviation for the Multiplicative model.}
157
+
158
+\item{d}{For Full models, a disminishing (negative) or increasing 
159
+  (positive) return as the peak is approached for multiplicative model.}
160
+  
161
+\item{o}{For Full models, mean value for the optimum model.}
162
+
163
+\item{O}{For Full models, standard deviation for the optimum model.}
130 164
 
131
-\item{model}{One of "RMF" (default), for Rough Mount Fuji, or "NK", for
132
-  Kauffman's NK model.}
165
+\item{p}{For Full models, the mean production value for each non 0
166
+  allele in the Optimum model component.}
167
+
168
+\item{P}{For Full models, the associated stdev (of non 0 alleles) in the
169
+Optimum model component.}
170
+
171
+
172
+
173
+\item{model}{One of "RMF" (default) for Rough Mount Fuji, "Additive" for
174
+ Additive model, "NK", for Kauffman's NK model, "Ising" for Ising model,
175
+ "Eggbox" for Eggbox model or "Full" for Full models.}
133 176
 } 
134 177
 
135 178
 
... ...
@@ -146,14 +189,56 @@ mean \code{mu} and standard deviation \code{sd}).}
146 189
   random variable (in this case, a normal deviate of mean \code{mu}
147 190
   and standard deviation \code{sd}).
148 191
 
149
-  Setting \eqn{c = 0} we obtain a House of Cards model. Setting \eqn{sd
150
-    = 0} fitness is given by the distance from the reference and if the
151
-    reference is the genotype with all positions mutated, then we have a
152
-    fully additive model (fitness increases linearly with the number of
153
-    positions mutated).
192
+  When using \code{model = "RMF"}, setting \eqn{c = 0} we obtain a House
193
+    of Cards model. Setting \eqn{sd = 0} fitness is given by the
194
+    distance from the reference and if the reference is the genotype
195
+    with all positions mutated, then we have a fully additive model
196
+    (fitness increases linearly with the number of positions mutated),
197
+    where all mutations have the same effect.
198
+
199
+  More flexible additive models can be used using \code{model =
200
+  "Additive"}. This model is like the Rough Mount Fuji model in Szendro
201
+  et al., 2013 or Franke et al., 2011, but in this case, each locus can
202
+  have different contributions to the fitness evaluation. This model is
203
+  also referred to as the "multiplicative" model in the literature as it
204
+  is additive in the log-scale (e.g., see Brouillet et al., 2015 or
205
+  Ferretti et al., 2016). The contribution of each mutated allele to the
206
+  log-fitness is a random deviate from a Normal distribution with
207
+  specified mean \code{mu} and standard deviation \code{sd}, and the
208
+  log-fitness of a genotype is the sum of the contributions of each
209
+  mutated allele. There is no "reference" genotype in the Additive
210
+  model.  There is no epistasis in the additve model because the effect
211
+  of a mutation in a locus does not depend on the genetic background, or
212
+  whether the rest of the loci are mutated or not.
213
+  
214
+
215
+  When using \code{model = "NK"} fitness is drawn from a uniform (0, 1)
216
+  distribution.
217
+  
218
+  
219
+  When using \code{model = "Ising"} for each pair of interacting loci, 
220
+  there is an associated cost if both alleles are not identical 
221
+  (and therefore 'compatible').
222
+  
223
+  
224
+  When using \code{model = "Eggbox"} each locus is either high or low fitness,
225
+  with a systematic change between each neighbor.
226
+  
227
+  
228
+  When using \code{model = "Full"}, the fitness is computed with different
229
+  parts of the previous models depending on the choosen parameters described 
230
+  above. 
231
+  
232
+  
233
+  For \code{model = "NK" | "Ising" | "Eggbox" | "Full"} the fitness
234
+  landscape is generated by directly calling the \code{fl_generate}
235
+  function of MAGELLAN
236
+  (\url{http://wwwabi.snv.jussieu.fr/public/Magellan/}). See details in
237
+  Ferretti et al. 2016, or Brouillet et al., 2015.
238
+  
154 239
 
155 240
   For OncoSimulR, we often want the wildtype to have a mean of
156
-  1. Reasonable settings are \code{mu = 1} and \code{wt_is_1 =
241
+  1. Reasonable settings when using RMF are \code{mu = 1} and \code{wt_is_1 =
157 242
   'subtract'} so that we simulate from a distribution centered in 1, and
158 243
   we make sure afterwards (via a simple shift) that the wildtype is
159 244
   actuall 1. The \code{sd} controls the standard deviation, with the
... ...
@@ -162,14 +247,6 @@ mean \code{mu} and standard deviation \code{sd}).}
162 247
   of the data can be large, specially if \code{g} (the number of genes)
163 248
   is large.
164 249
 
165
-
166
-  When using \code{model = "NK"}, the model used is Kauffman's NK model
167
-  (see details in Ferretti et al., or Brouillet et al., below), as
168
-  implemented in MAGELLAN
169
-  (\url{http://wwwabi.snv.jussieu.fr/public/Magellan/}). This fitness
170
-  landscape is generated by directly calling the \code{fl_generate}
171
-  function of MAGELLAN. Fitness is drawn from a uniform (0, 1)
172
-  distribution.
173 250
   
174 251
 } 
175 252
 
... ...
@@ -214,10 +291,12 @@ MAGELLAN web site: \url{http://wwwabi.snv.jussieu.fr/public/Magellan/}
214 291
 }
215 292
 
216 293
 \author{ Ramon Diaz-Uriarte for the RMF and general wrapping
217
-code. S. Brouillet, G. Achaz, S. Matuszewski, H. Annoni, and L. Ferreti
218
-for the MAGELLAN code.
219
-
220
-}
294
+  code. S. Brouillet, G. Achaz, S. Matuszewski, H. Annoni, and
295
+  L. Ferreti for the MAGELLAN code. Further contributions to the
296
+  additive model and to wrapping MAGELLAN code and documentation from
297
+  Guillermo Gorines Cordero, Ivan Lorca Alonso, Francisco Muñoz Lopez,
298
+  David Roncero Moroño, Alvaro Quevedo, Pablo Perez, Cristina Devesa,
299
+  Alejandro Herrador.}
221 300
 
222 301
 \seealso{
223 302
   
... ...
@@ -234,6 +313,7 @@ for the MAGELLAN code.
234 313
 ## Random fitness for four genes-genotypes,
235 314
 ## plotting and simulating an oncogenetic trajectory
236 315
 
316
+
237 317
 r1 <- rfitness(4)
238 318
 plot(r1)
239 319
 oncoSimulIndiv(allFitnessEffects(genotFitness = r1))
... ...
@@ -243,7 +323,26 @@ oncoSimulIndiv(allFitnessEffects(genotFitness = r1))
243 323
 rnk <- rfitness(5, K = 3, model = "NK")
244 324
 plot(rnk)
245 325
 oncoSimulIndiv(allFitnessEffects(genotFitness = rnk))
246
-}
247 326
 
327
+## Additive model
328
+radd <- rfitness(4, model = "Additive", mu = 0.2, sd = 0.5)
329
+plot(radd)
330
+
331
+## Eggbox model
332
+regg = rfitness(g=4,model="Eggbox", e = 2, E=2.4)
333
+plot(regg)
334
+
335
+
336
+## Ising model
337
+ris = rfitness(g=4,model="Ising", i = 0.002, I=2)
338
+plot(ris)
339
+
340
+
341
+## Full model
342
+rfull = rfitness(g=4, model="Full", i = 0.002, I=2, 
343
+                 K = 2, r = TRUE,
344
+                 p = 0.2, P = 0.3, o = 0.3, O = 1)
345
+plot(rfull)
346
+}
248 347
 \keyword{ datagen }
249 348
 
... ...
@@ -84,7 +84,10 @@ sh. See \code{\link{allFitnessEffects}}}
84 84
   input, as "restriction table" in  \code{\link{allFitnessEffects}}.
85 85
 
86 86
 }
87
-\author{Ramon Diaz-Uriarte}
87
+\author{Ramon Diaz-Uriarte. Code for transitive closure taken from the
88
+  nem package, whose authors are Holger Froehlich, Florian Markowetz,
89
+  Achim Tresch, Theresa Niederberger, Christian Bender, Matthias Maneck,
90
+  Claudio Lottaz, Tim Beissbarth}
88 91
 
89 92
 \examples{
90 93
 (a1 <- simOGraph(10))
... ...
@@ -57,6 +57,49 @@ test_that("Error if wrong arguments", {
57 57
                  fixed = TRUE)
58 58
 })
59 59
 
60
+test_that("Additive, Ising, Eggbox, Full exercising", {
61
+
62
+    expect_output(print(rfitness(4, model = "Ising")), "Fitness", fixed = TRUE)
63
+expect_output(print(rfitness(4, model = "Eggbox")), "Fitness", fixed = TRUE)
64
+expect_output(print(rfitness(4, model = "Full")), "Fitness", fixed = TRUE)
65
+expect_output(print(rfitness(4, model = "Additive")), "Fitness", fixed = TRUE)
66
+expect_output(print(rfitness(4, model = "Ising", i = 0.0002, I = 0.5,
67
+                             circular = TRUE)), "Fitness", fixed = TRUE)
68
+expect_output(print(rfitness(4, model = "Ising", i = 2, I = 0,
69
+                             circular = TRUE)), "Fitness", fixed = TRUE)
70
+expect_output(print(rfitness(4, model = "Ising", i = 2, I = -2,
71
+                             circular = TRUE)), "Fitness", fixed = TRUE)
72
+expect_output(print(rfitness(4, model = "Ising", i = 2, I = 0.5,
73
+                             circular = FALSE)), "Fitness", fixed = TRUE)
74
+expect_output(print(rfitness(4, model = "Eggbox", e = 2, E = 0)),
75
+              "Fitness", fixed = TRUE)
76
+expect_output(print(rfitness(4, model = "Eggbox", e = 2, E = 2.4)),
77
+              "Fitness", fixed = TRUE)
78
+expect_output(print(rfitness(4, model = "Eggbox", e = 0, E = 2.4)),
79
+              "Fitness", fixed = TRUE)
80
+expect_output(print(rfitness(4, model = "Full", i = 0.002, I = 2,
81
+                K = 2, r = TRUE,
82
+                p = 0.2, P = 0.3, o = 0.3, O = 1)), "Fitness", fixed = TRUE)
83
+expect_error(rfitness(4, model = "Full", K = 5), "It makes no sense",
84
+              fixed = TRUE)
85
+expect_output(print(rfitness(4, model = "Full", i = 0.002, I = 2,
86
+                K = 2, r = TRUE,
87
+                p = 0.2, P = 0.3, o = 0.3, O = 1,
88
+                s = 0.5, S=0, d = 0.1,
89
+                e = 1, E = 1.5,
90
+                H = 0.5)), "Fitness", fixed = TRUE)
91
+expect_output(print(rfitness(4, model = "Additive", mu = 0, sd = 1)),
92
+              "Fitness", fixed = TRUE)
93
+expect_output(print(rfitness(4, model = "Additive", mu = 1, sd = 0)),
94
+              "Fitness", fixed = TRUE)
95
+expect_output(print(rfitness(4, model = "Additive", mu = 3, sd = 1)),
96
+              "Fitness", fixed = TRUE)
97
+expect_output(print(rfitness(4, model = "Additive", mu = -4, sd = 10)),
98
+              "Fitness", fixed = TRUE)
99
+
100
+})
101
+
102
+
60 103
 cat(paste("\n Ending exercise-rfitness at", date(), "\n"))
61 104
 cat(paste("  Took ", round(difftime(Sys.time(), inittime, units = "secs"), 2), "\n\n"))
62 105
 rm(inittime)
... ...
@@ -8455,6 +8455,100 @@ plot(rnk)
8455 8455
 oncoSimulIndiv(allFitnessEffects(genotFitness = rnk))
8456 8456
 ```
8457 8457
 
8458
+
8459
+## Random fitness landscapes from an additive model {#additivemodel}
8460
+
8461
+This model evaluates fitness with different contributions of each allele,
8462
+which will be randomly generated.
8463
+
8464
+Given a number a genes by the user, the code uses rnorm to generate random
8465
+contribution for the mutated allele in each locus. Later, this
8466
+constributions will be used in the generation of the matrix that gives the
8467
+value of fitness for each combination of wild type/mutated alleles by
8468
+addition of the values for each locus and combination.
8469
+
8470
+```{r addex1}
8471
+radd <- rfitness(4, model = "Additive", mu = 0, sd = 0.5)
8472
+plot(radd)
8473
+```
8474
+
8475
+## Random fitness landscapes from Eggbox model {#eggboxmodel}
8476
+
8477
+You can also use Eggbox model [e.g., @ferretti_measuring_2016;
8478
+@brouillet_magellan:_2015], where each locus is either high or low fitness
8479
+(depending on the "e" parameter value), with a systematic change between
8480
+each neighbor. We call the function `fl_generate` from MAGELLAN
8481
+[@brouillet_magellan:_2015] to generate these landscapes.
8482
+
8483
+
8484
+<!-- "e" parameter: every other genotype is +/- e. -->
8485
+<!-- "E" parameter: add noise on the mean effect for eggbox. -->
8486
+
8487
+
8488
+
8489
+```{r eggex1}
8490
+regg <- rfitness(4, model = "Eggbox", e = 1, E = 0.5)
8491
+plot(regg)
8492
+```
8493
+
8494
+## Random fitness landscapes from Ising model {#isingmodel}
8495
+
8496
+In the Ising model [e.g., @ferretti_measuring_2016;
8497
+@brouillet_magellan:_2015], loci are arranged sequentially and each
8498
+locus interacts with its physical neighbors. For each pair of interacting
8499
+loci, there is a cost to (log)fitness if both alleles are not identical (and
8500
+therefore 'compatible'); in this case, the cost for incompatibility `i` is
8501
+applied. The last and the first loci will interact only if 'circular' is
8502
+set. The implementation of this model is decribedin
8503
+[@brouillet_magellan:_2015], and we use a call to MAGELLAN code to
8504
+generate the landscape.
8505
+
8506
+
8507
+
8508
+
8509
+```{r isingex1}
8510
+ris <- rfitness(g = 4, model = "Ising", i = 0.002, I = 0.45)
8511
+plot(ris)
8512
+```
8513
+
8514
+
8515
+## Random fitness landscapes from Full models {#fullmodel}
8516
+
8517
+MAGELLAN also offers the possibility to combine different models with their own
8518
+parameters in order to generate a Full model. The models combined are:
8519
+  
8520
+* House of cards: `H` is the number of interacting genes. 
8521
+* Multiplicative (what we call additive): `s` and `S` mean and SD for generating random fitnesses. `d` 
8522
+is a diminishing (negative) or increasing (positive) return as you approach the peak.
8523
+* Kauffman NK: each locus interacts with any other `K` loci that can be chosen
8524
+randomly pasing "r = TRUE" or among its neigbors.
8525
+* RMF: explained at \@ref(rmfmodel).
8526
+* Ising: `i` and `I` mean and SD for incompatibility. If `circular` option is
8527
+provided, the last and first alleles can interact (circular arrangement).
8528
+* Eggbox: In this model, each locus is considered as high or low fitness. From
8529
+one locus to another the fitness changes its sign so that one is benefficial
8530
+and its neigbor is detrimental.
8531
+`e` and `E` are fitness and noise for fitness.
8532
+* Optimum: there is an optimum fitness contribution defined by `o` $mu$ and 
8533
+`O` $sigma$ and every locus has a production `p`and `P` (also mean and sd 
8534
+                                                         respectively).
8535
+
8536
+All models can be taken into account for the fitness calculation.  With
8537
+default parameters, neither of Ising, Eggbox or Optimum contribute to
8538
+fitness lanadcape generation as all `i`, `e`, `o` and `p` all `== 0`.
8539
+Also, as all parameters refering to standard deviations have value `==
8540
+-1`, those are also have no effect unless changed. Further details can be
8541
+found in MAGELLAN's webpage
8542
+<http://wwwabi.snv.jussieu.fr/public/Magellan/> and
8543
+[@brouillet_magellan:_2015].
8544
+
8545
+```{r fullex1}
8546
+rnk <- rfitness(4, model = "Full", i = 0.5, I = 0.1,
8547
+                e = 2, s = 0.3, S = 0.08)
8548
+plot(rnk)
8549
+```
8550
+
8551
+
8458 8552
 ## Epistasis and fitness landscape statistics {#epistats}
8459 8553
 
8460 8554
 We can call MAGELLAN's [@brouillet_magellan:_2015] `fl_statistics` to
... ...
@@ -1,15 +1,15 @@
1 1
 \usepackage[%
2
-		shash={912503a},
3
-		lhash={912503a07aa95ca7ba70aa1022276cb71f1cab46},
2
+		shash={5f42713},
3
+		lhash={5f42713ca15e08bcae15e39a2c9e5c9f239b2e99},
4 4
 		authname={ramon diaz-uriarte (at Phelsuma)},
5 5
 		authemail={rdiaz02@gmail.com},
6
-		authsdate={2019-12-19},
7
-		authidate={2019-12-19 09:18:41 +0100},
8
-		authudate={1576743521},
6
+		authsdate={2020-01-30},
7
+		authidate={2020-01-30 14:08:28 +0100},
8
+		authudate={1580389708},
9 9
 		commname={ramon diaz-uriarte (at Phelsuma)},
10 10
 		commemail={rdiaz02@gmail.com},
11
-		commsdate={2019-12-19},
12
-		commidate={2019-12-19 09:18:41 +0100},
13
-		commudate={1576743521},
11
+		commsdate={2020-01-30},
12
+		commidate={2020-01-30 14:08:28 +0100},
13
+		commudate={1580389708},
14 14
 		refnames={ (HEAD -> master, origin/master, origin/HEAD)}
15 15
 	]{gitsetinfo}
16 16
\ No newline at end of file