R/lsd-beast-plot.R

In OJWatson/sims4: Epidemic generational simualator

LSD_BEAST_Plot <- function(num=c(4,5,6)){
  
  par(mfrow=c(length(num),1))
  
  for (j in 1:length(num)){
    
    load(paste("C:/Users/Oliver/Google Drive/Academic Work/Imperial/O15-12/BEAST/Data Collection/LSD",num[j],"/mu_1e-04_R0_2.4_Tg_50/casetree.RData",sep=""))
    
    dir=paste("C:/Users/Oliver/Google Drive/Academic Work/Imperial/O15-12/BEAST/Data Collection/LSD",num[j],"/",sep="")
    #dir.create(path=dir)
    network <- paste(network.path,tail(strsplit(dir,"/")[[1]],1),"\\",sep="")
    #dir.create(path=network)
    network.dir <- paste(network.path,tail(strsplit(dir,"/")[[1]],1),"\\","mu_",mu,"_R0_",R0,"_Tg_",Tg,"\\",sep="")
    #dir.create(path=network.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))
    m.res <- matrix(nrow = n.bs, ncol = length(sigma.reps))
    
    ## import beast
    network.root.dir=paste("LSD",num[j],"\\mu_1e-04_R0_2.4_Tg_50\\",sep="")
    local.root.dir=paste("LSD",num[j],"/mu_1e-04_R0_2.4_Tg_50/",sep="")
    beast.res <- Network_Pull_Tracer(network.root.dir = network.root.dir, local.root.dir = local.root.dir)
    root <- unique(beast.res$E_TMRCA)
    
    
    ## set up parameter loops
    for (i in 1:length(sigma.reps)){
      
      ## create dirs for results
      sub.dir <- paste(param.dir,"sigma",sigma.reps[i],"/",sep="")
      #dir.create(path=sub.dir)
      network.sub.dir <- paste(network.dir,"sigma",sigma.reps[i],"\\",sep="")
      #dir.create(path=network.sub.dir)
      
      ## write and execute LSD related script
      #sh.sub.dir <- paste("\"",gsub(pattern = "/",replacement = "\\",fixed=TRUE,sub.dir),sep="")
      #cmd <- paste(lsd.dir," -i ",sh.sub.dir,"bs_tree.txt\" -d ",sh.sub.dir,"bs_date.date\" -r a -n ",n.bs,sep="")
      #system
      
      ## 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]
      m.res[,i] <- as.numeric(substr(lsd.tb[,4], 1, nchar(as.character(lsd.tb[,4]))-1))
      
    }
    boxplot(res,xlab = "sigma",ylab = "TMRCA Error",ylim=c(-100, 100))
    par(new=T)
    for (i in 1:length(sigma.reps)){  
      points(x=c(i,i),y=c(root-beast.res$O_TMRCA_95L[i],root-beast.res$O_TMRCA_95H[i]), type = "l", col = "red",lwd = 2)
      points(x=i,y=root-beast.res$O_TMRCA_Med[i], col = "green",lwd = 2)  
      
    }
    
    abline(h=case.tree$Time[1],col="blue")
    par(new=F)
  }
  
}