#setwd("/Users/nanettecoetzer/Documents/Bioinformatics/MAIZE_project/eQTL_pipeline/April_2012_pipeline_scripts/July_2012/output_July")

d.all <- read.table("all_hotspots_out.txt", header=T, sep="\t")
d.cis <- read.table("cis_hotspots_out.txt", header=T, sep="\t")
d.trans <- read.table("trans_hotspots_out.txt", header=T, sep="\t")

t <- read.table("full_summary.txt", header=F, sep="\t")
thresh.all = as.numeric(paste(t[12,2]))
thresh.cis = as.numeric(paste(t[15,2]))
thresh.trans = as.numeric(paste(t[18,2]))

chr_ends <- c()
for (i in 1:d.all[dim(d.all)[1],2]) {
	k = d.all[d.all[,2]==i,]
	chr_ends <- c(chr_ends,k[dim(k)[1],1])
}
num.chr = length(chr_ends)

prev.i = 0
chr_ends.mid = c()
for (i in chr_ends) {
	n = floor((i-prev.i)/2)	
	chr_ends.mid = c(chr_ends.mid,prev.i+n)
	prev.i = i
}

chrvector = c()
for (j in 1:num.chr) {
	chrvector = c(chrvector,paste("chr",j,sep=""))
}

excess.all = d.all[(d.all[,9]=='*'),]
excess.cis = d.cis[(d.cis[,9]=='*'),]
excess.trans = d.trans[(d.trans[,9]=='*'),]

##########################  all eQTL  #########################
# plot eQTL per cM --> all

m = sort(d.all[,6],TRUE)[1]
s = sort(d.all[,6],TRUE)[2]
if (m > (s+100)) {
	max_eqtl_freq = s
	}
if (m <= (s+100)) {
	max_eqtl_freq = m
	}

pdf("genome_wide_EcM_all.pdf",width=12)

par(mfrow=c(1,1))
plot(cbind(d.all[,1],d.all[,6]),type="h",ylab="eQTL per cM",main="Frequency plot",xlab="Sliding window bins (4 to 5.9 cM each)",lwd=1.5, ylim = c(-(0.03*m), m), col="black")

lines(cbind(excess.all[,1],excess.all[,6]),col="red",type="h",lwd=1.5)

dotted = c(0,chr_ends[1:length(chr_ends)]+0.5)
abline(h=thresh.all)
abline(v=dotted,lty="dotted",col="grey2",lwd=2)

text(chr_ends.mid,rep(-(max_eqtl_freq*0.03),num.chr),labels=chrvector,cex=0.8)

dev.off()

pdf("chr_plots_EcM_all.pdf",width=12)

if (num.chr <= 6) {
	par(mfrow=c(2,3))
	}
if (num.chr>6 & num.chr<=9) {
	par(mfrow=c(3,3))
	}
if (num.chr>9 & num.chr<=12) {
	par(mfrow=c(3,4))
	}
if (num.chr>12 & num.chr<=15) {
	par(mfrow=c(3,5))
	}
if (num.chr>15 & num.chr<=20) {
	par(mfrow=c(4,5))
	}
prev_k <- 0
#m = max_eqtl_freq
for (k in chr_ends) {
	nr = which(k == chr_ends)
	plot(cbind(d.all[(prev_k+1):k,1],d.all[(prev_k+1):k,6]),type="l",ylab="eQTL per cM",main=paste("Chr ",nr),xlab="Sliding window bins",lwd=1.5, ylim = c(-(0.03*s), s), col="black")
	dotted = c(prev_k,k+1)
	abline(h=thresh.all)
	abline(v=dotted,lty="dotted",col="grey2",lwd=2)
	
	sub_excess = excess.all[excess.all[,2]==nr,]
	lines(cbind(sub_excess[,1], sub_excess[,6]),col="red",type="h",lwd=1.5)
	prev_k <- k
}

dev.off()


##########################  cis eQTL  #########################
# plot eQTL per cM --> cis

m = sort(d.cis[,6],TRUE)[1]
s = sort(d.cis[,6],TRUE)[2]
if (m > (s+100)) {
	max_eqtl_freq = s
	}
if (m <= (s+100)) {
	max_eqtl_freq = m
	}

pdf("genome_wide_EcM_cis.pdf",width=12)

