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)

cross_before_data2 <- cross_before
newmap_before_data2 <- newmap_before
plotmap_before_data2 <- ~plotMap(cross_before_data2,newmap_before_data2,shift=TRUE,main='',ylab='')
```

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

### Distribution of calculated distance between markers

```{r}
## BEFORE
#initialize variables
mark <- c()
chr <- c()
pos <- c()
place <- c()
previous <- c()
follow <- c()
kn_previous <- c()
kn_follow <- c()

#get information in newmap
for(i in names(newmap_before)){
  mark <- c(mark,names(newmap_before[[i]]))
  chr <- c(chr,rep(i,times=length(newmap_before[[i]])))
  pos <- c(pos,unname(newmap_before[[i]]))
  place <- c(place,"first",rep("middle",times=(length(newmap_before[[i]])-2)),"last")
  prev <- c(NA,unname(newmap_before[[i]])[1:length(newmap_before[[i]])-1])
  previous <- c(previous,prev)
  fol <- c(unname(newmap_before[[i]])[2:length(newmap_before[[i]])],NA)
  follow <- c(follow,fol)
}

annot <- annot_mini %>% filter(marker %in% mark)
kn_pos <- annot$cM_cox
kn_prev <- c(NA, annot[1:(nrow(annot) - 1), "cM_cox"])
kn_previous <- c(kn_previous, kn_prev)
kn_fol <- c(annot[2:nrow(annot), "cM_cox"], NA)
kn_follow <- c(kn_follow, kn_fol)

#create tab with positions
rec_ratios <- tibble(marker = mark,
                  chr = chr,
                  place = place,
                  pos = pos,
                  previous = previous,
                  prev_dif = pos-previous,
                  kn_pos = kn_pos,
                  kn_previous = kn_previous,
                  kn_prev_dif = kn_pos - kn_previous)

rec_ratios <- rec_ratios %>% mutate(kn_prev_dif = case_when(is.na(previous) == TRUE ~ NA_real_, T ~ kn_prev_dif))


rec_ratios <- rec_ratios %>% mutate(rat_prev = prev_dif/kn_prev_dif)


#remove if dist < 1cM
rec_ratios %<>% filter(!prev_dif<1 & !kn_prev_dif<1) 

#mean sd
rec_ratios %>% summarise(mean=mean(rat_prev,na.rm=TRUE),
                       sd=sd(rat_prev,na.rm=TRUE),
                       max=max(rat_prev,na.rm=TRUE)) 


rec_ratios_before_data2 <- rec_ratios %>% ggplot(aes(x=rat_prev)) +
  geom_histogram(binwidth = .1,alpha=0.4, position="identity",fill="#990000") +
  scale_x_log10(limits=c(0.2,1000)) +
  labs(x="Ratio between the calculated and the known distances",
       y="Marker pairs",
       fill="",
       title="Distances between adjacent markers") +
  geom_vline(xintercept = 5,linetype="dashed") +
  ggpubr::theme_classic2() +
  theme(plot.title = element_text(face="plain",size=14,hjust=0.5)) 

rec_ratios_before_data2

#percentage of markers in each group
rec_ratios %>% mutate(group_rat=case_when(rat_prev<5 ~ 0,
                                          rat_prev >=5 ~ 1)) %>%
  group_by(group_rat) %>% summarise(n=n(),p=n()/nrow(rec_ratios))

#meansd of peak <5
rec_ratios %>% filter(rat_prev<5) %>% summarise(mean=mean(rat_prev),
                                                sd=sd(rat_prev))

#max ratio
rec_ratios %>% summarise(max=max(rat_prev))

