# correlation between genetic distance and geographic distance

rm(list=ls())

library(tidyverse)
library(ggalt)
library(ggforce)
library(marmap)

# Calculate least-cost distance crossing the marine environment using bathymetry data of British Isles
Seas <- readGEBCO.bathy("./British_Isles_bathy.nc")
# ensure we have a maximum depth cutoff to ensure distances don't go across land. 
trans <- trans.mat(Seas, min.depth=20, max.depth = -200)
#  coordinates of marine outline plots for each river catchment
lat_long <- read.csv("./outlet_coordinates.csv")
row.names(lat_long) <- lat_long$..Outlet
lat_long <- lat_long[,-1]
lat_long <- lat_long[,c(2,1)]

dist1 <- lc.dist(trans, lat_long, res= "dist")
dist.mat <- as.matrix(dist1)
dist.mat <- data.frame(dist.mat)
write.csv(dist.mat, "./marine_disatnces_between_outlets_km.csv")


# calculate isolation_by_distance
Fst_geo <- read.csv("./pairwise_FST_geographic_distance.csv")

# get log values
Fst_geo$log_FST <- log(Fst_geo$pairwise_FST+1)
Fst_geo$log_distance <- log(Fst_geo$Total_distance+1)

summary(lm(Fst_geo$log_FST~Fst_geo$log_distance))

fill <- c("#E69F00","#009E73","#56B4E9")
lb1 <- paste("R^2 == 0.009")
ggplot(data=Fst_geo, aes(log_distance, log_FST))+
  geom_point(shape = 21, size=2,aes(fill= Comparison_type))+
  geom_mark_hull(concavity= 5, expand=0, radius=0, aes(fill=Comparison_type), show.legend = F)+
  geom_smooth(method= 'lm',colour="black")+
  labs(x="log(Geographic distance+1) ", y="log(Pairwise FST+1)", colour="Comparison Type")+  
  scale_fill_manual(values=fill, labels=c('Different HA', 'Same HA', "Same Lake"))+
  theme(panel.background = element_rect(fill = 'white', colour= 'black'), 
        axis.title.x = element_text(size=16, face= "bold"),
        axis.text.x = element_text(size=12),
        axis.title.y = element_text(size=16, face= "bold"),
        axis.text.y = element_text(size=12),
        legend.text = element_text(size=12),
        legend.key=element_blank(),
        legend.title= element_blank(),
        legend.position = "none")+
  annotate("text", label= "P < 0.001", x=0.5, y=0.5, size =5)+
  annotate("text", label= lb1, x=0.5, y=0.48, size =5, parse= T)


# Run comparisons for within hydrometric areas only

sameHA <- filter(Fst_geo, Comparison_type == "Same Hydrometric Area")
summary(lm(sameHA$log_FST~sameHA$log_distance))

lb1 <- paste("R^2 == 0.17")
ggplot(data=sameHA, aes(log_distance, log_FST))+
  geom_point(shape = 21, size=2,aes(fill= Comparison_type))+
  geom_smooth(method= 'lm',colour="black")+
  scale_fill_manual(values="#009E73")+
  labs(x="log(Geographic distance+1) ", y="log(Pairwise FST+1)", colour="Comparison Type")+  
  theme(panel.border = element_rect(colour = "black", fill= "NA"),  
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.text.x = element_text(size=12),
        axis.title.y = element_text(size=16, face= "bold"),
        axis.title.x = element_text(size=16, face= "bold"),
        axis.text.y = element_text(size=12),
        legend.position = "none")+
  annotate("text", label= "P < 0.001", x=2.3, y=0.4, size =5)+
  annotate("text", label= lb1, x=2.3, y=0.38, size =5, parse= T)



# Run comparisons for across hydrometric areas only

diffHA <- filter(Fst_geo, Comparison_type == "Different Hydrometric Area")
summary(lm(diffHA$log_FST~diffHA$log_distance))

lb1 <- paste("R^2 == 0.003")
ggplot(data=diffHA, aes(log_distance, log_FST))+
  geom_point(shape = 21, size=2,aes(fill= Comparison_type))+
  geom_smooth(method= 'lm',colour="black")+
  scale_fill_manual(values="#E69F00")+
  labs(x="log(Geographic distance+1) ", y="log(Pairwise FST+1)", colour="Comparison Type")+  
  theme(panel.border = element_rect(colour = "black", fill= "NA"),  
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.text.x = element_text(size=12),
        axis.title.y = element_text(size=16, face= "bold"),
        axis.title.x = element_text(size=16, face= "bold"),
        axis.text.y = element_text(size=12),
        legend.position = "none")+
  annotate("text", label= "P = 0.011", x=3.9, y=0.5, size =5)+
  annotate("text", label= lb1, x=3.9, y=0.48, size =5, parse= T)
