article_figures.Rmd

---
title: "article_figures.Rmd"
author: "Marie Bourdon"
date: "July 2021"
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(stuart)

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

# Data load and use of stuart functions

## genos dataframe

Data frame from stuart with genotypes of 176 F2 individuals and parental strains.
```{r genos}
data(genos)
summary(genos)
```
## phenos dataframe

Data frame from stuart with phenotypes of 176 F2 individuals for a quantitative trait.
```{r phenos}
data(phenos)
summary(phenos)
```
## annotation file

Annotation file from K.Broman: https://kbroman.org/MUGAarrays/mini_annotations.html
```{r annot}
annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv"))
summary(annot_mini)
```
## parental strains genotypes

2 genotypes tables will be used:

* One table with our genotypes of the strains used in the cross: "strains"

* One table with reference genotypes: "strns_ref"

```{r strains}
#our genotypes

#create tibble with individivual names
parental_strains <- tibble::tibble(parent1 = c("StrainsA_1","StrainsA_2"),
                                   parent2 = c("StrainsB_1","StrainsB_2"))


#create data frame with geno_strains
strains <- geno_strains(annot=annot_mini,geno=genos,
                        strn=parental_strains,cols=c("chr","cM_cox"))
rm(parental_strains)

#remove genotypes of parental strains from genos data frame
genos %<>% filter(!Sample.ID %in% c("StrainsA_1","StrainsA_2","StrainsB_1","StrainsB_2"))

#summary
summary(strains)
```

```{r strns_ref}
#reference genotypes
#load parental strains genotype data from Neogen
strns_ref <- read.csv("files/ref_genotypes.csv")

#join with annotation file from miniMUGA
strns_ref <- strns_ref %>% select(name,parent1,parent2) %>% right_join(annot_mini,.,by=c("marker"="name")) %>% select(c(marker,chr,cM_cox,parent1,parent2))

#summary
summary(strns_ref)
```

## Use of stuart's functions
```{r}
data(stuart_tab)
summary(stuart_tab)

tab2 <- mark_match(stuart_tab,ref=strns_ref)
tab2 <- mark_poly(tab2)
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1,homo1X=c(0.1,1),homo2X=c(0.1,1),heteroX=c(0.1,1))
tab2 <- mark_allele(tab2,ref=strns_ref,cross="F2",par1="parent1",par2="parent2")

data(stuart_newmap)
tab2 <- mark_estmap(tab=tab2,map=stuart_newmap,annot=annot_mini)
```

