R/create_params_colonyV2.R

In zakrobinson/RSibPurge: A convenient R-interface for COLONY Version 2.0.6.5

#' Create COLONY, Version 2.0.6.5 Input files
#'
#' Create input file for Colony, Version 2.0.6.5. Author Of COLONY: Jinliang Wang, Institute of Zoology, Zoological Society of London
#'
#'@param projectname a text string that describes the Colony run "project".
#'@param project_path a directory where input and output files will be stored.
#'@param offspring_alleles a data frame with individual IDs in the first column and alleles or genotypes following.
#'@param geno Logical if offspring_alleles are in genotype per column format versus allele per column format.
#'@param nCode if geno=T then how many digits are the genotypes coded in (i.e. 12 heterzygote has a nCode of 2)
#'@param rand_seed Seed for random number generator
#'@param update_allelefreq an indicator of value 1 (or 0) to instruct Colony to (or not to) update allele frequencies during the simulated annealing process in searching for the ML configuration.
#'@param Dio_Mono a value of either 2 or 1 to indicate dioecious species or monoecious species.
#'@param inbred a value of either 0 or 1 to indicate the absence or presence of inbreeding.
#'@param haplodiploid a value of either 0 or 1 to indicate diploid species or haplodiploid species.
#'@param monogamy_M an indicator value (0 or 1) to specify whether males are polygamous (=0) or monogamous (=1).
#'@param monogamy_F an indicator value (0 or 1) to specify whether females are polygamous (=0) or monogamous (=1).
#'@param clones an indicator value (1 or 0) to specify whether clones (or duplicated individuals) are to be inferred (=1) or not (=0).
#'@param FS_scaling an indicator value (1 or 0) to specify whether full sibship size is to be scaled (=1) or not (=0).
#'@param FS_prior_M give an indicator value of 0, 1, 2 or 3 to indicate no sibship prior, weak sibship prior, medium sibship prior or strong sibship prior for males.
#'@param FS_prior_F give an indicator value of 0, 1, 2 or 3 to indicate no sibship prior, weak sibship prior, medium sibship prior or strong sibship prior for females.
#'@param nRuns Provide a value for the number of replicate runs for the dataset (project)
#'@param RunLength Give a value of 1, 2, 3, 4 to indicate short, medium, long, very long run
#'@param likelihood Give an indicated value of 0, 1 or 2 to tell Colony to use the pairwise-likelihood score (PLS), full likelihood (FL), or the FL and PLS combined (FPLS) method.
#'@param precision Give an indicated value of 0/1/2/3 to use Low/Medium/High/VeryHigh precision in calculating the full likelihood.
#'@param alle_drop allele dropout rate.
#'@param fasle_alle other typing error rate.
#'@author Zak Robinson, Contact: zachary.robinson(at)umontana.com
#'@return returns path to parameter output file. Parameter output file is written to disk.
#'@export



