Hi Biostars community

I am working with RRBS data and currently I am trying the methylKit package to do pairwise comparisons between groups from a time course experiment.

I have encountered an issue with adding co-variates to the differential methylation analysis and if anyone knows what I am doing wrong I would be thankful if you could help me out.

From the original filtered and merged methylBase-object (contains all my 27 samples), I use reorganize() to select only two time points for the differential methylation test.

    meth =reorganize(x,sample.ids=c("1_1_2", "1_1_4", "1_2_2", "1_2_4", "10_1_2", "10_1_4", "10_2_2", "10_2_4"),treatment=c(0,0,0,0,1,1,1,1) )

covariates = data.frame(bisulfitelibs = c("d","d","d","d","a","c","a","b"))`

myDiff_withoutcovariates=calculateDiffMeth(meth,covariates = covariates,mc.cores = 4)

The output of getMethylDiff(myDiff, difference = 25, qvalue = 0.01)

is rather weird: I get p-values of 1 and 0.99, and the q-values are all 0.

When I do not use to co-variate argument, I get more " normal " p-value and q-values.

The batch "bisulfitelibs", which I use as covariates, come from the lab procedure. I split up all the samples and did bisulfite conversion 4 separate times, so I am thinking that this might have an impact on the data. Also, using the batch effect control feature in methylKit I found that the bisulfite processing might be relevant to add as a covariate.

Thanks a bunch for any suggestions of how to solve this.

Best, Line