rm(list=ls())
require(ggplot2)
require(reshape)
require(grid)
require(gridExtra)
library(dplyr)
library(tidyr)
library(ggplot2); theme_set(theme_bw())
library(windowscanr)

#load("~/Documents/Data/ddRAD/range_wide_r3/figure3_final_G_trial3.RData")
# used for cross G (f2ld)
# used to plot fst, dxy, pi and tajima's d ...

#### read the input table ####
cross1 <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)
names(cross1)
cross <- cross1[with(cross1,order(LG,Position)),]

### fst between####

bet <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/fst/hier/hier_fst_bet.tsv",sep="\t",header=T)

dim(bet)

ids <- c()
chrs <- c()
bps <- c()
fsts <- c()
for ( i in 1:nrow(cross) ) {
  print(i)
  contig <- cross[i,"Contig"]
  contig_pos <- cross[i,"Contig_Pos"]
  lg <- cross[i,"LG"]
  lg_pos <- cross[i,"Position"]
  index <- which(bet$CHROM %in% contig)
  count <- 1
  for (j in index) {
    id <- paste(bet[j,"CHROM"],bet[j,"POS"])
    ids <- append(ids,id)
    chrs <- append(chrs,lg)
    position <- lg_pos + count * 0.011 
    bps <- append(bps,position)
    fsts <- append(fsts,bet[j,"WEIR_AND_COCKERHAM_FST"])
    count <- count + 1
  }
}

fst.bet.table <- data.frame(id=ids,chr=chrs,bp=bps,fst=fsts)
dim(fst.bet.table)

### qtl table ####
table <- read.table("~/Dropbox/Xiaodong PhD thesis/Chapter 3 Population Genomics/Analysis/qtl_G.csv",sep=",",header=T)
table2 <- table[with(table,order(Chr,Pos)),]
table3 <- table2


### fst with species, from amova ####
wth <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/fst/hier/hier_fst_wth.txt",sep="\t",header=T)
dim(wth)

ids <- c()
chrs <- c()
bps <- c()
fsts <- c()
for ( i in 1:nrow(cross) ) {
  print(i)
  contig <- cross[i,"Contig"]
  contig_pos <- cross[i,"Contig_Pos"]
  lg <- cross[i,"LG"]
  lg_pos <- cross[i,"Position"]
  index <- which(wth$CHROM %in% contig)
  count <- 1
  for (j in index) {
    id <- paste(bet[j,"CHROM"],wth[j,"POS"])
    ids <- append(ids,id)
    chrs <- append(chrs,lg)
    position <- lg_pos + count * 0.011 
    bps <- append(bps,position)
    fsts <- append(fsts,wth[j,"WEIR_AND_COCKERHAM_FST"])
    count <- count + 1
  }
}

fst.wth.table <- data.frame(id=ids,chr=chrs,bp=bps,fst=fsts)
dim(fst.wth.table)


### fst with SD ####
setwd("/home/ivan/Documents/Data/ddRAD/range_wide_r3/fst/wthSpecies/sd/")
files<-list.files(pattern="*.fst")
pairs <- gsub(".weir.fst","",files)

fst <- c()
colnames <- c()
for (i in 1:length(pairs)) {
  j <- pairs[i]
  print(j)
  t <-read.table ( file=paste(j,"weir.fst",sep=".") ,header=T,sep="\t" ) 
  colnames <- append(colnames, paste(j,"fst",sep=".") )
  fst <- append(fst,t[,3])
}

fst.matrix <-matrix(fst,ncol=length(pairs),byrow = F)
colnames(fst.matrix) <- colnames
fst.matrix.2 <- cbind(t[,1:2],fst.matrix)
#na.miss <- apply(fst.matrix.2, 1, function(x) sum(is.na(x)))
#na.index <- which(na.miss >= 45)
#fst.matrix.3 <- fst.matrix.2[-na.index,]
fst.matrix.3 <- fst.matrix.2

fst.mean <- apply(fst.matrix.3[,3:ncol(fst.matrix.3)],1,function(x) mean(x,na.rm=T))
fst.median <- apply(fst.matrix.3[,3:ncol(fst.matrix.3)],1,function(x) median(x,na.rm=T))
fst.max <- apply(fst.matrix.3[,3:ncol(fst.matrix.3)],1,function(x) max(x,na.rm=T))

