---
title: "figure_1_samples_and_data"
author: "Jilong"
date: "3/3/2023"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(ggrepel)
library(ggimage)
library(treeio)
library(ggtree)
library(ape)
library(tidyverse)
library(grid)
library(ggpubr)
```

```{r}
pop <- c("S. lineatus","S. mimosarum","S. dumicola","S. dumicola","S. mimosarum","S. sarasinorum","S. tentoriicola","S. bicolor","S. sarasinorum","S. pacificus","S. africanus")
lat <- c(30.805311,-18.868499,-19.654697,-28.052735,-28.887923,27.007796,-33.458177,-15.93,7.458205,22.000555,-27.900512)
long <- c(34.825607,47.479210,17.327452,18.714993,30.042178,80,21.148091,25.336534,80.776628,79.603192,28.486073)
social <- c("subsocial","social","social","social","social","social","subsocial","subsocial","social","subsocial","subsocial")
data_type <- c("assembly","population","populaton","population","population","assembly","assembly","assembly","population","population","RNA")
sample_df <- data.frame(pop,long,lat,social,data_type)%>%
  mutate(sp = factor(pop,levels=c("S. pacificus","S. sarasinorum","S. bicolor","S. dumicola","S. tentoriicola","S. mimosarum","S. africanus","S. lineatus")))

world_map <- ne_countries(scale = "small", returnclass = "sf")

# Set the theme of our ggplot2 objects for this session
theme_set(theme_bw())

# The ggplot call is normal
DUM <- ggplot() + 
  geom_sf(data = world_map,fill=NA,color = "grey60") +
  geom_point(data = sample_df%>%filter(sp == "S. dumicola"),aes(y=lat,x=long))+
  geom_text(aes(x=1,y=-15,label = "Pairwise dN/dS:\n0.3366±0.0057\nSocial period:\n119kya(95%CI:72kya-170kya)"),hjust = 0)+
  scale_color_manual(values = c("#e34de0","#5b0c9c","#A1A7F5","#BB381C","#FAA390","#e3a60b","#34b1c9","#0D3C3E"))+
  #geom_label_repel(data = sample_df,aes(x=lat,y=long,label=rev(img)))+
  xlim(c(0, 30)) +
  ylim(c(-36, -11)) +
  xlab("longitude")+
  ylab("latitude")+
  ggtitle("S. dumicola")+
  theme(plot.title = element_text(face = "italic", family = "Times New Roman"),
        legend.title = element_text(), 
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())+
  guides(color=guide_legend(title="Species"))+
  NULL
DUM

MIM <- ggplot() + 
  geom_sf(data = world_map,fill=NA,color = "grey60") +
  geom_point(data = sample_df%>%filter(sp == "S. mimosarum"),aes(y=lat,x=long))+
  geom_text(aes(x=24,y=-15,label = "Pairwise dN/dS:\n0.3043±0.0028\nSocial period:\n279kya(95%CI:230kya-332kya)"),hjust = 0)+
  scale_color_manual(values = c("#e34de0","#5b0c9c","#A1A7F5","#BB381C","#FAA390","#e3a60b","#34b1c9","#0D3C3E"))+
  #geom_label_repel(data = sample_df,aes(x=lat,y=long,label=rev(img)))+
  xlim(c(23, 53)) +
  ylim(c(-36, -11)) +
  xlab("longitude")+
  ylab("latitude")+
  ggtitle("S.mimosarum")+
  theme(plot.title = element_text(face = "italic", family = "Times New Roman"),
        legend.title = element_text(), 
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())+
  guides(color=guide_legend(title="Species"))+
  NULL

MIM
SARA <- ggplot() + 
  geom_sf(data = world_map,fill=NA,color = "grey60") +
  geom_point(data = sample_df%>%filter(sp == "S. sarasinorum"),aes(y=lat,x=long))+
  geom_text(aes(x=61,y=26,label = "Pairwise dN/dS:\n0.2864±0.0035\nSocial period:\n260kya(95%CI:231kya-289kya)"),hjust = 0)+
  scale_color_manual(values = c("#e34de0","#5b0c9c","#A1A7F5","#BB381C","#FAA390","#e3a60b","#34b1c9","#0D3C3E"))+
  #geom_label_repel(data = sample_df,aes(x=lat,y=long,label=rev(img)))+
  xlim(c(60, 90)) +
  ylim(c(5, 30)) +
  xlab("longitude")+
  ylab("latitude")+
  ggtitle("S. sarasinorum")+
  theme(plot.title = element_text(face = "italic", family = "Times New Roman"),
        legend.title = element_text(), 
        legend.background = element_blank(),
        legend.key = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())+
  guides(color=guide_legend(title="Species"))+
  NULL

