Fig02.Rmd

---
title: "Fig 2 Generating a high-confidence set of deleterious variants"
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
---

# Fig2

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

## dftotal and dftotal_withAncestral

```{r}
dftotal <- read_tsv("data/005_SNPAnnoDB/geneSNPAnno.txt.gz") %>%
  mutate(Sub=str_sub(Transcript, 9, 9)) 

  # filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>%
  # mutate(IfRefisAnc = if_else(Ref == Ancestral,"Anc","Der")) %>% 
  # mutate(Sub=str_sub(Transcript, 9, 9)) %>% 

dftotal_ancestral <- dftotal %>% 
  filter(!is.na(Ancestral),Ancestral == Ref | Ancestral == Alt) 

### 现将基因组过滤，只保留有 ancestral 状态的，并且添加LoadGroup、IfRefisAnc、CommonOrRare 等分组信息
dftotal_withAncestral <- dftotal %>% 
  filter(Region=="CDS",Variant_type=="SYNONYMOUS"|Variant_type=="NONSYNONYMOUS")%>%
  filter(!is.na(Ancestral),Ancestral == Ref | Ancestral == Alt) %>%
  mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  mutate(IfRefisAnc = if_else(Ref == Ancestral,"Anc","Der")) %>% 
  mutate(CommonOrRare = if_else(Maf <= 0.05, "Rare","Common")) %>% 
  filter(!is.na(LoadGroup))
### 特别注意，这里根据MAF进行分组，按照维基百科的定义
### 2个filter 5个mutate 1个filter

```

#### =====================

## VariantsType count all region

```{r}
### 需要手动选择
factor1 <- "Region"
factor2 <- "Variant_type"
factor3 <- "Effect_VEP"
factor4 <- "Effect_snpEff"

df_variantsType1 <- dftotal %>% group_by(Region) %>% 
  # filter(!is.na(Region)) %>% 
  # summarise(n=n()) %>% 
  rename(Type=Region) %>% 
  count() %>% 
  mutate(fre=n/2672838)

df_variantsType2 <- dftotal %>% group_by(Variant_type) %>% 
  filter(!is.na(Variant_type)) %>% 
  # summarise(n=n()) %>% 
  rename(Type=Variant_type) %>% 
  count() %>% 
  mutate(fre=n/2672838)

df_variantsType3 <- dftotal %>% group_by(Effect_VEP) %>% 
  filter(!is.na(Effect_VEP)) %>% 
  summarise(n=n()) %>% 
  rename(Type=Effect_VEP)

print("Effect_VEP class:")
df_variantsType3$Type

df_variantsType4 <- dftotal %>% group_by(Effect_snpEff) %>% 
  filter(!is.na(Effect_snpEff)) %>% 
  summarise(n=n()) %>% 
  rename(Type=Effect_snpEff)

print("=======================")
print("Effect_SnpEff class:")
df_variantsType4$Type


colB <- brewer.pal(9,"Blues")[c(2,3,4,5,6,7)]

p <- df_variantsType1 %>% 
ggplot(aes(x=Type,y=n,fill=Type))+
  geom_bar(stat = 'identity',position = "stack")+
  # geom_text(aes(label=n),size=3,nudge_y = 1)+
  geom_text(aes(label=paste(n,"\n(",round(fre,4),")",sep = "")),size=3,
            position = position_dodge(0.9))+
  labs(x="", y="SNP count")+
  scale_fill_brewer(palette = "Set2")+
  # scale_fill_manual(values = colB)+
  # scale_fill_brewer(palette = "Blues")+
  theme_classic()+
  theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(legend.position = "none",plot.margin=unit(rep(1,4),'lines'),
          axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 10))

p

p <- df_variantsType2 %>% 
ggplot(aes(x=Type,y=n,fill=Type))+
  geom_bar(stat = 'identity',position = "stack")+
  # geom_text(aes(label=n),size=3,nudge_y = 1)+
  geom_text(aes(label=paste(n,"\n(",round(fre,4),")",sep = "")),size=3,
            position = position_dodge(0.9))+
  labs(x="", y="SNP count")+
  scale_fill_brewer(palette = "Set2")+
  # scale_fill_manual(values = colB)+
  # scale_fill_brewer(palette = "Blues")+
  theme_classic()+
  theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(legend.position = "none",plot.margin=unit(rep(1,4),'lines'),
          axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 10))

p

```

## 同义非同义区分亚基因组

```{r}
factor1 <- "Region"
factor2 <- "Variant_type"
factor3 <- "Effect_VEP"
factor4 <- "Effect_snpEff"

df_variantsType1 <- dftotal %>% group_by(Region) %>% 
  # filter(!is.na(Region)) %>% 
  summarise(n=n()) %>% 
  rename(Type=Region)

df_variantsType2 <- dftotal %>% group_by(Variant_type) %>% 
  filter(!is.na(Variant_type)) %>% 
  summarise(n=n()) %>% 
  rename(Type=Variant_type)


df_synNonsynbySub <- dftotal %>% 
  filter(!is.na(Variant_type)) %>% 
  group_by(Sub) %>% 
  count(Variant_type) %>% 
  mutate(fre=n/2672838) %>% 
  arrange(Variant_type)


### *****************************plot***********************************************
### 只画柱状图 dodge
p <- df_synNonsynbySub %>%
ggplot(aes(x=Sub,y=n,fill=Variant_type))+
  geom_bar(stat = 'identity',position = "dodge")+
  geom_text(aes(label=paste(n,"\n(",round(fre,4),")",sep = "")),size=3,
            position = position_dodge(0.9))+
  # geom_text(aes(label=paste(n,"\n(",round(fre,2),")",sep = "")),size=3,nudge_y = 1)+
  labs(x="", y="SNP count")+
  scale_fill_manual(values = brewer.pal(9,"Set2"))+
  theme_classic()+
  # theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(
    # legend.position = "none",
    plot.margin=unit(rep(1,4),'lines'),
    axis.line = element_line(size=0.4, colour = "black"),
    text = element_text(size = 10))

p

```




#### --------------------

## \*\*Fig2:variants count all region

```{r}
### "Intron" "UTR_3"  "UTR_5" 
df1 <- dftotal %>% group_by(Region) %>% 
  # filter(!is.na(Region)) %>% 
  summarise(n=n()) %>% 
  rename(Type=Region) %>% 
  filter(!Type=="CDS")

### "NONSYNONYMOUS" "STOP-GAIN"     "SYNONYMOUS" 
df2 <- dftotal %>% group_by(Variant_type) %>% 
  filter(!is.na(Variant_type)) %>% 
  summarise(n=n()) %>% 
  rename(Type=Variant_type) %>% 
  filter(!Type %in%  c("NONCODING","START-LOST","STOP-LOSS")) 

# expression(paste("Length-Freq of", italic("E. coruscans"), "by Gear")) 

### "Essential splice"
### 主要发生在内含子区域
df3 <- dftotal %>% 
  filter(Impact_VEP=="HIGH") %>%
  filter(str_detect(Effect_VEP,"splice_")) %>% 
  # filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>%
  count() %>% 
  mutate(Type="Essential splice") %>% 
  relocate(Type,.before = n)

### Intergenic
df4 <- tibble(Type="Intergenic",n=193285480)
  
df <- rbind(df1,df2,df3,df4)

factorA <- c("Intergenic","Intron","UTR_5","UTR_3","SYNONYMOUS","NONSYNONYMOUS","Essential splice","STOP-GAIN")
labelsA <- c("Intergenic","Intron","UTR_5","UTR_3","Synonymous","Nonsynonymous","Essential splice","Nonsense")
colB <- c(rep("gray70",4),"#004680","#e69f00","darkred","darkred")

colB <- rev(colB)

df$Type <- factor(df$Type,levels = rev(factorA),labels = rev(labelsA))

```

### \*plot

```{r}
library(ggbreak)
p <- df %>% 
  ggplot(aes(x= Type,y=n,fill=Type))+
  geom_bar(stat = 'identity',position = "stack")+
  # geom_text(aes(label=n),size=3,nudge_y = 1)+
  geom_hline(yintercept =c(400000,1400000), linetype="dashed", color = "gray")+
  geom_hline(yintercept =c(2000000,190000000), linetype="dashed", color = "gray")+
  scale_fill_manual(values = colB)+
  scale_y_break(c(400000, 1400000), scales = 0.5)+
  scale_y_break(c(2000000, 190000000), scales = 0.5)+
  scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  labs(x="",y="SNP count")+
  theme_classic()+
  coord_flip()+
  theme(axis.text.x= element_text(angle=90, vjust=0,hjust = 0))+
  theme(legend.position = "none",plot.margin=unit(rep(1,4),'lines'),
          axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 12))


p

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

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

```

#### --------------------

## GERP distribution

```{r}
gerpDir <- "data/006_Fig2/gerp/001_gerp_subset"

ChrRef <-rep(str_c(rep(1:7,each=3),rep(c("A","B","D"),7)),each=2)  
df_chrConvertion <- tibble(Chr=c(1:42),ChrRef=ChrRef,Sub=str_sub(ChrRef,2))

df_intergenic <- dir(gerpDir, full.names = T, pattern = ".txt.gz") %>% 
  map(read_tsv) %>% 
  bind_rows() %>% 
  filter(GeneFuture=="INTERGENIC") %>% 
  rename(Chr=ChrID) %>% 
  left_join(.,df_chrConvertion, by = "Chr")

df_region <- dftotal %>% 
  filter(Region %in%  c("Intron","UTR_5","UTR_3")) %>%
  select(Chr,Pos,GERP=Gerp_16way,GeneFuture=Region) %>% 
  left_join(.,df_chrConvertion, by = "Chr") 

df_variantsType <- dftotal %>% 
  filter(Variant_type %in%  c("NONSYNONYMOUS","SYNONYMOUS","STOP-GAIN")) %>%
  select(Chr,Pos,GERP=Gerp_16way,GeneFuture=Variant_type) %>% 
  left_join(.,df_chrConvertion, by = "Chr") 

df_splice <- dftotal %>% 
  filter(Impact_VEP=="HIGH") %>%
  filter(str_detect(Effect_VEP,"splice_")) %>% 
  select(Chr,Pos,GERP=Gerp_16way) %>% 
  mutate(GeneFuture="Essential splice") %>% 
  left_join(.,df_chrConvertion, by = "Chr") 

df <- rbind(df_intergenic,df_region,df_variantsType,df_splice)

factorA <- c("INTERGENIC","Intron","UTR_5","UTR_3","SYNONYMOUS","NONSYNONYMOUS","Essential splice","STOP-GAIN")
labelsA <- c("Intergenic","Intron","UTR_5","UTR_3","Synonymous","Nonsynonymous","Essential splice","Nonsense")
colB <- c(rep("gray70",4),"#004680","#e69f00","darkred","darkred")

df$GeneFuture <- factor(df$GeneFuture,levels = factorA,labels = labelsA)

dfgerpdis <- df %>% 
    slice_sample(prop = 0.3) 

```

### plot

```{r}
### histogram
p <- dfgerpdis %>%
  group_by(GeneFuture) %>% 
  ggplot(aes(x=GERP,fill=GeneFuture))+
  geom_histogram(aes(y=..ndensity..), alpha=0.95)+
  facet_grid(GeneFuture~., scales = "free")+
  # facet_grid(.~GeneFuture, scales = "free")+
  labs(x="GERP score",y="Density")+
  scale_y_continuous(breaks = seq(0,1,0.5))+
  scale_fill_manual(values = colB)+
  scale_color_manual(values = colB)+
  theme(
    strip.background = element_blank(),
    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(plot = p, "~/Documents/test.pdf", width = 3, height = 5)

```

### 求平均值、中位数
```{r}
df_ave <- dfgerpdis %>% 
  group_by(GeneFuture) %>%
  summarise(mean=mean(GERP,na.rm=TRUE))
  
```




## GERP bin  by syn nonsyn del proportion

```{r}
df <- dftotal %>% filter(Region=="CDS",Variant_type=="SYNONYMOUS"|Variant_type=="NONSYNONYMOUS") %>%
  # mutate(Sub=str_sub(Transcript, 9, 9)) %>% 
  mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(Derived_PolyPhen2_prob = ifelse(is.na(Derived_PolyPhen2_prob),100,Derived_PolyPhen2_prob)) %>% 
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>%
  filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>% 
  mutate(IfRefisAnc = if_else(Ref == Ancestral,"Anc","Der")) %>% 
  ### data deal
  filter(Gerp_16way > 0) %>%
  # mutate(bin=cut_interval(Gerp_16way,n=10)) %>% 
  mutate(bin=cut(Gerp_16way,breaks=c(0,0.3,0.6,0.9,1.2,1.49,2.1,2.4,2.7,3,3.3))) %>% 
  group_by(bin) %>% 
  count(LoadGroup) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(desc(bin)) 

dftest <- df %>% 
  filter(LoadGroup=="Nonsynonymous")
  

#### facet by sub
colB <- c("#d5311c","#e69f00","#004680")
p <- df %>% 
  ggplot(aes(x=bin,y=fre,group=LoadGroup,color=LoadGroup))+
  geom_line(aes(linetype=LoadGroup),size=1)+
  # scale_linetype_manual(values=c("dotdash","dotdash","dotdash"))+
  # scale_linetype_manual(values=c("dotdash","dotdash","dotdash"))+
  labs(y="Fraction in each category",x="GERP score")+
  geom_vline(xintercept = 5,linetype="dashed",color="blue")+
  scale_color_manual(values = colB)+
  scale_y_continuous(limits = c(0,1))+
  scale_x_discrete(labels = c("0","0.3","0.6","0.9","1.2","1.5","2.1","2.4","2.7","3","3.3"))+
  # coord_cartesian(xlim = c(14,20))+
  theme_classic()+
  # theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(plot.margin=unit(rep(1,4),'lines'),
        # legend.position = "none",
        legend.position = c(0.8,0.8),
        axis.line = element_line(size=0.4, colour = "black"),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 = 4)
# ggsave(p,filename = "~/Documents/test.pdf",height = 2.36,width = 2.36*4/3)
ggsave(p,filename = "~/Documents/test.pdf",height = 6.81,width = 9.5, units = "cm")


```

#### --------------------
## Deleterious count
### 1.区分亚基因组

