---
title: "MethylKit"
author: "Laura H Spencer"
date: "11/17/2022"
output: html_document
knit: (function(input, ...) {
    rmarkdown::render(
      input,
      output_file = paste0(
        xfun::sans_ext(input), '-', Sys.Date(), '.html'
      ),
      envir = globalenv()
    )
  })
---

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

### Load libraries

```{r message=FALSE, warning=FALSE, results=FALSE}
list.of.packages <- c("tidyverse", "reshape2", "here", "methylKit", "ggforce", "matrixStats", "Pstat", "viridis", "plotly", "readr", "GenomicRanges", "vegan", "factoextra", "PerformanceAnalytics", "corrplot", "janitor", "googlesheets4", "ggrepel", "clipr", "plotly") #add new libraries here 

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

# Load all libraries 
lapply(list.of.packages, FUN = function(X) {
  do.call("require", list(X)) 
})

# Load MethylKit
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("methylKit")
require(methylKit)

`%!in%` = Negate(`%in%`)
```

### Obtain session information

```{r}
sessionInfo()
```

### Import bismark alignment .bams into MethylKit. NOTE: This was done on Sedna! 

I ran this code interactively on Sedna 
`srun --pty /bin/bash` # First, I entered an interactive Sedna node
`module load R` # Then, load R
`R` #start R interactively 
Now proceed with the following code 

```{r}
# ## AGAIN- this code was run interactively on Sedna! 
# # Install and load programs
# 
# # Load MethylKit
# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# 
# BiocManager::install("methylKit")
# 
# require(methylKit)
# require(tidyverse)
# 
# # Load sample info dataframe 
# load(file="/home/lspencer/DuMOAR/R/sample.info")
# 
# # Triple check that treatments are in the order of sequence sample number (1-20)
# (sample.info %>% arrange(sample_seq))$sample # List out library prep sample number 
# (sample.info %>% arrange(sample_seq))$sample_seq # List out sequence sample 
# (sample.info %>% arrange(sample_seq))$high_or_low_co2 # List out treatment
# 
# # Read .bam files into Methylkit on Sedna 
# (file.list=list("/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH01_06_S1_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH01_14_S2_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH01_22_S3_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH01_38_S4_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH03_04_S5_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH03_15_S6_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH03_33_S7_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH05_01_S8_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH05_06_S9_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH05_21_S10_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH05_24_S11_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH07_06_S12_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH07_11_S13_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH07_24_S14_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH09_02_S15_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH09_13_S16_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH09_28_S17_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH10_01_S18_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH10_08_S19_R1_bismark_bt2_pe.deduplicated.sorted.bam",
#                "/scratch/lspencer/DuMOAR/aligned/trim-before-concat/CH10_11_S20_R1_bismark_bt2_pe.deduplicated.sorted.bam"))
# 
# #### The following command reads sorted BAM files and calls methylation percentage per base, and creates a methylRaw object for CpG methylation. It also assigns minimum coverage of 2x and the treatments (in this case, the CO2 treatment)
# 
# meth_obj = processBismarkAln(location = file.list, sample.id = list("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18", "19", "20"), 
#                           assembly = "PGA_assembly.fasta", read.context="CpG", mincov=2, 
#                           treatment = as.numeric((sample.info %>% arrange(sample_seq))$high_or_low_co2))
# 
# #### Save the MethylKit object; re-doing the previous step is memory/time intensive, so best to use the saved object moving forward. 
# 
# save(meth_obj, file = "/home/lspencer/DuMOAR/R/meth_obj")
```

```{r}
getwd()
```

#### Read in Bismark summary stats and some library prep stats 

```{r}
# Read in data from GoogleSheet 
sample.info <- read_sheet("https://docs.google.com/spreadsheets/d/1ym0XnYVts98tIUCn0kIaU6VuvqxzV7LoSx9RHwLdiIs/edit?usp=sharing", "DNAisolations_only") %>% 
  clean_names() %>% 
  mutate_at(vars(tank, treatment, high_or_low_co2, developmental_stage, dna_siolation_batch), as.factor) %>%
  
  # Join with bismark summary report 
  right_join(

  # Read in bismark summary report
    read_delim(file="../results/bismark/scaffold-only/bismark_summary_report.txt", delim="\t") %>% mutate(File = str_replace_all(File, '_R1_bismark_bt2_pe.bam', "")) %>%
  clean_names() %>% 
   separate(file, into = c("a", "b", "sample_seq"), remove = F, sep = "_") %>%
   mutate(sample_prep=paste(a, b, sep = "-")) %>%
    mutate(sample_seq=as.numeric(str_replace_all(sample_seq, "S", ""))) %>%
   dplyr::select(file, sample_prep, sample_seq, everything(), -a, -b) %>%
  mutate(perc_totalread_unique=unique_reads_remaining/total_reads, 
         CpGs_Meth=100*methylated_cp_gs/(methylated_cp_gs+unmethylated_cp_gs), 
         CHGs_Meth=100*methylated_chgs/(methylated_chgs+unmethylated_chgs), 
         CHHs_Meth=100*methylated_ch_hs/(methylated_ch_hs+unmethylated_ch_hs)),

  by = c("sample"="sample_prep")) %>%

# add mbd librar prep stats
  left_join(
  
  read_sheet("https://docs.google.com/spreadsheets/d/1ym0XnYVts98tIUCn0kIaU6VuvqxzV7LoSx9RHwLdiIs/edit?usp=sharing", "MBD_calcs") %>% 
  clean_names() %>% 
  dplyr::select(sample, total_recovery_ng, percent_recovery_percent, library_bioanlyzer_mean_fragment_size_bp,
         library_bioanalyzer_molarity_pmol_l, qubit_concentration_ng_u_l), 

  by = "sample")


save(sample.info, file="../data/sample.info")
load(file="../data/sample.info")
```

### Investigate possible factors influencing low %CpGmeth 

#### Correlation bubble plot 

```{r}
cor(sample.info %>% dplyr::select_if(is.numeric)) %>% corrplot(tl.cex=0.65) 
```

#### Look closer at factors that appear to correlate with CpGs_Meth to get correlation statistics  

```{r, warning=FALSE}
#trace("chart.Correlation", edit=T) #to edit function
# In line 17, change the cex value for what you want, i changed it to 2

chart.Correlation(sample.info %>% 
                    dplyr::select(CpGs_Meth, a260_280, mean_mw_highest_peak, aligned_reads, 
                           unaligned_reads, ambiguously_aligned_reads, no_genomic_sequence, 
                           unique_reads_remaining, total_cs, perc_totalread_unique, CHGs_Meth, CHHs_Meth), 
                  histogram=F, pch=19) 
```

Question asked by Zymo (creater of PicoMethyl kit)  - any correlation between the samples with abnormal CpG methylation and samples that had more yield from the MBD enrichment step?

Answer: no not really 

```{r}
chart.Correlation(sample.info %>% 
                    dplyr::select(CpGs_Meth, total_recovery_ng, percent_recovery_percent, 
                           library_bioanlyzer_mean_fragment_size_bp, library_bioanalyzer_molarity_pmol_l, 
                           qubit_concentration_ng_u_l), 
                  histogram=F, pch=19) 
# CpG methylation is not associated with MBD yield


# Use this code to generate 
# sample.info %>% dplyr::select(sample_seq, high_or_low_co2, CpGs_Meth, 
#                        total_recovery_ng, percent_recovery_percent, library_bioanlyzer_mean_fragment_size_bp, 
#                        library_bioanalyzer_molarity_pmol_l, qubit_concentration_ng_u_l) %>%
#   pivot_longer(cols = c("total_recovery_ng", "percent_recovery_percent", "library_bioanlyzer_mean_fragment_size_bp", 
#                        "library_bioanalyzer_molarity_pmol_l", "qubit_concentration_ng_u_l"), names_to = "MBD.stat") %>%
#   mutate(MBD.stat=as.factor(MBD.stat)) %>%
# 
# ggplot(aes(x=value, y=CpGs_Meth)) + 
#   geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
#   scale_color_manual(values=c("red", "blue")) +
#            geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
#   theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
# ggtitle("% CpG Methylation ~ MBD Stats") + 
#   facet_wrap(~MBD.stat, scales = "free")

# Weak correlation with MBD library fragment size, but not convincing
ggplot(sample.info, aes(x=library_bioanlyzer_mean_fragment_size_bp, y=CpGs_Meth)) +
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) +
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) +
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
    annotate("text", x = 500, y = 17, label = "r = 0.42", size=6) + 
  ggtitle("% CpG Methylation ~ library_bioanlyzer_mean_fragment_size_bp") 

# Low CpG methylation is associated with high CHH methylation 
ggplot(sample.info, aes(x=CHHs_Meth, y=CpGs_Meth)) +
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) +
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) +
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
    annotate("text", x = 33, y = 70, label = "r = -0.59", size=6) + 
  ggtitle("% CpG Methylation ~ CHHs_Meth") 
```

