[BioC] EdgeR multi-factor testing questions

[Yanzhu [guest]]

I’m currently using EdgeR to analyze RNASeq data. I have one question of whether I have chosen the correct method for analyzing my multi-factorial experiment.  
I have 3 factors: factor A with 16 levels (from level0 to level15), factor B with three levels (from level0 to leve2), and factor Sex with two levels (male and female).  I’m interested in testing all main effects: A, B and Sex, all two-way interaction terms:  A:B, A:Sex and B:Sex; and the three-way interaction term: A:B:Sex. 

My code:

group<-paste(A,B,Sex,sep=".")
design<-model.matrix(~A+B+Sex+A:B+A:Sex+B:Sex+A:B:Sex)
y<-DGEList(counts=readcounts,group=group)
y<-calcNormFactors(y,method="TMM")


y<-estimateGLMCommonDisp(y,design)
y<-estimateGLMTagwiseDisp(y,design)


fit<-glmFit(y,design)
colnames(fit)

colnames(fit)
 [1] "(Intercept)" "Alevel1" "Alevel2" "Alevel3" "Alevel4"        [6] "Alevel5" "Alevel6" "Alevel7" "Alevel8" "Alevel9"             
[11]"Alevel10" "Alevel11" "Alevel12" "Alevel13" "Alevel14"             
[16] "Alevel15" "Blevel1" "Blevel2" "Sexmale" "Alevel1:Blevel1"        

[21] "Alevel2:Blevel1" "Alevel3:Blevel1" "Alevel4:Blevel1" "Alevel5:Blevel1" "Alevel6:Blevel1"        

[26] "Alevel7:Blevel1" "Alevel8:Blevel1" "Alevel9:Blevel1" "Alevel10:Blevel1" "Alevel11:Blevel1"        

[31] "Alevel12:Blevel1" "Alevel13:Blevel1" "Alevel14:Blevel1" "Alevel15:Blevel1" "Alevel1:Blevel2"        

[36] "Alevel2:Blevel2"  "Alevel3:Blevel2" "Alevel4:Blevel2" "Alevel5:Blevel2" "Alevel6:Blevel2"        

[41] "Alevel7:Blevel2"  "Alevel8:Blevel2" "Alevel9:Blevel2"   "Alevel10:Blevel2" "Alevel11:Blevel2"        

[46] "Alevel12:Blevel2" "Alevel13:Blevel2" "Alevel14:Blevel2" "Alevel15:Blevel2" "Alevel1:Sexmale"  
   
[51] "Alevel2:Sexmale" "Alevel3:Sexmale" "Alevel4:Sexmale" "Alevel5:Sexmale" "Alevel6:Sexmale"    
 
[56] "Alevel7:Sexmale" "Alevel8:Sexmale" "Alevel9:Sexmale" "Alevel10:Sexmale" "Alevel11:Sexmale"
     
[61] "Alevel12:Sexmale" "Alevel13:Sexmale" "Alevel14:Sexmale" "Alevel15:Sexmale" "Blevel1:Sexmale"   
     
[66] "Blevel2:Sexmale" "Alevel1:Blevel1:Sexmale" "Alevel2:Blevel1:Sexmale" "Alevel3:Blevel1:Sexmale" "Alevel4:Blevel1:Sexmale"

[71] "Alevel5:Blevel1:Sexmale" "Alevel6:Blevel1:Sexmale" "Alevel7:Blevel1:Sexmale" "Alevel8:Blevel1:Sexmale" "Alevel9:Blevel1:Sexmale"

[76] "Alevel10:Blevel1:Sexmale" "Alevel11:Blevel1:Sexmale" "Alevel12:Blevel1:Sexmale" "Alevel13:Blevel1:Sexmale" "Alevel14:Blevel1:Sexmale"

[81] "Alevel15:Blevel1:Sexmale" "Alevel1:Blevel2:Sexmale" "Alevel2:Blevel2:Sexmale" "Alevel3:Blevel2:Sexmale" "Alevel4:Blevel2:Sexmale"

[86] "Alevel5:Blevel2:Sexmale" "Alevel6:Blevel2:Sexmale" "Alevel7:Blevel2:Sexmale" "Alevel8:Blevel2:Sexmale" "Alevel9:Blevel2:Sexmale"

[91] "Alevel10:Blevel2:Sexmale" "Alevel11:Blevel2:Sexmale" "Alevel12:Blevel2:Sexmale" "Alevel13:Blevel2:Sexmale" "Alevel14:Blevel2:Sexmale" 
[96] "Alevel15:Blevel2:Sexmale"

(1) To test the main effect, I use the following code:
lrt_A<-glmLRT(fit,coef=c(2:16))
lrt_B<-glmLRT(fit,coef=c(17:18))
lrt_Sex<-glmLRT(fit,coef=19)

(2) To test the two-way interaction terms, I use the following code:
lrt_AB<-glmLRT(fit,coef=c(20:49))

lrt_ASex<-glmLRT(fit,coef=c(50:64))

lrt_BSex<-glmLRT(fit,coef=c(65:66))

(3) To test the three-way interaction terms, I use the following code:
lrt_ABSex<-glmLRT(fit,coef=c(67:96))

My question: is my code correct for each of the testing above?



Thank you!










 -- 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] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] edgeR_3.2.4  limma_3.16.8

loaded via a namespace (and not attached):
[1] tools_3.0.1


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