```{r}
### GERP & PPH2 方法定义有害突变
df <- dftotal %>% 
  mutate(Sub=str_sub(Transcript, 9, 9)) %>% 
  filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>% 
  # filter(!is.na(DAF)) %>%
  mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>%
  group_by(Sub) %>% 
  count(LoadGroup) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(Sub,desc(LoadGroup))

### GERP 方法定义有害突变
df2 <- dftotal %>% 
  mutate(Sub=str_sub(Transcript, 9, 9)) %>% 
  filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>%
  filter(!is.na(DAF)) %>%
  # mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>% 
  mutate(LoadGroup=ifelse(Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=2.14, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>% 
  group_by(Sub) %>% 
  count(LoadGroup) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(Sub,desc(LoadGroup))


### *****************************plot***********************************************
colB <- c("#d5311c","#e69f00","#004680")
df$LoadGroup <- factor(df$LoadGroup,levels = c("Deleterious","Nonsynonymous","Synonymous"))

### 只画柱状图 dodge
p <- df %>%
ggplot(aes(x=Sub,y=n,fill=LoadGroup))+
  geom_bar(stat = 'identity',position = "dodge")+
  geom_text(aes(label=paste(n,"\n(",round(fre,2),")",sep = "")),size=3,
            position = position_dodge(0.9))+
  # geom_text(aes(label=paste(n,"\n(",round(fre,2),")",sep = "")),size=3,nudge_y = 1)+
  labs(x="", y="SNP count")+
  scale_fill_manual(values = colB)+
  theme_classic()+
  # theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(legend.position = "none",plot.margin=unit(rep(1,4),'lines'),
          axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 10))

p

### 只画柱状图 stack
p <- df %>% 
  ggplot(aes(x=Sub, y=n,fill=LoadGroup)) +
  geom_bar(stat="identity",alpha = 1)+  #,position=position_dodge()
  geom_text(aes(label=  paste(prettyNum(n,big.mark = ","),"(",round(fre,2),")")),position =  position_stack(),vjust=1.5, size=3) +
  scale_fill_manual(values=colB)+
  labs(y="SNP count",x="")+
  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 = 12)
        ,legend.position = 'none'
        )

p

### 只画柱状图 stack 百分比是1
p <- df %>% 
  ggplot(aes(x=Sub, y=n,fill=LoadGroup)) +
  geom_bar(stat="identity",alpha = 1,position = "fill")+  #,position=position_dodge()
  geom_text(aes(label=  paste(prettyNum(n,big.mark = ","),"(",round(fre,2),")")),position = position_fill(),vjust=1.5, size=3) +
  scale_fill_manual(values=colB)+
  labs(y="Proportion",x="")+
  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 = 12)
        ,legend.position = 'none'
        )

p

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

```

### 2.不区分亚基因组

```{r}
### 不区分亚基因组
### GERP & PPH2 方法定义有害突变
df <- dftotal %>% 
  mutate(Sub=str_sub(Transcript, 9, 9)) %>% 
  filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>%
  filter(!is.na(DAF)) %>%
  # mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  # mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>% 
  # filter(Sub=="B") %>%
  count(LoadGroup) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(desc(LoadGroup))


### GERP 方法定义有害突变
df2 <- dftotal %>% 
  mutate(Sub=str_sub(Transcript, 9, 9)) %>% 
  filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>%
  filter(!is.na(DAF)) %>%
  # mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>% 
  mutate(LoadGroup=ifelse(Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=2.14, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>% 
  # filter(Sub=="B") %>% 
  count(LoadGroup) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(desc(LoadGroup))


### *****************************plot***********************************************
colB <- c("#d5311c","#e69f00","#004680")
df$LoadGroup <- factor(df$LoadGroup,levels = c("Deleterious","Nonsynonymous","Synonymous"))

### 只画柱状图 dodge
p <- df %>%
ggplot(aes(x=LoadGroup,y=n,fill=LoadGroup))+
  geom_bar(stat = 'identity',position = "dodge")+
  geom_text(aes(label=paste(n,"\n(",round(fre,2),")",sep = "")),size=1,
            position = position_dodge(0.9))+
  # geom_text(aes(label=paste(n,"\n(",round(fre,2),")",sep = "")),size=3,nudge_y = 1)+
  labs(x="", y="SNP count")+
  scale_fill_manual(values = colB)+
  theme_classic()+
  theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(legend.position = "none",plot.margin=unit(rep(1,4),'lines'),
          axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 7))

p

ggsave(p,filename = "~/Documents/test.pdf",height = 1.725,width = 2.3)

```

### 3.不区分亚基因组- High impact/Del/Nonsyn

```{r}
### 不区分亚基因组
### GERP & PPH2 方法定义有害突变
df_highImpact <- dftotal %>% 
  filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>%
  filter(!is.na(DAF)) %>%
  filter(Impact_VEP == "HIGH") %>% 
  count() %>% 
  add_column(LoadGroup="High impact", .before = "n")

  
  
df_del_nonsyn_syn <- dftotal %>% 
  mutate(Sub=str_sub(Transcript, 9, 9)) %>% 
  filter(!is.na(Ancestral),Ancestral==Ref|Ancestral==Alt) %>%
  filter(!is.na(DAF)) %>%
  # mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  # mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>% 
  # filter(Sub=="B") %>%
  count(LoadGroup) 

df <- df_del_nonsyn_syn %>% 
  bind_rows(df_highImpact) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(desc(LoadGroup))

### *****************************plot***********************************************
colB <- c("darkred","#d5311c","#e69f00","#004680")
df$LoadGroup <- factor(df$LoadGroup,levels = c("High impact","Deleterious","Nonsynonymous","Synonymous"))

### 只画柱状图 dodge
p <- df %>%
ggplot(aes(x=LoadGroup,y=n,fill=LoadGroup))+
  geom_bar(stat = 'identity',position = "dodge")+
  geom_text(aes(label=paste(n,"\n(",round(fre,2),")",sep = "")),size=3,
            position = position_dodge(0.9))+
  # geom_text(aes(label=paste(n,"\n(",round(fre,2),")",sep = "")),size=3,nudge_y = 1)+
  labs(x="", y="SNP count")+
  scale_fill_manual(values = colB)+
  theme_classic()+
  theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(legend.position = "none",plot.margin=unit(rep(1,4),'lines'),
          axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 10))

p

ggsave(p,filename = "~/Documents/test.pdf",height = 2.9,width = 2.7)

```



#### ===== Gene level ================

## 韦恩图 gene common 和 gene rare 列表交集

```{r}
### https://rdrr.io/cran/VennDiagram/man/venn.diagram.html
### 含有 common 基因的基因列表
dfcommon <- dftotal_withAncestral %>% 
  group_by(Transcript) %>% 
  count(CommonOrRare) %>% 
  filter(CommonOrRare == "Common") %>% 
  dplyr::select(Transcript)

### 含有 rare 基因的基因列表
dfrare <- dftotal_withAncestral %>% 
  group_by(Transcript) %>% 
  count(CommonOrRare) %>% 
  filter(CommonOrRare == "Rare") %>%
  dplyr::select(Transcript)

library(VennDiagram)
#Make the plot
venn.diagram(
  x = list(dfcommon$Transcript,dfrare$Transcript), 
  main = "common and rare genes",
  category.names = c("common" , "rare" ),
  filename = "~/Documents/gene_statistic.png",
  output = T ,
  
  # Output features
  imagetype="png" ,
  # height = 480 , 
  # width = 480 , 
  # resolution = 300,
  # compression = "lzw",
  
  # Circles
  lwd = 3,
  col=c("#440154ff", '#21908dff'),
  fill = c(alpha("#440154ff",0.3), alpha('#21908dff',0.3)),
  
  # Numbers
  cex = 0.5,
  fontfamily = "sans",
  cat.cex = 0.3,
  cat.default.pos = "outer",
  ## cat.pos = c(-27, 27, 135),
  ## cat.dist = c(0.055, 0.055, 0.085),
  cat.pos = c(-27, 27),
  cat.dist = c(0.055, 0.055),
  cat.fontfamily = "sans",
  cat.col = c("#440154ff", '#21908dff')
  # rotation = 1
)
```

## Gene count summary

```{r}
### 1.计算多少个基因含有突变 47490
genenum <- length(unique(dftotal_withAncestral$Transcript))
print(paste("Gene num with SNP: ",genenum))

### 1.基因分类为 common rare commonAndRare

df_geneCategory <- dftotal_withAncestral %>% 
  dplyr::select(ID,Transcript,CommonOrRare) %>% 
  group_by(Transcript,CommonOrRare) %>% 
  count() %>% 
  pivot_wider(names_from = CommonOrRare,values_from = n,values_fill = 0) %>% 
  mutate(GeneGroup=ifelse(Common>0 & Rare >0 ,"WithCommonandRareSNPs",
                          ifelse(Common>0 & Rare ==0,"OnlyCommonSNPs",
                                 ifelse(Common==0 & Rare >0,"OnlyRareSNPs",NA))))

  
### 
df_geneCategory_num <- df_geneCategory %>% 
  group_by(GeneGroup) %>% count()

print(df_geneCategory_num)

p <- df_geneCategory_num %>% 
  ggplot(aes(x=GeneGroup, y=n,fill=GeneGroup)) +
  geom_bar(stat="identity",alpha = 1)+  #,position=position_dodge()
  geom_text(aes(label=prettyNum(n,big.mark = ",")),position =  position_stack(),vjust=1.5, size=3) +
  scale_fill_brewer(palette = "Set3")+
  labs(y="Gene count",x="")+
  # 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 = element_blank(),
        legend.position = 'none',
        text = element_text(size = 12))

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


### ======================================
### 2.基因分类为 common-del/nodel  rare-del/nodel commonAndRare-del/nodel

df_geneCategory2 <- dftotal_withAncestral %>% 
  dplyr::select(ID,Transcript,CommonOrRare,LoadGroup) %>% 
  mutate(IfDel=ifelse(LoadGroup=="Deleterious","Del","NoDel")) %>%
  group_by(Transcript,IfDel) %>% 
  count() %>% 
  pivot_wider(names_from = IfDel,values_from = n,values_fill = 0) %>% 
  mutate(DelGroup=ifelse(Del>0 ,"HasDeleteriousSNP","NoDeleteriousSNP")) %>% 
  left_join(df_geneCategory,by="Transcript")

df_geneCategory2_num <- df_geneCategory2 %>% 
  group_by(GeneGroup,DelGroup) %>% 
  count()
 
print(df_geneCategory2_num)

### ================== 画图 ====================
colB <- c("#d5311c","gray") ###一个是有害突变，一个是没有有害突变
p <- df_geneCategory2_num %>% 
  ggplot(aes(x=GeneGroup, y=n,fill=DelGroup)) +
  geom_bar(stat="identity",alpha = 1)+  #,position=position_dodge()
  geom_text(aes(label=prettyNum(n,big.mark = ",")),position =  position_stack(),vjust=1.5, size=3) +
  scale_fill_manual(values=colB)+
  labs(y="Gene count",x="")+
  # 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 = 12)
        # ,legend.position = 'none'
        )

p
ggsave(p,filename = "~/Documents/test.pdf",height = 3,width = 4)
ggsave(p,filename = "~/Documents/p.pdf",height = 6.33,width = 7.82,units = "cm")


```

## Gene MAF LoadGroup

![](figs/FigS_geneCount/misc/GeneCount_pipeline.png)

```{r}

## 每个基因至少出现3次，即每个基因至少含有Synonymous Nonsynonymous Deleterious
## 过滤后，每个基因按照rare的总数进行排序，由低到高，并进行间隔抽样，抽取100个试试看

### 每个基因含有的rare变异的分布范围
### 结论：合计有46093个基因含有稀有变异
df_common_rare <- dftotal_withAncestral %>%
  group_by(Transcript,CommonOrRare) %>% 
  count() %>% 
  ungroup() %>% 
  filter(CommonOrRare == "Rare") %>% 
  select(Transcript,RareVariantsCount=n)

print(str_c(nrow(df_common_rare)," genes which has rare SNPs."))
print(str_c("The range of rare SNPs is ",min(df_common_rare$RareVariantsCount)," - ",max(df_common_rare$RareVariantsCount)))

print(str_c("平均每个基因含有稀有变异的数目在",mean(df_common_rare$RareVariantsCount)))
### "平均每个基因含有稀有变异的数目在6.46545028529278"

df_stage1 <- dftotal_withAncestral %>% 
  group_by(Transcript,CommonOrRare,LoadGroup) %>% 
  count(LoadGroup) %>% 
  mutate(Sub = str_sub(Transcript,9,9)) %>% 
  pivot_wider(names_from = CommonOrRare, values_from  = n,values_fill = 0)

### 在几个基因中鉴定到了有害突变
df_delGeneNumber <- dftotal_withAncestral %>% 
  group_by(Transcript,LoadGroup) %>% 
  count(LoadGroup) %>% 
  filter(LoadGroup=="Deleterious")

print(str_c(nrow(df_delGeneNumber)," 个基因含有有害突变。","\n","平均每个基因含有有害突变的数目在", mean(df_delGeneNumber$n)))

df_delNoly1delSNP <- df_delGeneNumber %>% 
  filter(n==1)

print(str_c(nrow(df_delNoly1delSNP)," 个基因只含有1个有害突变。"))


df_filter <- df_stage1 %>% 
  group_by(Transcript) %>% 
  count() %>% 
  filter(n == 3) %>% 
  select(Transcript) %>% 
  left_join(.,df_common_rare,by="Transcript") %>% 
  filter(!is.na(RareVariantsCount)) %>% 
  ungroup()
  ### nuique 一下计数

df_sample <- df_filter %>% 
  distinct(RareVariantsCount,.keep_all=TRUE) %>% 
  arrange(RareVariantsCount)

df <- df_sample %>% 
  left_join(.,df_stage1,by="Transcript") %>% 
  rename(Gene = Transcript)

## 处理抽样基因的CDS长度
df_geneHC <- read_tsv("data/006_Fig2/wheat_v1.1_nonoverlap_addPos.txt") %>%
  select(Longest_transcript,CDSLength) %>% 
  rename(Transcript=Longest_transcript)

df_geneCDSlength <- df_sample %>% 
  left_join(.,df_geneHC,by="Transcript") %>% 
  rename(Gene = Transcript)

```

### plot
```{r}

colB <- c("#d5311c","#e69f00","#004680")
p <- ggplot()+
  geom_bar(data = df, mapping = aes(x=reorder(Gene,RareVariantsCount) ,y=-Rare,fill=LoadGroup),stat = 'identity')+
  geom_bar(data = df, mapping = aes(x=Gene,y=Common,fill = LoadGroup),stat = 'identity')+
  # geom_segment(x=0,xend= 100000,y=0,yend=0,size=1,color='white')+
  geom_line(data = df_geneCDSlength,mapping = aes(x=Gene,y= 7*CDSLength/1000,group = 1),color = "gray")+
  annotate("text",x=-Inf,y= 0 ,hjust=-0.3,vjust=4,label="MAF <= 5%",size=3)+
  annotate("text",x=-Inf,y= 0,hjust=-0.3,vjust=-4,label="MAF > 5%",size=3)+
  scale_y_continuous(limit = c(-200,120), breaks = seq(-200,100,50),labels = c(200,150,100,50,0,50,100)
                     ,sec.axis = sec_axis(~./7, name="CDS length (kb)",breaks = seq(0,15,5))
                     )+
  xlab("Gene rank (by number of rare coding SNVs)")+ylab("Number of variants")+
  scale_fill_manual(values = colB,name="")+
  scale_x_discrete(limits = df_sample$Transcript, breaks = c("TraesCS1A02G169600.1","TraesCS1A02G040600.1","TraesCS1A02G091500.1","TraesCS1B02G160800.1")
                   ,labels = c(1,25,50,75)
                   )+
  # theme(axis.text.x = element_blank(),axis.ticks.x = element_blank(),axis.line.x = 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 = 11))
p
# ggsave(p,filename = "figs/Fig2/misc/Gene_MAF_LoadGroup.pdf",height = 3,width = 6)
ggsave(p,filename = "~/Documents/test.pdf",height = 2,width = 3.6)

```

