Hi ALL,

I have questions with limma-voom RNA-Seq data analysis. My project has three factor: A (16 levels), B(three levels) and C(two levels), where A is random effect while B and C are fixed effects.

I would like to test the main effects of B and C, as well as A. To test the effect B, I used the following coding.

Group <- factor(paste(Design$A,Design$B,Design$C,Design$plateID,sep="."))
design <- model.matrix(~B+C)
y <- DGEList(counts=T,group=Group) ## T is the count data
                                                                                                                                                                      
y <- calcNormFactors(y,method="RLE")                  
v <- voom(y,design,plot=TRUE)

corfit <- duplicateCorrelation(v,design,block=Design$A)
corfit$consensus

fitRan <- lmFit(v,design,block=Design$A,correlation=corfit$consensus)
fitRan <- eBayes(fitRan)
topTable(fitRan,coef=c(2,3)) ## for testing the effect of factor B, which has three levels.

My questions:

1. Is the coding above correct for testing the main effect of B?

2. How to test the random effect A?

Thanks,

Yanzhu

Yes, the code is correct for testing for B. (Although usually a filtering step is needed. Have you filtered out very lowly expressed genes?)

You can't test for factor A. It doesn't really make sense to test for a random effect.