Dear Yanzhu,
The culprit is your code here:
glmFit(y,design,offset=offset )
whereby you are feeding your own offset to glmFit() and hence
overwriting
any normalization that has been done. Why are you doing this? If you
just remove the offset=offset argument and allow edgeR to work as
normal,
then the results will be correct.
Best wishes
Gordon
> Date: Tue, 28 Jan 2014 06:38:26 -0800 (PST)
> From: "Yanzhu [guest]" <guest at="" bioconductor.org="">
> To: bioconductor at r-project.org, mlinyzh at gmail.com
> Subject: [BioC] EdgeR multi-factors testing questions
>
>
> Hi Ryan,
>
> The P-values are exactly identical, please see the example of the
testing the three-way interaction (as below (II)):
>
> (I) I have checked the normalization factors for both methods,
> (i) TMM:
> head(y$sample$norm.factors)
> [1] 0.9404963 0.9056276 0.9928999 0.9961922 1.0109106 0.9115558
>
> (ii) UQ:
> head(y$sample$norm.factors)
> [1] 0.9883357 0.9625269 1.0414372 1.0143681 1.1134091 0.8992270
>
> (II) Example: the following is how I test the three-way interaction
term:
> (i) TMM:
> group<-paste(L,S,R,sep=".")
> design<-model.matrix(~L+R+S+L:R+L:S+R:S+L:R:S)
> y<-DGEList(counts=counts,group=group)
> y<-calcNormFactors(y,method="TMM")
> offset=log((y$sample)[,2])
> y<-estimateGLMCommonDisp(y,design)
> y<-estimateGLMTagwiseDisp(y,design)
>
> fiteTMM_LRS<-glmFit(y,design,offset=offset )
>
> ### testing the three-way interaction term L:R:S
> lrteTMM_LRS<-glmLRT(fiteTMM_LRS,coef=c(67:96))
>
> ### P-values:
> head((lrteTMM_LRS$table)[,33])
> [1] 1.233769e-30 5.648507e-30 4.254337e-06 9.304154e-05 8.504918e-01
1.075495e-41
>
> (ii) UQ:
> ##################### UQ
> group<-paste(L,S,R,sep=".")
> design<-model.matrix(~L+R+S+L:R+L:S+R:S+L:R:S)
> y<-DGEList(counts=counts,group=group)
> y<-calcNormFactors(y,method="upperquartile",p=0.75)
> offset=log((y$sample)[,2])
> y<-estimateGLMCommonDisp(y,design)
> y<-estimateGLMTagwiseDisp(y,design)
>
> fiteUQ_LRS<-glmFit(y,design,offset=offset )
>
> ### testing the three-way interaction term L:R:S
> lrteUQ_LRS<-glmLRT(fiteUQ_LRS,coef=c(67:96))
>
> ### P-values:
> head((lrteUQ_LRS$table)[,33])
> [1] 1.233769e-30 5.648507e-30 4.254337e-06 9.304154e-05 8.504918e-01
1.075495e-41
>
> The P-values from both normalization methods are exactly identical,
here I only show part of it.
>
>
> Thanks.
>
>
> Yanzhu
>
> -------------------------------------------------------
> I don't think it's surprising. TMM and upper quartile tend to give
very
> similar normalization factors. Are you getting *exactly* identical
> P-values, or only very similar ones? Also, you can check the
> normalization factors themselves to see how much of a difference
there
> is. If "d" is your DGEList object, you can get the normalization
> factors by typeing "d$samples$norm.factors.
>
> -Ryan
>
>
>
> -- output of sessionInfo():
>
> sessionInfo()
> R version 3.0.1 (2013-05-16)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252 LC_MONETARY=English_United States.1252
> [4] LC_NUMERIC=C LC_TIME=English_United
States.1252
>
> attached base packages:
> [1] parallel stats graphics grDevices utils datasets
methods base
>
> other attached packages:
> [1] DESeq_1.12.1 lattice_0.20-15 locfit_1.5-9.1
Biobase_2.20.1 BiocGenerics_0.6.0 edgeR_3.2.4 limma_3.16.8
>
> loaded via a namespace (and not attached):
> [1] annotate_1.38.0 AnnotationDbi_1.22.6 DBI_0.2-7
genefilter_1.42.0 geneplotter_1.38.0 grid_3.0.1
> [7] IRanges_1.18.4 RColorBrewer_1.0-5 RSQLite_0.11.4
splines_3.0.1 stats4_3.0.1 survival_2.37-4
> [13] XML_3.98-1.1 xtable_1.7-1
______________________________________________________________________
The information in this email is confidential and
intend...{{dropped:4}}