### 稀有变异数分布
```{r}
df_filter %>% 
  ggplot(aes(x=.$RareVariantsCount))+
  geom_histogram()
```
### 非同义突变数目分布
```{r}
df_highImpact <- dftotal %>% 
  filter(Impact_VEP=="HIGH") %>% 
  mutate(LoadGroup="HighImpact")

df1 <- dftotal_withAncestral %>% 
  select(-IfRefisAnc,CommonOrRare) %>% 
  bind_rows(df_highImpact) %>% 
  # filter(LoadGroup=="Nonsynonymous") %>% 
  group_by(LoadGroup) %>% 
  count(Transcript)
  # mutate(bin=cut_width(n,width = 3,boundary = 0)) %>% 
  # group_by(LoadGroup) %>% 
  # count(bin) %>% 
  # mutate(fre=n/sum(n))

colB <- c("#d5311c","darkred","#e69f00","#004680")
p <- df1 %>% 
  ggplot(aes(x=n,fill=LoadGroup,color=LoadGroup))+
  # geom_histogram(aes(y=..ndensity..), alpha=0.5,bins = 15)+
  geom_histogram(alpha=0.7,bins = 15)+
  # geom_density(alpha=.2, fill="#FF6666",)+
  # ggplot(aes(x=bin,y=fre,group=LoadGroup,fill=LoadGroup))+
  # geom_bar(stat = 'identity',position = "dodge")+
  labs(x="Number of derived SNPs in one gene", y="Gene number")+
  scale_fill_manual(values = colB)+
  scale_color_manual(values = colB)+
  facet_wrap(.~LoadGroup, scales = "free")+
  theme_classic()+
  theme(panel.grid = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(size=0.3, colour = "black"),
    legend.position = 'none',
    text = element_text(size = 12))

p

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

数目统计

```{r}
## 每个基因含有的变异数范围
df <- dftotal %>% filter(Region=="CDS",Variant_type=="SYNONYMOUS"|Variant_type=="NONSYNONYMOUS")%>% 
  mutate(Derived_SIFT = ifelse(is.na(Derived_SIFT),100,Derived_SIFT)) %>% 
  mutate(Gerp = ifelse(is.na(Gerp),-100,Gerp)) %>% 
  mutate(LoadGroup = ifelse(Variant_type=="NONSYNONYMOUS",ifelse(Derived_SIFT<0.05&Gerp>=1,"Deleterious","Nonsynonymous"),"Synonymous")) %>% 
  filter(!Ancestral == "NA") %>% 
  filter(Ancestral== Ref|Ancestral==Alt) %>% 
  mutate(IfRefisAnc = if_else(Ref == Ancestral,"Anc","Der")) %>% 
  mutate(CommonOrRare = if_else(Maf <= 0.05, "Rare","Common")) %>%  ### 特别注意，这里根据MAF进行分组，按照维基百科的定义
  count(Transcript) 

## 1-199 个变异

## 每个基因含有的common和rare变异的分布范围
df <- dftotal %>% filter(Region=="CDS",Variant_type=="SYNONYMOUS"|Variant_type=="NONSYNONYMOUS")%>% 
  mutate(Derived_SIFT = ifelse(is.na(Derived_SIFT),100,Derived_SIFT)) %>% 
  mutate(Gerp = ifelse(is.na(Gerp),-100,Gerp)) %>% 
  mutate(LoadGroup = ifelse(Variant_type=="NONSYNONYMOUS",ifelse(Derived_SIFT<0.05&Gerp>=1,"Deleterious","Nonsynonymous"),"Synonymous")) %>% 
  filter(!Ancestral == "NA") %>% 
  filter(Ancestral== Ref|Ancestral==Alt) %>% 
  mutate(IfRefisAnc = if_else(Ref == Ancestral,"Anc","Der")) %>% 
  mutate(CommonOrRare = if_else(Maf <= 0.05, "Rare","Common")) %>%  ### 特别注意，这里根据MAF进行分组，按照维基百科的定义
  group_by(CommonOrRare) %>% 
  count(Transcript) %>% 
  summarise(maxNum = max(n),minNum = min(n))



## 每个基因含有的不同loadGroup变异的分布范围
df <- dftotal %>% filter(Region=="CDS",Variant_type == "SYNONYMOUS"|Variant_type=="NONSYNONYMOUS")%>% 
  mutate(Derived_SIFT = ifelse(is.na(Derived_SIFT),100,Derived_SIFT)) %>% 
  mutate(Gerp = ifelse(is.na(Gerp),-100,Gerp)) %>% 
  mutate(LoadGroup = ifelse(Variant_type=="NONSYNONYMOUS",ifelse(Derived_SIFT<0.05&Gerp>=1,"Deleterious","Nonsynonymous"),"Synonymous")) %>% 
  filter(!Ancestral == "NA") %>% 
  filter(Ancestral== Ref|Ancestral==Alt) %>% 
  mutate(IfRefisAnc = if_else(Ref == Ancestral,"Anc","Der")) %>% 
  mutate(CommonOrRare = if_else(Maf <= 0.05, "Rare","Common")) %>%  ### 特别注意，这里根据MAF进行分组，按照维基百科的定义
  group_by(LoadGroup) %>% 
  count(Transcript) %>% 
  summarise(maxNum = max(n),minNum = min(n))

```

#### ======= Sample size==============
## 正文用图：CDF (scale_x)
```{r}
pfun_CDF_scale_x <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
### 找到每种分类下的最大值
aoFindMaxLine <- function(dfinput){
  maxtaxanum = max(dfinput$TaxaNum)
  dfout <- dfinput %>% filter(TaxaNum==maxtaxanum)
  return(dfout)
}

### 发现每个分类的最大值 "Group"         "TaxaNum"       "ObservedCount" "Loop"    
dfpoint <- df %>% group_by(Group) %>% 
  group_modify(~{.x %>% aoFindMaxLine()}) %>% ### 将每个分类的最大值合并入文件中
  select(Group,TotalNum = ObservedCount) %>% 
  distinct(.,Group,.keep_all = T) ###每种分类有多个循环loop，所以有重复，需要去除

p <- df %>% 
  left_join(.,dfpoint,by = "Group") %>% 
  filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  mutate(CDF = ObservedCount/TotalNum) %>% 
  ggplot(aes(x=TaxaNum,y=CDF))+
  stat_summary(geom = "errorbar"
               ,width = 0.1
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Number of sites observed\nto be variable")+
  scale_color_manual(values = colB)+
  # labs(x="Sample size",y="CDF")+
  labs(x="Sample size",y="Percent observed")+
  scale_x_log10(limits = c(1,maxSampleSize))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = 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 = file.path("~/Documents/",str_c(outname,"_CDF_scale_x.pdf")),height = 2.22,width = 1.8)
}

```

### p0函数调用
```{r}

df <- read_tsv("data/010_sampleSize2VariantsDiscovery/ABsub.txt")

pfun_CDF_scale_x(df,1026,"ABsub")

### -----------------------------------
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/Dsub.txt")
pfun_CDF_scale_x(df,1026,"Dsub")

dftest <- df %>% 
  group_by(TaxaNum,Group) %>% 
  summarise(mean=mean(ObservedCount))

```


#### ---------------------
## 分析- LoadGroup
### *methd1:count (no scale,scale_x,scale_xy)
```{r}
pfun_count <- function(df,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")

p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=ObservedCount))+
  stat_summary(geom = "errorbar"
               # ,width = 8   
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Number of sites observed\nto be variable")+
  scale_color_manual(values = colB)+
  scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = 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 = file.path("~/Documents/",str_c(outname,"_count.pdf")),height = 2.62,width = 3.5)
}

pfun_count_scale_x <- function(df,maxSampleSize,outname){
  
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
  p <- df %>% filter(Group %in% group1) %>% 
    mutate(Group=factor(Group,levels = group1)) %>%
    ggplot(aes(x=TaxaNum,y=ObservedCount))+
    stat_summary(geom = "errorbar"
               # ,width = 8   
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
    stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
    stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Number of sites observed\nto be variable")+
  # scale_color_manual(values = colB)+
  scale_x_log10(limits = c(1,maxSampleSize))+
  scale_color_manual(values = colB)+
  scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  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 = file.path("figs/Fig_sampleSize2VariantsDiscovery/misc2",str_c(outname,"_count_scale_x.pdf")),height = 2.22,width = 3.6)
ggsave(p,filename = file.path("~/Documents/",str_c(outname,"_count_scale_x.pdf")),height = 2.22,width = 3.5)
}

pfun_count_scale_xy <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=ObservedCount))+
  stat_summary(geom = "errorbar"
               # ,width = 8   
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Number of sites observed\nto be variable")+
  # scale_color_manual(values = colB)+
  scale_color_manual(values = colB)+
  scale_x_log10(limits = c(1,maxSampleSize))+
  scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  scale_y_log10()+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  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 = file.path("~/Documents/",str_c(outname,"_count_scale_xy.pdf")),height = 2.62,width = 3.5)
}

```

### p0 函数调用
```{r}
# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/ABsub.txt")

df <- read_tsv("data/010_sampleSize2VariantsDiscovery/ABsub.txt")
# pfun_count(df,"ABsub")
pfun_count_scale_x(df,1026,"ABsub")
# pfun_count_scale_xy(df,1026,"ABsub")

# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/Dsub.txt")
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/Dsub.txt")

# pfun_count(df,"Dsub")
# pfun_count_scale_x(df,1026,"Dsub")
# pfun_count_scale_xy(df,1026,"Dsub")

```


### *methd2:CDF (no scale,scale_x,scale_xy)
```{r}
pfun_CDF <- function(df,outname){

  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
### 找到每种分类下的最大值
aoFindMaxLine <- function(dfinput){
  maxtaxanum = max(dfinput$TaxaNum)
  dfout <- dfinput %>% filter(TaxaNum==maxtaxanum)
  return(dfout)
}

### 发现每个分类的最大值 "Group"         "TaxaNum"       "ObservedCount" "Loop"    
dfpoint <- df %>% group_by(Group) %>% 
  group_modify(~{.x %>% aoFindMaxLine()}) %>% ### 将每个分类的最大值合并入文件中
  select(Group,TotalNum = ObservedCount) %>% 
  distinct(.,Group,.keep_all = T) ###每种分类有多个循环loop，所以有重复，需要去除

p <- df %>% 
  left_join(.,dfpoint,by = "Group") %>% 
  filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  mutate(CDF = ObservedCount/TotalNum) %>% 
  ggplot(aes(x=TaxaNum,y=CDF))+
  stat_summary(geom = "errorbar"
               # ,width = 0.1
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Number of sites observed\nto be variable")+
  scale_color_manual(values = colB)+
  # labs(x="Sample size",y="CDF")+
  labs(x="Sample size",y="Percent observed")+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = 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 = file.path("~/Documents/",str_c(outname,"_CDF.pdf")),height = 2.62,width = 3.5)
}

pfun_CDF_scale_x <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
### 找到每种分类下的最大值
aoFindMaxLine <- function(dfinput){
  maxtaxanum = max(dfinput$TaxaNum)
  dfout <- dfinput %>% filter(TaxaNum==maxtaxanum)
  return(dfout)
}

### 发现每个分类的最大值 "Group"         "TaxaNum"       "ObservedCount" "Loop"    
dfpoint <- df %>% group_by(Group) %>% 
  group_modify(~{.x %>% aoFindMaxLine()}) %>% ### 将每个分类的最大值合并入文件中
  select(Group,TotalNum = ObservedCount) %>% 
  distinct(.,Group,.keep_all = T) ###每种分类有多个循环loop，所以有重复，需要去除

p <- df %>% 
  left_join(.,dfpoint,by = "Group") %>% 
  filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  mutate(CDF = ObservedCount/TotalNum) %>% 
  ggplot(aes(x=TaxaNum,y=CDF))+
  stat_summary(geom = "errorbar"
               ,width = 0.1
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Number of sites observed\nto be variable")+
  scale_color_manual(values = colB)+
  # labs(x="Sample size",y="CDF")+
  labs(x="Sample size",y="Percent observed")+
  scale_x_log10(limits = c(1,maxSampleSize))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = 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 = file.path("~/Documents/",str_c(outname,"_CDF_scale_x.pdf")),height = 2.22,width = 3.5)
}

pfun_CDF_scale_xy <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
### 找到每种分类下的最大值
aoFindMaxLine <- function(dfinput){
  maxtaxanum = max(dfinput$TaxaNum)
  dfout <- dfinput %>% filter(TaxaNum==maxtaxanum)
  return(dfout)
}

### 发现每个分类的最大值 "Group"         "TaxaNum"       "ObservedCount" "Loop"    
dfpoint <- df %>% group_by(Group) %>% 
  group_modify(~{.x %>% aoFindMaxLine()}) %>% ### 将每个分类的最大值合并入文件中
  select(Group,TotalNum = ObservedCount) %>% 
  distinct(.,Group,.keep_all = T) ###每种分类有多个循环loop，所以有重复，需要去除

p <- df %>% 
  left_join(.,dfpoint,by = "Group") %>% 
  filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  mutate(CDF = ObservedCount/TotalNum) %>% 
  ggplot(aes(x=TaxaNum,y=CDF))+
  stat_summary(geom = "errorbar"
               ,width = 0.1
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Number of sites observed\nto be variable")+
  scale_color_manual(values = colB)+
  # labs(x="Sample size",y="CDF")+
  labs(x="Sample size",y="Percent observed")+
  scale_x_log10(limits = c(1,maxSampleSize))+
  scale_y_log10()+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = element_blank(), 
      panel.background = element_blank(),
      axis.line = element_line(size=0.4, colour = "black"),
      legend.position = 'none',
      text = element_text(size = 12))


p
ggsave(p,filename = file.path("~/Documents/",str_c(outname,"_CDF_scale_xy.pdf")),height = 3,width = 4)
}

```


### p0 函数调用
```{r}
# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/ABsub.txt")
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/ABsub.txt")
# pfun_CDF(df,"ABsub")
pfun_CDF_scale_x(df,1026,"ABsub")
# pfun_CDF_scale_xy(df,1026,"ABsub")

# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/Dsub.txt")
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/Dsub.txt")

# pfun_CDF(df,"Dsub")
# pfun_CDF_scale_x(df,1026,"Dsub")
# pfun_CDF_scale_xy(df,1026,"Dsub")

```



### *methd3:Increasement (no scale,scale_x,scale_xy)
```{r}
pfun_Increasement <- function(df,maxSampleSize ,outname){ 

  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")

  ## 建立2个数据框，第一个是去除taxaNum为1的数据框 df_taxaNumNo1 ;第二个是去除taxaNum为最大值的数据框 df_taxaNumNomax
  df_taxaNumNo1 <- df %>% filter(TaxaNum != 1) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNomax <- df %>% filter(TaxaNum != maxSampleSize) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNo1$Increasement <- df_taxaNumNo1$ObservedCount - df_taxaNumNomax$ObservedCount
  dftotalnum <- df %>% filter(TaxaNum == maxSampleSize,Loop==1) %>% select(Group, ObservedCount)
  
  df <- df_taxaNumNo1 %>%left_join(dftotalnum,by="Group") %>% 
  mutate(percent = Increasement/ObservedCount.y)

p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=Increasement))+
  stat_summary(geom = "errorbar"
               ,width = 20   
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="New sites observed\nto be variable per additional sample")+
  scale_color_manual(values = colB)+
  scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = element_blank(),
      panel.background = element_blank(),
      axis.line = element_line(size=0.4, colour = "black"),
      legend.position = 'none',
      text = element_text(size = 12))

