-
Notifications
You must be signed in to change notification settings - Fork 1
/
MethDMRByMethyAnalysis.R
73 lines (68 loc) · 3.04 KB
/
MethDMRByMethyAnalysis.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# DMR analysis by methyAnalysis package
# It require the data analysis from idat
# Slide-window smoothing and Differential methylation test
# plotMethylationHeatmapByGene could be provided in the last
# Data: Oct/12/2016
# Email: [email protected]
source("http://bioconductor.org/biocLite.R")
biocLite("methyAnalysis")
biocLite("coMET")
library("methyAnalysis")
## load the example data
data(exampleMethyGenoSet)
## show MethyGenoSet class
slotNames(exampleMethyGenoSet)
MethyGenoSet(rowRanges, exprs, methylated, unmethylated, detection = NULL, pData = NULL,annotation = "", universe = NULL, assays=NULL)
options(width=50)
library(methyAnalysis)
data(exampleMethyGenoSet)
## show MethyGenoSet class
slotNames(exampleMethyGenoSet)
# showClass('MethyGenoSet')
## get chromosome location information
head(rowRanges(exampleMethyGenoSet))
## retrieve methylation data
dim(exprs(exampleMethyGenoSet))
## Sample information
colData(exampleMethyGenoSet)
###################################################
### code chunk number 3: Slide-window smoothing of the DNA methylation data
###################################################
methyGenoSet.sm <- smoothMethyData(exampleMethyGenoSet, winSize = 250)
## winsize is kept as an attribute
attr(methyGenoSet.sm, 'windowSize')
###################################################
### code chunk number 4: Differential methylation test
###################################################
## get sample type information
sampleType <- colData(exampleMethyGenoSet)$SampleType
## Do differential test
allResult <- detectDMR.slideWin(exampleMethyGenoSet, sampleType=sampleType)
head(allResult)
###################################################
### code chunk number 5: Identifying DMR
###################################################
## Identify the DMR (Differentially Methylated Region) by setting proper parameters.
## Here we just use default ones
allDMRInfo = identifySigDMR(allResult)
names(allDMRInfo)
###################################################
### code chunk number 6: Identifying DMR
###################################################
## Annotate significant DMR info
DMRInfo.ann <- annotateDMRInfo(allDMRInfo, 'TxDb.Hsapiens.UCSC.hg19.knownGene')
## output the DMR information
export.DMRInfo(DMRInfo.ann, savePrefix='testExample')
###################################################
### code chunk number 7: Export data for external tools (eval = FALSE)
###################################################
## ## output in IGV supported "gct" file
## export.methyGenoSet(exampleMethyGenoSet, file.format='gct', savePrefix='test')
## ## output in BigWig files
## export.methyGenoSet(exampleMethyGenoSet, file.format='bw', savePrefix='test')
###################################################
### code chunk number 9: heatmap_phenotype
###################################################
plotMethylationHeatmapByGene('6493', methyGenoSet=exampleMethyGenoSet, includeGeneBody=TRUE,
phenoData=colData(exampleMethyGenoSet), genomicFeature='TxDb.Hsapiens.UCSC.hg19.knownGene',
newPlot=FALSE)