[BioC] EdgeR multi-factors testing questions

[Yanzhu [guest]]

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        



--
Sent via the guest posting facility at bioconductor.org.