p
ggsave(p,filename = file.path("~/Documents/",str_c(outname,"_Increasement.pdf")),height = 3,width = 4)
}

pfun_Increasement_scale_x <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
  ## 建立2个数据框，第一个是去除taxaNum为1的数据框 df_taxaNumNo1 ;第二个是去除taxaNum为最大值的数据框 df_taxaNumNomax
  df_taxaNumNo1 <- df %>% filter(TaxaNum != 1) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNomax <- df %>% filter(TaxaNum != maxSampleSize) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNo1$Increasement <- df_taxaNumNo1$ObservedCount - df_taxaNumNomax$ObservedCount
  dftotalnum <- df %>% filter(TaxaNum == maxSampleSize,Loop==1) %>% select(Group, ObservedCount)

    df <- df_taxaNumNo1 %>%left_join(dftotalnum,by="Group") %>% 
  mutate(percent = Increasement/ObservedCount.y)
    
p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=Increasement))+
  stat_summary(geom = "errorbar"
               ,width = 40 
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="New sites observed\nto be variable per additional sample")+
  scale_color_manual(values = colB)+
  scale_x_log10(limits = c(100,maxSampleSize))+
  scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = 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 = file.path("~/Documents/",str_c(outname,"_Increasement_scale_x.pdf")),height = 2.22,width = 3.5)
}


pfun_Increasement_scale_xy <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
  ## 建立2个数据框，第一个是去除taxaNum为1的数据框 df_taxaNumNo1 ;第二个是去除taxaNum为最大值的数据框 df_taxaNumNomax
  df_taxaNumNo1 <- df %>% filter(TaxaNum != 1) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNomax <- df %>% filter(TaxaNum != maxSampleSize) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNo1$Increasement <- df_taxaNumNo1$ObservedCount - df_taxaNumNomax$ObservedCount
  dftotalnum <- df %>% filter(TaxaNum == maxSampleSize,Loop==1) %>% select(Group, ObservedCount)

    df <- df_taxaNumNo1 %>%left_join(dftotalnum,by="Group") %>% 
  mutate(percent = Increasement/ObservedCount.y)
    
p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=Increasement))+
  stat_summary(geom = "errorbar"
               ,width = 40   
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="New sites observed\nto be variable per additional sample")+
  scale_color_manual(values = colB)+
  scale_x_log10(limits = c(100,maxSampleSize))+
  scale_y_log10()+
  # scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = element_blank(),
      panel.background = element_blank(),
      axis.line = element_line(size=0.4, colour = "black"),
      legend.position = 'none',
      text = element_text(size = 12))

p
ggsave(p,filename = file.path("~/Documents/",str_c(outname,"_Increasement_scale_xy.pdf")),height = 3,width = 4)
}

```

### p0 函数调用
```{r}
# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/ABsub.txt")
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/ABsub.txt")

# pfun_Increasement(df,1026,"ABsub")
pfun_Increasement_scale_x(df,1026,"ABsub")
# pfun_Increasement_scale_xy(df,1026,"ABsub")

# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/Dsub.txt")
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/Dsub.txt")

# pfun_Increasement(df,850,"Dsub")
# pfun_Increasement_scale_x(df,850,"Dsub")
# pfun_Increasement_scale_xy(df,850,"Dsub")
```


### *methd4:Increasement CDF (no scale,scale_x,scale_xy)
```{r}
pfun_Increasement_CDF <- function(df,maxSampleSize ,outname){ 

  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
  ## 建立2个数据框，第一个是去除taxaNum为1的数据框 df_taxaNumNo1 ;第二个是去除taxaNum为最大值的数据框 df_taxaNumNomax
  df_taxaNumNo1 <- df %>% filter(TaxaNum != 1) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNomax <- df %>% filter(TaxaNum != maxSampleSize) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNo1$Increasement <- df_taxaNumNo1$ObservedCount - df_taxaNumNomax$ObservedCount
  dftotalnum <- df %>% filter(TaxaNum == maxSampleSize,Loop==1) %>% select(Group, ObservedCount)
  
  df <- df_taxaNumNo1 %>%left_join(dftotalnum,by="Group") %>% 
  mutate(percent = Increasement/ObservedCount.y)

p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=percent))+
  stat_summary(geom = "errorbar"
               ,width = 50   
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Percent of new sites observed\nto be variable per additional sample")+
  scale_color_manual(values = colB)+
  scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = element_blank(),
      panel.background = element_blank(),
      axis.line = element_line(size=0.4, colour = "black"),
      legend.position = 'none',
      text = element_text(size = 12))

p
ggsave(p,filename = file.path("~/Documents/",str_c(outname,"_Increasement_CDF.pdf")),height = 3,width = 4)
}

pfun_Increasement_CDF_scale_x <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
  ## 建立2个数据框，第一个是去除taxaNum为1的数据框 df_taxaNumNo1 ;第二个是去除taxaNum为最大值的数据框 df_taxaNumNomax
  df_taxaNumNo1 <- df %>% filter(TaxaNum != 1) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNomax <- df %>% filter(TaxaNum != maxSampleSize) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNo1$Increasement <- df_taxaNumNo1$ObservedCount - df_taxaNumNomax$ObservedCount
  dftotalnum <- df %>% filter(TaxaNum == maxSampleSize,Loop==1) %>% select(Group, ObservedCount)

    df <- df_taxaNumNo1 %>%left_join(dftotalnum,by="Group") %>% 
  mutate(percent = Increasement/ObservedCount.y)
    
p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=percent))+
  stat_summary(geom = "errorbar"
               # ,width = 50   
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Percent of new sites observed\nto be variable per additional sample")+
  scale_color_manual(values = colB)+
  scale_x_log10(limits = c(100,maxSampleSize))+
  # scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = 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 = file.path("~/Documents/",str_c(outname,"_Increasement_CDF_scale_x.pdf")),height = 2.22,width = 3.5)
}


pfun_Increasement_CDF_scale_xy <- function(df,maxSampleSize,outname){
  group1 <- c("Deleterious","Nonsynonymous","Synonymous")
  colB <- c("#d5311c","#e69f00","#004680")
  ## 建立2个数据框，第一个是去除taxaNum为1的数据框 df_taxaNumNo1 ;第二个是去除taxaNum为最大值的数据框 df_taxaNumNomax
  df_taxaNumNo1 <- df %>% filter(TaxaNum != 1) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNomax <- df %>% filter(TaxaNum != maxSampleSize) %>% arrange(TaxaNum,Loop,Group)
  df_taxaNumNo1$Increasement <- df_taxaNumNo1$ObservedCount - df_taxaNumNomax$ObservedCount
  dftotalnum <- df %>% filter(TaxaNum == maxSampleSize,Loop==1) %>% select(Group, ObservedCount)

    df <- df_taxaNumNo1 %>%left_join(dftotalnum,by="Group") %>% 
  mutate(percent = Increasement/ObservedCount.y)
    
p <- df %>% filter(Group %in% group1) %>% 
  mutate(Group=factor(Group,levels = group1)) %>%
  ggplot(aes(x=TaxaNum,y=percent))+
  stat_summary(geom = "errorbar"
               ,width = 50  
               ,fun.min = min
               ,fun.max = max
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
  stat_summary(fun = "mean"
               ,geom = "point"
               ,aes(color = Group)
               # ,alpha = 0.5
               )+ 
   stat_summary(fun = "mean"
               ,geom = "line"
               ,aes(color = Group)
               # ,alpha = 0.9
               )+
  labs(x="Sample size",y="Percent of new sites observed\nto be variable per additional sample")+
  scale_color_manual(values = colB)+
  scale_x_log10(limits = c(100,maxSampleSize))+
  scale_y_log10()+
  # scale_y_continuous(labels = scales::comma_format(big.mark = ','))+
  # facet_wrap(~Group,ncol = 3,scales = "free")+
  theme_classic()+
  theme(panel.grid = element_blank(),
      panel.background = element_blank(),
      axis.line = element_line(size=0.4, colour = "black"),
      legend.position = 'none',
      text = element_text(size = 12))

p
ggsave(p,filename = file.path("~/Documents/",str_c(outname,"_Increasement_CDF_scale_xy.pdf")),height = 3,width = 4)
}
```


### p0 函数调用
```{r}
# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/ABsub.txt")
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/ABsub.txt")

# pfun_Increasement_CDF(df,1026,"ABsub")
pfun_Increasement_CDF_scale_x(df,1026,"ABsub")
# pfun_Increasement_CDF_scale_xy(df,1026,"ABsub")

# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/010_sampleSize_to_variantsDiscovery/001_out_bySub/Dsub.txt")
df <- read_tsv("data/010_sampleSize2VariantsDiscovery/Dsub.txt")

# pfun_Increasement_CDF(df,850,"Dsub")
# pfun_Increasement_CDF_scale_x(df,850,"Dsub")
# pfun_Increasement_CDF_scale_xy(df,850,"Dsub")
``` 


#### ======================
## Allele age

### A:SNP count with allele age

```{r}
col <- c("#004680","#e69f00","#d5311c")
## probably damaging [0.957,1]
## possibly damaging [0.453,0.956]
dfAlleleAge <- dftotal %>% 
  filter(Region=="CDS",Variant_type=="SYNONYMOUS"|Variant_type=="NONSYNONYMOUS")%>%
  filter(!is.na(Ancestral),Ancestral == Ref | Ancestral == Alt) %>%
  mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>%
  filter(!is.na(AlleleAge_J), !is.na(AlleleAge_M), !is.na(AlleleAge_R)) %>%
  mutate(LoadGroup=factor(LoadGroup, levels = c("Synonymous","Nonsynonymous","Deleterious")))

df_db_summary <- dfAlleleAge %>% 
  group_by(LoadGroup,Sub) %>% 
  summarise(n=n())

q <- dfAlleleAge %>% 
  ggplot(aes(x=Sub, fill=LoadGroup))+
  geom_bar(position = position_dodge(width = 0.9))+
  geom_text(data = df_db_summary, aes(label=n, y=n), 
            position = position_dodge(width = 0.9), size=2.5)+
  scale_fill_manual(values = col)+
  labs(y="The num of SNP having allele age")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        # text = element_text(size = 14),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
```

### B: allele age distribution

```{r}
q <- dfAlleleAge %>% 
  mutate(Region=factor(Region, levels = c("UTR_5","Intron","CDS","UTR_3"))) %>% 
  ggplot(aes(x=AlleleAge_J, fill=LoadGroup))+
  geom_histogram(aes(y=..count..), color="white")+
  facet_grid(Sub~LoadGroup, scales = "free")+
  # scale_color_brewer(palette = "Set2")+
  # scale_fill_brewer(palette = "Set2")+
  scale_fill_manual(values = col)+
  scale_color_manual(values = col)+
  scale_x_continuous(trans = "log10",
                     breaks = scales::trans_breaks("log10", function(x) 10^x),
                     labels = scales::trans_format("log10",   scales::math_format(10^.x)))+
  labs(x="Mutation age (Years)")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        # axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 16),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
```

### C: allele age vs DAF

```{r}
q <- dfAlleleAge %>% 
  group_by(Sub,LoadGroup) %>% 
  slice_sample(prop = 0.2) %>%
  ggplot(aes(x=AlleleAge_J, y=DAF,color=LoadGroup, fill=LoadGroup))+
  geom_point(alpha=0.2)+
  facet_grid(Sub~LoadGroup)+
  scale_color_manual(values = col)+
  scale_fill_manual(values = col)+
  scale_x_continuous(trans = "log10",
                     # breaks = scales::trans_breaks("log10", function(x) 10^x),
                     labels = scales::trans_format("log10",   scales::math_format(10^.x)))+
  labs(x="Mutation age (Years)",y="Derived allele frequency")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        # axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 16),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
ggsave(plot = q, "~/Documents/temp.pdf", width = 6, height = 4)

```

### plot
```{r}
q <- dfAlleleAge %>% 
  group_by(Sub,LoadGroup) %>% 
  slice_sample(prop = 0.2) %>%
  ggplot(aes(x=DAF, y=AlleleAge_J,color=LoadGroup, fill=LoadGroup))+
  geom_point(alpha=0.2)+
  facet_grid(Sub~LoadGroup)+
  scale_color_manual(values = col)+
  scale_fill_manual(values = col)+
  scale_y_continuous(trans = "log10",
                     # breaks = scales::trans_breaks("log10", function(x) 10^x),
                     labels = scales::trans_format("log10",scales::math_format(10^.y)))+
  # labs(x="Mutation age (Years)",y="Derived allele frequency")+
  theme_classic()+
  theme(legend.title = element_blank(),
        legend.position = "none",
        legend.background = element_rect(fill="transparent"),
        strip.background = element_blank(),
        # axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 16),
        panel.border = element_rect(colour = "black",fill = "transparent"))
q
ggsave(plot = q, "~/Documents/temp.pdf", width = 6, height = 4)
```

### ** 密度图allele age vs DAF (2D plot)
- https://www.cnblogs.com/icydengyw/p/12309455.html
```{r}
p <- dfAlleleAge %>% 
  group_by(Sub,LoadGroup) %>% 
  slice_sample(prop = 0.5) %>%
  ggplot(aes(x=AlleleAge_J, y=DAF))+
  # stat_density_2d(aes(x = AlleleAge_J, y = DAF, fill = stat(nlevel)),
  #                 geom = "polygon", n = 500, bins = 10, contour = TRUE) +
  # geom_bin2d(aes(x = AlleleAge_J, y = DAF),bins = 10) +
  ### raster 方法
  stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
  scale_fill_distiller(palette=4, direction=1)+
  # scale_fill_fermenter(palette=4, direction=1)+
  # scale_fill_gradient()+
  ### hex 方法
  # geom_hex(bins = 50) +
  # scale_fill_gradient(
  #   low = "#9696F2", 
  #   high = "#0A0A3D")+

  # geom_point(alpha=0.5)+
  facet_grid(.~LoadGroup)+
  # scale_color_manual(values = col)+
  # scale_fill_manual(values = col)+
  # scale_fill_viridis_c(option = "A")+
  # scale_fill_continuous(type="viridis")+
  scale_x_continuous(trans = "log10",
                     expand = c(0, 0),
                     # breaks = scales::trans_breaks("log10", function(x) 10^x),
                     labels = scales::trans_format("log10",   scales::math_format(10^.x)))+
  scale_y_continuous(expand = c(0, 0))+
  labs(x="Derived allele age (Years)",y="Derived allele frequency")+
  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),
        panel.border = element_rect(colour = "black",fill = "transparent"))
p

ggsave(p, filename = "~/Documents/test.pdf", height = 2.22, width = 3.6)

```


### D: Del proportion on diff allele age 
```{r}
dabins <- c(0,1,2,4,8,16,32, 64, 128)
# ### 分亚基因组
# df <- dfAlleleAge %>% 
#   mutate(bin=cut(AlleleAge_J/1000,breaks=dabins)) %>%
#   filter(!is.na(LoadGroup), !is.na(bin)) %>% 
#   group_by(Sub,bin) %>% 
#   count(LoadGroup) %>% 
#   mutate(fre=n/sum(n)) 

