---
title: "03-methylkit"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Running DML with BAMs


https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/060321-oly/1_ATCACG_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam

/Volumes/block/wd/061421

```{bash}
pwd
```


```{bash, eval=FALSE}
wget -r \
--no-directories --no-parent \
-P /Volumes/block/wd/061421/ \
--quiet \
-A *deduplicated.sorted.bam \ https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/060321-oly/
```
/Volumes/block/wd/061421/1_ATCACG_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam

```{r}
bam.list=list('/Volumes/block/wd/061421/1_ATCACG_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam',
                 '/Volumes/block/wd/061421/2_CGATGT_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam',
                 '/Volumes/block/wd/061421/3_TTAGGC_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam',
                 '/Volumes/block/wd/061421/4_TGACCA_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam',
                 '/Volumes/block/wd/061421/5_ACAGTG_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam',
                 '/Volumes/block/wd/061421/6_GCCAAT_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam',
                 '/Volumes/block/wd/061421/7_CAGATC_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam',
                 '/Volumes/block/wd/061421/8_ACTTGA_L001_R1_001_trimmed_bismark_bt2.deduplicated.sorted.bam'
)
```


```{bash}
ls -lh /Volumes/block/wd/061421/*.sorted.bam
```




```{r}
library(methylKit)
```


```{r, eval = FALSE}
myobj_01 = processBismarkAln(location = sb.list, sample.id = list("1","2","3","4","5","6","7","8"), assembly = "v081", read.context="CpG", mincov=2, treatment = c(0,0,0,0,1,1,1,1))
```

```{r}
save(myobj_01, file = "../analyses/myobj_01")
```



```{r}
load("../analyses/myobj_01")
```



```{r}
filtered.myobj=filterByCoverage(myobj_01,lo.count=10,lo.perc=NULL,
                                      hi.count=1000,hi.perc=NULL)

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

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


PCASamples(meth_filter)

```



```{r}
myDiff=calculateDiffMeth(meth_filter,mc.cores=8)

```




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





```




```{r}
head(filtered.myobj)
```


















```{r}
myDiff_tile=calculateDiffMeth(meth_filter_tile,mc.cores=4)

```


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





```


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





```




```{r}
write.table(myDiff25p, file = "../analyses/myDiff25p.tab", sep = "\t")
```


