Fig03.Rmd

---
title: "Data analysis"
author: "Aoyue Bi, Daxing Xu"
date: "`r Sys.Date()`"
output: 
  html_document: # html_notebook, html_document
    toc: true
    toc_float:
      toc_collapsed: true
    toc_depth: 3
    number_sections: true
    theme: simplex # default, united, cerulean, journal, flatly, darkly, readable, spacelab, united, cosmo, lumen, paper, sandstone, simplex, and yeti
    highlight: haddock # default, tango, pygments, kate, monochrome, espresso, zenburn, haddock, and breezedark
    # bibliography: bibliography.json
---

# Fig3

```{r setup, echo =F, message = T, warning = F}
library(tidyverse)
# library(RColorBrewer)
# library(ggpubr)
# library(egg)
## eval 代码要不要运行
## echo 代码要不要输出
## include 图要不要输出
## warning 警告要不要输出
## message 默认Bin等信息要不要输出
```

#### =================
## Del load ratio
```{r}
choiceA <- "Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5"

### 先处理del文件
dfdel <- read_tsv("data/004_delCount/delRatio.txt.gz") 

dftaxaDB <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>% 
  select(Taxa=VcfID,GenomeType,GroupbyContinent)

factorA <- c("AT","WE","DE","FTT","LR_EU","CL","LR_EA")
labelsA <- c("Ae. tauschii", "Wild emmer","Domesticated emmer","Free-threshing tetraploids","Landrace_European","Cultivar","Landrace_EastAsian")

factorB <- c("Ratio_002_nonsynonymous","Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5","Ratio_007_VEP")
labelsB <- c("Nonsynonymous","Deleterious","High impact")

df <- dfdel %>% left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(DelCountGroup == "TotalDerivedSNPCount") %>% 
  select(-c("IfSelected","Group_refobj","DelCountGroup")) %>% 
  filter(GroupbyContinent %in% factorA) %>% 
  filter(Stand_LoadGroup %in% factorB) %>% 
  droplevels() %>% 
  mutate(GroupbyContinent= factor(GroupbyContinent, levels = factorA,labels = labelsA)) %>% 
  mutate(Stand_LoadGroup=factor(Stand_LoadGroup,levels = factorB, labels = labelsB))

colB <- c("#87cef9","#ffd702","#7f5701","#016699", "#fc6e6e","#9900ff","gray80")

```

### plot
```{r}
p <- df %>%
  ggplot(aes(x=Subgenome,y=Stand_delCount,fill=GroupbyContinent))+
  ylab("Mutation burden")+xlab("")+
  ### volin
  geom_violin(position = position_dodge(0.8),alpha=0.8,width=2) +
  # geom_boxplot(position = position_dodge(0.8), width=0.01, color="white", alpha=0.2,outlier.color = NA) +
  geom_boxplot(position = position_dodge(0.8), width=0.2, color="white",
                lwd=0.2, ###调节box的线条粗细
               alpha=0.2,outlier.color = NA) +  # geom_point(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.8),aes(color=GroupbyContinent),alpha=0.3)+
  stat_summary(fun=mean, geom="point", color="red",position = position_dodge(0.8),size=0.15)+
  # stat_summary(fun=mean, geom="point", color="red",position = position_dodge(0.8),size=0.5)+

  ### boxplot
  # geom_boxplot(position = position_dodge(0.8),outlier.color = NA, alpha = 0.8) +
  # geom_point(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.8),aes(color=GroupbyContinent),alpha=0.3)+
  # stat_summary(fun=mean, geom="point", color="red",position = position_dodge(0.8),size=0.5)+
  
  ### pointrange
  # stat_summary(aes(color=GroupbyContinent),fun.data=mean_sdl,position = position_dodge(0.8),geom="pointrange",size=0.3)+
  
  scale_color_manual(values = colB,name='')+
  scale_fill_manual(values = colB,name='')+
  facet_grid(Stand_LoadGroup~.,scales = "free")+
  # facet_wrap(Stand_LoadGroup~.,scales = "free", ncol = 1)+
  # ggtitle("Del total count under all genes in vmap 2.1-2021") +
  theme(strip.background = element_blank())+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.4)
    # ,legend.position = 'none'
    ,text = element_text(size = 20))
p

ggsave(p,filename = "~/Documents/test.pdf",height = 3,width = 5)
# ggsave(p,filename = "figs/Fig3/misc/whole_load.pdf",height = 3,width = 4)
ggsave(p,filename = "~/Documents/test.png",height = 16,width = 18,units = "cm")

### get legend
library("ggpubr")
ggp_legend <- get_legend(p)
ggsave(as_ggplot(ggp_legend),filename = "~/Documents/legend.pdf",height = 3,width = 2)
```

### plot 分成3个图
```{r}
p <- df %>%
  ggplot(aes(x=Subgenome,y=Stand_delCount,fill=GroupbyContinent))+
  ylab("Mutation burden")+xlab("")+
  ### volin
  geom_violin(position = position_dodge(1),
              alpha=0.8,
              lwd=0.25, ###调节小提琴周围的线条粗细
              # violinwidth=10,
              width=2) +
  geom_boxplot(position = position_dodge(1), 
               width=0.4, color="white",
               lwd=0.2, ###调节box的线条粗细
               alpha=0.2,outlier.color = NA) +
  stat_summary(fun=mean, geom="point", color="red",position = position_dodge(1),size=0.1)+
  scale_color_manual(values = colB,name='')+
  scale_fill_manual(values = colB,name='')+
  # facet_grid(Stand_LoadGroup~.,scales = "free")+
  facet_wrap(~Stand_LoadGroup,scales = "free", ncol = 3)+
  theme(strip.background = element_blank())+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.4),
    # legend.position = 'none',
    legend.position = 'bottom',
    text = element_text(size = 18))
p

ggsave(p,filename = "~/Documents/test.pdf",height = 2,width = 7.2)
ggsave(p,filename = "~/Documents/test.png",height = 12,width = 32,units = "cm")

### get legend
library("ggpubr")
ggp_legend <- get_legend(p)
ggsave(as_ggplot(ggp_legend),filename = "~/Documents/legend.pdf",height = 3,width = 2)
```



### *！显著性检验 nonsynonymous
- 在分析显著性差异前，先将数据框赋予一个新的变量，存储下来。
```{r}
dfdf <- df
```

### 求六倍体亚基因组突变负荷的均值
```{r}
df_ave_mutationLoad <- dfdf %>% 
  group_by(Subgenome,GenomeType, Stand_LoadGroup) %>% 
  summarise(mean=mean(Stand_delCount))

df_ave_mutationLoad
```

#### 1.先计算平均值 标准差之类的
```{r}
### 平均数，标准差，方差等统计
  
df <- dfdf %>% 
  filter(Stand_LoadGroup == "Nonsynonymous") %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent)) %>% 
  group_by(Subgenome,GroupbyContinent) %>% 
summarise(n=n(),mean=mean(Stand_delCount),median=median(Stand_delCount),max=max(Stand_delCount),min=min(Stand_delCount),sd=sd(Stand_delCount), se = sd/sqrt(n)) %>% 
  ### 为数据框加上一个行号
  ungroup %>% 
  mutate(id=row_number()) %>% 
  relocate(id,.after = GroupbyContinent) %>% 
  relocate(mean,.before = n)


df
write.table(df,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 



```

#### 2.写一个小函数，专门计算升高降低百分比
```{r}
### 写一个小函数，专门计算升高降低百分比
aochangePercent <- function(before, after){
  out <-  (after-before)/before
  print(paste("before",before,"\n","after",after,"out",out))
    return(out)
}

## we - de
print("we -de")
aochangePercent(df[1,4],df[2,4])
aochangePercent(df[7,4],df[8,4])
print("******************************")
## de-ftt
print("de -ftt")
aochangePercent(df[2,4],df[3,4])
aochangePercent(df[8,4],df[9,4])
print("******************************")
## ftt-lr_eu
print("ftt - lr_eu")
aochangePercent(df[3,4],df[4,4])
aochangePercent(df[9,4],df[10,4])
print("******************************")
## ae-lr_eu cl lr_ea
print("ae - lr_eu cl lr_ea")
aochangePercent(df[13,4],df[14,4])
aochangePercent(df[13,4],df[15,4])
aochangePercent(df[13,4],df[16,4])
print("******************************")
## lr_eu - cl
print("lr_eu - cl")
aochangePercent(df[4,4],df[5,4])
aochangePercent(df[10,4],df[11,4])
aochangePercent(df[14,4],df[15,4])
## lr_eu - lr_ea
print("lr_eu - lr_ea")
aochangePercent(df[4,4],df[6,4])
aochangePercent(df[10,4],df[12,4])
aochangePercent(df[14,4],df[16,4])
print("******************************")

```
[1] "we -de"
[1] "before 0.823689315191681 \n after 0.828781520401876 out 0.00618219165439815"
[1] "before 0.747479099423756 \n after 0.751070594198122 out 0.00480480962897143"
[1] "******************************"
[1] "de -ftt"
[1] "before 0.828781520401876 \n after 0.825297132659683 out -0.00420422953024223"
[1] "before 0.751070594198122 \n after 0.75282539292127 out 0.00233639652078369"
[1] "******************************"
[1] "ftt - lr_eu"
[1] "before 0.825297132659683 \n after 0.829414645816221 out 0.00498912814984454"
[1] "before 0.75282539292127 \n after 0.745761747417329 out -0.00938284703247111"
[1] "******************************"
[1] "ae - lr_eu cl lr_ea"
[1] "before 0.790304233630642 \n after 0.811286037203311 out 0.0265490208451479"
[1] "before 0.790304233630642 \n after 0.811762753030456 out 0.0271522263030698"
[1] "before 0.790304233630642 \n after 0.811550511456706 out 0.02688366950593"
[1] "******************************"
[1] "lr_eu - cl"
[1] "before 0.829414645816221 \n after 0.829154673300759 out -0.000313440951126173"
[1] "before 0.745761747417329 \n after 0.749658880085458 out 0.00522570738124463"
[1] "before 0.811286037203311 \n after 0.811762753030456 out 0.00058760511741107"
[1] "lr_eu - lr_ea"
[1] "before 0.829414645816221 \n after 0.829480080951301 out 7.88931512234021e-05"
[1] "before 0.745761747417329 \n after 0.74440363031215 out -0.00182111392798463"
[1] "before 0.811286037203311 \n after 0.811550511456706 out 0.000325993843437309"
[1] "******************************"

#### 3.正态分布检验和同质性检验
```{r}
### 先检验 slightly load

df <- dfdf %>% 
  filter(Stand_LoadGroup == "Nonsynonymous") %>% 
  filter(Subgenome == "A") %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 检验 A亚基因组中各个亚群是否为正态分布
factor1 <- levels(df$GroupbyContinent)
for (i in factor1){
  p <- with(df, shapiro.test(Stand_delCount[GroupbyContinent == i]))
  print(paste(p[2],i))
}


#### 再检验各组数据是否具有同质性,这里只进行不同倍性的同质性检验
dfhomogeneity <- dfdf %>% 
  filter(Stand_LoadGroup == "Nonsynonymous") %>% 
  filter(Subgenome == "A") %>% 
  filter(GenomeType == "AABB" | GenomeType == "AABBDD")
res.ftest <- var.test(Stand_delCount ~ GenomeType, data = dfhomogeneity)
res.ftest

dfhomogeneity <- dfdf %>% 
  filter(Stand_LoadGroup == "Nonsynonymous") %>% 
  filter(Subgenome == "D") %>% 
  filter(GenomeType == "DD" | GenomeType == "AABBDD")
res.ftest <- var.test(Stand_delCount ~ GenomeType, data = dfhomogeneity)
res.ftest

### 结论： 有显著性差别，不来自同一整体

### 因为其个数都大于30，可以用t检验
```

#### 4.显著性检验
```{r}
### 先检验 slightly load
  ### ************** A sub
df <- dfdf %>% filter(Stand_LoadGroup == "Nonsynonymous") %>% 
  filter(Subgenome == "A") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p

  ### ************** B sub
  df <- dfdf %>% filter(Stand_LoadGroup == "Nonsynonymous") %>% 
  filter(Subgenome == "B") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p
  
    ### ************** D sub
  df <- dfdf %>% filter(Stand_LoadGroup == "Nonsynonymous") %>% 
  filter(Subgenome == "D") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p
  
### 将最终结果存储在下边，方便查看
```

- A sub 的两两比较 slightly

Pairwise comparisons using t tests with non-pooled SD 

data:  df$Stand_delCount and df$GroupbyContinent 

      WE      DE    FTT   LR_EU CL   
DE    0.001   -     -     -     -    
FTT   0.603   0.058 -     -     -    
LR_EU 8.0e-05 0.677 0.018 -     -    
CL    4.1e-05 0.765 0.018 0.769 -    
LR_EA 3.1e-05 0.640 0.014 0.930 0.651

P value adjustment method: BH

- ***************************************
- B sub 的两两比较  slightly

Pairwise comparisons using t tests with non-pooled SD 

data:  df$Stand_delCount and df$GroupbyContinent 

      WE      DE      FTT     LR_EU   CL     
DE    0.00443 -       -       -       -      
FTT   0.00027 0.11342 -       -       -      
LR_EU 0.12833 1.3e-09 3.4e-09 -       -      
CL    0.05741 0.05741 0.00351 1.1e-08 -      
LR_EA 0.00734 1.0e-14 7.1e-12 0.03294 < 2e-16

P value adjustment method: BH 

- ***************************************
- D sub 的两两比较  slightly

Pairwise comparisons using t tests with non-pooled SD 

data:  df$Stand_delCount and df$GroupbyContinent 

      AT      LR_EU CL  
LR_EU 1.5e-15 -     -   
CL    1.5e-15 0.48  -   
LR_EA 1.5e-15 0.61  0.61

P value adjustment method: BH 


### *！显著性检验 Deleterious
- 在分析显著性差异前，先将数据框赋予一个新的变量，存储下来。
```{r}
dfdf <- df
```
#### 1.先计算平均值 标准差之类的
```{r}
### 平均数，标准差，方差等统计
  
df <- dfdf %>% 
  filter(Stand_LoadGroup == "Deleterious") %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent)) %>% 
  group_by(Subgenome,GroupbyContinent) %>% 
summarise(n=n(),mean=mean(Stand_delCount),median=median(Stand_delCount),max=max(Stand_delCount),min=min(Stand_delCount),sd=sd(Stand_delCount), se = sd/sqrt(n)) %>% 
  ### 为数据框加上一个行号
  ungroup %>% 
  mutate(id=row_number()) %>% 
  relocate(id,.after = GroupbyContinent) %>% 
  relocate(mean,.before = n)


df
write.table(df,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 



```

### 求六倍体亚基因组突变负荷的均值
```{r}
df_ave_mutationLoad <- dfdf %>% 
  filter(Stand_LoadGroup=="Deleterious") %>% 
  group_by(GenomeType, Stand_LoadGroup) %>% 
  summarise(mean=mean(Stand_delCount))

df_ave_mutationLoad

### 写一个小函数，专门计算升高降低百分比
aochangePercent <- function(before, after){
  out <-  (after-before)/before
  print(paste("before",before,"\n","after",after,"out",out))
    return(out)
}

## we - de
print("AABB - DD")
aochangePercent(df_ave_mutationLoad[1,3],df_ave_mutationLoad[3,3])
print("******************************")
```


#### 2.写一个小函数，专门计算升高降低百分比
```{r}
### 写一个小函数，专门计算升高降低百分比
aochangePercent <- function(before, after){
  out <-  (after-before)/before
  print(paste("before",before,"\n","after",after,"out",out))
    return(out)
}

## we - de
print("we -de")
aochangePercent(df[1,4],df[2,4])
aochangePercent(df[7,4],df[8,4])
print("******************************")
## de-ftt
print("de -ftt")
aochangePercent(df[2,4],df[3,4])
aochangePercent(df[8,4],df[9,4])
print("******************************")
## ftt-lr_eu
print("ftt - lr_eu")
aochangePercent(df[3,4],df[4,4])
aochangePercent(df[9,4],df[10,4])
print("******************************")
## ae-lr_eu cl lr_ea
print("ae - lr_eu cl lr_ea")
aochangePercent(df[13,4],df[14,4])
aochangePercent(df[13,4],df[15,4])
aochangePercent(df[13,4],df[16,4])
print("******************************")
## lr_eu - cl
print("lr_eu - cl")
aochangePercent(df[4,4],df[5,4])
aochangePercent(df[10,4],df[11,4])
aochangePercent(df[14,4],df[15,4])
## lr_eu - lr_ea
print("lr_eu - lr_ea")
aochangePercent(df[4,4],df[6,4])
aochangePercent(df[10,4],df[12,4])
aochangePercent(df[14,4],df[16,4])
print("******************************")

```
[1] "we -de"
[1] "before 0.0908550905233067 \n after 0.094014408031045 out 0.0347731479825877"
[1] "before 0.0815522541853451 \n after 0.0807457269330525 out -0.00988969906901194"
[1] "******************************"
[1] "de -ftt"
[1] "before 0.094014408031045 \n after 0.0926389755499388 out -0.0146300179931153"
[1] "before 0.0807457269330525 \n after 0.0797135387778713 out -0.0127831923048637"
[1] "******************************"
[1] "ftt - lr_eu"
[1] "before 0.0926389755499388 \n after 0.0936783866530534 out 0.011220019402679"
[1] "before 0.0797135387778713 \n after 0.0782731025090795 out -0.0180701583554789"
[1] "******************************"
[1] "ae - lr_eu cl lr_ea"
[1] "before 0.0756236944627057 \n after 0.0818045555100843 out 0.0817318049758429"
[1] "before 0.0756236944627057 \n after 0.0821116875782117 out 0.0857931255752855"
[1] "before 0.0756236944627057 \n after 0.0818179481296709 out 0.0819089005234984"
[1] "******************************"
[1] "lr_eu - cl"
[1] "before 0.0936783866530534 \n after 0.0941582161720334 out 0.00512209417906781"
[1] "before 0.0782731025090795 \n after 0.079570539003144 out 0.0165757642469073"
[1] "before 0.0818045555100843 \n after 0.0821116875782117 out 0.00375446167040755"
[1] "lr_eu - lr_ea"
[1] "before 0.0936783866530534 \n after 0.09292348129595 out -0.00805847948576665"
[1] "before 0.0782731025090795 \n after 0.077133318004367 out -0.0145616369886495"
[1] "before 0.0818045555100843 \n after 0.0818179481296709 out 0.000163714838410854"
[1] "******************************"

