EdgeR design-matrix design

[extended.wobble]

Dear List,

we are working on an RNA-Seq experiment including 3 different cultivars of the same plant species during a drought treatment.
The plants have been grown in two different growth chambers. When analysing the data we had to find out that the "effect" of the chambers on gene expression is actually stronger than the drought treatment effect. Hence, we want to remove this (technical) chamber effect from the analysis in order to extract only those genes that are differentially expressed in response to the drought treatment.

Unfortunately we are a bit stuck with correctly formulating the design matrix needed to accomplish this analysis using edgeR. We created a design matrix without the chamber factor according to the manual like this:

genotypes <- factor( c("cultivar1", "cultivar1", "cultivar1", "cultivar1", "cultivar2", "cultivar2", "cultivar2", "cultivar2") )
groups <- factor( c("drought", "drought", "watered", "watered", "drought", "drought", "watered", "watered") )
design <- model.matrix(~genotypes + groups)

view(design)
	(Intercept)	genotypes	groups
1	1		1		0
2	1		1		0
3	1		1		1
4	1		1		1
5	1		0		0
6	1		0		0
7	1		0		1
8	1		0		1

but the situation we are facing looks like this:

1	cultivar1 drought chamber1
2	cultivar1 drought chamber2
3	cultivar1 watered chamber1
4	cultivar1 watered chamber2
5	cultivar2 drought chamber1
6	cultivar2 drought chamber2
7	cultivar2 watered chamber1
8	cultivar2 watered chamber2

How do we correctly set up the design matrix to eliminate the chamber effect and is there a way to explicitly define the contrasts we want to extract (analogous to the contrast matrix function in limma?)

thanks for any help

[Gordon Smyth]

Just add chambers to the design matrix

You can adjust for possible chamber effects by including the factor in the design matrix:

design <- model.matrix(~genotypes+groups+chambers)
y <- estimateGLMCommonDisp(y, design)
fit <- glmFit(y,design)
results.groups <- glmLRT(y,fit,coef=2)
topTags(results.groups)

This gives a list of transcripts different between watered vs drought, adjusting for cultivar and chamber differences.

To get the watered vs drought effects for cultivar1 and cultivar2 separately, you can use

design <- model.matrix(~genotypes+genotypes:groups+chambers)

then

results <- glmLRT(y,fit,coef=3)

for cultivar1 and

results <- glmLRT(y,fit,coef=4)

for cultivar2.

Best wishes
Gordon

PS. I answer questions much more often on the Bioconductor list than on seqanswers.

[extended.wobble]

Thanks a lot for your help.

Is there a documentation or manual that would have contained this information?

cheers.

PS: I didn't want to post the same topic twice, but started here at seqanswers.

[extended.wobble]

Thanks again Gordon

your advice did already make things much clearer for us. One interesting
question we would like to answer is which genes are differentially expressed
between the two genotypes in response to drought treatment, as these would
be the genes that might be responsible for different levels of tolerance.

To extract these we put together the following design matrix, assuming that
coefficient 5 would contain these genes. Is that correct?

genotypes <- factor( c("cultivar1", "cultivar1", "cultivar1", "cultivar1", "cultivar2", "cultivar2", "cultivar2", "cultivar2") )
chamber <- factor( c("K1", "K3", "K1", "K3", "K1", "K3", "K1", "K3") )
groups <- factor( c("Drought", "Drought", "Watered", "Watered", "Drought", "Drought", "Watered", "Watered") )

design <- model.matrix(~genotypes+groups+chamber+genotypes:groups)

 
(Intercept) genotypesCultivar2 groupsWa chamberK3 genotypesCultivar2:groupsWa
Dr           1            0        0         0                     0
Dr           1            0        0         1                     0
Wa           1            0        1         0                     0
Wa           1            0        1         1                     0
Dr           1            1        0         0                     0
Dr           1            1        0         1                     0
Wa           1            1        1         0                     1
Wa           1            1        1         1                     1

 
 
finally <- glmLRT(DGE_List, glmfit, coef=5)

should contain the significant changes between the two cultivars
regarding their reaction to drought, with chamber effects eliminated.

Thanks a lot in advance for your time