Hi, I'm working with whole-genome EM-seq data from an experiment with two factors:

• Two genotypes (with contrasting response to a pathogen)
• Infected and non-infected

Only one tissue and one time analyzed, with two biological replicates for each treatment combination.

I would like to know:

• Which sites/regions are differentially methylated in response to infection independently of genotype. (Figure 1)
• Which sites/regions are diff. methylated in response to infection and between genotypes (Figure 2)
• Which sites/regions are diff. methylated in response to infection in a genotype-dependent manner. (Figure 3)
• Which sites/regions are differentially methylated between genotypes independently of infection. (Figure 4)

I followed the Bismark alignment and extraction pipeline, and now I'm following the methylKit package tutorials.

According to it, in cases like mine with two factors and two biological replicates, I can apply logistic regressions to test for differential methylation. In the calculateDiffMeth fuction, one of the factors should be "treatment" and the other should be "covariate". I assume that to get sites diff. methylated in response to infection (goals 1-3) I need to load the infection factor as "treatment" and the genotypes as "covariate".

BUT, the resulting object is simply a list of sites/regions and their pvalue, qvalue and methdiff, without specifying the significance of each component in the regression. Which means that, for example, I can't know whether the interaction between factors was significant (case Figure 3).

My question is,

• How do I access the full information of the tested logistic regressions?
• And if I can't, how do I work around it?

I understand that to get the sites diff. methylated due to genotype (goal 4) I should switch genotype as "treatment" and infection as "covariate". But would this somehow alter the correction I should apply to the p-values?

And again, how do I identify sites diff. methylated in response to infection depending of genotype?

I thought about testing diff. meth. against infection separately for each genotype and then contrast the resulting lists, but once again, would this somehow alter the correction I must apply to the p-values?

Here's an example of my code and my results:

> meth.CHG_promoters
methylBase object with 61192 rows
--------------
       chr start  end strand coverage1 numCs1 numTs1 coverage2 numCs2 numTs2 coverage3 numCs3 numTs3 coverage4 numCs4 numTs4 coverage5 numCs5 numTs5 coverage6 numCs6 numTs6 coverage7 numCs7 numTs7 coverage8 numCs8 numTs8
1 HanXRQCP    97 1597      -     41775    936  40839     40096   1047  39049     49664   1041  48623     40612    628  39984     69526   3333  66193     75164   1084  74080     53468   1418  52050     53950   2110  51840
2 HanXRQCP  1266 2766      +     46631    866  45765     42949   1025  41924     54426   1054  53372     45499    669  44830     75929   3657  72272     86452   1419  85033     60449   1657  58792     59809   2511  57298
3 HanXRQCP  1524 3024      -     29216    419  28797     25039    497  24542     33037    540  32497     28064    374  27690     44816   2021  42795     53065    812  52253     37224    934  36290     35293   1300  33993
4 HanXRQCP  1876 3376      -     29126    420  28706     24639    465  24174     32690    481  32209     27573    354  27219     43789   1889  41900     53152    783  52369     37005    883  36122     34954   1186  33768
5 HanXRQCP  4681 6181      -     34258    838  33420     28101    716  27385     37988    756  37232     33004    552  32452     49461   1798  47663     64667    700  63967     45144    883  44261     42012   1227  40785
6 HanXRQCP  5495 6995      -     43559    979  42580     37327    974  36353     49203   1055  48148     42015    782  41233     64553   3094  61459     84218   1236  82982     57475   1503  55972     56022   2263  53759
--------------
sample.ids: infected_var1_r1 infected_var1_r2 control_var1_r1 control_var1_r2 infected_var2_r1 infected_var2_r2 control_var2_r1 control_var2_r2 
destranded FALSE 
assembly: HanXRQ-2.0 
context: CHG 
treament: 1 1 0 0 1 1 0 0 
resolution: region 

> myDiff.CHG_promoters <- calculateDiffMeth(meth.CHG_promoters,
+ covariates=data.frame(genotype=c('var1','var1','var1','var1','var2','var2','var2','var2')),
+ overdispersion = "MN",adjust="BH",mc.cores=30)

> myDiff.CHG_promoters
methylDiff object with 61192 rows
--------------
       chr start  end strand    pvalue qvalue    meth.diff
1 HanXRQCP    97 1597      - 0.9161677      1  0.196036107
2 HanXRQCP  1266 2766      + 0.9603097      1  0.089608732
3 HanXRQCP  1524 3024      - 0.9686611      1  0.108272516
4 HanXRQCP  1876 3376      - 0.9622652      1  0.163918331
5 HanXRQCP  4681 6181      - 0.8713753      1  0.134653151
6 HanXRQCP  5495 6995      - 0.9290056      1 -0.001156726
--------------
sample.ids:  infected_var1_r1 infected_var1_r2 control_var1_r1 control_var1_r2 infected_var2_r1 infected_var2_r2 control_var2_r1 control_var2_r2  
destranded FALSE 
assembly: HanXRQ-2.0 
context: CHG 
treament: 1 1 0 0 1 1 0 0 
resolution: region

PS: I understand the three classic cytosine contexts (CpG, CHG and CHH) are analyzed independently from each other.