#### 3.正态分布检验和同质性检验
```{r}
### 先检验 slightly load

df <- dfdf %>% 
  filter(Stand_LoadGroup == "Deleterious") %>% 
  filter(Subgenome == "A") %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 检验 A亚基因组中各个亚群是否为正态分布
factor1 <- levels(df$GroupbyContinent)
for (i in factor1){
  p <- with(df, shapiro.test(Stand_delCount[GroupbyContinent == i]))
  print(paste(p[2],i))
}


#### 再检验各组数据是否具有同质性,这里只进行不同倍性的同质性检验
dfhomogeneity <- dfdf %>% 
  filter(Stand_LoadGroup == "Deleterious") %>% 
  filter(Subgenome == "A") %>% 
  filter(GenomeType == "AABB" | GenomeType == "AABBDD")
res.ftest <- var.test(Stand_delCount ~ GenomeType, data = dfhomogeneity)
res.ftest

dfhomogeneity <- dfdf %>% 
  filter(Stand_LoadGroup == "Deleterious") %>% 
  filter(Subgenome == "D") %>% 
  filter(GenomeType == "DD" | GenomeType == "AABBDD")
res.ftest <- var.test(Stand_delCount ~ GenomeType, data = dfhomogeneity)
res.ftest

### 结论： 有显著性差别，不来自同一整体

### 因为其个数都大于30，可以用t检验
```

#### 4.显著性检验
```{r}
### 再检验 strongly load

df <- dfdf %>% filter(Stand_LoadGroup == "Deleterious") %>% 
  filter(Subgenome == "A") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p

  ### ************** B sub
  df <- dfdf %>% filter(Stand_LoadGroup == "Deleterious") %>% 
  filter(Subgenome == "B") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p
  
    ### ************** D sub
  df <- dfdf %>% filter(Stand_LoadGroup == "Deleterious") %>% 
  filter(Subgenome == "D") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p
  
### 将最终结果存储在下边，方便查看
```


### *！显著性检验 High impact
- 在分析显著性差异前，先将数据框赋予一个新的变量，存储下来。
```{r}
dfdf <- df
```
#### 1.先计算平均值 标准差之类的
```{r}
### 平均数，标准差，方差等统计
  
df <- dfdf %>% 
  filter(Stand_LoadGroup == "High impact") %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent)) %>% 
  group_by(Subgenome,GroupbyContinent) %>% 
summarise(n=n(),mean=mean(Stand_delCount),median=median(Stand_delCount),max=max(Stand_delCount),min=min(Stand_delCount),sd=sd(Stand_delCount), se = sd/sqrt(n)) %>% 
  ### 为数据框加上一个行号
  ungroup %>% 
  mutate(id=row_number()) %>% 
  relocate(id,.after = GroupbyContinent) %>% 
  relocate(mean,.before = n)


df
write.table(df,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 



```

### 求六倍体亚基因组突变负荷的均值
```{r}
df_ave_mutationLoad <- dfdf %>% 
  filter(Stand_LoadGroup=="High impact") %>% 
  group_by(GenomeType, Stand_LoadGroup) %>% 
  summarise(mean=mean(Stand_delCount))

df_ave_mutationLoad

### 写一个小函数，专门计算升高降低百分比
aochangePercent <- function(before, after){
  out <-  (after-before)/before
  print(paste("before",before,"\n","after",after,"out",out))
    return(out)
}

## we - de
print("AABB - DD")
aochangePercent(df_ave_mutationLoad[1,3],df_ave_mutationLoad[3,3])
print("******************************")
```


#### 2.写一个小函数，专门计算升高降低百分比
```{r}
### 写一个小函数，专门计算升高降低百分比
aochangePercent <- function(before, after){
  out <-  (after-before)/before
  print(paste("before",before,"\n","after",after,"out",out))
    return(out)
}

## we - de
print("we -de")
aochangePercent(df[1,4],df[2,4])
aochangePercent(df[7,4],df[8,4])
print("******************************")
## de-ftt
print("de -ftt")
aochangePercent(df[2,4],df[3,4])
aochangePercent(df[8,4],df[9,4])
print("******************************")
## ftt-lr_eu
print("ftt - lr_eu")
aochangePercent(df[3,4],df[4,4])
aochangePercent(df[9,4],df[10,4])
print("******************************")
## ae-lr_eu cl lr_ea
print("ae - lr_eu cl lr_ea")
aochangePercent(df[13,4],df[14,4])
aochangePercent(df[13,4],df[15,4])
aochangePercent(df[13,4],df[16,4])
print("******************************")
## lr_eu - cl
print("lr_eu - cl")
aochangePercent(df[4,4],df[5,4])
aochangePercent(df[10,4],df[11,4])
aochangePercent(df[14,4],df[15,4])
## lr_eu - lr_ea
print("lr_eu - lr_ea")
aochangePercent(df[4,4],df[6,4])
aochangePercent(df[10,4],df[12,4])
aochangePercent(df[14,4],df[16,4])
print("******************************")

```
[1] "we -de"
[1] "before 0.0156743807247154 \n after 0.0145735966540833 out -0.0702282335720177"
[1] "before 0.0121845611753524 \n after 0.0116592793531028 out -0.043110442361445"
[1] "******************************"
[1] "de -ftt"
[1] "before 0.0145735966540833 \n after 0.0146838286234365 out 0.00756381365352542"
[1] "before 0.0116592793531028 \n after 0.0117577844411589 out 0.00844864292833312"
[1] "******************************"
[1] "ftt - lr_eu"
[1] "before 0.0146838286234365 \n after 0.0139760167149071 out -0.0482034983301021"
[1] "before 0.0117577844411589 \n after 0.0106316958287079 out -0.0957738779858082"
[1] "******************************"
[1] "ae - lr_eu cl lr_ea"
[1] "before 0.0119451271339075 \n after 0.0113120986563776 out -0.052994704069153"
[1] "before 0.0119451271339075 \n after 0.0115741681830399 out -0.0310552534693894"
[1] "before 0.0119451271339075 \n after 0.0116515381701452 out -0.0245781363790545"
[1] "******************************"
[1] "lr_eu - cl"
[1] "before 0.0139760167149071 \n after 0.0139085673019989 out -0.00482608273044301"
[1] "before 0.0106316958287079 \n after 0.0109144700398256 out 0.0265972819081376"
[1] "before 0.0113120986563776 \n after 0.0115741681830399 out 0.023167188920732"
[1] "lr_eu - lr_ea"
[1] "before 0.0139760167149071 \n after 0.0130999981539277 out -0.0626801311739285"
[1] "before 0.0106316958287079 \n after 0.010548266337176 out -0.00784724214049229"
[1] "before 0.0113120986563776 \n after 0.0116515381701452 out 0.0300067674512494"
[1] "******************************"

#### 3.正态分布检验和同质性检验
```{r}
### 先检验 slightly load

df <- dfdf %>% 
  filter(Stand_LoadGroup == "High impact") %>% 
  filter(Subgenome == "A") %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 检验 A亚基因组中各个亚群是否为正态分布
factor1 <- levels(df$GroupbyContinent)
for (i in factor1){
  p <- with(df, shapiro.test(Stand_delCount[GroupbyContinent == i]))
  print(paste(p[2],i))
}


#### 再检验各组数据是否具有同质性,这里只进行不同倍性的同质性检验
dfhomogeneity <- dfdf %>% 
  filter(Stand_LoadGroup == "High impact") %>% 
  filter(Subgenome == "A") %>% 
  filter(GenomeType == "AABB" | GenomeType == "AABBDD")
res.ftest <- var.test(Stand_delCount ~ GenomeType, data = dfhomogeneity)
res.ftest

dfhomogeneity <- dfdf %>% 
  filter(Stand_LoadGroup == "High impact") %>% 
  filter(Subgenome == "D") %>% 
  filter(GenomeType == "DD" | GenomeType == "AABBDD")
res.ftest <- var.test(Stand_delCount ~ GenomeType, data = dfhomogeneity)
res.ftest

### 结论： 有显著性差别，不来自同一整体

### 因为其个数都大于30，可以用t检验
```
[1] "0.00293826133074987 WE"
[1] "0.237271564834943 DE"
[1] "0.0132178713200821 FTT"
[1] "0.934000460463813 LR_EU"
[1] "1.18935349900269e-05 CL"
[1] "0.000164430891960339 LR_EA"

	F test to compare two variances

data:  Stand_delCount by GenomeType
F = 2.0808, num df = 192, denom df = 407, p-value = 8.665e-10
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 1.640292 2.668161
sample estimates:
ratio of variances 
          2.080804 


	F test to compare two variances

data:  Stand_delCount by GenomeType
F = 0.35394, num df = 407, denom df = 35, p-value = 1.021e-06
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.2049301 0.5513441
sample estimates:
ratio of variances 
         0.3539367 

#### 4.显著性检验
```{r}
### 再检验 strongly load

df <- dfdf %>% filter(Stand_LoadGroup == "High impact") %>% 
  filter(Subgenome == "A") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p

  ### ************** B sub
  df <- dfdf %>% filter(Stand_LoadGroup == "High impact") %>% 
  filter(Subgenome == "B") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p
  
    ### ************** D sub
  df <- dfdf %>% filter(Stand_LoadGroup == "High impact") %>% 
  filter(Subgenome == "D") %>%  #### 这里可以切换亚基因组
  mutate(GroupbyContinent = factor(GroupbyContinent))

#### 两两比较，默认是 two.side
  p <- pairwise.t.test(df$Stand_delCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH")
  p
  
### 将最终结果存储在下边，方便查看
```

Pairwise comparisons using t tests with non-pooled SD 

data:  df$Stand_delCount and df$GroupbyContinent 

      WE      DE      FTT     LR_EU   CL     
DE    8.6e-14 -       -       -       -      
FTT   2.1e-09 0.45    -       -       -      
LR_EU < 2e-16 1.0e-06 1.3e-06 -       -      
CL    < 2e-16 2.3e-09 3.1e-08 0.44    -      
LR_EA < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16

P value adjustment method: BH 

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$Stand_delCount and df$GroupbyContinent 

      WE      DE      FTT     LR_EU   CL     
DE    1.1e-07 -       -       -       -      
FTT   0.00049 0.38974 -       -       -      
LR_EU < 2e-16 < 2e-16 4.6e-16 -       -      
CL    < 2e-16 3.4e-16 2.0e-11 3.3e-05 -      
LR_EA < 2e-16 < 2e-16 < 2e-16 0.23699 1.7e-08

P value adjustment method: BH 

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$Stand_delCount and df$GroupbyContinent 

      AT      LR_EU   CL   
LR_EU 1.0e-05 -       -    
CL    0.004   1.0e-05 -    
LR_EA 0.019   2.1e-07 0.070

P value adjustment method: BH 


## Del count
```{r}
### 先处理del文件
dfdel <- read_tsv("data/004_delCount/delCount.txt.gz") 

dftaxaDB <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>% 
  select(Taxa=VcfID,GenomeType,GroupbyContinent)

factorA <- c("AT","WE","DE","FTT","LR_EU","CL","LR_EA")
factorB <- c("002_nonsynonymous","031_GERP16way_1.5andDerived_PolyPhen2_0.5","007_VEP")
labelsB <- c("Nonsynonymous","Deleterious","High impact")

df <- dfdel %>% left_join(.,dftaxaDB,by="Taxa") %>% 
  select(-c("IfSelected","Group_refobj","SiteCountWithMinDepth")) %>% 
  filter(GroupbyContinent %in% factorA) %>% 
  filter(VariantsGroup %in% factorB) %>% 
  droplevels() %>% 
  mutate(GroupbyContinent= factor(GroupbyContinent, levels = factorA)) %>% 
  mutate(VariantsGroup=factor(VariantsGroup,levels = factorB, labels = labelsB))

colB <- c("#87cef9","#ffd702","#7f5701","#016699", "#fc6e6e","#9900ff","gray70")
```


### plot
```{r}
p <- df %>%
  ggplot(aes(x=Subgenome,y=TotalDerivedSNPCount,fill=GroupbyContinent))+
  labs(y="No. derived SNPs",x="")+
  ### volin
  # geom_violin(position = position_dodge(0.8),alpha=0.8,width=2) +
  # geom_boxplot(position = position_dodge(0.8), width=0.01, color="white", alpha=0.2,outlier.color = NA) +
  # # geom_point(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.8),aes(color=GroupbyContinent),alpha=0.3)+
  # stat_summary(fun=mean, geom="point", color="red",position = position_dodge(0.8),size=0.5)+

  ### boxplot
  # geom_boxplot(position = position_dodge(0.8),outlier.color = NA, alpha = 0.8) +
  # geom_point(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.8),aes(color=GroupbyContinent),alpha=0.3)+
  # stat_summary(fun=mean, geom="point", color="red",position = position_dodge(0.8),size=0.5)+
  
  ### pointrange
  stat_summary(aes(color=GroupbyContinent),fun.data=mean_sdl,
               fun.args = list(mult = 1),
               position = position_dodge(0.8),geom="pointrange")+
  stat_summary(fun=mean, geom="point", color="white",position = position_dodge(0.8))+
  
  scale_color_manual(values = colB,name='')+
  scale_fill_manual(values = colB,name='')+
  facet_grid(VariantsGroup~.,scales = "free")+
  # facet_wrap(Stand_LoadGroup~.,scales = "free", ncol = 1)+
  # ggtitle("Del total count under all genes in vmap 2.1-2021") +
  theme(strip.background = element_blank())+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.4)
    ,legend.position = 'none'
    ,text = element_text(size = 9))
p

ggsave(p,filename = "~/Documents/test.pdf",height = 3,width = 4)
# ggsave(p,filename = "figs/Fig3/misc/whole_load_count.pdf",height = 3,width = 4)
```










```{r}
### 写一个小函数，专门计算升高降低百分比
aochangePercent <- function(before, after){
  out <-  (after-before)/before
  print(paste("before",before,"\n","after",after,"out",out))
    return(out)
}

## we - de
print("we -de")
aochangePercent(df[1,4],df[2,4])
aochangePercent(df[7,4],df[8,4])
print("******************************")
## de-ftt
print("de -ftt")
aochangePercent(df[2,4],df[3,4])
aochangePercent(df[8,4],df[9,4])
print("******************************")
## ftt-lr 
print("ftt - lr")
aochangePercent(df[3,4],df[4,4])
aochangePercent(df[9,4],df[10,4])
print("******************************")
## ae-lr cl
print("ae - lr cl")
aochangePercent(df[13,4],df[14,4])
aochangePercent(df[13,4],df[15,4])
print("******************************")
## lr - cl
print("lr - cl")
aochangePercent(df[4,4],df[5,4])
aochangePercent(df[10,4],df[11,4])
aochangePercent(df[14,4],df[15,4])
print("******************************")

```




#### =============
## XPCLR load ratio
```{r}
dfdelRatio <- read.delim("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/004_summary_XPCLR_top005/012_merge/delRatio.txt")

### 导入种质数据库,并选取特定的列
dftaxaDB <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>%
  select(Taxa=VcfID,GenomeType,GroupbyContinent,IfOnVmap2_S643)

df <- dfdelRatio %>% left_join(.,dftaxaDB,by="Taxa")
df$GroupbyContinent <- factor(df$GroupbyContinent, levels = c("WE","DE","FTT","LR_EU","CL"))
colB <- c("#ffd702","#7f5701","#016699", "#fc6e6e","#9900ff")

df$Stand_LoadGroup <- factor(df$Stand_LoadGroup,levels = c("Ratio_002_nonsynonymous","Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5","Ratio_007_VEP"), labels = c("Nonsynonymous","Deleterious","High impact"))

df$Group_refobj <- factor(df$Group_refobj, levels = c("wede","dedurum","lrcul"))

p <- df %>% 
  filter(IfSelected == 1) %>% 
  filter(DelCountGroup == "TotalDerivedSNPCount") %>% 
  ggplot(aes(x=Subgenome,y=Stand_delCount))+
  ylab("Mutation burden")+xlab("")+
  ### pointrange
  stat_summary(aes(color=GroupbyContinent),fun.data=mean_sdl,
               fun.args = list(mult = 1),
               position=position_dodge(0.8),geom="pointrange")+
  stat_summary(aes(group = GroupbyContinent), fun=mean, geom="point", color="white",position = position_dodge(0.8))+
  scale_color_manual(values = colB,name='')+
  scale_fill_manual(values = colB,name='')+
  facet_grid(Stand_LoadGroup~Group_refobj,scales = "free")+
  theme(strip.background = element_blank())+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7)
    ,legend.position = 'none'
    ,text = element_text(size = 10))
p

# ggsave(p,filename = "~/Documents/test.pdf",height = 3,width = 4)
ggsave(p,filename = "~/Documents/test.pdf",height =3,width = 3)

```

## XPCLR load count
```{r}
dfdelcount <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/004_summary_XPCLR_top005/012_merge/delCount_xpclr.txt")

### 导入种质数据库,并选取特定的列
dftaxaDB <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>%
  select(Taxa=VcfID,GenomeType,GroupbyContinent)

factorA <- c("002_nonsynonymous","031_GERP16way_1.5andDerived_PolyPhen2_0.5","007_VEP")
labelsA <- c("Nonsynonymous","Deleterious","High impact")

labelsB <- c("Wild emmer","Domesticated emmer", "Durum","Landrace","Cultivar")

df <- dfdelcount %>% left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(VariantsGroup %in% factorA)

df$GroupbyContinent <- factor(df$GroupbyContinent, levels = c("WE","DE","FTT","LR_EU","CL"),
                              labels = labelsB)
colB <- c("#ffd702","#7f5701","#016699", "#fc6e6e","#9900ff")
df$VariantsGroup <- factor(df$VariantsGroup,levels = factorA, labels = labelsA)
df$Group_refobj <- factor(df$Group_refobj, levels = c("wede","dedurum","lrcul"))

p <- df %>% 
  filter(IfSelected == 1) %>% 
  ggplot(aes(x=Subgenome,y=TotalDerivedSNPCount))+
  ylab("Derived allele count per accession")+xlab("")+
  ### pointrange
  stat_summary(aes(color=GroupbyContinent),fun.data=mean_sdl,
               fun.args = list(mult = 1),size=0.4,
               position=position_dodge(0.95),geom="pointrange")+
  stat_summary(aes(group = GroupbyContinent), fun=mean, geom="point", color="white",
               size=0.9,
               position = position_dodge(0.95))+
  scale_color_manual(values = colB,name='')+
  scale_fill_manual(values = colB,name='')+
  facet_grid(VariantsGroup~Group_refobj,scales = "free")+
  theme(strip.background = element_blank(),
        strip.text = element_text(size=6))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.5),
    # legend.position = 'none',
    legend.position = 'bottom',legend.title = element_blank(),legend.key = element_blank(),
    text = element_text(size = 9))
p

# ggsave(p,filename = "~/Documents/test.pdf",height = 3,width = 4)
# ggsave(p,filename = "~/Documents/test.pdf",height = 2.4,width = 2.4)
# ggsave(p,filename = "~/Documents/test.pdf",height = 3,width = 3)
ggsave(p,filename = "~/Documents/test.pdf",height = 2.8,width = 2.8)

library("ggpubr")
ggp_legend <- get_legend(p)
ggsave(as_ggplot(ggp_legend),filename = "~/Documents/legend.pdf",height = 3,width = 7.2)

```

### ！显著性检验

