This is my first question here. So sorry If I am not that clear:

I am using edgeR/limma to compare chromatin accessibility for three groups in two different condition (treated or not). Since the study has a matched design, I am using a nested model to account for inter-individual variability within the groups:

design2 = model.matrix(~ Block + Condition:Treatment)
dge <- DGEList(Matrix)
dge <- calcNormFactors(dge)
v <- voom(dge, design2)
fit <-lmFit(v, design2)
fit <- eBayes(fit)
res_fit <- decideTests(fit)
summary.TestResults(res_fit)

GroupA.Treatment GroupB.Treatment GroupC.Treatment
Down 1164 0 0
NotSig 35989 39917 39917
Up 2764 0 0

desgin2:

X.Intercept Block2 Block3 Block4…. GroupA.Treatment GroupB.Treatment GroupC.Treatment
1 1 0 0 1 0 0
1 1 0 0 0 0 0
1 0 1 0 1 0 0
1 0 1 0 0 0 0
1 0 0 1 0 1 0
1 0 0 1 0 0 0
1 0 0 0 0 1 0
1 0 0 0 0 0 0
1 0 0 1 0 0 1
1 0 0 1 0 0 0
1 0 0 0 0 0 1
1 0 0 0 0 0 0

I don’t know if GroupB and GroupC don’t respond to treatment or if due to the group condition the chromatin was already altered at the baseline. So, I tested that with the model:

design3 = model.matrix(~ Condition + Treatment, data=pheno)
v2 <- voom(dge, design3)
fit1 <-lmFit(v2, design3)
fit1 <- eBayes(fit1)
res_fit1 <- decideTests(fit1)
summary.TestResults(res_fit1)

desgin3:

X.Intercept GroupB GroupC Treatment
1 1 0 0
1 1 0 1
1 1 0 0
1 1 0 1
1 0 1 0
1 0 1 1
1 0 1 0
1 0 1 1

The inter individual variability seems higher than the treatment itself. Is there a way to account for inter individual variability at the baseline with the last model? I cannot add “Block” as a covariate due to collinearity?

Hi Vini, welcome to Bioconductor.

I have edited your question to format the code and R output better. You would find it very rewarding to become more familiar with markdown, which is used as the editor on this Support forum: https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet

It is hard to answer your question because the output you show doesn't match the R code you claim to have run. For example, the output shows summary(res_fit) to have three columns (GroupA.Treatment, GroupB.Treatment and GroupC.Treatment), but that is impossible, first because design2 has more than three columns and second because there is no Group variable in the call to model.matrix(). I suspect that wherever you have mentioned Condition in your post you really mean Group.

When you ask a question here, it is very important that you run all the code at once in a new R session so that you know the output matches the code. Copying disconnected code and output is a big no-no.

Another strange issue is that design2 and design3 have different numbers of rows whereas, according to your post, they should be the same.

Two more minor comments. You never need to use summary.TestResults(). Just use summary() and the correct method will be called automatically. Also, the correct use of voom assumes that you have filtered out genes that are not expressed at worthwhile levels. I don't see any filtering step in your code.

If you want any more help, you would need to show us the targets frame for your experiment, i.e., the three factors Block, Group and Treatment associated with your experiment. Just print the data.frame containing those variables. We don't need to see design matrices. You also need to correct the code and output so that it is consistent and truly reflects your experiment. At the moment, we can't see what the problem is. Why should it there be any problem including Block in the second design matrix given that you included it in the first?