stephane_plot

library(scales)
library(data.table)
setwd("~/projects/collaberations/stephane/")
rawdata<-as.data.frame(fread("test.txt",header = T,stringsAsFactors = F))
head(rawdata)

plot_vcf<-function(rawdata,start,end){
  if (start==0 & end==0){
    dat<-rawdata
  } else {
    dat<-subset(rawdata,pos>=start & pos<=end)
  }
  samples<-unique(dat$sample)
  if (length(samples)<2){stop("Not enough variants in window for one of the samples")}
  dat.s1<-subset(dat,sample==samples[1])
  dat.s2<-subset(dat,sample==samples[2])
  
  x.max<-max(dat$pos)
  x.min<-0
  y.max<-2
  y.min<-0
  plot(1,1,type="n",xlim=c(x.min,x.max),ylim=c(y.min,y.max),yaxt="n",ylab="",xlab="Genome Position")
  mtext(text = samples[1],side = 2,line = 1,at = 0.5)  
  mtext(text = samples[2],side = 2,line = 1,at = 1.5) 
  
  ybot.s1<-0
  ytop.s1<-ybot.s1+dat.s1$freq
  col.s1<-sapply(dat.s1$freq,function(x){rgb(x,0,1-x)})
  segments(x0 = dat.s1$pos, y0 = ybot.s1, x1 = dat.s1$pos, y1 = ytop.s1,col=col.s1)
  
  ybot.s2<-1
  ytop.s2<-ybot.s2+dat.s2$freq
  col.s2<-sapply(dat.s2$freq,function(x){rgb(x,0,1-x)})
  
  segments(x0 = dat.s2$pos, y0 = ybot.s2, x1 = dat.s2$pos, y1 = ytop.s2,col=col.s2)
}  

lowfreq<-subset(rawdata,freq<0.7)
#pdf("test.pdf")
plot_vcf(lowfreq,0,0)
#dev.off()