Mercurial > repos > rnateam > methylkit

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools/methylkit commit 14f0b39f64982773ef0367379b915f742eabcc1b
author rnateam
date Wed, 21 Dec 2016 17:30:57 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a8705df7c57f
1 library("methylKit")
2
3 file.list = list("input_test1.myCpG.txt", "input_test2.myCpG.txt", "input_control1.myCpG.txt", "input_control2.myCpG.txt")
4
5 myobj=methRead( file.list,
6 sample.id=list("test 1","test 2","control 1","control 2"),assembly="hg18",pipeline="amp",treatment=c(1,1,0,0))
7
8 pdf('output_statistics.pdf')
9 for (obj in myobj){
10 getMethylationStats(obj,plot=TRUE,both.strands=FALSE)
11 getCoverageStats(obj,plot=TRUE,both.strands=FALSE)
12 }
13 devname = dev.off()
14
15 # unite function
16 methidh = unite(myobj)
17
18 pdf("output_correlation.pdf")
19 getCorrelation(object = methidh, plot=TRUE, method = "pearson")
20 devname = dev.off()
21
22 # differential methylation
23 myDiff = calculateDiffMeth(methidh, overdispersion="none",
24 adjust="SLIM", effect="wmean", test="Chisq",
25 slim=FALSE, weighted.mean=FALSE)
26
27 bedgraph(myDiff, file.name="output_myDiff.bedgraph", col.name="meth.diff",
28 unmeth=FALSE, log.transform=FALSE, negative=FALSE, add.on="")
29
30 MethPerChr = diffMethPerChr(myDiff, plot=FALSE,
31 qvalue.cutoff=0.01,
32 meth.cutoff=25)
33 write.table(MethPerChr, sep="\t", row.names=FALSE, quote=FALSE, file="output_MethPerChr.tsv")
34
35 MethylDiff = getMethylDiff(myDiff, difference=25,
36 qvalue=0.01, type="all")
37 bedgraph(MethylDiff, file.name="output_MethylDiff.bedgraph", col.name="meth.diff",
38 unmeth=FALSE,log.transform=FALSE,negative=FALSE,add.on="")
39
40 pdf( "output_clustering.pdf" )
41 methClust = clusterSamples(methidh, dist="correlation", method="ward")
42 devname = dev.off()
43
44 pdf( "output_PCA.pdf" )
45 PCASamples(methidh)
46 devname = dev.off()
47
48 ## methSeg works for methylRaw or methylDiff with resolution region,
49 ## so methylBase has to be tiled before
50 tileRaw = tileMethylCounts(myobj[[1]])
51 tileBase = tileMethylCounts(methidh)
52 tileDiff = calculateDiffMeth(tileBase)
53
54 ## methseg generates Granges
55 segRaw = methSeg(tileRaw, diagnostic.plot = FALSE)
56 segDiff = methSeg(tileDiff, diagnostic.plot = FALSE)
57
58 ## and can be exported as BED
59 methSeg2bed(segments = segRaw, filename = "output_seg_raw.bed")
60 methSeg2bed(segments = segDiff, filename = "output_seg_diff.bed")