wth.SD <- cbind(fst.matrix.3[,1:2],fst.max)
inf.index <- which(wth.SD$fst.max == "-Inf")
wth.SD <- wth.SD[-inf.index,]

ids <- c()
chrs <- c()
bps <- c()
fsts <- c()
for ( i in 1:nrow(cross) ) {
  print(i)
  contig <- cross[i,"Contig"]
  contig_pos <- cross[i,"Contig_Pos"]
  lg <- cross[i,"LG"]
  lg_pos <- cross[i,"Position"]
  index <- which(wth.SD$CHROM %in% contig)
  count <- 1
  for (j in index) {
    id <- paste(wth.SD[j,"CHROM"],wth.SD[j,"POS"])
    ids <- append(ids,id)
    chrs <- append(chrs,lg)
    position <- lg_pos + count * 0.011 
    bps <- append(bps,position)
    fsts <- append(fsts,wth.SD[j,"fst.max"])
    count <- count + 1
  }
}

fst.wth.SD.table <- data.frame(id=ids,chr=chrs,bp=bps,fst=fsts)
head(fst.wth.SD.table)



### fst with SL####
setwd("/home/ivan/Documents/Data/ddRAD/range_wide_r3/fst/wthSpecies/sl/")

files<-list.files(pattern="*.fst")
pairs <- gsub(".weir.fst","",files)

fst <- c()
colnames <- c()
for (i in 1:length(pairs)) {
  j <- pairs[i]
  print(j)
  t <-read.table ( file=paste(j,"weir.fst",sep=".") ,header=T,sep="\t" ) 
  colnames <- append(colnames, paste(j,"fst",sep=".") )
  fst <- append(fst,t[,3])
}

fst.matrix <-matrix(fst,ncol=length(pairs),byrow = F)
colnames(fst.matrix) <- colnames
fst.matrix.2 <- cbind(t[,1:2],fst.matrix)
#na.index <- which(na.miss >= 27)
#na.miss <- apply(fst.matrix.2, 1, function(x) sum(is.na(x)))
#fst.matrix.3 <- fst.matrix.2[-na.index,]

fst.matrix.3 <- fst.matrix.2
fst.mean <- apply(fst.matrix.3[,3:ncol(fst.matrix.3)],1,function(x) mean(x,na.rm=T))
fst.median <- apply(fst.matrix.3[,3:ncol(fst.matrix.3)],1,function(x) median(x,na.rm=T))
fst.max <- apply(fst.matrix.3[,3:ncol(fst.matrix.3)],1,function(x) max(x,na.rm=T))

wth.SL <- cbind(fst.matrix.3[,1:2],fst.max)
inf.index <- which(wth.SL$fst.max == "-Inf")
wth.SL <- wth.SL[-inf.index,]



ids <- c()
chrs <- c()
bps <- c()
fsts <- c()
for ( i in 1:nrow(cross) ) {
  print(i)
  contig <- cross[i,"Contig"]
  contig_pos <- cross[i,"Contig_Pos"]
  lg <- cross[i,"LG"]
  lg_pos <- cross[i,"Position"]
  index <- which(wth.SL$CHROM %in% contig)
  count <- 1
  for (j in index) {
    id <- paste(wth.SL[j,"CHROM"],wth.SL[j,"POS"])
    ids <- append(ids,id)
    chrs <- append(chrs,lg)
    position <- lg_pos + count * 0.011 
    bps <- append(bps,position)
    fsts <- append(fsts,wth.SL[j,"fst.max"])
    count <- count + 1
  }
}

fst.wth.SL.table <- data.frame(id=ids,chr=chrs,bp=bps,fst=fsts)
head(fst.wth.SL.table)

### highest fst between SD populations ####


bet <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/fst/wthSpecies/sd/FO1_POL1.weir.fst",sep="\t",header=T)
dim(bet)


ids <- c()
chrs <- c()
bps <- c()
fsts <- c()
for ( i in 1:nrow(cross) ) {
  print(i)
  contig <- cross[i,"Contig"]
  contig_pos <- cross[i,"Contig_Pos"]
  lg <- cross[i,"LG"]
  lg_pos <- cross[i,"Position"]
  index <- which(bet$CHROM %in% contig)
  count <- 1
  for (j in index) {
    id <- paste(bet[j,"CHROM"],bet[j,"POS"])
    ids <- append(ids,id)
    chrs <- append(chrs,lg)
    position <- lg_pos + count * 0.011 
    bps <- append(bps,position)
    fsts <- append(fsts,bet[j,"WEIR_AND_COCKERHAM_FST"])
    count <- count + 1
  }
}

