R/msABC.commander.R

In gehara/PipeMaster: Building and Simulating Coalescent models.

# Sample parameters and generates msABC string
# @description this function sample parameters from priors and trasform them to coalescent scale to generate a msABC 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.msABC function. One mey want to run this function to check the msABC string.
#
msABC.commander<-function(model,use.alpha=use.alpha,arg){

  # 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<-model$loci

  # sample migrations if present and bind sampled parameters
  if(is.null(mig.pars)==F){
    mig.pars<-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<-parameters#,loci[,c(1,4)])

  ####### End of parameter sampling #######################################
  #########################################################################

  ####### Generate ms string ##############################################
  ####### Convertion 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
  # }
  Ne0 <- 100000
  ms.scalar <- 4*Ne0

  # rescale to inheritance scalar and transform size parameters to relative to Ne0
  size.pars[,4:5] <- as.numeric(size.pars[,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

  commands <- list(NULL,NULL)
  #### ms string command
  string <- PipeMaster:::ms.string.generator(model, size.pars, mig.pars,time.pars, use.alpha=use.alpha, scalar = as.numeric(loci[1,3]))
  #################################### theta and structure ###########################
  ######### generate -t and -I part of the command
  y <- paste(sum(as.numeric(model$I[1,4:ncol(model$I)])),1,paste(model$I[1,2:ncol(model$I)],collapse=" "),collapse=" ")
  ######### generate locfile part of the command
  loc.string <- paste("--frag-begin --finp .",arg,"locfile.txt --N ",Ne0," --frag-end",sep="")

  if(model$I[1, 3]=="1"){
    string <- gsub("-n 1","-en 0 1", string)
  }

  #### final command
  commands[[1]] <- paste(y,string,loc.string, collapse=" ")

  #### attach sampled parameters
  commands[[2]] <- t(parameters)
  return(commands)
  }