# q <- df %>% 
#   # filter(Sub=="A") %>% 
#   ggplot(aes(x=bin, y=fre, color=LoadGroup,group=LoadGroup))+
#   geom_point()+
#   geom_line()+
#   facet_grid(.~Sub)+
#   scale_color_manual(values = col)+
#   scale_x_discrete(label=dabins)+
#   labs(x="Derived allele age (Ka)",y="Deleterious SNPs (%)") +
#   theme_classic()+
#   theme(legend.title = element_blank(),
#         legend.position = "none",
#         legend.background = element_rect(fill="transparent"),
#         strip.background = element_blank(),
#         axis.line = element_blank(),
#         panel.border = element_rect(colour = "black",fill = "transparent"))
# q
# 
# ggsave(plot = q, "~/Documents/temp.pdf", width = 5, height = 3)

### 不分亚基因组

df <- dfAlleleAge %>% 
  mutate(bin=cut(AlleleAge_J/1000,breaks=dabins)) %>%
  filter(!is.na(LoadGroup), !is.na(bin)) %>% 
  group_by(bin) %>% 
  count(LoadGroup) %>% 
  mutate(fre=n/sum(n)) 

p <- df %>% 
  ggplot(aes(x=bin, y=fre, color=LoadGroup,group=LoadGroup))+
  geom_point(size=3)+
  geom_line(size=1)+
  scale_color_manual(values = col)+
  scale_x_discrete(label=dabins)+
  # labs(x="Derived allele age (Ka)",y="Deleterious SNPs (%)") +
  labs(x="Derived allele age (Ka)",y="Proportion of SNPs (%)") +
  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(plot = p, "~/Documents/test.pdf", height = 2.22,width = 1.8)
# ggsave(plot = p, "~/Documents/test.pdf", height = 2.9,width = 2.4)
ggsave(plot = p, "~/Documents/test.pdf", height = 2.38,width = 3.59)

```

### E: barplot: allele age in different types
```{r}

dfave <- dfAlleleAge %>% 
  select(Sub,LoadGroup,AlleleAge_J) %>% 
  group_by(LoadGroup,Sub) %>% 
  summarise(ave=mean(AlleleAge_J),sd=sd(AlleleAge_J),N=n(),se=sd/sqrt(N))

sum <- sum(dfave$N)

dfsum <- dfave %>% 
  group_by(LoadGroup) %>% 
  summarise(sum=sum(N))

#### 合计有 22220个位点能够计算出时间的信息，
#### 其中syn有10963个位点, nonsyn有9622个位点, del有1635个位点。

p <- dfave %>% 
  ggplot(aes(x=Sub,y=ave/1000,fill=LoadGroup))+
  geom_bar(stat = "identity",position = position_dodge(),alpha=0.9)+
  geom_errorbar(aes(ymin=(ave-se)/1000,ymax=(ave+se)/1000),width=0.2,color="black",
                position = position_dodge(.9))+
  scale_fill_manual(values = c("#004680","#e69f00","#d5311c"))+
  labs(x="", y="Average mutation age (x1,000)")+
  theme_classic()
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.4),
    # legend.position = 'none',
    legend.title = 'Variant type',
    text = element_text(size = 10))


p 

ggsave(plot = p, "~/Documents/test.pdf", height = 3,width = 5)

```

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

df <- dfAlleleAge %>% 
  select(Sub,LoadGroup,AlleleAge_J) %>% 
  mutate(LoadGroup = factor(LoadGroup))

df$LoadGroup <- factor(df$LoadGroup,levels = c("Deleterious","Nonsynonymous","Synonymous"))

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



#### =======================
## 1.Gene expression vs strongly burden
```{r}
df_burden1 <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")
df_expression <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneExpression.averageTPM.txt.gz")

df_burden1_expression <- df_burden1 %>% 
  left_join(y = df_expression,by = c("geneName"="Gene")) %>% 
  filter(!is.na(Subgenome))

temp <- df_burden1_expression %>% 
  # filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  group_by(SlightlyOrStrongly, Subgenome,bin) %>% 
  count() %>%   
  filter(n > 20, !is.na(bin))

colN1 <- brewer.pal(9,"Reds")[c(7)]
colA <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(13)
colC <- rev(colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(8))
colN2 <- brewer.pal(9,"Reds")[c(9)]

colB <- c(colN1,colA,colC,colN2)
# colB <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(40)
# xLabel <- expression(2^{-12},"",2^{-8},"",2^{-4},"",2^{0},"",2^{4},"",2^{8},"")
# xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})
xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})
p <- df_burden1_expression %>% 
  filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  filter(bin %in% temp$bin) %>%
  ggplot(aes(x=bin, y=mean_burden)) +
# stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
# geom_boxplot(aes(fill=bin),outlier.shape = NA) +
# geom_point(position = position_jitter(width = 0.2), size=0.1)+
# stat_summary(fun=median, geom="point", color="blue",
#              position = position_dodge(1),
#              size=0.1)+
  stat_summary(aes(color=bin),fun.data = mean_se,size=0.2)+
  scale_fill_manual(values = colB)+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  labs(x="TPM", y= "Gene burden")+
  # scale_y_continuous(trans = "log2",
  #                    breaks = trans_breaks("log2",function(x) 2^x),
  #                    labels = trans_format("log2",math_format(2^.x)))+
  scale_x_discrete(breaks=c("(-12,-11]","(-9,-8]","(-6,-5]","(-3,-2]","(0,1]","(3,4]","(6,7]","(9,10]"),
                   labels=expression(2^{-12},2^{-9},2^{-6},2^{-3},2^{0},2^{3},2^{6},2^{9}))+
  # scale_x_discrete(labels=xLabel)+
  scale_color_manual(values = colB)+
  coord_cartesian(ylim = c(0,0.04))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7),
    legend.position = 'none',
    text = element_text(size = 14))

p
ggsave(plot = p, "~/Documents//Burden1_vs_TPM.pdf",height = 2,width = 2)
ggsave(p,filename = "~/Documents/p.pdf",height = 6.33,width = 7.82,units = "cm")

```


## 1.去除异常值：Gene expression vs strongly burden
```{r}
df_burden1 <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")
df_expression <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneExpression.averageTPM.txt.gz")

df_burden1_expression <- df_burden1 %>% 
  left_join(y = df_expression,by = c("geneName"="Gene")) %>% 
  filter(!is.na(Subgenome))

temp <- df_burden1_expression %>% 
  # filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  group_by(SlightlyOrStrongly, Subgenome,bin) %>% 
  count() %>%   
  filter(n > 20, !is.na(bin))

# colN1 <- brewer.pal(9,"Reds")[c(7)]
colA <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(13)
colC <- rev(colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(8))
# colN2 <- brewer.pal(9,"Reds")[c(9)]

colB <- c(colN1,colA,colC,colN2)
xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})
p <- df_burden1_expression %>% 
  filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  filter(bin %in% temp$bin) %>%
  ggplot(aes(x=bin, y=mean_burden)) +
# stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
# geom_boxplot(aes(fill=bin),outlier.shape = NA) +
# geom_point(position = position_jitter(width = 0.2), size=0.1)+
# stat_summary(fun=median, geom="point", color="blue",
#              position = position_dodge(1),
#              size=0.1)+
  stat_summary(aes(color=bin),fun.data = mean_se,size=0.2)+
  scale_fill_manual(values = colB)+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  labs(x="TPM", y= "Gene burden")+
  # scale_y_continuous(trans = "log2",
  #                    breaks = trans_breaks("log2",function(x) 2^x),
  #                    labels = trans_format("log2",math_format(2^.x)))+
  scale_x_discrete(breaks=c("(-12,-11]","(-9,-8]","(-6,-5]","(-3,-2]","(0,1]","(3,4]","(6,7]","(9,10]"),
                   labels=expression(2^{-12},2^{-9},2^{-6},2^{-3},2^{0},2^{3},2^{6},2^{9}))+
  # scale_x_discrete(labels=xLabel)+
  scale_color_manual(values = colB)+
  coord_cartesian(ylim = c(0,0.02))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7),
    legend.position = 'none',
    text = element_text(size = 10.38))

p
# ggsave(plot = p, "~/Documents/Burden1_vs_TPM.pdf",height = 2,width = 2)
ggsave(plot = p, "~/Documents/Burden1_vs_TPM.pdf",height = 2.14,width = 2.35)

```

### 每个 bin 的基因数量
```{r}
df_burden1 <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")
df_expression <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneExpression.averageTPM.txt.gz")

df_burden1_expression <- df_burden1 %>% 
  left_join(y = df_expression,by = c("geneName"="Gene")) %>% 
  filter(!is.na(Subgenome)) %>% 
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0))

temp <- df_burden1_expression %>% 
  # filter(SlightlyOrStrongly=="StronglyBurden") %>%
  # group_by(SlightlyOrStrongly, Subgenome,bin) %>% 
  group_by(SlightlyOrStrongly,bin) %>% 
  count() %>%   
  filter(n > 20, !is.na(bin))

colN1 <- brewer.pal(9,"Reds")[c(7)]
colA <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(13)
colC <- rev(colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(8))
colN2 <- brewer.pal(9,"Reds")[c(9)]
colB <- c(colN1,colA,colC,colN2)

xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})

p <- temp %>% 
  filter(SlightlyOrStrongly=="StronglyBurden") %>%
  ggplot(aes(x=bin, y=n)) +
  geom_bar(stat = 'identity',position = "dodge")+
  theme_classic()
  # theme_void()
  
p
ggsave(plot = p, "~/Documents//Burden1_vs_TPM.pdf",height = 2,width = 2)
ggsave(p,filename = "~/Documents/p.pdf",height = 6.33,width = 7.82,units = "cm")
```


## 1.2 新加内容：哪些功能类型的基因突变负荷较低？哪些功能类型的基因突变负荷较高？
答：从GO富集分析的结果来看，
```{r}
df_burden1 <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")
df_expression <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneExpression.averageTPM.txt.gz")

df_burden1_expression <- df_burden1 %>% 
  left_join(y = df_expression,by = c("geneName"="Gene")) %>% 
  filter(!is.na(Subgenome))

temp <- df_burden1_expression %>% 
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  group_by(SlightlyOrStrongly, Subgenome,bin) %>% 
  count() %>%   
  filter(n > 20, !is.na(bin))

df_geneWithHighExpression <- df_burden1_expression %>% 
  filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  filter(bin %in% temp$bin) %>% 
  ### 当TPM_mean_log > 6时，检测到 1506 个基因是高表达的
  ### 当TPM_mean_log > 8时，检测到 203 个基因是高表达的
  ### 当TPM_mean_log > 7且小于9时，检测到 552 个基因是高表达的
  ### 当TPM_mean_log > 9时，检测到 60 个基因是高表达的，这60个基因的负荷有点高，需要去掉的
  # filter(TPM_mean_log > 8,TPM_mean_log < 9) %>%
  # filter(TPM_mean_log > 7, TPM_mean_log < 9) %>%
  filter(TPM_mean_log > 0, TPM_mean_log < 1) %>%
  select(geneName)

```


### p
```{r}
##第一步：准备基因列表文件
### 库文件df2,获取高表达基因的最长转录本
df2 <- read_tsv("data/006_Fig2/wheat_v1.1_nonoverlap_addPos.txt") %>%
  select(Gene,Transcript=Longest_transcript)
  
### gene -> transcript -> version01G
df3 <- df_geneWithHighExpression %>% 
  left_join(df2,by=c("geneName" = "Gene")) %>% 
  select(Transcript) %>% 
  mutate(Transcript_version01=str_replace(Transcript,"02G","01G")) %>% 
  select(Transcript_version01)

### get gene list and write out
outfileHighExpress <- "data/006_Fig2/geneList/geneList_withHighExpress.txt"
write_tsv(df3,outfileHighExpress,col_names = F)

##第二步：开始GO分析画图
library("clusterProfiler")
library(RColorBrewer)
library(tidyverse)
### ①导入基因列表
gene <- read.table(outfileHighExpress,header = F,sep="\t")
gene <- as.factor(gene$V1)
# head(gene)

### ②导入注释文件
term2gene <- read.table("/Users/Aoyue/project/wheatVMapII/003_dataAnalysis/005_vcf/019_popGen/104_XPCLR/007_GO/001_input/TERM2GENE_v1__HCgenes_v1.0_repr.TEcleaned.txt",header=F,sep="\t")
term2name <- read.table("/Users/Aoyue/project/wheatVMapII/003_dataAnalysis/005_vcf/019_popGen/104_XPCLR/007_GO/001_input/TERM2name_v1__HCgenes_v1.0_repr.TEcleaned.txt",header=F,sep="\t")
head(term2gene)
head(term2name)

### ③富集分析
x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.05)
# x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 1, pAdjustMethod = "BH", qvalueCutoff = 1)

### ④绘制条形图-软件包自带的，不能省略这一步。一定要看一下结果，再做自己的图
barplot(x)+ scale_fill_gradient(low="#7acdc9",high="#d59b4f")
barplot(x,showCategory = 15)
dotplot(x)+ scale_color_gradient(low="#7acdc9",high="#d59b4f")

### ⑤输出结果，自行画图
# library(topGO)
# plotGOgraph(x)
# 设置结果输出文件
ouf <- c("data/006_Fig2/geneList/GO/001_GO_geneList_withHighExpress.txt")
write.csv(x,str_replace(ouf,".txt",".csv"),row.names =FALSE)
write_tsv(as_tibble(x),path = ouf)

### ⑥自定义画图，可自行美化及分析
dataFile <- ouf
df <- read_tsv(dataFile) %>%
  arrange(p.adjust) ###重点，这里按照显著性排序的！！！

df$Description <- factor(df$Description,levels = df$Description)
########################## 自定义画图 气泡图
# 
# p <- df %>% 
#   filter(pvalue<0.04) %>% 
#   ### 由于显著富集的通路很多，取前20展示
#   top_n(20) %>% 
#   ggplot(aes(x=GeneRatio,y= Description,col=p.adjust,size=Count))+
#   geom_point()+theme_bw()+ylab("")+
#   # facet_grid(.~Ontology)+
#   scale_color_gradient(low="#7acdc9",high="#d59b4f")
# 
# p
# 
# ggsave(p,filename = "~/Documents/GO_dotplot.pdf",height =4,width = 5)
# ggsave(p,filename = "~/Documents/GO_dotplot.png",height = 6.33,width = 7.82,units = "cm")

########################## 自定义画图 柱状图
#### 这个图画出来以后，要从下往上看，最下面的是最重要的，一定要注意！！！
p <- df %>% 
  filter(pvalue < 0.04) %>% 
  ### 由于显著富集的通路很多，取前20展示
  slice_head(n = 15) %>% 
  ggplot(aes(x=Description,y=Count,fill=p.adjust))+
  geom_bar(stat = 'identity')+theme_bw()+xlab("")+
  scale_fill_gradient(low="#7acdc9",high="#d59b4f")+
  # facet_grid(.~Sub)+
  coord_flip()+
  theme(text = element_text(size = 14))

