I am trying to use RnBeads to analyze DNA-methylation data, but I already have a matrix of beta-values. I cannot find a way to import data in that format and launch a differential analysis. Has anyone encountered the same problem?
If you have normalised Beta values, and you trust the normalisation method (i.e. there's no way to get the IDAT files), then you can simply read that into R, convert to M values:
m <- log2(beta/(1 - beta))
Then perform a differential methylation analysis using Limma
I would actually prefer using the beta values rather than the transformed values (that is what I do within COHCAP), since there can (rarely) be results that actually show a noticable percent methylation difference that is lost in the transformation. You probably also want to keep a |delta-beta| filter, to avoid results with small differences that are close to 0 and 1 (which is probably the more frequent issue).
For RnBeads, I think you can write R scripts for more specific steps, but I think they kind of expect you to analyze data with normalization performed within RnBeads. I think the last time I saw something that was totally off was SWAN normalization with an EPIC array (but one weird result for one project isn't a good predictor of introduction of possible artifacts in other projects), but I also don't think that was with RnBeads (although I have to apologize that I don't know for sure).
Nevertheless, I'll also add an answer about possible ways to start analysis in RnBeads with the raw data as a separate answer.
This answer is three years old, and I'm struggling to see the point of this comment other than a plug for your package? The point of using M values is that Limma fits a homoskedastic model, and M values are homoskedastic. Beta values are heteroskedastic mostly in the <0.2 & >0.8 regions [1], and thus using them to fit a model in Limma would break a fundamental assumption, and result in more false positives. That's the rationale for my suggestion, and like any quantitative analysis, a magnitude of change threshold is sensible. Rnbeads uses Limma under the hood here.
I didn't notice the answer was from this long ago - sorry about that. I also appreciate you adding links for the discussion. However, I think the answer may be useful to some people.
I think having multiple freely available methods is a good thing. In other words, I don't want anybody to use any one strategy and immediately take those results to publication (which is part of what I was communicating in one of my COHCAP comments). So, I certainly don't want them to only use COHCAP: I am suggesting that either using COHCAP or their own code for analysis of beta-values may be worth as considering as another possibility.
For example, with BS-Seq, I would usually test both COHCAP and methylKit, and sometimes I would even tell people I thought methylKit was a better fit for their specific project.
While I am familiar with many people using M-values, I think there is something important about looking for candidates to further characterize (although I apologize that I may not be explaning myself precisely). For example, if you can use the default settings in COHCAP, I think the methylated and unmethylated thresholds can have more of an effect than the p-value / FDR filter. However, I think that we agree that using limma with M-values may further increase sensitivity (beyond using limma with beta values, which would probably be more sensitive than COHCAP). So, I am OK with using M-values, if you have a |delta-beta| filter.
Also, I think it sounds like we might agree that there can be value in applying a |delta-beta| filter of 0.2 (and/or using the methylated / unmethylated thresholds, such as in COHCAP, where the default would look for one group with a beta value > 0.7 and the other group with a beta value < 0.3), even if the delta-beta and p-values are calculated separately.
For some reason, I thought RnBeads had a delta-beta threshold, but it looks like I remembered that wrong (or am just having difficulty finding it). However, it looks like they provide beta values for visualization in the resulting reports, so that you could add extra filters outside of RnBeads (in that example, I believe mean.diff is equivalent to delta-beta).
As a warning, this code that I wrote may be slightly outdated.
However, if you have the raw .idats (which I think is important: if possible, I think testing different normalizations can give you a better feel for the robustness of your results), I think this code using a benchmark (described here) with some methods may also help as a starting point?
http://sourceforge.net/projects/cohcap/files/Protocol_Exchange_Example.zip
Otherwise, I think RnBeads has some good tutorials:
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
What kind of analysis do you want to perform? Differential Methylation?
Yes, sorry I didn't specify that.
Hi Nick,
I am analysing WGBS data, could you please guide me how to get matrix of normalised beta-value ?
Please do not refresh old threads with new questions. First use google and the search function and if this does not help, open a new question.
Thanks for your polite reply. I tried your second suggestion and opened a new thread A: WGCNA analysis for bisulfite data. and waited for 8 weeks. I refreshed the old thread when I couldn't get help there.
If your question does not get responses please consider adding more background information and the efforts you took to solve your problem to your question. That motivates people to invest time and effort into your question. You can also bump a question back to the top of the page by just editing it without making changes. Sometimes questions have simply been missed by those who have the expertise to answer.
Yes - that is a good point.
While I apologize for not seeing this before, I believe RnBeads is designed for Illumina Arrays (not BS-Seq).
Strictly speaking, I would say BS-Seq has a "percentage methylation" value (usually represented between 0 and 100%) and "beta" (between 0 and 1) is more frequently used for Illumina Methylation Arrays. If you have a minimum coverage value (such as at least 10x or 20x reads), I think using COHCAP can be useful, and some other methods that allow that input may be OK (but I have only done that with targeted BS-Seq data or RRBS data: a WGBS data table would be much larger).
Also, if you use annotations for region analysis, you'll be skipping a lot of WGBS sites (although I believe RnBeads has a tiling window option for analysis, not just pre-defined CpG Islands and/or promoter regions).