Luis,

Your design has 4 groups,which means that the inter-group ANOVA-like 
test has 4 - 1 = 3 degrees of freedom. Hence, the ANOVA-like test is 
performed by testing 3 coefficients/contrasts equal to zero. The 
procedures I have described will yield will test where the null 
hypothesis is all group means being equal, and the alternative 
hypothesis is that any two or more group means differ significantly from 
each other. It is true that between 4 groups, there are 6 possible 
pairwise comparisons to be made. If you want individual significance 
rankings for each of the 6 possible contrasts, you can do that too, by 
making 6 separate calls to glmLRT, each with a single coefficient or 
contrast. Doing so will give you an appropriate rank order of 
significance for each contrast, but the exact FDR values will probably 
be too liberal, since there is no correction for the fact that you are 
performing 6 tests on each gene instead of just 1. Ideally, you would 
use a post-hoc testing method to assess individual contrasts within 
genes that are called significant by their FDR in the ANOVA-like test. 
However, I'm not sure how to perform such tests with the negative 
binomial distribution. You might consider using the limma-voom method 
instead of edgeR, since that will perform an actual ANOVA based on the 
normal distribution, so standard ANOVA post-hoc analysis is possible. If 
anyone has a take on how to do rigorous post-hoc testing on negative 
binomial GLMs, I'd be interested in hearing about it.

If you just want the logFC values for the missing 3 contrasts, then 
simply note that they can all be written in terms of the X vs 10 
contrasts. For example, the 20 vs 40 contrast would be the 10 vs 40 
contrast minus the 10 vs 20 contrast.

-Ryan

