[BioC] edgeR: GLM for multi-factor and mulit-level designs

[Dorota Herman]

thanks a lot for your respond

best wishes
Dorota


Ryan C. Thompson wrote:
> I think you have both tests correct. Each coefficient in the model 
> matrix from conds is an "X vs A" coefficient. So the "condsB" 
> coefficient is the coefficient for B vs A. That means that you can 
> test B vs A by simply testing that coefficient equal to zero, as you 
> have done. Note that you could also make a contrast for the same test 
> like this:
>
> B_A<-makeContrasts(condsB,levels=design)
>
> Essentially, you are testing the null hypothesis that the specified 
> coefficient is zero for each gene.
>
> which would make a vector with just one nonzero element. When you test 
> the contrast of condsC-condsD, you are subtracting "D vs A" from "C vs 
> A", which is the same as adding "C vs A" and "A vs D". Thus this 
> difference gives "C vs D", which is what you want. Condition A is not 
> being ignored; it simply cancels out of the contrast in question. 
> Instead of testing whether a single coefficient is equal to zero, you 
> are testing whether the difference of two coefficients is zero.
>
> By the way, note the order: "condsC - condsD" will give the same 
> statistics as "condsD - condsC", but will give opposite fold changes, 
> so choose the one that makes the most sense to you (typically 
> "treatment - control" for simple experiments).
>
>
> On Tue 27 Nov 2012 07:18:20 AM PST, Dorota Herman wrote:
>>
>> Hello everyone,
>>
>> I have been playing with the GLM approach for RNA-seq data in DESeq
>> and edgeR but I am fairly new in DE analyses. I am interested in
>> pairwise comparisons in multi-factor multi-level designs. My question
>> concerns my understanding of an application of the glmLRT function
>>
>> #My code is
>>>
>>> countsTable <- read.delim(file)
>>> header <-
>>
>> c('A_1','A_2','A_3','B_1','B_2','B_3','C_1','C_2','C_3','D_1','D_2','D_3') 
>>
>>
>>>
>>> names(countsTable) <- header
>>> conds <- factor(c('A','A','A’,'B','B','B','C','C','C','D','D','D'))
>>> Ex<-factor(c('exper1', 'exper2', 'exper3', 'exper2', 'exper3', 
>>> 'exper4',
>>
>> 'exper1', 'exper3', 'exper4', 'exper2', 'exper3', 'exper4'))
>>>
>>> group <- conds
>>> dge <- DGEList(counts=countsTable,group=group)
>>> dge <- calcNormFactors(dge, method='TMM')
>>> design <- model.matrix(~Ex+conds)
>>> rownames(design)<-colnames(dge)
>>> dge <- estimateGLMCommonDisp(dge,design)
>>> dge <- estimateGLMTrendedDisp(dge, design)
>>> dge <- estimateGLMTagwiseDisp(dge, design)
>>> fit <- glmFit(dge, design)
>>
>>
>> #my design looks like:
>>>
>>> design
>>
>> (Intercept) Eexper2 Eexper3 Eexper4 condsB condsC condsD
>> A_1 1 0 0 0 0 0 0
>> A_2 1 1 0 0 0 0 0
>> A_3 1 0 1 0 0 0 0
>> B_1 1 1 0 0 1 0 0
>> B_2 1 0 1 0 1 0 0
>> B_3 1 0 0 1 1 0 0
>> C_1 1 0 0 0 0 1 0
>> C_2 1 0 1 0 0 1 0
>> C_3 1 0 0 1 0 1 0
>> D_1 1 1 0 0 0 0 1
>> D_2 1 0 1 0 0 0 1
>> D_3 1 0 0 1 0 0 1
>> attr(,"assign")
>> [1] 0 1 1 1 2 2 2
>> attr(,"contrasts")
>> attr(,"contrasts")$E
>> [1] "contr.treatment"
>> attr(,"contrasts")$conds
>> [1] "contr.treatment"
>>
>> I understand that R function rewrote the model matrix because of the
>> identifiability problem for parameter estimations. However it causes
>> my confusion in further usage of that design for the pairwise
>> comparisons.
>>
>> In a case when I want to obtain differentially expressed genes between
>> A and B, I understand I should use the function:
>>>
>>> lrt <- glmLRT(fit,coef="condsB")
>>
>> Is it correct?
>>
>> In a case when I want to obtain differentially expressed genes between
>> C and D (*without taking into account A*), are these calling functions
>> correct?
>>>
>>> C_D<-makeContrasts(condsC-condsD,levels=design)
>>> lrt <- glmLRT(fit,contrast=C_D)
>>
>>
>> Does it mean that glmLRT function takes into account first conds (A)
>> when we use ‘coef’ parameter and discard it when we use ‘contrast’
>> parameter? Or it means that the second analysis, between C and D takes
>> into account differential expression with A too?
>>
>> I hope my explanation of the question is not too confusing.
>> Best wishes
>> Dorota


-- 
==================================================================
Dorota Herman, PhD 
VIB Department of Plant Systems Biology, Ghent University
Technologiepark 927
9052 Gent, Belgium
Tel: +32 (0)9 3313692
Email:dorota.herman at psb.vib-ugent.be
Web: http://www.psb.ugent.be