Correlations between CpGs_Meth and the following:  
- a260_280: 0.52* - meh     
- mean_mw_highest_peak: 0.51* - not predictive, but low CpG meth samples have lower mean MW.    
- aligned_reads: 0.66** - low CpG meth = fewer aligned reads    
- unaligned_reads: -0.77***   
- ambiguously_aligned_reads: -0.60**   
- unique_reads_remaining: 0.73***  
- total_cs: 0.72***  
- perc_totalread_unique: 0.82***  
- CHGs_Meth: -0.38  
- CHHs_Meth: -0.59**  

```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=a260_280, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 1.5, y = 75, label = "r = 0.52", size=6) + 
ggtitle("% CpG Methylation ~ a260_280")
```

```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=mean_mw_highest_peak, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 10000, y = 100, label = "r = 0.51", size=6) + 
ggtitle("% CpG Methylation ~ Mean MW (highest peak)")
```

```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=aligned_reads, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 6e6, y = 85, label = "r = 0.66", size=6) + 
ggtitle("% CpG Methylation ~ # of aligned reads")
```
```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=unaligned_reads, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 1e7, y = 85, label = "r = -0.77", size=6) + 
ggtitle("% CpG Methylation ~ # unaligned reads")
```

```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=ambiguously_aligned_reads, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 5e6, y = 85, label = "r = -0.60", size=6) + 
ggtitle("% CpG Methylation ~ # ambiguously aligned reads")
```
```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=no_genomic_sequence, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  #annotate("text", x = 285, y = 75, label = "r = -0.93", size=6) + 
ggtitle("% CpG Methylation ~ # unaligned- no genomic sequence")
```


No. of reads not aligning due to no genomic sequence corrlates very well with % methylation, but the actual # of reads is super low across all samples (<300!). What does this mean?! 


```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=unique_reads_remaining, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 3.5e6, y = 95, label = "r = 0.73", size=6) + 
ggtitle("% CpG Methylation ~ # unique aligned reads")
```

```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=total_cs, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 6e7, y = 90, label = "r = 0.72", size=6) + 
ggtitle("% CpG Methylation ~ total # Cs")
```
```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=perc_totalread_unique, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = .15, y = 85, label = "r = 0.82", size=6) + 
ggtitle("% CpG Methylation ~ % total reads uniquely aligned")
```

```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=CHGs_Meth, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 2.25, y = 70, label = "r = -0.38", size=6) + 
ggtitle("% CpG Methylation ~ % CHG Methylation")
```

```{r, warning=FALSE,message=FALSE}
ggplot(sample.info, aes(x=CHHs_Meth, y=CpGs_Meth)) + 
  geom_point(aes(label=sample_seq, color=high_or_low_co2), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  annotate("text", x = 30, y = 75, label = "r = -0.59", size=6) + 
ggtitle("% CpG Methylation ~ % CHH Methylation")
```


### Load methylKit object that was created on Sedna 

```{r}
load("../results/methylkit/meth_obj")
```

#### What does the resulting methylkit object look like? 

```{r}
# Check out data for sample #1 
head(meth_obj[[1]])
```

### Check the basic stats about the methylation data - coverage and percent methylation.

#### First look at % CpG methylation for each sampe 

```{r}
for (i in 1:20) {
getMethylationStats(meth_obj[[i]],plot=T,both.strands=TRUE)
     mtext(paste("Sample", i, sep=" "), side=3, line = -6)
}
```
#### Now look at coverage for each sample 

```{r}
for (i in 1:20) {
 getCoverageStats(meth_obj[[i]],plot=TRUE,both.strands=TRUE) 
   mtext(paste("Sample", i, sep=" "), side=3, line = -6)
}
```

### Generate CpG methylation histograms after filtering for coverage

#### Filter out loci where coverage is less than 5x or greater than 100x. 

```{r}
meth_5x=filterByCoverage(meth_obj,lo.count=5,lo.perc=NULL,
                                      hi.count=100,hi.perc=NULL)
save(meth_5x, file="../results/methylkit/meth_5x")
```

#### Now look at CpG methylation after filtering for 5x coverage 

```{r}
for (i in 1:20) {
getMethylationStats(meth_5x[[i]],plot=T,both.strands=TRUE)
  mtext(paste("Sample", i, sep=" "), side=3, line = -6)
}
```

#### Filter out loci where coverage is less than 10x or greater than 100x. 

```{r}
meth_10x=filterByCoverage(meth_obj,lo.count=10,lo.perc=NULL,
                                      hi.count=100,hi.perc=NULL)
save(meth_10x, file="../results/methylkit/meth_10x")
```

#### Now look at CpG methylation after filtering for 10x coverage 

```{r}
for (i in 1:20) {
getMethylationStats(meth_10x[[i]],plot=T,both.strands=TRUE)
  mtext(paste("Sample", i, sep=" "), side=3, line = -6)
}
```

### Unite methylation data into one object for methylkit functions

#### column bind all samples, and destrand

```{r}
meth_unite=methylKit::unite(meth_obj, destrand=TRUE, min.per.group=1L)
save(meth_unite, file = "../results/methylkit/meth_unite") #save object to file
#load(file = "../results/methylkit/meth_unite") #load object from file if needed
```

### Cluster analysis 

```{r}
clusterSamples(meth_unite, dist="correlation", method="ward", plot=TRUE)
```

```{r}
PCASamples(meth_unite, screeplot = T)
PCASamples(meth_unite)
```

#### Reshape methylation to inspect % methylation at various filtering thresholds 
 
```{r}
# Here - use all methylation data
#load(file = "../results/methylkit/meth_unite") #save object to file

# NOTE: since the sample IDs in the meth_unite() object are sequential from 1:20, the sample IDs in the column names containing data (coverage, numCs, numTs) are correct in this instance. 

meth_unite_reshaped <- meth_unite %>% 
  melt(id=c("chr", "start", "end", "strand"), value.name = "count") %>%
  drop_na(count) %>%
  mutate(stat=gsub('[0-9]+', '', variable)) %>%
  mutate(sample=gsub('[a-zA-Z]+', '', variable)) %>%
  dplyr::select(-variable) %>%
  pivot_wider(names_from = stat, values_from = count) %>%
  mutate(percMeth = 100*(numCs/coverage)) %>%
  mutate(sample=as.numeric(sample))

save(meth_unite_reshaped, file = "../results/methylkit/meth_unite_reshaped")
#load(file = "../results/methylkit/meth_unite_reshaped")  
```

```{r, cache=TRUE, eval=FALSE}
# # Here - use a smaller data set, already filtered for min coverage 5x 
# meth_unite_5x <- methylKit::unite(meth_5x, destrand=TRUE, min.per.group=1L)
# save(meth_unite_5x, file = "../results/methylkit/meth_unite_5x") #save object to file
# load(file = "../results/methylkit/meth_unite_5x") #save object to file
# 
# meth_unite_5x_reshaped <- meth_unite_5x %>%
#   melt(id=c("chr", "start", "end", "strand"), value.name = "count") %>%
#   drop_na(count) %>%
#   mutate(stat=gsub('[0-9]+', '', variable)) %>%
#   mutate(sample=gsub('[a-zA-Z]+', '', variable)) %>%
#   dplyr::select(-variable) %>%
#   pivot_wider(names_from = stat, values_from = count) %>%
#   mutate(percMeth = 100*(numCs/coverage)) %>%
#   mutate(sample=as.numeric(sample))
# 
# save(meth_unite_5x_reshaped, file = "meth_unite_5x_reshaped")
# save(meth_unite_5x_reshaped, file = "../results/methylkit/meth_unite_5x_reshaped")
```

#### No coverage threshold, what is mean % methylated and how many loci are there?

```{r, warning=F}
meth_unite_reshaped %>% 
  group_by(sample) %>%
  summarize(mean = mean(percMeth), nloci = n()) %>%
  left_join(sample.info[c("treatment", "sample_seq", "high_or_low_co2")], by=c("sample"="sample_seq")) %>%
  ggplot(aes(x=nloci, y=mean, label=sample)) + geom_point(aes(colour=high_or_low_co2)) + 
  geom_label_repel(min.segment.length=0.25,aes(colour=high_or_low_co2)) + theme_minimal() +
  ggtitle("% methylation ~ # loci, no coverage minimum") +
  scale_color_manual(values = c("red", "blue")) + 
  geom_smooth(method = "loess", se=T, color="black", linetype="dotted", fill="gray80")
```

#### 5x coverage, what is mean % methylated and how many loci are there?

```{r, warning=F}
meth_unite_reshaped %>% 
  group_by(sample) %>%
  filter(coverage>=5) %>%
  summarize(mean = mean(percMeth), nloci = n()) %>%
  left_join(sample.info[c("treatment", "sample_seq", "high_or_low_co2")], by=c("sample"="sample_seq")) %>%
  ggplot(aes(x=nloci
             , y=mean, label=sample)) + geom_point(aes(colour=high_or_low_co2)) + 
  geom_label_repel(min.segment.length=0.25,aes(colour=high_or_low_co2)) + theme_minimal() +
  ggtitle("% methylation ~ # loci, 5x coverage minimum") +
  scale_color_manual(values = c("red", "blue")) + 
  geom_smooth(method = "loess", se=T, color="black", linetype="dotted", fill="gray80")
```

