---
title: "43-intron-exon-methylation"
output: html_document
---

Lets see if we can get methylation information for intron and exon component of each gene

should get feature tracks

http://eagle.fish.washington.edu/Cvirg_tracks/C_virginica-3.0_intron.bed
http://eagle.fish.washington.edu/Cvirg_tracks/C_virginica-3.0_Gnomon_exon.bed

```{bash}
cd ../genome-features
wget http://eagle.fish.washington.edu/Cvirg_tracks/C_virginica-3.0_intron.bed
```
```{bash}
cd ../genome-features
wget http://eagle.fish.washington.edu/Cvirg_tracks/C_virginica-3.0_Gnomon_exon.bed
```

```{bash}
head ../genome-features/C_virginica-3.0_Gnomon_exon.bed
```
```{bash}
head ../genome-features/C_virginica-3.0_intron.bed
```

Can we get and exon and intron bed with gene ID info tied to it?

```{bash}
#test

NAME=5F

/home/shared/bedtools2/bin/intersectBed \
-wb \
-a ../genome-features/C_virginica-3.0_Gnomon_exon.bed \
-b ../genome-features/C_virginica-3.0_Gnomon_genes.bed \
| awk -v name=$NAME -v OFS="\t" '{ print $1, $2, $3, $7, name}' \
| head

```
```{bash}
/home/shared/bedtools2/bin/intersectBed \
-wb \
-a ../genome-features/C_virginica-3.0_Gnomon_exon.bed \
-b ../genome-features/C_virginica-3.0_Gnomon_genes.bed \
| awk -v OFS="\t" '{ print $1, $2, $3, $7}' \
| head
```
```{bash}
/home/shared/bedtools2/bin/intersectBed \
-wb \
-a ../genome-features/C_virginica-3.0_Gnomon_exon.bed \
-b ../genome-features/C_virginica-3.0_Gnomon_genes.bed \
| awk -v OFS="\t" '{ print $1, $2, $3, $7}' \
> ../genome-features/C_virginica-3.0_Gnomon_exon-geneID.bed
```

```{bash}
/home/shared/bedtools2/bin/intersectBed \
-wb \
-a ../genome-features/C_virginica-3.0_intron.bed \
-b ../genome-features/C_virginica-3.0_Gnomon_genes.bed \
| awk -v OFS="\t" '{ print $1, $2, $3, $7}' \
> ../genome-features/C_virginica-3.0_Gnomon_intron-geneID.bed
```






```{bash}
 cd ../data/big/
 
 FILES=$(ls *bedgraph)
 
 cd -
 
 for file in ${FILES}
 do
    NAME=$(echo ${file} | awk -F "_" '{print $1}')
    echo ${NAME}
   /home/shared/bedtools2/bin/intersectBed \
   -wb \
   -a ../data/big/${NAME}_R1_val_1_10x.bedgraph \
   -b ../genome-features/C_virginica-3.0_Gnomon_exon-geneID.bed \
   | awk -v name=$NAME -v OFS="\t" '{ print $0, name}' \
   > ../output/43-exon-intron-methylation/${NAME}_mExon.out
 done  


```




```{bash}
head ../output/43-exon-intron-methylation/36F_mExon.out
```


```{bash}
 cd ../data/big/
 
 FILES=$(ls *bedgraph)
 
 cd -
 
 for file in ${FILES}
 do
    NAME=$(echo ${file} | awk -F "_" '{print $1}')
    echo ${NAME}
   /home/shared/bedtools2/bin/intersectBed \
   -wb \
   -a ../data/big/${NAME}_R1_val_1_10x.bedgraph \
   -b ../genome-features/C_virginica-3.0_Gnomon_intron-geneID.bed \
   | awk -v name=$NAME -v OFS="\t" '{ print $0, name}' \
   > ../output/43-exon-intron-methylation/${NAME}_mIntron.out
 done  


```


```{bash}
head -2 ../output/43-exon-intron-methylation/36F_mIntron.out
```


```{bash}
cat ../output/43-exon-intron-methylation/*_mExon.out > ../output/43-exon-intron-methylation/exon-meth_all-samples.out
```

```{bash}
cat ../output/43-exon-intron-methylation/*_mIntron.out > ../output/43-exon-intron-methylation/intron-meth_all-samples.out
```














```{r}
exon_meth <- read.delim("../output/43-exon-intron-methylation/exon-meth_all-samples.out", header = FALSE)
```

```{r}
intron_meth <- read.delim("../output/43-exon-intron-methylation/intron-meth_all-samples.out", header = FALSE)
```




```{r}
library(tidyverse)
```


```{r}
exon_meth
```






```{r}
ggplot() +
  geom_histogram(data = exon_meth, aes(x = V4), fill = "blue", alpha = 0.2) 
```


```{r}
em <- exon_meth %>%
   mutate(art = paste(V8, V9, sep = '_')) %>%
   group_by(art) %>%
   summarize(avg_meth = mean(V4))


```


```{r}
em
```


```{r}
int <- intron_meth %>%
   mutate(art = paste(V8, V9, sep = '_')) %>%
   group_by(art) %>%
   summarize(avg_meth = mean(V4))


```





```{r}
exint <- inner_join(em, int, by = "art") %>%
   separate(art, into = c("gene", "sample"), sep = "_") %>%
   separate(sample, into = c("number", "sex"), sep = -1)

```

```{r}
exint
```
```{r}
sumgene <- inner_join(em, int, by = "art") %>%
   separate(art, into = c("gene", "sample"), sep = "_") %>%
   separate(sample, into = c("number", "sex"), sep = -1) %>%
   group_by(gene) %>%
   summarize(avg_meth_exon = mean(avg_methy.x), avg_meth_intron = mean(avg_meth.y))

```



```{r}
library("RColorBrewer")
```

```{r}
ggplot(exint, aes(x = avg_meth.x, y = avg_meth.y)) +
   geom_point(alpha = 0.6, size = .2)
   
```



```{r, fig.width=20,fig.height=10}

ggplot(exint, aes(x = avg_meth.x, y = avg_meth.y, color = sex)) +
  geom_point(alpha = 1, size = .2) +
  facet_wrap(~number) +
  #geom_smooth(method = "lm", color = "black")
  stat_smooth(method = "glm", color = "black")
```
```{r, fig.width=20,fig.height=10}

ggplot(exint, aes(x = avg_meth.x, y = avg_meth.y)) +
  geom_point(alpha = 1, size = .2) +
  facet_wrap(~sex)
```




```{r}
library("ggExtra")
```


```{r}

p <- ggplot(exint, aes(x = avg_meth.x, y = avg_meth.y)) +
  geom_point(alpha = 0.2, size = 1) +
  facet_wrap(~sex) 
ggMarginal(p, type = "histogram")
```


```{r}
ggplot(exint) +
  geom_histogram(x = avg_meth.x) +
  geom_histogram(x = avg_meth.y) +
  facet_wrap(~sex)
```










#Creating exon and intron methylation file where libary is row

```{r}
em %>%
  separate(art, into = c("gene", "sample_id"), sep = '_') %>%
  spread(value = avg_meth, key = gene)
```



```{r}
em %>%
  separate(art, into = c("gene", "sample_id"), sep = '_') %>%
  spread(value = avg_meth, key = gene) %>%
  write.csv(., "../output/43-exon-intron-methylation/exon-methylation.csv")
```

```{r}
int %>%
  separate(art, into = c("gene", "sample_id"), sep = '_') %>%
  spread(value = avg_meth, key = gene) %>%
  write.csv(., "../output/43-exon-intron-methylation/intron-methylation.csv")
```




