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)