[adegenet-commits] r1020 - in pkg: R inst

[noreply at r-forge.r-project.org]

Author: jombart
Date: 2012-07-11 18:09:53 +0200 (Wed, 11 Jul 2012)
New Revision: 1020

Removed:
   pkg/inst/doc/
Modified:
   pkg/R/haploGen.R
Log:
- removed useless inst/doc
- fixed substitution process in haplogen -> now different types of mutations.


Modified: pkg/R/haploGen.R
===================================================================
--- pkg/R/haploGen.R	2012-07-10 17:52:45 UTC (rev 1019)
+++ pkg/R/haploGen.R	2012-07-11 16:09:53 UTC (rev 1020)
@@ -8,10 +8,10 @@
 ## mean.gen.time, sd.gen.time: average time for transmission and its standard deviation (normal dist)
 ## mean.repro, sd.repro: average number of transmissions and its standard deviation (normal dist)
 ##
-haploGen <- function(seq.length=10000, mu=0.0001, t.max=20,
-                     gen.time=function(){round(rnorm(1,5,1))},
-                     repro=function(){round(rnorm(1,2,1))}, max.nb.haplo=1e3,
-                     geo.sim=FALSE, grid.size=5, lambda.xy=0.5,
+haploGen <- function(seq.length=1e4, mu.transi=1e-4, mu.transv=mu.transi/2, t.max=20,
+                     gen.time=function(){1+rpois(1,0.5)},
+                     repro=function(){rpois(1,1.5)}, max.nb.haplo=1e3,
+                     geo.sim=FALSE, grid.size=10, lambda.xy=0.5,
                      mat.connect=NULL,
                      ini.n=1, ini.xy=NULL){
 
@@ -19,55 +19,27 @@
     if(!require(ape)) stop("The ape package is required.")
 
 
-    ## ## HACK TO FIX APE'S BUG ##
-    ## env <- environment(as.list.DNAbin)
-    ## as.list.DNAbin.new <- function (x, ...){
-    ##     if (is.list(x))
-    ##         return(x)
-    ##     if (is.null(dim(x)))
-    ##         obj <- list(x)
-    ##     else {
-    ##         n <- nrow(x)
-    ##         obj <- vector("list", n)
-    ##         for (i in 1:n) obj[[i]] <- x[i, ]
-    ##         names(obj) <- rownames(x)
-    ##     }
-    ##     class(obj) <- "DNAbin"
-    ##     obj
-    ## }
-
-    ## as.list.DNAbin <- as.list.DNAbin.new
-    ## unlockBinding("as.list.DNAbin", env)
-    ## assignInNamespace("as.list.DNAbin", as.list.DNAbin.new, ns="ape", envir=env)
-    ## assign("as.list.DNAbin", as.list.DNAbin.new, envir=env)
-    ## lockBinding("as.list.DNAbin", env)
-
-    ## END HACK ##
-
-
-
     ## HANDLE ARGUMENTS ##
-    if(is.numeric(mu)){
-        mu.val <- mu
-        mu <- function(){return(mu.val)}
-    }
-
+    ## if numeric vector, taken as proba for t=0,1,2,...(length-1)
     if(is.numeric(gen.time)){
-        gen.time.val <- gen.time
-        gen.time <- function(){gen.time.val}
+        gen.time.val <- gen.time/sum(gen.time)
+        gen.time <- function(){sample(0:(length(repro.val)-1), size=1, prob=gen.time.val)}
     }
 
+    ## if numeric vector, taken as proba for ndesc=0,1,2,...(length-1)
     if(is.numeric(repro)){
-        repro.val <- repro
-        repro <- function(){repro.val}
+        repro.val <- repro/sum(repro)
+        repro <- function(){sample(0:(length(repro.val)-1),size=1,prob=repro.val)}
     }
 
 
+
     ## GENERAL VARIABLES ##
     NUCL <- as.DNAbin(c("a","t","c","g"))
+    TRANSISET <- list('a'=as.DNAbin('g'), 'g'=as.DNAbin('a'), 'c'=as.DNAbin('t'), 't'=as.DNAbin('c'))
+    TRANSVSET <- list('a'=as.DNAbin(c('c','t')), 'g'=as.DNAbin(c('c','t')), 'c'=as.DNAbin(c('a','g')), 't'=as.DNAbin(c('a','g')))
     res <- list(seq=as.matrix(as.DNAbin(character(0))), dates=integer(), ances=character())
     toExpand <- logical()
-    ##mu <- mu/365 # mutation rate by day
     myGrid <- matrix(1:grid.size^2, ncol=grid.size, nrow=grid.size)
 
 
@@ -80,21 +52,42 @@
         return(res)
     }
 
-    ## create substitutions for defined SNPs
+    ## create substitutions for defined SNPs - no longer used
     substi <- function(snp){
-        res <- sapply(1:length(snp), function(e) sample(setdiff(NUCL,e),1)) # ! sapply does not work on DNAbin vectors directly
+        res <- sapply(1:length(snp), function(i) sample(setdiff(NUCL,snp[i]),1)) # ! sapply does not work on DNAbin vectors directly
         class(res) <- "DNAbin"
         return(res)
     }
 
+    ## create transitions for defined SNPs
+    transi <- function(snp){
+        res <- unlist(TRANSISET[as.character(snp)])
+        class(res) <- "DNAbin"
+        return(res)
+    }
+
+    ## create transversions for defined SNPs
+    transv <- function(snp){
+        res <- sapply(TRANSVSET[as.character(snp)],sample,1)
+        class(res) <- "DNAbin"
+        return(res)
+    }
+
     ## duplicate a sequence (including possible mutations)
     seq.dupli <- function(seq, T){ # T is the number of time units between ancestor and decendent
-        ## toChange <- as.logical(rbinom(n=seq.length, size=1, prob=mu())) # can be faster
-        temp <- rbinom(n=1, size=seq.length*T, prob=mu()) # total number of mutations
-        toChange <- sample(1:seq.length, size=temp, replace=FALSE)
-        if(sum(toChange)>0) {
-            seq[toChange] <- substi(seq[toChange])
+        ## transitions ##
+        n.transi <- rbinom(n=1, size=seq.length*T, prob=mu.transi) # total number of transitions
+        if(n.transi>0) {
+            idx <- sample(1:seq.length, size=n.transi, replace=FALSE)
+            seq[idx] <- transi(seq[idx])
         }
+
+        ## transversions ##
+        n.transv <- rbinom(n=1, size=seq.length*T, prob=mu.transv) # total number of transitions
+        if(n.transv>0) {
+            idx <- sample(1:seq.length, size=n.transv, replace=FALSE)
+            seq[idx] <- transv(seq[idx])
+        }
         return(seq)
     }
 
@@ -107,7 +100,7 @@
     ## how many days before duplication occurs ?
     time.dupli <- function(){
         ##res <- round(rnorm(1, mean=mean.gen.time, sd=sd.gen.time))
-        res <- gen.time()
+        res <- round(gen.time()) # force integers
         res[res<0] <- 0
         return(res)
     }