SD.max.table <- data.frame(id=ids,chr=chrs,bp=bps,fst=fsts)
dim(SD.max.table)

### highest fst between SL populations ####

bet <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/fst/wthSpecies/sl/ESP1_RUS1.weir.fst",sep="\t",header=T)
dim(bet)


ids <- c()
chrs <- c()
bps <- c()
fsts <- c()
for ( i in 1:nrow(cross) ) {
  print(i)
  contig <- cross[i,"Contig"]
  contig_pos <- cross[i,"Contig_Pos"]
  lg <- cross[i,"LG"]
  lg_pos <- cross[i,"Position"]
  index <- which(bet$CHROM %in% contig)
  count <- 1
  for (j in index) {
    id <- paste(bet[j,"CHROM"],bet[j,"POS"])
    ids <- append(ids,id)
    chrs <- append(chrs,lg)
    position <- lg_pos + count * 0.011 
    bps <- append(bps,position)
    fsts <- append(fsts,bet[j,"WEIR_AND_COCKERHAM_FST"])
    count <- count + 1
  }
}

SL.max.table <- data.frame(id=ids,chr=chrs,bp=bps,fst=fsts)
dim(SL.max.table)

### pi wth SD####
pi.sd <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/pi/sd/sd.sites.pi",sep="\t",header=T)
contig_length <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/pi/interval/length.txt",sep="\t",header=T)
names(pi.sd)


index <-contig_length$Chr %in% unique(pi.sd$CHROM)
contig_length.subset <- contig_length[index,]
names(contig_length.subset) <- c("contig","length")

contig.sum <- c()
pi.sum <- c()
for (count in 1:nrow(contig_length.subset)) {
  
  i<- contig_length.subset[count,1]
  
  index <- which(pi.sd$CHROM %in% i)
  pi <- sum(pi.sd[index,"PI"],na.rm=T)  
  
  
  pi2 <- pi / contig_length.subset[count,2]
  pi.sum <- append(pi.sum,pi2)
  contig.sum <- append(contig.sum,as.character(i))
}

pi.table.sd.all <- data.frame(contigs=contig.sum,pi=pi.sum)
head(pi.table.sd.all)
G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)
G.pi.sd.table <- merge(G,pi.table.sd.all, by.x="Contig",by.y="contigs")
G.pi.sd.table.2 <-G.pi.sd.table[ with(G.pi.sd.table,order(LG,Position)),]
G.pi.sd.table.final <- G.pi.sd.table.2[,c(1,4,5,14)]
names(G.pi.sd.table.final) <- c("contig","chr","bp","fst")
head(G.pi.sd.table.final)

### pi wth SL ####
pi.sl <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/pi/sl/sl.sites.pi",sep="\t",header=T)
contig_length <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/pi/interval/length.txt",sep="\t",header=T)
names(pi.sl)


index <-contig_length$Chr %in% unique(pi.sl$CHROM)
contig_length.subset <- contig_length[index,]
names(contig_length.subset) <- c("contig","length")
contig.sum <- c()
pi.sum <- c()
for (count in 1:nrow(contig_length.subset)) {
  
  i<- contig_length.subset[count,1]
  
  index <- which(pi.sl$CHROM %in% i)
  pi <- sum(pi.sl[index,"PI"],na.rm=T)  
  
  
  pi2 <- pi / contig_length.subset[count,2]
  pi.sum <- append(pi.sum,pi2)
  contig.sum <- append(contig.sum,as.character(i))
}

pi.table.sl.all <- data.frame(contigs=contig.sum,pi=pi.sum)
dim(pi.table.sl.all)
head(pi.table.sl.all)

G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old//G_map_new.csv",sep=",",header=T)
G.pi.sl.table <- merge(G,pi.table.sl.all, by.x="Contig",by.y="contigs")
G.pi.sl.table.2 <- G.pi.sl.table[ with(G.pi.sl.table,order(LG,Position)),]

G.pi.sl.table.final <- G.pi.sl.table.2[,c(1,4,5,14)]
names(G.pi.sl.table.final) <- c("contig","chr","bp","fst")

