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

file.list <- list('/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/methylKit_GilesCode/M1_methylkit.txt','/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/methylKit_GilesCode/M3_methylkit.txt',
                  '/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/methylKit_GilesCode/T1D3_methylkit.txt', '/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/methylKit_GilesCode/T1D5_methylkit.txt',
                  '/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/methylKit_GilesCode/T3D3_methylkit.txt', '/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/methylKit_GilesCode/T3D5_methylkit.txt')

myobj=read(file.list,sample.id=list("M1","M3","T1D3","T1D5","T3D3","T3D5"),assembly="v9",treatment=c(1,1,1,1,1,1))

meth<-unite(myobj)

head(meth)

nrow(meth)

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

PCA<-PCASamples(meth)

getCorrelation(meth,plot=T)

PM<-percMethylation(meth)

write.csv(PM, file="/Volumes/web/Mollusk/bs_larvae_exp/AdapterTrimmed_FilteredFiles/methylKit_GilesCode/MethylationValuesMatrix.csv")