article_figures.Rmd

---
title: "article_figures.Rmd"
author: "Marie Bourdon"
date: "June 2022"
output: html_document
---

# Goal and raw data

The goal of this script is to produce figure for the stuart package manuscript.

This scripts uses data from the package, and other files found in the "files" directory.

This script uses the qtl_plot() ggplot based function to plot QTL mapping results. This function is in the script "QTL_plot.R".

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(magrittr)
library(qtl)
library(cowplot)
library(grid)
library(gridExtra)
library(gridGraphics)
library(scales)

library(stuart)

source("files/QTL_plot.R")
source("files/find_linked_markers.R")
```

# Load R objects from datasets

```{r}
load("data1/data1.rda")
load("data2/data2.rda")
load("data3/data3.rda")
annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv"))
```


# Graphs

## Grid mark_prop

### Graph missing genotypes

```{r graph_NA}
na_plot <- tab2 %>% mutate(prop_NA=n_NA/176) %>% ggplot(aes(x=prop_NA)) +
  geom_histogram() +
  scale_y_log10(expand = c(.01, .01)) +
  scale_x_continuous(expand = c(.005, .005)) +
  theme_classic() +
  labs(title="Proportion of missing genotypes",
       x="Proportion of NA",y="Number of markers") +
  theme(
    aspect.ratio=0.8,
    plot.title = element_text(hjust = 0.4,size=14))

na_plot
```

Proportions of markers with more than 75% of missing genotypes:
```{r prop_missing}
tab2 %>% mutate(prop_NA=n_NA/176) %>% filter(prop_NA > 0.50) %>% summarise(p=n()/count(tab2)%>%pull())
```


### Graph proportion of genotypes

```{r graph_prop}
prop_plot <- tab2 %>% filter(n_NA<88) %>% filter(!chr %in% c("M","X","Y")) %>%
  ggplot(aes(x=n_HM1/(n_HM1+n_HM2+n_HT),y=n_HM2/(n_HM1+n_HM2+n_HT),color=as.factor(exclude_prop))) +
  geom_point() +
  scale_color_manual(values=c("#66bd63","#b2182b"),labels = c("Retained", "Excluded")) +
  geom_hline(yintercept = 0.1,linetype="dashed",size=.3) +
  geom_vline(xintercept = 0.1,linetype="dashed",size=.3) +
  geom_abline(intercept = 0.9, slope=-1,linetype="dashed",size=.3) +
  scale_x_continuous(expand = c(.005, .005)) +
  scale_y_continuous(expand = c(.005, .005)) +
  labs(title="Exclusion of markers by mark_prop()",
       x="Proportion of A/A individuals",
       y="Proportion of B/B individuals",
       color="Exclusion") +
  theme_classic() +
  theme(aspect.ratio=0.8,
        legend.position=c(0.8,0.8),
        legend.title = element_blank()) +
  theme(plot.title = element_text(hjust = 0.4,size=14)) +
  geom_segment(x = 0.2, y = 0.78,
               xend = 0.25, yend = 0.72,
               color = "#1d91c0",
               arrow = arrow(type="closed",length=unit(4,"points"))) +
  geom_segment(x = 0.39, y = 0.58,
               xend = 0.44, yend = 0.52,
               color = "#1d91c0",
               arrow = arrow(type="closed",length=unit(4,"points")))
prop_plot
```

### Grid

Figure with graphs proportion of NA and exclusion depending on genotype proportions


```{r grid_na_prop}
gridnaprop <- plot_grid(na_plot,prop_plot,
                  ncol=1,
                  labels = c('A', 'B'), 
                  label_size = 20,
                  byrow = TRUE
                  )

gridnaprop

ggsave("figures/fig2.pdf",gridnaprop,width=5,height=8)
ggsave("figures/fig2.png",gridnaprop,width=5,height=8)


