
# Control Focus

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

## sample metadata



|Sample.ID|OldSample.ID|Treatment|Sex|TreatmentN|Parent.ID|
|---------|------------|---------|---|----------|---------|
|12M      |S12M        |Exposed  |M  |3         |EM05     |
|13M      |S13M        |Control  |M  |1         |CM04     |
|16F      |S16F        |Control  |F  |2         |CF05     |
|19F      |S19F        |Control  |F  |2         |CF08     |
|22F      |S22F        |Exposed  |F  |4         |EF02     |
|23M      |S23M        |Exposed  |M  |3         |EM04     |
|29F      |S29F        |Exposed  |F  |4         |EF07     |
|31M      |S31M        |Exposed  |M  |3         |EM06     |
|35F      |S35F        |Exposed  |F  |4         |EF08     |
|36F      |S36F        |Exposed  |F  |4         |EF05     |
|39F      |S39F        |Control  |F  |2         |CF06     |
|3F       |S3F         |Exposed  |F  |4         |EF06     |
|41F      |S41F        |Exposed  |F  |4         |EF03     |
|44F      |S44F        |Control  |F  |2         |CF03     |
|48M      |S48M        |Exposed  |M  |3         |EM03     |
|50F      |S50F        |Exposed  |F  |4         |EF01     |
|52F      |S52F        |Control  |F  |2         |CF07     |
|53F      |S53F        |Control  |F  |2         |CF02     |
|54F      |S54F        |Control  |F  |2         |CF01     |
|59M      |S59M        |Exposed  |M  |3         |EM01     |
|64M      |S64M        |Control  |M  |1         |CM05     |
|6M       |S6M         |Control  |M  |1         |CM02     |
|76F      |S76F        |Control  |F  |2         |CF04     |
|77F      |S77F        |Exposed  |F  |4         |EF04     |
|7M       |S7M         |Control  |M  |1         |CM01     |
|9M       |S9M         |Exposed  |M  |3         |EM02     |





File list control only..

```{r include=FALSE}
file.list_c=list('../bg_data/13M_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/16F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/19F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/39F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/44F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/52F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/53F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/54F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/64M_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/6M_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/76F_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam',
                '../bg_data/7M_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam'
)
```


```{r eval=FALSE, include=FALSE}
myobj_c = processBismarkAln(location = file.list_c,
  sample.id = list("13M","16F","19F","39F","44F","52F","53F","54F", "64M", "6M", "76F", "7M"),
  assembly = "cv",
  read.context="CpG",
  mincov=2,
  treatment = c(0,1,1,1,1,1,1,1,0,0,1,0))
```


```{r include=FALSE}
save(myobj_c, file = "../analyses/myobj_c")

```

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


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

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

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


PCASamples(meth_filter)

```

```{r include=FALSE}
myDiff_c=calculateDiffMeth(meth_filter,mc.cores=36)

```



```{r include=FALSE}
save(myDiff_c, file = "../analyses/myDiff_c")

```

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

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

```{r}
# get all differentially methylated bases
myDiff_c75p=getMethylDiff(myDiff_c,difference=75,qvalue=0.01)
```



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


```{r}
myDiff_c75p.tab <- read.csv(file = "../analyses/myDiff_c75p.tab", sep = "\t")
```




```{r}
head(myDiff_c75p)
```


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

```{r}
sex_dml <- dplyr::select(myDiff_c75p.tab, chr, start, end, meth.diff) %>% 
  mutate(start = start -1) %>%
  mutate(end = end +1) %>%
  mutate_if(is.numeric, as.integer) %>% 
  mutate(TYPE ="sex_DML") %>% 
  select(chr, start, end, TYPE, meth.diff) 

```


```{r}
write.table(sex_dml, file = "../analyses/sex_dml.bed", sep = "\t", row.names = FALSE,col.names = FALSE, quote = FALSE)

write.table(sex_dml, file = "../igv/sex_dml.bed", sep = "\t", row.names = FALSE,col.names = FALSE, quote = FALSE)
```