I want to run some epigenome wide association studies based on city and pollutant type, could you please suggest if I can get DNA methylation sites (CpG) association based on pollutant type for each city. like below I am trying to find Pb toxic elements association site in four different cities which has varying level of Pb. we do have 10 toxic element so not sure should I run model with one toxic element at a time or should include all in this model?

#model matrix

var<-model.matrix(~log2Pb*City + sex + Age +CD8T +CD4T +NK + Bcell +Mono, data=data1)

fit<-lmFit(mval,var)


fit2<-eBayes(fit,trend=TRUE, robust=TRUE)

probe<-topTable(fit2,adjust="BH",coef=???,num=Inf)

try to run interaction model using limma but not sure if this is the right approach and what to mention in coef????

Thanks for your all help!

The interaction term should be the last one in your design matrix, so

topTable(fit2, ncol(var), Inf, p.value = 0.05)

Will give you the CpGs with an FDR<0.05 for the interaction term.

There are many ways to specify the interaction coefficients, but maybe the easiest here would be:

topTable(fit2, coef=grep(":",colnames(var)), num=Inf)

To see how this works, just type colnames(var) or colnames(fit2). You will easily see the interaction coefficients. Then you specify them to topTable either by name or by number.