rm(mark,chr,pos,place,previous,follow,kn_previous,kn_follow,fol,kn_fol,kn_pos,kn_prev,i,prev)
```

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

cross_after_data2 <- cross_after
newmap_after_data2 <- newmap_after
plotmap_after_data2 <- ~plotMap(cross_after_data2,newmap_after_data2,shift=TRUE,main="")


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

### map length

```{r}
#calc map
length_calc <- newmap_after[["1"]][[length(newmap_after[["1"]])]] + newmap_after[["2"]][[length(newmap_after[["2"]])]] +
  newmap_after[["3"]][[length(newmap_after[["3"]])]] + newmap_after[["4"]][[length(newmap_after[["4"]])]] +
  newmap_after[["5"]][[length(newmap_after[["5"]])]] + newmap_after[["6"]][[length(newmap_after[["6"]])]] +
  newmap_after[["7"]][[length(newmap_after[["7"]])]] + newmap_after[["8"]][[length(newmap_after[["8"]])]] +
  newmap_after[["9"]][[length(newmap_after[["9"]])]] + newmap_after[["10"]][[length(newmap_after[["10"]])]] +
  newmap_after[["11"]][[length(newmap_after[["11"]])]] + newmap_after[["12"]][[length(newmap_after[["12"]])]] +
  newmap_after[["13"]][[length(newmap_after[["13"]])]] + newmap_after[["14"]][[length(newmap_after[["14"]])]] +
  newmap_after[["15"]][[length(newmap_after[["15"]])]] + newmap_after[["16"]][[length(newmap_after[["16"]])]] +
  newmap_after[["17"]][[length(newmap_after[["17"]])]] + newmap_after[["18"]][[length(newmap_after[["18"]])]] +
  newmap_after[["19"]][[length(newmap_after[["19"]])]] + newmap_after[["X"]][[length(newmap_after[["X"]])]]

#known map
lengh_known <- cross_after[["geno"]][["1"]][["map"]][[length(cross_after[["geno"]][["1"]][["map"]])]] +
  cross_after[["geno"]][["2"]][["map"]][[length(cross_after[["geno"]][["2"]][["map"]])]] +
  cross_after[["geno"]][["3"]][["map"]][[length(cross_after[["geno"]][["3"]][["map"]])]] + 
  cross_after[["geno"]][["4"]][["map"]][[length(cross_after[["geno"]][["4"]][["map"]])]] +
  cross_after[["geno"]][["5"]][["map"]][[length(cross_after[["geno"]][["5"]][["map"]])]] + 
  cross_after[["geno"]][["6"]][["map"]][[length(cross_after[["geno"]][["6"]][["map"]])]] +
  cross_after[["geno"]][["7"]][["map"]][[length(cross_after[["geno"]][["7"]][["map"]])]] + 
  cross_after[["geno"]][["8"]][["map"]][[length(cross_after[["geno"]][["8"]][["map"]])]] +
  cross_after[["geno"]][["9"]][["map"]][[length(cross_after[["geno"]][["9"]][["map"]])]] + 
  cross_after[["geno"]][["10"]][["map"]][[length(cross_after[["geno"]][["10"]][["map"]])]] +
  cross_after[["geno"]][["11"]][["map"]][[length(cross_after[["geno"]][["11"]][["map"]])]] + 
  cross_after[["geno"]][["12"]][["map"]][[length(cross_after[["geno"]][["12"]][["map"]])]] +
  cross_after[["geno"]][["13"]][["map"]][[length(cross_after[["geno"]][["13"]][["map"]])]] + 
  cross_after[["geno"]][["14"]][["map"]][[length(cross_after[["geno"]][["14"]][["map"]])]] +
  cross_after[["geno"]][["15"]][["map"]][[length(cross_after[["geno"]][["15"]][["map"]])]] + 
  cross_after[["geno"]][["16"]][["map"]][[length(cross_after[["geno"]][["16"]][["map"]])]] +
  cross_after[["geno"]][["17"]][["map"]][[length(cross_after[["geno"]][["17"]][["map"]])]] + 
  cross_after[["geno"]][["18"]][["map"]][[length(cross_after[["geno"]][["18"]][["map"]])]] +
  cross_after[["geno"]][["19"]][["map"]][[length(cross_after[["geno"]][["19"]][["map"]])]] + 
  cross_after[["geno"]][["X"]][["map"]][[length(cross_after[["geno"]][["X"]][["map"]])]]

