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.