methylkit.R

#########
#Script running methylkit on the gold standard/aging data set
#########


library(methylKit)
library(bsseq)
library(BiocParallel)
library(here)
library(data.table)
library(stringr)

mc_cores <- 30
register(MulticoreParam(mc_cores))
set.seed(0)

# Define the directory of output
output_dir <- here::here('results/Methylkit_aging')
if (!dir.exist(output_dir)) dir.create(output_dir, recursive = TRUE)

###############################################################################
## Get Data
###############################################################################

control_fname  <- here::here('data/EGAZ00001016574_90_cpg.txt.gz')
case_fname  <- here::here('data/EGAZ00001016575_90_cpg.txt.gz')

convert_data_format <- function(fname) {
	dt <- fread(fname)
	colnames <- names(dt)
	dt[, total_reads := get(colnames[8]) + get(colnames[9])]
	dt[, strand := '+']
	fwrite(dt, str_replace(fname, '.txt.gz', '.methylkit.gz'), sep="\t")
}

convert_data_format(control_fname)
convert_data_format(case_fname)

file.list <- list(case_fname, control_fname)
file.list <- as.list(str_replace(file.list, '.txt.gz', '.methylkit.gz'))

myobj <- methRead(file.list,
	sample.id=list("case1", "control1"),
	assembly="hg38", treatment=c(1, 0),
	context="CpG", mincov = 10,
	pipeline=list(fraction = TRUE, chr.col = 1, start.col = 2,
		end.col = 3, coverage.col = 12, freqC.col = 7, strand.col = 13)
)

meth <- unite(myobj, destrand=FALSE, min.per.group=1L)

myDiff <- calculateDiffMeth(meth, mc.cores = mc_cores)

# get hyper methylated bases
myDiffHyper <- getMethylDiff(myDiff, qvalue=0.01, type="hyper")

# get hypo methylated bases
myDiffHypo <- getMethylDiff(myDiff, qvalue=0.01, type="hypo")

# get all differentially methylated bases
myDiffAll <- getMethylDiff(myDiff, qvalue=0.01)