R/lsd-bootstrap.R

In OJWatson/sims4: Epidemic generational simualator

## fucntion declaration

LSD_Bootstrap <- function(t.exp=TRUE,N=1000,n.tips=50,n.bs=50,sigma.reps=seq(from=1,to=(N/n.tips)-0.5,by=0.5),
                          fasta.file="Data/1000.fasta",mu=0.0001,R0=2.4,Tg=50,
                          dir="C:/Users/Oliver/Google Drive/Academic Work/Imperial/O15-12/BEAST/Data Collection/LSD/",
                          lsd.dir="C:\\Users\\Oliver\\Desktop\\lsd"){
  
  ## HANDLE VARIABLES ##
  ##########################################################
  if(all(sigma.reps>=1)==FALSE) stop("sigma reps contains values less than 1")
  
  ## CREATE RESULTS STORAGE ##
  ##########################################################
  dir.create(path=dir)
  param.dir <- paste(dir,"mu_",mu,"_R0_",R0,"_Tg_",Tg,"/",sep="")
  dir.create(path=param.dir)
  res <- matrix(nrow = n.bs, ncol = length(sigma.reps))
  
  ## MAIN LOOP ##
  ##########################################################
  
  ## set up parameter loops
  for (i in 1:length(sigma.reps)){
    
    sub.dir <- paste(param.dir,"sigma",sigma.reps[i],"/",sep="")
    dir.create(path=sub.dir)
    
    if(i==1){
      ## simulate tree
      if(t.exp==TRUE){
        case.tree <- Time_Sim(fasta.file,mu=mu,R0=R0,Tg=Tg)
        save(case.tree,file=paste(param.dir,"casetree.RData",sep=""))
      } else {
        case.tree <- Generation_Sim(fasta.file=,mu=mu,R0=R0,Tg=Tg)
      }
      
      ## create the full tree for labelling reasons
      ft <- Case_To_Phylo(case.tree, full=TRUE)
    }
    
    ## set up bootstrap loops
    Optimised_Tip_Samples(case.tree=case.tree,file=sub.dir,n.tips=n.tips,
                          sigma=sigma.reps[i],bs=TRUE,n.bs=n.bs)
    
    ## create necessary LSD files
    ## create related trees
    for(j in 1:n.bs){
      nex <- ape::read.nexus.data(paste(sub.dir,"bs",j,".nex",sep=""))
      if(j==1){ 
        tree <- ape::write.tree(ape::njs(dist.dna(as.DNAbin.list(nex),model = "F84")),file=paste(sub.dir,"bs_tree.txt",sep=""))
        df <- matrix(nrow=N,ncol=2)
        df[1:N,1] <- c(ft$node.label,ft$tip.label)
        df[1:N,2] <- sapply(strsplit(c(ft$node.label,ft$tip.label), "_"), "[[", 3)
        write.table(df[,2],file=paste(sub.dir,"bs_date.date",sep=""),quote = FALSE,na = "",row.names = df[,1],col.names = N)
      } else {
      tree <- write.tree(njs(dist.dna(as.DNAbin.list(nex),model = "F84")),file=paste(sub.dir,"bs_tree.txt",sep=""),append=TRUE)
      }
    }
    
    ## write and execute LSD related script
    sh.sub.dir <- paste("\"",gsub(pattern = "/",replacement = "\\",fixed=TRUE,sub.dir),sep="")
    
    if(i==1){
      cmd.list <- list(cmd = paste(lsd.dir," -i ",sh.sub.dir,"bs_tree.txt\" -d ",sh.sub.dir,"bs_date.date\" -r a -n ",50,sep=""))
      cmd <- paste(lsd.dir," -i ",sh.sub.dir,"bs_tree.txt\" -d ",sh.sub.dir,"bs_date.date\" -c -r a -n ",50,sep="")
    } else {
      cmd <- paste(cmd,"\n",lsd.dir," -i ",sh.sub.dir,"bs_tree.txt\" -d ",sh.sub.dir,"bs_date.date\" -c -r a -n ",50,sep="")
      cmd.list$cmd[i] <- paste(lsd.dir," -i ",sh.sub.dir,"bs_tree.txt\" -d ",sh.sub.dir,"bs_date.date\" -r a -n ",50,sep="")
    }
    ## colate and plotbootstrap error against sigma
    
  
    
    ## OLD COMMAND CODE THAT KEPT THROWING ERROR 5 CHECK ##
    #############################################################
    #cmd <- paste(lsd.dir," -i ",sh.sub.dir,"bs_tree.txt\" -d ",sh.sub.dir,"bs_date.date\" -c -r as -n ",n.bs,sep="")
    #system(cmd)
    ## colate and plotbootstrap error against sigma
    
    #lsd.tb <- read.table(file=paste(sub.dir,"bs_tree.txt.result",sep=""),skip = 20,nrows = n.bs)
    #res[,i] <- lsd.tb[,6]
    
  }
  
  ## BIG UGLY HACK AS THE PROGRAM KEEPS TIMING OUT OR SOMETHING ##
  ################################################################
  fileConn2<-file("C:/Users/Oliver/Desktop/lsd.bat")
  writeLines(cmd, fileConn2)
  close(fileConn2)
  system("C:\\Users\\Oliver\\Desktop\\lsd.bat")
  
  for (i in 1:length(sigma.reps)){
    sub.dir <- paste(param.dir,"sigma",sigma.reps[i],"/",sep="")
    info = file.info(paste(sub.dir,"bs_tree.txt.result",sep=""))
    if (info$size == 0 | is.na(info$size)){
      system(cmd.list$cmd[i])
    }
    lsd.tb <- read.table(file=paste(sub.dir,"bs_tree.txt.result",sep=""),skip = 20,nrows = n.bs)
    res[,i] <- lsd.tb[,6]
  }
  
  boxplot(res)
  return(res)
}