#### 10x coverage, what is mean % methylated and how many loci are there?

```{r, warning=F}
meth_unite_reshaped %>% 
  group_by(sample) %>%
  filter(coverage>=10) %>%
  summarize(mean = mean(percMeth), nloci = n()) %>% 
  left_join(sample.info[c("treatment", "sample_seq", "high_or_low_co2")], by=c("sample"="sample_seq")) %>%
  ggplot(aes(x=nloci, y=mean, label=sample)) + geom_point(aes(colour=high_or_low_co2)) + 
  geom_label_repel(label.size=0.1, min.segment.length = 0.5, aes(colour=high_or_low_co2)) + theme_minimal() +
  ggtitle("% methylation ~ # loci, 10x coverage minimum") + xlab("# Loci") + ylab("Mean % CpG Methylation") +
  scale_color_manual(values = c("red", "blue")) + 
  geom_smooth(method = "loess", se=T, color="black", linetype="dotted", fill="gray80")
```

#### What is mean % methylated with 15x coverage?

```{r, warning=F}
meth_unite_reshaped %>% 
  group_by(sample) %>%
  filter(coverage>=15) %>%
  summarize(mean = mean(percMeth), nloci = n()) %>%
  left_join(sample.info[c("treatment", "sample_seq", "high_or_low_co2")], by=c("sample"="sample_seq")) %>%
  ggplot(aes(x=nloci, y=mean, label=sample)) + geom_point(aes(colour=high_or_low_co2)) + 
  geom_label_repel(min.segment.length=0.25,aes(colour=high_or_low_co2)) + theme_minimal() +
  ggtitle("% methylation ~ # loci, 15x coverage minimum") +
  scale_color_manual(values = c("red", "blue")) + 
  geom_smooth(method = "loess", se=T, color="black", linetype="dotted", fill="gray80")
```
#### Is % methylation a function of coverage within a sample? 

The following samples have low CpG methylation: 5, 6, 8 (kinda), 9, 13, 14, 16, 19. Check out each sample individually 

Weird samples have lots of loci with 0% methylation 

```{r}
meth_unite_reshaped %>%
  filter(coverage>=5) %>%
  ggplot() + geom_point(aes(x=coverage, y=percMeth), size=0.05) +
  theme_minimal() + 
  facet_wrap(~sample)
```

**Bad samples:** all have mean % methylation <70%
control pco2: 5, 6, 13, 14
high pco2: 9, 16, 19, 8 (kinda)

**Good samples:** all have mean % methylation >= 70% and nloci meeting 10x > 3e5
control pCO2: 1, 2, 3, 4, 7, 12
high pCO2:    10, 11, 15, 17, 18, 20

At 10x coverage, here's the mean % methylation and number of loci per sample 

```{r}
meth_unite_reshaped %>% 
  group_by(sample) %>%
  filter(coverage>=10) %>%
  #filter(sample %in%  c("1", "2", "3", "4", "7", "10", "11", "12", "15", "17", "18", "20")) %>%
  #filter(sample %!in%  c("1", "2", "3", "4", "7", "10", "11", "12", "15", "17", "18", "20")) %>%
  summarize(mean = mean(percMeth), nloci = n()) %>% 
  arrange(sample) #%>% dplyr::select(nloci) %>% write_clip()
```


```{r}
# Summary stats by "good" and "bad" samples  

# Bad samples
sample.info %>% filter(sample_seq %!in%  c("1", "2", "3", "4", "7", "10", "11", "12", "15", "17", "18", "20")) %>%
  summarise(mean.aligned=mean(perc_totalread_unique), sd.aligned=sd(perc_totalread_unique),
            min.aligned=min(perc_totalread_unique), max.aligned=max(perc_totalread_unique),
            mean.reads=mean(total_reads), sd.reads=sd(total_reads),
            mean.CpGs=mean(total_cs), sd.CpGs=sd(total_cs))

# Good samples
sample.info %>% filter(sample_seq %in%  c("1", "2", "3", "4", "7", "10", "11", "12", "15", "17", "18", "20")) %>%
  summarise(mean.aligned=mean(perc_totalread_unique), sd.aligned=sd(perc_totalread_unique),
            min.aligned=min(perc_totalread_unique), max.aligned=max(perc_totalread_unique),
            mean.reads=mean(total_reads), sd.reads=sd(total_reads),
            mean.CpGs=mean(total_cs), sd.CpGs=sd(total_cs))

# sample.info %>% arrange(sample_seq) %>% dplyr::select(sample_seq, sample) %>% 
#   filter(sample_seq %!in%  c("1", "2", "3", "4", "7", "10", "11", "12", "15", "17", "18", "20")) %>% write_clip()
```

### Run analyses without the low-CpG methylation samples 
Good samples: 
control pCO2: 1, 2, 3, 4, 7, 12
high pCO2:    10, 11, 15, 17, 18, 20

#### Generate PCA without weird samples (those with very low meth rates)

```{r}
meth_unite_nolows <- reorganize(methylObj = meth_unite, 
                                sample.ids = c("1", "2", "3", "4", "7", "10", "11", "12", "15", "17", "18", "20"),
                                treatment = as.numeric((sample.info %>% 
                                                          arrange(sample_seq))$high_or_low_co2)[c(1:4, 7, 10, 11, 12, 15, 17, 18, 20)])
PCASamples(meth_unite_nolows)
```

```{r}
(sample.info %>% arrange(sample_seq))[c(1:4, 7, 10, 11, 12, 15, 17, 18, 20),]
```
```{r}
meth_obj_nolow <- reorganize(methylObj = meth_obj, sample.ids = c("1", "2", "3", "4", "7", "10", "11", "12", "15", "17", "18", "20"),
                             treatment = as.numeric((sample.info %>% 
                                                          arrange(sample_seq))$high_or_low_co2)[c(1:4, 7, 10, 11, 12, 15, 17, 18, 20)])

meth_5x=filterByCoverage(meth_obj_nolow,
                         lo.count=5,lo.perc=NULL,
                                      hi.count=100,hi.perc=NULL)

meth_unite_5x=methylKit::unite(meth_5x, destrand=TRUE, min.per.group=5L)
# For some reason the sample id's aren't accurately assigned to each column. 


PCA.5x <- PCASamples(meth_unite_5x, obj.return = TRUE) #Run a PCA analysis on percent methylation for all samples. methylKit uses prcomp to create the PCA matrix
summary(PCA.5x)  
#PC1=11.6% PC2=10.6% 


meth_10x=filterByCoverage(meth_obj_nolow,
                         lo.count=10,lo.perc=NULL,
                                      hi.count=100,hi.perc=NULL)

meth_unite_10x=methylKit::unite(meth_10x, destrand=TRUE)  # by not setting min.per.group manually, all samples must meet threshold for locus to be retained

PCA.10x <- PCASamples(meth_unite_10x, obj.return = TRUE) #Run a PCA analysis on percent methylation for all samples. methylKit uses prcomp to create the PCA matrix
summary(PCA.10x)  
#PC1=12.0% PC2=10.4% 

save(meth_10x, file="../results/methylkit/meth_10x")
save(meth_unite_10x, file="../results/methylkit/meth_unite_10x")
```
Important note!  The column names in the meth_unite_10x object are NOT the sample IDs. 

THE UNITE() FUNCTION ASSIGNS COLUMN NAMES USING THE CONSECUTIVE NUMBERS 1-13, WHILE THE SAMPLE ID'S ARE 1,2,3,4,7,10,11,12,15,17,18,20.  

```{r}
getSampleID(meth_unite_10x) #correct sample IDs
colnames(meth_unite_10x) # columns containing data do NOT have sample iDs, but are rather the n'th item in the vector of sample IDs 
```


PCA analysis and figure 

```{r}
plotCustomization <- (sample.info %>% arrange(sample_seq))[c(1:4, 7, 10, 11, 12, 15, 17, 18, 20),c("high_or_low_co2", "sample_seq", "tank")] %>%
  mutate(color=case_when(high_or_low_co2 == "H" ~ "firebrick3", TRUE ~ "dodgerblue3"),
         symbol=case_when(high_or_low_co2 == "H" ~ 16, TRUE ~ 17))

# Save size: width=650 height=600
par(mar = c(5, 5, 1, 1)) #Specify inner and outer margins
PCA.figure <- ordiplot(PCA.10x, choices = c(1, 2), type = "none", display = "sites", cex = 0.4, xlab = "", ylab = "", xaxt = "n", yaxt = "n") #Use ordiplot to create base biplot. Do not add any points
points(PCA.figure, "sites", col = as.character(plotCustomization$color), pch = plotCustomization$symbol, cex = 2.2)  #use points
#text(PCA.figure$sites, col = as.character(plotCustomization$color), labels=plotCustomization$sample_seq, cex=1.5) #use sample number 

#Add multiple white boxes on top of the default black box to manually change the color
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
#ordiellipse(PCA.filtered, plotCustomization$population, show.groups = "HC", col = "firebrick3") #Add confidence ellipse around the samples in elevated pCO2
#ordiellipse(PCA.filtered, plotCustomization$population, show.groups = "SS", col = "dodgerblue3") #Add confidence ellipse around the samples in ambient pCO2
axis(side =  1, labels = TRUE, col = "grey80", cex.axis = 1.7) #Add x-axis
mtext(side = 1, text = "PC 1 (12.0%)", line = 3, cex = 1.5) #Add x-axis label
axis(side =  2, labels = TRUE, col = "grey80", cex.axis = 1.7) #Add y-axis
mtext(side = 2, text = "PC 2 (10.4%)", line = 3, cex = 1.5) #Add y-axis label
legend("topleft", 
       pch = c(16, 17), 
       legend = c("High pCO2", "Control pCO2"), 
       col = c(c("firebrick3", "dodgerblue3")), 
       cex = 1.3, bty = "n") #Add a legend with information about ambient and elevated samples
```

