Hi,

I have some 450k data for ~1000 samples. I want to use ComBat to adjust for 2 known batches. I would like to adjust for Age_Group and Sentrix_ID (cf. example below). My variable of interest is Sample_Group. I don't really know how to proceed (run ComBat twice? which is yes accordingly to old posts) and what to put in the model.matrix. I was thinking to do:

First:

batch <- pheno$Age_Group
mod <- model.matrix(~as.factor(Sample_Group) + as.factor(Sentrix_ID) + as.factor(Sentrix_Position), pheno)

And then:

batch <- pheno$Sentrix_ID
mod <- model.matrix(~as.factor(Sample_Group) + as.factor(Age_Group) + as.factor(Sentrix_Position), pheno)

>data_table()

>pheno_table()

Thank you in advance for your assistance.

Hi, 

Its actually unclear to me what to recommend to you at this point. First off, I must say I have never tried a two-pass application of ComBat in this way, so I have no idea if it would work. Maybe it would be fine, but I don't know what potential problems you might run into.

Is your sample size large enough to allow you to combine the two variables? This would probably be best if you can do it (I see six samples above, but I'm assuming this is just a small sample of the complete dataset?)  

Also, are you sure that there are effects due to both of these covariates in your data? Sometimes what you think should be the batch variable isn't really important. Do they cluster by these variables? As PCA and side-by-side sample boxplots are also very informative. 

Finally, and probably most important, it seems to me that Age group is something you would like to include in a downstream model, not necessarily remove in the batch adjustment stage. Why do you want to remove it as opposed to just account for it at the differential expression stage?

Thanks!

Evan

Hi Evan,

Thanks for your feedbacks.

Is your sample size large enough to allow you to combine the two variables? This would probably be best if you can do it (I see six samples above, but I'm assuming this is just a small sample of the complete dataset?)

My sample size is ~1500 samples. By combining the 2 batch variables do you mean something like that:

Also, are you sure that there are effects due to both of these covariates in your data? Sometimes what you think should be the batch variable isn't really important. Do they cluster by these variables? As PCA and side-by-side sample boxplots are also very informative.

PCA plot shows clear effect related to Sentrix_ID.

I can't be sure for the Age_Group variable. However I'm working on a infantile disease, and a lot of my control samples are "old" patients (compare to my patients samples which are mostly young). It has been shown that methylation changes with age and in that regards I need to correct for this "potential" biological batch effect.

Finally, and probably most important, it seems to me that Age group is something you would like to include in a downstream model, not necessarily remove in the batch adjustment stage. Why do you want to remove it as opposed to just account for it at the differential expression stage?

Unfortunately I'm not really an expert in doing models :S. We are using the final adjusted data_table to compare methylation at targeted regions in our patients versus controls. In addition it is "easier" for us to have an adjusted data_table for meta analysis purposes.

Thanks!!

So according to this post C: ComBat: 3 adjustment variables & continuous adjustment variables I was thinking to do a 3 steps adjustment (I added Sentrix_Position as batch effect too). Does it look right (ie. model.matrix)?

#1 - correct for Sentrix_ID effect
bat.1<- ComBat(dat=datMval,
  batch=pheno$Sentrix_ID,
  mod=model.matrix(~as.factor(Sample_Group) + as.factor(Sentrix_Position) + as.factor(Age), data=pheno))

#2 - correct for Sentrix_Position effect
bat.2<- ComBat(dat=bat.1,
  batch=pheno$Sentrix_Position,
  mod=model.matrix(~as.factor(Sample_Group) + as.factor(Age), data=pheno))

#3 - correct for age effect
bat.3<- ComBat(dat=bat.2,
  batch=pheno$Age,
  mod=model.matrix(~as.factor(Sample_Group), data=pheno))

Thank you!