rm(na_plot,prop_plot)
```

## Graph before after


### Grid before



```{r grid_before}
grid_before <- plot_grid(as_grob(plotmap_before),rec_ratios_before,pheno_before_annot,pheno_before_zoom,ncol=2,labels=c("A","B","C","D"),label_size=20)
grid_before
ggsave("figures/fig1prev.pdf",grid_before,width=14,height = 10)
```


### Grid after



```{r grid fig3, fig.height = 7, fig.width = 13, fig.align = "center"}
grid_after <- plot_grid(as_grob(plotmap_after),rec_ratios_after,pheno_after_plot,pheno_after_zoom,ncol=2,labels=c("A","B","C","D"),label_size=20)
grid_after
ggsave("figures/fig3prev.pdf",grid_after,width=14,height = 10)
```

# Differences between real and reference genotypes for parental strains

```{r dif_number}
dif <- full_join(strains,strns_ref,by=c("marker","chr","cM_cox")) %>% 
  mutate(dif=case_when((!parent1.x%in%c("N","H") & 
                          !parent1.y%in%c("N","H") & 
                          parent1.x!=parent1.y) ~ 1,
                       (!parent2.x%in%c("N","H") &
                          !parent2.y%in%c("N","H") & 
                          parent2.x!=parent2.y) ~ 1, 
                       T~0))

dif %>% filter(dif==1) %>% count()

rm(dif)
```

## Table alleles different between parental strains and F2s

```{r allele}
#investigation of the role of mark_allele function
#prove that some marker with non corresponding alleles between parents and F2

#keep only markers that are exlcuded with mark_allele
allele <- tab2 %>% 
  filter(exclude_allele==1&exclude_poly==0&exclude_prop==0&exclude_na==0&exclude_estmap==0)
strains_allele <- strains %>% filter(marker %in% allele$marker)

#number of markers not excluded by other functoins
tab2 %>% 
  filter(exclude_poly==0&exclude_prop==0&exclude_na==0&exclude_estmap==0) %>% nrow()

#join with strains genotypes to have parental strains
allele <- left_join(allele,strains_allele,by=c("marker"="marker"))


#most of markers excluded with mark_allele that were not excluded with other functions have N/H as genotype for parents
#keep only those with non missing/heterozygous genotypes
allele %<>% filter(parent1 != "N" & parent2 != "N")
allele %<>% select(marker,parent1,parent2,allele_1,allele_2)


print(allele)

print(xtable::xtable(allele, type = "latex"), file = "tables/tab_alleles.tex",include.rownames=FALSE)
rm(allele,strains_allele)
```

## Differences between real and reference genotypes for parental strains

```{r dif_table}
dif <- full_join(strains,strns_ref,by=c("marker","chr","cM_cox")) %>% 
  mutate(dif=case_when((!parent1.x%in%c("N","H") & 
                          !parent1.y%in%c("N","H") & 
                          parent1.x!=parent1.y) ~ 1,
                       (!parent2.x%in%c("N","H") &
                          !parent2.y%in%c("N","H") & 
                          parent2.x!=parent2.y) ~ 1, 
                       T~0))

dif %>% filter(dif==1) %>% count()
```

## Pheno data format

```{r pheno}
format_pheno <- phenos[1:6,]
write_csv(format_pheno,file="figures/tableS1.csv")
print(xtable::xtable(format_pheno, type = "latex"), file = "tables/tab_alleles.tex",include.rownames=FALSE)
```


## Markers with unmatching alleles that lead to distortion of genetic map

```{r}
tab2 %>% filter(exclude_allele==1 & exclude_prop==0) %>% select(marker:n_NA) %>% left_join(.,strains) %>% filter(!parent1 %in% c("N","H") & !parent2 %in% c("N","H"))

allele_rec <- tab2 %>% filter(marker %in% c("gJAX00038569","mUNC21540855","gUNC21555204","gUNC21596600")) %>% 
  select(marker:n_NA) %>% left_join(.,strains,by = c("marker", "chr", "cM_cox"))

# gJAX00038569
map_1 <- newmap_before[["12"]][c("gUNC21523346","gJAX00038569","gUNCHS034222")] #50cM on each side
pos_1 <- cross_before[["geno"]][["12"]][["map"]][c("gUNC21523346","gJAX00038569","gUNCHS034222")]

# mUNC21540855
map_2 <- newmap_before[["12"]][c("gUNCHS034222","mUNC21540855","SFJ123443466")] #50cM on each side
pos_2 <- cross_before[["geno"]][["12"]][["map"]][c("gUNCHS034222","mUNC21540855","SFJ123443466")]

