---
title: "Compute and Plot Linkage Decay"
author: "Samuel Decker"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "j:/camelina/camelina_r/DISSERTATION")
```

## Load the required libraries

# most of this code was taken from this tutorial:
https://www.biostars.org/p/300381/#300423

```{R Library Loading, include = FALSE}
library(dplyr)
library(stringr)
library(ggplot2)
```

## Load the data

```{R Load the Data}
#### load the workspace image.
#load("Compute and plot Linkage Decay.RData")

#### or load the data from scratch:
# read.delim is not fast, but it does work for this case.  
all_markers_summary <- read.delim("all_markers_LD.ld.summary", sep="", header=T, check.names = F, stringsAsFactors = F)
#ld_summary <- summary(all_markers_summary)

# add column names to the summary data
colnames(all_markers_summary) <- c("distance", "r2")
```

## Save the environment
```{R save the environment}
save.image("Compute and plot Linkage Decay.RData")
```

## Compute the linkage decay

```{R Compute Linkage Decay}
#### The tutorial sets the intervals to 10kb. I want to do 5, I am trying to maximize the granularity of this analysis
all_markers_summary$distance_interval <- cut(all_markers_summary$distance, 
                                         breaks=seq(from=min(all_markers_summary)-1,
                                                    to=max(all_markers_summary$distance)+1,
                                                    by=5000))

#### compute the mean and/or median r2 within blocks
LD_data <- all_markers_summary %>% group_by(distance_interval) %>% summarise(mean=mean(r2),median=median(r2))

#### helper step to get the mid points of the distance intervals for ploting
LD_data <- LD_data %>% mutate(start=as.integer(str_extract(str_replace_all(distance_interval,"[\\(\\)\\[\\]]",""),"^[0-9-e+.]+")),
                        end=as.integer(str_extract(str_replace_all(distance_interval,"[\\(\\)\\[\\]]",""),"[0-9-e+.]+$")),
                        mid=start+((end-start)/2))
```

## Plot the Linkage Decay

```{R plot linkage decay}
#### plot the data:
LD_decayplot <- ggplot()+
  geom_point(data=LD_data,aes(x=start,y=mean),size=1.0,colour="black")+
  geom_line(data=LD_data,aes(x=start,y=mean),linewidth=1.3,alpha=0.5,colour="black")+
  labs(x="Distance Between Markers (kilobases)",y=expression(LD~(r^{2})))+
  scale_x_continuous(breaks=c(1*10^4,5*10^4,1*10^5,1.5*10^5,2*10^5,2.5*10^5),labels=c("10","50","100","150","200","250"))+
  ggtitle("Linkage Decay") + theme(axis.text=element_text(size=14), axis.title=element_text(size=20), 
                                  plot.title=element_text(size=24), panel.background=element_rect("lightblue", "darkblue"),
                                  axis.ticks=element_line(linewidth=1.5))
LD_decayplot + ggtitle("Camelina sativa Linkage Decay") +
  theme(plot.title = element_text(lineheight = 1.5, face="bold"))
```

## Save the Final environment

```{R save the environment again}
#save.image("Compute and plot Linkage Decay.RData")
```