library(methylKit) library(readr) setwd("~/Documents/Geoduck/methratio-out") #file.list2 <- list.files(path = ".", pattern = "*.txt") #for(i in 1:length(file.list2)) { # temp.df <- read_delim(file.list2[i], "\t", escape_double = FALSE, trim_ws = TRUE) # temp.df2 <- temp.df[temp.df$context == 'CG', ] # temp.df2$strand <- gsub("-", "R", temp.df2$strand) # temp.df2$strand <- gsub("\\+", "F", temp.df2$strand) # write.table(temp.df2, paste0("CPG_", file.list2[i]), sep = "\t", col.names = FALSE) #} setwd("~/Documents/Geoduck/methratio-out") file.names <- list("CPG_methratio-out-EPI-103_S27_L005.sam.txt", "CPG_methratio-out-EPI-104_S28_L005.sam.txt", "CPG_methratio-out-EPI-111_S29_L005.sam.txt", "CPG_methratio-out-EPI-113_S30_L005.sam.txt", "CPG_methratio-out-EPI-119_S31_L005.sam.txt", "CPG_methratio-out-EPI-120_S32_L005.sam.txt", "CPG_methratio-out-EPI-127_S33_L005.sam.txt", "CPG_methratio-out-EPI-128_S34_L005.sam.txt", "CPG_methratio-out-EPI-135_S35_L005.sam.txt", "CPG_methratio-out-EPI-136_S36_L005.sam.txt", "CPG_methratio-out-EPI-143_S37_L005.sam.txt", "CPG_methratio-out-EPI-145_S38_L005.sam.txt") sample.names <- list("EPI-103", "EPI-104", "EPI-111", "EPI-113", "EPI-119", "EPI-120", "EPI-127", "EPI-128", "EPI-135", "EPI-136", "EPI-143", "EPI-145") treatment <- c(1, 1, 2, 2, 0, 0, 1, 1, 0, 0, 2, 2) methObj <- methRead(location = file.names, sample.id = sample.names, treatment = treatment, assembly = "10K", mincov = 3, header = TRUE, context = "CpG", resolution = "base", pipeline = list(fraction = TRUE, chr.col = 2, start.col = 3, end.col = 3, coverage.col = 7, strand.col = 4, freqC.col = 6)) filteredObj <- filterByCoverage(methObj, lo.count = 3, low.perc = NULL, hi.count = NULL, hi.perc = 95) for(i in 1:12) { # png(file = paste0("~/Documents/Geoduck/methratio-output/", sample.names[i], "-methStats.png")) getMethylationStats(filteredObj[[i]], plot = T, both.strands = F) dev.copy(jpeg, filename = paste0(sample.names[i], "methstats")) dev.off() } for(i in 1:12) { #png(file = paste0("~/Documents/Geoduck/methratio-output/", sample.names[i], "-covStats.png")) getCoverageStats(filteredObj[[i]], plot = T, both.strands = F) dev.copy(jpeg, filename = paste0(sample.names[i], "covstats")) dev.off() } meth <- unite(filteredObj, destrand = F) #png(file = "~/Documents/Geoduck/methratio-output/Dendrogram.png") CpGDendro <- clusterSamples(meth, dist = "correlation", method = "ward", plot = T) dev.copy(jpeg, filename = paste0(sample.names[i], "dendrogram")) dev.off() #png(file = "~/Documents/Geoduck/methratio-output/pca.png") PCA <- PCASamples(meth, scale = T, center = T) dev.copy(jpeg, filename = paste0(sample.names[i], "pca")) dev.off()