# gUNC21555204
map_3 <- newmap_before[["12"]][c("SFJ123443466","gUNC21555204","gJAX00340618")] #30-60cM on each side
pos_3 <- cross_before[["geno"]][["12"]][["map"]][c("SFJ123443466","gUNC21555204","gJAX00340618")]

# gUNC21596600
map_4 <- newmap_before[["12"]][c("JAX00341204","gUNC21596600","UNC21601697")] #30-60cM on each side
pos_4 <- cross_before[["geno"]][["12"]][["map"]][c("JAX00341204","gUNC21596600","UNC21601697")]


df <- tibble(gJAX00038569 = c(map_1,pos_1),
       mUNC21540855 = c(map_2,pos_2),
       gUNC21555204 = c(map_3,pos_3),
       gUNC21596600 = c(map_4,pos_4))
df <- as_tibble(cbind(nms = names(df), t(df))) %>% rename(marker = nms,
                                                          rec_before = V2,
                                                          rec = V3,
                                                          rec_after = V4,
                                                          pos_before = V5,
                                                          pos = V6,
                                                          pos_after = V7)

df %<>% mutate(rec_before = as.numeric(rec) - as.numeric(rec_before),
              rec_after = as.numeric(rec_after) - as.numeric(rec),
              pos_before = as.numeric(pos) - as.numeric(pos_before),
              pos_after = as.numeric(pos_after) - as.numeric(pos)) %>%
  select(marker,rec_before,rec_after,pos_before,pos_after)

allele_rec <- full_join(allele_rec,df,by="marker") %>% select(marker:parent2,rec_before,rec_after,pos_before,pos_after)
allele_rec


rm(map_1,map_2,map_3,map_4,pos_1,pos_2,pos_3,pos_4,df,allele_rec)
```

## All 3 datasets

### Phenotype distribution

```{r}
grid_pheno <- plot_grid(pheno_data1,pheno_data2,pheno_data3,ncol=2,labels=c("A","B","C","D"),label_size=20)
grid_pheno

ggsave("figures/figS1.pdf",grid_pheno,width=14,height = 10)
```

### Grid with all 3 data sets

```{r}
colnames(peak1_tab) <- c("marker","chr","pos","allele\n1","allele\n2",
                 "n\nHM1","n\nHM2","n\nHT","n\nNA")
colnames(peak3_tab) <- c("marker","chr","pos","allele\n1","allele\n2",
                 "n\nHM1","n\nHM2","n\nHT","n\nNA")


narrow_grid <- ggdraw() +
  draw_plot(pheno_before_annot_data2,0,.7,.5,.3) +
  draw_plot(pheno_before_annot_data3,.5,.7,.5,.3) +
  draw_plot(peak1,0,.35,.3,.3) +
  draw_plot(tableGrob(peak1_tab,rows=NULL),.3,.35,.7,.3) +
  draw_plot(peak3,0,0,.3,.3) +
  draw_plot(tableGrob(peak3_tab,rows=NULL),.3,0,.7,.3) +
  draw_plot_label(c("A","B","C","D","E","F"),
                  c(0,.5,0,.3,0,.3),c(1,1,.65,.65,.3,.3))

narrow_grid
ggsave(narrow_grid,file="figures/figS3prev.pdf",width=10,height=14)
```


## Dataset sup

### Map before and rec ratio distrib

```{r}
grid_mapbef <- plot_grid(plotmap_before_data2,plotmap_before_data3,rec_ratios_before_data2,rec_ratios_before_data3,ncol=2,labels=c("A","B","C","D"),label_size=20)

grid_mapbef
ggsave("figures/figS2prev.pdf",grid_mapbef ,width=14,height = 10)
```

### Map after and rec ratio distrib

```{r}
grid_mapaf <- plot_grid(plotmap_after_data2,plotmap_after_data3,rec_ratios_after_data2,rec_ratios_after_data3,ncol=2,labels=c("A","B","C","D"),label_size=20)

ggsave("figures/figS4prev.pdf",grid_mapaf,width=14,height = 10)
```

### Scans

```{r}
grid_scans <- plot_grid(pheno_after_plot_data2,pheno_after_plot_data3,ncol=2,labels=c("A","B"),label_size=20)