p

ggsave(p,filename = "~/Documents/GO_barplot.pdf",height =4,width = 7)
ggsave(p,filename = "~/Documents/GO_barplot.png",height = 10.16,width =17.78,units = "cm")

### 第三步：想分清楚是 BP CC MF 中具体哪个过程？
############################## 加入 Ontology ####################################
df_GOdb <- read_tsv("data/009_GO/000_DB/GDdb.txt") %>% 
  distinct(GOID,.keep_all = T) %>% 
  select(GOID,Description_byGOdb=Description,Ontology)

df_BP_CC_MF <- df %>% 
  left_join(df_GOdb,by=c("ID"="GOID")) %>%
    arrange(p.adjust)


########################## 柱状图
p <- df_BP_CC_MF %>%
  filter(pvalue < 0.04) %>%
  ### 由于显著富集的通路很多，取前20展示
  slice_head(n = 15) %>% 
  ggplot(aes(x=Description,y=Count,fill=p.adjust))+
  geom_bar(stat = 'identity')+theme_bw()+xlab("")+
  scale_fill_gradient(low="#7acdc9",high="#d59b4f")+
  facet_grid(.~ Ontology)+
  coord_flip()+
  theme(text = element_text(size = 14))

p

## save height 4 width 7.5
ggsave(p,filename = "~/Documents/GO_barplot_ont.pdf",height =4,width = 7)
ggsave(p,filename = "~/Documents/GO_barplot_ont.png",height = 10.16,width =18.59,units = "cm")

######################### 气泡图

# p <- df_BP_CC_MF %>%
#   ### 由于显著富集的通路很多，取前20展示
#   top_n(15) %>%
#   ggplot(aes(x=GeneRatio,y= Description,col=p.adjust,size=Count))+
#   geom_point()+theme_bw()+ylab("")+
#   scale_color_gradient(low="#7acdc9",high="#d59b4f")+
#   facet_grid(.~Ontology)
# 
# p
# 
# ggsave(p,filename = "~/Documents/GO_dotplot_ont.png",height = 10.16,width =17.78,units = "cm")

```

## P相关性系数
```{r}
df_burden1 <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")
df_expression <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneExpression.averageTPM.txt.gz")

df_burden1_expression <- df_burden1 %>% 
  left_join(y = df_expression,by = c("geneName"="Gene")) %>% 
  filter(!is.na(Subgenome))

temp <- df_burden1_expression %>% 
  # filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  group_by(SlightlyOrStrongly, Subgenome,bin) %>% 
  count() %>%   
  filter(n > 20, !is.na(bin))

colN1 <- brewer.pal(9,"Reds")[c(7)]
colA <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(13)
colC <- rev(colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(8))
colN2 <- brewer.pal(9,"Reds")[c(9)]

colB <- c(colN1,colA,colC,colN2)
# colB <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(40)
# xLabel <- expression(2^{-12},"",2^{-8},"",2^{-4},"",2^{0},"",2^{4},"",2^{8},"")
# xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})
xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})
p <- df_burden1_expression %>% 
  filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>%
  filter(TPM_mean_log > 4) %>% 
  filter(mean_burden < 0.5) %>% 
  # mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  # filter(bin %in% temp$bin) %>%
  ggplot(aes(x=TPM_mean_log, y=mean_burden)) +
# stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
# geom_boxplot(aes(fill=bin),outlier.shape = NA) +
# geom_point(position = position_jitter(width = 0.2), size=0.1)+
# stat_summary(fun=median, geom="point", color="blue",
#              position = position_dodge(1),
#              size=0.1)+
  # stat_summary(aes(color=bin),fun.data = mean_se,size=0.2)+
  geom_point(aplha=0.3,size=2,color="red")+
  stat_cor(
    aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
    label.x.npc = 0.3,label.y.npc = 0.01,size= 2)+
  # scale_fill_manual(values = colB)+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  labs(x="TPM", y= "Gene burden")+
  # scale_y_continuous(trans = "log2",
  #                    breaks = trans_breaks("log2",function(x) 2^x),
  #                    labels = trans_format("log2",math_format(2^.x)))+
  # scale_x_discrete(breaks=c("(-12,-11]","(-9,-8]","(-6,-5]","(-3,-2]","(0,1]","(3,4]","(6,7]","(9,10]"),
                   # labels=expression(2^{-12},2^{-9},2^{-6},2^{-3},2^{0},2^{3},2^{6},2^{9}))+
  # scale_x_discrete(labels=xLabel)+
  # scale_color_manual(values = colB)+
  coord_cartesian(ylim = c(0,0.04))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7),
    legend.position = 'none',
    text = element_text(size = 9))

p
ggsave(plot = p, "~/Documents//Burden1_vs_TPM.pdf",height = 2,width = 2)


```

```{r}
dff <- df_burden1_expression %>% 
  filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  filter(!is.na(mean_burden), is.finite(mean_burden),mean_burden!=0)
```


## 2.Gene expression vs slightly burden
```{r}
df_burden1 <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")
df_expression <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneExpression.averageTPM.txt.gz")

df_burden1_expression <- df_burden1 %>% 
  left_join(y = df_expression,by = c("geneName"="Gene")) %>% 
  filter(!is.na(Subgenome))

temp <- df_burden1_expression %>% 
  # filter(SlightlyOrStrongly=="StronglyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  group_by(SlightlyOrStrongly, Subgenome,bin) %>% 
  count() %>%   
  filter(n > 20, !is.na(bin))

colN1 <- brewer.pal(9,"Reds")[c(7)]
colA <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(13)
colC <- rev(colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(8))
colN2 <- brewer.pal(9,"Reds")[c(9)]

colB <- c(colN1,colA,colC,colN2)
# colB <- colorRampPalette(brewer.pal(9,"Reds")[c(2,4,5,6,7,8,9)])(40)
# xLabel <- expression(2^{-12},"",2^{-8},"",2^{-4},"",2^{0},"",2^{4},"",2^{8},"")
# xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})
xLabel <- expression("",2^{-12},"","",2^{-9},"","",2^{-6},"","",2^{-3},"","",2^{0},"","",2^{3},"","",2^{6},"","",2^{9})
p <- df_burden1_expression %>% 
  filter(SlightlyOrStrongly=="SlightlyBurden") %>%
  mutate(TPM_mean_log = log2(TPM_mean)) %>% 
  mutate(bin=cut_width(TPM_mean_log,width = 1,boundary = 0)) %>% 
  filter(bin %in% temp$bin) %>%
  ggplot(aes(x=bin, y=mean_burden)) +
# stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
# geom_boxplot(aes(fill=bin),outlier.shape = NA) +
# geom_point(position = position_jitter(width = 0.2), size=0.1)+
# stat_summary(fun=median, geom="point", color="blue",
#              position = position_dodge(1),
#              size=0.1)+
  stat_summary(aes(color=bin),fun.data = mean_se,size=0.2)+
  scale_fill_manual(values = colB)+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  labs(x="TPM", y= "Slightly gene burden")+
  # scale_y_continuous(trans = "log2",
  #                    breaks = trans_breaks("log2",function(x) 2^x),
  #                    labels = trans_format("log2",math_format(2^.x)))+
  scale_x_discrete(breaks=c("(-12,-11]","(-9,-8]","(-6,-5]","(-3,-2]","(0,1]","(3,4]","(6,7]","(9,10]"),
                   labels=expression(2^{-12},2^{-9},2^{-6},2^{-3},2^{0},2^{3},2^{6},2^{9}))+
  scale_x_discrete(labels=xLabel)+
  scale_color_manual(values = colB)+
  # coord_cartesian(ylim = c(0,0.04))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7),
    legend.position = 'none',
    text = element_text(size = 9))

p
ggsave(plot = p, "~/Documents//Burden1_vs_TPM.pdf",height = 2,width = 2)

```


## 3.Tissue expression VS mutation burden
- caption: Deleterious mutation of genes expressed across tissues.
conclusion: Deleterious mutations are depleted in genes expressed in grain and anther.

```{r message=FALSE, warning=FALSE}
### Step1: 将29万个基因在35个组织中的表达量过滤，只保留高置信度的最长转录本;将Load的结果也合并到HC数据库中
### 第一个数据框
df_expression <- read_tsv("data/011_geneExpression/expressionByTissue.txt.gz")

### 第二个数据框
dataFile <- c("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")

df_load <- read_tsv(dataFile) %>%
  rename(Gene = geneName)

dataFile <- c("data/011_geneExpression/geneHC.txt.gz")
df_geneHC <- read_tsv(dataFile) %>%
  rename(Transcripts = Longest_transcript) %>%
  left_join(.,df_expression,by="Transcripts")  ###合并第1个数据框

## expression 文件与load文件合并
  df_expression_load <- df_load %>%
    left_join(.,df_geneHC,by="Gene")  ###合并第2个数据框

# write.table(df_geneHC,file="~/Documents/test.txt",quote=F,col.name=T,row.names=F,sep = "\t")
  
  
  ### Step2： 转换表格，将组织表达 -> 长表
df_expression_load_long <- df_expression_load %>%
  pivot_longer(.,cols = "Internode #2":"third leaf sheath", names_to = "Tissue",values_to = "TPM") %>%
  filter(TPM < 0.5)

### 计算load和组织的关系
dfstatistic <- df_expression_load_long %>%
  # pivot_longer(.,cols = FreSyn:RatioHGDelVSSyn,names_to = "LoadType",values_to = "LoadValue") %>%
  filter(!is.na(mean_burden), is.finite(mean_burden)) %>%
  group_by(SlightlyOrStrongly,Tissue) %>% ### 计算不同分类下，不同组织load的平均值 标准差
  summarise(ave = mean(mean_burden), sd = sd(mean_burden), n= n(), se = sd/sqrt(n))

# write.table(dfstatistic,file="~/Documents/barplot.txt",quote=F,col.name=T,row.names=F,sep = "\t")

```

接上- 开始画图
### strongly burden by individual
```{r}
### 将多个组织的命名方式都添加进去，便于排序画图使用
dataFile <- c("data/011_geneExpression/TissueHashMap.txt")
dfTissueHashMap <- read_tsv(dataFile) %>%
  rename(Tissue = Tissue_inOri)

AoCase <- "StronglyBurden"

df <- dfstatistic %>%
  left_join(.,dfTissueHashMap,by="Tissue") %>%
  filter(SlightlyOrStrongly == AoCase) %>%
  mutate(High_level_tissue=factor(High_level_tissue,levels = c("roots","leaves/shoots","spike","grain"))) %>%
  arrange(High_level_tissue,-ave) ## 注意排序！！！

### 画图的时候根据数据框的顺序画出来
### 35 个组织的顺序定义
df$Tissue <- factor(df$Tissue,levels = df$Tissue)

colB <- c("#9a784c","#749b64","#e6b700","#b3b563") ## 4个组织的颜色

p <- df %>%
  ggplot(aes(x= Tissue, y=ave,fill = High_level_tissue)) +
  geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=ave-se, ymax=ave+se), width=.1,
                # color="gray",
                 position=position_dodge(.9))+
  # labs(x="",y= "Strongly burden in CDS")+
  labs(x="",y= "Mutation burden")+
  geom_hline(yintercept = 0.0095, color = "gray",linetype="dashed")+
  theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  scale_fill_manual(values = colB)+
  theme(plot.margin=unit(rep(1,4),'lines'))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7)
    ,legend.position = 'none'
    ,text = element_text(size = 8))

p

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

### 查看每个组织有多少个基因含有burden值
q <- df %>%
  ggplot(aes(x= Tissue, y=n,fill = High_level_tissue)) +
   geom_bar(stat="identity", position=position_dodge()) +
  labs(x="",y= "Gene number with strongly burden")+
  geom_hline(yintercept = 0.0318, color = "gray",linetype="dashed")+
  theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  scale_fill_manual(values = colB)+
  theme(plot.margin=unit(rep(1,4),'lines'))+
  theme(
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7)
    ,legend.position = 'none'
    ,text = element_text(size = 12))

q

```

### slightly burden by individual
```{r}
### 将多个组织的命名方式都添加进去，便于排序画图使用
dataFile <- c("data/011_geneExpression/TissueHashMap.txt")
dfTissueHashMap <- read.delim(dataFile) %>%
  rename(Tissue = Tissue_inOri)

AoCase <- "SlightlyBurden"

df <- dfstatistic %>%
  left_join(.,dfTissueHashMap,by="Tissue") %>%
  filter(SlightlyOrStrongly == AoCase) %>%
  mutate(High_level_tissue=factor(High_level_tissue,levels = c("roots","leaves/shoots","spike","grain"))) %>%
  arrange(High_level_tissue,-ave) ## 注意排序！！！

### 画图的时候根据数据框的顺序画出来
### 35 个组织的顺序定义
df$Tissue <- factor(df$Tissue,levels = df$Tissue)

colB <- c("#9a784c","#749b64","#e6b700","#b3b563") ## 4个组织的颜色

p <- df %>%
  ggplot(aes(x= Tissue, y=ave,fill = High_level_tissue)) +
   geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=ave-se, ymax=ave+se), width=.2,
                 position=position_dodge(.9))+
  labs(x="",y= "Slightly burden in CDS")+
  geom_hline(yintercept = 0.14, color = "gray",linetype="dashed")+
  theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  scale_fill_manual(values = colB)+
  theme(plot.margin=unit(rep(1,4),'lines'))+
  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/slightlyBurden_in_Tissue.pdf",height = 5,width = 7)

### 查看每个组织有多少个基因含有burden值
p <- df %>%
  ggplot(aes(x= Tissue, y=n,fill = High_level_tissue)) +
   geom_bar(stat="identity", position=position_dodge()) +
  labs(x="",y= "Gene number with strongly burden")+
  geom_hline(yintercept = 0.0318, color = "gray",linetype="dashed")+
  theme(axis.text.x = element_text(vjust = 1, hjust = 1, angle = 45))+
  scale_fill_manual(values = colB)+
  theme(plot.margin=unit(rep(1,4),'lines'))+
  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


```

## 4.Number of tissues vs mutation burden
```{r}
### 首先找到每个基因组织表达的宽度
df_tissueBreath <- read_tsv("data/011_geneExpression/002_tissueBreath/002_Azhurnaya_tissue_breadth_100Kgenet_1TPM_miniVersion.txt")

df_tissueBreath %>% 
  ggplot(aes(x=ExpressionBreadth))+
  geom_histogram()+
  theme_classic()

### 其次找到每个基因的load值
df_burden1 <- read_tsv("data/fromDaxing/Burden_Expression/004_burdenMean_vs_expression/GeneBurden.averageIndividualBurden.txt.gz")

df_burden1 %>% 
  ggplot(aes(x= log10(mean_burden)))+
  facet_wrap(~SlightlyOrStrongly)+
  geom_histogram()+
  theme_classic()

### 再次将load和组织表达宽度进行相关性分析
df_burden1_tissueBreath <- df_burden1 %>% 
  left_join(y = df_tissueBreath,by = c("geneName"="Gene"))

