Fig01.Rmd

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

# Fig1

```{r}
setwd("/Users/Aoyue/project/wheatVMapII/003_dataAnalysis/zz_rProject/R4VMap2Figs")
```

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

## Fig1|World map
### 检查Taxa来自多少个国家
```{r}
df <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt")
dfcountry <- df %>% filter(!is.na(orgCty)) %>%  group_by(orgCty) %>% count()
numCountry <- nrow(dfcountry)
print(str_c("The sample in VMap2 project comes from ",numCountry," countries"))

### 按照亚群分类
df <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt")
dfcountry <- df %>% 
  filter(!is.na(orgCty)) %>% 
  group_by(Subspecies_by6_TreeValidated,orgCty) %>% count() %>% 
  group_by(Subspecies_by6_TreeValidated) %>% count()

numCountry <- nrow(dfcountry)
print(str_c("The sample in VMap2 project comes from ",numCountry," countries"))

```
### map group by subspecies- 地图白色背景
- 726 个点是有经纬度信息的
```{r}
library(maps)
library(maptools)
df <- read_tsv("data/001_TaxaDB/WheatVMap2_GermplasmInfo.txt")

unique(df$Subspecies_by6_TreeValidated)
factorA <- c("Strangulata","Wild_emmer","Domesticated_emmer","Free_threshing_tetraploids","OtherTetraploids","Landrace","Cultivar","OtherHexaploids")
colB <- c("#87cef9","#ffd702","#7f5701","#016699","#00f3ff", "#fc6e6e","#9900ff","#fe63c2")
df$Subspecies_by6_TreeValidated <- factor(df$Subspecies_by6_TreeValidated,levels = factorA)

df <- df %>% 
  arrange(desc(Subspecies_by6_TreeValidated)) %>% 
  filter(!is.na(Longitude))

mp<-NULL #定义一个空的地图
# 注意： borders函数中的colour是地图国家边界线的颜色,fill是国家的填充颜色
# mapworld<-borders("world",colour = "#dedfe0",fill="#dedfe0") #绘制基本地图
# fcfcfc
mapworld<-borders("world",colour = "#f6f6f6",fill="#f6f6f6") #绘制基本地图， 白色背景
mpvoid<-ggplot()+mapworld + ylim(-55,90)+theme_void()
# mpvoid
mp<-mpvoid+
  geom_point(data=df,
             aes(x=df$Longitude,y=df$Latitude,fill = df$Subspecies_by6_TreeValidated),
             shape=21,color="gray",stroke = 0.08,
             alpha=0.75,size=2.5)+
  scale_fill_manual(values = colB)+
  theme(legend.position = 'none')
mp
# ggsave(mp,filename = "~/Documents/test.pdf", height = 4,width = 10)
ggsave(mp,filename = "~/Documents/test.pdf", height = 2.76 ,width = 6.89)
# ggsave(mp,filename = "~/Documents/test.png", height = 2.76 ,width = 6.89)

```

#### ====
## Fig.1\|PopNum barplot

```{r}
Subspecies <- c("Ae.tauschii","Wild emmer","Domesticated emmer","Free-threshing tetraploids",
                "Other tetraploid","Landrace","Cultivar","Other hexaploid")
n <- c(36,91,52,50,19,306,223,285)
id <- rep(1,8)
df <- data.frame(id=id,Subspecies=Subspecies,n=n)
head(df)
df
colB <- c('#87cef9','#fe63c2','#9900ff','#fc6e6e','#00f3ff','#016699','#7f5701','#ffd702')
df$Subspecies <- factor(df$Subspecies,
                        levels = c( "Ae.tauschii","Other hexaploid","Cultivar","Landrace", "Other tetraploid",
                                   "Free-threshing tetraploids", "Domesticated emmer","Wild emmer"))
p <- ggplot(df,aes(x=id,y=n,fill=Subspecies))+
  geom_bar(stat="identity",alpha=0.95)+
  scale_fill_manual(values = colB)+
  coord_flip()+
  theme_void()+
  theme(legend.position = 'none')
p
ggsave(p,filename = "~/Documents/test.pdf",height = 0.5,width = 4)

```

### pie plot
```{r}
# Create Data
data <- data.frame(
  group=c('AABBDD','AABB','DD'),
  value=c(814,212,36)
)

# Compute the position of labels
data <- data %>% 
  arrange(desc(group)) %>%
  mutate(prop = value / sum(data$value) *100) %>%
  mutate(ypos = cumsum(prop)- 0.5*prop )

# Basic piechart
p <- ggplot(data, aes(x="", y=prop, fill=group)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void() + 
  theme(legend.position="none") +
  geom_text(aes(y = ypos, label = group), color = "black", size=3) +
  scale_fill_manual(values = c("#ffd702",'#fc6e6e','#87cef9'))
p  
ggsave(p,filename = "/Users/Aoyue/Documents/test.pdf",height = 2.57,width = 2.57)
```


#### =================
## LD
### Data import
```{r}
table_ploidy <- read.delim("data/012_LD/LD.ploid.txt.gz")
table_subspecies <- read.delim("data/012_LD/LD.subspecies.txt.gz")
table_subgenome <- read.delim("data/012_LD/LD.subgenome.txt.gz")
```

