I am (maybe foolishly) trying to carry out a 5-way analysis in limma on some 450K methylation data.
I have carried out unsupervised clustering on methylation data and derived group calls for samples. I have then carried out this analysis to try and define probes which significantly define each cluster - that is probes which are differentially methylated between the clusters.
Let groupcall
be some vector of cluster calls from 1 to 5 and Mval
an M-value matrix
f <- factor(groupcall, levels = c("Group1", "Group2", "Group3", "Group4", "Group5"))
design <- model.matrix(~0+f)
colnames(design) <- c("Group1", "Group2", "Group3", "Group4", "Group5")
fit <- lmFit(Mval, design)
contrast.matrix <- makeContrasts(Group2-Group1,
Group3-Group2,
Group4-Group3,
Group5-Group4,
Group5-Group1,
Group5-Group2,
Group5-Group3,
Group4-Group1,
Group4-Group2,
Group3-Group1,
levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
results.bh <- decideTests(fit2, method = "global", adjust.method = "BH", p.value = 0.001)
Question 1: Is this a suitable way to approach this problem in terms of model design and fitting of model? I am ideally looking for probes that describe a cluster.
Having obtained results.bh
I see that it is a table of 1,0,-1 denoting significance
TestResults matrix
Contrasts
Group2 - Group1 Group3 - Group2 Group4 - Group3 Group5 - Group4 Group5 - Group1 Group5 - Group2 Group5 - Group3 Group4 - Group1 Group4 - Group2 Group3 - Group1
cg00009292 1 0 1 0 1 1 1 1 0 0
cg00044796 1 -1 1 -1 -1 -1 0 0 -1 -1
cg00064255 0 -1 1 0 0 -1 0 0 0 0
cg00069860 0 -1 0 1 1 1 1 0 0 0
cg00172812 1 -1 0 -1 -1 -1 0 0 -1 0
Question 2: How do I then go about extracting the cluster-specific probes?
Any guidance you can offer is appreciated.
Have you tried chAmp?
Thanks for replying. ChAmp can only handle differential methylation analysis between two groups. This is why I'm using limma itself rather than dmpFinder() in minfi or any of the other functions available which either require simpler models or don't implement my model properly.
Hello, any update about your problem? I would like to try the same approach with limma >:
Hi @yura.grabovska any update about your problem? I am interested to use the similar approach, kindly help.