Generating Day 10 genomic feature files.

I needed to make some bedgraphs for Steven, and since the code I’d used previously was buried in a larger section, I decided to break it out into it’s own little script.

Housekeeping stuff, loading needed packages, setting working directory, and getting all of the CG output files from Bismark into a vector to iterate over.

library(readr)
library(data.table)

setwd("~/Documents/Genewiz_hdd/Day10/CGoutput/")

files <- list.files(".", pattern = "*.output")

print(files)

Loops over the file names reading in the data, removing all rows that contain an NA, saving location data (I subtract 1 because IGV is apparently 0 indexed, so when you don’t, everything is shifted one to the right), combining the two strand data for methylation and total coverage, and then setting a 5x coverage requirement to be included.

I then save the file as a bedgraph, and use sed -i to append the minimum header required to count as a bedgraph type=bedGraph.

for(i in 1:length(files))   {
  
  temp <- read_delim(paste0("~/Documents/Genewiz_hdd/Day10/CGoutput/", files[i]), 
    "\t", escape_double = FALSE, trim_ws = TRUE, na = ".")
  
  temp <- na.omit(temp)
  
  temp.2 <- as.data.table(cbind(temp$`#CHROM`))
  temp.2$sloc <-temp$POS - 1
  temp.2$eloc <- temp$POS - 1
  
  temp.2$meth <- temp$`Watson METH` + temp$`Crick METH`
  temp.2$totcov <- temp$`Watson COVERAGE` + temp$`Crick COVERAGE`  
  temp.3 <- temp.2[temp.2$totcov > 5,]    

  temp.final <- as.data.table(cbind(temp.3$V1, temp.3$sloc, temp.3$eloc, (temp.3$meth / temp.3$totcov)))  
  
  
  filename <- paste0(substr(files[i], 1, 7), "_percmeth.bedgraph")
  
  write.table(temp.final, file = filename, quote = FALSE, sep = "\t", col.names = FALSE, row.names = FALSE)
  system(paste0("sed -i 'itype=bedGraph' ", filename))
}