- 在分析显著性差异前，先将数据框赋予一个新的变量，存储下来。
```{r}
dfdf <- df
```

#### 1.先计算平均值 标准差之类的

```{r}
### 平均数，标准差，方差等统计
  df <- dfdf %>% 
  select(Taxa,Subgenome,VariantsGroup,Group_refobj,
         TotalDerivedSNPCount,GenomeType,GroupbyContinent) %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent)) %>% 
  group_by(VariantsGroup,Group_refobj,Subgenome,GroupbyContinent) %>% 
summarise(n=n(),mean=mean(TotalDerivedSNPCount),median=median(TotalDerivedSNPCount),max=max(TotalDerivedSNPCount),min=min(TotalDerivedSNPCount),sd=sd(TotalDerivedSNPCount), se = sd/sqrt(n)) %>% 
  ### 为数据框加上一个行号
  ungroup %>% 
  mutate(id=row_number()) %>% 
  relocate(id,.after = GroupbyContinent) %>% 
  relocate(mean,.before = n)


df
write.table(df,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 

#### 注意：该数据框用于下文计算升高降低比例使用，不要随意改动！！！

```

#### 2.写一个小函数，专门计算升高降低百分比

```{r}
aochangePercent <- function(before, after){
  out <-  (after-before)/before
  print(paste("before",before,"\n","after",after,"out",out))
    return(out)
}

## Nonsynonymous
print("************* Nonsynonymous *********start")
### we-de
print("we -de")
aochangePercent(df[1,6],df[2,6])
aochangePercent(df[3,6],df[4,6])

### de-durum
print("de-durum")
aochangePercent(df[5,6],df[6,6])
aochangePercent(df[7,6],df[8,6])

### lr-cul
print("lr-cul")
aochangePercent(df[9,6],df[10,6])
aochangePercent(df[11,6],df[12,6])
aochangePercent(df[13,6],df[14,6])
print("************* Nonsynonymous *********end")



## Deleterious
print("************* Deleterious *********start")
### we-de
print("we -de")
aochangePercent(df[15,6],df[16,6])
aochangePercent(df[17,6],df[18,6])

### de-durum
print("de-durum")
aochangePercent(df[19,6],df[20,6])
aochangePercent(df[21,6],df[22,6])

### lr-cul
print("lr-cul")
aochangePercent(df[23,6],df[24,6])
aochangePercent(df[25,6],df[26,6])
aochangePercent(df[27,6],df[28,6])
print("************* Deleterious *********end")


## High impact
print("************* High impact *********start")
### we-de
print("we -de")
aochangePercent(df[29,6],df[30,6])
aochangePercent(df[31,6],df[32,6])

### de-durum
print("de-durum")
aochangePercent(df[33,6],df[34,6])
aochangePercent(df[35,6],df[36,6])

### lr-cul
print("lr-cul")
aochangePercent(df[37,6],df[38,6])
aochangePercent(df[39,6],df[40,6])
aochangePercent(df[41,6],df[42,6])
print("************* High impact *********end")
```

[1] "************* Nonsynonymous *********start"
[1] "we -de"
[1] "before 331.445054945055 \n after 320.153846153846 out -0.0340666080930989"
[1] "before 236.978021978022 \n after 230.951923076923 out -0.0254289357755623"
[1] "de-durum"
[1] "before 347.576923076923 \n after 319.527777777778 out -0.0806991012258248"
[1] "before 259.788461538462 \n after 238.5 out -0.0819453697534977"
[1] "lr-cul"
[1] "before 294.646551724138 \n after 290.103603603604 out -0.0154182972580232"
[1] "before 244.362068965517 \n after 247.873873873874 out 0.0143713176239811"
[1] "before 274.931034482759 \n after 268.824324324324 out -0.0222117891125792"
[1] "************* Nonsynonymous *********end"
[1] "************* Deleterious *********start"
[1] "we -de"
[1] "before 44.4725274725275 \n after 42.5673076923077 out -0.0428403755868545"
[1] "before 35.8516483516484 \n after 40.3173076923077 out 0.12455938697318"
[1] "de-durum"
[1] "before 46.3076923076923 \n after 43.8055555555556 out -0.054032853451458"
[1] "before 36.625 \n after 31.4166666666667 out -0.142207053469852"
[1] "lr-cul"
[1] "before 45.8965517241379 \n after 44.9864864864865 out -0.0198286190021727"
[1] "before 34.9051724137931 \n after 36.3738738738739 out 0.042076900313502"
[1] "before 38.198275862069 \n after 36.981981981982 out -0.0318415910832973"
[1] "************* Deleterious *********end"
[1] "************* High impact *********start"
[1] "we -de"
[1] "before 5.03846153846154 \n after 4.64423076923077 out -0.0782442748091603"
[1] "before 3.22527472527473 \n after 2.71153846153846 out -0.159284497444634"
[1] "de-durum"
[1] "before 1.68269230769231 \n after 1.66666666666667 out -0.00952380952380949"
[1] "before 3.31730769230769 \n after 2.86111111111111 out -0.137520128824477"
[1] "lr-cul"
[1] "before 2.54310344827586 \n after 3.4009009009009 out 0.337303405100015"
[1] "before 2.08620689655172 \n after 2.63513513513514 out 0.26312262675899"
[1] "before 6.56896551724138 \n after 6.26126126126126 out -0.0468421177082593"
[1] "************* High impact *********end"

#### 3.显著性检验

```{r}
# p <- dfdf %>% 
#   filter(Subgenome != "D") %>% 
#   mutate(GroupbyContinent = factor(GroupbyContinent)) %>% 
#   group_by(VariantsGroup,Group_refobj,Subgenome) %>% 
#   group_map(~{pairwise.t.test(.x$TotalDerivedSNPCount,.x$GroupbyContinent
#                                  ,paired = F,pool.sd = FALSE
#                                  ,var.equal=F,p.adjust.method = "BH")
#     
#     })
         
### 第一步：由于自己技术有限，不能在数据框分组中套用显著性检验函数，故一一进行检验
### 需要进行21次显著性检验

### ***********************第二步：单次测试 成功
# df <- dfdf %>% 
#   filter(VariantsGroup=="Nonsynonymous") %>% 
#   filter(Group_refobj=="wede") %>% 
#   filter(Subgenome=="A")
# 
# #### 两两比较，默认是 two.side
# p <- pairwise.t.test(df$TotalDerivedSNPCount,df$GroupbyContinent
#                        ,paired = F,pool.sd = FALSE
#                        ,var.equal = F
#                        ,p.adjust.method = "BH")
# p
 
### ***********************第三步：函数测试 成功
fun_getSig <- function(df){
  p <- pairwise.t.test(df$TotalDerivedSNPCount,df$GroupbyContinent
                       ,paired = F,pool.sd = FALSE
                       ,var.equal = F
                       ,p.adjust.method = "BH") 
  p

}

# df <- dfdf %>% 
#   filter(VariantsGroup=="Nonsynonymous") %>% 
#   filter(Group_refobj=="wede") %>% 
#   filter(Subgenome=="A") %>% 
#   fun_getSig
  
#################### wede
### ***********************函数调用，正式运行
### 1
print("Nonsynonymous ---- wede ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="Nonsynonymous") %>% 
  filter(Group_refobj=="wede") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 2
print("Nonsynonymous ---- wede ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="Nonsynonymous") %>% 
  filter(Group_refobj=="wede") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

### 3
print("Deleterious ---- wede ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="Deleterious") %>% 
  filter(Group_refobj=="wede") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 4
print("Deleterious ---- wede ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="Deleterious") %>% 
  filter(Group_refobj=="wede") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

### 5
print("High impact ---- wede ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="High impact") %>% 
  filter(Group_refobj=="wede") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 6
print("High impact ---- wede ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="High impact") %>% 
  filter(Group_refobj=="wede") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

#################### dedurum
  ### ***********************函数调用，正式运行
### 7
print("Nonsynonymous ---- dedurum ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="Nonsynonymous") %>% 
  filter(Group_refobj=="dedurum") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 8
print("Nonsynonymous ---- dedurum ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="Nonsynonymous") %>% 
  filter(Group_refobj=="dedurum") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

### 9
print("Deleterious ---- dedurum ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="Deleterious") %>% 
  filter(Group_refobj=="dedurum") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 10
print("Deleterious ---- dedurum ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="Deleterious") %>% 
  filter(Group_refobj=="dedurum") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

### 11
print("High impact ---- dedurum ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="High impact") %>% 
  filter(Group_refobj=="dedurum") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 12
print("High impact ---- dedurum ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="High impact") %>% 
  filter(Group_refobj=="dedurum") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

#################### lrcul
  ### ***********************函数调用，正式运行
### 13
print("Nonsynonymous ---- lrcul ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="Nonsynonymous") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 14
print("Nonsynonymous ---- lrcul ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="Nonsynonymous") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

### 15
print("Nonsynonymous ---- lrcul ---- D")
df <- dfdf %>% 
  filter(VariantsGroup=="Nonsynonymous") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="D") %>% 
  fun_getSig
df

### 16
print("Deleterious ---- lrcul ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="Deleterious") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 17
print("Deleterious ---- lrcul ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="Deleterious") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

### 18
print("Deleterious ---- lrcul ---- D")
df <- dfdf %>% 
  filter(VariantsGroup=="Deleterious") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="D") %>% 
  fun_getSig
df


### 19
print("High impact ---- lrcul ---- A")
df <- dfdf %>% 
  filter(VariantsGroup=="High impact") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="A") %>% 
  fun_getSig
df

### 20
print("High impact ---- lrcul ---- B")
df <- dfdf %>% 
  filter(VariantsGroup=="High impact") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="B") %>% 
  fun_getSig
df

### 21
print("High impact ---- lrcul ---- D")
df <- dfdf %>% 
  filter(VariantsGroup=="High impact") %>% 
  filter(Group_refobj=="lrcul") %>% 
  filter(Subgenome=="D") %>% 
  fun_getSig
df

```



[1] "Nonsynonymous ---- wede ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   WE   
DE 0.041

P value adjustment method: BH 
[1] "Nonsynonymous ---- wede ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   WE  
DE 0.12

P value adjustment method: BH 
[1] "Deleterious ---- wede ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   WE   
DE 0.026

P value adjustment method: BH 
[1] "Deleterious ---- wede ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   WE     
DE 3.5e-07

P value adjustment method: BH 

[1] "High impact ---- wede ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   WE   
DE 0.045

P value adjustment method: BH 
[1] "High impact ---- wede ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   WE    
DE 0.0015

P value adjustment method: BH 
[1] "Nonsynonymous ---- dedurum ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

    DE   
FTT 0.021

P value adjustment method: BH 
[1] "Nonsynonymous ---- dedurum ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

    DE   
FTT 0.019

P value adjustment method: BH 

[1] "Deleterious ---- dedurum ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

    DE  
FTT 0.17

P value adjustment method: BH 
[1] "Deleterious ---- dedurum ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

    DE     
FTT 0.00032

P value adjustment method: BH 
[1] "High impact ---- dedurum ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

    DE  
FTT 0.95

P value adjustment method: BH 
[1] "High impact ---- dedurum ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

    DE   
FTT 0.028

P value adjustment method: BH 

[1] "Nonsynonymous ---- lrcul ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU
CL 0.04 

P value adjustment method: BH 
[1] "Nonsynonymous ---- lrcul ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU
CL 0.083

P value adjustment method: BH 
[1] "Nonsynonymous ---- lrcul ---- D"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU  
CL 0.00086

P value adjustment method: BH 
[1] "Deleterious ---- lrcul ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU
CL 0.094

P value adjustment method: BH 

[1] "Deleterious ---- lrcul ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU
CL 0.014

P value adjustment method: BH 
[1] "Deleterious ---- lrcul ---- D"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU  
CL 0.00052

P value adjustment method: BH 
[1] "High impact ---- lrcul ---- A"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU
CL 1e-13

P value adjustment method: BH 
[1] "High impact ---- lrcul ---- B"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU  
CL 0.00023

P value adjustment method: BH 

[1] "High impact ---- lrcul ---- D"

	Pairwise comparisons using t tests with non-pooled SD 

data:  df$TotalDerivedSNPCount and df$GroupbyContinent 

   LR_EU
CL 0.028

P value adjustment method: BH 

#### 4.其他类型的平均值统计
```{r}
### 平均数，标准差，方差等统计
  df <- dfdf %>% 
  select(Taxa,Subgenome,VariantsGroup,Group_refobj,
         TotalDerivedSNPCount,GenomeType,GroupbyContinent) %>% 
  mutate(GroupbyContinent = factor(GroupbyContinent)) %>% 
  group_by(VariantsGroup,Group_refobj,GroupbyContinent) %>% 
summarise(n=n(),mean=mean(TotalDerivedSNPCount),median=median(TotalDerivedSNPCount),max=max(TotalDerivedSNPCount),min=min(TotalDerivedSNPCount),sd=sd(TotalDerivedSNPCount), se = sd/sqrt(n)) %>% 
  ### 为数据框加上一个行号
  ungroup %>% 
  mutate(id=row_number()) %>% 
  relocate(id,.after = GroupbyContinent) %>% 
  relocate(mean,.before = n)


df
write.table(df,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 
```


#### ==============
## DelLoad only Hexaploid
```{r}
choiceA <- "Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5"

### 先处理del文件
dfdel <- read_tsv("data/004_delCount/delRatio.txt.gz") 

dftaxaDB <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>% 
  select(Taxa=VcfID,GenomeType,GroupbyContinent)

factorA <- c("LR_EU","CL","LR_EA")

# factorB <- c("Ratio_002_nonsynonymous","Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5","Ratio_007_VEP")
# labelsB <- c("Nonsynonymous","Deleterious","High impact")

factorB <- c("Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5")
labelsB <- c("Deleterious")

df <- dfdel %>% left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(DelCountGroup == "TotalDerivedSNPCount") %>% 
  select(-c("IfSelected","Group_refobj","DelCountGroup")) %>% 
  filter(GroupbyContinent %in% factorA) %>% 
  filter(Stand_LoadGroup %in% factorB) %>% 
  droplevels() %>% 
  mutate(GroupbyContinent= factor(GroupbyContinent, levels = factorA)) %>% 
  mutate(Stand_LoadGroup=factor(Stand_LoadGroup,levels = factorB, labels = labelsB))

colB <- c( "#fc6e6e","#9900ff","gray90")

```

### plot
```{r}
p <- df %>%
  ggplot(aes(x=Subgenome,y=Stand_delCount,fill=GroupbyContinent))+
  ylab("Mutation burden")+xlab("")+
  ### volin
  # geom_violin(position = position_dodge(0.8),alpha=0.8,width=2) +
  # geom_boxplot(position = position_dodge(0.8), width=0.02, color="white", alpha=0.2) +
  # geom_point(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.8),aes(color=GroupbyContinent),alpha=0.3)+
  # stat_summary(fun=mean, geom="point", color="red",position = position_dodge(0.8),size=0.5)+

  ### boxplot
  geom_boxplot(position = position_dodge(0.8),outlier.color = NA, alpha = 0.8) +
  geom_point(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.8),aes(color=GroupbyContinent),alpha=0.3)+
  stat_summary(fun=mean, geom="point", color="red",position = position_dodge(0.8),size=0.5)+
  
  ### pointrange
  # stat_summary(aes(color=GroupbyContinent),fun.data=mean_sdl,position = position_dodge(0.8),geom="pointrange",size=0.3)+
  
  scale_color_manual(values = colB,name='')+
  scale_fill_manual(values = colB,name='')+
  facet_wrap(Stand_LoadGroup~.,scales = "free", ncol = 1)+
  # ggtitle("Del total count under all genes in vmap 2.1-2021") +
  theme(strip.background = element_blank())+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7)
    ,legend.position = 'none'
    ,text = element_text(size = 12))
p

# ggsave(p,filename = "~/Documents/test.pdf",height = 3,width = 4)
ggsave(p,filename = "figs/Fig3/misc/whole_load.pdf",height = 3,width = 4)
```


#### ===============
## Introgression
#### ==============
## Data important
```{r}
# table <- read_tsv("data/fromDaxing/fd/006_introgression_burden/001_popFdIntervalBurden/introgressionDonorBurden.txt.gz") %>% 
#   mutate(Sub=substring(Chr, 2,2)) %>% 
#   mutate(P2=factor(P2,levels = c("LR_AM","LR_AF","LR_CSA","LR_WA","LR_EU","CL","LR_EA"))) %>% 
#   mutate(P3=factor(P3,levels = c("WE","DE","FTT","AT"),labels = c("WE","DE","FTT","AT"))) 
# 
# df <- table %>% 
#   group_by(Chr,Start,End,P3) %>%
#   summarise(across(contains("introgressed"),.fns = sum, na.rm = T))
# 
# df1 <- df %>% 
#   mutate(IntrogressedSynTotalCount=ifelse(NumAccessionIntrogressed==0,NA,IntrogressedSynTotalCount),
#          IntrogressedNonTotalCount=ifelse(NumAccessionIntrogressed==0,NA,IntrogressedNonTotalCount),
#          IntrogressedDelTotalCount=ifelse(NumAccessionIntrogressed==0,NA,IntrogressedDelTotalCount)) %>% 
#   mutate(NonIntrogressedSynTotalCount=ifelse(NumAccessionNonIntrogressed==0,NA,NonIntrogressedSynTotalCount),
#          NonIntrogressedNonTotalCount=ifelse(NumAccessionNonIntrogressed==0,NA,NonIntrogressedNonTotalCount),
#          NonIntrogressedDelTotalCount=ifelse(NumAccessionNonIntrogressed==0,NA,NonIntrogressedDelTotalCount))
# 
# write_tsv(df1,"/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/RProject/R4VMap2FigsS1062/data/fromDaxing/fd/006_introgression_burden/001_popFdIntervalBurden/002_introgressionDonorBurden.txt.gz")



table <- read_tsv("data/fromDaxing/fd/006_introgression_burden/001_popFdIntervalBurden/002_introgressionDonorBurden.txt.gz") %>% 
  mutate(Sub=str_sub(Chr,2,2)) %>% 
  mutate(Sub=factor(Sub, levels = c("A","B","D")), SubID=str_sub(Chr,1,1)) %>% 
  mutate(IntrogressedWeaklyLoad=(IntrogressedNonTotalCount+IntrogressedDelTotalCount)/IntrogressedSynTotalCount) %>%
  mutate(IntrogressedStronglyLoad=(IntrogressedDelTotalCount)/IntrogressedSynTotalCount) %>% 
  mutate(NonIntrogressedWeaklyLoad=(NonIntrogressedNonTotalCount+NonIntrogressedDelTotalCount)/ NonIntrogressedSynTotalCount) %>% 
  mutate(NonIntrogressedStronglyLoad=(NonIntrogressedDelTotalCount)/NonIntrogressedSynTotalCount)

### 结果里有 NA  NaN Inf
dfMaxMin <- table %>% 
  pivot_longer(cols = IntrogressedWeaklyLoad:NonIntrogressedStronglyLoad,names_to = "Group",values_to = "Load") %>% 
  filter(!is.infinite(Load), !is.na(Load)) %>% 
  group_by(Sub,P3,Group) %>% 
  summarise(max=max(Load))

df <- table %>% 
  pivot_longer(cols = IntrogressedWeaklyLoad:NonIntrogressedStronglyLoad,names_to = "Group",values_to = "Load") %>% 
  left_join(y = dfMaxMin,by = c("Sub","P3","Group")) %>% 
  mutate(Load=ifelse(Load>max, max, Load)) %>%  #### 对缺失值的处理
  select(-c("max")) %>% 
  pivot_wider(names_from = "Group", values_from = "Load") %>% 
  ## Normalization
  mutate(WeaklyLoadNew2=(IntrogressedWeaklyLoad-NonIntrogressedWeaklyLoad)/(IntrogressedWeaklyLoad+NonIntrogressedWeaklyLoad),
         StronglyLoadNew2=(IntrogressedStronglyLoad-NonIntrogressedStronglyLoad)/(IntrogressedStronglyLoad+NonIntrogressedStronglyLoad)) %>% 
  pivot_longer(cols = c(WeaklyLoadNew2,StronglyLoadNew2),names_to = "SlightlyOrStrongly",values_to = "IntrogressionBurdenIndex") %>% 
  mutate(SlightlyOrStrongly=factor(SlightlyOrStrongly, levels = c("WeaklyLoadNew2","StronglyLoadNew2"), labels = c("Nonsynonymous","Deleterious")))

df$P3 <- factor(df$P3,levels = c("AT","WE","DE","FTT"))
colB <- c("#87cef9","#ffd702","#7f5701","#016699")
```


