Comparing conditions were altered when including adjustCovariate for some covariates

[MinhTTNguyen]

Hello,

I'm reaching out for help with my EM-Seq data analysis using dmrseq to identify DMRs between two puberty epoch groups (Late vs. Early) within the same sex (M: male). I have three samples per group. Since my dataset is large (human genome), I created a bsseq object for each of the 24 human chromosomes and ran dmrseq on each unsmoothed bsseq object using the following command:

regions <- dmrseq(bs=dat, 
                               cutoff=0.05, 
                               maxPerms=10,
                               testCovariate="group",
                               adjustCovariate=c("pH","Bisulfite_conversion"))

However, I noticed an issue in the output displayed in RStudio:

Assuming the test covariate group is a factor.
Condition: 3_LateM1_EarlyM vs. 1_EarlyM3_LateM

The expected condition should be "3_LateM vs. 1_EarlyM", but when I include the adjustCovariate argument, dmrseq() alters the condition formatting unexpectedly.

Interestingly, when I exclude adjustCovariate=c("pH", "Bisulfite_conversion"), the function runs correctly. Moreover, when adjusting for some other covariates, dmrseq() behaves as expected, but for certain covariates, the "condition" variable gets altered.

I assigned "group" as factor in the pData().

Since I don't see any issues with the input data or metadata, I'm uncertain why dmrseq() modifies the "condition" variable when adjusting for certain covariates. Could you please help me with some explanation about this issue?

Thank you very much in advance!

[kdkorthauer]

Hi @MinhTTNguyen,

Thanks for reaching out! If I understand correctly, you have 6 samples total, evenly divided into two groups (early vs late puberty epoch). All samples are Male, and you are interested in comparing the two groups, adjusting for covariates "pH" and "Bisulfite_conversion". So that I can figure out what is causing this issue, it would help to know more about the adjustment covariates. Could you simply print out the relevant columns of pData for me (group, pH, and Bisulfite_conversion) for all six rows? Thanks!

Best,
Keegan

[MinhTTNguyen]

Hi Keegan,

Thanks so much for your response. Yes, I have 6 samples from 2 groups (Early vs Late puberty epoch), and 3 samples per group. All of the samples are Male. I want to adjust the 2 covariates "pH", and "Bisulfite_conversion" when running dmrseq() function.

Here is my metadata table:
> pData(dat)
DataFrame with 6 rows and 11 columns
group Bisulfite_conversion pmi pH RIN CAUSE_OF_DEATH REPOSITORY Substances hL2.3 excitatory inhibitory

X5387 1_EarlyM 86.2 13.0 6.31 7.5 Drowning Maryland 0 0.109958 0.589455 0.300587
S.469 3_LateM 97.9 37.5 6.56 0.0 Suicide Douglas 0 0.300633 0.473090 0.226276
X1908 1_EarlyM 99.4 13.0 NA 0.0 Hanging Maryland 0 0.225814 0.537684 0.236502
S.139 3_LateM 99.4 35.0 6.59 0.0 Accident Douglas 0 0.266940 0.505645 0.227415
S.135 3_LateM 99.0 68.5 5.94 0.0 Natural Douglas Nil 0.139144 0.544718 0.316138
X1790 1_EarlyM 99.4 18.0 5.54 8.2 Multiple Injuries Maryland 0 0.206619 0.523388 0.269993

To give you more information, "dat" is the bsseq object that stores the methylation data of the 6 samples. Here is more information about "dat":
> dat
An object of type 'BSseq' with
538846 methylation loci
6 samples
has not been smoothed
All assays are in-memory

Also, here is my code when running dmrseq:

regions <- dmrseq(
  bs=dat,
  cutoff = 0.05,
  maxPerms=2, # I will change maxPerms to 10 when running the analysis, "maxPerms=2" is only for testing 
  testCovariate="group",
  adjustCovariate=c("pH","Bisulfite_conversion"))

This is what was printed out on the screen when running dmrseq:

Assuming the test covariate group is a factor.
Condition: 3_LateM1_EarlyM vs 1_EarlyM3_LateM
Adjusting for covariate (s): Bisulfite_conversion, pH
Parallelizing using 62 workers/cores (backend: BiocParallel:MulticoreParam).
Computing on 1 chromosome(s) at a time.

