---
title: "Re-running ZR1394 samples in BSMap"
output: html_notebook
---





```{r}

setwd("~/Documents/Bismark-Mapping-Eff/Galore")
file.names.list <- list.files("~/Documents/Bismark-Mapping-Eff/Galore/", pattern = "*.fq")


for(i in 1:length(file.names.list))   {
  
  system(paste0("/home/shared/bsmap/BSMAP-master/bsmap -a ~/Documents/Bismark-Mapping-Eff/Galore/", file.names.list[i], " -d ~/Documents/Bismark-Mapping-Eff/Bowtie1/Ostrea_lurida-Scaff-10k.fa -o ~/Documents/Bismark-Mapping-Eff/Galore/", file.names.list[i], ".sam"))
  
}


```

Now to run everything through Methratio, the BSMap methylation extractor. I tried running a BSMap generated file through bismark methylation extractor, but it generated an error due to unexpected characters in the file. 

```{r}

meth.files.list <- list.files(path = "~/Documents/Bismark-Mapping-Eff/Galore/", pattern = "*.sam")

for(i in 1:length(meth.files.list))   {
  
  system(paste0("python /home/shared/bsmap/BSMAP-master/methratio.py -u -z -d ~/Documents/Bismark-Mapping-Eff/Bowtie2/Ostrea_lurida-Scaff-10k.fa -o ~/Documents/Bismark-Mapping-Eff/Galore/methratio-out/methratio-out-", meth.files.list[i], ".txt ~/Documents/Bismark-Mapping-Eff/Galore/", meth.files.list[i]))
}

```

Now we move on to Methylkit, a native R package.

```{r}
library(readr)
setwd("~/Documents/Bismark-Mapping-Eff/Galore/methratio-out")

file.list2 <- list.files(path = ".", pattern = "*.txt")

for(i in 1:length(file.list2))   {
  temp.df <- read_delim(file.list2[i], "\t", escape_double = FALSE, trim_ws = TRUE)
  temp.df2 <- temp.df[temp.df$context == 'CG', ]
  temp.df2$strand <- gsub("-", "R", temp.df2$strand)
  temp.df2$strand <- gsub("\\+", "F", temp.df2$strand)
  write.table(temp.df2, paste0("CPG_", file.list2[i]), sep = "\t", col.names = FALSE)
  
}

```

```{r}
file.list3 <- list("CPG_methratio-out-zr1394_1_trimmed.fq.sam.txt", "CPG_methratio-out-zr1394_2_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_3_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_4_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_5_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_6_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_7_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_8_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_9_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_10_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_11_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_12_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_13_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_14_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_15_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_16_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_17_trimmed.fq.sam.txt",
"CPG_methratio-out-zr1394_18_trimmed.fq.sam.txt")

sample.id <- list("hc1_2B",
"hc1_4B",
"hc2_15B",
"hc2_17",
"hc3_1",
"hc3_5",
"hc3_7",
"hc3_10",
"hc3_11",
"ss2_9B",
"ss2_14B",
"ss2_18B",
"ss3_3B",
"ss3_14B",
"ss3_15B",
"ss3_16B",
"ss3_20",
"ss3_18")

treatment <- c(rep(0, 9), rep(1,9))
setwd("~/Documents/Bismark-Mapping-Eff/Galore/methratio-out")
testObj <- methRead(file.list3, sample.id = sample.id, treatment = treatment, assembly = "10K", mincov = 3, header = TRUE, context = "CpG", resolution = "base", pipeline = list(fraction = TRUE, chr.col = 2, start.col = 3, end.col = 3, coverage.col = 7, strand.col = 4, freqC.col = 6))

filteredObj <- filterByCoverage(testObj, lo.count = 3, 
    low.perc = NULL, hi.count = NULL, hi.perc = 95)  
```

```{r}

for(i in 1:18)   {
  getMethylationStats(filteredObj[[i]], plot = T, both.strands = F)
}

```

```{r}

for(i in 1:18)   {
  getCoverageStats(filteredObj[[i]], plot = T, both.strands = F)
}

```

```{r}

meth <- unite(filteredObj, destrand = F)

getCorrelation(meth, plot = F)

CpGDendro <- clusterSamples(meth, dist = "correlation", method = "ward", plot = T)

PCA <- PCASamples(meth, scale = T, center = T)
```