---
title: "Meth_Compare_QiaCompare"
output: html_document
---

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


#library
```{r}
library(methylKit)
library(data.table)
library(dplyr)
setwd("~/GitHub/Meth_Compare")
```


##~~~~~~~~~~~~~~
####FULL DATASET -04/13/20
##~~~~~~~~~~~~~~
#get Mcap dedup the files from gannet
```{bash}
wget -r -l1 --no-parent -A "*merged_CpG_evidence.cov" \
https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/ --directory-prefix="../data/"

```

#get Mcap nodedup the files from gannet
```{bash}
wget -r -l1 --no-parent -A "*merged_CpG_evidence.cov" \
https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/ --directory-prefix="../data/"

```

#get Pact dedup the files from gannet
```{bash}
wget -r -l1 --no-parent -A "*merged_CpG_evidence.cov" \
https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/ --directory-prefix="../data/"

```



#get Pact nodedup the files from gannet
```{bash}
wget -r -l1 --no-parent -A "*merged_CpG_evidence.cov" \
https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/ --directory-prefix="../data/"

```

#make md5 for the cov files
```{bash}
find ~/GitHub/Meth_Compare/data \
-type f -exec md5sum {} + \
> ~/GitHub/Meth_Compare/data/all_emu_cov_files_md5sum.txt
```

#get the MD5 from Gannet
```{bash}
wget  \
https://gannet.fish.washington.edu/metacarcinus/FROGER_meth_compare/20200410/all_031520-TG-bs_files_GANNET_md5sum.txt --directory-prefix="../data/"
```

#compare MD5 of cov files on emu to GANNET - join by MD5 and make sure all cov files on emu match to gannet files
```{r}
emu_md5<-fread("~/GitHub/Meth_Compare/data/all_emu_cov_files_md5sum.txt", header=F)
d5<-fread("~/GitHub/Meth_Compare/data/all_031520-TG-bs_files_GANNET_md5sum.txt", header=F)
join_md5<-left_join(emu_md5,gannet_md5, by = c("V1" = "V1"))
View(join_md5)

```

#get Pact_C1 dedup the files from gannet
```{bash}
wget -r -l1 --no-parent -A "*merged_CpG_evidence.cov" \
https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/032320-Pact-C1/Pact_C1/dedup/ --directory-prefix="../data/"

```

#get Pact_C1 dedup the files from gannet
```{bash}
wget -r -l1 --no-parent -A "*merged_CpG_evidence.cov" \
https://gannet.fish.washington.edu/seashell/bu-mox/scrubbed/032320-Pact-C1/Pact_C1/nodedup/ --directory-prefix="../data/"

```



##~~~~~~~~~~~~~~
####MCap Plots and Analyses ~start here
##~~~~~~~~~~~~~~
#make a list of cov files
```{r}
library(methylKit)
MC_all_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth10_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth11_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth12_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth16_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth17_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth18_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth13_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth14_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth15_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```

#The following command reads coverage files
```{r}
myobj_MC_all_TG=methRead(MC_all_TG.list, sample.id = list("Mcap_W_10","Mcap_W_11","Mcap_W_12","Mcap_M_16","Mcap_M_17","Mcap_M_18","Mcap_R_13","Mcap_R_14","Mcap_R_15"), assembly = "Map_genome", treatment = c(0,0,0,0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.MC_all_TG=filterByCoverage(myobj_MC_all_TG,lo.count=10)

#################UNITE##############################
meth_MC_all_TG<-unite(filtered.myobj.MC_all_TG)
nrow(meth_MC_all_TG)
head(meth_MC_all_TG)
jpeg("../analyses/MethCompare_correlation_MC_all_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_MC_all_TG,plot=TRUE)
dev.off()
```

#make a PCA
```{r}
PCASamples(meth_MC_all_TG)
```



##WGBS v MBD
#make a list of cov files
```{r}
library(methylKit)
MC_WvM_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth10_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth11_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth12_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth16_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth17_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth18_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov'
)

```


