---
title: "01-methykit"
output: html_document
---

# Done on RAVEN


```{r load packages, echo=TRUE, message=FALSE, warning=FALSE}
library(methylKit)
library(tidyverse)
```

Bismark code
```

# Directories and programs
bismark_dir="/gscratch/srlab/programs/Bismark-0.22.3"
bowtie2_dir="/gscratch/srlab/programs/bowtie2-2.3.4.1-linux-x86_64/"
samtools="/gscratch/srlab/programs/samtools-1.9/samtools"
reads_dir="/gscratch/srlab/sr320/data/cg-big/"
fastqc="/gscratch/srlab/programs/fastqc_v0.11.9/fastqc"
genome_folder="/gscratch/srlab/sr320/data/Cgig-genome/roslin_M/"


# 
# 
# ${bismark_dir}/bismark_genome_preparation \
# --verbose \
# --parallel 28 \
# --path_to_aligner ${bowtie2_dir} \
# ${genome_folder}




find ${reads_dir}*_R1.fq.gz \
| xargs basename -s _R1.fq.gz | xargs -I{} ${bismark_dir}/bismark \
--path_to_bowtie ${bowtie2_dir} \
-genome ${genome_folder} \
-p 4 \
-score_min L,0,-0.6 \
--non_directional \
-1 ${reads_dir}{}_R1.fq.gz \
-2 ${reads_dir}{}_R2.fq.gz

find *.bam | \
xargs basename -s .bam | \
xargs -I{} ${bismark_dir}/deduplicate_bismark \
--bam \
--paired \
{}.bam



${bismark_dir}/bismark_methylation_extractor \
--bedGraph \
--counts \
--comprehensive \
--merge_non_CpG \
--multicore 28 \
--buffer_size 75% \
*deduplicated.bam

# Sort files for methylkit and IGV

find *deduplicated.bam | \
xargs basename -s .bam | \
xargs -I{} ${samtools} \
sort --threads 28 {}.bam \
-o {}.sorted.bam

# Index sorted files for IGV
# The "-@ 16" below specifies number of CPU threads to use.

find *.sorted.bam | \
xargs basename -s .sorted.bam | \
xargs -I{} ${samtools} \
index -@ 28 {}.sorted.bam
```


## Code to wget Bismark files down from Roberts lab server

```{bash}
wget -r \
--no-check-certificate \
--no-directories --no-parent \
-P ../data/ \
-A "[0-9]*sorted.bam" \
https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/061021-big/
```

_Samples IDS_
"0501", "0502", "0503", "3501", "3502", "3503", ----SUS
"5201", "5202", "5203", "5901", "5902", "5903" <---- RES
                        
                        

```{r include=FALSE}
file.list=list('../data/0501_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/0502_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/0503_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/3501_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/3502_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/3503_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/5201_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/5202_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/5203_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/5901_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/5902_R1_bismark_bt2_pe.deduplicated.sorted.bam',
                 '../data/5903_R1_bismark_bt2_pe.deduplicated.sorted.bam'
)
```
                        
```{r eval=FALSE, include=FALSE}
myobj = processBismarkAln(location = file.list, sample.id = list("0501", "0502", "0503", "3501", "3502", "3503","5201", "5202", "5203", "5901", "5902", "5903"), assembly = "ROSLIN", read.context="CpG", mincov=2, treatment = c(0,0,0,0,0,0,1,1,1,1,1,1))
```                      
          
```{r, eval = FALSE}
save(myobj, file = "../analyses/myobj")

```


```{r include=FALSE}
load("../analyses/myobj")
```



          
```{r message=FALSE, warning=FALSE}
filtered.myobj=filterByCoverage(myobj,lo.count=5,lo.perc=NULL,
                                      hi.count=NULL,hi.perc=99)

meth_filter=unite(filtered.myobj, min.per.group=5L, destrand=TRUE)

clusterSamples(meth_filter, dist="correlation", method="ward", plot=TRUE)


PCASamples(meth_filter)

```    


          
```{r include=FALSE}
myDiff=calculateDiffMeth(meth_filter,mc.cores=24)

```

```{r}
# get hyper methylated bases
myDiff.hyper=getMethylDiff(myDiff,difference=50,qvalue=0.01,type="hyper")
#
# get hypo methylated bases
myDiff.hypo=getMethylDiff(myDiff,difference=50,qvalue=0.01,type="hypo")
#
#
# get all differentially methylated bases
myDiff=getMethylDiff(myDiff,difference=50,qvalue=0.01)

```

```{r}
write.table(myDiff, file = "../analyses/myDiff.tab", sep = "\t")
```
        
        
There were `r nrow(myDiff)` differentially methylated loci iddentifed. 
                        