data2.Rmd

---
title: "data2.Rmd"
author: "Marie Bourdon"
date: "01/06/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(stuart)
library(magrittr)
library(readr)
library(stringr)
library(qtl)
library(ggplot2)
source("../files/QTL_plot.R")
```

## Load
```{r}
genos <- read_csv("geno_data2.csv",show_col_types = FALSE) #genotypes of F2s
annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv")) #annotation file for miniMUGA
phenos <- read_csv("pheno_data2.csv",show_col_types = FALSE) #phenotypes of F2s
parents <- read_csv("parents_data2.csv",show_col_types = FALSE) #genotypes of parental strains (genotyped with F2s)
strns_ref <- read_csv("ref_geno_data2.csv",show_col_types = FALSE) #reference genotypes of parental strains
```

```{r,cache=TRUE}
tab <- tab_mark(genos,annot_mini,"cM_cox")
```

## Before: creation of Rqtl csv file 

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

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

# create rqtl csv file
cross_before <- write_rqtl(geno=genos,pheno=phenos,tab=tab_before,ref=strns_ref,par1="StrainA",par2="StrainB",prefix="F2-",pos="cM_cox",path="cluster/cross_before.csv")

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


load("cluster/newmap_before.rda")
plotMap(cross_before,newmap_before,shift=TRUE)
```

### Before: plot genome scan

```{r before_scan}
# load rda with perm
load("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
```

## create file for parental strains genotyped

```{r}
#our genotypes

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


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

#summary
summary(strains)
```


## Use of stuart's functions

```{r}
tab2 <- mark_match(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="StrainA",par2="StrainB")
```



### estmap

```{r}
# create rqtl csv file
write_rqtl(geno=genos,pheno=phenos,tab=tab2,ref=strains,par1="StrainA",par2="StrainB",prefix="F2-",pos="cM_cox",path="cluster/cross_after.csv")

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


load("cluster/newmap_after.rda")
plotMap(cross_after,newmap_after,shift=TRUE)

tab2 <- mark_estmap(tab=tab2,map=newmap_after,annot=annot_mini) #0 marker removed
```


```{r after_scan}
# load rda with perm
load("cluster/after_1000p.rda")

# calculate thresholds
threshold_after <- summary(after_1000p,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_after, step=2.0, off.end=0.0, 
                                 error.prob=1.0E-4, map.function="haldane", stepwidth="fixed")


pheno_after <- scanone(cross=cross_after, 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) +
    theme(legend.position = "none") +
    ggtitle("")
pheno_after_plot
```

## Number of markers kept after each function

```{r barplot}
# TO DO WHEN ESTMAP OK
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,functions_df)
```