source("/Volumes/web/scaphapoda/Yanouk/Methylkitforpesticides.R") load("/Users/srlab/Downloads/methylKit/data/methylKit.RData") source("/Volumes/web/scaphapoda/Yanouk/Methylkitforpesticides.R") 2*3 2+5989855 load("/Users/srlab/Downloads/methylKit/data/methylKit.RData") install.packages("data.table")# source("http://www.bioconductor.org/biocLite.R")# biocLite("GenomicRanges") source("/Volumes/web/scaphapoda/Yanouk/Methylkitforpesticides.R") install.packages("data.table")# source("http://www.bioconductor.org/biocLite.R")# biocLite("GenomicRanges") go<-read("Mix54methratio_final.txt",sample.id="test",assembly="mm9",header=TRUE, context="CpG", resolution="base",# pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,# coverage.col=4,strand.col=3,freqC.col=5 )# ) go<-read("/Volumes/web/scaphapoda/Yanouk/DecData/Mix54methratio_final.txt",sample.id="test",assembly="mm9",header=TRUE, context="CpG", resolution="base",# pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,# coverage.col=4,strand.col=3,freqC.col=5 )# ) go<-read("TJGR_GonadPE_BS_v9_10x_format.txt",sample.id="test",assembly="mm9",header=TRUE, context="CpG", resolution="base",# pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,# coverage.col=4,strand.col=3,freqC.col=5 ) go<-read("TJGR_GonadPE_BS_v9_10x_format.txt",sample.id="test",assembly="mm9",header=TRUE, context="CpG", resolution="base",# pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,# coverage.col=4,strand.col=3,freqC.col=5 )# ) source("/Volumes/web/scaphapoda/Yanouk/Methylkitforpesticides.R") go<-read("//Volumes/web/scaphapoda/Yanouk/DecData/mix54formethylkit.txt",sample.id="test",assembly="mm9",header=TRUE, context="CpG", resolution="base",# pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,# coverage.col=4,strand.col=3,freqC.col=5 )# ) go<-read("/Volumes/web/scaphapoda/Yanouk/DecData/mix54formethylkit.txt",sample.id="test",assembly="mm9",header=TRUE, context="CpG", resolution="base",# pipeline=list(fraction=TRUE,chr.col=1,start.col=2,end.col=2,# coverage.col=4,strand.col=3,freqC.col=5 )# ) library(methylKit)# # install.packages("data.table")# source("http://www.bioconductor.org/biocLite.R")# biocLite("GenomicRanges")# library(data.table)# library(GenomicRanges) source("/Volumes/web/scaphapoda/Yanouk/Methylkitforpesticides.R") file.list <- list("Control37formethylkit", "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("Control","Pesticides"),assembly="v9",treatment=c(1,0)) file.list <- list("Control37formethylkit.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("Control","Pesticides"),assembly="v9",treatment=c(1,0)) getMethylationStats(myobj[[1]],plot=F,both.strands=F)# getMethylationStats(myobj[[2]],plot=T,both.strands=F) meth<-unite(myobj,destrand=FALSE)# head(meth)# nrow(meth) getCorrelation(meth,plot=T) pdf(file="graph_M3_T3D3_gonad.pdf")# getCorrelation(meth,plot=T)# dev.off() source("/Volumes/web/scaphapoda/Yanouk/Methylkitforpesticides.R") source("/Volumes/web/scaphapoda/Yanouk/Methylkitforpesticides.R") file.list <- list("Control37formethylkit.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("Control","Pesticides"),assembly="v9",treatment=c(1,0)) getMethylationStats(myobj[[1]],plot=F,both.strands=F)# getMethylationStats(myobj[[2]],plot=T,both.strands=F) meth<-unite(myobj,destrand=FALSE)# head(meth)# nrow(meth getCorrelation(meth,plot=T) meth<-unite(myobj,destrand=FALSE)# head(meth)# nrow(meth) getCorrelation(meth,plot=T) pdf(file="graph_M3_T3D3_gonad.pdf")# getCorrelation(meth,plot=T)# dev.off() myDiff<-calculateDiffMeth(meth)# nrow(myDiff)# head(myDiff)# write.csv(myDiff,file="diffmeth") myDiff25p.hypo <- get.methylDiff(myDiff, difference = 25,# qvalue = 0.01, type = "hypo")# nrow(myDiff25p.hypo)# write.csv(myDiff25p.hypo,file="diffmethylation") file.list <- list("gonad_methylkit", "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") file.list <- list("gonad_methylkit", "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)) 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") myDiff25p.hypo <- get.methylDiff(myDiff, difference = 20,# qvalue = 0.01, type = "hypo")# nrow(myDiff25p.hypo)# write.csv(myDiff25p.hypo,file="diffmethylation") get hypo methylated bases# myDiff25p.hypo <- get.methylDiff(myDiff, difference = 10,# qvalue = 0.01, type = "hypo")# nrow(myDiff25p.hypo)# write.csv(myDiff25p.hypo,file="diffmethylation")