library(methylKit)
library(data.table)
library(GenomicRanges)
library(knitr)

file.list <- list('/Volumes/web/Mollusk/bs_larvae_exp/Filtered_Larvae_Files/methylkit/methylkit_out_M1b.txt',
                  '/Volumes/web/Mollusk/bs_larvae_exp/Filtered_Larvae_Files/methylkit/methylkit_out_M3b.txt','/Volumes/web/Mollusk/bs_larvae_exp/Filtered_Larvae_Files/methylkit/methylkit_out_T1D3b.txt',
                  '/Volumes/web/Mollusk/bs_larvae_exp/Filtered_Larvae_Files/methylkit/methylkit_out_T1D5b.txt','/Volumes/web/Mollusk/bs_larvae_exp/Filtered_Larvae_Files/methylkit/methylkit_out_T3D3b.txt',
                  '/Volumes/web/Mollusk/bs_larvae_exp/Filtered_Larvae_Files/methylkit/methylkit_out_T3D5b.txt')

myobj<-read(file.list,pipeline=list(fraction=TRUE,header=T, chr.col=1,start.col=2,end.col=3,strand.col=4, coverage.col=5,freqC.col=6, freqT.col=7),
            sample.id=list("M1","M3","T1D3","T1D5","T3D3","T3D5"),assembly="v9",treatment=c(1,0,0,0,0,0))

meth<-unite(myobj,destrand=FALSE)

head(meth)

nrow(meth)

hc<- clusterSamples(meth, dist="correlation", method="ward", plot=T)

PCA<-PCASamples(meth)
PCAsamples(meth)

getCorrelation(meth,plot=T)

myDiff=calculateDiffMeth(meth)
write.csv(myDiff,file="/Volumes/web/Mollusk/bs_larvae_exp/methylKit/myDiff.csv")
# get hyper methylated bases
myDiff25p.hyper=get.methylDiff(myDiff,difference=25,qvalue=0.01,type="hyper")
# get hypo methylated bases
myDiff25p.hypo=get.methylDiff(myDiff,difference=25,qvalue=0.01,type="hypo")
# get all differentially methylated bases
myDiff25p=get.methylDiff(myDiff,difference=25,qvalue=0.01)