temp <- df_burden1_tissueBreath %>% 
  ###手动分割bin
  # mutate(bin=cut(ExpressionBreadth,breaks = c(0,2,4,6,10,15,20,23))) %>% 
  mutate(bin=cut_width(ExpressionBreadth,width = 2,boundary = 0)) %>%
  filter(!is.na(bin)) %>% 
  group_by(SlightlyOrStrongly,bin) %>% 
  count() 

### 颜色配置
colA <- colorRampPalette(brewer.pal(9,"Reds")[c(4,5,6,7,8,9)])(6)
colC <- rev(colorRampPalette(brewer.pal(9,"Reds")[c(4,5,6,7,8,9)])(5))
colB <- c(colA,colC)

p <- df_burden1_tissueBreath %>% 
  filter(SlightlyOrStrongly=="StronglyBurden") %>%
  # filter(mean_burden!=0,is.finite(mean_burden)) %>%
  # filter(is.finite(mean_burden)) %>
  # mutate(bin=cut(ExpressionBreadth,breaks = c(0,2,4,6,10,15,20,23))) %>% 
  mutate(bin=cut_width(ExpressionBreadth,width = 2,boundary = 0)) %>%
  filter(!is.na(bin)) %>% 
  ggplot(aes(x=bin, y=mean_burden,color=bin)) +
  stat_summary(fun.data = mean_se,size=0.8)+
  scale_fill_manual(values = colB)+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  labs(x="Tissue breadth of gene expression", y= "Strongly mutation load per gene")+
  scale_x_discrete(breaks=
                     c("[0,2]","(6,8]","(12,14]","(18,20]"),
                   labels=c("2","8","14","20"))+
  scale_color_manual(values = colB)+
  # coord_cartesian(ylim = c(0,0.09))+
  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(plot = p, "~/Documents/strongly_vs_tissueBreath.pdf",height = 3,width = 4)


### 颜色配置需要修改一下，根据变化的趋势
colC <- rev(colorRampPalette(brewer.pal(9,"Reds")[c(4,5,6,7,8,9)])(11))
colB <- c(colC)
p <- df_burden1_tissueBreath %>% 
  filter(SlightlyOrStrongly=="SlightlyBurden") %>%
  # filter(mean_burden!=0,is.finite(mean_burden)) %>%
  # filter(is.finite(mean_burden)) %>
  # mutate(bin=cut(ExpressionBreadth,breaks = c(0,2,4,6,10,15,20,23))) %>% 
  mutate(bin=cut_width(ExpressionBreadth,width = 2,boundary = 0)) %>%
  filter(!is.na(bin)) %>% 
  ggplot(aes(x=bin, y=mean_burden,color=bin)) +
  stat_summary(fun.data = mean_se,size=0.8)+
  scale_fill_manual(values = colB)+
  # facet_grid(SlightlyOrStrongly~.,scales = "free")+
  labs(x="Tissue breadth of gene expression", y= "Slightly mutation load per gene")+
  scale_x_discrete(breaks=
                     c("[0,2]","(6,8]","(12,14]","(18,20]"),
                   labels=c("2","8","14","20"))+
  scale_color_manual(values = colB)+
  # coord_cartesian(ylim = c(0,0.09))+
  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(plot = p, "~/Documents/slightly_vs_tissueBreath.pdf",height = 3,width = 4)



# p <- df_burden1_tissueBreath %>% 
#   filter(SlightlyOrStrongly=="StronglyBurden") %>%
#   filter(mean_burden!=0,is.finite(mean_burden)) %>%
#   # filter(is.finite(mean_burden)) %>% 
#   # mutate(bin=cut(ExpressionBreadth,breaks = c(0,2,4,6,10,15,20,23))) %>% 
#   mutate(bin=cut_width(ExpressionBreadth,width = 2,boundary = 0)) %>%
#   filter(!is.na(bin)) %>% 
#   ggplot(aes(x=bin, y=mean_burden,color=bin)) +
# # stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
# # geom_boxplot(aes(fill=bin),outlier.shape = NA) +
# # geom_point(position = position_jitter(width = 0.2), size=0.1)+
# # stat_summary(fun=median, geom="point", color="blue",
# #              position = position_dodge(1),
# #              size=0.1)+
#   # stat_summary(fun.data = mean_se,size=0.2)+
#   geom_violin(aes(fill=bin),position = position_dodge(0.8)) +
#   stat_summary(fun=mean, geom="point", color="black",position = position_dodge(0.8),size=0.5)+
#   scale_fill_manual(values = colB)+
#   # facet_grid(SlightlyOrStrongly~.,scales = "free")+
#   labs(x="Tissue breadth of gene expression", y= "Strongly mutation load per gene")+
#   scale_x_discrete(breaks=
#                      c("[0,2]","(6,8]","(12,14]","(18,20]"),
#                    labels=c("2","8","14","20"))+
#   scale_color_manual(values = colB)+
#   coord_cartesian(ylim = c(0,0.09))+
#   theme(
#     panel.background = element_blank(),
#     panel.border = element_rect(colour = "black", fill=NA, size=0.7),
#     legend.position = 'none',
#     text = element_text(size = 9))
# 
# p
# ggsave(plot = p, "~/Documents/Burden1_vs_TPM.pdf",height = 2,width = 2)


# ### 过滤方式不一样
# p <- df_burden1_tissueBreath %>% 
#   filter(SlightlyOrStrongly=="StronglyBurden") %>%
#   # filter(mean_burden!=0,is.finite(mean_burden)) %>%
#   filter(is.finite(mean_burden)) %>%
#   # mutate(bin=cut(ExpressionBreadth,breaks = c(0,2,4,6,10,15,20,23))) %>% 
#   mutate(bin=cut_width(ExpressionBreadth,width = 2,boundary = 0)) %>%
#   filter(!is.na(bin)) %>% 
#   ggplot(aes(x=bin, y=mean_burden,color=bin)) +
# # stat_boxplot(geom = "errorbar",width=0.05,position = position_dodge(1))+ # add error bar
# # geom_boxplot(aes(fill=bin),outlier.shape = NA) +
# # geom_point(position = position_jitter(width = 0.2), size=0.1)+
# # stat_summary(fun=median, geom="point", color="blue",
# #              position = position_dodge(1),
# #              size=0.1)+
#   # stat_summary(fun.data = mean_se,size=0.2)+
#   geom_violin(aes(fill=bin),position = position_dodge(0.8)) +
#   stat_summary(fun=mean, geom="point", color="black",position = position_dodge(0.8),size=0.5)+
#   scale_fill_manual(values = colB)+
#   # facet_grid(SlightlyOrStrongly~.,scales = "free")+
#   labs(x="Tissue breadth of gene expression", y= "Strongly mutation load per gene")+
#   scale_x_discrete(breaks=
#                      c("[0,2]","(6,8]","(12,14]","(18,20]"),
#                    labels=c("2","8","14","20"))+
#   scale_color_manual(values = colB)+
#   coord_cartesian(ylim = c(0,0.09))+
#   theme(
#     panel.background = element_blank(),
#     panel.border = element_rect(colour = "black", fill=NA, size=0.7),
#     legend.position = 'none',
#     text = element_text(size = 9))
# 
# p
# ggsave(plot = p, "~/Documents/Burden1_vs_TPM2.pdf",height = 2,width = 2)

```

#### ========================
## del ratio whole genome scan
```{r}
### step1:根据原来数据框，求出 2方面值，1：单位CDS的结果，2：有害突变和同义突变的比
# df <- read_tsv("data/006_Fig2/genomeScan/002_delVSsynOnChr_10000000window1000000step_addEffectiveCDSLength_addRecom.txt")

# df <- read_tsv("/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/009_genomeScan_delvcSyn/002/001_delVSsynOnChr_20000000window5000000step_addEffectiveCDSLength.txt")

# df1 <- df %>%
#   mutate(SynFre=.$a001_synonymous/.$CDSLength,
#          NonsynFre=.$a002_nonsynonymous/.$CDSLength,
#          DelFre=.$a003_nonsynGERPandDerivedPPH2/.$CDSLength,
#          PPH2Fre=.$a004_nonsynDerivedPPH2/.$CDSLength,
#          GERPFre=.$a005_GERP_16way/.$CDSLength,
#          StopGainFre=.$a006_StopGain/.$CDSLength,
#          ### ratio
#          NonsynSynRatio=.$a002_nonsynonymous/.$a001_synonymous,
#          DelSynRatio=.$a003_nonsynGERPandDerivedPPH2/.$a001_synonymous,
#          PPH2SynRatio=.$a004_nonsynDerivedPPH2/.$a001_synonymous,
#          GERPSynRatio=.$a005_GERP_16way/.$a001_synonymous,
#          StopGainSynRatio=.$a006_StopGain/.$a001_synonymous
#          )

# write_tsv(df1,"data/006_Fig2/genomeScan/003_delVSsynOnChr_10000000window1000000step_addEffectiveCDSLength_addRecom.txt")

# write_tsv(df1,"/Users/Aoyue/project/wheatVMap2_1000/002_dataAnalysis/004_annoDB/009_genomeScan_delvcSyn/002/001_delVSsynOnChr_20000000window5000000step_addEffectiveCDSLength.txt")

### step2:根据step1的结果，求出 windowScan 滑窗后的数据框
# library(tidyverse)
# library(windowscanr)
# df1 <- read_tsv("data/006_Fig2/genomeScan/003_delVSsynOnChr_10000000window1000000step_addEffectiveCDSLength_addRecom.txt") 
# 
# df <- df1 %>% 
#   filter(!is.na(.$DelSynRatio),!is.na(.$RecombinationRate)) 
# 
# ### check del/syn 最大值和最小值
# max(df$DelSynRatio)
# min(df$DelSynRatio)
# max(df$RecombinationRate)
# min(df$RecombinationRate)
# 
# pos_win <- winScan(x = df, 
#                    groups = "CHROM", 
#                    position = "BIN_START_scale" , 
#                    values = c("NonsynSynRatio", 
#                               "DelSynRatio",
#                               "PPH2SynRatio",
#                               "GERPSynRatio",
#                               "StopGainSynRatio",
#                               "RecombinationRate"), 
#                    win_size = 5,
#                    win_step = 2,
#                    funs = c("mean", "sd"))
# 
# head(pos_win)
# df <- pos_win
# df3 <- mutate(df,Sub=substring(df$CHROM,2,2))
# df <- df3
# 
# write_tsv(df,"data/006_Fig2/genomeScan/004_delVSsynOnChr_windowScanbyRpackage.txt")
```
## plot
```{r}
df <- read_tsv("data/006_Fig2/genomeScan/004_delVSsynOnChr_windowScanbyRpackage.txt")

# rec <- data.frame(x1=c(0,9,28,50,80), x2=c(9,28,50,80,100), y1=c(-Inf,-Inf,-Inf,-Inf,-Inf), y2=c(Inf,Inf,Inf,Inf,Inf), region=c('R1','R2a','C','R2b','R3'))
colB <- c("#fffcd4","#1e3c84","#3ab4c4","#3ab4c4","#1e3c84")
p <- df %>% 
  ggplot()+
  # geom_rect(data=rec, mapping=aes(xmin=x1, xmax=x2, ymin=y1, ymax=y2,fill=region))+
  annotate(geom = "rect", xmin = 28.14,xmax = 53.08, ymin = -Inf, ymax = Inf, fill = "#fffcd4",alpha = 1)+
  geom_line(mapping=aes(x=df$win_start,y=df$DelSynRatio_mean,group=CHROM),col='#F4A460')+
  stat_smooth(mapping=aes(x=df$win_start,y=df$DelSynRatio_mean),col = "#d5311c",method = "loess",formula = y~x,span = 0.04,se=F,size=1) +
  # geom_point()+
  geom_line(mapping=aes(x=df$win_start,y=df$RecombinationRate_mean/2,group=CHROM),col='#56B4E9')+
  stat_smooth(df, mapping=aes(x=df$win_start,y=df$RecombinationRate_mean/2),col = "blue",method = "loess",formula = y~x,span = 0.04,se=F,size=1) +
    # geom_point(df_centromere,mapping = aes(x=x,y=y),color='black',size=0.5)+
  xlab("Scale position")+ylab("Del SNPs/syn SNPs") +
  theme(legend.position = 'none')+
  scale_fill_manual(values = colB)+
  scale_y_continuous(expand = c(0, 0),sec.axis = sec_axis(~.*2, name="Recombination rate"))+
  scale_x_continuous(expand = c(0, 0))+
  theme(
    # plot.margin=unit(rep(1,4),'lines'),
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7)
    ,legend.position = 'none'
    ,text = element_text(size = 18))
p
# ggsave(p,filename = "~/Documents/test.pdf",height = 2,width = 3.6)
ggsave(p,filename = "~/Documents/test.png",height = 7.62,width = 17.98,units = "cm")

```


## nonsyn ratio
```{r}
df <- read_tsv("data/006_Fig2/genomeScan/004_delVSsynOnChr_windowScanbyRpackage.txt")

# rec <- data.frame(x1=c(0,9,28,50,80), x2=c(9,28,50,80,100), y1=c(-Inf,-Inf,-Inf,-Inf,-Inf), y2=c(Inf,Inf,Inf,Inf,Inf), region=c('R1','R2a','C','R2b','R3'))
colB <- c("#fffcd4","#1e3c84","#3ab4c4","#3ab4c4","#1e3c84")
p <- df %>% 
  ggplot()+
  # geom_rect(data=rec, mapping=aes(xmin=x1, xmax=x2, ymin=y1, ymax=y2,fill=region))+
  annotate(geom = "rect", xmin = 28.14,xmax = 53.08, ymin = -Inf, ymax = Inf, fill = "#fffcd4",alpha = 1)+
  geom_line(mapping=aes(x=df$win_start,y=df$NonsynSynRatio_mean,group=CHROM),col='#F4A460')+
  stat_smooth(mapping=aes(x=df$win_start,y=df$NonsynSynRatio_mean),col = "#d5311c",method = "loess",formula = y~x,span = 0.04,se=F,size=1) +
  # geom_point()+
  geom_line(mapping=aes(x=df$win_start,y=df$RecombinationRate_mean*2,group=CHROM),col='#56B4E9')+
  stat_smooth(df, mapping=aes(x=df$win_start,y=df$RecombinationRate_mean*2),col = "blue",method = "loess",formula = y~x,span = 0.04,se=F,size=1) +
    # geom_point(df_centromere,mapping = aes(x=x,y=y),color='black',size=0.5)+
  xlab("Scale position")+ylab("Nonsynonymous SNPs/syn SNPs") +
  theme(legend.position = 'none')+
  scale_fill_manual(values = colB)+
  scale_y_continuous(expand = c(0, 0),sec.axis = sec_axis(~./2, name="Recombination rate"))+
  scale_x_continuous(expand = c(0, 0))+
  theme(
    # plot.margin=unit(rep(1,4),'lines'),
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7)
    ,legend.position = 'none'
    ,text = element_text(size = 9))
p
ggsave(p,filename = "~/Documents/test.pdf",height = 2,width = 3.6)
```

## VEP high impact
```{r}
df <- read_tsv("data/006_Fig2/genomeScan/004_delVSsynOnChr_windowScanbyRpackage.txt")