#ratio 
length_rat <- length_calc/lengh_known
length_rat
```

```{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_data2 <- 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("Dataset 2: genome scan") +
  theme(plot.title = element_text(face="plain",size=14,hjust=0.5)) 
pheno_after_plot_data2
```

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

### Distribution of calculated distance between markers

```{r}
## AFTER
#initialize variables
mark <- c()
chr <- c()
pos <- c()
place <- c()
previous <- c()
follow <- c()
kn_previous <- c()
kn_follow <- c()

#get information in newmap
for(i in names(newmap_after)){
  mark <- c(mark,names(newmap_after[[i]]))
  chr <- c(chr,rep(i,times=length(newmap_after[[i]])))
  pos <- c(pos,unname(newmap_after[[i]]))
  place <- c(place,"first",rep("middle",times=(length(newmap_after[[i]])-2)),"last")
  prev <- c(NA,unname(newmap_after[[i]])[1:length(newmap_after[[i]])-1])
  previous <- c(previous,prev)
  fol <- c(unname(newmap_after[[i]])[2:length(newmap_after[[i]])],NA)
  follow <- c(follow,fol)
}

annot <- annot_mini %>% filter(marker %in% mark)
kn_pos <- annot$cM_cox
kn_prev <- c(NA, annot[1:(nrow(annot) - 1), "cM_cox"])
kn_previous <- c(kn_previous, kn_prev)
kn_fol <- c(annot[2:nrow(annot), "cM_cox"], NA)
kn_follow <- c(kn_follow, kn_fol)

#create tab with positions
rec_ratios <- tibble(marker = mark,
                  chr = chr,
                  place = place,
                  pos = pos,
                  previous = previous,
                  prev_dif = pos-previous,
                  kn_pos = kn_pos,
                  kn_previous = kn_previous,
                  kn_prev_dif = kn_pos - kn_previous)

rec_ratios <- rec_ratios %>% mutate(kn_prev_dif = case_when(is.na(previous) == TRUE ~ NA_real_, T ~ kn_prev_dif))

rec_ratios <- rec_ratios %>% mutate(rat_prev = prev_dif/kn_prev_dif)

#remove if dist < 1cM
rec_ratios %<>% filter(!prev_dif<1 & !kn_prev_dif<1) 

#mean sd
rec_ratios %>% summarise(mean=mean(rat_prev,na.rm=TRUE),
                       sd=sd(rat_prev,na.rm=TRUE)) 


rec_ratios_after_data2 <- rec_ratios %>% ggplot(aes(x=rat_prev)) +
  geom_histogram(binwidth = .1,alpha=0.4, position="identity",fill="#2171b5") +
  scale_x_log10(limits=c(0.2,1000)) +
  labs(x="Ratio between the calculated and the known distances",
       y="Marker pairs",
       fill="",
       title="Distances between adjacent markers") +
  ggpubr::theme_classic2() +
  theme(plot.title = element_text(face="plain",size=14,hjust=0.5)) 

rec_ratios_after_data2

rm(mark,chr,pos,place,previous,follow,kn_previous,kn_follow,fol,kn_fol,kn_pos,kn_prev,i,prev)
```

## Narrow peaks

Investigation of high lod score peaks 


```{r before_ann}
chrs <- c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "X")
ann_dat_text<-data.frame(
    chr=factor(chrs,
               levels=chrs),
    lod=c(rep(NA,10),8,rep(NA,9)),
    label=c(rep(NA,10),"p3",rep(NA,9)),
    x=c(rep(NA,10),22,rep(NA,9))
)

