[adegenet-commits] r486 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 17 18:34:41 CET 2009
Author: jombart
Date: 2009-11-17 18:34:41 +0100 (Tue, 17 Nov 2009)
New Revision: 486
Modified:
pkg/R/haploPop.R
Log:
not yet finished, have to be able to reuse an haploPop when tracking diversity (otherwise works)
Modified: pkg/R/haploPop.R
===================================================================
--- pkg/R/haploPop.R 2009-11-17 16:10:24 UTC (rev 485)
+++ pkg/R/haploPop.R 2009-11-17 17:34:41 UTC (rev 486)
@@ -468,3 +468,251 @@
return(invisible(tre))
} # end plot.haploPop
+
+
+
+
+
+
+
+
+
+
+############
+## haploPopDiv
+############
+haploPopDiv <- function(n.steps=20, ini.obj=NULL, haplo.length=1e6, mu=1e-5, n.snp.ini=1,
+ birth.func=function(){ sample(0:3, 1, prob=c(.05, .45, .35, .15))},
+ max.pop.size=function(){1e4}, max.nb.pop=30, ini.pop.size=10, regen=FALSE,
+ p.new.pop=function(){1e-4}, kill.func=function(age){age>1},
+ quiet=FALSE, clean.haplo=FALSE) {
+
+
+ ## SOME CHECKS
+ ## if(is.numeric(ini.pop.size)){
+ ## ini.pop.size.val <- ini.pop.size
+ ## ini.pop.size <- function(){ini.pop.size.val}
+ ## }
+
+ if(is.numeric(max.pop.size)){
+ max.pop.size.val <- max.pop.size
+ max.pop.size <- function(){max.pop.size.val}
+ }
+
+ if(is.numeric(p.new.pop)){
+ p.new.pop.val <- p.new.pop
+ p.new.pop <- function(){p.new.pop.val}
+ }
+
+ if(is.numeric(birth.func)){
+ birth.func.val <- birth.func[1]
+ birth.func <- function(){birth.func.val}
+ }
+
+ if(is.numeric(kill.func)){
+ kill.func.val <- kill.func[1]
+ kill.func <- function(age){age>kill.func.val}
+ }
+
+
+ ## GLOBAL VARIABLES ##
+ SNP.POOL <- 1:haplo.length
+ vecS <- 1 # will be redefined later, but needed for evolveOnePop definition
+
+ ## AUXILIARY FUNCTIONS ##
+ createMutations <- function(N){ # L:genome length; N: pop size
+ nb.mutations <- sum(rbinom(N, size=haplo.length, prob=mu))
+ return( sample(SNP.POOL, size=nb.mutations, replace=TRUE) )
+ }
+
+ assignMutations <- function(myPop, mutations){ # mypop: list of `haplotypes'; mutations: vector of SNPs
+ if(length(mutations)==0 | length(myPop)==0) return(myPop)
+ id <- sample(1:length(myPop), size=length(mutations), replace=TRUE)
+ mutations <- split(mutations, id)
+ myPop[as.integer(names(mutations))] <- mapply(c, myPop[as.integer(names(mutations))], mutations, SIMPLIFY=FALSE)
+ return(myPop)
+ }
+
+ if(!regen){
+ ## VERSION FOR NO REGENERATION OF SUSCEPTIBLES
+ evolveOnePop <- function(myPop, myS, myAge){ # myPop: pop to evolve; myS: nb of susceptible in the pop; myAge: vector of ages
+ ## kill 'em bastards (= old strains)
+ myAge <- myAge + 1
+ toKill <- kill.func(myAge)
+ myPop[toKill] <- NULL
+ myAge <- myAge[!toKill]
+
+ ## generate new strains for new generation
+ sampSize <- round(min( length(myPop)*birth.func(), myS)) # number of strains for next step
+ if(sampSize<1){ # if no sample
+ return(list(pop=myPop, S=myS, age=myAge))
+ }
+ newGen <- myPop[sample(1:length(myPop), sampSize, replace=TRUE)] # sample strains for new generations
+ newGen <- assignMutations(newGen, createMutations(sampSize)) # mutate strains
+ newAge <- rep(0, sampSize) # new ages for newborns
+
+ ## merge old and new generation
+ myPop <- c(myPop,newGen)
+ myAge <- c(myAge, newAge)
+
+ ## possibly create one or more new pop
+ if((length(listPop) < max.nb.pop) & (p.new.pop()>0)) { # total number of pop. limitation
+ nbNewPop <- rbinom(1, length(myPop), prob=p.new.pop())
+ } else {
+ nbNewPop <- 0
+ }
+ if(nbNewPop>0){
+ ## newPop <- sample(listPop, size=nbNewPop, replace=TRUE) # wrong
+ newPop <- lapply(sample(myPop, size=nbNewPop, replace=TRUE), as.list)
+ listPop <<- c(listPop, newPop)
+ vecS <<- c(vecS, replicate(nbNewPop, max.pop.size()) )
+ listAges <<- c(listAges, replicate(nbNewPop, 0, simplify=FALSE) )
+ } # end new pop
+ return(list(pop=myPop, S=myS-sampSize, age=myAge))
+ } # end no regen version
+ } else { ## REGEN VERSION
+ evolveOnePop <- function(myPop, myS, myAge){ # myPop: pop to evolve; myS: nb of susceptible in the pop; myAge: vector of ages
+ ## kill 'em bastards (= old strains)
+ myAge <- myAge + 1
+ toKill <- kill.func(myAge)
+ myPop[toKill] <- NULL
+ myAge <- myAge[!toKill]
+ myS <- max.pop.size() ## DIFFERENCE between the two versions of the function
+
+ ## generate new strains for new generation
+ sampSize <- round(min( length(myPop)*birth.func(), myS)) # number of strains for next step
+ if(sampSize<1){ # if no sample
+ return(list(pop=myPop, S=myS, age=myAge))
+ }
+ newGen <- myPop[sample(1:length(myPop), sampSize, replace=TRUE)] # sample strains for new generations
+ newGen <- assignMutations(newGen, createMutations(sampSize)) # mutate strains
+ newAge <- rep(0, sampSize) # new ages for newborns
+
+ ## merge old and new generation
+ myPop <- c(myPop,newGen)
+ myAge <- c(myAge, newAge)
+
+ ## possibly create one or more new pop
+ if((length(listPop) < max.nb.pop) & (p.new.pop()>0)) { # total number of pop. limitation
+ nbNewPop <- rbinom(1, length(myPop), prob=p.new.pop())
+ } else {
+ nbNewPop <- 0
+ }
+ if(nbNewPop>0){
+ ## newPop <- sample(listPop, size=nbNewPop, replace=TRUE) # wrong
+ newPop <- lapply(sample(myPop, size=nbNewPop, replace=TRUE), as.list)
+ listPop <<- c(listPop, newPop)
+ vecS <<- c(vecS, replicate(nbNewPop, max.pop.size()) )
+ listAges <<- c(listAges, replicate(nbNewPop, 0, simplify=FALSE) )
+ } # end new pop
+ return(list(pop=myPop, S=myS, age=myAge)) ## DIFFERENCE between the two versions of the function
+ } # end no regen version
+ } ## end evolveOnePop (both versions)
+
+
+
+ ## INITIATE SIMULATIONS ##
+ ## INITIALIZE FROM SCRATCH
+ if(is.null(ini.obj)){
+ vecS <- max.pop.size() - n.snp.ini # susceptibles
+ haplo.ini <- sample(SNP.POOL, n.snp.ini, replace=TRUE)
+ ANCES <- haplo.ini
+ listPop <- list()
+ listPop[[1]] <- lapply(1:ini.pop.size, function(i) haplo.ini) # contains only one population of identical clones to start with
+ listAges <- list() # will contain vectors of ages of haplotypes (a time of appearance, age=0)
+ listAges[[1]] <- rep(0, ini.pop.size)
+ } else { ## INITIALIZE WITH PROVIDED OBJECT
+ if(!inherits(ini.obj, "haploPop")) stop("x is not a haploPopDiv object")
+ vecS <- ini.obj$S
+ ANCES <- attr(ini.obj, "ances")
+ listPop <- ini.obj$pop
+ listAges <- ini.obj$ages
+ }
+
+ res <- list(tab=list(), popSize=integer())
+ res$tab[[1]] <- table(unlist(listPop))
+ res$popSize[1] <- sum(sapply(listPop, length))
+
+ ## MAKE SIMULATIONS ##
+
+ ## evolve all populations
+ i <- 1L
+ if(!quiet){
+ cat("\nSimulating populations of haplotypes through time: \n")
+ }
+ ##while((sum(vecS)>0) & (i<(n.steps+1))){ # evolve all generations
+ while(i<(n.steps+1)){ # evolve all generations
+ i <- i + 1L # update iterator
+ if(!quiet){
+ catStep <- max(round(n.steps/100), 10)
+ cat(ifelse((i %% catStep)==0, paste(" ...", i), ""))
+ }
+
+
+ ## make populations evolve of one generation
+ ##idx <- which(vecS>0) # make sure that new pop won't evolve this time ! leads to not dying
+ idx <- 1:length(listPop) # make sure that new pop won't evolve this time
+ if(length(idx)>0){
+ for(j in idx){
+ temp <- evolveOnePop(listPop[[j]], vecS[j], listAges[[j]])
+ listPop[[j]] <- temp$pop
+ vecS[j] <- temp$S
+ listAges[[j]] <- temp$age
+ }
+
+ }
+
+ ## ## purge non-susceptible pop
+ ## listPop <- listPop[vecS>0]
+ ## vecS <- vecS[vecS>0]
+
+ ## purge empty populations
+ toKeep <- sapply(listPop, length)>0
+ listPop <- listPop[toKeep]
+ vecS <- vecS[toKeep]
+ listAges <- listAges[toKeep]
+
+ res$tab[[i]] <- table(unlist(listPop))
+ res$popSize[i] <- sum(sapply(listPop, length))
+
+ ## stop if all pop go extinct
+ if(length(listPop)==0L){
+ cat("\n All populations went extinct at time",i,"\n")
+
+ return(res)
+ }
+
+
+ ## FOR DEBUGGING
+ ## cat("\n=== ",i," ===")
+ ## cat("\nlistPop")
+ ## print(listPop)
+ ## cat("\nvecS")
+ ## print(vecS)
+ ## cat("\nlistAges")
+ ## print(listAges)
+ ## END DEBUGGING
+ } # end while
+
+ if(!quiet){
+ cat("\n... done! \n")
+ }
+
+ ## END OF SIMULATIONS ##
+
+ ## STORE HAPLOPOP OBJECT
+ obj <- list(pop=listPop, ages=listAges, S=vecS)
+ class(obj) <- "haploPop"
+ obj$call <- match.call()
+ attr(obj,"ances") <- ANCES # ancestral genotype
+
+ cat("\nStored haploPop object in 'last.haploPop'\n")
+ assign("last.haploPop", obj, envir= .GlobalEnv)
+
+
+ ## RETURN RES
+ return(res)
+
+} # end haploPopDiv
+
+
More information about the adegenet-commits
mailing list