# 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() +
  theme_classic() +
  labs(title="Proportion of missing genotyped",
       x="Proportion of NA",y="Number of markers") +
  theme(
    aspect.ratio=0.8,
    plot.title = element_text(hjust = 0.4,face="bold",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) +
  #geom_text(aes(label=ifelse(exclude_prop=="1",SNP.Name,'')),hjust=0, vjust=0,size=2) +
  labs(title="Exclusion of markers with mark_prop()",
       x="Proportion of homozygous individuals AA",
       y="Proportion of homozygous individuals BB",
       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,face="bold",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}
grid <- plot_grid(na_plot,prop_plot,
                  ncol=1,
                  labels = c('A', 'B'), 
                  label_size = 20,
                  byrow = TRUE
                  )

grid

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

rm(na_plot,prop_plot)
```

## Graph before after

### Before: creation of Rqtl csv file 

```{r cross_before}
# filter at minima: remove non polymorphic and NA>0.5 
tab_before <- mark_na(stuart_tab)
tab_before <- mark_poly(tab_before)

# create rqtl csv file
write_rqtl(geno=genos,pheno=phenos,tab=tab_before,ref=strns_ref,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path="files/cluster/cross_before.csv")
```

### Before: newmap and permutation

These objects were produced on our cluster with the following script: /files/cluster/before_after.R

### Before: plot estimated map

```{r before_map}
# import cross
cross_before <- read.cross(format="csv",file="files/cluster/cross_before.csv",
                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

# load rda with newmap
load("files/cluster/newmap_before.rda")

# plot
plotMap(cross_before,newmap_before,shift=TRUE)
plotmap_before <- ~plotMap(cross_before,newmap_before,shift=TRUE,main="Before stuart", ylab='')
```

### Before: plot genome scan

```{r before_scan}
# load rda with perm
load("files/cluster/before_1000p.rda")
# load("files/before_1000p.rda")


# calculate thresholds
threshold_before <- summary(before_1000p,alpha=c(0.05,0.1,0.63)) #donne lod score pour risque 0.05, 0.1, 0.63

# scanone
cross_before <- calc.genoprob(cross_before, step=2.0, off.end=0.0, 
                                 error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")


pheno_before <- scanone(cross=cross_before, chr=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X", "Y"), pheno.col="Pheno", model="normal", method="em")
summary(pheno_before)

# Plot
pheno_before_plot <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_before[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom") +
    theme(legend.position = "none") +
    ggtitle("")
pheno_before_plot
```


```{r before_scan}
pheno_before_zoom <- qtl_plot(pheno_before,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_before[1:3]),
         ylab="LOD score",
         chrs = "13",
         size=0.6) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position (cM)") +
    ggtitle("")
pheno_before_zoom
```

### Before figure

```{r grid_before}
grid_before <- plot_grid(as_grob(plotmap_before),pheno_before_plot,pheno_before_zoom,ncol=1,labels=c("A","B","C"),label_size=20)

ggsave("figures/fig1.pdf",grid_before,width=7,height = 15)
ggsave("figures/fig1.png",grid_before,width=7,height = 15)
```


### After: creation of Rqtl csv file 

```{r cross_after}
# filter with stuart functions: use the good data for parental strains (strains df)
tab2 <- mark_match(stuart_tab,ref=strains)
tab2 <- mark_poly(tab2)
tab2 <- mark_na(tab2)
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1,homo1X=c(0.1,1),homo2X=c(0.1,1),heteroX=c(0.1,1))
tab2 <- mark_allele(tab2,ref=strains,cross="F2",par1="parent1",par2="parent2")


# create rqtl csv file
write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path="files/cluster/cross_after.csv")
```

### After: newmap 

These objects were produced on our cluster with the following script: /files/cluster/before_after.R

### After: plot estimated map

```{r after_map}
# import cross
cross_after <- read.cross(format="csv",file="files/cluster/cross_after.csv",
                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

# load rda with newmap
load("files/cluster/newmap_after.rda")

# plot
plotMap(cross_after,newmap_after,shift=TRUE)
plotmap_after <- ~plotMap(cross_after,newmap_after,shift=TRUE,main="After stuart")
```
### Remove last problematic markers with mark_estmap

```{r after_estmap}
#filter with mark_estmap
tab2 <- mark_estmap(tab2,newmap_after,annot_mini)

# create rqtl csv file
write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox",path="files/cluster2/cross_after2.csv")
```

### After: plot estimated map 2

```{r after_map2}
# import cross
cross_after2 <- read.cross(format="csv",file="files/cluster2/cross_after2.csv",
                              genotypes=c("0","1","2"),na.strings=c("NA"), convertXdata=TRUE)

# load rda with newmap
load("files/cluster2/newmap_after2.rda")

# plot
plotMap(cross_after2,newmap_after2,shift=TRUE)
plotmap_after <- ~plotMap(cross_after2,newmap_after2,shift=TRUE,main="After stuart")
```

### After: plot genome scan

```{r after_scan}
# load rda with perm
load("files/cluster2/after_1000p2.rda")

# calculate thresholds
threshold_after <- summary(after_1000p2,alpha=c(0.05,0.1,0.63)) #donne lod score pour risque 0.05, 0.1, 0.63

# scanone
cross_after <- calc.genoprob(cross_after2, step=2.0, off.end=0.0, 
                                 error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")


pheno_after <- scanone(cross=cross_after2, chr=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X", "Y"), pheno.col="Pheno", model="normal", method="em")
summary(pheno_after)

# Plot
pheno_after_plot <- qtl_plot(pheno_after,lod=data.frame(group = c("alpha=0.05", "alpha=0.1","alpha=0.63"),
                 lod = threshold_after[1:3]),
         ylab="LOD score",
         title="QTL mapping",
         x.label = element_blank(),
         size=0.6,
         strip.pos="bottom") +
    theme(legend.position = "none") +
    ggtitle("")
pheno_after_plot
```
### Grid after

```{r grid fig1, fig.height = 7, fig.width = 13, fig.align = "center"}
grid_after <- plot_grid(as_grob(plotmap_after),pheno_after_plot,ncol=1,labels=c("A","B"),label_size=20)

ggsave("figures/fig3.pdf",grid_after,width=7,height = 10)
ggsave("figures/fig3.png",grid_after,width=7,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()
```

## 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)
```

```{r dif_table}
table_dif <- dif %>% filter(dif==1) %>% select(marker,parent1_ref=parent1.y,parent1_geno=parent1.x,parent2_ref=parent2.y,parent2_geno=parent2.x) %>% head()
knitr::kable(table_dif)
```

## Number of markers kept after each function

```{r barplot}
none <- tab2 %>% nrow()
match <- tab2 %>% filter(exclude_match==0) %>% nrow()
allele <- tab2 %>% filter(exclude_match==0&exclude_allele==0) %>% nrow()
naf <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0) %>% nrow()
poly <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0&exclude_poly==0) %>% nrow()
prop <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0&exclude_poly==0&exclude_prop==0) %>% nrow()
estmap <- tab2 %>% filter(exclude_match==0&exclude_allele==0&exclude_na==0&exclude_poly==0&exclude_prop==0&exclude_estmap==0) %>% nrow()

functions_df <- tibble(fct=c("none","match","allele","na","poly","prop","estmap"),
                       markers=c(none,match,allele,naf,poly,prop,estmap))

functions_plot <- functions_df %>% ggplot(aes(x=markers,y=fct)) +
  geom_bar(stat="identity",width=0.6) +
  geom_text(aes(label=markers), hjust=1.3, color="white", size=3.5) +
  scale_y_discrete(limits=c("estmap","prop","poly", "na", "allele","match","none")) +
  theme(aspect.ratio=0.7) +
  labs(title="Number of markers kept after each step",
       x="Number of markers",
       y="Function used") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.4,face="bold",size=14))

functions_plot
rm(none,allele,match,poly,prop,barplot_df)
```

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