##  load ratio (Introgressed/Nonintrogressed) IntrogressionBurdenIndex
```{r}

p <- df %>% 
  ggplot(aes(x=Sub, y=IntrogressionBurdenIndex, fill=P3))+
  # geom_boxplot()+
  # geom_lv(position = position_dodge(0.8),width.method = "height",color="black",conf=0.9,k=4)+
  # stat_summary(geom = "point",fun = mean,size=1, color="black",
  #              position = position_dodge(0.8))+
  geom_hline(yintercept = 0, linetype=2, color="grey")+
  stat_summary(geom = "pointrange", fun.data = mean_cl_boot, aes(color=P3), 
               position = position_dodge(0.8))+
  labs(y="Introgressed VS\nnonintrogressed burden",x="Subgenome")+
  # scale_fill_brewer(palette = "Set2")+
  # scale_color_brewer(palette = "Set2")+
  scale_fill_manual(values = colB)+
  scale_color_manual(values = colB)+
  # scale_y_log10(breaks = scales::trans_breaks("log2", function(x) 2^x),
  #               labels = scales::trans_format("log2", scales::math_format(2^.x)))+
  # coord_cartesian(xlim = c(0,1.0))+
  facet_grid(SlightlyOrStrongly~.)+
  theme(
    # axis.title.x = element_blank(),
    strip.background = element_blank(),
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.4),
    legend.position = 'none',
    text = element_text(size = 9))
p

# ggsave(p,filename = "~/Documents/test.pdf",height = 2.4,width = 2.4)
# ggsave(p,filename = "~/Documents/test.pdf",height = 2.4,width = 3.5)
ggsave(p,filename = "~/Documents/test.pdf",height = 2.38,width = 3.59)


dftest <- df %>% 
  filter(!is.na(IntrogressionBurdenIndex)) %>% 
  group_by(P3,SlightlyOrStrongly) %>% 
  count()

dftest2 <- df %>% 
  # filter(is.na(IntrogressionBurdenIndex)) %>% 
  group_by(P3,SlightlyOrStrongly) %>% 
  count()

```

## whole genome distribution
```{r}
dfall <- df
```

```{r}
library(slider)
df <- dfall %>% 
  filter(!is.na(IntrogressionBurdenIndex)) %>% 
  group_by(Chr,P3) %>% 
  mutate(bin=slide(Start,mean,.before = 50,.after = 50) %>% unlist(),
         binBurdenIndex=slide(IntrogressionBurdenIndex, mean,.before = 50,.after = 50) %>% unlist())

#### 设置每条染色体的点
data1 <- data.frame(RefChr = paste(1:7,rep("A",7),sep = ""),x=c(213,327,319,266,254,286,362),y=rep(0,7))
data2 <- data.frame(RefChr = paste(1:7,rep("B",7),sep = ""),x=c(241,348,347,304,201,325,309),y=rep(0,7))
data3 <- data.frame(RefChr = paste(1:7,rep("D",7),sep = ""),x=c(170,268,240,185,187,215,339),y=rep(0,7))
data <- rbind(data1,data2,data3) %>% 
  mutate(Sub=str_sub(RefChr,2,2), SubID=str_sub(RefChr,1,1)) %>% 
  filter(RefChr %in% c("1A","1B","1D"))

p <- df %>%
  group_by(SlightlyOrStrongly) %>% 
  # filter(Chr %in% c("1A","1B","1D")) %>% 
  group_map(~{ggplot(.x,aes(x=bin/1000000,y=binBurdenIndex,color=P3))+
  # geom_point()+
  geom_line(size=0.4)+
  geom_hline(yintercept = 0, linetype=2, color="grey",size=0.4)+
  facet_grid(SubID~Sub)+
  geom_point(aes(x,y),color = "blue",size=0.5,data = data)+
  scale_fill_manual(values = colB)+
  scale_color_manual(values = colB)+
  scale_y_continuous(breaks = seq(-0.3,0.1,by=0.3))+
  # geom_smooth(method = "loess",span=0.2)+
  labs(x="Genomic position (Mb)", y="Introgressed VS 
       nonIntrogressed burden")+
  theme_classic()+
  theme(
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        legend.position = "none",
        legend.title = element_blank(),
        # axis.title.x = element_blank(),
        axis.line = element_line(size = 0.3),
        text = element_text(size = 9),
        panel.border = element_rect(colour = "transparent",fill = "transparent"))})
p
ggsave("~/Documents/temp.pdf",p[[1]],width =4.73 , height = 2.31)


p <- df %>%
  filter(P3=="AT") %>%
  ggplot(aes(x=bin/1000000,y=binBurdenIndex,color=SlightlyOrStrongly))+
  # geom_point()+
  geom_line()+
  geom_point(aes(x,y),color = "blue",size=1,data = data)+
  geom_hline(yintercept = 0, linetype=2, color="grey")+
  facet_grid(SubID~Sub)+
  scale_fill_manual(values = colB)+
  scale_color_manual(values = colB)+
  labs(x="Genomic position (Mb)")+
  # geom_smooth(method = "loess",span=0.2)+
  theme_classic()+
  theme(
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        # legend.position = "none", legend.title = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        # text = element_text(size = 14),
        panel.border = element_rect(colour = "black",fill = "transparent"))

p
ggsave(p,filename = "~/Documents/test.pdf",height = 7.2,width = 7.2)
  
```
#### ------------ load vs PGI vs timing
- 从这里开始一直运行，后续才能得到 load 和 PGI 以及 inrogression timing 的关系
## Calculate Relative timing (R script, using Version1)

```{r}
taxaInfo <- "data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt"
chrconvention <- "data/001_TaxaDB/chrConvertionRule_byAoyue.txt"
df_taxaInfo <- read_tsv(taxaInfo,show_col_types = FALSE) %>% 
  select(VcfID,GroupbyContinent,IntrogressionID)


library(tidyverse)
# args <- commandArgs(trailingOnly = T)

pop2indi_fd_dxy_dir <- "data/fromDaxing/015_introgression_pop2Indi/003_fd_dxy_testData/003_pop2indi_fd_dxy_v1"
outFile <- ""

# pop2indi_fd_dxy_dir <- "005_fd_dxy_data_LR_EU_AF_AM/"

df_fd_dxy_AB <- dir(pop2indi_fd_dxy_dir,full.names = T) %>% 
    keep(~str_sub(basename(.x),5,5) %in% c("A","B")) %>% 
    map(~read_tsv(.x,show_col_types = FALSE,col_types = cols(.default = "n",scaffold="c",DonorPredictedByMindxy="c",
                                                             P2="c",P3="c")) %>% 
            mutate(P2=str_c(str_split(basename(.x),"_")[[1]][[3]]))) %>% 
    bind_rows() %>% 
    mutate(DonorPredictedByMindxy_fd = ifelse(is.na(fd),"NONE",ifelse(fd > fd_mean & DonorPredictedByMindxy!="NONE", DonorPredictedByMindxy, "NONE"))) %>% 
    mutate(RelativeTiming_AC=Min_dxy/dxy_FTT_IndianDwarf) %>% 
    select(scaffold:sites,P2,DonorPredictedByMindxy_fd,mean_dxy_outgroup,Min_dxy,RNDMin,RelativeTiming_AC)

df_fd_dxy_D <- dir(pop2indi_fd_dxy_dir,full.names = T) %>% 
    keep(~str_sub(basename(.x),5,5) %in% c("D")) %>% 
    map(~read_tsv(.x,show_col_types = FALSE,col_types = cols(.default = "n",scaffold="c",DonorPredictedByMindxy="c",
                                                             P2="c",P3="c")) %>% 
            mutate(P2=str_c(str_split(basename(.x),"_")[[1]][[3]]))) %>% 
    bind_rows() %>% 
    mutate(DonorPredictedByMindxy_fd = ifelse(is.na(fd),"NONE",ifelse(fd > fd_mean & DonorPredictedByMindxy!="NONE", DonorPredictedByMindxy, "NONE"))) %>% 
    mutate(RelativeTiming_AC=Min_dxy/dxy_AT_Indian_dwarf) %>% 
    select(scaffold:sites,P2,DonorPredictedByMindxy_fd,mean_dxy_outgroup,Min_dxy,RNDMin,RelativeTiming_AC)

df_fd_dxy <- bind_rows(df_fd_dxy_AB, df_fd_dxy_D) 

df <- df_fd_dxy %>% 
    mutate(IfILS=ifelse(RelativeTiming_AC >= 0.95, 1,0)) %>% 
    mutate(DonorPredictedByMindxy_fd_ILS=ifelse(DonorPredictedByMindxy_fd!="NONE" & IfILS==0, DonorPredictedByMindxy_fd, "NONE")) %>% 
    mutate(DonorPredictedByMindxy_fd_ILS=factor(DonorPredictedByMindxy_fd_ILS,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
    select(scaffold:sites,P2,RelativeTiming_AC,DonorPredictedByMindxy_fd_ILS) %>% 
    filter(!DonorPredictedByMindxy_fd_ILS %in% c("NONE")) %>% 
    mutate(Sub=str_sub(scaffold,2,2)) %>% 
    group_by(P2,DonorPredictedByMindxy_fd_ILS,Sub) %>% 
    summarise(mean_RelativeTiming_AC=mean(RelativeTiming_AC,na.rm=T),
              sd=sd(RelativeTiming_AC,na.rm=T),
              WindowNum=n())

head(df)
# write_tsv(df, outFile)
```

### Distribution of Relative timing

```{r}
df_taxaInfo <- read_tsv(taxaInfo) %>% 
  select(IntrogressionID,GroupbyContinent)

df1 <- read_tsv("data/fromDaxing/015_introgression_pop2Indi/009_relativingTiming/RelativeTiming_AC_LR_EA_CL.txt.gz") %>% 
  mutate(P2=str_sub(P2,1,6)) %>% 
  mutate(P2=str_c(str_sub(P2,1,2),str_sub(P2,3,6),sep = "_")) %>% 
  left_join(y = df_taxaInfo,by = c("P2"="IntrogressionID"))

df2 <- read_tsv("data/fromDaxing/015_introgression_pop2Indi/009_relativingTiming/RelativeTiming_AC_LR_EU_AF_AM.txt.gz") %>% 
  mutate(P2=str_c(str_sub(P2,1,2),str_sub(P2,3,6),sep = "_")) %>% 
  left_join(y = df_taxaInfo,by = c("P2"="IntrogressionID"))

df <- bind_rows(df1,df2)
dfAB <- df %>% 
  filter(Sub %in% c("A","B")) %>% 
  group_by(P2,DonorPredictedByMindxy_fd_ILS) %>% 
  mutate(proportion=WindowNum/sum(WindowNum),Timing=mean_RelativeTiming_AC*proportion) %>% 
  summarise(Time=sum(Timing))

dfD <- df %>% 
  filter(Sub %in% c("D")) %>% 
  select(P2,DonorPredictedByMindxy_fd_ILS,mean_RelativeTiming_AC) %>% 
  rename(Time=mean_RelativeTiming_AC)
  
dfTiming <- bind_rows(dfAB,dfD) %>% 
  left_join(y = df_taxaInfo,by = c("P2"="IntrogressionID")) %>% 
  filter(GroupbyContinent %in% c("LR_EU","LR_EA","CL"))

# factorA <- c("AT","WE","DE","FTT")
# labelsA <- c("Strangulata","Wild emmer","Domesticated emmer","Free-threshing tetraploids")
# colB <- c("#87cef9","#ffd702","#7f5701","#016699")

factorA <- c("LR_EU","CL","LR_EA")
labelsA <- c("Landrace_EU","Cultivar","Landrace_EA")
colB <- c("#fc6e6e","#9900ff","gray90")


## AB
q <- dfTiming %>% 
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(DonorPredictedByMindxy_fd_ILS=factor(DonorPredictedByMindxy_fd_ILS,levels = c("WE","AT","DE","FTT"))) %>%
  filter(!DonorPredictedByMindxy_fd_ILS %in% c("AT")) %>%
  ggplot(aes(x=Time, fill=GroupbyContinent, color=GroupbyContinent))+
  # geom_histogram(alpha=0.6,position = position_identity())+
  geom_density(alpha=0.6)+
  facet_grid(DonorPredictedByMindxy_fd_ILS~.)+
  scale_fill_manual(values = colB)+
  scale_color_manual(values = colB)+
  labs(x="Relative timing of introgression")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        # text = element_text(size = 16),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
ggsave("~/Documents/temp.png",q,width = 4,height = 3)


## D
q <- dfTiming %>% 
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(DonorPredictedByMindxy_fd_ILS=factor(DonorPredictedByMindxy_fd_ILS,levels = c("WE","AT","DE","FTT"))) %>%
  filter(DonorPredictedByMindxy_fd_ILS %in% c("AT")) %>%
  ggplot(aes(x=Time, fill=GroupbyContinent, color=GroupbyContinent))+
  # geom_histogram(alpha=0.6,position = position_identity())+
  geom_density(alpha=0.6)+
  facet_grid(DonorPredictedByMindxy_fd_ILS~.)+
  scale_fill_manual(values = colB)+
  scale_color_manual(values = colB)+
  labs(x="Relative timing of introgression")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        # text = element_text(size = 16),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q

```

## Introgression time PGI VS mutation burden
## Scatter burden vs Timing vs PGI

**Burden, WildAncestryPerTaxa, Time** **不区分亚基因组**
- WE_Burden是AB亚基因组burden之和
- WE_WildAncestryPerTaxa是AB亚基因组WildAncestryPerTaxa之和
- WE_Time是AB基因组所有window的均值

### Tidy

```{r}
taxaInfo <- "data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt"
chrconvention <- "data/001_TaxaDB/chrConvertionRule_byAoyue.txt"
df_taxaInfo <- read_tsv(taxaInfo,show_col_types = FALSE) %>% 
  select(VcfID,GroupbyContinent,IntrogressionID)

genomeSize <- 14066280851

factorA <- c("LR_EU","CL","LR_EA")
labelsA <- c("Landrace_EU","Cultivar","Landrace_EA")
colB <- c("#fc6e6e","#9900ff","gray90")

df_chrconvention <- read_tsv(chrconvention,show_col_types = FALSE) %>% select(ChrID,OriChrName) %>% 
  mutate(Chr=str_remove(OriChrName,"chr"), Sub=str_sub(Chr,2,2)) %>% 
  select(-OriChrName) 


df_Time_PGI <- read_tsv("data/fromDaxing/015_introgression_pop2Indi/007_PGI/PGI.txt.gz") %>% 
  left_join(y = df_taxaInfo,by = c("P2"="IntrogressionID")) %>% 
  mutate(Genomesize=genomeSize) %>% 
  mutate(WildAncestryPerTaxa=Genomesize*PGI/1000000) %>% 
  group_by(GroupbyContinent,P2,DonorPredictedByMindxy_fd_ILS) %>% 
  summarise(TotalWildAncestryPerTaxa=sum(WildAncestryPerTaxa)) %>% 
  filter(!DonorPredictedByMindxy_fd_ILS %in% c("NONE")) %>%
  # filter(GroupbyContinent %in% c("LR_EU","LR_EA","CL")) %>%
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = factorA)) %>% 
  left_join(y = dfTiming,by = c("P2","DonorPredictedByMindxy_fd_ILS","GroupbyContinent"))


df_Burden_Time_PGI <- dir("data/fromDaxing/015_introgression_pop2Indi/008_fd_burden/002_individualBased/",full.names = T) %>% 
  map(read_tsv,show_col_types = FALSE) %>% 
  bind_rows() %>% 
  left_join(y = df_chrconvention, by = "ChrID") %>% 
  group_by(Taxa,Donor) %>%
  summarise(across(starts_with("num"),sum)) %>%
  mutate(Slightly=(numDerivedInNonsyn+0.5*numHeterInNonsyn)/(numDerivedInSyn+0.5*numHeterInSyn)) %>% 
  mutate(Strongly=(numDerivedInHGDeleterious+0.5*numHeterInHGDeleterious)/(numDerivedInSyn+0.5*numHeterInSyn)) %>% 
  select(-starts_with("num")) %>% 
  pivot_longer(cols = Slightly:Strongly, names_to = "SlightlyOrStrongly", values_to = "Burden") %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  left_join(y = df_taxaInfo, by = c("Taxa"="VcfID")) %>% 
  # ABD_0862 ABD_0863 is NA
  filter(!is.na(GroupbyContinent)) %>% 
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_AM","LR_AF","LR_CSA","LR_WA","LR_EU","CL","LR_EA"))) %>% 
  # filter(GroupbyContinent  %in% c("LR_EU","LR_EA","CL")) %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  left_join(y = df_Time_PGI,by = c("IntrogressionID"="P2","Donor"="DonorPredictedByMindxy_fd_ILS","GroupbyContinent"))


```

### Scatter {.tabset}

#### PGI vs Timing