SARA
```

```{r}
library(ggpubr)
figure_1<- ggarrange(SARA,DUM,MIM,ncol = 3,labels = c("(A)","(B)","(C)"),align = "hv")
ggsave("population_pair.png",figure_1,width = 3000,height = 1000, units = "px")


```
```{r}
boostrap_autosomes <- read_tsv("auto_500_single_pac_clean.tab")
dnds_x <- read_tsv("x_100_single_pac_clean.tab")
auto_df <- boostrap_autosomes%>%
  mutate(branch = factor(nodename,levels = c("LIN","AFR","AFR_MIM","MIM","LIN_AFR_MIM_TENT_DUM_BI_PAC_SARA","TENT","TENT_DUM","DUM","TENT_DUM_BI_PAC_SARA","BI","BI_PAC_SARA","SARA","PAC_SARA","PAC")))%>%
  filter(branch!= "LIN_AFR_MIM_TENT_DUM_BI_PAC_SARA")%>%
  mutate(w = dN/dS)%>%
  group_by(branch)%>%
  summarise(
    w_median = mean(w),
    w_se = sd(w)/sqrt(500),
    w_q025 = w_median - 3*w_se,
    w_q975 = w_median + 3*w_se
  )%>%
  select(branch,w_median,w_q025,w_q975)%>%
  mutate(chrom = "Autosomes")

x_df <- dnds_x%>%
  mutate(branch = factor(nodename,levels = c("LIN","AFR","AFR_MIM","MIM","LIN_AFR_MIM_TENT_DUM_BI_PAC_SARA","TENT","TENT_DUM","DUM","TENT_DUM_BI_PAC_SARA","BI","BI_PAC_SARA","SARA","PAC_SARA","PAC")))%>%
  filter(branch!= "LIN_AFR_MIM_TENT_DUM_BI_PAC_SARA")%>%
  mutate(w = dN/dS)%>%
  group_by(branch)%>%
  summarise(
    w_median = mean(w),
    w_se = sd(w)/sqrt(100),
    w_q025 = w_median - 3*w_se,
    w_q975 = w_median + 3*w_se
  )%>%
  select(branch,w_median,w_q025,w_q975)%>%
  mutate(chrom = "X chromosomes")

dnds_df <- rbind(auto_df,x_df)

dnds_plot <- dnds_df %>%
  ggplot(aes(y=branch,color = branch),)+
  geom_point(aes(x=w_median,shape = chrom,size =1))+
  geom_segment(aes(x = w_q025,xend=w_q975,yend = branch))+
  geom_text(aes(x=w_median,label=round(w_median,3)),hjust = 0,vjust = -0.5)+
  xlab("dN/dS")+
  #scale_color_manual(values = c("#0D3C3E","#ffa200","#7a8e91","#e3a60b","#FAA390","#7a8e91","#BB381C","#7a8e91","#A1A7F5","#7a8e91","#5b0c9c","#7a8e91","#9b3fd1"))+
  scale_color_manual(values = c("#415CE5","#415CE5","#7a8e91","#D23333","#415CE5","#7a8e91","#D23333","#7a8e91","#415CE5","#7a8e91","#D23333","#7a8e91","#415CE5"))+
  guides(color=guide_legend(title="Chromosomes"))+
  scale_x_continuous(breaks = seq(0, 0.25, by = 0.02)) +
  theme(axis.text.y=element_blank(),
        axis.title.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        #legend.position = c(0.8,0.9),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        legend.position = "none",
        axis.line.x = element_line(color = "black", size = 0.5))
dnds_plot
  
```


```{r}
library(ggplot2)

# Vector of branch names
name_vec <- factor(c("LIN", "AFR", "AFR_MIM", "MIM","MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"), 
              levels  = c("LIN", "AFR", "AFR_MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"))
