Aggregating methylation ratio by loci

Steven wanted to start looking at the methylome of the geoduck, so the script below is for aggregating all samples percent methylation information into a single (large) file. Since Emu is currently being strangled by PECAN, I’m doing this on my laptop with files I either had from other projects, or grabbed from Owl.

First the Day 10 stuff. Set my working directory and grab all the file names.

setwd("~/Documents/RobertsLab/DMRGraphs")
list.files(path = ".")

 [1] "DiffMethRegions.bed"           "DMR-Box-Plots.pdf"             "EPI-103-percmeth.bedgraph.txt" "EPI-104-percmeth.bedgraph.txt"
 [5] "EPI-111-percmeth.bedgraph.txt" "EPI-113-percmeth.bedgraph.txt" "EPI-119-percmeth.bedgraph.txt" "EPI-120-percmeth.bedgraph.txt"
 [9] "EPI-127-percmeth.bedgraph.txt" "EPI-128-percmeth.bedgraph.txt" "EPI-135-percmeth.bedgraph.txt" "EPI-136-percmeth.bedgraph.txt"
[13] "EPI-143-percmeth.bedgraph.txt" "EPI-145-percmeth.bedgraph.txt" "graphs.nb.html"                "graphs.Rmd"                   
[17] "rsconnect"                    

Start my fence with the first file, read it in and save it to my methCounts data table.

EPI_103_percmeth_bedgraph <- read_delim("~/Documents/RobertsLab/DMRGraphs/EPI-103-percmeth.bedgraph.txt", 
    "\t", escape_double = FALSE, col_names = FALSE, 
    col_types = cols(X4 = col_double()), 
    trim_ws = TRUE)
EPI_103_percmeth_bedgraph <- as.data.table(EPI_103_percmeth_bedgraph)
head(EPI_103_percmeth_bedgraph)

print(EPI_103_percmeth_bedgraph[7,])

Now I iterate over all the other files in the Day10 directory, merging them with the previous version of methCounts. I do two different merges, one that includes all loci, and one that includes only loci shared between samples. Minimum coverage to be included is 5x.

methCounts.all <- methCounts
methCounts.short <- methCounts


for(i in 2:length(sample.files))   {
  
  temp <- read_delim(paste0("~/Documents/RobertsLab/DMRGraphs/", sample.files[i]), 
    "\t", escape_double = FALSE, col_names = FALSE, 
    col_types = cols(X4 = col_double()), 
    trim_ws = TRUE)
  
  temp2 <- as.data.table(paste0(temp$X1, "-", temp$X2))
  temp2$X2 <- temp$X4
  colnames(temp2)[1] <- "loc"
  colnames(temp2)[2] <- substr(sample.files[i], 1, 7)
  
  methCounts.all <- merge(methCounts.all, temp2, by = "loc", all = TRUE)
  methCounts.short <- merge(methCounts.short, temp2, by = "loc")
}

Next I make a column of row means, for average methylation across sample animals.

Next I write the files to disk.

write.table(methCounts.all, file = "Day10-methcounts-all.csv", quote = FALSE, sep = ",", col.names = TRUE, row.names = FALSE)

write.table(methCounts.short, file = "Day10-methcounts-short.csv", quote = FALSE, sep = ",", col.names = TRUE, row.names = FALSE)

I do the same as above, but with Day 135 files.

Then Day 145 files.

Then I merge all the files into a single file and write it to disk.