library(data.table)
small_indels<-read.table("~/gdr/assoc/pcaFiltered/indels/summary_stats/summary.indelReport.txt")
bed<-read.table("genes.bed")
big_dels<-read.table("big_indels.txt")
if(!exists("snps")){snps<-as.data.frame(fread("snps.window.txt"))}

plotLocus<-function(gene,window=5000){
	gene_start <- bed[match(gene,bed$V4),2]
	gene_end <- bed[match(gene,bed$V4),3]
	lw <- gene_start-window
	rw <- gene_end+window
	temp_bed <- subset(bed,V2>lw & V3<rw)
	xmax <- max(temp_bed$V3)
	xmin <- min(temp_bed$V2)
	xlims <- c(xmin,xmax)
	ylims <- c(0,3)
	plot(0,yaxt="n",type="n",xlim=xlims,ylim=ylims,ylab="Large Deletions",xlab="Genome position (Mb)")
	arrow_coords <- t(apply(temp_bed,1,function(x){if (x[5]=="+"){c(as.numeric(x[2]),as.numeric(x[3]))} else {c(as.numeric(x[3]),as.numeric(x[2]))}}))
	arrows(as.numeric(arrow_coords[,1]),c(0,0.5),arrow_coords[,2],c(0,0.5),lwd=5,col="blue",length=0.1)
	xmid <- (temp_bed$V2+temp_bed$V3)/2
	text(xmid,c(0.25,0.75),temp_bed$V4)

	temp_snps <- subset(snps, V1>lw & V1<rw)
	max_snps <- max(temp_snps$V2)
	points(scale(temp_snps$V2,F,max_snps)+2 ~ temp_snps$V1,type="l")
	points(scale(temp_snps$V3,F,max_snps)+2 ~ temp_snps$V1,type="l",col="dark green")
	points(small_indels$V2,rep(2,dim(small_indels)[1]),col="red",pch=4)
	abline(h=c(1,2))
	axis(side=2,at=scale(seq(0,round(max_snps,-1),50),F,max_snps)+2,labels=seq(0,round(max_snps,-1),50))
	mtext(text="SNPs/KB",side=2,line=3,at=2.5)
	mtext(text="Genes",side=2,line=3,at=0.5)
	temp_big_dels<-subset(big_dels,V3>lw-10000 & V4<rw+10000)
	num_dels <- dim(temp_big_dels)[1]
	if (num_dels>0){for (i in 1:num_dels){
			yh <- scale(i,F,num_dels+1)+1
			segments(temp_big_dels[i,3],yh,temp_big_dels[i,4],yh,col="dark green")
		}
	}
}

#plotLocus("pncA")