R/ms.commander.forplot.R

In gehara/PipeMaster: Building and Simulating Coalescent models.

# 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.commander.forplot<-function(model,use.alpha=use.alpha){

  # empty parameter vector
  parameters<-vector()

  # bind Ne, mig and Time priors
  size.pars<- data.frame(rbind(model$flags$n,model$flags$en$size))

  mig.pars<-data.frame(rbind(model$flags$m,model$flags$em$size))
  if(nrow(mig.pars)==0) mig.pars<-NULL

  time.pars<-data.frame(rbind(model$flags$ej,model$flags$en$time,model$flags$em$time))
  if(nrow(time.pars)==0) time.pars<-NULL

  options(warn=-1)
  size.pars[,4:5] <- t(apply(size.pars[,4:5],1,as.numeric))
  if(size.pars[,6]=="rnorm"){
    size.pars[,4:5] <- size.pars[,4]
    size.pars[,6] <- "runif"
  } else {size.pars[,4:5] <- apply(size.pars[,4:5],1, mean)}
  size.pars <- as.matrix(size.pars)

  if(is.null(mig.pars)==F){
  mig.pars[,4:5] <- t(apply(mig.pars[,4:5],1, as.numeric))
  if(mig.pars[,6]=="rnorm"){
    mig.pars[,4:5] <- mig.pars[,4]
    mig.pars[,6] <- "runif"
  } else {mig.pars[,4:5] <- apply(mig.pars[,4:5],1, mean)}
  mig.pars<-as.matrix(mig.pars)
  }

  if(is.null(time.pars)==F){
  time.pars[,4:5] <- t(apply(time.pars[,4:5],1,as.numeric))
  options(warn=-1)
  if(time.pars[,6]=="rnorm"){
    time.pars[,4:5] <- time.pars[,4]
    time.pars[,6] <- "runif"
  } else { time.pars[,4:5] <- apply(time.pars[,4:5],1,mean)}
  time.pars<-as.matrix(time.pars)
  }
  options(warn=0)
  # sample Ne, div.time and mutation rate
  size.pars <- 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<-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<-sample.pars(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<-rbind(parameters,loci[,c(1,4)])

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

  ####### Generate ms string ##############################################
  ####### Convertion to coalescent scale #####################################

  # generate coalescent scalar. Arbitrary value

  #### if single population
    Ne0<-1000000
    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<-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(parameters)
  return(commands)
  }