### dxy ####
dxy <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/dxy/Dxy.tsv",sep="\t",header=F)
contig_length <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/pi/interval/length.txt",sep="\t",header=T)
names(dxy) <- c("CHROM", "POS" ,"DXY"  )


index <- contig_length$Chr %in% unique(dxy$CHROM)
contig_length.subset <- contig_length[index,]
names(contig_length.subset) <- c("contig","length")

contig.sum <- c()
dxy.sum <- c()
for (count in 1:nrow(contig_length.subset)) {
  
  i<- contig_length.subset[count,1]
  
  index <- which(dxy$CHROM %in% i)
  pi <- sum(dxy[index,"DXY"],na.rm=T)  
  
  
  pi2 <- pi / contig_length.subset[count,2]
  dxy.sum <- append(dxy.sum,pi2)
  contig.sum <- append(contig.sum,as.character(i))
}

dxy.table.all <- data.frame(contigs=contig.sum,dxy=dxy.sum)
head(dxy.table.all)


G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)
G.dxy.table <- merge(G,dxy.table.all, by.x="Contig",by.y="contigs")
G.dxy.table[ with(G.dxy.table,order(LG,Position)),]
G.dxy.table.2 <- G.dxy.table[ with(G.dxy.table,order(LG,Position)),]


G.dxy.table.final <- G.dxy.table.2[,c(1,4,5,14)]
names(G.dxy.table.final) <- c("contig","chr","bp","fst")
head(G.dxy.table.final)

### hfst contig based ####
hfstC.table.all <- read.table("~/Documents/Data/ddRAD/range_wide_r3/fst/ContigBased/contigFst.csv",sep=",",header=T)
head(hfstC.table.all)
( threshold <- quantile(hfstC.table.all$fst,0.95,na.rm=T) )
G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)

G.hfstC.table <- merge(G,hfstC.table.all, by.x="Contig",by.y="chr")
G.hfstC.table[ with(G.hfstC.table,order(LG,Position)),]
G.hfstC.table.2 <- G.hfstC.table[ with(G.hfstC.table,order(LG,Position)),]


hfstC.table <- G.hfstC.table.2[,c(1,4,5,14)]
names(hfstC.table) <- c("contig","chr","bp","fst")
head(hfstC.table)

( island_contigs <- hfstC.table[which(hfstC.table$fst >= threshold),"contig"] )


## highest population within sd
hfstC.SD.table.all <- read.table("~/Documents/Data/ddRAD/range_wide_r3/fst/wthSpecies/contigFst_highest_pop/contigFst_highestSD.tsv",sep="\t",header=T)
head(hfstC.SD.table.all)

G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)

G.hfstC.SD.table <- merge(G,hfstC.SD.table.all, by.x="Contig",by.y="chr")
G.hfstC.SD.table[ with(G.hfstC.SD.table,order(LG,Position)),]
G.hfstC.SD.table.2 <- G.hfstC.SD.table[ with(G.hfstC.SD.table,order(LG,Position)),]


hfstC.SD.table <- G.hfstC.SD.table.2[,c(1,4,5,14)]
names(hfstC.SD.table) <- c("contig","chr","bp","fst")
head(hfstC.SD.table)

## highest population within sl
hfstC.SL.table.all <- read.table("~/Documents/Data/ddRAD/range_wide_r3/fst/wthSpecies/contigFst_highest_pop/contigFst_highestSL.tsv",sep="\t",header=T)
head(hfstC.SL.table.all)

G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)

G.hfstC.SL.table <- merge(G,hfstC.SL.table.all, by.x="Contig",by.y="chr")
G.hfstC.SL.table[ with(G.hfstC.SL.table,order(LG,Position)),]
G.hfstC.SL.table.2 <- G.hfstC.SL.table[ with(G.hfstC.SL.table,order(LG,Position)),]


hfstC.SL.table <- G.hfstC.SL.table.2[,c(1,4,5,14)]
names(hfstC.SL.table) <- c("contig","chr","bp","fst")
head(hfstC.SL.table)

# with species, between pop
hfstC_wth.all.table <- read.table("~/Documents/Data/ddRAD/range_wide_r3/fst/ContigBased/contigFst_wthspecies.csv",sep=",",header=T)
head(hfstC_wth.all.table)
W <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/W_map_new.csv",sep=",",header=T)
W.hfstC_wth.table <- merge(W,hfstC_wth.all.table, by.x="Contig",by.y="chr")
W.hfstC_wth.table[ with(W.hfstC_wth.table,order(LG,Position)),]
W.hfstC_wth.table.2 <- W.hfstC_wth.table[ with(W.hfstC_wth.table,order(LG,Position)),]
hfstC_wth.table <- W.hfstC_wth.table.2[,c(4,5,6)]
names(hfstC_wth.table) <- c("chr","bp","fst")
head(hfstC_wth.table)