### Save PC scores and R object for possible integration with other data

NOTE: this PCA was run with _all_ methylated loci data after filtering for 10x across all retained samples. 

```{r}
save(PCA.10x, file="../results/methylkit/PCA.10x")
write_delim(data.frame(PCA.10x$x) %>% tibble::rownames_to_column("sample"), "../results/methylkit/PC-scores-methylation-10x.tab",  delim = '\t', col_names = TRUE)
```

### Quick check to confirm that removed samples should definitely be removed: 

Filter for 10x across ALL samples (this is including the low-CpG methylation samples) to assess clustering. 
Yes, the weird samples still are separated along PC1, which explains the vast majority of variation (68%), from the other samples. Removal of those samples is confirmed. 

```{r}
meth_10x_allsamples=filterByCoverage(meth_obj,
                         lo.count=10,lo.perc=NULL,
                                      hi.count=100,hi.perc=NULL)

meth_unite_10x_allsamples=methylKit::unite(meth_10x_allsamples, destrand=TRUE, min.per.group=10L) #require all 10 samples to have 10x
PCA.10x_allsamples <- PCASamples(meth_unite_10x_allsamples, obj.return = TRUE) #Run a PCA analysis on percent methylation for all samples. methylKit uses prcomp to create the PCA matrix
summary(PCA.10x_allsamples)  
#PC1=68.0% PC2=6.8% 
```

```{r}
plotCustomization_allsamples <- (sample.info %>% arrange(sample_seq))[,c("high_or_low_co2", "sample_seq", "tank")] %>%
  mutate(color=case_when(high_or_low_co2 == "H" ~ "firebrick3", TRUE ~ "dodgerblue3"),
         symbol=case_when(high_or_low_co2 == "H" ~ 16, TRUE ~ 17))

# Save size: width=650 height=600
par(mar = c(5, 5, 1, 1)) #Specify inner and outer margins
PCA.figure_allsamples <- ordiplot(PCA.10x_allsamples, choices = c(1, 2), type = "none", display = "sites", cex = 0.4, xlab = "", ylab = "", xaxt = "n", yaxt = "n") #Use ordiplot to create base biplot. Do not add any points
#points(PCA.figure_allsamples, "sites", col = as.character(plotCustomization_allsamples$color), pch = plotCustomization_allsamples$symbol, cex = 2.5)  #use points
text(PCA.figure_allsamples$sites, col = as.character(plotCustomization_allsamples$color), labels=plotCustomization_allsamples$sample_seq, cex=1.3) #use tank number 

#Add multiple white boxes on top of the default black box to manually change the color
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
box(col = "white")
#ordiellipse(PCA.filtered, plotCustomization$population, show.groups = "HC", col = "firebrick3") #Add confidence ellipse around the samples in elevated pCO2
#ordiellipse(PCA.filtered, plotCustomization$population, show.groups = "SS", col = "dodgerblue3") #Add confidence ellipse around the samples in ambient pCO2
axis(side =  1, labels = TRUE, col = "grey80", cex.axis = 1.7) #Add x-axis
mtext(side = 1, text = "PC 1 (60%)", line = 3, cex = 1.5) #Add x-axis label
axis(side =  2, labels = TRUE, col = "grey80", cex.axis = 1.7) #Add y-axis
mtext(side = 2, text = "PC 2 (9.5%)", line = 3, cex = 1.5) #Add y-axis label
legend("bottomright", 
       pch = c(16, 17), 
       legend = c("High pCO2", "Low pCO2"), 
       col = c(c("firebrick3", "dodgerblue3")), 
       cex = 1.3, bty = "n") #Add a legend with information about ambient and elevated samples
```

```{r, warning=FALSE}
getCorrelation(meth_unite_10x_allsamples,plot=TRUE)
```

### Hierarchical Clustering using filtered methylation data

The function clusters samples using hclust function and various distance metrics derived from percent methylation per base or per region for each sample.

```{r}
clusterSamples(meth_unite_10x, dist="correlation", method="ward", plot=TRUE)
```


```{r}
getData(meth_unite_5x) %>% nrow() #340,394 loci @ 5x coverage
getData(meth_unite_10x) %>% nrow() #146,761 loci @ 10x coverage
```

Percent methylation matrix (rows=region/base, columns=sample) can be extracted from methylBase object by using percMethylation function. This can be useful for downstream analyses. 

Here I create % methylation matrices from the filtered object, and use it to do another cluster analysis 

```{r}
perc.meth=percMethylation(meth_unite_10x, rowids=T)
hist(perc.meth, col="gray50" )

# # Save % methylation df to object and .tab file 
# save(perc.meth, file = "../analyses/methylation-filtered/R-objects/perc.meth") #save object to file 
# #load(file = "../analyses/methylation-filtered/R-objects/perc.meth") #load object if needed
# write.table((as.data.frame(perc.meth) %>% tibble::rownames_to_column("contig")), file = here::here("analyses", "methylation-filtered", "percent-methylation-filtered.tab"), sep = '\t', na = "NA", row.names = FALSE, col.names = TRUE)

# What percentage of loci have ZERO methylation? 
100*table(perc.meth==0)[2]/table(perc.meth==0)[1] 
# 1.8% of loci unmethylated, averaged across all samples 

# Mean % methylation across all samples 
mean(perc.meth, na.rm=TRUE) 
#79.2%

sample.info %>% dplyr::select(high_or_low_co2, sample_seq) %>% arrange(sample_seq)
(perc.meth %>% as.matrix())[,c("1","2","3","4","7","12")] %>% head() # low pCO2 
(perc.meth %>% as.matrix())[,c("10","11","15", "17", "18", "20")] %>% head() # high pCO2 

# Mean % methylation for each treatment 
mean((perc.meth %>% as.matrix())[,c("1","2","3","4","7","12")], na.rm=TRUE) # Low CO2 exposed - 79.5%
mean((perc.meth %>% as.matrix())[,c("10","11","15", "17", "18", "20")], na.rm=TRUE) # High CO2 exposed - 78.9% 

perc.meth.summ <- data.frame(sample=1:12, methylated=1:12, coverage=1:12, mean.perc=1:12)
for (i in 1:12) {
  perc.meth.summ[i,"sample"] <- getSampleID(meth_unite_10x)[i]
  perc.meth.summ[i,"methylated"] <-mean(as.vector(getData(meth_unite_10x)[,grep(c("numCs"), colnames(meth_unite_10x))][i][,1]), na.rm=T)
  perc.meth.summ[i,"coverage"] <- mean(as.vector(getData(meth_unite_10x)[,grep(c("coverage"), colnames(meth_unite_10x))][i][,1]), na.rm=T)
  perc.meth.summ[i,"mean.perc"] <- mean(as.vector(getData(meth_unite_10x)[,grep(c("numCs"), colnames(meth_unite_10x))][i][,1]) / 
                                     as.vector(getData(meth_unite_10x)[,grep(c("coverage"), colnames(meth_unite_10x))][i][,1]), na.rm=T)
}
mean(perc.meth.summ[perc.meth.summ$sample %in% c("1","2","3","4","7","12"),"mean.perc"]) #Loe pCO2 treatment = 79.5%
mean(perc.meth.summ[perc.meth.summ$sample %in% c("10","11","15", "17", "18", "20"),"mean.perc"]) #78.9%

paste("No. of loci in filtered dataset:", nrow(perc.meth)) #146,761 loci 
paste("No. of samples:", ncol(perc.meth)) #12 sampes 

# % CpG Methylation by sample after filtering for 10x coverage (i.e. uniting)
perc.meth.summ[,c("sample", "mean.perc")] %>% mutate(mean.perc=round(100*mean.perc, digits = 1)) %>% View()
```

```{r}
myDiff=calculateDiffMeth(meth_unite_10x)
```