pheno_before_annot_data2 <- pheno_before_plot +  geom_text(
    # the new dataframe for annotating text
    data = ann_dat_text,
    mapping = aes(x = x,
                  y = lod,
                  label = label,
                  color="blue")
  ) +
  ggtitle("Dataset 2: genome scan") +
  theme(plot.title = element_text(face="plain",size=14,hjust=0.5)) 
pheno_before_annot_data2

rm(ann_dat_text)
rm(chrs)
```


### Peak 3
grid_scans
```{r peak3_zoom}
peak3 <- 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",
         size=0.6,
         strip.pos="bottom",
         chr="11",
         rug=TRUE) +
    theme(legend.position = "none",
          strip.background = element_blank(),
          strip.text.x = element_blank()) +
    xlab("Position on chromosome 11 (cM)") +
    xlim(c(11,31)) +
    ggtitle("Peak p3") +
    theme(plot.title = element_text(face="plain",size=14,hjust=0.5)) 
peak3
```

1 peak on chromosome 11 on 1 pseudomarker : c11.loc18, postionned between SNT111392585 and mJAX00308021.

Here are the infos on genotype counts for these markers:

```{r summary_geno_peak3}
peak3_tab <- tab_before %>% filter(marker %in% c("SNT111392585","mJAX00308021")) %>% select(marker:n_NA)
peak3_tab 
```

For SNT111392585, all individuals except 1 are homozygous so this marker should be removed. The proportions for mJAX00308021 seem correct.

Graph:

```{r geno_plot_peak7}
phenotypes <- cross_before[["pheno"]]
map <- cross_before[["geno"]][["11"]][["map"]] 
map <- tibble(marker=names(map),pos=map)
genotypes <- cross_before[["geno"]][["11"]][["data"]]
genotypes <- as_tibble(genotypes)
phenogeno <- cbind(phenotypes,genotypes)
phenogeno %<>% pivot_longer(mbackupUNC110000218:gUNC20538837,names_to="marker",values_to="genotype")
pgmap <- full_join(phenogeno,map,by="marker")

geno_plot3 <- pgmap %>% filter(pos > 15 & pos < 25) %>%
  filter(id %in% sample(phenotypes$id,10)) %>%
  ggplot(aes(x=pos,y=as.factor(id))) +
  geom_point(aes(color=as.factor(genotype))) +
  coord_cartesian(ylim = c(1, 10), expand = TRUE, clip = "off") +
  annotate(geom="text",y=-1,size=3,
           x = map %>% filter(pos > 15 & pos < 25) %>% pull(pos),
           label = map %>% filter(pos > 15 & pos < 25) %>% pull(marker),
           angle=90) +
  labs(x="Position (cM)",y="Individual",color="Genotype") +
  theme_bw() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"),
        axis.title.x = element_text(margin = margin(t = 50)))
geno_plot3

rm(pgmap,phenotypes,map,genotypes,phenogeno)
```

The two markers before the pseudomarker have an excess of homozygous.



## Phenotype distributions

```{r pheno_distrib}
pheno_data2 <- phenos %>% ggplot(aes(x=Pheno)) +
  geom_histogram(binwidth=0.1) +
  ggpubr::theme_classic2() +
  labs(y="F2 individuals", x="Quantitative phenotype",title="Dataset 2") +
  theme(plot.title = element_text(hjust=0.5,size=14)) 
pheno_data2
```

## Missing genotypes

```{r}
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
```

## Genotype proportions

```{r}
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) +
  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))

prop_plot
```


```{r save_narrow}
tab2_data2 <- tab2
save(pheno_before_annot_data2,pheno_data2,peak3,peak3_tab,tab2_data2,
     plotmap_before_data2,plotmap_after_data2,
     rec_ratios_before_data2,rec_ratios_after_data2,
     cross_before_data2,newmap_before_data2,cross_after_data2,newmap_after_data2,
     pheno_after_plot_data2,
     file="data2.rda")
rm(tab2_data2)
```