### SD Tajima's D ####
sd.tajimaD <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/TajimaD/SD.tajimaD.out",sep="\t",header=F)
head(sd.tajimaD)
sd.tajiD <- sd.tajimaD[,c(1,7)]
names(sd.tajiD) <- c("contigs","TajimaD")
head(sd.tajiD)

G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)
G <- G[,c(1:5)]
G.sd.tajiD.table <- merge(G,sd.tajiD, by.x="Contig",by.y="contigs")
G.sd.tajiD.table[ with(G.sd.tajiD.table,order(LG,Position)),]
G.sd.tajiD.table.2 <- G.sd.tajiD.table[ with(G.sd.tajiD.table,order(LG,Position)),]


G.sd.tajiD.table.final <- G.sd.tajiD.table.2[,c(1,4,5,6)]
names(G.sd.tajiD.table.final) <- c("contig","chr","bp","fst")
head(G.sd.tajiD.table.final)



### SL Tajima's D ####
sl.tajimaD <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_r3/TajimaD/SL.tajimaD.out",sep="\t",header=F)
head(sl.tajimaD)
sl.tajiD <- sl.tajimaD[,c(1,7)]
names(sl.tajiD) <- c("contigs","TajimaD")
head(sl.tajiD)

G <- read.table("/home/ivan/Documents/Data/ddRAD/range_wide_old/G_map_new.csv",sep=",",header=T)
G<- G[,c(1:5)]
G.sl.tajiD.table <- merge(G,sl.tajiD, by.x="Contig",by.y="contigs")
G.sl.tajiD.table[ with(G.sl.tajiD.table,order(LG,Position)),]
G.sl.tajiD.table.2 <- G.sl.tajiD.table[ with(G.sl.tajiD.table,order(LG,Position)),]


G.sl.tajiD.table.final <- G.sl.tajiD.table.2[,c(1,4,5,6)]
names(G.sl.tajiD.table.final) <- c("contig","chr","bp","fst")
head(G.sl.tajiD.table.final)



#### plot ####
#save.image("~/Dropbox/Xiaodong PhD thesis/Chapter 3 Population Genomics/Analysis/plot_code/Fig3_new_W.RData")

### fst between with qtl as a reference for position ###
data=fst.bet.table
chr="chr"
pos="bp"
value="fst"

chr.index <- which(names(data)==chr)
pos.index <- which(names(data)==pos)
value.index <- which(names(data)==value)

dataset <- data[,c(chr.index,pos.index,value.index)]
names(dataset) <- c("chrom","opos","value")

gpos <- rep(NA, dim(dataset)[1])
gpos[1] <- 0
accumulate <- 0

global.lg.min <- c()
global.lg.max <- c()
for ( i in unique(dataset$chrom) ) {
  space <- 5
  ### deal with SNP dataset
  index <- which(dataset$chrom == i)
  lg.table <- dataset[index,]
  lg.min <- min(lg.table$opos,na.rm=T)
  global.lg.min <- append(global.lg.min,lg.min)
  lg.max <- max(lg.table$opos,na.rm = T)
  global.lg.max <- append(global.lg.max,lg.max)
  (lg.interval <- lg.max - lg.min)
  pos <- lg.table$opos - lg.min + accumulate
  gpos[index] <- pos
  ### deal with SNP dataset
  
  ### deal with QTL
  index2 <- which(table2$Chr ==i)
  lg.table.2 <- table2[index2,]
  table3[index2,c("Pos","Lower","Upper")] <- table2[index2,c("Pos","Lower","Upper")] - lg.min + accumulate
  ### deal with QTL
  
  accumulate <- accumulate + lg.interval + space
  
}

dataset$gpos <- gpos
dataset$group <- dataset$chrom %% 2

global.lg.min 
global.lg.max

( posmin <- tapply(dataset$gpos,dataset$chr, min) );
( posmax <- tapply(dataset$gpos,dataset$chr, max) );
posmean <- (posmin + posmax) /2 
pos.lg <- posmax - posmin