```{r}
# get  deferentially methylated bases using 12% difference and save to files - 12,971 loci 
myDiff12p=getMethylDiff(myDiff,difference=12,qvalue=0.01)
paste("Number diff. methylated loci with 12% difference: ", getData(myDiff12p) %>% nrow())

# get all differentially methylated bases - 1,427 loci
myDiff25p=getMethylDiff(myDiff,difference=25,qvalue=0.01)
paste("Number diff. methylated loci with 25% difference: ", getData(myDiff25p) %>% nrow())

# Extract hypo and hyper methylated loci. Here "type" must think that the low pCO2 is the experimental treatment, since hypo and hyper are switched. 
# get hypo methylated bases - 843 loci
myDiff25p.hyper=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hyper") 
paste("Number loci hypo-methylated in high pCO2 with 25% difference: ", getData(myDiff25p.hyper) %>% nrow())

# get hypo methylated bases - 593 loci 
myDiff25p.hypo=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hypo")
paste("Number loci hyper-methylated in high pCO2 with 25% difference: ", getData(myDiff25p.hypo) %>% nrow())
  
#write.table(myDiff12p, file = "../analyses/DMLs/myDiff12p.tab", sep = "\t", col.names = TRUE)
#write_delim(getData(myDiff12p) %>%  mutate(start = start-1, end = end+1) %>% dplyr::select(chr, start, end, meth.diff),
#            "../analyses/DMLs/dml12.bed",  delim = '\t', col_names = TRUE) 
```

```{r}
write.table(myDiff25p, file = "../results/methylkit/myDiff25p.tab", sep = "\t", col.names = TRUE)
save(myDiff25p, file="../results/methylkit/myDiff25p")
#load(file="../analyses/DMLs/myDiff25p") #load if needed
```

### Take the DMLs to a bed

```{r}
# Write out bedfile for DMLs 
dml25 <- myDiff25p %>% getData() %>% 
  mutate(start = start-1, end = end+1) %>% 
  dplyr::select(chr, start, end, meth.diff)

write_delim(dml25, "../results/methylkit/dml25.bed",  delim = '\t', col_names = TRUE) 
write_delim(dml25, "../results/methylkit/dml25_forIGV.bed",  delim = '\t', col_names = FALSE) #no column names for IGV
save(dml25, file="../results/methylkit/dml25")
#load(file="../results/methylkit/dml25") #load dml25 if needed

# Write out bedfile for all loci fed into DML analysis, to use as background for enrichment 
mydiff.all <-  mutate(getData(myDiff), start = start -1, end = end + 1) %>% dplyr::select(chr, start, end, meth.diff) %>% mutate_if(is.numeric, as.integer)
write_delim(mydiff.all, "../results/methylkit/mydiff-all.bed",  delim = '\t', col_names = TRUE)
save(mydiff.all, file="../results/methylkit/mydiff.all")
#load(file="../results/methylkit/mydiff.all")
```

```{r}
paste("No. of DMLs that had higher methylation in high pCO2 = ", dml25 %>%filter(meth.diff<0)%>%nrow(), " (", 100*round(dml25 %>%filter(meth.diff<0)%>%nrow()/nrow(dml25), digits = 2), "%)", sep="") 
paste("No. of DMLs that had lower methylation in high pCO2= ", dml25 %>%filter(meth.diff>0)%>%nrow(), " (", 100*round(dml25 %>%filter(meth.diff>0)%>%nrow()/nrow(dml25), digits = 2), "%)", sep="") 

# Save DML bed files with hyper- and hypo-methylated loci in high pCO2 treatment 
# meth.diff < 0 = loci where methylation rate is HIGHER in high pCO2 treatment 

write_delim(dml25 %>%filter(meth.diff<0), "../results/methylkit/dml25-hypermeth_in_OA.bed",  delim = '\t', col_names = TRUE) 
write_delim(dml25 %>%filter(meth.diff>0), "../results/methylkit/dml25-hypometh_in_OA.bed",  delim = '\t', col_names = TRUE) 
```

Plot mean % methylation for loci by those that hypo- and hyper-methylated in OA-exposed crab 

```{r}
meth_unite_reshaped %>%
  filter(coverage>=10) %>%
  filter(sample %in% c(1,2,3,4,7,12, #low pCO3
                       8,10,11,15,17,18,20)) %>% #high pCO2
  mutate(locus=paste(chr, start-1, sep=".")) %>%  # Create locus column; subtract 1 from the "start" column to match dml25
  filter(locus %in% (dml25 %>%filter(meth.diff<0) %>% mutate(locus=paste(chr, start, sep=".")))$locus) %>% #keep loci with meth.diff < 0, higher meth in high pCO2
  left_join(sample.info[,c("high_or_low_co2", "sample_seq")], by=c("sample"="sample_seq")) %>%
  group_by(high_or_low_co2, locus) %>%
  summarise(mean = mean(percMeth)) %>%
  ggplot(aes(x=high_or_low_co2, y=mean)) + geom_violin() + geom_jitter() + 
  theme_minimal() + xlab("pCO2 treatment") +
  ggtitle("Loci with higher methylation in OA (n=593)")


#dml25 %>%filter(meth.diff>0) %>% nrow()

meth_unite_reshaped %>%
  filter(coverage>=10) %>%
  filter(sample %in% c(1,2,3,4,7,12, #low pCO3
                     8,10,11,15,17,18,20)) %>% #high pCO2
  mutate(locus=paste(chr, start-1, sep=".")) %>%  # Create locus column; subtract 1 from the "start" column to match dml25
  filter(locus %in% (dml25 %>%filter(meth.diff>0) %>% mutate(locus=paste(chr, start, sep=".")))$locus) %>% #keep loci with meth.diff > 0, lower meth in high pCO2
  left_join(sample.info[,c("high_or_low_co2", "sample_seq")], by=c("sample"="sample_seq")) %>%
  group_by(high_or_low_co2, locus) %>%
  summarise(mean = mean(percMeth)) %>%
  ggplot(aes(x=high_or_low_co2, y=mean)) + geom_violin() + geom_jitter() + 
  theme_minimal() + xlab("pCO2 treatment") +
  ggtitle("Loci with lower methylation in OA (n=834)")
```
###  Subset the filtered methylation count dataframe for only those that are differentially methylated. Save object to file. 

```{r}
dml25_counts <- selectByOverlap(meth_unite_10x,as(myDiff25p,"GRanges"))
class(dml25_counts) #class should be methylBase
perc.methDMLs= percMethylation(dml25_counts, rowids=T) #create a percent methylated object for DMLs

# save count and percent objects to be used in subsequent notebook 
save(dml25_counts, file = "../results/methylkit/dml25_counts")
save(perc.methDMLs, file="../results/methylkit/perc.methDMLs")
```
## permANOVA and dispersion test, % methylation by treatment 

```{r}
perc.meth.d<-vegdist(data.frame(t(perc.meth)),"euclidean", na.rm = TRUE)
meth.perm <- adonis(data.frame(t(perc.meth)) ~ high_or_low_co2, data=plotCustomization,
                          permutations=1000, method='euclidean', na.rm = TRUE) #yes, global population difference

meth.perm$aov.tab # 

hist(meth.perm$f.perms, main="Histogram of pseudo-F values under the null model",
     xlab="pseudo=F values", col="gray", xlim=c(0.5,2.5))
abline(v=1.9743, col="red") # indicates F-value is significant.

source("../references/biostats.R")
pca.eigenval(PCA.10x) #The Proporation of Variance = %variance explained
screeplot(PCA.10x, bstick=TRUE) #that's weird- PC's 1-5 aren't sign. but the higher PCs are. 

# # check sign. of PCs - can't get this to work, though.
# ordi.monte(PCA.10x, ord="pca", dim=2, perm = 1000, scale = F, plot = F)

anova(meth.bdp<-betadisper(perc.meth.d,plotCustomization$high_or_low_co2)) #no sign. difference in dispersion, which is good
boxplot(meth.bdp)

# # perform PCoA just for fun on distance matrix 
# meth.pcoa<-cmdscale(perc.meth.d, eig=TRUE, add=T)
# ordiplot(meth.pcoa, choices = c(1, 2), type="text", display="sites", xlab="PC 1", ylab="PC 2")
```

### Reshape filtered data to create a new column with sample info, and calculate %methylation for each locus/sample 

```{r}
# Get sample ID vector from methylBase object (created from unite() function)
samples <- getSampleID(meth_unite_10x)

meth_unite_10x_reshaped <- melt(data=meth_unite_10x, id=c("chr", "start", "end", "strand"), value.name = "count") %>%
 tidyr::separate(col = "variable" , into = c("variable", "sample_temp"), sep = "(?<=[a-zA-Z])\\s*(?=[0-9])") %>%
  dcast(chr+start+end+strand+sample_temp ~  variable) %>%

  #Correct sample IDs extracted from unite() function 
  mutate(sample=case_when(
    sample_temp== 1 ~ samples[1],
    sample_temp== 2 ~ samples[2],
    sample_temp== 3 ~ samples[3],
    sample_temp== 4 ~ samples[4],
    sample_temp== 5 ~ samples[5],
    sample_temp== 6 ~ samples[6],
    sample_temp== 7 ~ samples[7],
    sample_temp== 8 ~ samples[8],
    sample_temp== 9 ~ samples[9],
    sample_temp== 10 ~ samples[10],
    sample_temp== 11 ~ samples[11],
    sample_temp== 12 ~ samples[12],
    sample_temp== 13 ~ samples[13])) %>%
  dplyr::select(-sample_temp) %>%
  
  mutate(percMeth = 100*(numCs/coverage)) %>% 
  left_join(sample.info[c("sample_seq", "high_or_low_co2")] %>% 
              mutate(sample_seq=as.character(sample_seq)), by=c("sample"="sample_seq")) %>%
  dplyr::rename("treatment"="high_or_low_co2")

save(meth_unite_10x_reshaped, file="../results/methylkit/meth_unite_10x_reshaped")
#load(file="../results/methylkit/meth_unite_10x_reshaped")

# save to tab file if needed
#readr::write_delim(meth_unite_10x_reshaped, "../results/methylkit/meth_unite_10x_reshaped.tab"),  delim = '\t', col_names = FALSE)
```