```{r}
## PGI and timing
library(scales)
factorA <- c("AT","WE","DE","FTT")
labelsA <- c("Strangulata","Wild emmer","Domesticated emmer","Free-threshing tetraploids")
colB <- c("#87cef9","#ffd702","#7f5701","#016699","grey")

q <- df_Burden_Time_PGI %>% 
    filter(GroupbyContinent  %in% c("LR_EU","LR_EA","CL")) %>%
    filter(SlightlyOrStrongly=="Slightly") %>%
  # filter(!IntrogressionID %in% c("CL_0056","CL_0187","LR_0101","CL_0025","LR_0040","LR_0096")) %>%
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  ggplot(aes(y=TotalWildAncestryPerTaxa,x=Time,label=GroupbyContinent,color=Donor,shape=GroupbyContinent))+
  geom_point(alpha=0.6)+
    # facet_grid(.~GroupbyContinent,scales = "free")+
    scale_y_continuous(name = "Wild ancestry per individual (Mb)",
                     trans = log10_trans(),
                     # limits = c(0, 4000),
                     breaks = c(1,10,100, 1000,10000),
                     labels = trans_format("log10", math_format(10^.x)))+
  scale_x_continuous(breaks = c(0.3,0.4,0.5,0.6))+
    annotation_logticks(sides = "l",short = unit(0.5,units = "mm"),
                      mid = unit(1,units = "mm"),long = unit(2,units = "mm"))+
  scale_color_manual(values = colB)+
  labs(x="Relative timing of introgression")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 9),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
# ggsave("~/Documents/temp.png",q,width = 4,height = 3)
ggsave(q,filename = "~/Documents/test.pdf",height = 2.4,width = 3.5)
# ggsave(q,filename = "~/Documents/test.pdf",height = 3,width = 4)

# ggplotly(q)

q <- df_Burden_Time_PGI %>% 
    filter(GroupbyContinent  %in% c("LR_EU","LR_EA","CL")) %>%
  filter(SlightlyOrStrongly=="Slightly") %>%
  # filter(!IntrogressionID %in% c("CL_0056","CL_0187","LR_0101","CL_0025","LR_0040","LR_0096")) %>%
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  ggplot(aes(y=TotalWildAncestryPerTaxa,x=Time,label=GroupbyContinent,color=Donor,shape=GroupbyContinent))+
  geom_point(alpha=0.6)+
  facet_grid(.~GroupbyContinent,scales = "free")+
  scale_y_continuous(name = "Wild ancestry per individual (Mb)",
                     trans = log10_trans(),
                     # limits = c(0, 4000),
                     breaks = c(1,10,100, 1000,10000),
                     labels = trans_format("log10", math_format(10^.x)))+
  scale_x_continuous(breaks = c(0.3,0.4,0.5,0.6))+
  annotation_logticks(sides = "l",short = unit(0.5,units = "mm"),
                      mid = unit(1,units = "mm"),long = unit(2,units = "mm"))+
  scale_color_manual(values = colB)+
  labs(x="Relative timing of introgression")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        # text = element_text(size = 16),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
ggsave("~/Documents/temp.png",q,width = 4,height = 3)

```


#### Burden vs PGI

```{r}
## burden and PGI

factorA <- c("AT","WE","DE","FTT")
labelsA <- c("Strangulata","Wild emmer","Domesticated emmer","Free-threshing tetraploids")
colB <- c("#87cef9","#ffd702","#7f5701","#016699","grey")


q <- df_Burden_Time_PGI %>% 
    filter(GroupbyContinent  %in% c("LR_EU","LR_EA","CL")) %>%
  # filter(SlightlyOrStrongly=="Strongly") %>%
  # filter(!IntrogressionID %in% c("CL_0056","CL_0187","LR_0101","CL_0025","LR_0040","LR_0096")) %>%
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  ggplot(aes(y=Burden,x=TotalWildAncestryPerTaxa,label=IntrogressionID,color=Donor))+
  geom_point(alpha=0.4)+
  # facet_grid(.~GroupbyContinent,scales = "free")+
  facet_grid(SlightlyOrStrongly~.,scales = "free")+
  scale_x_continuous(name = "Wild ancestry per individual (Mb)",
                     trans = log10_trans(), 
                     # limits = c(0, 4000), 
                     breaks = c(1,10,100, 1000,10000),
                     labels = trans_format("log10", math_format(10^.x)))+
  annotation_logticks(sides = "b",short = unit(0.5,units = "mm"),
                      mid = unit(1,units = "mm"),long = unit(2,units = "mm"))+
  # coord_cartesian(clip = "off")+
  scale_color_manual(values = colB)+
  labs(y="Mutation burden")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 9),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
# ggsave("~/Documents/temp.png",q,width = 4,height = 3)
ggsave(q,filename = "~/Documents/test.pdf",height = 2.4,width = 3.5)
# ggsave(q,filename = "~/Documents/test.pdf",height = 3,width = 4)


q <- df_Burden_Time_PGI %>% 
    filter(GroupbyContinent  %in% c("LR_EU","LR_EA","CL")) %>%
  # filter(SlightlyOrStrongly=="Strongly") %>%
  # filter(!IntrogressionID %in% c("CL_0056","CL_0187","LR_0101","CL_0025","LR_0040","LR_0096")) %>%
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  ggplot(aes(y=Burden,x=TotalWildAncestryPerTaxa,label=IntrogressionID,color=Donor))+
  geom_point(alpha=0.4)+
  facet_grid(SlightlyOrStrongly~GroupbyContinent,scales = "free")+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  scale_x_continuous(name = "Wild ancestry per individual (Mb)",
                     trans = log10_trans(), 
                     # limits = c(0, 4000), 
                     breaks = c(1,10,100, 1000,10000),
                     labels = trans_format("log10", math_format(10^.x)))+
    annotation_logticks(sides = "b",short = unit(0.5,units = "mm"),
                      mid = unit(1,units = "mm"),long = unit(2,units = "mm"))+
  # coord_cartesian(clip = "off")+
  scale_color_manual(values = colB)+
  labs(y="Mutation burden")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        # text = element_text(size = 16),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q

```

#### Burden vs Timing

```{r}
## Burden and timing
q <- df_Burden_Time_PGI %>% 
    filter(GroupbyContinent  %in% c("LR_EU","LR_EA","CL")) %>%
  # filter(SlightlyOrStrongly=="Slightly") %>%
  # filter(!IntrogressionID %in% c("CL_0056","CL_0187","LR_0101","CL_0025","LR_0040","LR_0096")) %>%
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  ggplot(aes(y=Burden,x=Time,label=IntrogressionID,color=Donor))+
  geom_point(alpha=0.6)+
  # facet_grid(.~GroupbyContinent,scales = "free")+
  facet_grid(SlightlyOrStrongly~.,scales = "free")+
  scale_color_manual(values = colB)+
  labs(y="Mutation burden",x="Relative timing of introgression")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 9),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
ggsave(q,filename = "~/Documents/test.pdf",height = 2.4,width = 3.5)

# ggsave(q,filename = "~/Documents/test.pdf",height = 3,width = 4)


q <- df_Burden_Time_PGI %>% 
    filter(GroupbyContinent  %in% c("LR_EU","LR_EA","CL")) %>%
  # filter(SlightlyOrStrongly=="Slightly") %>%
  filter(!IntrogressionID %in% c("CL_0056","CL_0187","LR_0101","CL_0025","LR_0040","LR_0096")) %>%
  mutate(GroupbyContinent=factor(GroupbyContinent,levels = c("LR_EU","CL","LR_EA"))) %>% 
  mutate(Donor=factor(Donor,levels = c("AT","WE","DE","FTT","NONE"))) %>% 
  ggplot(aes(y=Burden,x=Time,label=IntrogressionID,color=Donor))+
  geom_point(alpha=0.6)+
  facet_grid(SlightlyOrStrongly~GroupbyContinent,scales = "free")+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  scale_color_manual(values = colB)+
  labs(y="Mutation burden",x="Relative timing of introgression")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.line = element_blank(),
        # text = element_text(size = 16),
        # axis.title.x = element_blank(),
        # axis.text.x = element_blank(),
        # axis.ticks.x = element_blank(),
        # axis.title = element_text(size = 8),
        # strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q

ggsave("~/Documents/temp.png",q,width = 4,height = 3)

```



#### ===============
## ***Extended Data***
#### ===============
## FigS:XPCLR positive control
### Important data
```{r}
dataFile <- c("data/020_xpclr/ab_WE_DE_log2_0.0001_200_100000.clrgs.txt.gz")
df_wede <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/ab_DE_Durum_log2_0.0001_200_100000.clrgs.txt.gz")
df_dedm <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/abd_log2_0.0001_200_100000.clrgs.txt.gz")
df_lrcl <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

df_wede <- df_wede %>% mutate(group="WE vs DE", subID=substr(ChrRef,1,1), sub=substr(ChrRef,2,2))
df_dedm <- df_dedm %>% mutate(group="DE vs Durum", subID=substr(ChrRef,1,1), sub=substr(ChrRef,2,2))
df_lrcl <- df_lrcl %>% mutate(group="LR vs CL", subID=substr(ChrRef,1,1), sub=substr(ChrRef,2,2))

df <- rbind(df_wede,df_dedm,df_lrcl)
df$group <- factor(df$group, levels = c("WE vs DE","DE vs Durum","LR vs CL"))

```

### Start to plot - 循环
- 一一验证比较慢，故采用循环的方式，瞬间出来
```{r}
df_positive <- read_tsv("data/020_xpclr/GeneList_forPositiveControl.txt")

aoFactor <- c("Btr-A1","Btr-B1","Q","Tg","Rht-B1","Rht-D1","Ppd-D1","Vrn-A1","Vrn-B1","Vrn-D1","Fhb1","Sm1","Lr10","Lr21","Lr34-Yr18-Sr57-Pm38-Ltn1-Sb1")
for (i in aoFactor) {
  dfgeneinfo <- df_positive %>% filter(GeneName==i)
  
  groupA <- dfgeneinfo$Group
  chrRefA <- dfgeneinfo$ChrRef
  genePosA <- dfgeneinfo$GenePos
  centromerePosA <- dfgeneinfo$CentromerePos
  xlabA <- dfgeneinfo$XLab
  
   p <- df %>% 
    filter(ChrRef == chrRefA, group == groupA) %>%
     ggplot(aes(x=PosRef/1000000, y=XPCLR_100score,col=group))+
     geom_line(size=0.5)+
     geom_point(aes(x=centromerePosA,y=0),color='grey')+ ### 标记着丝粒
    # geom_hline(yintercept = 19.48245)+
    ### 标记蓝色点
    geom_point(aes(x=genePosA,y=0),color='blue')+ ###标记基因
    # geom_point(aes(x=gene2,y=0),color='blue')+
    labs(x=xlabA,y="XP-CLR score")+
    facet_wrap(~group,scales = "free_x",ncol = 1)+
    # scale_x_continuous(breaks = seq(0,600,by = 300))+
    # scale_x_continuous(limits = c(0,50), breaks = seq(0,25,by = 5))+
    scale_y_continuous(breaks = seq(0,40,by = 25))+
    # scale_color_manual(values = colB)+
    scale_color_brewer(palette = "Set2")+
    theme_classic()+
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()
      ,panel.background = element_blank()
      ,axis.line = element_line(size=0.3, colour = "black")
      ,legend.position = 'none'
      ,text = element_text(size = 12))+
    theme(strip.background = element_rect(fill = "lightgray",color = "white"))

  p
  
  filename <- str_c("~/Documents/",i,".png",sep = "")


# ggsave(p,filename = "~/Documents/test.pdf",height = 2,width = 5)
ggsave(p,filename = filename,height = 2,width = 5) ###由于线条过多，故保存成图片
}

```

#### ======
## xpclr-whole genome
## ***XPCLR score
## WE vs. DE manhattan plot for add signal dot

```{r}
################################# ##
#                               #
#  XPCLR genome line -- manhattan 100kb 50 kb 
#                               #
################################# ##
library(tidyverse)
# ### log2 snp method whole-genome SNP 3%
dataFile <- c("data/020_xpclr/ab_WE_DE_log2_0.0001_200_100000.clrgs.txt.gz")
df_wede <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/ab_DE_Durum_log2_0.0001_200_100000.clrgs.txt.gz")
df_dedm <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/abd_log2_0.0001_200_100000.clrgs.txt.gz")
df_lrcl <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

### ******************************************************************************
### choice: 1.df 2.color 3.threshold  
### ******************************************************************************
### 1. data
df <- df_wede
# df <- df_dedm
# df <- df_lrcl
### 2.color
test <- "wede"
# test <- "dedm"
# test <- "lrcl"
colB <- NULL
labelB <- NULL

# colB1 <- c("#ffd702","#7f5701") ### we de
# colB2 <- c("#7f5701","#016699") ### de durum
# colB3 <- c("#fc6e6e","#9900ff")  ### lr cl
## use subgenome color
colB1 <- c("#fd8582","#967bce") ### we de
colB2 <- c("#fd8582","#967bce") ### de durum
colB3 <- c("#fd8582","#967bce","#4bcdc6")  ### lr cl

### if else failed to use
if(test == "wede"){
  colB <- colB1
  labelB <- "WE vs. DE"
} 
if(test == "dedm"){
  colB <- colB2
  labelB <- "DE vs. Durum"
}
if(test == "lrcl"){
  colB <- colB3
  labelB <- "Landrace vs. Cultivar"
}

### 3. 
# top <- 0.01
top <- 0.05

# head(df)
# nrow(df)
### step2: 修改列名，使之固定
names(df)[names(df)=="XPCLR_100score"] = "Ave"
names(df)[names(df)=="CHROM"] = "Chr"
names(df)[names(df)=="ChrRef"] = "CHROM"
names(df)[names(df)=="PosRef"] = "BIN_START"
max(df$Ave)
min(df$Ave)
### step3: 去除NA和无限小的值
df <- dplyr::filter(df,!is.na(Ave) & is.finite(Ave))
nrow(df)
head(df)
max(df$Ave) # 48.50477
min(df$Ave) # 0
### step4: 获取TOP K 截距
### step4a: 排序
df_sorted <- df %>% 
  dplyr::arrange(-Ave) 
head(df_sorted)
### step4b: 获取top1% 的值，即水平截距
threshod <- df_sorted[top*nrow(df_sorted),8]
threshod

# write.table(d,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 
### step5: 建立累积坐标系统和x轴标签位置
don <- df %>%
  # Compute chromosome size
  group_by(CHROM) %>%
  summarise(chr_len=max(BIN_START)) %>%  ##每条染色体的最大值赋值给 chr_len
  # Calculate cumulative position of each chromosome
  mutate(tot=cumsum(as.numeric(chr_len))-chr_len) %>%  ##所有之和减去该行的染色体对应的长度，如：1B起始为1A长度，1D起始为1A 1B之和的长度
  select(-chr_len) %>% ##去除 chr_len这一列
  # Add this info to the initial dataset
  left_join(df, ., by=c("CHROM"="CHROM")) %>%  ## 将原数据框和don数据框进行合并
  # Add a cumulative position of each SNP
  arrange(CHROM, BIN_START) %>% ##根据CHROM BIN_START 进行排序
  mutate( BPcum=BIN_START+tot)  ##增加新的一列，改变每行点的位置

### 记录增量
dfadd <- df %>%
  # Compute chromosome size
  group_by(CHROM) %>%
  summarise(chr_len=max(BIN_START)) %>%  ##每条染色体的最大值赋值给 chr_len
  # Calculate cumulative position of each chromosome
  mutate(tot=cumsum(as.numeric(chr_len))-chr_len) %>%  ##所有之和减去该行的染色体对应的长度，如：1B起始为1A长度，1D起始为1A 1B之和的长度
  select(-chr_len)  ##去除 chr_len这一列
# dfadd  

## x轴标签位置：求出每条染色体在坐标轴上的中间位置
axisdf = don %>% group_by(CHROM) %>% summarize(center=( max(BPcum) + min(BPcum) ) / 2 )

### 建立函数关系，添加驯化信号的位点
### 建立函数关系，添加驯化信号的位点
### 建立函数关系，添加驯化信号的位点
geneList <- read.table("data/020_xpclr/gene.txt",sep = "\t",header = T)
geneList
geneList$Gene <- as.vector(geneList$Gene)
geneList$Chromosome <- as.vector(geneList$Chromosome)
geneList$Start <- as.numeric(geneList$Start)  ## 为了可以进行预测
geneList$End <- as.numeric(geneList$End)

### 写一个函数，返回基因的标签，以及在图上的range位置
aoGetGenepos <- function(geneRow){
  gene <- geneList[geneRow,] ### get gene's dataframe, only has one line
  labelgene <- gene[1,1]### get gene's label
  chrGene <- as.vector(gene[1,2]) ### get gene's chr
  posSignal <- gene[,3:4] ### get gene's pos, a dataframe
  posOnchr <- c(posSignal[1,1],posSignal[1,2]) ### change to x 
  ### select chr, only one line
  dfadd_one <- dfadd %>% 
    filter(CHROM == chrGene) 
  # dfadd_one
  posSignalCum <- posOnchr + as.numeric(dfadd_one[1,2])
  # geneSummary <- list( labelgene,as.numeric(posSignalCum[1]),as.numeric(posSignalCum[2]))
  geneSummary <- list( labelgene,posSignalCum[1],posSignalCum[2])
  
  return(geneSummary)
}
 
#### btr1 
btra1 <- aoGetGenepos(5)
btrb1 <- aoGetGenepos(6)
btra2 <- aoGetGenepos(8)
btrb2 <- aoGetGenepos(9)

### step6: 颜色设置并画图
### 开始作图

### sample
don <- sample_frac(don,0.5)
p <- ggplot(don, aes(x=BPcum, y=Ave)) +  ##新的横坐标，Y轴的XPCLR值不变
  # Show all points
  # geom_point( aes(color=as.factor(CHROM)), alpha=0.2, size=0.5) +  ## alpha=0.5 first test
  geom_line(aes(color=as.factor(CHROM)), alpha=0.9, size=0.2) +  ## alpha=0.5 first test
  annotate("text",x=Inf,y= threshod ,hjust= 2,vjust= -0.3,label="top 5%", size = 3)+
  geom_hline(yintercept = threshod, linetype= "dashed", color = "gray",size=0.2)+
  annotate("text",x=0,y= Inf ,hjust= -0.2,vjust= 1.1,label= labelB, size = 2)+ ## add group annotation
  scale_color_manual(values = rep(colB, 14 )) + ##合计42条染色体。故21*2=42
  # custom X axis: 每条染色体对应一个中间位置，完美！
  scale_x_continuous(expand = c(0,0), label = axisdf$CHROM, breaks= axisdf$center) +
  scale_y_continuous(expand = c(0,1),limits = c(-6,54)) +
  # scale_y_continuous(expand = c(0, 0),limits = c(-1,52) ) +     # remove space between plot area and x axis
  
  ## add arrow to represent the selection signal 
  # geom_vline(xintercept = unlist(btra1[2]),color="black",size=0.2)+
  # geom_segment(aes(x=unlist(btra1[2]) ,xend = unlist(btra1[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(btra1[2]) ,xend = unlist(btra1[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") + 
  annotate(geom = "text",x = unlist(btra1[3]),y=0.2,label= unlist(btra1[1]),fontface = 'italic',hjust = -0.2,size=1.4,color="blue") +

  # geom_segment(aes(x=unlist(btrb1[2]) ,xend = unlist(btrb1[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(btrb1[2]),xend = unlist(btrb1[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(btrb1[3]),y=0.2,label=unlist(btrb1[1]),fontface = 'italic',hjust = -0.2,size=1.4,color="blue") +
  
  # annotate(geom = "segment",x=unlist(btra2[2]),xend = unlist(btra2[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(btra2[3]),y=-4,label=unlist(btra2[1]),fontface = 'italic',hjust = -0.3,size=1.4,color="blue") +
  # geom_vline(xintercept = unlist(btra2[2]),color="red",size=0.5)+
  
  # annotate(geom = "segment",x=unlist(btrb2[2]),xend = unlist(btrb2[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(btrb2[3]),y=-4,label=unlist(btrb2[1]),fontface = 'italic',hjust = -0.3,size=1.4,color="blue") +
  # geom_vline(xintercept = unlist(btrb2[2]),color="red",size=0.5)+
  
   # Custom the theme:
  xlab("")+ylab("XP-CLR")+
  theme(plot.margin=unit(rep(1,4),'lines'))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 9),legend.position = 'none')

p
# ggsave(p,filename = "/Users/Aoyue/Documents/test.png",width = 7,height = 7*3/16)
ggsave(p,filename = "/Users/Aoyue/Documents/test.pdf",width = 7,height = 7*3/16)

```