sp <- "ggplot(dataset, aes(x=gpos,y=value,colour=as.factor(group))) +geom_point(size=0.5,shape=20)" 
for ( i in unique(dataset$chrom) ) {
  j <- i %%2 + 3
  subset <- dataset[ which(dataset$chrom == i ),]
  na.index <- which(is.na(subset$value))
  if (length(na.index)>0) {
    subset <- subset[-na.index,]
  }
  mod.loess <- loess(subset$value ~ subset$gpos,span=.5)
  loess.predict <- predict(mod.loess)
  #lines(loess.predict, x=subset$gpos, col="red",lwd=2)  
  assign(paste("loess_predict_df",i,sep="_"),  data.frame(gpos=subset$gpos,value=loess.predict,valuelb=rep(0, length(loess.predict)),
                                                          valueub = loess.predict,group=rep(j,length(subset$gpos))) )
  
}
# plot(dataset$value, x=dataset$gpos, main="Loess Smoothing and Prediction",pch=20)
# lines(loess.predict, x=subset$gpos, col="red",lwd=2)


gl<-paste( "geom_line(size=.6,data=loess_predict_df_", 1:12, ")",sep="",collapse="+") # line width set here
gl2 <- paste( "geom_ribbon(aes(ymin=valuelb, ymax=valueub),data=loess_predict_df_", 1:12, ")",sep="",collapse="+")

panel <- paste(sp , gl,gl2,sep="+")
plot.final <- eval(parse(text=panel))  + 
  scale_color_manual(values=c("#E69F00", "#56B4E9","red","black")) + 
  theme(text=element_text(family="Times New Roman", face="bold", size=12)) 

###
plot1 <- plot.final + labs(y="fst bet.") +  theme(legend.position="none") +
  #geom_hline(yintercept = quantile(bet$WEIR_AND_COCKERHAM_FST,0.95,na.rm = T),linetype=3,col="black")+
  geom_vline(xintercept=posmean,linetype=2,col="gray",lwd=0.5)+
  
  theme(
    plot.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    
    axis.title.x=element_blank(),
    axis.text.x=element_blank(),
    axis.ticks.x=element_blank(),
    axis.text=element_text(size=12)
  ) +
  scale_y_continuous(limits = c(-0.05,1)) +
  xlim(0, sum(global.lg.max-global.lg.min)+space*11)

global.lg.max
global.lg.min
( global.x.label <- posmean )


#### new QTL plot
table4 <- table3[with(table3, order(Type,QTL)), ]

df <- data.frame()
raw.plot <- ggplot(df) + geom_point() + xlim(0, sum(global.lg.max-global.lg.min)+space*11) + ylim(3, 20)

cols=c("red","green","blue")
cols=c("red","green","gold","blue")


y_value <- c()
annotate_pos <- c()
y <- 0
old.type <- "new"
old.qtl <- "qtl"
for (i in 1:nrow(table4)) {
  qtl <- gsub(".[0-9]","",table4[i,"QTL"])
  print(qtl)
  if (table4[i,"Type"] == old.type && qtl == old.qtl) {
    y <- y
  } else if ( table4[i,"Type"] == old.type && qtl != old.qtl) {
    y <- y + 1
  } else if (table4[i,"Type"] != old.type) {
    y <- y + 3
    annotate_pos <- append(annotate_pos,y)
  }
  y_value <- append(y_value,y)
  old.type = table4[i,"Type"]
  old.qtl = qtl
}

table5<- cbind(table4,y_value)

qtl.plot <- raw.plot + geom_segment(data=table5,aes(x=Lower,xend=Upper,y=y_value,yend=y_value) ,col=cols[as.numeric(table5$Type)],size=1) +
             geom_point(data=table5,aes(x=Pos,y=y_value),size=1,col=cols[as.numeric(table5$Type)])  + 
             geom_vline(xintercept=posmean,linetype=2,col="gray",lwd=0.5)+ 
              geom_text(aes(y=annotate_pos,x=rep(0,length(annotate_pos))),label=substring(levels(table5$Type),1,3) )+
              theme(text=element_text(family="Times New Roman", face="bold", size=12),
                    plot.background = element_blank(),
                    panel.grid.major = element_blank(),
                    panel.grid.minor = element_blank(),
                    axis.title = element_blank(),
                    axis.ticks = element_blank(),
                    axis.text = element_blank() )
                    