create_params_colonyV2<-function(projectname="test1",project_path=getwd(), offspring_alleles,geno=T,nCode=NULL, rand_seed=1234,update_allelefreq=0,Dio_Mono=2,inbred=0,haplodiploid = 0, monogamy_M=0, monogamy_F=0, clones=0, FS_scaling=1,FS_prior=0,FS_prior_M=NA,FS_prior_F=NA,nRuns=1,RunLength=2,likelihood=1,precision=3, alle_drop=0,false_alle=0.005){

## input and output filename generation
project_path<-gsub(x = project_path,pattern = "(.+)/$",replacement = "\\1")
out_inputfile_name=paste0(project_path,"/",projectname,"_Colony_params",".txt")
out_outputfile=paste0(project_path,"/",projectname)

offspring_alleles<-apply(offspring_alleles,2,as.character)

if(geno==T){
  offspring_alleles<-cbind(offspring_alleles[,1],genotypes2alleles(offspring_alleles[,-1],ncode = nCode))
   }



fix_off_alleles<-function(xx){
    return(format(x=xx,width = max(nchar(xx))+2,justify = "right"))
  }

offspring_alleles[is.na(offspring_alleles)] <- 0
  on.exit(expr = suppressWarnings(sink()))
  sink(file = out_inputfile_name,append = F)
  cat("'",projectname,"'","     ! Dataset Name",sep="")
  cat("\n")
  cat("'",out_outputfile,"'","    ! Output File Name",sep="")
  cat("\n")
  cat(as.character(nrow(offspring_alleles)),"    ! Number of offspring genotypes",sep="")
  cat("\n")
  cat(as.character((ncol(offspring_alleles)-1)/2),"    ! Number of loci",sep="")
  cat("\n")
  cat(as.character(rand_seed),"    ! Seed for random number generator",sep="")
  cat("\n")
  cat(as.character(update_allelefreq),"     ! 0/1=Not updating/updating allele frequency",sep="")
  cat("\n")
  cat(as.character(Dio_Mono),"    ! 2/1=Dioecious/Monoecious species",sep="")
  cat("\n")
  cat(as.character(inbred),"     ! 0/1=No inbreeding/inbreeding",sep="")
  cat("\n")
  cat(as.character(haplodiploid),"     ! 0/1=Diploid species/HaploDiploid species",sep="")
  cat("\n")
  cat(as.character(monogamy_M),"  ",as.character(monogamy_F),"     ! 0/1=Polygamy/Monogamy for males & females",sep="")
  cat("\n")
  cat(as.character(clones),"     ! 0/1=Clone inference =No/Yes",sep="")
  cat("\n")
  cat(as.character(FS_scaling),"    ! 0/1=Full sibship size scaling =No/Yes",sep="")
  cat("\n")

  if(FS_prior==0){
  cat(as.character(FS_prior),"     ! 0,1,2,3=No,weak,medium,strong sibship size prior; mean paternal & meteral sibship size")
  }
  else{
  cat(paste(as.character(FS_prior),as.character(FS_prior_M),as.character(FS_prior_F)),"      ! 0,1,2,3=No,weak,medium,strong sibship size prior; mean paternal & meteral sibship size")
  }

  cat("\n")
  cat("0","      ! 0/1=Unknown/Known population allele frequency")
  cat("\n")
  cat(as.character(nRuns),"     ! Number of runs")
  cat("\n")
  cat(as.character(RunLength),"     ! Length of run")
  cat("\n")
  cat("0","      ! 0/1=Monitor method by Iterate#/Time in second")
  cat("\n")
  cat("100000","     ! Monitor interval in Iterate# / in seconds")
  cat("\n")
  cat("0","     ! non-Windows version")
  cat("\n")
  cat(as.character(likelihood),"     ! 0/1/2=PairLikelihood score/Fulllikelihood/FPLS")
  cat("\n")
  cat(as.character(precision),"      ! 1/2/3=low/medium/high Precision for Fulllikelihood")
  cat("\n\n")
  cat(format(paste("Loc", seq(1,(ncol(offspring_alleles)-1)/2,1),sep = ""),width = 8,justify = "right"),"        !Marker names\n")
  cat(format(rep(0,(ncol(offspring_alleles)-1)/2),width = 8,justify = "right"),"        !Marker types, 0/1 = codominant/dominant\n")
  cat(format(rep(alle_drop,(ncol(offspring_alleles)-1)/2),nsmall = 4,width = 8,justify = "right",scientific = F),"        !Allelic dropout rate\n")
  cat(format(rep(false_alle,(ncol(offspring_alleles)-1)/2),nsmall= 4,width = 8,justify = "right",scientific = F),"        !False allele rate\n\n\n")
  tp<-apply(offspring_alleles, MARGIN = 2,FUN = fix_off_alleles)
  apply(tp,MARGIN = 1,FUN = function(x){cat(x,"\n",sep = "")})
  cat("\n\n")
  cat("0  0     !prob. of dad/mum included in the candidates\n0  0       !numbers of candiadte males & females\n\n\n0  0          !#known fater-offspring dyads, paternity exclusion threshold\n\n0  0          !#known moter-offspring dyads, maternity exclusion threshold\n\n\n0             !#known paternal sibship with unknown fathers\n\n\n\n0             !#known maternal sibship with unknown mothers\n\n\n0             !#known paternity exclusions\n\n\n0             !#known maternity exclusions\n\n\n\n0             !#known paternal sibship exclusions\n\n\n0             !#known maternal sibship exclusions\n")
  sink()
  return(out_inputfile_name)

}