## Check out overall distribution of filtered methylation data by  treatment 

Any trends? No, similar distribution of % CpG methyation 

```{r}
meth_unite_10x_reshaped %>% 
  ggplot(aes(x= treatment, y=percMeth, fill= treatment)) +
  geom_violin(width=1, alpha=0.8) + theme_minimal() +
  ylab("% CpG Methylation") + theme(axis.title.x = element_blank(), legend.position = "none",
                                    text = element_text(size=16)) +
  scale_x_discrete(labels=c("H" = "High pCO2", "L" = "Control pCO2")) + 
  scale_fill_manual(values=c("firebrick3", "dodgerblue3"))
```

View distributions in density plot 

Looks as though there is slightly more loci with higher methylation rates in the High pCO2 treatment (red distribution is shifted to the right a bit)

```{r}
#ggplotly(
meth_unite_10x_reshaped %>% 
  ggplot(aes(x=percMeth, fill= treatment)) +
  geom_density(alpha=0.5) + 
  scale_fill_manual(values=c("firebrick3", "dodgerblue3"), 
                    labels=c("H" = "High pCO2", "L" = "Control pCO2"), name=NULL) +
  #xlim(90,100) + #use to find location on x axis of second highest peak
  theme_minimal() +
  ylab("Density") + xlab("% CpG Methylation") +
  theme(text = element_text(size=18), legend.position = "top") #)
```

Find location of peaks in each treatment in the above figure. Used this as a reference: http://ianmadd.github.io/pages/PeakDensityDistribution.html

Answer:  
- Largest peak for high pCO2 treatment is at 99.98% 
- Largest peak for low pCO2 treatment (control) is 94.56% 

```{r}
## Find largest peak 

# Find largest peak for High pCO2
density(subset(meth_unite_10x_reshaped,  treatment=="H" & percMeth!="NA")$percMeth) # Y maximum is 4.780e-02

# Find where that maximum lies (this gives the index #)
which.max(density(subset(meth_unite_10x_reshaped,  treatment=="H" & percMeth!="NA")$percMeth)$y) #answer = index #497

# Find the x value where the largest Y value is
density(subset(meth_unite_10x_reshaped,  treatment=="H" & percMeth!="NA")$percMeth)$x[497] #H max (4.895e-02) is @ 99.90%

##--------------

# Find largest peak for Low pCO2, i.e. control
density(subset(meth_unite_10x_reshaped,  treatment=="L" & percMeth!="NA")$percMeth) # Y maximum is 4.663e-02  

which.max(density(subset(meth_unite_10x_reshaped,  treatment=="L" & percMeth!="NA")$percMeth)$y) #answer = 472

density(subset(meth_unite_10x_reshaped,  treatment=="L" & percMeth!="NA")$percMeth)$x[472] #L max (0.0469560  ) is 94.52%

## --------------

## Find 2nd largest peak. NOTE: I looked at the ggplotly density figure above and noted that the 2nd peaks end before ~98%, so I subset the percMeth data to retain only those <98%, then find the peak 

# Find 2nd largest peak for High pCO2
density(subset(meth_unite_10x_reshaped,  treatment=="H" & percMeth!="NA" & percMeth<98)$percMeth) # Y maximum is .04997

# Find where that maximum lies (this gives the index #)
which.max(density(subset(meth_unite_10x_reshaped,  treatment=="H" & percMeth!="NA" & percMeth<98)$percMeth)$y) #answer = index #478

# Find the x value where the largest Y value is
density(subset(meth_unite_10x_reshaped,  treatment=="H" & percMeth!="NA" & percMeth<98)$percMeth)$x[478] #H max (5.281e-02) is @ 94.40%

#---
# Find 2nd largest peak for Low pCO2
density(subset(meth_unite_10x_reshaped,  treatment=="L" & percMeth!="NA" & percMeth>95)$percMeth) # Y maximum is 1.1686

# Find where that maximum lies (this gives the index #)
which.max(density(subset(meth_unite_10x_reshaped,  treatment=="L" & percMeth!="NA" & percMeth>95)$percMeth)$y) #answer = index #475

# Find the x value where the largest Y value is
density(subset(meth_unite_10x_reshaped,  treatment=="L" & percMeth!="NA" & percMeth>95)$percMeth)$x[475] #L max (4.953e-02) is @ 100.00%

# Very similar 2nd peak CpG methylation - High pCO2 = 94.24%, and Low pCO2 = 94.61% (only slightly higher in low pCO2 treatment)

#----- Double check # of loci 
meth_unite_10x_reshaped %>%
  filter( treatment=="H" & percMeth!="NA") %>%
  dplyr::select(chr, start) %>% unique() %>% nrow() #146,761

meth_unite_10x_reshaped %>%
  filter( treatment=="L" & percMeth!="NA") %>%
  dplyr::select(chr, start) %>% unique() %>% nrow() #146,761 (good they are the same)

#----- average % meth by pop

# High pCO2 mean = 78.92%, median = 88.27%
meth_unite_10x_reshaped %>%
  filter( treatment=="H" & percMeth!="NA") %>%
  dplyr::select(percMeth) %>% summary()

# Low (control) pCO2 mean = 79.54%, median = 88.71%
meth_unite_10x_reshaped %>%
  filter( treatment=="L" & percMeth!="NA") %>%
  dplyr::select(percMeth) %>% summary()
```

Takeaway = high and low (control) pCO2 have similar %CpG distributions on average, but the distribution of %meth varies a bit - both treatments' % methylation peak around 94% and 100%, but in the high pCO2 treatment the 100% peak is higher, and the 94% peak is lower. So, high pCO2 has more instances of 100% methylation. 

slightly higher methylation rates in the high pCO2 reared crab. 


### Calculations of % methylation by  treatment 

```{r}
# mean % methylation across samples  
dml25_counts_unique <- unique(getData(dml25_counts)[,c("chr", "start")]) %>%
  mutate(chr_start=paste(chr, start, sep="_"))

#Calculate ratio between mean % methylation in High pCO2 : Low pCO2 
DML.ratios <- meth_unite_10x_reshaped %>% 
  mutate(chr_start=paste(chr, start, sep="_")) %>%
  filter(chr_start %in% dml25_counts_unique$chr_start) %>% 
  group_by( treatment, chr, start) %>%
  summarise(median_percentMeth = median(percMeth, na.rm = TRUE)) %>% #mean_percentMeth = mean(percMeth, na.rm = TRUE), 
  dcast(chr + start ~  treatment) %>% 
  mutate(ratio_H.L = H/L,  #in this ratio column, values <1 = High pCO2 hypomethylated, values >1 = Low pCO2 hypomethylated 
         chr_start=paste(chr, start, sep="_"),
         diff_H.L = H-L)   #this is the difference between High pCO2 & Low pCO2 (H - L) 

nrow(DML.ratios)==nrow(dml25_counts_unique) #confirm same # loci; should=TRUE
save(DML.ratios, file="../results/methylkit/DML.ratios")
#load(file="../analyses/methykit/DML.ratios")

# For all loci, calculate mean and sd percent methylation across samples by  treatment 
DML.calcs <- meth_unite_10x_reshaped %>% 
  group_by( treatment, chr, start) %>% 
  dplyr::summarise(
    mean_percMeth = mean(percMeth, na.rm=TRUE),
    sd_percMeth=sd(percMeth, na.rm=TRUE),
    n())  %>% 
  filter(chr %in% c(DML.ratios$chr) & 
           start %in% c(DML.ratios$start))
```
Interesting! DMLs on average have lower % methylation in the high pCO2 treatment 

```{r}
DML.ratios %>% 
pivot_longer(cols = c(H, L),
   names_to = "treatment",
   values_to = "mean_meth") %>%
  mutate(treatment=as.factor(treatment)) %>%
  ggplot(aes(x=treatment, y=mean_meth, fill=treatment)) +
  geom_violin(alpha=0.2) + geom_boxplot(alpha=0.4) +
  scale_fill_manual(values=c("firebrick3", "dodgerblue3"), 
                    labels=c("H" = "High pCO2", "L" = "Control pCO2"), name=NULL) + 
theme_minimal() + ggtitle("Methylation rates of DMLs") + 
  ylab("Methylation (mean within treatment)") + xlab(NULL)


DML.ratios %>% 
pivot_longer(cols = c(H, L),
   names_to = "treatment",
   values_to = "mean_meth") %>%
  mutate(treatment=as.factor(treatment)) %>%
  group_by(treatment) %>%
  summarise(mean=mean(mean_meth), 
             median=median(mean_meth),)
```

