I am doing differential expression on a fairly large RNA-Seq dataset, where multiple organs are obtained from the same mice. I would like to model the effects of individual mice using limma's random effects via duplicateCorrelation
. As library sizes are very consistent across all samples (and I might want to play around with genas
), I was thinking of using limma-trend.
I've seen some posts on using voom
and voomWithArrayWeights
together with duplicateCorrelation
, but nothing on limma-trend.
Given that a standard limma-trend pipeline goes as follows (starting from a count matrix EM
and design mod
):
dge <- DGEList(EM) dge <- calcNormFactors(dge) v <- cpm(dge, log=TRUE, prior.count=3) fit <- lmFit(v, design=mod) eb <- eBayes(fit, trend=TRUE, robust=TRUE) decideTests(eb) topTable(eb) ...
How would this be modified to accommodate duplicateCorrelation
and/or arrayWeights
?
+1 but, regarding the last sentence, there's no need to iterate arrayWeights() and duplicateCorrelation(). arrayWeights doesn't even take correlation as an argument.
Also, I don't recommend that you iteration duplicateCorrelation() and voom(). I don't do that myself. Sure you can, but I don't see any real advantage in doing so.
In line with the original question, which sequence of function calls would you use in that case for combining duplicateCorrelation with voom/voomWithArrayWeights?