### R code from vignette source 'vignettes/MEDIPS/inst/doc/MEDIPS.Rnw' ################################################### ### code chunk number 1: MEDIPS.Rnw:61-63 (eval = FALSE) ################################################### ## packageDescription("BSgenome") ## packageDescription("gtools") ################################################### ### code chunk number 2: MEDIPS.Rnw:69-72 (eval = FALSE) ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("BSgenome") ## biocLite("gtools") ################################################### ### code chunk number 3: MEDIPS.Rnw:82-84 (eval = FALSE) ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("MEDIPS") ################################################### ### code chunk number 4: MEDIPS.Rnw:124-125 ################################################### library("BSgenome") ################################################### ### code chunk number 5: MEDIPS.Rnw:128-129 ################################################### available.genomes() ################################################### ### code chunk number 6: MEDIPS.Rnw:135-137 (eval = FALSE) ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("BSgenome.Hsapiens.UCSC.hg19") ################################################### ### code chunk number 7: MEDIPS.Rnw:151-152 ################################################### library(MEDIPS) ################################################### ### code chunk number 8: MEDIPS.Rnw:158-159 ################################################### library(BSgenome.Hsapiens.UCSC.hg19) ################################################### ### code chunk number 9: MEDIPS.Rnw:171-172 ################################################### file=system.file("extdata", "MeDIP_hESCs_chr22.txt", package="MEDIPS") ################################################### ### code chunk number 10: MEDIPS.Rnw:178-179 ################################################### CONTROL.SET=MEDIPS.readAlignedSequences(BSgenome="BSgenome.Hsapiens.UCSC.hg19", file=file, numrows=672866) ################################################### ### code chunk number 11: MEDIPS.Rnw:187-188 ################################################### CONTROL.SET ################################################### ### code chunk number 12: MEDIPS.Rnw:217-218 ################################################### CONTROL.SET=MEDIPS.genomeVector(data=CONTROL.SET, bin_size=50, extend=400) ################################################### ### code chunk number 13: MEDIPS.Rnw:236-237 (eval = FALSE) ################################################### ## MEDIPS.exportWIG(file="output.rpm.control.WIG", data=CONTROL.SET, raw=T, descr="hESCs.rpm") ################################################### ### code chunk number 14: MEDIPS.Rnw:256-257 ################################################### sr.control=MEDIPS.saturationAnalysis(data=CONTROL.SET, bin_size=50, extend=400, no_iterations=10, no_random_iterations=1) ################################################### ### code chunk number 15: MEDIPS.Rnw:278-279 ################################################### sr.control ################################################### ### code chunk number 16: MEDIPS.Rnw:289-290 ################################################### MEDIPS.plotSaturation(sr.control) ################################################### ### code chunk number 17: MEDIPS.Rnw:326-327 ################################################### CONTROL.SET=MEDIPS.getPositions(data=CONTROL.SET, pattern="CG") ################################################### ### code chunk number 18: MEDIPS.Rnw:348-349 ################################################### cr.control=MEDIPS.coverageAnalysis(data=CONTROL.SET, extend=400, no_iterations=10) ################################################### ### code chunk number 19: MEDIPS.Rnw:360-361 ################################################### cr.control ################################################### ### code chunk number 20: MEDIPS.Rnw:372-373 ################################################### MEDIPS.plotCoverage(cr.control) ################################################### ### code chunk number 21: MEDIPS.Rnw:405-406 ################################################### er.control=MEDIPS.CpGenrich(data=CONTROL.SET) ################################################### ### code chunk number 22: MEDIPS.Rnw:412-413 ################################################### er.control ################################################### ### code chunk number 23: MEDIPS.Rnw:441-442 ################################################### CONTROL.SET=MEDIPS.couplingVector(data=CONTROL.SET, fragmentLength=700, func="count") ################################################### ### code chunk number 24: MEDIPS.Rnw:482-483 (eval = FALSE) ################################################### ## MEDIPS.exportWIG(file="PatternDensity.WIG", data=CONTROL.SET, pattern.density=TRUE, descr="Pattern.density") ################################################### ### code chunk number 25: MEDIPS.Rnw:499-500 ################################################### CONTROL.SET=MEDIPS.calibrationCurve(data=CONTROL.SET) ################################################### ### code chunk number 26: MEDIPS.Rnw:512-515 ################################################### png("CalibrationPlotCONTROL.png") MEDIPS.plotCalibrationPlot(data=CONTROL.SET, linearFit=T, plot_chr="chr22") dev.off() ################################################### ### code chunk number 27: MEDIPS.Rnw:518-519 (eval = FALSE) ################################################### ## MEDIPS.plotCalibrationPlot(data=CONTROL.SET, linearFit=T, plot_chr="chr22") ################################################### ### code chunk number 28: MEDIPS.Rnw:645-646 ################################################### CONTROL.SET=MEDIPS.normalize(data=CONTROL.SET) ################################################### ### code chunk number 29: MEDIPS.Rnw:667-668 (eval = FALSE) ################################################### ## MEDIPS.exportWIG(file="output.rms.control.WIG", data=CONTROL.SET, raw=F, descr="hESCs.rms") ################################################### ### code chunk number 30: MEDIPS.Rnw:674-675 ################################################### CONTROL.SET ################################################### ### code chunk number 31: MEDIPS.Rnw:711-713 ################################################### file=system.file("extdata", "hg19.chr22.txt", package="MEDIPS") promoter = MEDIPS.methylProfiling(data1 = CONTROL.SET, ROI_file = file, math = mean, select = 2) ################################################### ### code chunk number 32: MEDIPS.Rnw:835-836 ################################################### hist(promoter$coupling, breaks=100, main="Promoter CpG densities", xlab="CpG coupling factors") ################################################### ### code chunk number 33: MEDIPS.Rnw:842-843 ################################################### hist(promoter$rpm_A[promoter$rpm_A!=0], breaks=100, main="RPM signals", xlab="reads/bin") ################################################### ### code chunk number 34: MEDIPS.Rnw:850-851 ################################################### hist(promoter$ams_A[promoter$ams_A!=0], breaks=100, main="AMS signals", xlab="absolute methylation score (ams)") ################################################### ### code chunk number 35: MEDIPS.Rnw:865-866 (eval = FALSE) ################################################### ## frames.frame500.step250=MEDIPS.methylProfiling(data1=CONTROL.SET, frame_size=500, step=250, math=mean, select=2) ################################################### ### code chunk number 36: MEDIPS.Rnw:873-874 (eval = FALSE) ################################################### ## write.table(frames.frame500.step250, file="frames.chr22.meth.txt", sep="\t", quote=F, col.names=T, row.names=F) ################################################### ### code chunk number 37: MEDIPS.Rnw:882-883 (eval = FALSE) ################################################### ## frames.frame500.step250=read.table(file="frames.chr22.meth.txt", header=T) ################################################### ### code chunk number 38: MEDIPS.Rnw:901-903 ################################################### file=system.file("extdata", "MeDIP_DE_chr22.txt", package="MEDIPS") TREAT.SET=MEDIPS.readAlignedSequences(BSgenome="BSgenome.Hsapiens.UCSC.hg19", file=file, numrows=863054) ################################################### ### code chunk number 39: MEDIPS.Rnw:906-907 ################################################### TREAT.SET=MEDIPS.genomeVector(data=TREAT.SET, bin_size=50, extend=400) ################################################### ### code chunk number 40: MEDIPS.Rnw:909-910 ################################################### TREAT.SET=MEDIPS.getPositions(data=TREAT.SET, pattern="CG") ################################################### ### code chunk number 41: MEDIPS.Rnw:912-913 ################################################### TREAT.SET=MEDIPS.couplingVector(data=TREAT.SET, fragmentLength=700, func="count") ################################################### ### code chunk number 42: MEDIPS.Rnw:915-916 ################################################### TREAT.SET=MEDIPS.calibrationCurve(data=TREAT.SET) ################################################### ### code chunk number 43: MEDIPS.Rnw:918-919 ################################################### TREAT.SET=MEDIPS.normalize(data=TREAT.SET) ################################################### ### code chunk number 44: MEDIPS.Rnw:930-932 ################################################### file=system.file("extdata", "Input_StemCells_chr22.txt", package="MEDIPS") INPUT.SET=MEDIPS.readAlignedSequences(BSgenome="BSgenome.Hsapiens.UCSC.hg19", file=file, numrows=352756) ################################################### ### code chunk number 45: MEDIPS.Rnw:935-936 ################################################### INPUT.SET=MEDIPS.genomeVector(data=INPUT.SET, bin_size=50, extend=400) ################################################### ### code chunk number 46: MEDIPS.Rnw:947-948 (eval = FALSE) ################################################### ## MEDIPS.exportWIG(file="output.rpm.input.WIG", data=INPUT.SET, raw=T, descr="INPUT.rpm") ################################################### ### code chunk number 47: MEDIPS.Rnw:951-953 (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) ################################################### ### code chunk number 48: MEDIPS.Rnw:959-960 (eval = FALSE) ################################################### ## INPUT.SET=MEDIPS.getPositions(data=INPUT.SET, pattern="CG") ################################################### ### code chunk number 49: MEDIPS.Rnw:963-965 (eval = FALSE) ################################################### ## cr.input=MEDIPS.coverageAnalysis(data=INPUT.SET, extend=400, no_iterations=10) ## MEDIPS.plotCoverage(cr.input) ################################################### ### code chunk number 50: MEDIPS.Rnw:968-969 (eval = FALSE) ################################################### ## er.input=MEDIPS.CpGenrich(data=INPUT.SET) ################################################### ### code chunk number 51: MEDIPS.Rnw:974-975 (eval = FALSE) ################################################### ## INPUT.SET=MEDIPS.couplingVector(data=INPUT.SET, fragmentLength=700, func="count") ################################################### ### code chunk number 52: MEDIPS.Rnw:978-979 (eval = FALSE) ################################################### ## INPUT.SET=MEDIPS.calibrationCurve(data=INPUT.SET) ################################################### ### code chunk number 53: MEDIPS.Rnw:982-985 (eval = FALSE) ################################################### ## png("CalibrationPlotINPUT.png") ## MEDIPS.plotCalibrationPlot(data=INPUT.SET, linearFit=F, plot_chr="chr22") ## dev.off() ################################################### ### code chunk number 54: MEDIPS.Rnw:1009-1010 ################################################### diff.meth=MEDIPS.methylProfiling(data1=CONTROL.SET, data2=TREAT.SET, input=INPUT.SET, frame_size=500, select=2) ################################################### ### code chunk number 55: MEDIPS.Rnw:1091-1092 (eval = FALSE) ################################################### ## write.table(diff.meth, file="diff.meth.txt", sep="\t", quote=F, col.names=T, row.names=F) ################################################### ### code chunk number 56: MEDIPS.Rnw:1097-1098 (eval = FALSE) ################################################### ## diff.meth=read.table(file="diff.meth.txt", header=T) ################################################### ### code chunk number 57: MEDIPS.Rnw:1139-1140 ################################################### diff.meth.control=MEDIPS.selectSignificants(frames=diff.meth, control=T, input=T, up=2, p.value=0.0001, quant=0.99) ################################################### ### code chunk number 58: MEDIPS.Rnw:1176-1177 (eval = FALSE) ################################################### ## diff.meth.control.merged=MEDIPS.mergeFrames(frames=diff.meth.control) ################################################### ### code chunk number 59: MEDIPS.Rnw:1191-1192 (eval = FALSE) ################################################### ## write.table(diff.meth.control, file = "DiffMethyl.Up.hESCs.bed", sep = "\t", quote = F, row.names = F, col.names = F) ################################################### ### code chunk number 60: MEDIPS.Rnw:1201-1203 ################################################### file=system.file("extdata", "hg19.chr22.txt", package="MEDIPS") diff.meth.control.annotated=MEDIPS.annotate(diff.meth.control, anno=file) ################################################### ### code chunk number 61: MEDIPS.Rnw:1224-1225 ################################################### length(unique(diff.meth.control.annotated[,4])) ################################################### ### code chunk number 62: MEDIPS.Rnw:1233-1234 ################################################### diff.meth.treat=MEDIPS.selectSignificants(frames=diff.meth, control=F, input=T, up=2, down=0.5, p.value=0.0001, quant=0.99) ################################################### ### code chunk number 63: MEDIPS.Rnw:1236-1237 (eval = FALSE) ################################################### ## write.table(diff.meth.treat, file = "DiffMethyl.Up.DE.bed", sep = "\t", quote = F, row.names = F, col.names = F) ################################################### ### code chunk number 64: MEDIPS.Rnw:1240-1242 ################################################### file=system.file("extdata", "hg19.chr22.txt", package="MEDIPS") diff.meth.treat.annotated=MEDIPS.annotate(diff.meth.treat, anno=file) ################################################### ### code chunk number 65: MEDIPS.Rnw:1245-1246 ################################################### length(unique(diff.meth.treat.annotated[,4]))