library(methylKit)

library(data.table)

library(GenomicRanges)

file.list <- list('/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/3xCoverage/methylKit/M1.txt','/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/3xCoverage/methylKit/M3.txt')

myobj=read(file.list,sample.id=list("M1","M3"),assembly="v9",treatment=c(1,0))

meth<-unite(myobj)

nrow(meth)

# calculate differential methylation p-values and q-values

myDiff=calculateDiffMeth(meth)

# get differentially methylated regions with 25% difference and qvalue<0.01

myDiff25p=get.methylDiff(myDiff,difference=25, qvalue=0.01)

nrow(myDiff25p)

write.csv(myDiff25p, file="/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/3xCoverage/methylKit/DMRs/myDiff25p_M1M3.csv")