h_line_x <- c(0,4,1,4,0,4,2,4,1,3,2,4,3,4)
h_line_xend <- c(5,5,4,5,1,5,4,5,2,5,3,5,4,5)
h_line_y <- factor(c("LIN", "AFR", "AFR_MIM", "MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"),levels = c("LIN", "AFR", "AFR_MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"))
h_line_yend <- factor(c("LIN", "AFR", "AFR_MIM", "MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"),levels = c("LIN", "AFR", "AFR_MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"))
hline_df <- data.frame(
  x= h_line_x,
  xend = h_line_xend,
  y = h_line_y,
  yend = h_line_yend
)
v_line_x <- c(0,1,2,3,4,4,4)
v_line_xend <- c(0,1,2,3,4,4,4)
v_line_y <- c("LIN","AFR_MIM","TENT_DUM","PAC_SARA","PAC","TENT","AFR")
v_line_yend <- c("MIM","TENT_DUM_BI_PAC_SARA","BI_PAC_SARA","BI","SARA","DUM","MIM")
vline_df <- data.frame(
  x = v_line_x,
  xend = v_line_xend,
  y = factor(v_line_y,levels  = c("LIN", "AFR", "AFR_MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC")),
  yend = factor(v_line_yend, levels = c("LIN", "AFR", "AFR_MIM", "MIM", 
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"))
)

#text_vec <- c("S. lineatus", "S. africanus", "A_M", "S. mimosarum",
#              "S. tentoriicola", "T_D", "S. dumicola", "T_D_B_P_S", "S. bicolor", "B_P_S",
#              "S. sarasinorum", "P_S", "S. pacificus")
text_vec <- c("S. lineatus", "S. africanus", "", "S. mimosarum",
              "S. tentoriicola", "", "S. dumicola", "", "S. bicolor", "",
              "S. sarasinorum", "", "S. pacificus")
color_vec <- c("N", "N", "N", "S",
              "N", "N", "S", "N", "N", "N",
              "S", "N", "N")
x_vec <- c(rep(5.5,13))
y_vec <- factor(c("LIN", "AFR", "AFR_MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"),levels = c("LIN", "AFR", "AFR_MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"))
text_df <- data.frame(
  x= x_vec,
  y=y_vec,
  social = color_vec,
  text = text_vec
)

point_x <- c(5,5,4,5,5,4,5,2,5,3,5,4,5)
point_y <- factor(c("LIN", "AFR", "AFR_MIM", "MIM", 
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"),levels = c("LIN", "AFR", "AFR_MIM", "MIM",
              "TENT", "TENT_DUM", "DUM", "TENT_DUM_BI_PAC_SARA", "BI", "BI_PAC_SARA",
              "SARA", "PAC_SARA", "PAC"))
point_df <- data.frame(x = point_x,y=point_y)
test_df <- data.frame(branch = name_vec,x = h_line_x,xend= h_line_xend)
tree_test <- ggplot()+
  geom_segment(data = hline_df, aes(x=x,xend=xend,y=y,yend = yend),show.legend = FALSE)+
  geom_segment(data = vline_df,aes(x=x,xend=xend,y=y,yend = yend),show.legend = FALSE)+
  geom_point(data = point_df,aes(x=x,y=y))+
  geom_text(data = text_df,aes(x=x,y=y,label=text,color = social),hjust = 0,
            fontface = "italic",family = "Times New Roman",show.legend  = F)+
  xlim(0,7)+
  scale_color_manual(values = c("black","red"))+
  theme(axis.text.y=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.y=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_blank())


```
```{r}
legend_p <- dnds_df %>%
  mutate(CHROM = ifelse(chrom == "Autosomes",chrom,"X"))%>%
  ggplot(aes(y=branch))+
  geom_point(aes(x=w_median,shape = CHROM))
legend_p
legend <- get_legend(legend_p)
```

```{r}
dnds_branch<-ggarrange(tree_test,dnds_plot,ncol = 2,widths = c(1,1.8),labels = c("(D)",""),align = "hv")
dnds_branch
figure_2 <- ggarrange(figure_1,dnds_branch,ncol = 1,heights = c(1,1.2),align = "hv")
ggsave("branch_dnds.png",figure_2,width = 3600,height = 2400, units = "px")
```


