Bioconductor: DMRcate

[Tim Peters]

Hello epigeneticists, bioinformaticians, and other interested parties,

I've written a small Bioconductor package called DMRcate that analyses 450k data. It extracts and ranks the most differentially methylated regions de novo (i.e. agnostic to its annotation) from a pair of labelled phenotypes, using kernel density modelling. All you need as input is your matrix of M-values with probe IDs as rownames, and a limma style model.matrix denoting phenotype. There are a number of parameters you can experiment with, including kernel bandwidth p-value cutoffs, beta-value cutoffs and the like. You can find it here: http://bioconductor.org/packages/2.14/bioc/html/DMRcate.html

Extra functionality includes:

• Extraction of most variable and most hypermethylated regions from non-labelled samples
  
• Optional pre-filtering step allowing removal of probes potentially confounded by SNPs and cross-hybridisation
• BedGraph generation
  
• Plotting of individual DMRs (by beta value)

Any questions or feedback are more than welcome. If you like you can contact me at <Tim.Peters at csiro.au>.

Regards,

Tim

Dr. Tim Peters

Postdoctoral Fellow - Bioinformatics and Statistics

CSIRO Computational Informatics

Ph: +61 2 9325 3266 

Riverside Life Sciences Centre, 11 Julius Avenue, North Ryde, NSW 2113, Australia

[Tim Peters]

A further update:

A few users have been asking me about the correct way to set up limma using unpaired data for DMRcate. Firstly, the design matrix should not have an intercept. Say we have 8 unpaired samples constituted of 4 in each phenotype where the matrix of M-values myUnpairedMs has phenotype 1 values in columns 1,3,5 and 7 and phenotype 2 in columns 2, 4, 6, and 8. Then the correct procedure would be:

> type <- factor(rep(c("1","2"), times=4))

> design <- model.matrix(~0+type)

> design

  type1 type2

1     1     0

2     0     1

3     1     0

4     0     1

5     1     0

6     0     1

7     1     0

8     0     1

attr(,"assign")

[1] 1 1

attr(,"contrasts")

attr(,"contrasts")$type

[1] "contr.treatment"

> myannotation <- annotate(myUnpairedMs, analysis.type = "differential", design = design, coef = 1, diff.metric = "FC", paired = FALSE)

[Tim Peters]

One caveat: if you have ggplot 2 loaded in your namespace you will get this error when you call annotate():

Error in geom$default_pos : $ operator is invalid for atomic vectors

So for the time being, use DMRcate:::annotate(), and I will put a more permanent fix in the next version...

Tim

[David Martino]

Thanks Tim,

if this were a normal limma design, I would imagine you would need to set up a contrast matrix to estimate the coefficients for the difference between type 1 and type 2 samples if you didn't have an intercept term. It seems not to be the case with DMRcate? In your function annotate with analysis.type='differential', coefficient 1 in the example here would estimate the mean log-methylation for the type 1 samples would it not? Thus to draw the difference between type 1 and type 2 you would either need a contrast matrix, or fit an intercept term by which the first coefficient would be the log methylation estimates for type 1 samples and the second coefficient would be the difference between type 1 and type 2, as in the examples on p40 of the limma manual.

have I got this wrong?

Thanks

David.

[Tim Peters]

Hi David,

The DMRcate pipeline does the contrasting itself internally for unpaired analyses. The contrast matrix is created within the annotate() function, specifically at line 43 of the source. The next version of DMRcate will have a bit more explicit user control for multi-level factors and the like. But for the time being it's easy enough to hack a new "weights" vector into your "annot" object with the weights you want.

Regards,

Tim

[Benjamin Green]

Hi Tim,

I am working on creating an unpaired DMR analysis using human 450K data. I had the program up and running, but noticed your group post about running unpaired data. I have removed the intercept from my design matrix, and have created a limma-style contrast matrix to read into the cpg.annotate function, but am continuously getting the error: "Error: coef %in% colnames(cont.matrix) is not TRUE". I am not sure what exactly the issue is and was wondering if you had an idea.

Thanks!

Ben

[Rachel Edgar]

Hi Tim,

We are interested in using DMRcate on our methylation data. Seems like a great package. We were able to install DMRcate fine in R version  3.1.0 perfectly fine and it was very easy to run!  We have a few servers which will be able to run DMRcate on bigger datasets but we have been getting the error

Warning message:

package ‘DMRcate’ is not available (for R version 3.1.1)