### For ALL loci, calculate difference between mean % methylation in H : L for divergent bar plot 

```{r}
#Calculate ratio between mean % methylation in High pCO2 : Low pCO2 (i.e. control) 
all.loci.ratios <- meth_unite_10x_reshaped %>% 
  mutate(chr_start=paste(chr, start, sep="_")) %>%
  group_by( treatment, chr, start) %>%
  summarise(mean_percentMeth = mean(percMeth, na.rm = TRUE)) %>% #could instead use ... allCs_percent = 100*(sum(numCs)/sum(coverage)), 
  dcast(chr + start ~  treatment) %>% 
  mutate(ratio_H.L = H/L,  #in this ratio column, values <1 = HC hypomethylated, values >1 = SS hypomethylated 
         chr_start=paste(chr, start, sep="_"),
         diff_H.L = H-L)   #this is the difference between H & L (H - L) 

# Any majuor differences in global methylatin among the  treatments? Here are summarys tats of % methylation across all shared loci: 
summary(all.loci.ratios$H) 
summary(all.loci.ratios$L)

# As determined in previous code chunks % meth is slightly higher in the high pCO2 treatment group, while % meth is slightly lower in DMLs in high pCO2 group 

# Do %methylation (all loci) distributions differ among  treatment using all loci? 
#ggplotly(
all.loci.ratios %>% 
pivot_longer(cols = c(H, L),
   names_to = "treatment",
   values_to = "mean_meth") %>%
  ggplot(aes(x=mean_meth, fill= treatment)) +
  geom_density(alpha=0.2) +
    scale_fill_manual(values=c("firebrick3", "dodgerblue3"), 
                    labels=c("H" = "High pCO2", "L" = "Control pCO2"), name=NULL) + 
theme_minimal() #)

# Using DMLs
DML.ratios %>% 
pivot_longer(cols = c(H, L),
   names_to = "treatment",
   values_to = "mean_meth") %>%
  ggplot(aes(x=mean_meth, fill= treatment)) +
  geom_density(alpha=0.2) +
    scale_fill_manual(values=c("firebrick3", "dodgerblue3"), 
                    labels=c("H" = "High pCO2", "L" = "Control pCO2"), name=NULL) + 
theme_minimal() + xlab("Mean methylation (within treatment)") + 
  ggtitle("% Methylation distribution for DMLs by treatment")
```


### Heatmap

```{r}
#load("../analyses/methylation-filtered/R-objects/meth_unite_10x_reshaped")
meth_unite_10x_reshaped %>% dplyr::select(sample, treatment) %>% distinct()

DMLs_heatmap <- 
  meth_unite_10x_reshaped %>% 
    left_join(sample.info[c("sample_seq", "tank")] %>%  mutate(sample_seq=as.character(sample_seq)), by=c("sample"="sample_seq")) %>% 
    arrange(treatment, tank) %>%
  mutate(sample=factor(sample, ordered = T, levels = c(3,4,7,1,2,12,20,17,8,10,18,11,15)), 
         chr_start=paste(chr, start, sep="_")) %>%
  filter(
           start %in% dml25_counts_unique$start & 
           chr %in% DML.ratios$chr)

#ggplotly(
ggplot(DMLs_heatmap, aes(sample, chr_start, fill= percMeth)) + xlab("Sample") + geom_tile(na.rm = T) +
  scale_y_discrete(limits=(DML.ratios[order(DML.ratios$ratio_H.L),]$chr_start)) + 
  #scale_fill_viridis(discrete=FALSE) 
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) + 
  scale_fill_distiller(palette = "YlGnBu", direction = 1) #)
  #scale_fill_gradient(low="gray5", high="white")

sample.info[c("high_or_low_co2", "sample_seq", "tank")]
```
### Alignment to mitochondrial DNA to estimate conversion rate 

```{r}
# Join sample info with bismark summary report 

(sample.info.mito <- sample.info %>% 
   dplyr::select(sample, tank, treatment, high_or_low_co2, sample_seq, developmental_stage, 
                 dna_siolation_batch, perc_totalread_unique, CpGs_Meth, CHGs_Meth, CHHs_Meth) %>% 
  
  right_join(

  # Read in bismark summary report
    read_delim(file="../results/bismark/mito/bismark_summary_report.txt", delim="\t") %>% mutate(File = str_replace_all(File, '_R1_bismark_bt2_pe.bam', "")) %>%
  clean_names() %>% 
   separate(file, into = c("a", "b", "sample_seq"), remove = F, sep = "_") %>%
   mutate(sample_prep=paste(a, b, sep = "-")) %>%
    mutate(sample_seq=as.numeric(str_replace_all(sample_seq, "S", ""))) %>%
   #dplyr::select(file, sample_prep, sample_seq, everything(), -a, -b) %>%
  mutate(perc_totalread_unique_mito=unique_reads_remaining/total_reads, 
         CpGs_Meth_mito=100*methylated_cp_gs/(methylated_cp_gs+unmethylated_cp_gs), 
         CHGs_Meth_mito=100*methylated_chgs/(methylated_chgs+unmethylated_chgs), 
         CHHs_Meth_mito=100*methylated_ch_hs/(methylated_ch_hs+unmethylated_ch_hs)) %>% 
    dplyr::select(sample_prep, perc_totalread_unique_mito, CpGs_Meth_mito, CHGs_Meth_mito, CHHs_Meth_mito),

  by = c("sample"="sample_prep")) %>%
  mutate(included=case_when(
    sample_seq== 1 ~ "yes",
    sample_seq== 2 ~ "yes",
    sample_seq== 3 ~ "yes",     
    sample_seq== 4 ~ "yes",
    sample_seq== 7 ~ "yes",
    sample_seq== 12 ~ "yes",
    sample_seq== 10 ~ "yes",
    sample_seq== 11 ~ "yes",
    sample_seq== 15 ~ "yes",
    sample_seq== 17 ~ "yes",
    sample_seq== 18 ~ "yes",
    sample_seq== 20 ~ "yes",
    TRUE ~ "no")))


# **Good samples:** all have mean % methylation >= 70% and nloci meeting 10x > 3e5
# control pCO2: 1, 2, 3, 4, 7, 12
# high pCO2:    10, 11, 15, 17, 18, 20

ggplot(sample.info.mito, aes(x=CpGs_Meth_mito, y=CpGs_Meth)) + 
  geom_point(aes(shape=high_or_low_co2, color=included), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           #geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  #annotate("text", x = 2.25, y = 70, label = "r = -0.38", size=6) + 
ggtitle("% CpG Methylation full genome ~ mitochondria") + 
  xlab("Mitochindria CpG meth %") + ylab("Full Genome CpG % meth")

ggplot(sample.info.mito, aes(x=perc_totalread_unique_mito, y=CpGs_Meth_mito)) + 
  geom_point(aes(shape=high_or_low_co2, color=included), size=3) + 
  scale_color_manual(values=c("red", "blue")) +
           #geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  #annotate("text", x = 2.25, y = 70, label = "r = -0.38", size=6) + 
ggtitle("mitochondria % CpG methylation ~ alignment rate") + 
  xlab("Mitochindria alignment rate") + ylab("Mitochindria CpG meth %")

ggplot(sample.info.mito, aes(x=included, y=CpGs_Meth_mito, color=included)) + 
  geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.15, aes(shape=high_or_low_co2)) + 
  theme_minimal() + scale_color_manual(values=c("red", "blue")) #+
           #geom_text(aes(label=sample_seq, color=high_or_low_co2), vjust = -0.7) + 
  #theme_minimal() + geom_smooth(method = "lm", se=T, color="black", linetype="dotted", fill="gray80") +
  #annotate("text", x = 2.25, y = 70, label = "r = -0.38", size=6) + 
#ggtitle("% CpG Methylation full genome ~ mitochondria") + 
#  xlab("Mitochindria CpG meth %") + ylab("Full Genome CpG % meth")

summary(sample.info.mito$CpGs_Meth_mito)
```

Create table of alignment rate and methylation rates in full genome and mitochondrial contexts to assess sample integrity 

