Dear bioconductor community,

I have recently analyzed two EPIC methylation array datasets using the same pipeline in R. Both of the datasets have the same variables and the comparison was PET Amyloid negative vs PET amyloid positive samples in order to find differentially methylated CpG sites. Both of the datasets are from EPIC array, but different batch(batch 0 and batch 1).

I've already tried to correct the batch effect using Harman batch correction

shifted_betas <- shiftBetas(betas=getBeta(mSetFunnormFlt), shiftBy=1e-4)
#mSetFunnormFlt is GenomicRatioSet which is the result of probe QC and Funnorm normalization
shifted_ms <- beta2m(shifted_betas)
plot(density(shifted_ms, 0.05), main="Shifted M-values, shiftBy = 1e-4",
     cex.main=0.7)
shifted_ms
methHarman <- harman(shifted_ms, expt=targets_pr$Sample_Group,
                     batch=targets_pr$batch, limit=0.65)
#Sample_Group was divided into P and N (PET amyloid positive and Negative)

ms_hm <- reconstructData(methHarman)

fit <- lmFit(ms_hm, design)

# design matrix was created by using the code below
# design <- model.matrix(~0+Sample_Group+Age+Sex+Center+Smoking+Bcell+CD4T+CD8T+Mono+NK+Neu, data=targets_pr_0)

After batch correction, I could not get any significant differentially methylated CpG site (PET Amyloid positive vs PET Amyloid negative)

When I did methylation analysis with same pipeline separately (batch0 and batch1), I found 2,756 significant DMCs from batch 0 but 0 significant DMCs from batch 1.

I thought that case-control unbalance may cause something wrong.

My question is whether it is possible and applicable in someway to correct the batch effect in unbalanced case-control dataset?

Thank you in advance.

Best,

Yujin Kim.

You should not apply batch correction functions like harman before a limma analysis. Instead, just add batch to the limma linear model.

The unbalanced batches do not cause any problem for limma except that the power to detect changes will be lower for batch1.

In this case, it sounds as if the data for batch1 behaves differently from batch0. Perhaps you should use the analysis of batch0 only?