```{r eval=FALSE, include=FALSE}
### 最初添加pos 设计代码
## 选择 btr2, 某一行依旧为数据框
btrb1 <- geneList[6,]
btrb1
gene <- btrb1
gene
labelBtrb1 <- gene[1,1] 
chrGene <- as.vector(gene[1,2])  
posSignal <- gene[,3:4]
posOnplot <- data.frame(x=c(posSignal[1,1],posSignal[1,2])) 
### select chr
donSignal <- don %>% 
  filter(CHROM == chrGene) 
donSignal
y <- donSignal$BPcum
x <- donSignal$BIN_START
### build regression
relation <- lm( y ~ x )
relation
# summary(relation)
### input data frame
posSignalCum <-  predict(relation,posOnplot)
posSignalCum


  ## add arrow to represent the selection signal 
  # geom_segment(data=genes, mapping=aes(x=x, y= -5, xend=xend, yend=0), arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="red") + 
  # annotate(geom = "segment",x=posSignalCum[1],xend = posSignalCum[2],y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="red") + 
  # annotate(geom = "text",x = posSignalCum[2],y=-3,label=labelBtrb1,fontface = 'italic',hjust = -0.5,size=1.2,color="red") 

```

## DE vs. Durum manhattan plot for add signal dot

```{r}
################################# ##
#                               #
#  XPCLR genome line -- manhattan 100kb 50 kb 
#                               #
################################# ##
library(tidyverse)

dataFile <- c("data/020_xpclr/ab_WE_DE_log2_0.0001_200_100000.clrgs.txt.gz")
df_wede <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/ab_DE_Durum_log2_0.0001_200_100000.clrgs.txt.gz")
df_dedm <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/abd_log2_0.0001_200_100000.clrgs.txt.gz")
df_lrcl <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

# df <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

### ******************************************************************************
### choice: 1.df 2.color 3.threshold  
### ******************************************************************************
### 1. data
# df <- df_wede
df <- df_dedm
# df <- df_lrcl
### 2.color
# test <- "wede"
test <- "dedm"
# test <- "lrcl"
colB <- NULL
labelB <- NULL

# colB1 <- c("#ffd702","#7f5701") ### we de
# colB2 <- c("#7f5701","#016699") ### de durum
# colB3 <- c("#fc6e6e","#9900ff")  ### lr cl
## use subgenome color
colB1 <- c("#fd8582","#967bce") ### we de
colB2 <- c("#fd8582","#967bce") ### de durum
colB3 <- c("#fd8582","#967bce","#4bcdc6")  ### lr cl

### if else failed to use
if(test == "wede"){
  colB <- colB1
  labelB <- "WE vs. DE"
} 
if(test == "dedm"){
  colB <- colB2
  labelB <- "DE vs. Durum"
}
if(test == "lrcl"){
  colB <- colB3
  labelB <- "Landrace vs. Cultivar"
}

### 3. 
 # top <- 0.01
top <- 0.05

# head(df)
# nrow(df)
### step2: 修改列名，使之固定
names(df)[names(df)=="XPCLR_100score"] = "Ave"
names(df)[names(df)=="CHROM"] = "Chr"
names(df)[names(df)=="ChrRef"] = "CHROM"
names(df)[names(df)=="PosRef"] = "BIN_START"
max(df$Ave)
min(df$Ave)
### step3: 去除NA和无限小的值
df <- dplyr::filter(df,!is.na(Ave) & is.finite(Ave))
nrow(df)
head(df)
max(df$Ave) # 48.50477
min(df$Ave) # 0
### step4: 获取TOP K 截距
### step4a: 排序
df_sorted <- df %>% 
  dplyr::arrange(-Ave) 
head(df_sorted)
### step4b: 获取top1% 的值，即水平截距
threshod <- df_sorted[top*nrow(df_sorted),8]
 print(paste(threshod," threshod ", top))

# write.table(d,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 
### step5: 建立累积坐标系统和x轴标签位置
don <- df %>%
  # Compute chromosome size
  group_by(CHROM) %>%
  summarise(chr_len=max(BIN_START)) %>%  ##每条染色体的最大值赋值给 chr_len
  # Calculate cumulative position of each chromosome
  mutate(tot=cumsum(as.numeric(chr_len))-chr_len) %>%  ##所有之和减去该行的染色体对应的长度，如：1B起始为1A长度，1D起始为1A 1B之和的长度
  select(-chr_len) %>% ##去除 chr_len这一列
  # Add this info to the initial dataset
  left_join(df, ., by=c("CHROM"="CHROM")) %>%  ## 将原数据框和don数据框进行合并
  # Add a cumulative position of each SNP
  arrange(CHROM, BIN_START) %>% ##根据CHROM BIN_START 进行排序
  mutate( BPcum=BIN_START+tot)  ##增加新的一列，改变每行点的位置

### 记录增量
dfadd <- df %>%
  # Compute chromosome size
  group_by(CHROM) %>%
  summarise(chr_len=max(BIN_START)) %>%  ##每条染色体的最大值赋值给 chr_len
  # Calculate cumulative position of each chromosome
  mutate(tot=cumsum(as.numeric(chr_len))-chr_len) %>%  ##所有之和减去该行的染色体对应的长度，如：1B起始为1A长度，1D起始为1A 1B之和的长度
  select(-chr_len)  ##去除 chr_len这一列
dfadd  
  
# don
# write.table(don,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t")

## x轴标签位置：求出每条染色体在坐标轴上的中间位置
axisdf = don %>% group_by(CHROM) %>% summarize(center=( max(BPcum) + min(BPcum) ) / 2 )

### 建立函数关系，添加驯化信号的位点
### 建立函数关系，添加驯化信号的位点
### 建立函数关系，添加驯化信号的位点
geneList <- read.table("data/020_xpclr/gene.txt",sep = "\t",header = T)
geneList
geneList$Gene <- as.vector(geneList$Gene)
geneList$Chromosome <- as.vector(geneList$Chromosome)
geneList$Start <- as.numeric(geneList$Start)  ## 为了可以进行预测
geneList$End <- as.numeric(geneList$End)

### 写一个函数，返回基因的标签，以及在图上的range位置
aoGetGenepos <- function(geneRow){
  gene <- geneList[geneRow,] ### get gene's dataframe, only has one line
  labelgene <- gene[1,1]### get gene's label
  chrGene <- as.vector(gene[1,2]) ### get gene's chr
  posSignal <- gene[,3:4] ### get gene's pos, a dataframe
  posOnchr <- c(posSignal[1,1],posSignal[1,2]) ### change to x 
  ### select chr, only one line
  dfadd_one <- dfadd %>% 
    filter(CHROM == chrGene) 
  # dfadd_one
  posSignalCum <- posOnchr + as.numeric(dfadd_one[1,2])
  # geneSummary <- list( labelgene,as.numeric(posSignalCum[1]),as.numeric(posSignalCum[2]))
  geneSummary <- list( labelgene,posSignalCum[1],posSignalCum[2])
  
  return(geneSummary)
}
 
#### btr1 
q <- aoGetGenepos(1)
q2 <- aoGetGenepos(2)
tg <- aoGetGenepos(4)

### step6: 颜色设置并画图
### 开始作图

### sample
don <- sample_frac(don,0.5)
p <- ggplot(don, aes(x=BPcum, y=Ave)) +  ##新的横坐标，Y轴的XPCLR值不变
  # Show all points
  # geom_point( aes(color=as.factor(CHROM)), alpha=0.2, size=0.5) +  ## alpha=0.5 first test
  geom_line(aes(color=as.factor(CHROM)), alpha=0.9, size=0.2) +  ## alpha=0.5 first test
  annotate("text",x=Inf,y= threshod ,hjust= 2,vjust= -0.3,label="top 5%", size = 3)+
  geom_hline(yintercept = threshod, linetype= "dashed", color = "gray",size=0.2)+
  annotate("text",x=0,y= Inf ,hjust= -0.2,vjust= 1.1,label= labelB, size = 2)+ ## add group annotation
  scale_color_manual(values = rep(colB, 14 )) + ##合计42条染色体。故21*2=42
  # custom X axis: 每条染色体对应一个中间位置，完美！
  scale_x_continuous(expand = c(0,0), label = axisdf$CHROM, breaks= axisdf$center) +
  scale_y_continuous(expand = c(0,1),limits = c(-6,54)) +
  # scale_y_continuous(expand = c(0, 0),limits = c(-1,52) ) +     # remove space between plot area and x axis
  
  ## add arrow to represent the selection signal 
  # geom_vline(xintercept = 5838260571,color="black",size=0.2)+
  # geom_segment(aes(x=6490914897,xend = 6490914897,y=0,yend = Inf) , size=0.2, color="blue") +
  # annotate(geom = "segment",x=6490914897 ,xend = 6490914897,y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") + 
  # annotate(geom = "text",x = 6490914897,y=-3,label= unlist(q[1]),fontface = 'italic',hjust = -0.4,size=1.2,color="blue") +
  # 
  # geom_segment(aes(x=7202289667 ,xend = 7202289667,y=0,yend = Inf) , size=0.2, color="blue") +
  # annotate(geom = "segment",x=7202289667,xend = 7202289667,y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  # annotate(geom = "text",x = 7202289667,y=-3,label=unlist(q2[1]),fontface = 'italic',hjust = -0.4,size=1.2,color="blue") +
  # # 
  # geom_segment(aes(x=2076798318 ,xend = 2076798318,y=0,yend = Inf) , size=0.2, color="blue") +
  # annotate(geom = "segment",x=2076798318,xend = 2080198318,y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  # annotate(geom = "text",x = 2080198318,y=-3,label=unlist(tg[1]),fontface = 'italic',hjust = -0.4,size=1.2,color="blue") +

  
  # geom_segment(aes(x=unlist(q[2]) ,xend = unlist(q[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(q[2]) ,xend = unlist(q[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") + 
  annotate(geom = "text",x = unlist(q[3]),y=-3,label= unlist(q[1]),fontface = 'italic',hjust = -0.4,size=1.4,color="blue") +
  
  # geom_segment(aes(x=unlist(q2[2]) ,xend = unlist(q2[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(q2[2]),xend = unlist(q2[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(q2[3]),y=-3,label=unlist(q2[1]),fontface = 'italic',hjust = -0.4,size=1.4,color="blue") +
  
  # geom_segment(aes(x=unlist(tg[2]) ,xend = unlist(tg[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(tg[2]),xend = unlist(tg[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(tg[3]),y=-3,label=unlist(tg[1]),fontface = 'italic',hjust = -0.4,size=1.4,color="blue") +
  
   # Custom the theme:
  xlab("")+ylab("XP-CLR")+
  theme(plot.margin=unit(rep(1,4),'lines'))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 9),legend.position = 'none')

p
# ggsave(p,filename = "/Users/Aoyue/Documents/test.png",width = 7,height = 7*3/16)
ggsave(p,filename = "/Users/Aoyue/Documents/test.pdf",width = 7,height = 7*3/16)


```

## LR vs. Cultivar manhattan plot for add signal dot

```{r}
################################# ##
#                               #
#  XPCLR genome line -- manhattan 100kb 50 kb 
#                               #
################################# ##
library(tidyverse)

# ### log2 snp method whole-genome SNP 3%
dataFile <- c("data/020_xpclr/ab_WE_DE_log2_0.0001_200_100000.clrgs.txt.gz")
df_wede <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/ab_DE_Durum_log2_0.0001_200_100000.clrgs.txt.gz")
df_dedm <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

dataFile <- c("data/020_xpclr/abd_log2_0.0001_200_100000.clrgs.txt.gz")
df_lrcl <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

# df <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

### ******************************************************************************
### choice: 1.df 2.color 3.threshold  
### ******************************************************************************
### 1. data
# df <- df_wede
# df <- df_dedm
df <- df_lrcl
### 2.color
# test <- "wede"
# test <- "dedm"
test <- "lrcl"
colB <- NULL
labelB <- NULL

# colB1 <- c("#ffd702","#7f5701") ### we de
# colB2 <- c("#7f5701","#016699") ### de durum
# colB3 <- c("#fc6e6e","#9900ff")  ### lr cl
## use subgenome color
colB1 <- c("#fd8582","#967bce") ### we de
colB2 <- c("#fd8582","#967bce") ### de durum
colB3 <- c("#fd8582","#967bce","#4bcdc6")  ### lr cl

### if else failed to use
if(test == "wede"){
  colB <- colB1
  labelB <- "WE vs. DE"
} 
if(test == "dedm"){
  colB <- colB2
  labelB <- "DE vs. Durum"
}
if(test == "lrcl"){
  colB <- colB3
  labelB <- "Landrace vs. Cultivar"
}

### 3. 
# top <- 0.01
top <- 0.05

# head(df)
# nrow(df)
### step2: 修改列名，使之固定
names(df)[names(df)=="XPCLR_100score"] = "Ave"
names(df)[names(df)=="CHROM"] = "Chr"
names(df)[names(df)=="ChrRef"] = "CHROM"
names(df)[names(df)=="PosRef"] = "BIN_START"
max(df$Ave)
min(df$Ave)
### step3: 去除NA和无限小的值
df <- dplyr::filter(df,!is.na(Ave) & is.finite(Ave))
nrow(df)
head(df)
max(df$Ave) # 48.50477
min(df$Ave) # 0
### step4: 获取TOP K 截距
### step4a: 排序
df_sorted <- df %>% 
  dplyr::arrange(-Ave) 
head(df_sorted)
### step4b: 获取top1% 的值，即水平截距
threshod <- df_sorted[top*nrow(df_sorted),8]
 print(paste(threshod," threshod ", top))

# write.table(d,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t") 
### step5: 建立累积坐标系统和x轴标签位置
don <- df %>%
  # Compute chromosome size
  group_by(CHROM) %>%
  summarise(chr_len=max(BIN_START)) %>%  ##每条染色体的最大值赋值给 chr_len
  # Calculate cumulative position of each chromosome
  mutate(tot=cumsum(as.numeric(chr_len))-chr_len) %>%  ##所有之和减去该行的染色体对应的长度，如：1B起始为1A长度，1D起始为1A 1B之和的长度
  select(-chr_len) %>% ##去除 chr_len这一列
  # Add this info to the initial dataset
  left_join(df, ., by=c("CHROM"="CHROM")) %>%  ## 将原数据框和don数据框进行合并
  # Add a cumulative position of each SNP
  arrange(CHROM, BIN_START) %>% ##根据CHROM BIN_START 进行排序
  mutate( BPcum=BIN_START+tot)  ##增加新的一列，改变每行点的位置

### 记录增量
dfadd <- df %>%
  # Compute chromosome size
  group_by(CHROM) %>%
  summarise(chr_len=max(BIN_START)) %>%  ##每条染色体的最大值赋值给 chr_len
  # Calculate cumulative position of each chromosome
  mutate(tot=cumsum(as.numeric(chr_len))-chr_len) %>%  ##所有之和减去该行的染色体对应的长度，如：1B起始为1A长度，1D起始为1A 1B之和的长度
  select(-chr_len)  ##去除 chr_len这一列
dfadd  
  
# don
# write.table(don,file="/Users/Aoyue/Documents/value.txt",quote=F,col.name=T,row.names=F,sep = "\t")

## x轴标签位置：求出每条染色体在坐标轴上的中间位置
axisdf = don %>% group_by(CHROM) %>% summarize(center=( max(BPcum) + min(BPcum) ) / 2 )

### 建立函数关系，添加驯化信号的位点
### 建立函数关系，添加驯化信号的位点
### 建立函数关系，添加驯化信号的位点
geneList <- read.table("data/020_xpclr/gene.txt",sep = "\t",header = T)
geneList
geneList$Gene <- as.vector(geneList$Gene)
geneList$Chromosome <- as.vector(geneList$Chromosome)
geneList$Start <- as.numeric(geneList$Start)  ## 为了可以进行预测
geneList$End <- as.numeric(geneList$End)

### 写一个函数，返回基因的标签，以及在图上的range位置
aoGetGenepos <- function(geneRow){
  gene <- geneList[geneRow,] ### get gene's dataframe, only has one line
  labelgene <- gene[1,1]### get gene's label
  chrGene <- as.vector(gene[1,2]) ### get gene's chr
  posSignal <- gene[,3:4] ### get gene's pos, a dataframe
  posOnchr <- c(posSignal[1,1],posSignal[1,2]) ### change to x 
  ### select chr, only one line
  dfadd_one <- dfadd %>% 
    filter(CHROM == chrGene) 
  # dfadd_one
  posSignalCum <- posOnchr + as.numeric(dfadd_one[1,2])
  # geneSummary <- list( labelgene,as.numeric(posSignalCum[1]),as.numeric(posSignalCum[2]))
  geneSummary <- list( labelgene,posSignalCum[1],posSignalCum[2])
  
  return(geneSummary)
}
 
#### btr1 
q <- aoGetGenepos(11)
q2 <- aoGetGenepos(12)
tg <- aoGetGenepos(13)
sig4 <- aoGetGenepos(20)
sig5 <- aoGetGenepos(21)
sig6 <- aoGetGenepos(22)
sig7 <- aoGetGenepos(23)
sig8 <- aoGetGenepos(24)


### step6: 颜色设置并画图
### 开始作图

### sample
# don <- sample_frac(don,0.1)
p <- ggplot(don, aes(x=BPcum, y=Ave)) +  ##新的横坐标，Y轴的XPCLR值不变
  # Show all points
  # geom_point( aes(color=as.factor(CHROM)), alpha=0.2, size=0.5) +  ## alpha=0.5 first test
  geom_line(aes(color=as.factor(CHROM)), alpha=0.9, size=0.2) +  ## alpha=0.5 first test
  annotate("text",x=Inf,y= threshod ,hjust= 2,vjust= -0.3,label="top 5%", size = 3)+
  geom_hline(yintercept = threshod, linetype= "dashed", color = "gray",size=0.2)+
  annotate("text",x=0,y= Inf ,hjust= -0.2,vjust= 1.1,label= labelB, size = 2)+ ## add group annotation
  scale_color_manual(values = rep(colB, 14 )) + ##合计42条染色体。故21*2=42
  # custom X axis: 每条染色体对应一个中间位置，完美！
  scale_x_continuous(expand = c(0,0), label = axisdf$CHROM, breaks= axisdf$center) +
  scale_y_continuous(expand = c(0,1),limits = c(-6,54)) +
  # scale_y_continuous(expand = c(0, 0),limits = c(-1,52)) +     # remove space between plot area and x axis
  
  ## add arrow to represent the selection signal 
  # geom_vline(xintercept = 5838260571,color="black",size=0.2)+
  #geom_segment(aes(x=unlist(q[2]) ,xend = unlist(q[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(q[2]) ,xend = unlist(q[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") + 
  annotate(geom = "text",x = unlist(q[3]),y=-3,label= unlist(q[1]),fontface = 'italic',hjust = -0.1,size=1.4,color="blue") +
  
  #geom_segment(aes(x=unlist(q2[2]) ,xend = unlist(q2[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(q2[2]),xend = unlist(q2[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(q2[3]),y=-3,label=unlist(q2[1]),fontface = 'italic',hjust = -0.1,size=1.4,color="blue") +
  
  #geom_segment(aes(x=unlist(tg[2]) ,xend = unlist(tg[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(tg[2]),xend = unlist(tg[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(tg[3]),y=-3,label=unlist(tg[1]),fontface = 'italic',hjust = -0.1,size=1.4,color="blue") +
  
  # geom_segment(aes(x=unlist(sig4[2]) ,xend = unlist(sig4[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(sig4[2]),xend = unlist(sig4[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(sig4[3]),y=-3,label=unlist(sig4[1]),fontface = 'italic',hjust = -0.1,size=1.4,color="blue") +
  
  # geom_segment(aes(x=unlist(sig5[2]) ,xend = unlist(sig5[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(sig5[2]),xend = unlist(sig5[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(sig5[3]),y=-3,label=unlist(sig5[1]),fontface = 'italic',hjust = -0.1,size=1.4,color="blue") +
  
  # geom_segment(aes(x=unlist(sig6[2]) ,xend = unlist(sig6[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(sig6[2]),xend = unlist(sig6[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(sig6[3]),y=-3,label=unlist(sig6[1]),fontface = 'italic',hjust = 1.1,size=1.4,color="blue") +
  
  # geom_segment(aes(x=unlist(sig7[2]) ,xend = unlist(sig7[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(sig7[2]),xend = unlist(sig7[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(sig7[3]),y=-3,label=unlist(sig7[1]),fontface = 'italic',hjust = -0.1,size=1.4,color="blue") +
  
  # geom_segment(aes(x=unlist(sig8[2]) ,xend = unlist(sig8[2]),y=0,yend = Inf) , size=0.2, color="blue") +
  annotate(geom = "segment",x=unlist(sig8[2]),xend = unlist(sig8[3]),y=-5,yend = 0, arrow=arrow(length = unit(0.08, "cm"), type = "closed"), size=0.1, color="blue") +
  annotate(geom = "text",x = unlist(sig8[3]),y=-3,label=unlist(sig8[1]),fontface = 'italic',hjust = -0.1,size=1.4,color="blue") +
  
   # Custom the theme:
  xlab("")+ylab("XP-CLR")+
  theme(plot.margin=unit(rep(1,4),'lines'))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 9),legend.position = 'none')

p
# ggsave(p,filename = "/Users/Aoyue/Documents/test.png",width = 7,height = 7*3/16)

ggsave(p,filename = "/Users/Aoyue/Documents/test.pdf",width = 7,height = 7*3/16)

```

