Hi I'd like to use DMRcate to identify DMRs between treatment and control samples at two different time points and then compare the DMRs at these two time points to see how many of them are common and unique. To see whether DMR at time point 1 is DMR or not at time point 2, I first need to figure out whether the corresponding region is analyzable or not at time point 2. How should I do this in DMRcate?

Based on the manual, the data used has biological replicates, I'm wondering whether DMRcate still works well if there is no biological replicate available. Thanks.

Hi xie186,

Yes, you absolutely can do this by hypothesising an interaction effect between treatment and timepoint. So say you have a minimally fully blocked experiment with factors treat and timepoint, with 2x2x2 design, then some boilerplate code might look like this:

> treat <- c(rep("Treatment", each=4), rep("Control", each=4))
> treat
[1] "Treatment" "Treatment" "Treatment" "Treatment" "Control"   "Control"   "Control"   "Control"  
> timepoint <- rep(c("T1", "T2"), times=4)
> timepoint
[1] "T1" "T2" "T1" "T2" "T1" "T2" "T1" "T2"

Then you would design your model matrix like so:

> design <- model.matrix(~treat*timepoint)
> design
  (Intercept) treatTreatment timepointT2 treatTreatment:timepointT2
1           1              1           0                          0
2           1              1           1                          1
3           1              1           0                          0
4           1              1           1                          1
5           1              0           0                          0
6           1              0           1                          0
7           1              0           0                          0
8           1              0           1                          0
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$treat
[1] "contr.treatment"

attr(,"contrasts")$timepoint
[1] "contr.treatment"

And then transform to a methModelMatrix using edgeR:

> methdesign <- edgeR::modelMatrixMeth(design)
> methdesign
   Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 Sample7 Sample8 (Intercept) treatTreatment timepointT2 treatTreatment:timepointT2
1        1       0       0       0       0       0       0       0           1              1           0                          0
2        1       0       0       0       0       0       0       0           0              0           0                          0
3        0       1       0       0       0       0       0       0           1              1           1                          1
4        0       1       0       0       0       0       0       0           0              0           0                          0
5        0       0       1       0       0       0       0       0           1              1           0                          0
6        0       0       1       0       0       0       0       0           0              0           0                          0
7        0       0       0       1       0       0       0       0           1              1           1                          1
8        0       0       0       1       0       0       0       0           0              0           0                          0
9        0       0       0       0       1       0       0       0           1              0           0                          0
10       0       0       0       0       1       0       0       0           0              0           0                          0
11       0       0       0       0       0       1       0       0           1              0           1                          0
12       0       0       0       0       0       1       0       0           0              0           0                          0
13       0       0       0       0       0       0       1       0           1              0           0                          0
14       0       0       0       0       0       0       1       0           0              0           0                          0
15       0       0       0       0       0       0       0       1           1              0           1                          0
16       0       0       0       0       0       0       0       1           0              0           0                          0

Given your bsseq object, you would call sequencing.annotate() like this, using coef=12 because the 12th column of methdesign is the interaction effect you are interested in:

seq_annot <- sequencing.annotate(my_bsseq, methdesign, all.cov = TRUE, coef = 12, fdr=0.05)

As for replicates, yes they are necessary. Above is a minimal example for a second-order effect model, which means you need minimum 8 WGBS samples designed in the structure described above. If you don't have at least n=2 for every subgroup of treatment x timepoint, your experiment will be incompletely blocked, and you will need to remove an entire level from treat or timepoint to have your coefficients estimated correctly.

Best, Tim