diff test-data/test.r @ 0:a8705df7c57f draft default tip

planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/methylkit commit 14f0b39f64982773ef0367379b915f742eabcc1b

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/test.r	Wed Dec 21 17:30:57 2016 -0500
@@ -0,0 +1,60 @@
+library("methylKit")
+
+file.list = list("input_test1.myCpG.txt", "input_test2.myCpG.txt", "input_control1.myCpG.txt", "input_control2.myCpG.txt")
+
+myobj=methRead( file.list,
+                sample.id=list("test 1","test 2","control 1","control 2"),assembly="hg18",pipeline="amp",treatment=c(1,1,0,0))
+
+pdf('output_statistics.pdf')
+for (obj in myobj){
+  getMethylationStats(obj,plot=TRUE,both.strands=FALSE)
+  getCoverageStats(obj,plot=TRUE,both.strands=FALSE)
+}
+devname = dev.off()
+
+# unite function
+methidh = unite(myobj)
+
+pdf("output_correlation.pdf")
+getCorrelation(object = methidh, plot=TRUE, method = "pearson")
+devname = dev.off()
+
+# differential methylation
+myDiff = calculateDiffMeth(methidh, overdispersion="none",
+                           adjust="SLIM", effect="wmean", test="Chisq",
+                           slim=FALSE, weighted.mean=FALSE)
+
+bedgraph(myDiff, file.name="output_myDiff.bedgraph", col.name="meth.diff",
+         unmeth=FALSE, log.transform=FALSE, negative=FALSE, add.on="")
+
+MethPerChr = diffMethPerChr(myDiff, plot=FALSE,
+                            qvalue.cutoff=0.01,
+                            meth.cutoff=25)
+write.table(MethPerChr, sep="\t", row.names=FALSE, quote=FALSE, file="output_MethPerChr.tsv")
+
+MethylDiff = getMethylDiff(myDiff, difference=25,
+                               qvalue=0.01, type="all")
+bedgraph(MethylDiff, file.name="output_MethylDiff.bedgraph", col.name="meth.diff",
+         unmeth=FALSE,log.transform=FALSE,negative=FALSE,add.on="")
+
+pdf( "output_clustering.pdf" )
+methClust = clusterSamples(methidh, dist="correlation", method="ward")
+devname = dev.off()
+
+pdf( "output_PCA.pdf" )
+PCASamples(methidh)
+devname = dev.off()
+
+## methSeg works for methylRaw or methylDiff with resolution region,
+## so methylBase has to be tiled before
+tileRaw = tileMethylCounts(myobj[[1]])
+tileBase = tileMethylCounts(methidh)
+tileDiff = calculateDiffMeth(tileBase)
+
+## methseg generates Granges
+segRaw = methSeg(tileRaw, diagnostic.plot = FALSE)
+segDiff = methSeg(tileDiff, diagnostic.plot = FALSE)
+
+## and can be exported as BED
+methSeg2bed(segments = segRaw, filename = "output_seg_raw.bed")
+methSeg2bed(segments = segDiff, filename = "output_seg_diff.bed")