This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Issue when using design matrix taking batch effect into account - EdgeR RNA-seq

My groups are Diet : ND or WD and Treatment : CT or KO.

They have been extracted at date_1 (all ND) or date_2 (all WD).

I have done this design matrix :

design.matrix <- model.matrix( ~ 0 + diet + treatment + diet:treatment)

design.matrix
   dietND dietWD treatmentKO dietWD:treatmentKO
1       1      0           0                  0
2       1      0           0                  0
3       1      0           0                  0
4       1      0           0                  0
5       1      0           0                  0
6       1      0           0                  0
7       1      0           1                  0
8       1      0           1                  0
9       1      0           1                  0
10      1      0           1                  0
11      1      0           1                  0
12      1      0           1                  0
13      0      1           0                  0
14      0      1           0                  0
15      0      1           0                  0
16      0      1           0                  0
17      0      1           0                  0
18      0      1           1                  1
19      0      1           1                  1
20      0      1           1                  1
21      0      1           1                  1
22      0      1           1                  1
attr(,"assign")
[1] 1 1 2 3
attr(,"contrasts")
attr(,"contrasts")$diet
[1] "contr.treatment"

attr(,"contrasts")$treatment
[1] "contr.treatment"

Trying to see without the batch effect by including the extract dates in a new design matrix :

table(dgef$samples$group, block_extract)
        block_extract
                 date_1          date_2
  ND_CT               6               0
  ND_KO               6               0
  WD_CT               0               5
  WD_KO               0               5

designBlock <- model.matrix( ~ 0 + diet + treatment + diet:treatment + block_extract)
designBlock
 dietND dietWD treatmentKO block_extractdate_18/02/2025 dietWD:treatmentKO
1       1      0           0                            0                  0
2       1      0           0                            0                  0
3       1      0           0                            0                  0
4       1      0           0                            0                  0
5       1      0           0                            0                  0
6       1      0           0                            0                  0
7       1      0           1                            0                  0
8       1      0           1                            0                  0
9       1      0           1                            0                  0
10      1      0           1                            0                  0
11      1      0           1                            0                  0
12      1      0           1                            0                  0
13      0      1           0                            1                  0
14      0      1           0                            1                  0
15      0      1           0                            1                  0
16      0      1           0                            1                  0
17      0      1           0                            1                  0
18      0      1           1                            1                  1
19      0      1           1                            1                  1
20      0      1           1                            1                  1
21      0      1           1                            1                  1
22      0      1           1                            1                  1
attr(,"assign")
[1] 1 1 2 3 4
attr(,"contrasts")
attr(,"contrasts")$diet
[1] "contr.treatment"

attr(,"contrasts")$treatment
[1] "contr.treatment"

attr(,"contrasts")$block_extract
[1] "contr.treatment"

But when I continue with the designBlock:

  logCPM_RMVexpe <- removeBatchEffect(
   cpm(
     dgef,
     log = T,
     prior = 1,
     norm = T
   ),
   batch = block_extract,
   design = designBlock,
   group = dgef$samples$group
 )

I have this message :

Coefficients not estimable: batch1.
Partial NA coefficients for 16443 probe(s)

When I try to estimate the dispersion :

dgef_block$design.matrix <- designBlock
dgef_block <- estimateDisp(dgef_block, designBlock, robust = T)

Error in glmFit.default(sely, design, offset = seloffset, dispersion = 0.05,  : 
Design matrix not of full rank.  The following coefficients not estimable:
block_extractdate_1

I don't really get why can't I use designBlock nor what the error means. I am not sure either if I can continue without estimating the dispersion, as I am following an existing code. Is it linked to the fact that each diet is batched separately?

rna-seq edger design-matrix
• 271 views • link
updated by Ram Ram • written by maxime • maxime

1 answer

You bock is exactly the same thing as your treatment - all your dietND samples are in one block and all your dietWD are in the other. Thus, there is no way for edgeR to know whether any difference between these samples is because of the diet or the extraction date.

This is called confounding. We say that your diet effect is completely confounded by extraction date effect. You can estimate that effect, but you can't mathemetically assign it to either diet or extraction. The only thing you can really do here is decide whether its more biologically plausible that any difference is due to diet or extraction day.

In future its worth considering such effects when designing expriments - for example, rather than extracting ND on one day and WD on the other, extract replicates 1-3 of ND and 1-3 of WD on day 1 and replicates 4-6 of ND and 4-5 of WD on day 2. But you'd also have to be careful to balance your KO and WT across the two days as well.

The art of doing this in such a way as to allow effective estimation of the batch effect is called "blocking", and is part of formal exprimanetal design theory.

• 0 views • link
written by i.sudbery i.sudbery

Log in to answer this question.

More posts like this