Entering edit mode
Luis alberto Martinez lopez
▴
10
@luis-alberto-martinez-lopez-5817
Last seen 10.4 years ago
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]]