### We calculated pairwise LD (r-square) using vcftools
#module load vcftools
#vcftools --gzvcf   input.vcf.gz    --geno-r2 --out   output.ld

## load required R packages
library(dplyr)
library(stringr)
library(ggplot2)
library(gridExtra)

## judaicum sub-population

data_jd <- read.delim("judaicum.LD.estimation.txt",sep="",header=F,check.names=F,stringsAsFactors=F)
colnames(data_jd) <- c("dist","rsq")
data_jd$distc <- cut(data_jd$dist,breaks=seq(from=min(data_jd$dist)-1,to=max(data_jd$dist)+1,by=100000))

dfr1 <- data_jd %>% group_by(distc) %>% summarise(mean=mean(rsq),median=median(rsq))

dfr1 <- dfr1 %>% mutate(start=as.integer(str_extract(str_replace_all(distc,"[\\(\\)\\[\\]]",""),"^[0-9-e+.]+")),
                        end=as.integer(str_extract(str_replace_all(distc,"[\\(\\)\\[\\]]",""),"[0-9-e+.]+$")),
                        mid=start+((end-start)/2))

P1 <- ggplot()+
  geom_point(data=dfr1,aes(x=start,y=mean),size=0.4,colour="grey20")+
  #geom_line(data=dfr1,aes(x=start,y=mean),size=0.3,alpha=0.5,colour="grey40")+
  labs(x="Distance (Megabases)",y=expression(LD~(r^{2})))+
  scale_x_continuous(labels=function(x)x/1000000, breaks = seq(0,900000000, by = 100000000),expand = c(0.01, 0)) +
  #scale_x_continuous(breaks=c(0,2*10^6,4*10^6,6*10^6,8*10^6),labels=c("0","2","4","6","8"))+
  theme_bw()

## Southern Levant sub-population

data_SL <- read.delim("Southern.Levant_Population.LD.estimation.txt",sep="",header=F,check.names=F,stringsAsFactors=F)
colnames(data_SL) <- c("dist","rsq")
data_SL$distc <- cut(data_SL$dist,breaks=seq(from=min(data_SL$dist)-1,to=max(data_SL$dist)+1,by=100000))

dfr2 <- data_SL %>% group_by(distc) %>% summarise(mean=mean(rsq),median=median(rsq))

dfr2 <- dfr2 %>% mutate(start=as.integer(str_extract(str_replace_all(distc,"[\\(\\)\\[\\]]",""),"^[0-9-e+.]+")),
                        end=as.integer(str_extract(str_replace_all(distc,"[\\(\\)\\[\\]]",""),"[0-9-e+.]+$")),
                        mid=start+((end-start)/2))

P2 <- ggplot()+
  geom_point(data=dfr2,aes(x=start,y=mean),size=0.4,colour="navy")+
  #geom_line(data=dfr1,aes(x=start,y=mean),size=0.3,alpha=0.5,colour="grey40")+
  labs(x="Distance (Megabases)",y=expression(LD~(r^{2})))+
  scale_x_continuous(labels=function(x)x/1000000, breaks = seq(0,900000000, by = 100000000),expand = c(0.01, 0)) +
  #scale_x_continuous(breaks=c(0,2*10^6,4*10^6,6*10^6,8*10^6),labels=c("0","2","4","6","8"))+
  theme_bw()


## Northern sub-population

data_NP <- read.delim("Northen_Population.LD.estimation.txt",sep="",header=F,check.names=F,stringsAsFactors=F)
colnames(data_NP) <- c("dist","rsq")
data_NP$distc <- cut(data_NP$dist,breaks=seq(from=min(data_NP$dist)-1,to=max(data_NP$dist)+1,by=100000))

dfr3 <- data_NP %>% group_by(distc) %>% summarise(mean=mean(rsq),median=median(rsq))

dfr3 <- dfr3 %>% mutate(start=as.integer(str_extract(str_replace_all(distc,"[\\(\\)\\[\\]]",""),"^[0-9-e+.]+")),
                        end=as.integer(str_extract(str_replace_all(distc,"[\\(\\)\\[\\]]",""),"[0-9-e+.]+$")),
                        mid=start+((end-start)/2))

## Combine all graph

P3 <- ggplot()+
  geom_point(data=dfr3,aes(x=start,y=mean),size=0.4,colour="red")+
  #geom_line(data=dfr1,aes(x=start,y=mean),size=0.3,alpha=0.5,colour="grey40")+
  labs(x="Distance (Megabases)",y=expression(LD~(r^{2})))+
  scale_x_continuous(labels=function(x)x/1000000, breaks = seq(0,900000000, by = 100000000),expand = c(0.01, 0)) +
  #scale_x_continuous(breaks=c(0,2*10^6,4*10^6,6*10^6,8*10^6),labels=c("0","2","4","6","8"))+
  theme_bw()

#P_all <- P1 + P2 +P3 

P_all <- grid.arrange(P1, P2, P3, ncol = 1)




#### merging plots together using another approach


library(ggplot2)

# Combine the three data frames into a single data frame
dfr1 <- bind_rows(
  mutate(dfr1, group = "judaicum"),
  mutate(dfr2, group = "Southern Levant"),
  mutate(dfr3, group = "Northern Population")
)

# Add a column for facetting based on the group
dfr1$facet_group <- factor(dfr1$group, levels = c("judaicum", "Southern Levant", "Northern Population"))

# Plot the data
P100 <- ggplot(dfr1, aes(x = start, y = mean, color = group)) +
  geom_point(size = 0.4) +
  labs(x = "Distance (Mb)", y = expression(LD ~ (r^{2}))) +
  scale_x_continuous(labels = function(x) x / 1000000, breaks = seq(0, 900000000, by = 100000000), expand = c(0.01, 0)) +
  theme_bw() +
  facet_wrap(~ facet_group, ncol = 1) +
  scale_color_manual(values = c("gray20", "red", "blue")) +
  labs(color = "Group") +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(face = "bold"),
    axis.text.y = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    strip.text = element_text(face = "bold")
  )

ggsave("Supplelmentary_FigS4_LD-decay.pub.pdf", P100, width = 11.5, height = 11.5, units = "in", limitsize = F) 