#The following command reads coverage files
```{r}
myobj_MC_WvM_TG=methRead(MC_WvM_TG.list, sample.id = list("Mcap_W_10","Mcap_W_11","Mcap_W_12","Mcap_R_16","Mcap_R_17","Mcap_R_18"), assembly = "Map_genome", treatment = c(0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x for all compared
```{r}

#################FILTERING###########################
filtered.myobj.MC_WvM_TG=filterByCoverage(myobj_MC_WvM_TG,lo.count=10)

#################UNITE##############################
meth_MC_WvM_TG<-unite(filtered.myobj.MC_WvM_TG)
nrow(meth_MC_WvM_TG)
head(meth_MC_WvM_TG)
jpeg("../analyses/MethCompare_correlation_MC_WvM_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_MC_WvM_TG,plot=TRUE)
dev.off()
```
#make a PCA
```{r}
PCASamples(meth_MC_WvM_TG)
```

###meth diff
```{r}
#calculate diffs
myDiff_WvM_TG<-calculateDiffMeth(meth_MC_WvM_TG,num.cores = 8,weighted.mean=FALSE,slim=TRUE)
head(myDiff_WvM_TG)
#sig diffs
myDiff_WvM_TG_20p <- getMethylDiff(myDiff_WvM_TG, difference = 50, qvalue = 0.01)
nrow(myDiff_WvM_TG_20p)

```

##WGBS v RRBS
#make a list of cov files
```{r}
library(methylKit)
MC_WvR_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth10_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth11_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth12_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth13_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth14_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth15_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```

#The following command reads coverage files
```{r}
myobj_MC_WvR_TG=methRead(MC_WvR_TG.list, sample.id = list("Mcap_W_10","Mcap_W_11","Mcap_W_12","Mcap_R_13","Mcap_R_14","Mcap_R_15"), assembly = "Map_genome", treatment = c(0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```


#filtering to 10x for all compared
```{r}

#################FILTERING###########################
filtered.myobj.MC_WvR_TG=filterByCoverage(myobj_MC_WvR_TG,lo.count=10)

#################UNITE##############################
meth_MC_WvR_TG<-unite(filtered.myobj.MC_WvR_TG)
nrow(meth_MC_WvR_TG)
head(meth_MC_WvR_TG)
jpeg("../analyses/MethCompare_correlation_MC_WvR_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_MC_WvR_TG,plot=TRUE)
dev.off()
```

#make a PCA
```{r}
PCASamples(meth_MC_WvR_TG)
```

###meth diff
```{r}
#calculate diffs
myDiff_WvR_TG<-calculateDiffMeth(meth_MC_WvR_TG,num.cores = 8,weighted.mean=FALSE,slim=TRUE)
head(myDiff_WvR_TG)
#sig diffs
myDiff_WvR_TG_50p <- getMethylDiff(myDiff_WvR_TG, difference = 50, qvalue = 0.01)
nrow(myDiff_WvR_TG_50p)

```

##MBD v RRBS
#make a list of cov files
```{r}
library(methylKit)
MC_MvR_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth16_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth17_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/dedup/Meth18_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth13_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth14_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Mcap_tg/nodedup/Meth15_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```

#The following command reads coverage files
```{r}
myobj_MC_MvR_TG=methRead(MC_MvR_TG.list, sample.id = list("Mcap_M_16","Mcap_M_17","Mcap_M_18","Mcap_R_13","Mcap_R_14","Mcap_R_15"), assembly = "Map_genome", treatment = c(0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```


#filtering to 10x for all compared
```{r}

#################FILTERING###########################
filtered.myobj.MC_MvR_TG=filterByCoverage(myobj_MC_MvR_TG,lo.count=10)

#################UNITE##############################
meth_MC_MvR_TG<-unite(filtered.myobj.MC_MvR_TG)
nrow(meth_MC_MvR_TG)
head(meth_MC_MvR_TG)
jpeg("../analyses/MethCompare_correlation_MC_MvR_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_MC_MvR_TG,plot=TRUE)
dev.off()
```

#make a PCA
```{r}
PCASamples(meth_MC_MvR_TG)
```

###meth diff
```{r}
#calculate diffs
myDiff_MvR_TG<-calculateDiffMeth(meth_MC_MvR_TG,num.cores = 8,weighted.mean=FALSE,slim=TRUE)
head(myDiff_MvR_TG)

#sig diffs
myDiff_MvR_TG_50p <- getMethylDiff(myDiff_MvR_TG, difference = 50, qvalue = 0.01)
nrow(myDiff_MvR_TG_50p)


```


##~~~~~~~~~~~~~~~~
####FULL DATASET TRIMGALORE _Pact
##~~~~~~~~~~~~~~

#Pact - all methods
#make a list of cov files
```{r}
library(methylKit)
Pact_all_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth1_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth2_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth3_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth7_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth8_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth9_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth4_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth5_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth6_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```


#The following command reads coverage files
```{r}
myobj_Pact_all_TG=methRead(Pact_all_TG.list, sample.id = list("Pact_W_1","Pact_W_2","Pact_W_3","Pact_M_7","Pact_M_8","Pact_M_9","Pact_R_4","Pact_R_5","Pact_R_6"), assembly = "Pact_genome", treatment = c(0,0,0,0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.Pact_all_TG=filterByCoverage(myobj_Pact_all_TG,lo.count=30)

#################UNITE##############################
meth_Pact_all_TG<-unite(filtered.myobj.Pact_all_TG)
nrow(meth_Pact_all_TG)
head(meth_Pact_all_TG)
jpeg("../analyses/MethCompare_correlation_Pact_all_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_Pact_all_TG,plot=TRUE)
dev.off()
```

#Pact - WvM methods
#make a list of cov files
```{r}
library(methylKit)
Pact_WvM_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth1_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth2_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth3_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth7_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth8_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth9_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```


#The following command reads coverage files
```{r}
myobj_Pact_WvM_TG=methRead(Pact_WvM_TG.list, sample.id = list("Pact_W_1","Pact_W_2","Pact_W_3","Pact_M_7","Pact_M_8","Pact_M_9"), assembly = "Pact_genome", treatment = c(0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.Pact_WvM_TG=filterByCoverage(myobj_Pact_WvM_TG,lo.count=10)

#################UNITE##############################
meth_Pact_WvM_TG<-unite(filtered.myobj.Pact_WvM_TG)
nrow(meth_Pact_WvM_TG)
head(meth_Pact_WvM_TG)
jpeg("../analyses/MethCompare_correlation_Pact_WvM_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_Pact_WvM_TG,plot=TRUE)
dev.off()
```

###meth diff
```{r}
#calculate diffs
myDiff_Pact_WvM_TG<-calculateDiffMeth(meth_Pact_WvM_TG,num.cores = 8,weighted.mean=FALSE,slim=TRUE)
head(myDiff_Pact_WvM_TG)

#sig diffs
myDiff_Pact_WvM_TG_50p <- getMethylDiff(myDiff_Pact_WvM_TG, difference = 50, qvalue = 0.01)
nrow(myDiff_Pact_WvM_TG_50p)
tail(myDiff_Pact_WvM_TG)

```


#Pact - WvR
#make a list of cov files
```{r}
library(methylKit)
Pact_WvR_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth1_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth2_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth3_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth4_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth5_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth6_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```

#####START UPDATING SCRIPT HERE -- once confirming Meth4 file#########
#The following command reads coverage files
```{r}
myobj_Pact_WvR_TG=methRead(Pact_WvR_TG.list, sample.id = list("Pact_W_1","Pact_W_2","Pact_W_3","Pact_R_4","Pact_R_5","Pact_R_6"), assembly = "Pact_genome", treatment = c(0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.Pact_WvR_TG=filterByCoverage(myobj_Pact_WvR_TG,lo.count=10)

#################UNITE##############################
meth_Pact_WvR_TG<-unite(filtered.myobj.Pact_WvR_TG)
nrow(meth_Pact_WvR_TG)
head(meth_Pact_WvR_TG)
jpeg("../analyses/MethCompare_correlation_Pact_WvR_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_Pact_WvR_TG,plot=TRUE)
dev.off()
```
###meth diff
```{r}
#calculate diffs
myDiff_Pact_WvR_TG<-calculateDiffMeth(meth_Pact_WvR_TG,num.cores = 8,weighted.mean=FALSE,slim=TRUE)
head(myDiff_Pact_WvR_TG)

#sig diffs
myDiff_Pact_WvR_TG_50p <- getMethylDiff(myDiff_Pact_WvR_TG, difference = 50, qvalue = 0.01)
nrow(myDiff_Pact_WvR_TG_50p)
tail(myDiff_Pact_WvR_TG)

```



#Pact - M v R methods
#make a list of cov files
```{r}
library(methylKit)
Pact_MvR_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth7_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth8_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth9_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth4_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth5_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/nodedup/Meth6_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```


#The following command reads coverage files
```{r}
myobj_Pact_MvR_TG=methRead(Pact_MvR_TG.list, sample.id = list("Pact_M_7","Pact_M_8","Pact_M_9","Pact_R_4","Pact_R_5","Pact_R_6"), assembly = "Pact_genome", treatment = c(0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.Pact_MvR_TG=filterByCoverage(myobj_Pact_MvR_TG,lo.count=10)

#################UNITE##############################
meth_Pact_MvR_TG<-unite(filtered.myobj.Pact_MvR_TG)
nrow(meth_Pact_MvR_TG)
head(meth_Pact_MvR_TG)
jpeg("../analyses/MethCompare_correlation_Pact_MvR_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_Pact_MvR_TG,plot=TRUE)
dev.off()
```

###meth diff
```{r}
#calculate diffs
myDiff_Pact_MvR_TG<-calculateDiffMeth(meth_Pact_MvR_TG,num.cores = 8,weighted.mean=FALSE,slim=TRUE)
head(myDiff_Pact_MvR_TG)

#sig diffs
myDiff_Pact_MvR_TG_50p <- getMethylDiff(myDiff_Pact_MvR_TG, difference = 50, qvalue = 0.01)
nrow(myDiff_Pact_MvR_TG_50p)
tail(myDiff_Pact_MvR_TG)

```

#####Testing stuff --trying within methods to see if individuals are limiting, try this for Pact..
```{r}
library(methylKit)
Pact_W_TG.list=list('../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth1_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth2_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov',
'../data/gannet.fish.washington.edu/seashell/bu-mox/scrubbed/031520-TG-bs/Pact_tg/dedup/Meth3_R1_001_val_1_bismark_bt2_pe..CpG_report.merged_CpG_evidence.cov')

```


#The following command reads coverage files
```{r}
myobj_Pact_W_TG=methRead(Pact_W_TG.list, sample.id = list("Pact_W_1","Pact_W_2","Pact_W_3"), assembly = "Pact_genome", treatment = c(0,0,0), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.Pact_W_TG=filterByCoverage(myobj_Pact_W_TG,lo.count=10)

#################UNITE##############################
meth_Pact_W_TG<-unite(filtered.myobj.Pact_W_TG)
nrow(meth_Pact_W_TG)
head(meth_Pact_all_TG)
```

```{r}
#meth_MC_all_TG -- testing a pretty PCA
library(reshape2)
library(ggplot2)
pcaobj<-PCASamples(meth_MC_all_TG, screeplot = F, adj.lim = c(0.4,0.2),
                    scale = F, center = T, comp = c(1,2), transpose = T, sd.filter = F,
                   sd.threshold = 0, filterByQuantile = F, obj.return = T)

summary(pcaobj)
pcscores<-as.data.frame(pcaobj$x[,1:2])
head(pcscores)
newColNames <- c("species", "method", "sample")
newCols <-colsplit(row.names(pcscores), "_", newColNames)
pcscores_id<-cbind(pcscores, newCols)

jpeg("../analyses/MethCompare_PCA_Mcap_all_041520.jpg", width = 500, height = 500)
ggplot(pcscores_id, aes(PC1, PC2)) +
  geom_point(aes(col=method), size = 4) +
  xlab("PC1 (51%)") + ylab("PC2 (21%)") +
  ggtitle("PCA Mcap %methylation across 3 methods")
dev.off()
```


```{r}
#meth_Pact_all_TG -- testing a pretty PCA
library(reshape2)
library(ggplot2)
pcaobj<-PCASamples(meth_Pact_all_TG, screeplot = F, adj.lim = c(0.4,0.2),
                    scale = F, center = T, comp = c(1,2), transpose = T, sd.filter = F,
                   sd.threshold = 0, filterByQuantile = F, obj.return = T)

summary(pcaobj)
pcscores<-as.data.frame(pcaobj$x[,1:2])
head(pcscores)
newColNames <- c("species", "method", "sample")
newCols <-colsplit(row.names(pcscores), "_", newColNames)
pcscores_id<-cbind(pcscores, newCols)
head(pcscores_id)

jpeg("../analyses/MethCompare_PCA_Pact_all_041520.jpg", width = 500, height = 500)
ggplot(pcscores_id, aes(PC1, PC2)) +
  geom_point(aes(col=method), size = 4) +
  xlab("PC1 (89%)") + ylab("PC2 (6%)") +
  ggtitle("PCA Pact %methylation across 3 methods")
dev.off()
```



```{r}
#meth_MC_all_TG -- %meth table
Mcap_perc<-percMethylation(meth_MC_all_TG,rowids = T)
colMeans(Mcap_perc, na.rm = T, dims = 1)

```

```{r}
#meth_Pact_all_TG -- %meth table
Pact_perc<-percMethylation(meth_Pact_all_TG,rowids = T)
colMeans(Pact_perc, na.rm = T, dims = 1)

```


#Pact - all methods
#make a list of cov files
```{r}
library(methylKit)
Pact_pan_TG.list=list('../analyses/pan-genome-cov/Meth1_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth2_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth3_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth7_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth8_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth9_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth4_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth5_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth6_R1_001_Pact-viaPan-merged.cov')

```


#The following command reads coverage files
```{r}
myobj_Pact_pan_TG=methRead(Pact_pan_TG.list, sample.id = list("Pact_W_1","Pact_W_2","Pact_W_3","Pact_M_7","Pact_M_8","Pact_M_9","Pact_R_4","Pact_R_5","Pact_R_6"), assembly = "Pact_genome", treatment = c(0,0,0,0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.Pact_pan_TG=filterByCoverage(myobj_Pact_pan_TG,lo.count=1)

#################UNITE##############################
meth_Pact_pan_TG<-unite(filtered.myobj.Pact_pan_TG)
nrow(meth_Pact_pan_TG)
head(meth_Pact_pan_TG)
jpeg("../analyses/MethCompare_correlation_Pact_pan_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_Pact_pan_TG,plot=TRUE)
dev.off()
```

#Pact - all methods
#make a list of cov files
```{r}
library(methylKit)
Pact_panMvR_TG.list=list('../analyses/pan-genome-cov/Meth7_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth8_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth9_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth4_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth5_R1_001_Pact-viaPan-merged.cov',
'../analyses/pan-genome-cov/Meth6_R1_001_Pact-viaPan-merged.cov')

```


#The following command reads coverage files
```{r}
myobj_Pact_panMvR_TG=methRead(Pact_panMvR_TG.list, sample.id = list("Pact_M_7","Pact_M_8","Pact_M_9","Pact_R_4","Pact_R_5","Pact_R_6"), assembly = "Pact_genome", treatment = c(0,0,0,1,1,1), context = "CpG", pipeline = "bismarkCoverage")
```

#filtering to 10x  - unite - correlation plot
```{r}

#################FILTERING###########################
filtered.myobj.Pact_panMvR_TG=filterByCoverage(myobj_Pact_panMvR_TG,lo.count=10)

#################UNITE##############################
meth_Pact_panMvR_TG<-unite(filtered.myobj.Pact_panMvR_TG)
nrow(meth_Pact_panMvR_TG)
head(meth_Pact_pan_TG)
jpeg("../analyses/MethCompare_correlation_Pact_pan_TG_041420.jpg", width = 1000, height = 600)
getCorrelation(meth_Pact_panMvR_TG,plot=TRUE)
dev.off()
```