I was wondering if anyone you know of has had this issue and if you have a quick fix. Other then reverting to R 3.1.0  which would be a bit of a head ache.

Any help you can offer will be appreciated!

Thanks

Rachel Edgar

Research Assistant, Kobor Lab

[Roderick Slieker]

Hi Rachel,

The error you get is not caused by the package. You have to update bioconductor to a newer version.

source("http://bioconductor.org/biocLite.R")
biocLite("BiocUpgrade")     ## R version 2.15 or later

Once all packages are updated you will be able to install it.

Good luck,

Roderick

Op vrijdag 14 november 2014 19:36:15 UTC+1 schreef Rachel Edgar:

[Rachel Edgar]

Thanks for the reply Roderick. My bioconductor was up to date but you were right I had out of date packages. After I ran:

biocLite()

BioC_mirror: http://bioconductor.org

Using Bioconductor version 3.0 (BiocInstaller 1.16.1), R version 3.1.1.

Old packages: 'DESeq2', 'edgeR', 'VariantAnnotation'

Update all/some/none? [a/s/n]: a

DMRcate installed perfectly.

Thanks for your help!

Rachel

[Tim Peters]

Hi Ben,

Thanks for your post, I ought to be checking this forum more regularly than I have recently been doing. The newer versions of DMRcate (since August) have changed the required approach for unpaired data because the earlier version wasn't so intuitive. If your design matrix doesn't have an intercept, you must now pass a contrast matrix to cpg.annotate as well, like you would do for using contrasts.fit() in limma. For simple 2 group data this is as easy as specifying cont.matrix=c(1, -1), coef=1 as arguments to cpg.annotate() (along with your nx2 design matrix).

You can now also pass more complex contrast matrices to cpg.annotate(), but it will only take 1 coefficient at a time. 

Good Luck,

Tim

[Tim Peters]

Thank you Roderick! And glad it worked for you, Rachel.

[Methylnick]

Hi Tim,

I have not been working on methylation for a little while but upon getting back into it I have come across your library, it is a nice package.

I was wondering if you have managed to publish a paper on it?

[Tim Peters]

Hi Methylnick (and others),

A provisional version of the DMRcate paper is now available here: http://www.epigeneticsandchromatin.com/content/8/1/6/abstract. 

Any questions or feedback are more than welcome, either here or by email.

Regards,

Tim

[Methylnick]

That is great, many thanks Tim for the heads up and congrats on publishing the work. 

Nick.

[Smeeta Shrestha]

 Hi !

Why do I get a "X" prepended to my column header. I am using DMRcate and use column header as "55-Normal" and "55-Tumor", but the name gets changes to "X55.Normal"  " X55.Tumor".

Any help would be really appreciated.

Thank you

Smeeta Shrestha
Post doc
Centre for Cellular and Molecular Biology
India.

[Roderick Slieker]

R does not allow column names to start with a number and will always place "X" in front of it...

Roderick

Op Mon Feb 02 2015 at 06:11:45 schreef Smeeta Shrestha <smeetas...@gmail.com>:

--
You received this message because you are subscribed to the Google Groups "Epigenomics forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to epigenomicsfor...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Smeeta Shrestha]

Thank you Roderick for the information. Managed to run DMRcate with small sample data set.

smeeta

[Smeeta Shrestha]

Hi !
Thank you for earlier assistance. I am having another problem with DMRcate.

I think I should put down the study background a bit .

1. I have total of 24 individuals in my analysis.
2. 12 are case which as per the DMRcate demo I put as "Tumor"
3. 12 are control which as per the DMRcase demo I put as  "Normal".
4 My matrix contains 462,813 probes (rows) and 24 individuals column and 1 column for Probe ID.
5. I have labelled the sample heads but not the probe id head. ( attachment 1)
6. i am using the DMR script and get this error.