# rec <- data.frame(x1=c(0,9,28,50,80), x2=c(9,28,50,80,100), y1=c(-Inf,-Inf,-Inf,-Inf,-Inf), y2=c(Inf,Inf,Inf,Inf,Inf), region=c('R1','R2a','C','R2b','R3'))
colB <- c("#fffcd4","#1e3c84","#3ab4c4","#3ab4c4","#1e3c84")
p <- df %>% 
  ggplot()+
  # geom_rect(data=rec, mapping=aes(xmin=x1, xmax=x2, ymin=y1, ymax=y2,fill=region))+
  annotate(geom = "rect", xmin = 28.14,xmax = 53.08, ymin = -Inf, ymax = Inf, fill = "#fffcd4",alpha = 1)+
  geom_line(mapping=aes(x=df$win_start,y=df$StopGainSynRatio_mean,group=CHROM),col='#F4A460')+
  stat_smooth(mapping=aes(x=df$win_start,y=df$StopGainSynRatio_mean),col = "#d5311c",method = "loess",formula = y~x,span = 0.04,se=F,size=1) +
  # geom_point()+
  geom_line(mapping=aes(x=df$win_start,y=df$RecombinationRate_mean/9,group=CHROM),col='#56B4E9')+
  stat_smooth(df, mapping=aes(x=df$win_start,y=df$RecombinationRate_mean/9),col = "blue",method = "loess",formula = y~x,span = 0.04,se=F,size=1) +
    # geom_point(df_centromere,mapping = aes(x=x,y=y),color='black',size=0.5)+
  xlab("Scale position")+ylab("High impact SNPs/syn SNPs") +
  theme(legend.position = 'none')+
  scale_fill_manual(values = colB)+
  scale_y_continuous(expand = c(0, 0),sec.axis = sec_axis(~.*9, name="Recombination rate"))+
  scale_x_continuous(expand = c(0, 0))+
  theme(
    # plot.margin=unit(rep(1,4),'lines'),
    panel.background = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=0.7)
    ,legend.position = 'none'
    ,text = element_text(size = 9))
p
ggsave(p,filename = "~/Documents/test.pdf",height = 2,width = 3.6)
```




#### ============
## Extended data
#### ============
## SNP count in gene region

## 查看 coding 区有哪些类型及个数
```{r}
df_test <- dftotal %>% 
  filter(Region=="CDS") %>% 
  group_by(Variant_type) %>% 
  count()
```

## FigS:SNP count in gene region
UTR           Synonymous    Stop-loss     Stop-gain     Start-lost    Nonsynonymous
以上几种类型的计数
```{r}
### Goal: 将 Region列和Variant_type列柔和在一起，成为新的一列，即将UTR的信息加入分析
### 第一步，Region 列，是 UTR_3 UTR_5的，新列就命名为UTR，不是的，则按照Variant_type类型添加

dfUTR <- dftotal %>% 
  filter(Region == "UTR_3" | Region == "UTR_5") %>% 
  mutate(VariantGroup = "UTR")

dfnotUTR <- dftotal %>% 
  filter(Region == "CDS") %>% 
  mutate(VariantGroup = Variant_type)

dfnew <- rbind(dfUTR,dfnotUTR)

df_variantGroup<- dfnew %>% 
  filter(!is.na(VariantGroup)) %>%
  group_by(VariantGroup) %>% 
  summarise(n=n()) %>% 
  ungroup() %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(desc(VariantGroup))

# colB <- c("#a9a9a9","#e69f00","#871f12","#d5311c","#5a150c","#004680")
colB <- c("#e69f00","#871f12","#d5311c","#5a150c","#004680","#a9a9a9")

df_variantGroup$VariantGroup <- factor(df_variantGroup$VariantGroup,levels = c("NONSYNONYMOUS","START-LOST","STOP-GAIN", "STOP-LOSS","SYNONYMOUS","UTR"), labels = c("Nonsynonymous","Start-lost", "Stop-gain","Stop-loss","Synonymous","UTR"))
p <- df_variantGroup %>%
ggplot(aes(x=VariantGroup,y=n,fill=VariantGroup))+
  geom_bar(stat = 'identity',position = "stack")+
  # geom_text(aes(label=n),size=3,nudge_y = 1)+
    geom_text(aes(label=paste(n,"\n(",round(fre,3),")",sep = "")),size=3,
            position = position_dodge(0.9))+
  labs(x="", y="SNP count")+
  # scale_fill_brewer(palette = "Set2")+
  scale_fill_manual(values = colB)+
  # scale_fill_brewer(palette = "Blues")+
  theme_classic()+
  theme(axis.text.x= element_text(angle=45, vjust=1,hjust = 1))+
  theme(legend.position = "none",plot.margin=unit(rep(1,4),'lines'),
          axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 12))
p

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

```

## FigS: AABBDD DAF
#### AABBDD: keep all sites in segregation
AABBDD AABB DD 三个群体都必须有分离
```{r}
###step1: add AAF to GeneDB
dfAAF <- read_tsv("data/015_AAF/002_aaf_byPloidy/002_AAF_byPloidy.txt.gz")

df <- dftotal %>% 
  # mutate(Sub=str_sub(Transcript, 9, 9)) %>%
  ### 定义load
  mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  filter(!is.na(LoadGroup)) %>% 
  filter(!Ancestral == "NA") %>% 
  filter(Ancestral==Ref|Ancestral==Alt) %>% 
  ### 将AAF数据代入其中
  left_join(.,dfAAF,by=c("Chr","Pos")) %>% 
  mutate(DAF_42 = ifelse(Ref==Ancestral,AAF_AABB_DD,1-AAF_AABB_DD)) %>% 
  mutate(DAF_06 = ifelse(Ref==Ancestral,AAF_AABBDD,1-AAF_AABBDD))

#### ======================================================================
### 必须在六倍体群体分离
df06 <- df %>% 
  mutate(Sub=factor(Sub,levels = c("A","B","D"),labels = c("A subgenome","B subgenome","D subgenome"))) %>% 
  select(Sub,LoadGroup,DAF_42,DAF_06) %>% 
  ### data deal
  filter(!is.na(DAF_06)) %>%  ###*****
  mutate(bin=cut_width(DAF_06,width = 0.05,boundary = 0)) %>% 
  group_by(Sub,LoadGroup) %>% 
  count(bin) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(Sub,desc(LoadGroup)) 
 
#### facet by sub
colB <- rep(c("#d5311c","#e69f00","#004680"),2)
q <- df06 %>% 
  ggplot(aes(x=bin,y=fre,group=LoadGroup,fill=LoadGroup))+
  geom_bar(stat = 'identity',position = "dodge")+
  labs(y="Proportion",x="Derived allele frequency")+
  # scale_color_brewer(palette = "Set2")+
  scale_fill_manual(values = colB)+
  # scale_fill_brewer(palette = "Blues")+
  scale_y_continuous(limits = c(0,0.855),breaks = seq(0,0.8,0.2))+
  scale_x_discrete(labels=c("0","","","","0.2","","","","0.4","","","","0.6","","","","0.8","","","",""))+
  facet_grid(Sub~.)+
  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 = 11))
q

ggsave(q,filename = "~/Documents/DAF_06.pdf",height = 3.17,width = 3.17)


#### ======================================================================
### 必须在四二倍体群体分离
df42 <- df %>% 
  mutate(Sub=factor(Sub,levels = c("A","B","D"),labels = c("A subgenome","B subgenome","D subgenome"))) %>% 
  select(Sub,LoadGroup,DAF_42,DAF_06) %>% 
  ### data deal
  filter(!is.na(DAF_42)) %>%  ###*****
  mutate(bin=cut_width(DAF_42,width = 0.05,boundary = 0)) %>% 
  group_by(Sub,LoadGroup) %>% 
  count(bin) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(Sub,desc(LoadGroup)) 
 
#### facet by sub
colB <- rep(c("#d5311c","#e69f00","#004680"),2)
q <- df42 %>% 
  ggplot(aes(x=bin,y=fre,group=LoadGroup,fill=LoadGroup))+
  geom_bar(stat = 'identity',position = "dodge")+
  labs(y="Proportion",x="Derived allele frequency")+
  # scale_color_brewer(palette = "Set2")+
  scale_fill_manual(values = colB)+
  # scale_fill_brewer(palette = "Blues")+
  scale_y_continuous(limits = c(0,0.855),breaks = seq(0,0.8,0.2))+
  scale_x_discrete(labels=c("0","","","","0.2","","","","0.4","","","","0.6","","","","0.8","","","",""))+
  facet_grid(Sub~.)+
  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 = 11))
q

ggsave(q,filename = "~/Documents/DAF_42.pdf",height = 3.17,width = 3.17)

```

### 显著性检验
colnames(df)
 [1] "ID"                           "Chr"                         
 [3] "Pos"                          "Ref"                         
 [5] "Alt"                          "Major"                       
 [7] "Minor"                        "Maf"                         
 [9] "Transcript"                   "Ancestral"                   
[11] "DAF"                          "Gerp"                        
[13] "Region"                       "Variant_type"                
[15] "Alt_SIFT"                     "Ref_SIFT"                    
[17] "Derived_SIFT"                 "Effect_VEP"                  
[19] "Impact_VEP"                   "Effect_snpEff"               
[21] "Impact_snpEff"                "Rank_Total"                  
[23] "Gerp_16way"                   "AMINO_POS"                   
[25] "UniProtKB_AC"                 "LIST_S2"                     
[27] "PhyloP"                       "Amino_Size"                  
[29] "PhyloP_RefMask"               "AlleleAge_J"                 
[31] "AlleleAge_M"                  "AlleleAge_R"                 
[33] "Alt_PolyPhen2_prediction"     "Alt_PolyPhen2_class"         
[35] "Alt_PolyPhen2_prob"           "Ref_PolyPhen2_prediction"    
[37] "Ref_PolyPhen2_class"          "Ref_PolyPhen2_prob"          
[39] "Derived_PolyPhen2_prediction" "Derived_PolyPhen2_class"     
[41] "Derived_PolyPhen2_prob"       "Sub"                         
[43] "LoadGroup"                    "AAF_AABB_DD"                 
[45] "AAF_AABBDD"                   "DAF_42"                      
[47] "DAF_06"                      
```{r}
df_sig06 <- df %>% 
  select(LoadGroup, DAF_06, Sub) %>% 
  filter(!is.na(DAF_06), DAF_06 < 0.05)

df_sig06$Sub <- as.factor(df_sig06$Sub)
df_sig06$LoadGroup <- as.factor(df_sig06$LoadGroup)

df_sig06_Asub <- df_sig06 %>% filter(Sub=="A")
df_sig06_Bsub <- df_sig06 %>% filter(Sub=="B")
df_sig06_Dsub <- df_sig06 %>% filter(Sub=="D")

pairwise.wilcox.test(df_sig06_Asub$DAF_06,df_sig06_Asub$LoadGroup,p.adjust.method = "BH")
print("***************A sub completed**************************")

pairwise.wilcox.test(df_sig06_Bsub$DAF_06,df_sig06_Bsub$LoadGroup,p.adjust.method = "BH")
print("***************B sub completed**************************")

pairwise.wilcox.test(df_sig06_Dsub$DAF_06,df_sig06_Dsub$LoadGroup,p.adjust.method = "BH")
print("***************D sub completed**************************")

pairwise.wilcox.test(df_sig06$DAF_06,df_sig06$LoadGroup,p.adjust.method = "BH")
print("*************** no Sub completed*********************")


```


## DAF DAF DAF DAF DAF
### !only daf 坐标轴改变,AB 合并在一起

```{r}
df <- dftotal %>% 
  mutate(Sub=str_sub(Transcript, 9, 9)) %>%
  filter(Region=="CDS",Variant_type=="SYNONYMOUS"|Variant_type=="NONSYNONYMOUS") %>% 
  filter(!is.na(Ancestral),Ancestral == Ref | Ancestral == Alt) %>%
  mutate(Gerp_16way=ifelse(is.na(Gerp_16way),-100,Gerp_16way)) %>%
  mutate(Derived_PolyPhen2_prob=ifelse(is.na(Derived_PolyPhen2_prob),-100,Derived_PolyPhen2_prob)) %>%
  mutate(LoadGroup=ifelse(Derived_PolyPhen2_prob>=0.5 & Variant_type %in% c("NONSYNONYMOUS") & Gerp_16way >=1.5, "Deleterious", ifelse(Variant_type=="SYNONYMOUS", "Synonymous", ifelse(Variant_type=="NONSYNONYMOUS", "Nonsynonymous", NA)))) %>% 
  mutate(IfRefisAnc = if_else(Ref == Ancestral,"Anc","Der")) %>% 
  ### data deal
  filter(!is.na(DAF)) %>% 
  mutate(bin=cut_width(DAF,width = 0.05,boundary = 0)) %>% 
  mutate(mergedSub = if_else(Sub == "D","D","AB")) %>% 
  group_by(mergedSub,LoadGroup) %>% 
  count(bin) %>% 
  mutate(fre=n/sum(n)) %>% 
  arrange(mergedSub,desc(LoadGroup)) %>% 
  ungroup() %>% 
  mutate(binStart = round(as.numeric(str_sub(str_split_fixed(bin,",",2)[,1],2)),2))
 
#### facet by sub
colB <- rep(c("#d5311c","#e69f00","#004680"),2)
q <- df %>% 
  filter(mergedSub=="AB") %>% 
  ggplot(aes(x=as.character(binStart) ,y=fre,group=LoadGroup,fill=LoadGroup))+
  geom_bar(stat = 'identity',position = "dodge")+
  labs(y="Proportion",x="Derived allele frequency")+
  # scale_color_brewer(palette = "Set2")+
  scale_fill_manual(values = colB)+
  # scale_y_continuous(limits = c(0,0.88),breaks = seq(0,0.8,0.2))+
  scale_x_discrete(labels=c("0","","","","0.2","","","","0.4","","","","0.6","","","","0.8","","","",""))+
  # scale_fill_brewer(palette = "Blues")+
  # facet_grid(mergedSub~.)+
  theme_classic()+
  # theme(axis.text.x= element_text(angle=90, vjust=1,hjust = 1))+
  theme(plot.margin=unit(rep(1,4),'lines'),
        legend.title = element_blank(),legend.position = "none",
        axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 10))+
    theme(strip.text = element_blank())


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

### D sub
q <- df %>% 
  filter(mergedSub=="D") %>% 
  ggplot(aes(x=as.character(binStart) ,y=fre,group=LoadGroup,fill=LoadGroup))+
  geom_bar(stat = 'identity',position = "dodge")+
  labs(y="Proportion",x="Derived allele frequency")+
  # scale_color_brewer(palette = "Set2")+
  scale_fill_manual(values = colB)+
  # scale_y_continuous(limits = c(0,0.88),breaks = seq(0,0.8,0.2))+
  scale_x_discrete(labels=c("0","","","","0.2","","","","0.4","","","","0.6","","","","0.8","","","",""))+
  # scale_fill_brewer(palette = "Blues")+
  # facet_grid(mergedSub~.)+
  theme_classic()+
  # theme(axis.text.x= element_text(angle=90, vjust=1,hjust = 1))+
  theme(plot.margin=unit(rep(1,4),'lines'),
        legend.title = element_blank(),legend.position = "none",
        axis.line = element_line(size=0.4, colour = "black"),text = element_text(size = 10))+
    theme(strip.text = element_blank())


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

```