###manhattan plots in qqplot2
library(ggplot2)
setwd("~/ncbi/zm/XPCLR/")
read.table("LW.gene.xpclr.position.txt",header=F,sep=',')->IL
head(IL)
#V1       V2 V3     V4
#1 GRMZM2G059865 1.159389  1   4854
#2 GRMZM5G888250 1.159389  1   9882

ILp<-data.frame(ID=1:dim(IL)[1], chr=factor(IL[,3]), coord=IL[,4], p=IL[,2])

head(ILp)
#ID chr  coord        p
#1  1   1   4854 1.159389
#2  2   1   9882 1.159389

##Position is defined as position on chromosome. Need to add largest position from last marker on previous chr so that markers are ordered correctly along x-axis
chrNum=10
for(i in 1:chrNum){
  ndx<-which(ILp[,2]==i)
  lstMrk<-max(ILp[ndx,3])
  if (i < chrNum) ndx2 <- which(ILp[, 2]==i+1)
  if (i < chrNum) ILp[ndx2, 3] <- ILp[ndx2,3] + lstMrk
}
##Use this little loop to find the midposition of every chromosome so nice chr labels can be added to the plot
bpMidVec <- vector(length=chrNum)
for (i in 1:chrNum){
  ndx <- which(ILp[, 2]==i)
  posSub <- ILp[ndx, 3]
  bpMidVec[i] <- ((max(posSub) - min(posSub))/2) + min(posSub)
}
##Use qplot function in ggplot2 to create Manhattan plot
cbPalette <- c("darkgrey", "#000000","#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","red")

tiff(filename = "LW.tif",width = 36,height = 18,units ="cm",compression="lzw",bg="white",res=300)
ggplot(ILp) + geom_point(aes(x=coord, y=p,colour=chr),show.legend = FALSE) + scale_color_manual(values=cbPalette)+#brewer(type="seq",palette="Set5")+
  scale_x_continuous(labels=as.character(1:chrNum),breaks=bpMidVec)+theme_bw(base_size=15)+labs(x="Chromosomes",y="XP-CLR score of Landraces and Wild")+ylim(0,152)
dev.off()