par(mfrow=c(1,1))
plot(cbind(d.cis[,1],d.cis[,6]),type="h",ylab="eQTL per cM",main="Frequency plot",xlab="Sliding window bins (4 to 5.9 cM each)",lwd=1.5, ylim = c(-(0.03*m), m), col="black")

lines(cbind(excess.cis[,1],excess.cis[,6]),col="red",type="h",lwd=1.5)

dotted = c(0,chr_ends[1:length(chr_ends)]+0.5)
abline(h=thresh.cis)
abline(v=dotted,lty="dotted",col="grey2",lwd=2)

text(chr_ends.mid,rep(-(max_eqtl_freq*0.03),num.chr),labels=chrvector,cex=0.8)

dev.off()

pdf("chr_plots_EcM_cis.pdf",width=12)

if (num.chr <= 6) {
	par(mfrow=c(2,3))
	}
if (num.chr>6 & num.chr<=9) {
	par(mfrow=c(3,3))
	}
if (num.chr>9 & num.chr<=12) {
	par(mfrow=c(3,4))
	}
if (num.chr>12 & num.chr<=15) {
	par(mfrow=c(3,5))
	}
if (num.chr>15 & num.chr<=20) {
	par(mfrow=c(4,5))
	}
prev_k <- 0
#m = max_eqtl_freq
for (k in chr_ends) {
	nr = which(k == chr_ends)
	plot(cbind(d.cis[(prev_k+1):k,1],d.cis[(prev_k+1):k,6]),type="l",ylab="eQTL per cM",main=paste("Chr ",nr),xlab="Sliding window bins",lwd=1.5, ylim = c(-(0.03*s), s), col="black")
	dotted = c(prev_k,k+1)
	abline(h=thresh.cis)
	abline(v=dotted,lty="dotted",col="grey2",lwd=2)
	
	sub_excess = excess.cis[excess.cis[,2]==nr,]
	lines(cbind(sub_excess[,1], sub_excess[,6]),col="red",type="h",lwd=1.5)
	prev_k <- k
}

dev.off()

##########################  trans eQTL  #########################
# plot eQTL per cM --> trans

m = sort(d.trans[,6],TRUE)[1]
s = sort(d.trans[,6],TRUE)[2]
if (m > (s+100)) {
	max_eqtl_freq = s
	}
if (m <= (s+100)) {
	max_eqtl_freq = m
	}

pdf("genome_wide_EcM_trans.pdf",width=12)

par(mfrow=c(1,1))
plot(cbind(d.trans[,1],d.trans[,6]),type="h",ylab="eQTL per cM",main="Frequency plot",xlab="Sliding window bins (4 to 5.9 cM each)",lwd=1.5, ylim = c(-(0.03*m), m), col="black")

lines(cbind(excess.trans[,1],excess.trans[,6]),col="red",type="h",lwd=1.5)

dotted = c(0,chr_ends[1:length(chr_ends)]+0.5)
abline(h=thresh.trans)
abline(v=dotted,lty="dotted",col="grey2",lwd=2)

text(chr_ends.mid,rep(-(max_eqtl_freq*0.03),num.chr),labels=chrvector,cex=0.8)

dev.off()

pdf("chr_plots_EcM_trans.pdf",width=12)

if (num.chr <= 6) {
	par(mfrow=c(2,3))
	}
if (num.chr>6 & num.chr<=9) {
	par(mfrow=c(3,3))
	}
if (num.chr>9 & num.chr<=12) {
	par(mfrow=c(3,4))
	}
if (num.chr>12 & num.chr<=15) {
	par(mfrow=c(3,5))
	}
if (num.chr>15 & num.chr<=20) {
	par(mfrow=c(4,5))
	}
prev_k <- 0
#m = max_eqtl_freq
for (k in chr_ends) {
	nr = which(k == chr_ends)
	plot(cbind(d.trans[(prev_k+1):k,1],d.trans[(prev_k+1):k,6]),type="l",ylab="eQTL per cM",main=paste("Chr ",nr),xlab="Sliding window bins",lwd=1.5, ylim = c(-(0.03*s), s), col="black")
	dotted = c(prev_k,k+1)
	abline(h=thresh.trans)
	abline(v=dotted,lty="dotted",col="grey2",lwd=2)
	
	sub_excess = excess.trans[excess.trans[,2]==nr,]
	lines(cbind(sub_excess[,1], sub_excess[,6]),col="red",type="h",lwd=1.5)
	prev_k <- k
}

dev.off()