ManhattanPlot6 <- function(data,contig="contig",chr="chr",pos="bp",value="fst",line=F,x.label=F,highlight=F,hlref=F,...) {
  if (highlight==T &&  hlref== F) {
    stop("Highlighting is enabled but no highlight reference is not given!")
  }
  contig.index <- which(names(data)==contig)
  chr.index <- which(names(data)==chr)
  pos.index <- which(names(data)==pos)
  value.index <- which(names(data)==value)
  
  dataset <- data[,c(contig.index,chr.index,pos.index,value.index)]
  names(dataset) <- c("contig","chrom","opos","value")
  
  gpos <- rep(NA, dim(dataset)[1])
  gpos[1] <- 0
  accumulate <- 0
  space <- 5
  for ( i in unique(dataset$chrom) ) {
    index <- which(dataset$chrom == i)
    lg.table <- dataset[index,]
    lg.min <- global.lg.min[i] 
    lg.max <- global.lg.max[i]
    lg.interval <- lg.max - lg.min
    pos <- lg.table$opos - lg.min + accumulate
    gpos[index] <- pos
    accumulate <- accumulate + lg.interval + space
  }
  dataset$gpos <- gpos
  dataset$group <- dataset$chrom %% 2
  
  posmin <- tapply(dataset$gpos,dataset$chr, min);
  posmax <- tapply(dataset$gpos,dataset$chr, max);
  posmean <- (posmin + posmax) /2
  pos.lg <- posmax - posmin
  
  if (highlight==T) {
    sp <- "ggplot(dataset, aes(x=gpos,y=value,colour=as.factor(group)))   +geom_point(...) +geom_point(data=dataset[which(dataset$contig %in% hlref),], aes(x=gpos,y=value),colour=\"red\",size=0.7)" 
  } else {
    sp <- "ggplot(dataset, aes(x=gpos,y=value,colour=as.factor(group)))  + geom_point(...)" 
  }
  ### loess curve ###
  for ( i in unique(dataset$chrom) ) {
    j <- i %%2 + 3
    subset <- dataset[ which(dataset$chrom == i ),]
    na.index <- which(is.na(subset$value))
    if (length(na.index) > 0) {
      subset <- subset[-na.index,] 
    }
    mod.loess <- loess(subset$value ~ subset$gpos,span=.5)
    loess.predict <- predict(mod.loess)
    #lines(loess.predict, x=subset$gpos, col="red",lwd=2)  
    assign(paste("loess_predict_df",i,sep="_"),  data.frame(gpos=subset$gpos,value=loess.predict,valuelb=rep(0, length(loess.predict)),
                                                            valueub = loess.predict,group=rep(j,length(subset$gpos))) )
    
  }
  # plot(dataset$value, x=dataset$gpos, main="Loess Smoothing and Prediction",pch=20)
  # lines(loess.predict, x=subset$gpos, col="red",lwd=2)
  
  
  gl<-paste( "geom_line(size=0.6,data=loess_predict_df_", 1:12, ")",sep="",collapse="+") # set line width here
  gl2 <- paste( "geom_ribbon(aes(ymin=valuelb, ymax=valueub),data=loess_predict_df_", 1:12, ")",sep="",collapse="+")
  ### end of loess curve ###
  
  ### window ###
  pos_win <- winScan(x = dataset, 
                     
                     position = "gpos", 
                     values = c("value"),
                     win_size = 3,
                     win_step = 3,
                     funs = c("function(x) mean(x,na.rm=T)")
  )
  names(pos_win)[ncol(pos_win)] <- "value"
  gw <- "geom_line(data=pos_win,aes(x=win_mid,y=value),colour=\"black\")"
  ### end of window ###
  
  if (line=="loess") {
    panel <- paste(sp , gl,gl2,sep="+")
  } else if (line=="window") {
    panel <- paste(sp,gw,sep="+")
  } else {
    panel <- sp
  }
  
  if (x.label==T) {	
    final.plot <- eval(parse(text=panel))  + 
      scale_color_manual( values=c("#E69F00", "#56B4E9","red","black") ) + 
      scale_x_continuous(breaks=c(global.x.label),labels=c("X/Y",1:11) ) +
      geom_vline(xintercept=posmean,linetype=2,col="gray",lwd=0.5)+
      xlab("Linkage groups") +
      theme(text=element_text(family="Times New Roman", face="bold", size=12),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            axis.text=element_text(size=12)
      ) 
  } else {
    final.plot <- eval(parse(text=panel))  + 
      scale_color_manual( values=c("#E69F00", "#56B4E9","red","black") ) + 
      geom_vline(xintercept=posmean,linetype=2,col="gray",lwd=0.5)+
      theme(text=element_text(family="Times New Roman", face="bold", size=12),
            plot.background = element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            
            axis.title.x=element_blank(),
            axis.text.x=element_blank(),
            axis.ticks.x=element_blank(),
            axis.text=element_text(size=12)
      ) 
  }
  
  return(final.plot)
  cat(pos.lg)
}

