# Sample parameters and generates ms string
# @description this function sample parameters from priors and trasform them to coalescent scale to generate a ms command from a model object generated by the Model Builder, main.menu() function.
# @param model A model object
# @param use.alpha Logical.If TRUE the most recent population size change will be exponential. If FALSE sudden demographic changes.
# @return a list with msABC command and sampled parameters.
# @note This function is used internally for the sim.snp.sumstat function. One my want to run this function to check the ms string.
#
ms.commander2<-function(model,use.alpha=use.alpha){
# empty parameter vector
parameters<-vector()
# bind Ne, mig and Time priors
size.pars<-rbind(model$flags$n,model$flags$en$size)
mig.pars<-rbind(model$flags$m,model$flags$em$size)
time.pars<-rbind(model$flags$ej,model$flags$en$time,model$flags$em$time)
# sample Ne, div.time and mutation rate
size.pars <- PipeMaster:::sample.w.cond(par.matrix=size.pars,cond.matrix=model$conds$size.matrix)
# bind Ne sampled parameters
parameters <- rbind(parameters, size.pars[,c(1,4)])
if(is.null(time.pars)==F){
time.pars <- PipeMaster:::sample.w.cond(par.matrix=time.pars,cond.matrix=model$conds$time.matrix)
# bind sampled time parameters
parameters <- rbind(parameters,time.pars[,c(1,4)])
}
loci <- PipeMaster:::sample.pars(model$loci)
# sample migrations if present and bind sampled parameters
if(is.null(mig.pars)==F){
mig.pars <- PipeMaster:::sample.w.cond(par.matrix=mig.pars,cond.matrix=model$conds$mig.matrix)
#bind sampled migration parameters
parameters<-rbind(parameters,mig.pars[,c(1,4)])
}
#### bind sampled mutation rate
parameters <- rbind(parameters,loci[,c(1,4)])
####### End of parameter sampling #######################################
#########################################################################
####### Generate ms string ##############################################
####### Conversion to coalescent scale #####################################
# generate coalescent scalar. Arbitrary value
#### if single population
if(model$I[1,3]=="1"){
Ne0<-as.numeric(size.pars[1,4])
ms.scalar<-4*Ne0
} else {
Ne0<-mean(as.numeric(model$flags$n[,4:5]))
ms.scalar<-4*Ne0
}
#### bind scaled theta per gene (4Ne0*m*pb)
loci <- cbind(loci, ms.scalar * as.numeric(loci[,4]) * as.numeric(loci[,2]))
#### convertion of time to coalescent scale
time.pars[,4:5] <- as.numeric(time.pars[,4])/ms.scalar
# rescale to inheritance scalar and transform size parameters to relative to Ne0
size.pars[,4:5] <- as.numeric(size.pars[,4])/Ne0
commands<-list(NULL)
for(u in 1:nrow(loci)){
string <- PipeMaster:::ms.string.generator(model,size.pars,mig.pars,time.pars,use.alpha=use.alpha, scalar = as.numeric(loci[u,3]))
#################################### theta and structure ###########################
######### generate -t and -I part of the command
if(model$I[1,3]!="1"){
y<-paste("-t",loci[u,7],paste(model$I[u,2:ncol(model$I)],collapse=" "),collapse=" ")
} else {
y<-paste("-t",loci[u,7],collapse=" ")
}
commands[[u]]<-paste(y,string, collapse=" ")
}
#### attach sampled parameters
commands[[nrow(loci)+1]]<-t(rbind(parameters,c("scalar",ms.scalar)))
return(commands)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.