ggsave("figures/figS5.pdf",grid_scans,width=14,height = 5)
```


## Comparision of non mendelian/wrong alleles markers between the 4 datasets


```{r}
colnames(tab2_data2) <- c(colnames(tab2_data2)[1:3],paste0(colnames(tab2_data2)[4:ncol(tab2_data2)],".2"))
colnames(tab2_data3) <- c(colnames(tab2_data3)[1:3],paste0(colnames(tab2_data3)[4:ncol(tab2_data3)],".3"))

tab_alldata <- full_join(tab2,tab2_data2) %>% full_join(.,tab2_data3)

#keep all markers with non mendelian prop in at least 1 dataset
prop_rem <- tab_alldata %>% filter((exclude_poly==0&exclude_prop==1) | (exclude_poly==0&exclude_allele==1) |
                       (exclude_poly.2==0&exclude_prop.2==1) | (exclude_poly.2==0&exclude_allele.2==1) |
                       (exclude_poly.3==0&exclude_prop.3==1) | (exclude_poly.3==0&exclude_allele.3==1))

nrow(prop_rem)

#make sure marker is prop excluded in at least 1 dataset
prop_rem %<>% mutate(exclpr_data1=case_when(exclude_poly==0&exclude_prop==1 ~ 1, T ~ 0)) %>%
              mutate(exclpr_data2=case_when(exclude_poly.2==0&exclude_prop.2==1 ~ 1, T ~ 0)) %>%
              mutate(exclpr_data3=case_when(exclude_poly.3==0&exclude_prop.3==1 ~ 1, T ~ 0)) %>%
  mutate(sompr=exclpr_data1+exclpr_data2+exclpr_data3) 

#note if marker retained in a data set
prop_rem %<>% mutate(excl_data1=case_when((exclude_match+exclude_poly+exclude_na+exclude_prop+exclude_allele+exclude_estmap)==0 ~ 0, T ~ 1)) %>%
             mutate(excl_data2=case_when((exclude_match.2+exclude_poly.2+exclude_na.2+exclude_prop.2+exclude_allele.2+exclude_estmap.2)==0 ~ 0, T ~ 1)) %>%
             mutate(excl_data3=case_when((exclude_match.3+exclude_poly.3+exclude_na.3+exclude_prop.3+exclude_allele.3+exclude_estmap.3)==0 ~ 0, T ~ 1)) %>%
  mutate(som=excl_data1+excl_data2+excl_data3) 

#number of markers ok in a least 1 dataset
prop_rem %>% filter(som<3) %>% nrow()

#number of markers always removed
prop_rem %<>% mutate(notpoly=case_when(som==3 & (exclude_poly+exclude_poly.2+exclude_poly.3==2) ~ 0,
                                       som==3 & is.na(exclude_poly.2)==TRUE & (exclude_poly+exclude_poly.3==1) ~ 0,
                                       som==3 & (exclude_poly+exclude_poly.2+exclude_poly.3==0) ~ 2,
                                       som==3 & is.na(exclude_poly.2)==TRUE & (exclude_poly+exclude_poly.3==0) ~ 2,
                                       som==3 & (exclude_poly+exclude_poly.2+exclude_poly.3<2) & (exclude_poly+exclude_poly.2+exclude_poly.3>0) ~ 1,
                                       som==3 & is.na(exclude_poly.2)==TRUE &  (exclude_poly+exclude_poly.3<1) & (exclude_poly+exclude_poly.3>0) ~ 1))

prop_rem %>% group_by(notpoly) %>% summarise(n=n())
```

## Are markers excluded by mark_estmap okay in other crosses ? 

```{r}
#keep all markers with non mendelian prop in at least 1 dataset
map_rem <- tab_alldata %>% filter(exclude_estmap==1 |
                       exclude_estmap.2==1 |
                       exclude_estmap.3==1)

map_rem %<>% mutate(som=case_when(exclude_match+exclude_poly+exclude_na+exclude_prop+exclude_allele+exclude_estmap==0 ~ 0,
                                 exclude_match.2+exclude_poly.2+exclude_na.2+exclude_prop.2+exclude_allele.2+exclude_estmap.2==0 ~ 0,
                                 exclude_match.3+exclude_poly.3+exclude_na.3+exclude_prop.3+exclude_allele.3+exclude_estmap.3==0 ~ 0,))
```

5/9 are okay in another cross.