### Top 5%
```{r}
# ### log2 snp method whole-genome SNP 3%


dataFile <- c("data/020_xpclr/ab_WE_DE_log2_0.0001_200_100000.clrgs.txt.gz")
df_wede <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

### 找到 top5 的值域是多少，并且输出bins数目  
df <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)
dftest <- df %>% slice_max(XPCLR_100score,prop = 0.05)
dffinalline <- dftest %>% slice(n())
regions <- nrow(dftest)
print(paste("Threshod is",round(dffinalline[1,8],3) ,"Regions total",regions))
```

```{r}
dataFile <- c("data/020_xpclr/ab_DE_Durum_log2_0.0001_200_100000.clrgs.txt.gz")
df_dedm <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

### 找到 top5 的值域是多少，并且输出bins数目  
df <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)
dftest <- df %>% slice_max(XPCLR_100score,prop = 0.05)
dffinalline <- dftest %>% slice(n())
regions <- nrow(dftest)
print(paste("Threshod is",round(dffinalline[1,8],3) ,"Regions total",regions))
```

```{r}
dataFile <- c("data/020_xpclr/abd_log2_0.0001_200_100000.clrgs.txt.gz")
df_lrcl <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)

### 找到 top5 的值域是多少，并且输出bins数目  
df <- read.table(dataFile, header=TRUE, sep="\t",row.names=NULL)
dftest <- df %>% slice_max(XPCLR_100score,prop = 0.05)
dffinalline <- dftest %>% slice(n())
regions <- nrow(dftest)
print(paste("Threshod is",round(dffinalline[1,8],3) ,"Regions total",regions))
```


### Top 5% 和 Top 1% 求值域;求受选择区域的CDS长度及占比
- a:ind top K and its threshod
- b:Get slective regions num
- c:Get proportion of whole-genome region gene region and cds region
- d:Get gene list and GO analysis

Table:
Group   TopK    Threshod    SelectiveSweep    GeneNum   Len_wholegenome   Len_Gene    Len_CDS   Prop_whole    Prop_Gene    Prop_CDS

#### Get nonoverlap trans
```{r}
# dataFile <- c("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/001_geneTable/wheat_v1.1_nonoverlap_addPos.txt.gz")
# df <- read.delim(dataFile) %>% 
#   filter(Is_Unique_gene.1.0. == 1) %>% 
#   select(Longest_transcript)
# 
# write_tsv(df,file = "/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/006_nonoverlapGene/nonoverlapGene.txt")
```
#### 函数构建及运行
```{r}
AoGetXPCLRsummary <- function(Group,TopK,bin,dataFile){
  library(bedtoolsr)
  # bin <- 100000
  ### 注意不加 chl mt and chrUn
  genomeSize_aabbdd <- 14066280851
  ### genomeSize_aabb <- 9652829943
  ### genomeSize_dd <- 3951074735 
  ### 注意不加 chr0
  geneSize <- 370393700
  cdsSize <- 130348274
  
  ### @para1 Group
  # Group <- "wede"
  ### @para2 TopK
  # TopK <- top1
  
  # ### log2 snp method whole-genome SNP 3%
  ### 找到 topK 的值域是多少，并且输出bins数目
  df <- read.delim(dataFile) %>% slice_max(XPCLR_100score, prop = TopK)
  dffinalline <- df %>% slice(n()) ### 最后一行的内容,为了找到值域的确切值
  ### @para3 Threshod
  Threshod <- round(dffinalline[1, 8], 3)
  ### @para4 SelectiveSweep
  SelectiveSweep <- nrow(df)
  
  print(paste("Threshod is", Threshod, "Regions total", SelectiveSweep))
  
  ### ============================================================
  ### step1:make bed
  dfquery <- df %>% select(CHROM, POS) %>% mutate(POS_end = POS + bin)
  dfgff3gene <-
    read.delim(
      "/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/007_bed_onWholeGenome/GFF3_gene_nonoverlap.bed"
    )
  dfgff3cds <-
    read.delim(
      "/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/007_bed_onWholeGenome/GFF3_CDS_nonoverlap.bed"
    )
  ### step2:calculate the length on whole genome, gene region, and cds region
  ### a::genomesize AABBDD : 14066280851 AABB : 9652829943  DD : 3951074735  do not add chl and mt
  dfwhole <- dfquery %>% mutate(Length = POS_end - POS) %>%
    summarise(sum_whole = sum(Length)) %>%
    mutate(prop_whole = sum_whole / genomeSize_aabbdd)
  
  dataFile <-
    c(
      "/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/000_chrDB/chrConvertionRule_byAoyue.txt"
    )
  dfchrDB <- read.delim(dataFile) %>% select(CHROM = ChrID, Subgenome)
  
  dfsub <- dfquery %>% mutate(Length = POS_end - POS) %>%
    left_join(.,dfchrDB,by="CHROM") %>% 
    group_by(Subgenome) %>% 
    summarise(sum_whole = sum(Length)) %>%
    mutate(prop_whole = sum_whole / genomeSize_aabbdd) %>% 
    mutate(Group = Group, TopK = TopK)
  
  ### b::目的：保留选择区和 gene 的交集，并且得到交集所在的基因
  dfintersect_gene <-
    bt.intersect(dfgff3gene, dfquery) %>% arrange(V1, V2) %>% mutate(Subtract = V3 - V2)
  ### c::目的：保留选择区和 cds 的交集，并且得到交集所在的基因
  dfintersect_cds <-
    bt.intersect(dfgff3cds, dfquery) %>% arrange(V1, V2) %>% mutate(Subtract = V3 - V2)
  
  ### @para5 GeneNum
  dfgeneSelected <- dfintersect_gene %>% group_by(V4) %>% count() %>% distinct(., V4) %>% mutate(Group = Group, TopK = TopK)

  GeneNum <- dfintersect_gene %>% group_by(V4) %>% count() %>% distinct(., V4) %>% nrow()
  ### @para6 Len_wholegenome
  Len_wholegenome <- dfwhole$sum_whole
  ### @para7 Len_Gene
  Len_Gene <- sum(dfintersect_gene$Subtract)
  ### @para8 Len_CDS
  Len_CDS <- sum(dfintersect_cds$Subtract)
  ### @para9 Prop_whole
  Prop_whole <- dfwhole$prop_whole
  ### @para10 Prop_Gene
  Prop_Gene <- Len_Gene / geneSize
  ### @para11 Prop_CDS
  Prop_CDS <- Len_CDS / cdsSize
  
  dffinal <-
    data.frame(
      Group = Group,
      TopK = TopK,
      Threshod = Threshod,
      SelectiveSweep = SelectiveSweep,
      GeneNum = GeneNum,
      Len_wholegenome = Len_wholegenome,
      Len_Gene = Len_Gene,
      Len_CDS = Len_CDS,
      Prop_whole = Prop_whole,
      Prop_Gene = Prop_Gene,
      Prop_CDS = Prop_CDS
    )
  
  dfout <- list(dffinal,dfsub,dfgeneSelected)  ### 测试失败
  # return(dffinal)
  # return(dfgeneSelected)
  # return(dfsub)
  return(dfout)
}


dataFile1 <-c("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/004_output/ab_WE_DE_log2_0.0001_600_100000.clrgs.txt.gz")
dataFile2 <- c("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/004_output/ab_DE_Durum_log2_0.0001_600_100000.clrgs.txt.gz")
dataFile3 <- c("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/004_output/abd_log2_0.0001_600_100000.clrgs.txt.gz")
 
Group <- c("wede","wede","dedurum","dedurum","lrcl","lrcl")
TopK <- c(0.05,0.01,0.05,0.01,0.05,0.01)
bin <- rep(100000,6)
dataFile <- c(dataFile1,dataFile1,dataFile2,dataFile2,dataFile3,dataFile3) 
dfpara <- data.frame(Group=Group,TopK=TopK,bin=bin,dataFile=dataFile,stringsAsFactors = F)

dfm <- dfpara %>%
  # group_modify(~{.x %>% AoGetXPCLRsummary()})  %>% 
  rowwise() %>% 
  summarise(sum = list(AoGetXPCLRsummary(Group,TopK,bin,dataFile))) %>% 
  unnest()

### output   
# filename <- c("~/Documents/test.txt")
# write_tsv(dfm,file = filename)
```


#### 单个文件测试
```{r}

Group <- dfpara$Group[1]
TopK <- dfpara$TopK[1]
bin <- dfpara$bin[1]
dataFile <- dfpara$dataFile[1]

library(bedtoolsr)
  # bin <- 100000
  ### 注意不加 chl mt and chrUn
  genomeSize_aabbdd <- 14066280851
  ### genomeSize_aabb <- 9652829943
  ### genomeSize_dd <- 3951074735 
  ### 注意不加 chr0
  geneSize <- 370393700
  cdsSize <- 130348274
  
  ### @para1 Group
  # Group <- "wede"
  ### @para2 TopK
  # TopK <- top1
  
  # ### log2 snp method whole-genome SNP 3%
  ### 找到 topK 的值域是多少，并且输出bins数目
  df <- read.delim(dataFile) %>% slice_max(XPCLR_100score, prop = TopK)
  dffinalline <- df %>% slice(n()) ### 最后一行的内容,为了找到值域的确切值
  ### @para3 Threshod
  Threshod <- round(dffinalline[1, 8], 3)
  ### @para4 SelectiveSweep
  SelectiveSweep <- nrow(df)
  
  print(paste("Threshod is", Threshod, "Regions total", SelectiveSweep))
  
  ### ============================================================
  ### step1:make bed
  dfquery <- df %>% select(CHROM, POS) %>% mutate(POS_end = POS + bin)
  dfgff3gene <-
    read.delim(
      "/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/007_bed_onWholeGenome/GFF3_gene_nonoverlap.bed"
    )
  dfgff3cds <-
    read.delim(
      "/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/009_xpclr/007_bed_onWholeGenome/GFF3_CDS_nonoverlap.bed"
    )
  ### step2:calculate the length on whole genome, gene region, and cds region
  ### a::genomesize AABBDD : 14066280851 AABB : 9652829943  DD : 3951074735  do not add chl and mt
  dfwhole <- dfquery %>% mutate(Length = POS_end - POS) %>%
    summarise(sum_whole = sum(Length)) %>%
    mutate(prop_whole = sum_whole / genomeSize_aabbdd)
  
  dataFile <-
    c(
      "/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/000_chrDB/chrConvertionRule_byAoyue.txt"
    )
  dfchrDB <- read.delim(dataFile) %>% select(CHROM = ChrID, Subgenome)
  
  dfsub <- dfquery %>% mutate(Length = POS_end - POS) %>%
    left_join(.,dfchrDB,by="CHROM") %>% 
    group_by(Subgenome) %>% 
    summarise(sum_whole = sum(Length)) %>%
    mutate(prop_whole = sum_whole / genomeSize_aabbdd)
  
  ### b::目的：保留选择区和 gene 的交集，并且得到交集所在的基因
  dfintersect_gene <-
    bt.intersect(dfgff3gene, dfquery) %>% arrange(V1, V2) %>% mutate(Subtract = V3 - V2)
  
  dfintersect_gene_bySub <- dfintersect_gene %>% 
    mutate()
  
  ### c::目的：保留选择区和 cds 的交集，并且得到交集所在的基因
  dfintersect_cds <-
    bt.intersect(dfgff3cds, dfquery) %>% arrange(V1, V2) %>% mutate(Subtract = V3 - V2)
  
  ### @para5 GeneNum
  dfgeneSelected <- dfintersect_gene %>% group_by(V4) %>% count() %>% distinct(., V4) %>% mutate(Group = Group, TopK = TopK)

  GeneNum <- dfintersect_gene %>% group_by(V4) %>% count() %>% distinct(., V4) %>% nrow()
  ### @para6 Len_wholegenome
  Len_wholegenome <- dfwhole$sum_whole
  ### @para7 Len_Gene
  Len_Gene <- sum(dfintersect_gene$Subtract)
  ### @para8 Len_CDS
  Len_CDS <- sum(dfintersect_cds$Subtract)
  ### @para9 Prop_whole
  Prop_whole <- dfwhole$prop_whole
  ### @para10 Prop_Gene
  Prop_Gene <- Len_Gene / geneSize
  ### @para11 Prop_CDS
  Prop_CDS <- Len_CDS / cdsSize
  
  dffinal <-
    data.frame(
      Group = Group,
      TopK = TopK,
      Threshod = Threshod,
      SelectiveSweep = SelectiveSweep,
      GeneNum = GeneNum,
      Len_wholegenome = Len_wholegenome,
      Len_Gene = Len_Gene,
      Len_CDS = Len_CDS,
      Prop_whole = Prop_whole,
      Prop_Gene = Prop_Gene,
      Prop_CDS = Prop_CDS
    )
  
  
```


####============
## 不同年代 cultivar 的变化
### 0.cultivar 年代分布
```{r}
df2 <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>%
  filter(IfOnVmap2_S643 != 1) %>% 
  filter(GroupbyContinent=="CL") %>% 
  filter(Group == "CAAS") %>% 
  filter(!is.na(ReleasedIntroducedYear)) %>% 
  mutate(ReleasedIntroducedYear = as.numeric(ReleasedIntroducedYear)) %>%   
  mutate(bin=cut(ReleasedIntroducedYear,breaks=c(1920,1950,1980,2020))) %>% 
  mutate(GroupTemp =  str_split_fixed(bin,",",n=2)[,1]) %>% 
  mutate(GroupTemp =as.numeric(str_sub(GroupTemp,2))) %>% 
  select(-bin) 

### 用年代和分组IMC MCC 画个图，看看国外引进来的在哪个年代
p <- df2 %>%
  group_by(ReleasedIntroducedYear,GroupByOriPublication_Abbre,GroupTemp) %>%
  count() %>% 
  ggplot(aes(x=as.character(ReleasedIntroducedYear),y=n,fill=GroupByOriPublication_Abbre))+
  geom_bar(stat="identity",alpha = 1)+  #,position=position_dodge()
  geom_text(aes(label=  paste(prettyNum(n,big.mark = ","))),position =  position_stack(),vjust=1.5, size=3,color="black")+
  labs(y="Count",x="")+
  guides(fill=guide_legend(title="Group"))+
  theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  theme(panel.background = element_rect(size = 0.7, colour = 'black',fill = 'white'),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),text = element_text(size = 10),
        legend.position = c(0.1,0.6))

p
ggsave(p,filename = "~/Documents/bar_YearAsCharacter.pdf",height = 3,width = 7)

### 用年代和分组IMC MCC 画个图，看看国外引进来的在哪个年代
p <- df2 %>%
  group_by(ReleasedIntroducedYear,GroupTemp) %>%
  count() %>% 
  ggplot(aes(x=ReleasedIntroducedYear,y=n,fill=ReleasedIntroducedYear))+
  geom_bar(stat="identity",alpha = 1)+  #,position=position_dodge()
  geom_text(aes(label=  paste(prettyNum(n,big.mark = ","))),position =  position_stack(),vjust=1.5, size=3,color="black") +
  labs(y="Count",x="")+
  scale_fill_distiller(type = "seq", palette = "Spectral", name="Introduced year",direction = 1)+
  # theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  theme(panel.background = element_rect(size = 0.7, colour = 'black',fill = 'white'),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),text = element_text(size = 10),
        legend.position = c(0.1,0.6))
        

p
ggsave(p,filename = "~/Documents/bar_YearAsNumeric.pdf",height = 3,width = 7)

### 用年代和分组IMC MCC 画个图，看看国外引进来的在哪个年代
p <- df2 %>%
  group_by(GroupTemp) %>% 
  count() %>% 
  ggplot(aes(x=as.character(GroupTemp),y=n,fill=GroupTemp))+
  geom_bar(stat="identity",alpha = 1)+  #,position=position_dodge()
  geom_text(aes(label=  paste(prettyNum(n,big.mark = ","))),position =  position_stack(),vjust=1.5, size=5,color="orange") +
  labs(y="Count",x="")+
  scale_fill_distiller(type = "seq", palette = "BuGn", name="Introduced year",direction = 1)+
  # theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  theme(panel.background = element_rect(size = 0.7, colour = 'black',fill = 'white'),panel.grid.major = element_blank(), panel.grid.minor = element_blank(),text = element_text(size = 13),
        legend.position = c(0.2,0.6))
        

p
ggsave(p,filename = "~/Documents/bar_bin.pdf",height = 3,width = 5)
```

