RUVseq design matrix for first pass of edger
1
0
Entering edit mode
tonja.r ▴ 80
@tonjar-7565
Last seen 7.4 years ago
United Kingdom

RUVg uses the first call of edgeR to find the factors of unwanted variance. The design matrix is constructed with the possible batch effect variable:

 x <- as.factor(rep(c("Ctl", "Trt"), each=2)) 
  design = model.matrix(~x, data=pData(set))
  y <- DGEList(counts=counts(set), group=x)


If I know that I have a batch effect between first replicates (Ctl1, Trt1) and second replicates (Clt2, Trt2). Would it make sense to introduce it already in the design matrix and perform a first call of edgeR and  then RUVg?

ruvseq • 1.7k views
ADD COMMENT
1
Entering edit mode
davide risso ▴ 950
@davide-risso-5075
Last seen 3 days ago
University of Padova

Hi Tonja,

given your example, I understand that you only have two replicates per class: two treated and two controls. If that's the case, then you cannot include both the known batch effect and a factor of RUV because you don't have enough degrees of freedom in the model.

In fact, what do you mean here by batch? The design that you are describing seems a paired design to me, i.e., you have two _pairs_ of treated, control experiments. You may want to take advantage of that in your analysis, by specifying a paired design in edgeR (see chapter 3.4 of the edgeR User's guide).

A minimal example would be:

expt <- factor(c(1, 1, 2, 2))
treatment <- factor(c("t", "c", "t", "c"))
design <- model.matrix(~expt+treatment)

Note that if you do that, you need to estimate three regression coefficients: for the intercept, for experiment (batch), and for treatment. Then you need to estimate the variance as well. Hence, you don't have any degrees of freedom left (4 observations, 4 parameters). And it is not possible to add any additional factor in the design, i.e., you cannot use RUV in this case.

Obviously, you could drop the expt variable, and use one factor from RUV instead, but it seems to be a better strategy in this case to make use of the information that the samples are paired.

 

ADD COMMENT
0
Entering edit mode

I guess I did not express myself well enough. The minimal example for edgeR analysis is clear for me and you got it right: Ctl1 and Trt1 were processed together, and Ctl2 and Trt2 together: so the bach effect would be between the first pair of treated/control and the second pair.


 If that's the case, then you cannot include both the known batch effect and a factor of RUV because you don't have enough degrees of freedom in the model.
I do not want to do this. In RUVg one performs the first pass (step A) of edgeR analysis to identify the most non-differential genes. The design matrix for this first pass is (as suggested in RUVseq manual)

x <- as.factor(rep(c("Ctl", "Trt"), each=2))
design = model.matrix(~x, data=pData(set))

one perform first pass of edgeR analysis , then the most non-differential genes are introduced to RUVg ​then one obtains the factos from RUVseq and inputs them into a new design matrix for the downstream analysis with edgeR (step B).
My question was if it makes sense to modify design matrix to account for batch effect already in step A: in the first pass of edgeR in RUVseq to identify the most non-differential genes.

Meaning, following:
step A:

expt <- factor(c(1, 1, 2, 2))
treatment <- factor(c("t", "c", "t", "c"))
design <- model.matrix(~expt+treatment)

Perform DE analysis with edgeR to find most non differential genes and input them into RUVg. Get factors (W_1) from RUVg and now the design matrix would look like following for the downstream analysis with edgeR (step B):

design <- model.matrix(~W_1+treatment)

 


 

ADD REPLY
1
Entering edit mode

The answer to your question depends on the design of the experiment and on what you want the factors of unwanted variation to capture.

*In general* (i.e., with an enough number of observations) you could have a step A that looks like this:

design <- model.matrix(~expt+treatment)

and a step B that looks like this:

design <- model.matrix(~expt+treatment+W)

This would mean that W captures _additional_ variation that is not captured by batch and treatment. With only two observations per class you cannot fit this model, hence you have to remove either W or expt. 

In this case, I'm not sure it would help  to use expt in the design of step A, because you want W to capture the batch effect as well. So I would just use the simpler step A:

design <- model.matrix(~treatment)
ADD REPLY

Login before adding your answer.

Traffic: 632 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6