RADMeth DMR in IGV: Looks…


Looks like I may have spoken too soon. I loaded up the RADMeth produced DMRs in IGV and the DMRs didn’t even remotely line up to sites with methylation in the Day 10 samples. This may be an artifact of the merging process, or something being wonky with the model. Off the top of my head, potential oddities for the model may include a lack of nearby neighbors for the p-value adjusting step, or some oddity in trying to calculate DMRs from a very sparse set of loci.

I’m going to try to figure out how to load up either the regression or the adjustment steps in IGV and see if that shines any light on the problem.


IGV session XML file

RADMeth output for Day 10:…

RADMeth output for Day 10:

Over the last weekend I ran RADMeth on the day 10 samples I had converted and got some results! This was a little surprising time wise, as the RADMeth manual suggests you run the program on a cluster with “a few hundred” nodes available. Maybe that was just wishful thinking on their part, or something else? I’m not honestly sure.

The program doesn’t really have any screen output, so I’ll just paste the little bit of code and the output here for brevity’s sake. The merging and whatnot was done in the notebook found here

First, we run the regression portion of the program with

radmeth regression -factor trt designmatrix6.txt Day10-merged.meth > cpgs.bed

That outputs a cpgs.bed file that looks like:


The first four columns are location, strand, and context information, and the 6th-9th columns are total coverage and methylated counts at that location. The 5th column is the interesting column, representing the value obtained from a log-likelihood ratio test. The MethPipe developers suggest not using this p-value directly, and instead adjusting based upon neighboring values using the adjust argument to radmeth below.

radmeth adjust -o cpgs.adj cpgs.bed

which results in:


The first 5 and last 4 columns are the same as above, with columns 6 and 7 representing modified p-values based on neighboring sites, and an FDR corrected value (I need to look in to what method they use for FDR correction.

Next, the developers suggest we merge the loci in to regions, I went ahead and did this, though we may be more interested in the loci themselves, so both datasets will be avaiable. Merging is done via

radmeth merge -o cpgs.output cpgs.adj

which gives us:


Hey, that’s not an empty file, which is surprising given that MACAU + q-value FDR correction resulted in no significant DMRs. Interesting. I’ll throw those in to IGV and update this shortly. Columns are scaffold, start/stop location, context, number of methylated cytosines in the region, and log-likelihood ratio value.

Output files will be located in owl once they’re finished uploading.

RADMeth uses a beta-binomial based regression model so we preserve the ability to directly model count data, but lose the ability to control for some biological covariates and model genetic relatedness based covariance which will be something to consider.