---
title: "10- Methylkit"
author:Steven Roberts
date: "`r format(Sys.time(), '%d %B, %Y')`" 
output: 
  github_document:
    toc: true
    toc_depth: 3
    number_sections: true
    html_preview: true
  html_document:
    theme: readable
    highlight: zenburn
    toc: true
    toc_float: true
    number_sections: true
    code_folding: show
    code_download: true
---

```{r setup, include=FALSE, }
library(knitr)
library(tidyverse)
library(kableExtra)
library(DESeq2)
library(pheatmap)
library(RColorBrewer)
library(data.table)
library(DT)
library(Biostrings)
library(methylKit)
knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = FALSE,         # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  fig.width = 6,       # Set plot width in inches
  fig.height = 4,      # Set plot height in inches
  fig.align = "center" # Align plots to the center
)
```

# list of sorted bams

```{r}
analysisFiles <- list(
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/105M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/106M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/107M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/109M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/239M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/240M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/241M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/242M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/269M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/270M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/271M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/272M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/69M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/70M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/71M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/72M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/78M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/79M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/80M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/81M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/92M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/94M_pe.deduplicated.sorted.bam",
  "/home/shared/16TB_HDD_01/sr320/github/project-mytilus-methylation/output/09-meth-quant/95M_pe.deduplicated.sorted.bam"
) #Put all .bam files into a list for analysis.
```


# setting treatment

```{r}
treatment <- list("HC", "HC", "HC", "HC", "AP", "AP", "AP", "AP", "BS", "BS", "BS", "BS", "EB", "EB", "EB", "EB", "SA", "SA", "SA", "SA", "SC", "SC", "SC")
```

```{r}
treatment <- unlist(treatment)
```

# setting sample IDs

```{r}
sample.IDs <- list("105M", "106M", "107M", "109M", "239M", "240M", "241M", "242M", "269M", "270M", "271M", "272M", "69M", "70M", "71M", "72M", "78M", "79M", "80M", "81M", "92M", "94M", "95M") #Create list of sample IDs
```


```{r}
myobj_all = processBismarkAln(location = analysisFiles,
  sample.id = sample.IDs,
  assembly = "Mt",
  read.context="CpG",
  mincov=2,
   treatment = c("HC", "HC", "HC", "HC", "AP", "AP", "AP", "AP", "BS", "BS", "BS", "BS", "EB", "EB", "EB", "EB", "SA", "SA", "SA", "SA", "SC", "SC", "SC"))
```

or 

```{r, eval = FALSE}
myobj_23 = processBismarkAln(location = analysisFiles, sample.id = list("105M", "106M", "107M", "109M", "239M", "240M", "241M", "242M", "269M", "270M", "271M", "272M", "69M", "70M", "71M", "72M", "78M", "79M", "80M", "81M", "92M", "94M", "95M"), assembly = "Mt", read.context="CpG", mincov=2, treatment = c(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1))
```




```{r, eval = FALSE}

library(methylKit)

myobj_23 = processBismarkAln(location = analysisFiles, sample.id = list("105M", "106M", "107M", "109M", "239M", "240M", "241M", "242M", "269M", "270M", "271M", "272M", "69M", "70M", "71M", "72M", "78M", "79M", "80M", "81M", "92M", "94M", "95M"), assembly = "Mt", read.context="CpG", mincov=2, treatment = c("HC", "HC", "HC", "HC", "AP", "AP", "AP", "AP", "BS", "BS", "BS", "BS", "EB", "EB", "EB", "EB", "SA", "SA", "SA", "SA", "SC", "SC", "SC"))
```


```{r}
library(doParallel)

# Set up parallel backend
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

# Run the parallelized code
myobj_all <- foreach(idx = seq_along(analysisFiles), .packages = "methylKit") %dopar% {
  processBismarkAln(
    location = analysisFiles[idx],
    sample.id = sample.IDs[idx],
    assembly = "cv",
    read.context = "CpG",
    mincov = 2,
    treatment = treatment[idx]
  )
}

# Stop the cluster
stopCluster(cl)
```