Detecting candidate regions with coefficient larger than 0.05 in magnitude.
...Chromosome chr22: Smoothed (0.24 min). Error in (function (cond) :
error in evaluating the argument 'args' in selecting a method for function 'do.call': BiocParallel errors
62 remote errors, element index: 1, 68, 135, 202, 269, 336, ...
4062 unevaluated and other errors
first remote error:
Error in data.frame(g.fac = rep(pDat[, colnames(design)[coeff[1]]], each = length(ix)), : arguments imply differing number of rows: 30, 25

The function seems to work well when I don't include adjustCovariate=c("pH","Bisulfite_conversion")

> regions <- dmrseq(
+   bs=dat,
+   cutoff = 0.05, #Matt: cutoff = 0.02
+   maxPerms=2, 
+   testCovariate="group")

Assuming the test covariate group is a factor.
Condition: 3_LateM vs 1_EarlyM
Parallelizing using 62 workers/cores (backend: BiocParallel:MulticoreParam).
Computing on 1 chromosome(s) at a time.

Detecting candidate regions with coefficient larger than 0.05 in magnitude.
...Chromosome chr22: Smoothed (0.3 min). 4215 regions scored (0.41 min).

• 4215 candidates detected
  Performing balanced permutations of condition across samples to generate a null distribution of region test statistics

Beginning permutation 1
...Chromosome chr22: Smoothed (0.24 min). 3667 regions scored (0.46 min).

• 1 out of 2 permutations completed (3667 null candidates)

Beginning permutation 2
...Chromosome chr22: Smoothed (0.24 min). 3447 regions scored (0.46 min).

• 2 out of 2 permutations completed (3447 null candidates)
  There were 50 or more warnings (use warnings() to see the first 50)

Your help is very much appreciated!

Regards,

Minh

[kdkorthauer]

Hi Minh,

As a quick check, can you please see if you get the same error if you adjust only for Bisulfite_conversion? I noticed that the pH is missing (NA) for one of the samples, so this could be related to the issue (but I need to dig further)... You won't be able to adjust for a covariate that has missing values, unfortunately.

[MinhTTNguyen]

Hi Keegan,

Thanks for your comment. I excluded the "pH" covariate in the adjustCovariate argument and it seemed that it worked!!! (Although I did not find any DMRs passing qvalue<0.05. But I should expand my dmrseq analysis to other chromosomes.)

> regions <- dmrseq(
+   bs=dat,
+   cutoff = 0.05, 
+   maxPerms=2, 
+   testCovariate="group",
+   adjustCovariate=c("Bisulfite_conversion"))

Assuming the test covariate group is a factor.
Condition: 3_LateM vs 1_EarlyM
Adjusting for covariate (s): Bisulfite_conversion
Parallelizing using 62 workers/cores (backend: BiocParallel:MulticoreParam).
Computing on 1 chromosome(s) at a time.

Detecting candidate regions with coefficient larger than 0.05 in magnitude.
...Chromosome chr22: Smoothed (0.24 min). 4215 regions scored (0.51 min).

• 4215 candidates detected
  Performing balanced permutations of condition across samples to generate a null distribution of region test statistics

Beginning permutation 1
...Chromosome chr22: Smoothed (0.23 min). 3667 regions scored (0.46 min).

• 1 out of 2 permutations completed (3667 null candidates)

Beginning permutation 2
...Chromosome chr22: Smoothed (0.25 min). 3447 regions scored (0.45 min).

• 2 out of 2 permutations completed (3447 null candidates)

Should we have any warning or error message showing that the metadata does not have valid values to process the dmrseq function instead of having the altered conditions? I am just thinking that it may be useful for users to identify where the issue comes from :-).

Again, thanks very much for your help! I really appreciate it!

[kdkorthauer]

Thank you for checking that. Yes, definitely - thank you for the suggestion! When I get a chance, I will add a flag to check for NAs in covariates that will throw a more informative error message.

Thanks again for reporting this issue!

[MinhTTNguyen]

Thanks very much Keegan for helping me solve the issue!

Minh