> patient <- factor(sub("-.*", "", colnames(myMs)))
> type <- factor(sub(".*-", "", colnames(myMs)))
> design <- model.matrix(~patient + type)
> myannotation <- cpg.annotate(myMs.noSNPs, analysis.type="differential",
+     design=design, coef=1)
Error in data.frame(ID = rownames(object), weights = weights, CHR = RSanno$chr,  :
  arguments imply differing number of rows: 419134, 419070

Please any assistance would be greatly appreciated.

regards
smeeta

[Pooja]

Hi Tim,

 

I recently started working with DMRs and I think DMRcate is a very nice package to work with. I have managed to learn the steps to find DMRs for the paired and unpaired data having categorical phenotypes. However in my case, I have a continuous phenotype (for example: lets say I want to know the effect of age on methylation). Could you help me know how can I work with continuous phenotypes?

 

> patient [1:5]

[1] X8261  X8164  X7342  X10318 X3949

700 Levels: X10080 X10108 X10115 X10118 X10130 X10136 X10138 X10146 ... X9485

 

> type [1:5]

[1] 0.99 1.22 1.27 1.55 1.05

 

How do I design and annotate for this?

 

Thanks a ton for your help!

 

Best,

Pooja

[Tim Peters]

Hi Pooja,

Thanks! DMRcate can definitely work with continuous response variables. All you have to do is set up the design matrix like you would for a linear regression in limma i.e.

design <- model.matrix(~type)

and then set coef=2 in the call to cpg.annotate().

Be warned that the $maxbetafc in the output matrix is less interpretable; it's a slope rather than a beta fold-change and therefore contingent on the range of your response variable. My advice is to look at the distribution of slopes and manually filter if you think this is important.

Cheers,
Tim

[Michael]

Hi Peter

I am having trouble with the annotation step. There is something going wrong with the matrix/df dimensions, but I cannot guess where the problem lies.

I am following the vignette example and do apply the close-to-SNP filtering. In any case, when I do the annotation step it seems to work in a first instance, but then there is an error: the output looks like this:

myann <- cpg.annotate(betas.filtered.log, analysis.type="differential", ,design=design.m, coef=dim(design.m)[2])
Your contrast returned 169081 individually significant probes. We recommend the default setting of pcutoff in dmrcate().
Error in data.frame(ID = rownames(object), weights = weights, CHR = RSanno$chr,  : 
  arguments imply differing number of rows: 438091, 438026

I tried the same with the data without applying the filtering SNP filtering step:
myann <- cpg.annotate(betas.log, analysis.type="differential",design=design.m, coef=dim(design.m)[2])
Your contrast returned 187097 individually significant probes. We recommend the default setting of pcutoff in dmrcate().
Error in data.frame(ID = rownames(object), weights = weights, CHR = RSanno$chr,  : 
  arguments imply differing number of rows: 484940, 484875

Both times, there is an offset of 65 rows in whatever objects the cpg.annotate() is trying put into a data.frame. Where could these 65 entrie come from and how can I prevent this?


FIY:
dim(betas.filtered.log)
[1] 438091     11
dim(betas.log)
[1] 484940     11

Any help would be appreciated

[Tim Peters]

Hi Michael,

My guess is that you've left the SNP probes in your beta matrix. The 450k array has 65 genotyping probes that are used for quality control for potential sample mix-ups. Their IDs are in the form "rsXXXXXX" instead of "cgXXXXXX". DMRcate has no use for them and so you should remove them from the matrix.

In addition, cpg.annotate() should only take M-values, never beta values (otherwise the beta fold change calculations will be wrong). 

Regards,

Tim

[Michael]

Hi Tim

Thanks for the quick answer! Your guess is spot-on right! I was breaking my head about that all day long altough a had read about the 65 control SNPs before. So removing those was the solution to the problems.

About the M-values: I was using the logged values of the beta's. This is what you call the M-values right? Maybe the variable name should be changed.

Just a little follow-up question: At what step would you recommend to remove the control SNPs: Before or after normalization?

Greetings

Michael

[Tim Triche, Jr.]

M-values are log2(M/U) with a small offset added to each which avoids positive and negative infinities/NaNs.  See Du et al 2010, http://www.biomedcentral.com/1471-2105/11/587, for discussion.  log(beta) is a nonsense quantity, do not use. It is arguably unwise to use m-values, too, unless you know exactly what you're doing and believe that fold-changes at ~15x coverage equivalence are a wise representation of biological changes. The latter statement is itself debatable. (M-values are terrific for performing downstream normalization due to their sensitivity!) I have been hemming and hawing for years as to whether they're as good for reproducible analysis, and should finally put up a manuscript on the subject in the near future, now that it's out-of-fashion and can be discussed rationally in the context of "will my results most likely be reproduced?"

I'd suggest funnorm for normalization, then rmSNPandCH() on the resulting GenomicRatioSet... it seems to be a good general-purpose approach.

Dr. Peters: this reminds me that I need to send you some patches.

Best,

another Tim

[Michael]

Hi Tim

Thanks for your explanations and literature hint! It helped understanding the logit (not log) transformation, which I was using as I was following the DMRcate vignette.

Just out of curiosity - why are you worried about reproducibility when using M-values? 

Best,

Michael

[Cameron Lassiter]

Hello Dr Peters,

I have been attempting to use your DMRcate package with using human 450K data and a continuous outcome variable. I am having difficulty interpreting my DMR finding, specifically the P -values returned. In my data set of 36 arrays, there are no significant ( lowest FDR is 0.3) CpG sites. The power is not very high so I call dmrcate() with pcutoff at 0.1. Doing this I get ~5,000 DMRs with p-values above 0.05 (low minpval ~ 10^-12, and means ~10^-5) . Some of these are in fact single probe regions, suggesting that something is off in my analysis. Hopefully you can help clarify what is happening here.

Many thanks,

Cameron Lassiter

P.S. Your package has worked well for me with categorical variables.

[Tim Peters]

Hi Cameron,

Since DMRcate gathers information locally from other probes, its post-smoothing per-CpG site test for effect is more sensitive than if probes are considered individually. So, to guard against this potential Type I error, I set the default in dmrcate() to include as many probes in DMR constituents as would be found if they were considered individually via limma (which is fairly conservative). 

Single probe "regions" are fine; it just means that only one probe's post-smoothed FDR has poked its head "above" the threshold you set in pcutoff. Either that or it's located in a genomically sparse section of the array, probewise. I usually filter the singletons out manually but made the decision to leave them in the results since they may be of interest to users. Feel free to set your own lower bound on how many probes you feel constitute a DMR.

pcutoff=0.1 is very generous in my opinion; and since the cutoff is made pre-aggregation then you'll get some DMRs with long stretches of negligible DM. I'd tighten it up considerably.

Best,

Tim

[Casey S]

Hi Tim,

The use of limma "under the hood" is what initially drew me to using DMRcate... it's a package I've used heavily for gene expression analysis. I have DNA methylation data derived from human peripheral whole blood, so investigating the methylation profiles of the various cellular components within is of interest. In gene expression, I've used limma and interaction designs (e.g. model.matrix(~ cell_prop_matrix * condition) and considering the various contrasts: cell1:treated, cell1:treated, etc.) to do so. I tried doing the same with DMRcate, but this results in lmFit throwing a... fit.

Is such interaction models possible in DMRcate?

Thanks!

[Tim Triche, Jr.]

If you felt like posting a reproducible example either here or to support.bioconductor.org, it is far more likely that someone would be able to help you

[yourela3]

Hi Tim,

I'm doing a sibling paired analysis and have  been  using cpg.annotate to find DMRs. I keep getting "Your contrast returned no individually significant probes. Set pcutoff manually in dmrcate() to return DMRs, but be warned there is an increased risk of Type I errors."

I've changed the pcutoff but nothing seems to work. I've attached my script below. I've converted the output of wateRmelon nasen method to M-values and I have included sibling status and diagnosis in the design matrix (whit an intercept) and have also added in surrogate variables. Any ideas on what could be the issue?

Thanks,

Alan

myannotation <- cpg.annotate(nasen.Mval, analysis.type="differential",contrasts=FALSE, design=modSv, pcutoff=0.05, coef=77)

[Tim Triche, Jr.]

Do you find significant probes if you perform the same fit with limma/lmFit and eBayes?

All that cpg.annotate is really doing is the following:

fit <- lmFit(object, design, ...)

#...

fit <- eBayes(fit)

tt <- topTable(fit, coef=coef, number = nrow(object))

nsig <- sum(tt$adj.P.Val < 0.05)

If nsig < 0 it's going to complain.  (You could argue that this 0.05 cutoff is a bug, since it's static... perhaps that's worth patching.)

[Tim Peters]

Hi Alan,

Following on from Tim Jr.,

Firstly, pcutoff is obsolete for cpg.annotate(), so it won't make a difference if you change it there. What's happening is that your contrast isn't returning any probes with an limma-derived fdr < 0.05; so it means that your effect is relatively subtle. myannotation should still be a legitimate object to pass to dmrcate(, but if you want to see where in the genome your contrast has the greatest effect, you'll need to specify pcutoff using dmrcate(). These are a different set of fdrs because they are derived from the kernel smoother and hence are likely to be generally smaller than myannotation$indfdr. Start with, say, 1e-4 and work your way upwards.

Cheers,

Tim

[Tim Peters]

Hi Casey,

Like Dr. Triche said, a reproducible example would help enormously. DMRcate should be able to take any design matrix/contrast/term valid for limma.

Cheers,

Tim

[yourela3]

Thanks Tim and Tim, 

I had completely overlooked the fact you can specify the pcutoff in dmrcate() and was changing it using cpg.annotate(). 

Alan

[Kristian Almstrup]

Dear Tim ao.

Is there any easy way to add legends to a DMR.plot?

I know that I can use legend() from the plot function after DMR.plot but then I have to redefine the legend incl. colours etc.

Thanks for a nice package.

Cheers,

Kristian

[Tim Peters]

Hi Kristian,

Not really, I'm afraid. Automated legend additions are left out for the same reason as they are in most R packages; good aesthetic and avoiding obscuring parts of the plot are difficult without visualising the plot first. What you can do, since DMR.plot() has an ellipsis, is expand the space by adding xlim and ylim arguments and wrap it up with your own custom call to legend(), utilising the extra space.

Best,

Tim

[Alder Olivia]

Hi Smeeta, 

I'm getting the same error, did you ever find a solution? 

> myannotation_oa <- cpg.annotate(myMs, analysis.type="differential",

+     design=design, contrasts=FALSE, coef=39)

>Error in data.frame(ID = rownames(object), weights = weights, CHR = RSanno$chr,  : 

  arguments imply differing number of rows: 26486, 25014

Many thanks 

Olivia 

[Alder Olivia]

Figured it out. 

[yourela3]

Hi everyone,

I'm trying to use dmrcate() on the replication phase of my study  but have encountered some issues. I've only taken  probes forward from my discovery phase with an unadjusted p-value<0.05 leaving me with approx 23,000 probes.  I encounter the error pasted below. This doesn't occur when I run analysis with all the probes that were remaining in my discovery phase after pre-processing. Any ideas? 

Error in approx(x = x, y = i, xout = xout, rule = 2) :

  need at least two non-NA values to interpolate

Thanks,

Alan

[Tim Peters]

Hi Alan,

Please don't pre-filter your probes on effect size before running DMRcate; this will bias the method. DMRcate is set up to receive as many valid CpG sites as possible (significantly differential or not) and will process them under-the-hood. The more CpG sites there are, the better the results will be. 

Cheers,

Tim

[Pär Larsson]

Hi Alder,

I'm having the same issue, may I ask how you got past it?

-Pär

[Tim Peters]

A significantly updated version of DMRcate is now available. I have created a new thread detailing the changes:

https://groups.google.com/forum/#!topic/epigenomicsforum/G3yO92pEqzQ

Tim 

[liyufan...@gmail.com]

Hello Michael ,

    I have met the same problems. But I don't know how to remove the SNP.Can help me?

在 2015年5月13日星期三 UTC+8下午11:14:36，Michael写道：

[liyufan...@gmail.com]

Hello Tim Peters,

    I have met the same problems. But I don't know how to remove the SNP. Can help me?Can you give me some hints?

在 2015年5月14日星期四 UTC+8下午4:12:31，Tim Peters写道：

[liyufan...@gmail.com]

Hello Peters,

I have used package called DMRcate to  analyse 450k data. I want to find the gene which associates the DMRs.And I met some questions.The outputs include gene_assoc, group, hg19coord, no.probes, minpval, meanpval and maxbetafc,when  I used the old version to analyse 450k data.

But when I update the package,the outputs became coord, no.cpgs, minfdr, Stouffer, maxbetafc,meanbetafc.There is no result of "gene_assoc".I want to find the gene names associating "coord",can you help me ? Can you please tel me how to associate the gene by using the newest DMRcate packages.

There follow the output results of the newest DMRcate.

 coord no.cpgs minfdr Stouffer  maxbetafc  meanbetafc
63999  chr6:33156164-33181870     265      0        0 -0.5008031 -0.02648790
63997  chr6:33128825-33149777     150      0        0  0.4176126  0.08611966
63917  chr6:32144195-32161004     128      0        0 -0.2574513 -0.03184096
63914  chr6:32114490-32123701     124      0        0 -0.4377015 -0.06195576
63889  chr6:31935801-31940855     101      0        0 -0.1555205 -0.02401999
12564 chr11:31817810-31841980     100      0        0 -0.4611059 -0.17113506

There follow the output results of the older DMRcate. 
gene_assoc                  group        hg19coord no.probes minpval    meanpval  maxbetafc
LOC100132215,OTX1 Body,TSS200,TSS1500 chr2:63273684-63287288     88  0 1.53809e-08 -0.1045803

Thank you !

在 2014年3月17日星期一 UTC+8上午8:38:59，Tim Peters写道：