################################################### ### chunk number 1: eval=FALSE ################################################### ## packageDescription("BSgenome") ## packageDescription("gtools") ################################################### ### chunk number 2: eval=FALSE ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("BSgenome") ## biocLite("gtools") ################################################### ### chunk number 3: eval=FALSE ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("MEDIPS") ################################################### ### chunk number 4: ################################################### library("BSgenome") ################################################### ### chunk number 5: ################################################### available.genomes() ################################################### ### chunk number 6: eval=FALSE ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("BSgenome.Hsapiens.UCSC.hg19") ################################################### ### chunk number 7: ################################################### library(MEDIPS) ################################################### ### chunk number 8: ################################################### library(BSgenome.Hsapiens.UCSC.hg19) ################################################### ### chunk number 9: ################################################### file=system.file("extdata", "MeDIP_hESCs_chr22.txt", package="MEDIPS") ################################################### ### chunk number 10: ################################################### CONTROL.SET=MEDIPS.readAlignedSequences(BSgenome="BSgenome.Hsapiens.UCSC.hg19", file=file, numrows=672866) ################################################### ### chunk number 11: ################################################### CONTROL.SET ################################################### ### chunk number 12: ################################################### CONTROL.SET=MEDIPS.genomeVector(data=CONTROL.SET, bin_size=50, extend=400) ################################################### ### chunk number 13: eval=FALSE ################################################### ## MEDIPS.exportWIG(file="output.rpm.control.WIG", data=CONTROL.SET, raw=T, descr="hESCs.rpm") ################################################### ### chunk number 14: ################################################### sr.control=MEDIPS.saturationAnalysis(data=CONTROL.SET, bin_size=50, extend=400, no_iterations=10, no_random_iterations=1) ################################################### ### chunk number 15: ################################################### sr.control ################################################### ### chunk number 16: ################################################### MEDIPS.plotSaturation(sr.control) ################################################### ### chunk number 17: ################################################### CONTROL.SET=MEDIPS.getPositions(data=CONTROL.SET, pattern="CG") ################################################### ### chunk number 18: ################################################### cr.control=MEDIPS.coverageAnalysis(data=CONTROL.SET, extend=400, no_iterations=10) ################################################### ### chunk number 19: ################################################### cr.control ################################################### ### chunk number 20: ################################################### MEDIPS.plotCoverage(cr.control) ################################################### ### chunk number 21: ################################################### er.control=MEDIPS.CpGenrich(data=CONTROL.SET) ################################################### ### chunk number 22: ################################################### er.control ################################################### ### chunk number 23: ################################################### CONTROL.SET=MEDIPS.couplingVector(data=CONTROL.SET, fragmentLength=700, func="count") ################################################### ### chunk number 24: eval=FALSE ################################################### ## MEDIPS.exportWIG(file="PatternDensity.WIG", data=CONTROL.SET, pattern.density=TRUE, descr="Pattern.density") ################################################### ### chunk number 25: ################################################### CONTROL.SET=MEDIPS.calibrationCurve(data=CONTROL.SET) ################################################### ### chunk number 26: ################################################### png("CalibrationPlotCONTROL.png") MEDIPS.plotCalibrationPlot(data=CONTROL.SET, linearFit=T, plot_chr="chr22") dev.off() ################################################### ### chunk number 27: eval=FALSE ################################################### ## MEDIPS.plotCalibrationPlot(data=CONTROL.SET, linearFit=T, plot_chr="chr22") ################################################### ### chunk number 28: ################################################### CONTROL.SET=MEDIPS.normalize(data=CONTROL.SET) ################################################### ### chunk number 29: eval=FALSE ################################################### ## MEDIPS.exportWIG(file="output.rms.control.WIG", data=CONTROL.SET, raw=F, descr="hESCs.rms") ################################################### ### chunk number 30: ################################################### CONTROL.SET ################################################### ### chunk number 31: ################################################### file=system.file("extdata", "hg19.chr22.txt", package="MEDIPS") promoter = MEDIPS.methylProfiling(data1 = CONTROL.SET, ROI_file = file, math = mean, select = 2) ################################################### ### chunk number 32: ################################################### hist(promoter$coupling, breaks=100, main="Promoter CpG densities", xlab="CpG coupling factors") ################################################### ### chunk number 33: ################################################### hist(promoter$rpm_A[promoter$rpm_A!=0], breaks=100, main="RPM signals", xlab="reads/bin") ################################################### ### chunk number 34: ################################################### hist(promoter$ams_A[promoter$ams_A!=0], breaks=100, main="AMS signals", xlab="absolute methylation score (ams)") ################################################### ### chunk number 35: eval=FALSE ################################################### ## frames.frame500.step250=MEDIPS.methylProfiling(data1=CONTROL.SET, frame_size=500, step=250, math=mean, select=2) ################################################### ### chunk number 36: eval=FALSE ################################################### ## write.table(frames.frame500.step250, file="frames.chr22.meth.txt", sep="\t", quote=F, col.names=T, row.names=F) ################################################### ### chunk number 37: eval=FALSE ################################################### ## frames.frame500.step250=read.table(file="frames.chr22.meth.txt", header=T) ################################################### ### chunk number 38: ################################################### file=system.file("extdata", "MeDIP_DE_chr22.txt", package="MEDIPS") TREAT.SET=MEDIPS.readAlignedSequences(BSgenome="BSgenome.Hsapiens.UCSC.hg19", file=file, numrows=863054) ################################################### ### chunk number 39: ################################################### TREAT.SET=MEDIPS.genomeVector(data=TREAT.SET, bin_size=50, extend=400) ################################################### ### chunk number 40: ################################################### TREAT.SET=MEDIPS.getPositions(data=TREAT.SET, pattern="CG") ################################################### ### chunk number 41: ################################################### TREAT.SET=MEDIPS.couplingVector(data=TREAT.SET, fragmentLength=700, func="count") ################################################### ### chunk number 42: ################################################### TREAT.SET=MEDIPS.calibrationCurve(data=TREAT.SET) ################################################### ### chunk number 43: ################################################### TREAT.SET=MEDIPS.normalize(data=TREAT.SET) ################################################### ### chunk number 44: ################################################### file=system.file("extdata", "Input_StemCells_chr22.txt", package="MEDIPS") INPUT.SET=MEDIPS.readAlignedSequences(BSgenome="BSgenome.Hsapiens.UCSC.hg19", file=file, numrows=352756) ################################################### ### chunk number 45: ################################################### INPUT.SET=MEDIPS.genomeVector(data=INPUT.SET, bin_size=50, extend=400) ################################################### ### chunk number 46: eval=FALSE ################################################### ## MEDIPS.exportWIG(file="output.rpm.input.WIG", data=INPUT.SET, raw=T, descr="INPUT.rpm") ################################################### ### chunk number 47: eval=FALSE ################################################### ## sr.input=MEDIPS.saturationAnalysis(data=INPUT.SET, bin_size=50, extend=400, no_iterations=10, no_random_iterations=1) ## MEDIPS.plotSaturation(sr.input) ################################################### ### chunk number 48: eval=FALSE ################################################### ## INPUT.SET=MEDIPS.getPositions(data=INPUT.SET, pattern="CG") ################################################### ### chunk number 49: eval=FALSE ################################################### ## cr.input=MEDIPS.coverageAnalysis(data=INPUT.SET, extend=400, no_iterations=10) ## MEDIPS.plotCoverage(cr.input) ################################################### ### chunk number 50: eval=FALSE ################################################### ## er.input=MEDIPS.CpGenrich(data=INPUT.SET) ################################################### ### chunk number 51: eval=FALSE ################################################### ## INPUT.SET=MEDIPS.couplingVector(data=INPUT.SET, fragmentLength=700, func="count") ################################################### ### chunk number 52: eval=FALSE ################################################### ## INPUT.SET=MEDIPS.calibrationCurve(data=INPUT.SET) ################################################### ### chunk number 53: eval=FALSE ################################################### ## png("CalibrationPlotINPUT.png") ## MEDIPS.plotCalibrationPlot(data=INPUT.SET, linearFit=F, plot_chr="chr22") ## dev.off() ################################################### ### chunk number 54: ################################################### diff.meth=MEDIPS.methylProfiling(data1=CONTROL.SET, data2=TREAT.SET, input=INPUT.SET, frame_size=500, select=2) ################################################### ### chunk number 55: eval=FALSE ################################################### ## write.table(diff.meth, file="diff.meth.txt", sep="\t", quote=F, col.names=T, row.names=F) ################################################### ### chunk number 56: eval=FALSE ################################################### ## diff.meth=read.table(file="diff.meth.txt", header=T) ################################################### ### chunk number 57: ################################################### diff.meth.control=MEDIPS.selectSignificants(frames=diff.meth, control=T, input=T, up=2, p.value=0.0001, quant=0.99) ################################################### ### chunk number 58: eval=FALSE ################################################### ## diff.meth.control.merged=MEDIPS.mergeFrames(frames=diff.meth.control) ################################################### ### chunk number 59: eval=FALSE ################################################### ## write.table(diff.meth.control, file = "DiffMethyl.Up.hESCs.bed", sep = "\t", quote = F, row.names = F, col.names = F) ################################################### ### chunk number 60: ################################################### file=system.file("extdata", "hg19.chr22.txt", package="MEDIPS") diff.meth.control.annotated=MEDIPS.annotate(diff.meth.control, anno=file) ################################################### ### chunk number 61: ################################################### length(unique(diff.meth.control.annotated[,4])) ################################################### ### chunk number 62: ################################################### diff.meth.treat=MEDIPS.selectSignificants(frames=diff.meth, control=F, input=T, up=2, down=0.5, p.value=0.0001, quant=0.99) ################################################### ### chunk number 63: eval=FALSE ################################################### ## write.table(diff.meth.treat, file = "DiffMethyl.Up.DE.bed", sep = "\t", quote = F, row.names = F, col.names = F) ################################################### ### chunk number 64: ################################################### file=system.file("extdata", "hg19.chr22.txt", package="MEDIPS") diff.meth.treat.annotated=MEDIPS.annotate(diff.meth.treat, anno=file) ################################################### ### chunk number 65: ################################################### length(unique(diff.meth.treat.annotated[,4]))