###




###

plot2 <- ManhattanPlot5(data=fst.wth.SD.table,chr="chr",pos="bp",value="fst",line="loess",x.label=F,size=.5,shape=20)+labs(y="fST S.d.") +  
  theme(legend.position="none")+ scale_y_continuous(limits = c(-0.05,1))

plot3 <- ManhattanPlot5(data=fst.wth.SL.table,chr="chr",pos="bp",value="fst",line="loess",x.label=F,size=.5,shape=20)+labs(y="fST S.l.") +  
  theme(legend.position="none")+ scale_y_continuous(limits = c(-0.05,0.6))


plot4 <- ManhattanPlot6(data=G.pi.sd.table.final,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+
  labs(y="pi for S.d.") +  theme(legend.position="none")+ scale_y_continuous(limits = c(0.00,0.016))

plot5 <- ManhattanPlot6(data=G.pi.sl.table.final,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+
  labs(y="pi for S.l.") +  theme(legend.position="none")+ scale_y_continuous(limits = c(0.00,0.016))

plot6<- ManhattanPlot6(data=G.dxy.table.final,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+
  labs(y="dXY bet.") +  theme(legend.position="none") + scale_y_continuous(limits = c(0,0.016)) 

plot7<-ManhattanPlot6(data=G.sd.tajiD.table.final,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+
  labs(y="Tajima's D sd") +  
  theme(legend.position="none") + scale_y_continuous(limits = c(-2,5))

plot8<-ManhattanPlot6(data=G.sl.tajiD.table.final,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+ 
  labs(y="Tajima's D  sl") +  
  theme(legend.position="none") + scale_y_continuous(limits = c(-2,5))

plot9<-ManhattanPlot5(data=gc.table,chr="chr",pos="bp",value="fst",line="loess",x.label=T,size=.5)+ labs(y="GC content") +  
  theme(legend.position="none") 

plot10<- ManhattanPlot6(data=hfstC.table,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+ 
  labs(y="contig based Fst") +  theme(legend.position="none") + scale_y_continuous(limits = c(0,1)) +
  geom_hline(yintercept=threshold,col="red",lty=2)

plot15<- ManhattanPlot5(data=hfstC_wth.table,chr="chr",pos="bp",value="fst",line="loess",x.label=F,size=.5)+ labs(y="contig based Fwth") +  
  theme(legend.position="none") + scale_y_continuous(limits = c(0,1))

plot12 <-ManhattanPlot6(data=hfstC.SD.table,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+
  labs(y="fST highest S.d.") +  theme(legend.position="none")+ scale_y_continuous(limits = c(-0.05,1))

plot13 <-ManhattanPlot6(data=hfstC.SL.table,contig="contig",chr="chr",pos="bp",value="fst",line="loess",highlight=T,hlref=island_contigs,x.label=F,size=0.7)+
  labs(y="fST highest S.l.") +  theme(legend.position="none")+ scale_y_continuous(limits = c(-0.05,1))


## main plot
g0 <- ggplotGrob(qtl.plot)
g10 <- ggplotGrob(plot10)
g6 <- ggplotGrob(plot6)
g4 <- ggplotGrob(plot4)
g5 <- ggplotGrob(plot5)
g12 <- ggplotGrob(plot12)
g13 <- ggplotGrob(plot13)

g <- gtable:::rbind_gtable(g0, g10, size="last")
g <- gtable:::rbind_gtable(g, g6, "first")
g <- gtable:::rbind_gtable(g, g4, "first")
g <- gtable:::rbind_gtable(g, g5, "first")
g <- gtable:::rbind_gtable(g, g12, "first")
g <- gtable:::rbind_gtable(g, g13, "first")
grid.newpage()
grid.draw(g)