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.