On 03/12/2013 06:42 PM, Luis alberto Martinez lopez wrote:
> thank you Ryan for your clear answer:
> but in your answer again only considered the comparison versus the 10 
> group(the first)
> i dont't have a control group so i'm looking for a way of doing an 
> anova like test in edgeR that compares 10 vs 20 ,20 vs 40, 40 vs 60,10 
> vs 40, 10 vs 60,20 vs 60 ie all the posibilities at the same time and 
> that give me a list with the genes more differential in any condition,
> thanks in advance,
> Luis
> *De:* Ryan C. Thompson <rct@thompsonclan.org>
> *Para:* Luis alberto Martinez lopez <lewisluis2013@yahoo.com>
> *CC:* "bioconductor@r-project.org" <bioconductor@r-project.org>
> *Enviado:* Lunes, 11 de marzo, 2013 22:46:19
> *Asunto:* Re: [BioC] anova like test in edgeR
>
> Luis,
>
> Note that your experimental design is the same whether or not you 
> include an intercept in the design matrix. If you include an 
> intercept, then the intercept represents the level in the D10 group, 
> and the other three coefficients represent the differences between 
> D20/40/60 and D10. So with an intercept, you would specify coef=2:4 to 
> do an ANOVA-like test for any differences among all four groups.
>
> If you choose not to use an intercept term, then the design matrix is 
> just a different parametrization of the same design. Instead of 
> representing differences, each coefficient simply represents the 
> expression level in a group. So you can still perform the same test, 
> you simply have to specify the contrasts (differences) explicitly:
>
> > anova.contrasts <- makeContrasts(contrasts=c("grpD20-grpD10", 
> "grpD40-grpD10", "grpD60-grpD10"), levels=design)
> > print(anovalike.contrasts)
>         Contrasts
> Levels   grpD20-grpD10 grpD40-grpD10 grpD60-grpD10
>   grpD10            -1            -1            -1
>   grpD20             1             0             0
>   grpD40             0             1             0
>   grpD60             0             0             1
> > glmLRT(fit2, contrast=anova.contrasts)
>
> Hope this helps,
>
> -Ryan Thompson
>
>
> On 03/11/2013 05:18 PM, Luis alberto Martinez lopez wrote:
>> Hello everyone:
>>
>> a wish to know if there is one way of doing an anova like test in edgeR, comparing all the possibilities between the groups defined in the experimental design.
>>
>> in the users guide the authors explain one way anova test, but the comparison is done only against the intercept group, is there  a way of doing a general test comparing all the posibilities and not only against the intercept?
>>
>> also i'd like to know what does the following differential expression analysis does because the output is different when i include an intercept, i can't deduce what comparison are being doing....
>>
>>> grp<-c("D10","D10","D20","D20","D40","D40","D60","D60")
>>>    dge = DGEList(counts=general.m, group=grp)
>>>    dge = calcNormFactors(dge)
>>> design <- model.matrix(~0+grp, data=dge$samples)
>> Warning message:
>> In model.matrix.default(~0 + grp, data = dge$samples) :
>>    variable 'grp' converted to a factor
>>>    design
>>         grpD10 grpD20 grpD40 grpD60
>> i10re1      1      0      0      0
>> i10re2      1      0      0      0
>> i20re1      0      1      0      0
>> i20re2      0      1      0      0
>> i40re1      0      0      1      0
>> i40re2      0      0      1      0
>> i60re1      0      0      0      1
>> i60re2      0      0      0      1
>> attr(,"assign")
>> [1] 1 1 1 1
>> attr(,"contrasts")
>> attr(,"contrasts")$grp
>> [1] "contr.treatment"
>>
>>> dge <- estimateGLMCommonDisp(dge,design)
>>> dge <- estimateGLMTrendedDisp(dge,design)
>> Loading required package: splines
>>>    dge <- estimateGLMTagwiseDisp(dge,design)
>>>    fit2 <- glmFit(dge, design)
>>> lrt.anova_like <- glmLRT(fit2, coef=2:4)
>>> lrt.anova_like
>> An object of class "DGELRT"
>> $coefficients
>>                grpD10    grpD20    grpD40    grpD60
>> comp101_c0 -12.55010 -16.12974 -13.27081 -16.12974
>> comp103_c0 -13.85721 -13.35290 -12.61155 -13.96824
>> comp120_c0 -12.33347 -13.98877 -16.12974 -16.12974
>> comp179_c0 -12.55135 -13.98566 -16.12974 -13.96827
>> comp192_c0 -13.85721 -13.98569 -13.27548 -12.67228
>> 23484 more rows ...
>>
>> $fitted.values
>>                i10re1    i10re2    i20re1    i20re2   i40re1   i40re2    i60re1
>> comp101_c0 2.0108643 1.9941877 0.0000000 0.0000000 0.841309 1.168177 0.0000000
>> comp103_c0 0.5020816 0.4979177 1.1140306 0.8859780 1.674705 2.325366 0.5414430
>> comp120_c0 2.5112747 2.4904480 0.5552870 0.4416145 0.000000 0.000000 0.0000000
>> comp179_c0 2.0083306 1.9916749 0.5570245 0.4429963 0.000000 0.000000 0.5414282
>> comp192_c0 0.5020816 0.4979177 0.5570036 0.4429797 0.837345 1.162673 2.1657497
>>                i60re2
>> comp101_c0 0.0000000
>> comp103_c0 0.4585719
>> comp120_c0 0.0000000
>> comp179_c0 0.4585593
>> comp192_c0 1.8342685
>> 23484 more rows ...
>>
>> $deviance
>> comp101_c0 comp103_c0 comp120_c0 comp179_c0 comp192_c0
>>    4.386443   3.070288   2.848548   3.917710   2.629182
>> 23484 more elements ...
>>
>> $df.residual
>> [1] 4 4 4 4 4
>> 23484 more elements ...
>>
>> $abundance
>> [1] -13.63381 -13.35715 -13.64240 -13.64482 -13.35716
>> 23484 more elements ...
>>
>> $design
>>         grpD10 grpD20 grpD40 grpD60
>> i10re1      1      0      0      0
>> i10re2      1      0      0      0
>> i20re1      0      1      0      0
>> i20re2      0      1      0      0
>> i40re1      0      0      1      0
>> i40re2      0      0      1      0
>> i60re1      0      0      0      1
>> i60re2      0      0      0      1
>> attr(,"assign")
>> [1] 1 1 1 1
>> attr(,"contrasts")
>> attr(,"contrasts")$grp
>> [1] "contr.treatment"
>>
>>
>> $offset
>>           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
>> [1,] 13.27698 13.26865 13.52514 13.29609 13.15729 13.48553 13.47707 13.31095
>> [2,] 13.27698 13.26865 13.52514 13.29609 13.15729 13.48553 13.47707 13.31095
>> [3,] 13.27698 13.26865 13.52514 13.29609 13.15729 13.48553 13.47707 13.31095
>> [4,] 13.27698 13.26865 13.52514 13.29609 13.15729 13.48553 13.47707 13.31095
>> [5,] 13.27698 13.26865 13.52514 13.29609 13.15729 13.48553 13.47707 13.31095
>> 23484 more rows ...
>>
>> $dispersion
>> [1] 0.2161256308 0.0003318302 0.0634519684 0.0003270568 0.0003318016
>> 23484 more elements ...
>>
>> $method
>> [1] "oneway"
>>
>> $samples
>>         group lib.size norm.factors
>> i10re1   D10   594988    0.9808651
>> i10re2   D10   590850    0.9795431
>> i20re1   D20   723182    1.0343009
>> i20re2   D20   567524    1.0481805
>> i40re1   D40   452242    1.1449003
>> i40re2   D40   671656    1.0703965
>> i60re1   D60   801682    0.8892299
>> i60re2   D60   685351    0.8809633
>>
>> $table
>>             logFC.grpD20 logFC.grpD40 logFC.grpD60    logCPM          LR
>> comp101_c0    -23.27030    -19.14573    -23.27030 0.2621324    647.0115
>> comp103_c0    -19.26416    -18.19462    -20.15192 0.6612798 194043.3654
>> comp120_c0    -20.18153    -23.27030    -23.27030 0.2497509   1999.3276
>> comp179_c0    -20.17704    -23.27030    -20.15196 0.2462508 196431.3868
>> comp192_c0    -20.17709    -19.15248    -18.28224 0.6612611 194057.4372
>>                    PValue
>> comp101_c0 6.475824e-140
>> comp103_c0  0.000000e+00
>> comp120_c0  0.000000e+00
>> comp179_c0  0.000000e+00
>> comp192_c0  0.000000e+00
>> 23484 more rows ...
>>
>> $comparison
>> [1] "grpD20" "grpD40" "grpD60"
>>
>> $df.test
>> [1] 3 3 3 3 3
>> 23484 more elements ...
>>
>>
>>> top<-topTags(lrt.anova_like)
>>> top
>> Coefficient:  grpD20 grpD40 grpD60
>>              logFC.grpD20 logFC.grpD40 logFC.grpD60    logCPM         LR PValue
>> comp628_c0      -23.2703    -23.27030    -19.23749 0.2462499 196434.924      0
>> comp2607_c0     -23.2703    -20.07173    -18.68223 0.4686227 195146.814      0
>> comp2858_c0     -23.2703    -23.27030    -16.85200 1.1206869 192585.615      0
>> comp2885_c0     -23.2703    -18.60303    -17.96733 0.6593429   4326.878      0
>> comp2904_c0     -23.2703    -23.27030    -16.21482 1.8318474   2065.059      0
>> comp2912_c0     -23.2703    -17.88139    -20.15196 1.1207038 192646.480      0
>> comp2933_c0     -23.2703    -18.59551    -18.68219 0.4686259 195111.455      0
>> comp2950_c0     -23.2703    -19.16586    -17.97899 0.4607410   1715.721      0
>> comp2999_c0     -23.2703    -23.27030    -23.27030 1.2524134   3853.539      0
>> comp3028_c0     -23.2703    -19.15248    -17.30519 0.9831770 192775.649      0
>>              FDR
>> comp628_c0    0
>> comp2607_c0   0
>> comp2858_c0   0
>> comp2885_c0   0
>> comp2904_c0   0
>> comp2912_c0   0
>> comp2933_c0   0
>> comp2950_c0   0
>> comp2999_c0   0
>> comp3028_c0   0
>>> lrt.anova_like <- glmLRT(fit2, coef=1:4)
>> Error in mglmLevenberg(y, design = design, dispersion = dispersion, offset = offset,  :
>>    BLAS/LAPACK routine 'DGEMM ' gave error code -13
>>
>>
>>
>>
>>
>> The DGELRT object says that is comparing
>>
>> $comparison
>> [1] "grpD20" "grpD40" "grpD60"
>>
>> but what exactly this means?
>>
>> Luis
>> 	[[alternative HTML version deleted]]
>>
>>
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor@r-project.org  <mailto:Bioconductor@r-project.org>
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
>


	[[alternative HTML version deleted]]