```{r}
sample.info %>% dplyr::select(sample_seq, sample, high_or_low_co2, treatment, tank, developmental_stage, 
                              dna_siolation_batch, perc_totalread_unique, CpGs_Meth, CHGs_Meth, CHHs_Meth) %>%
  rename("perc_totalread_unique"="perc_totalread_unique_genome", 
         "CpGs_Meth"="CpGs_Meth_genome", "CHGs_Meth"="CHGs_Meth_genome", "CHHs_Meth"="CHHs_Meth_genome") %>% 
  mutate(perc_totalread_unique_genome= perc_totalread_unique_genome*100)  %>% 
  right_join(

  # Read in bismark summary report
    read_delim(file="../results/bismark/mito/bismark_summary_report.txt", delim="\t") %>% mutate(File = str_replace_all(File, '_R1_bismark_bt2_pe.bam', "")) %>%
  clean_names() %>% 
   separate(file, into = c("a", "b", "sample_seq"), remove = F, sep = "_") %>%
   mutate(sample_prep=paste(a, b, sep = "-")) %>%
    mutate(sample_seq=as.numeric(str_replace_all(sample_seq, "S", ""))) %>%
   #dplyr::select(file, sample_prep, sample_seq, everything(), -a, -b) %>%
  mutate(perc_totalread_unique_mito=100*unique_reads_remaining/total_reads, 
         CpGs_Meth_mito=100*methylated_cp_gs/(methylated_cp_gs+unmethylated_cp_gs), 
         CHGs_Meth_mito=100*methylated_chgs/(methylated_chgs+unmethylated_chgs), 
         CHHs_Meth_mito=100*methylated_ch_hs/(methylated_ch_hs+unmethylated_ch_hs)) %>% 
    dplyr::select(sample_prep, perc_totalread_unique_mito, CpGs_Meth_mito, CHGs_Meth_mito, CHHs_Meth_mito),

  by = c("sample"="sample_prep")) %>%
  mutate_at(vars(perc_totalread_unique_genome, CpGs_Meth_genome, CHGs_Meth_genome, CHHs_Meth_genome,
                 CpGs_Meth_mito, CHGs_Meth_mito, CHHs_Meth_mito), funs(round(., 1))) %>%
  mutate(perc_totalread_unique_mito=round(perc_totalread_unique_mito, digits = 2)) %>%
  mutate(included=case_when(
    sample_seq== 1 ~ "yes",
    sample_seq== 2 ~ "yes",
    sample_seq== 3 ~ "yes",     
    sample_seq== 4 ~ "yes",
    sample_seq== 7 ~ "yes",
    sample_seq== 12 ~ "yes",
    sample_seq== 10 ~ "yes",
    sample_seq== 11 ~ "yes",
    sample_seq== 15 ~ "yes",
    sample_seq== 17 ~ "yes",
    sample_seq== 18 ~ "yes",
    sample_seq== 20 ~ "yes",
    TRUE ~ "no")) %>%
  arrange(desc(included), desc(high_or_low_co2), sample_seq) %>% 
  write.csv(file="../results/alignment-meth-summary.csv", row.names=FALSE) 

```


## BONEYARD 

Divergent Bar plot of DMLs 

```{r}
# DML.ratios %>% 
#            mutate(hypermethylated=ifelse(diff_H.L>0, "H", "L")) %>%
#   ggplot(aes(x = chr_start, y = diff_H.L, fill=hypermethylated)) +
#   geom_bar(stat = "identity", width=0.41) +
#     scale_x_discrete(limits=(DML.ratios[order(DML.ratios$diff_H.L),]$chr_start)) + 
#   scale_y_continuous(limits=c(-80,80), breaks=c(-75, -50, -25, 0),
#                        labels=c("-75"="75%", "-50"="50%", "-25"="25%", "0"="0%"),
#                      sec.axis = sec_axis(trans=~.*1, breaks=c(75, 25, 50, 0),
#                        labels=c("-75"="75%", "-50"="50%", "-25"="25%", "0"="0%")))  +
#     theme(axis.text.y=element_blank(), axis.ticks.y = element_blank(), 
#           legend.position = "none", text = element_text(size=14, color="gray30"),
#           panel.background = element_blank(),
#           panel.border = element_rect(colour = "gray30", fill=NA, size=.4)) +
#     #ggtitle("DMLs, showing % methylation difference among  treatments") +
#     scale_fill_manual(values=c("firebrick3", "dodgerblue3")) +
#     #annotate(geom="text", x=-20, y=-65, label="Loci hypermethylated in\nSouth Sound", color="dodgerblue3") +
#     #annotate(geom="text", x=260, y=62, label="Loci hypermethylated in\nHood Canal", color="firebrick3") +
#   ylab("% methylation difference among  treatments") + xlab("Locus") + coord_flip()
# ggsave(filename = "../analyses/DMLs/DMLs-diff-barplot.pdf", width=7, height=5)
```

### Extract correlation matrix (pearson) and distance matrix (manhattan) from the % methylated matrix

FYI I checked that the correlation matrix in `getCorrelation` is derived from the percent methylated matrix. It is! 

```{r}
# cor.pears <- cor(perc.meth, method="pearson") # create pearson correlation matrix. 
# dist.manhat <- dist(t(perc.meth), method = "manhattan", diag = T, upper = F)
# dist.manhat.df <- reshape2::melt(as.matrix(dist.manhat), varnames = c("row", "col"))
# 
# cor.pears.DMLs <- cor(perc.methDMLs, method="pearson") 
# dist.manhat.DMLs <- dist(t(perc.methDMLs), method = "manhattan", diag = T, upper = F) 
# dist.manhat.DMLs.df <- reshape2::melt(as.matrix(dist.manhat.DMLs), varnames = c("row", "col"))
```

### Read in sample number key, merge with distance matrix dataframes

```{r}
#key <- read.csv(here::here("data", "sample-key.csv"))

# dist.manhat.df <- merge(x=merge(x=dist.manhat.df, y=key, by.x="row", by.y="MBD.FILENAME"), y=key, by.x="col", by.y="MBD.FILENAME")
# colnames(dist.manhat.df) <- c("SeqNum.row", "SeqNum.col", "dist.manh", "SampNum.row", "SampNum.col")
# write.csv(dist.manhat.df, file=here::here("analyses", "methylation-filtered", "dist.manhat.csv"), row.names=FALSE) 
# 
# dist.manhat.DMLs.df <- merge(x=merge(x=dist.manhat.DMLs.df, y=key, by.x="row", by.y="MBD.FILENAME"), y=key, by.x="col", by.y="MBD.FILENAME")
# colnames(dist.manhat.DMLs.df) <- c("SeqNum.row", "SeqNum.col", "dist.manh", "SampNum.row", "SampNum.col")
# write.csv(dist.manhat.DMLs.df, file=here::here("analyses", "methylation-filtered", "dist.manhat.DMLs.csv"), row.names=FALSE)
```


### Principal Component Analyses, DMLs only 

```{r, fig.show='hide'}
# load("../analyses/DMLs/R-objects/dml25_counts")
# 
# PCA.DMLs <- PCASamples(dml25_counts, obj.return = TRUE, sd.filter = T) #Run a PCA analysis on percent methylation for all samples. methylKit uses prcomp to create the PCA matrix
# nrow(PCA.DMLs$rotation)
# 
# summary(PCA.DMLs) #Look at summary statistics. The first PC explains 36.3% of variation, the second PC explains 8.5% of variation.  # new analysis/filtering AND >5x replace with NA: PC1=36.5% PC2=9.3%

```

### Save DML PC scores and R object for integration with other data 

```{r}
# save(PCA.DMLs, file="../analyses/DMLs/R-objects/PCA.DMLs")
# write_delim(data.frame(PCA.DMLs$x) %>% tibble::rownames_to_column("sample"), "../analyses/DMLs/PC-scores-DMLs.tab",  delim = '\t', col_names = TRUE)
```

### PCA Figure, DMLs only 

```{r}
# par(mar = c(5, 5, 1, 1)) #Specify inner and outer margins
# PCA.figure <- ordiplot(PCA.DMLs, choices = c(1, 2), type = "none", display = "sites", cex = 0.5, xlab = "", ylab = "", xaxt = "n", yaxt = "n") #Use ordiplot to create base biplot. Do not add any points
# points(PCA.figure, "sites", col = as.character(plotCustomization$color), pch = plotCustomization$symbol, cex = 2.5) #Add each sample. Darker samples are ambient, lighter samples are elevated pCO2
# #Add multiple white boxes on top of the default black box to manually change the color
# box(col = "white")
# box(col = "white")
# box(col = "white")
# box(col = "white")
# box(col = "white")
# box(col = "white")
# box(col = "white")
# box(col = "white")
# box(col = "white")
# box(col = "white")
# #ordiellipse(PCA.DMLs, plotCustomization$ treatment, show.groups = "HC", col = "firebrick3") #Add confidence ellipse around the samples in elevated pCO2
# #ordiellipse(PCA.DMLs, plotCustomization$ treatment, show.groups = "SS", col = "dodgerblue3") #Add confidence ellipse around the samples in ambient pCO2
# axis(side =  1, labels = TRUE, col = "grey80", cex.axis = 1.7) #Add x-axis
# mtext(side = 1, text = "PC 1 (36.5%)", line = 3, cex = 1.5) #Add x-axis label
# axis(side =  2, labels = TRUE, col = "grey80", cex.axis = 1.7) #Add y-axis
# mtext(side = 2, text = "PC 2 (9.3%)", line = 3, cex = 1.5) #Add y-axis label
# legend("topright", 
#        pch = c(16, 17), 
#        legend = c("Hood Canal", "South Sound"), 
#        col = c(c("firebrick3", "dodgerblue3")), 
#        cex = 1.4, bty = "n") #Add a legend with information about ambient and elevated samples
# 
# save(PCA.figure, file="../analyses/DMLs/R-objects/PCA.figure")
```


OTHER POSSIBLE ANALYSIS OPTIONS
- randomize which samples are "Treatment" and run 10 permutations to see now many loci are DMLs then 
- Try to pull SNPs from RNASeq data to identify pop structure 
- Research cool software for integrating meth+rnaseq. AI options for integration? 
