###############GENERATING A METHYLLIST FOR DIFFERENTIAL METH ANALYSIS#####
#going to try to read in as a list since I can't follow the code to do it separately (***NOTE:! labeling samples as treatment or control here will use averages when generating differential methylation - need to look at this more closely before using these differentially methylated regions for any downstream analysis - it won't matter for correlations though***)
file.list <- list("gonad_methylkit.txt", "mix54formethylkit.txt")
myobj<-read( file.list,pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,
coverage.col=4,strand.col=3,freqC.col=5 ),
sample.id=list("Gonad","Pesticides"),assembly="v9",treatment=c(1,0))
#QC - gives histograms and a few stats
getMethylationStats(myobj[[1]],plot=F,both.strands=F)
getMethylationStats(myobj[[2]],plot=T,both.strands=F)
#unite the samples to create a list of CpG covered by all samples then check it and count the number of rows (or loci) that have data for all samples
meth<-unite(myobj,destrand=FALSE)
head(meth)
nrow(meth)
#check correlation between samples
getCorrelation(meth,plot=T)
#this makes same graph, but saves it higher resolution!
pdf(file="graph_M3_T3D3_gonad.pdf")
getCorrelation(meth,plot=T)
dev.off()
#calculating differentially methylated bases
myDiff<-calculateDiffMeth(meth)
nrow(myDiff)
head(myDiff)
write.csv(myDiff,file="diffmeth")
# # get hypo methylated bases
myDiff25p.hypo <- get.methylDiff(myDiff, difference = 25,
qvalue = 0.01, type = "hypo")
nrow(myDiff25p.hypo)
write.csv(myDiff25p.hypo,file="diffmethylation")