Dear Community,

I have a question about the voom-limma of RNA-Seq data analysis. I have multi-factor design with three factors: A (16 levels), B (2 levels) and C (three levels), and A is the random effect.

My question is how to treat A as random effect in the model fitting in limma. I am not quite clear about it even after referring to the limma user manual.

My code are given as below:

design<-model.matrix(~A+B+C+A:B+A:C+B:C+A:B:C)

Group<-factor(paste(Design$A,Design$B,Design$C,sep="."))

y<-DGEList(counts=T,group=Group)               ####T is the RNASeq count data

y<-calcNormFactors(y)                  
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(seq(2,16,1)))

My questions:

1. I want to fit the multi-factor ANOVA model with treating factor A as random effect, Are the codes above correct?

2. As I treat factor A as random effect, I want to the variance explained by factor A? Is there any option in limma report the variation explained by the random effect?

3. I want to test the factor A effect, is it correct that I use  topTable(fitRan,coef=c(seq(2,16,1)),number=Inf)? From the second column to 16th column of coefficients matrix are for factor A's level.

The answer to all your questions is no. You are trying to fit a model with random interactions (A.B, A.C, A.B.C), but limma does not do that, nor does any other Bioconductor package.

Why do you want to do that? It is really a good idea? I have never seen an RNA-seq dataset for which that would be sensible.

If your intention is simply to treat A as a random factor, then you need to remove A and its interactions from the design matrix, because you can't treat a factor as fixed and random at the same time.