### 1.划分 Bin，del Ratio: cultivar 和 landrace load ratio 变化，区分亚基因组
```{r}
dfdelRatio <- read_tsv("data/004_delCount/delRatio.txt.gz")
dftaxaDB <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>% select(Taxa=VcfID,GroupbyContinent,ReleasedIntroducedYear,Group) %>% 
  mutate(ReleasedIntroducedYear = str_replace(ReleasedIntroducedYear,"s",""))

### 有年代信息的cultivar
df1 <- dfdelRatio %>% 
  left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(GroupbyContinent == "CL") %>% 
  filter(!is.na(ReleasedIntroducedYear)) %>% 
  filter(DelCountGroup=="TotalDerivedSNPCount") %>% 
  ### 划分 Bin
  mutate(ReleasedIntroducedYear = as.numeric(ReleasedIntroducedYear)) %>%   
  # mutate(bin=cut_width(ReleasedIntroducedYear,width = 30,boundary = 0,closed = "left")) %>%
  mutate(bin=cut(ReleasedIntroducedYear,breaks=c(1920,1960,1980,2000,2020))) 
  # mutate(binString = str_split_fixed(bin,",",n=2)[,1]) %>% 
  # mutate(binString =as.numeric(str_sub(binString,2)))
  # arrange(.,by=binString) 

### LR_EU 及无年代信息的 cultivar
df2 <- dfdelRatio %>% 
  left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(Group != "CAAS") %>% 
  filter(GroupbyContinent %in%  c("CL","LR_EU")) %>% 
  filter(DelCountGroup=="TotalDerivedSNPCount") %>% 
  ### 划分 Bin
  mutate(bin=GroupbyContinent) 

### 设置因子顺序及标签
mygroup <- c("Ratio_002_nonsynonymous","Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5","Ratio_007_VEP")
mygrouplabelA <- c("Nonsynonymous","Deleterious","High impact")

mygroup2 <- c("LR_EU","(1.92e+03,1.96e+03]","(1.96e+03,1.98e+03]","(1.98e+03,2e+03]","(2e+03,2.02e+03]","CL")
mygroup2labelA <- c("LR","1920s-1960s","1960s-1980s","1980s-2000s","2000s-2020s","CL_timeUnknown")

df <- rbind(df1,df2) %>% 
  filter(Stand_LoadGroup %in% mygroup ) %>% 
  droplevels() %>% 
  mutate(Stand_LoadGroup=factor(Stand_LoadGroup,levels = mygroup,labels = mygrouplabelA)) %>% 
  mutate(bin=factor(bin,levels = mygroup2,labels = mygroup2labelA))

### plot
library(RColorBrewer)
# colB <- colorRampPalette(brewer.pal(9,"Reds")[c(3,7,9)])(5)
colB <- c("#fc6e6e",brewer.pal(9,"Blues")[c(4,6,7,8)],"#9900ff") 
  
p <- df %>% 
  ggplot(aes(x=Subgenome, y= Stand_delCount)) +
  #### pointrange
  stat_summary(aes(color=bin),fun.data=mean_sdl,
               fun.args = list(mult = 1),
               position=position_dodge(1),geom="pointrange")+
  stat_summary(aes(group = bin), fun=mean, geom="point", color="white",position = position_dodge(1))+
  scale_fill_manual(values = colB,name='')+  
  scale_color_manual(values = colB,name='')+
  # scale_x_discrete(labels=c("1920s-1950s","1950s-1980s","1980s-2020s"))+
  labs(x="",y="Mutation load")+
  facet_wrap(~Stand_LoadGroup,scales = "free")+
  # theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  theme_classic()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(size=0.4, colour = "black"),
        # legend.position = 'none',
        text = element_text(size = 10))
p
ggsave(p,filename = "~/Documents/test.pdf",height = 2.5,width = 7)
```

## 2.划分 Bin，del Count: cultivar 和 landrace load ratio 变化，区分亚基因组
```{r}
dfdel <- read_tsv("data/004_delCount/delCount.txt.gz")
dftaxaDB <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>% select(Taxa=VcfID,GroupbyContinent,ReleasedIntroducedYear,Group) %>% 
  mutate(ReleasedIntroducedYear = str_replace(ReleasedIntroducedYear,"s",""))

### 有年代信息的cultivar
df1 <- dfdel %>% 
  left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(GroupbyContinent == "CL") %>% 
  filter(!is.na(ReleasedIntroducedYear)) %>% 
  ### 划分 Bin
  mutate(ReleasedIntroducedYear = as.numeric(ReleasedIntroducedYear)) %>%   
  # mutate(bin=cut_width(ReleasedIntroducedYear,width = 30,boundary = 0,closed = "left")) %>%
  mutate(bin=cut(ReleasedIntroducedYear,breaks=c(1920,1960,1980,2000,2020))) 
  # mutate(binString = str_split_fixed(bin,",",n=2)[,1]) %>% 
  # mutate(binString =as.numeric(str_sub(binString,2)))
  # arrange(.,by=binString) 

### LR_EU 及无年代信息的 cultivar
df2 <- dfdel %>% 
  left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(Group != "CAAS") %>% 
  filter(GroupbyContinent %in%  c("CL","LR_EU")) %>% 
  ### 划分 Bin
  mutate(bin=GroupbyContinent) 

### 设置因子顺序及标签
mygroup <- c("002_nonsynonymous","031_GERP16way_1.5andDerived_PolyPhen2_0.5","007_VEP")
mygrouplabelA <- c("Nonsynonymous","Deleterious","High impact")

mygroup2 <- c("LR_EU","(1.92e+03,1.96e+03]","(1.96e+03,1.98e+03]","(1.98e+03,2e+03]","(2e+03,2.02e+03]","CL")
mygroup2labelA <- c("LR","1920s-1960s","1960s-1980s","1980s-2000s","2000s-2020s","CL_timeUnknown")

df <- rbind(df1,df2) %>% 
  filter(VariantsGroup %in% mygroup) %>% 
  droplevels() %>% 
  mutate(VariantsGroup=factor(VariantsGroup,levels = mygroup,labels = mygrouplabelA)) %>% 
  mutate(bin=factor(bin,levels = mygroup2,labels = mygroup2labelA))

### plot
library(RColorBrewer)
# colB <- colorRampPalette(brewer.pal(9,"Reds")[c(3,7,9)])(5)
colB <- c("#fc6e6e",brewer.pal(9,"Blues")[c(4,6,7,8)],"#9900ff") 
  
p <- df %>% 
  ggplot(aes(x=Subgenome, y= TotalDerivedSNPCount)) +
  #### pointrange
  stat_summary(aes(color=bin),fun.data=mean_sdl,
               fun.args = list(mult = 1),
               position=position_dodge(1),geom="pointrange")+
  stat_summary(aes(group = bin), fun=mean, geom="point", color="white",position = position_dodge(1))+
  scale_fill_manual(values = colB,name='')+  
  scale_color_manual(values = colB,name='')+
  # scale_x_discrete(labels=c("1920s-1950s","1950s-1980s","1980s-2020s"))+
  labs(x="",y="Derived allele count per individual")+
  facet_wrap(~VariantsGroup,scales = "free")+
  # theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  theme_classic()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(size=0.4, colour = "black"),
        # legend.position = 'none',
        text = element_text(size = 10))
p
ggsave(p,filename = "~/Documents/test.pdf",height = 2.5,width = 7)
```


### 求平均值
colnames(df)
 [1] "Taxa"                   "Subgenome"              "VariantsGroup"          "IfSelected"            
 [5] "Group_refobj"           "TotalDerivedSNPCount"   "HomoDerivedSNPCount"    "HeterDerivedSNPCount"  
 [9] "SiteCountWithMinDepth"  "GroupbyContinent"       "ReleasedIntroducedYear" "Group"       
[13] "bin" 
```{r}
df_ave <- df %>% 
  filter(bin !="LR") %>% 
  group_by(VariantsGroup,Subgenome) %>% 
  summarise(mean = mean(TotalDerivedSNPCount))

df_ave2 <- df_ave %>% 
  group_by(VariantsGroup) %>% 
  summarise(sum = sum(mean))

```






## 3.whole genome & selection region - ratio
```{r}
############################## choice

### 建立画图函数，然后调用画图
fun_plotDelLoad_Boxplot <- function(para1,Prefix){

### 外部定义变量后，内部函数可以调用，也可以重新命名。因此，choiceA1_top5，choiceA1_top1，choiceA3_wholegenome这3个变量可以不用在这里再次写出
# choiceA1_top5 <- c("data/023_del_XPCLR/top5/delCount.txt.gz")
# choiceA2_top1 <- c("data/023_del_XPCLR/top1/delCount.txt.gz")
# choiceA3_wholegenome <- c("data/004_delCount/delCount.txt.gz")

### 导入种质数据库,并选取特定的列
dftaxaDB <- read.delim("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>% select(Taxa=VcfID,GroupbyContinent,ReleasedIntroducedYear) %>% 
  mutate(ReleasedIntroducedYear = str_replace(ReleasedIntroducedYear,"s",""))

### topK
if(para1 %in% c(choiceA1_top5,choiceA2_top1)){
  dfdelRatio <- read.delim(para1) %>% 
  left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(IfSelected==1) %>% select(-IfSelected)
}else if(para1 %in% c(choiceA3_wholegenome)){
  dfdelRatio <- read.delim(para1) %>% 
  left_join(.,dftaxaDB,by="Taxa")
}
 
mygroup <- c("Ratio_002_nonsynonymous","Ratio_031_GERP16way_1.5andDerived_PolyPhen2_0.5","Ratio_007_VEP")
mygrouplabelA <- c("Nonsynonymous","Deleterious","High impact")

df <- dfdelRatio %>% 
  filter(GroupbyContinent == "CL") %>% select(-GroupbyContinent) %>%
  filter(!is.na(ReleasedIntroducedYear)) %>% 
  filter(DelCountGroup=="TotalDerivedSNPCount") %>%  ###选择该种类型的ratio值
  #### 可以选择哪一列作为有害突变的计数 ###
  filter(Stand_LoadGroup %in% mygroup) %>% 
  mutate(Stand_LoadGroup=factor(Stand_LoadGroup,levels = mygroup,labels = mygrouplabelA)) %>% 
  ### 不区分亚基因组
  group_by(Taxa,ReleasedIntroducedYear,Stand_LoadGroup) %>% 
  summarise(across(.cols = c(Stand_delCount), .fns = mean)) %>% 
  ### 划分 Bin
  mutate(ReleasedIntroducedYear = as.numeric(ReleasedIntroducedYear)) %>%   
  mutate(bin=cut(ReleasedIntroducedYear,breaks=c(1920,1960,1980,2000,2020))) %>%
  filter(!is.na(bin))

### plot
library(RColorBrewer)
colB <- colorRampPalette(brewer.pal(9,"Reds")[c(3,5,7,9)])(5)

p <- df %>% 
  ggplot(aes(x=bin, y= Stand_delCount)) +
  stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
  geom_boxplot(aes(fill=bin),alpha = 0.2,outlier.color = NA) +
  geom_point(position = position_jitterdodge(0.6),alpha = 0.8,aes(color=bin),size=2)+
  stat_summary(fun=mean, geom="point", color="yellow",position = position_dodge(0.8),size=2)+ ### 画出平均值
  scale_fill_manual(values = colB,name='')+
  scale_color_manual(values = colB,name='')+  
  scale_x_discrete(labels=c("1920s-1960s","1960s-1980s","1980s-2000s","2000s-2020s"))+
  labs(x="",y="Mutation load per accession")+
  facet_wrap(~Stand_LoadGroup,scales = "free",ncol = 1)+
  theme(strip.background = element_blank())+
  theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 30))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7),
    legend.position = 'none',
    text = element_text(size = 10))
p
# ggsave(p,filename = "~/Documents/test.pdf",height = 6,width = 4)

ggsave(p,filename = str_c("~/Documents/",Prefix,"_delLoad_changeByYears.pdf"),height = 6,width = 4)
return(p)
}


######  开始调用函数
choiceA1_top5 <- c("data/023_del_XPCLR/delRatio_top0.05.txt.gz")
choiceA2_top1 <- c("data/023_del_XPCLR/delRatio_top0.01.txt.gz")
choiceA3_wholegenome <- c("data/004_delCount/delRatio.txt.gz")
df1 <- fun_plotDelLoad_Boxplot(choiceA3_wholegenome,"whole_genome")
df2 <- fun_plotDelLoad_Boxplot(choiceA1_top5,"top5")
df3 <- fun_plotDelLoad_Boxplot(choiceA2_top1,"top1")
df1
df2
df3

###### 合并所画的图形

### M2
q <- ggpubr::ggarrange(plotlist=list(df1,df2,df3),
                       labels = c("Whole-genome","Top 5% of selection regions","Top 1% of selection regions"),
                       font.label = list(size = 7),
                       ncol = 3,nrow = 1)
cowplot::save_plot(q,filename = "~/Documents/test.pdf",base_height = 7.2,base_width = 7.2)

```


## 4.whole genome & selection region - count
```{r}
############################## choice

### 建立画图函数，然后调用画图
fun_plotDelLoad_Boxplot <- function(para1,Prefix){

### 外部定义变量后，内部函数可以调用，也可以重新命名。因此，choiceA1_top5，choiceA1_top1，choiceA3_wholegenome这3个变量可以不用在这里再次写出
# choiceA1_top5 <- c("data/023_del_XPCLR/top5/delCount.txt.gz")
# choiceA2_top1 <- c("data/023_del_XPCLR/top1/delCount.txt.gz")
# choiceA3_wholegenome <- c("data/004_delCount/delCount.txt.gz")

### 导入种质数据库,并选取特定的列
dftaxaDB <- read.delim("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt") %>% select(Taxa=VcfID,GroupbyContinent,IfOnVmap2_S643,ReleasedIntroducedYear) %>% 
  mutate(ReleasedIntroducedYear = str_replace(ReleasedIntroducedYear,"s",""))

### topK
if(para1 %in% c(choiceA1_top5,choiceA2_top1)){
  dfdelCount <- read.delim(para1) %>% 
  left_join(.,dftaxaDB,by="Taxa") %>% 
  filter(IfSelected==1) %>% select(-IfSelected)
}else if(para1 %in% c(choiceA3_wholegenome)){
  dfdelCount <- read.delim(para1) %>% 
  left_join(.,dftaxaDB,by="Taxa")
}
 
mygroup <- c("002_nonsynonymous","031_GERP16way_1.5andDerived_PolyPhen2_0.5","007_VEP")
mygrouplabelA <- c("Nonsynonymous","Deleterious","High impact")

df <- dfdelCount %>% 
  filter(GroupbyContinent == "CL") %>% 
  filter(!is.na(ReleasedIntroducedYear)) %>% 
  #### 可以选择哪一列作为有害突变的计数 ###
  filter(VariantsGroup %in% mygroup) %>% 
  mutate(VariantsGroup=factor(VariantsGroup,levels = mygroup,labels = mygrouplabelA)) %>% 
  ### 不区分亚基因组
  group_by(Taxa,ReleasedIntroducedYear,VariantsGroup) %>% 
  summarise(across(.cols = c(TotalDerivedSNPCount:SiteCountWithMinDepth), .fns = sum)) %>% 
  ### 划分 Bin
  mutate(ReleasedIntroducedYear = as.numeric(ReleasedIntroducedYear)) %>%   
  mutate(bin=cut(ReleasedIntroducedYear,breaks=c(1920,1960,1980,2000,2020))) %>%
  filter(!is.na(bin))

### plot
library(RColorBrewer)
colB <- colorRampPalette(brewer.pal(9,"Reds")[c(3,5,7,9)])(5)

p <- df %>% 
  ggplot(aes(x=bin, y= TotalDerivedSNPCount)) +
  stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
  geom_boxplot(aes(fill=bin),alpha = 0.2,outlier.color = NA) +
  geom_point(position = position_jitterdodge(0.6),alpha = 0.8,aes(color=bin),size=2)+
  stat_summary(fun=mean, geom="point", color="yellow",position = position_dodge(0.8),size=2)+ ### 画出平均值
  scale_fill_manual(values = colB,name='')+
  scale_color_manual(values = colB,name='')+  
  scale_x_discrete(labels=c("1920s-1960s","1960s-1980s","1980s-2000s","2000s-2020s"))+
  labs(x="",y="No. derived SNPs per accession")+
  facet_wrap(~VariantsGroup,scales = "free",ncol = 1)+
  theme(strip.background = element_blank())+
  theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 30))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7),
    legend.position = 'none',
    text = element_text(size = 10))
p
# ggsave(p,filename = "~/Documents/test.pdf",height = 6,width = 4)

ggsave(p,filename = str_c("~/Documents/",Prefix,"_delLoad_changeByYears.pdf"),height = 6,width = 4)
return(p)
}


######  开始调用函数
choiceA1_top5 <- c("data/023_del_XPCLR/delCount_top0.05.txt.gz")
choiceA2_top1 <- c("data/023_del_XPCLR/delCount_top0.01.txt.gz")
choiceA3_wholegenome <- c("data/004_delCount/delCount.txt.gz")
df1 <- fun_plotDelLoad_Boxplot(choiceA3_wholegenome,"whole_genome")
df2 <- fun_plotDelLoad_Boxplot(choiceA1_top5,"top5")
df3 <- fun_plotDelLoad_Boxplot(choiceA2_top1,"top1")
df1
df2
df3

###### 合并所画的图形

### M2
q <- ggpubr::ggarrange(plotlist=list(df1,df2,df3),
                       labels = c("Whole-genome","Top 5% of selection regions","Top 1% of selection regions"),
                       font.label = list(size = 7),
                       ncol = 3,nrow = 1)
cowplot::save_plot(q,filename = "~/Documents/test.pdf",base_height = 7.2,base_width = 7.2)
```