Hi,

So I'm trying to test differential gene expression on case/control study across all time points when accounting for the within subject variability. I know the full and reduced models between which I want to test (if I could use LRT) with but I'm a bit unsure how to apply it to the limma+voom pipeline presented in Limma manual.

Models for LRT testing:

Full: ~condition + time + condition:subject.nested + condition:time (model.matrix with 36 coefs)

Reduced: ~condition + time + condition:subject.nested (model.matrix with 34 coefs)

The code I've been using in LIMMA+VOOM pipeline:

dge <- DGEList(counts=countdata, samples = coldata, group = coldata$CASE_CONTROL)

keep <- rowSums(cpm(dge)>1) >= 2

dge <- dge[keep, , keep.lib.sizes=FALSE]

dge <- calcNormFactors(dge)

v <- voom(dge, design, plot=TRUE)

fit <- lmFit(v, design)
fit <- eBayes(fit, robust = TRUE)

I know I can output the results with topTable but I'm not quite sure how to define the coef's in proper way so that I could get the answers I want --> which genes are differentially expressed "between" case/control groups across any of the time points when accounting for the within subject variability. I've tried resLIMMAfilt3536 <- topTable(fit, coef=35:36) - but not sure if it the way to achieve wanted answers. What kind of statistical backgroud topTable uses to extract the results from fit? 

Really would appreciate the help.

Setting coef=35:36 will perform a moderated F-test against the null hypothesis that the 35th and 36th coefficients are both equal to zero. I'm guessing that these coefficients correspond to the condition-specific time effect in your full model; so rejections will identify genes that exhibit a time effect in either condition. Of course, this is dependent on correctly interpreting the coefficients, which is not always simple in a complex parametrization.