### only hexaploid by subgenome
```{r}
colB <- c("#fd8582","#967bce","#4bcdc6") 

p <- table_subgenome %>% 
  mutate(Ploidy = factor(Ploidy,levels = c("tetraploid","hexaploid","diploid"),labels = c("AABB","AABBDD","DD"))) %>% 
  mutate(Sub=factor(Sub,levels = c("AA","BB","DD"),labels = c("A","B","D"))) %>% 
  filter(Ploidy == "AABBDD") %>%
  slice_sample(prop = .1) %>%
  ggplot(aes(x=window_kb/1000,y=meanOfR2, color=Sub, fill=Sub))+
  # geom_point(alpha=0.05)+
  # geom_smooth()+
  stat_smooth(method = "gam", formula = y ~ s(x,k=110),span = 0.2, size = 1,se=T) +
  scale_color_manual(values = colB)+
  labs(x="Pairwise distance (Mb)", y=expression("Mean"~italic(r)^2))+
  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 = 8))

p

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

# ggsave(p,filename = "~/Documents/test.pdf",height = 1.6 ,width = 1.76)
```
### claculate the r2
```{r}
dfnew <- table_subgenome %>% filter(Ploidy=="hexaploid")

```

```{r}
df1 <- dfnew %>% 
  filter(Sub == "AA") %>% 
  filter(meanOfR2 > 0.49, meanOfR2 < 0.5)

### 7662.08 kb

df2 <- dfnew %>% 
  filter(Sub == "BB") %>% 
  filter(meanOfR2 > 0.49, meanOfR2 < 0.5)

### 460.2692 kb

df3 <- dfnew %>% 
  filter(Sub == "DD") %>% 
  filter(meanOfR2 > 0.49, meanOfR2 < 0.5)

### 4404.953 kb

mean(df1$window_kb)
mean(df2$window_kb)
mean(df3$window_kb)


```

### by subspecies
```{r}
# unique(table_subspecies$Group)
factorA <- c("wildEmmer","domesticatedEmmer","freeThreshTetraploid","landrace","cultivar","tauschii")
labelsA <- c("WE","DE","FTT","LR","CL","AT")
colB <- c("#ffd702","#7f5701","#016699", "#fc6e6e","#9900ff","#87cef9")

p <- table_subspecies %>% 
  mutate(group=factor(Group,levels=factorA,labels=labelsA)) %>% 
  group_by(group) %>% 
  slice_sample(prop = .1) %>%
  ggplot(aes(x=window_kb/1000,y=meanOfR2, color=group, fill=group))+
  # geom_point(alpha=0.05)+
  # geom_smooth(se = T)+
  stat_smooth(method = "gam", formula = y ~ s(x,k=110),span = 0.2, size = 1,se=T) +
  # geom_line(size=0.0001,alpha=0.9)+
  scale_color_manual(values = colB)+
  labs(x="Pairwise distance (Mb)", y=expression("Mean"~italic(r)^2))+
  scale_x_continuous(limits = c(0,20))+
  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 = 8))

p

# ggsave(p,filename = "~/Documents/test.pdf",height = 1.6 ,width = 1.76)
ggsave(p,filename = "~/Documents/test.pdf",height = 2.4 ,width = 3.23)

```

### claculate the r2

```{r}
df1 <- table_subspecies %>% 
  filter(Group=="landrace") %>% 
  filter(meanOfR2 > 0.49, meanOfR2 < 0.5)


df2 <- table_subspecies %>% 
  filter(Group=="cultivar") %>% 
  filter(meanOfR2 > 0.49, meanOfR2 < 0.5)


df3 <- table_subspecies %>% 
  filter(Group=="wildEmmer") %>% 
  filter(meanOfR2 > 0.19, meanOfR2 < 0.2)


df4 <- table_subspecies %>% 
  filter(Group=="domesticatedEmmer") %>% 
  filter(meanOfR2 > 0.49, meanOfR2 < 0.5)

df5 <- table_subspecies %>% 
  filter(Group=="freeThreshTetraploid") %>% 
  filter(meanOfR2 > 0.49, meanOfR2 < 0.5)

df6 <- table_subspecies %>% 
  filter(Group=="tauschii") %>% 
  filter(meanOfR2 > 0.19, meanOfR2 < 0.2)


mean(df1$window_kb)
mean(df2$window_kb)
mean(df3$window_kb)
mean(df4$window_kb)
mean(df5$window_kb)
mean(df6$window_kb)


```

```{r}
# unique(table_subspecies$Group)
factorA <- c("wildEmmer","domesticatedEmmer","freeThreshTetraploid","landrace","cultivar","tauschii")
labelsA <- c("WE","DE","FTT","LR","CL","AT")
colB <- c("#ffd702","#7f5701","#016699", "#fc6e6e","#9900ff","#87cef9")

p <- table_subspecies %>% 
  mutate(group=factor(Group,levels=factorA,labels=labelsA)) %>% 
  group_by(group) %>% 
  slice_sample(prop = .1) %>%
  filter(group=="LR") %>% 
  ggplot(aes(x=window_kb/1000,y=meanOfR2, color=group, fill=group))+
  # geom_point(alpha=0.05)+
  # geom_smooth(se = T)+
  stat_smooth(method = "gam", formula = y ~ s(x,k=110),span = 0.2, size = 1,se=T) +
  # geom_line(size=0.0001,alpha=0.9)+
  scale_color_manual(values = colB)+
  labs(x="Pairwise distance (Mb)", y=expression("Mean"~italic(r)^2))+
  scale_x_continuous(limits = c(0,20))+
  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 = 8))

p

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