Question: limma-trend with random effects and/or array weights
0
2.1 years ago by
maltethodberg140
Sweden
maltethodberg140 wrote:

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

modified 2.1 years ago by Gordon Smyth38k • written 2.1 years ago by maltethodberg140
Answer: limma-trend with random effects and/or array weights
4
2.1 years ago by
Gordon Smyth38k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth38k wrote:

The reason why there aren't any posts is because there's no special issue, everything just works as usual:

aw <- arrayWeights(v, design)
d <- duplicateCorrelation(v, design, block=mouse, weights=aw)
fit <- lmFit(v, design=mod, weights=aw, block=mouse, correlation=d\$consensus)
eb <- eBayes(fit, trend=TRUE, robust=TRUE)


Note also that you should filter low expression genes after running cpm() but before running any of the above.

ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by Gordon Smyth38k
Answer: limma-trend with random effects and/or array weights
3
2.1 years ago by
The Scripps Research Institute, La Jolla, CA
Ryan C. Thompson7.4k wrote:

The issues that occur when mixing voom with arrayWeights and/or duplicateCorrelation happen because the voom step comes before the arrayWeights or duplicateCorrelation steps, but in order to be fully accurate, the voom step needs information from those other steps. The solution implemented in voomWithQualityWeights is essentially to alternate running voom and arrayWeights until convergence. (In practice, convergence is pretty much guaranteed after 2 runs of each, so that's what voomWithQualityWeights actually does instead of checking any convergence criterion.) I've written a helper function to do something similar with voom and duplicateCorrelation here: https://github.com/DarwinAwardWinner/CD4-csaw/blob/e6de96bcefcb3fcd928ea972dc4bc2d511e6ddbd/scripts/utilities.R#L392-L445. Note that I have not thoroughly investigated the convergence behavior of this function, but I'm pretty sure it should converge. If you use it, you should audit the code to make sure it's doing what you expect.

However, if you are using limma-trend instead of voom, I don't think there is any issue. The eBayes step comes after both arrayWeights and duplicateCorrelation, so there is no chicken-and-egg problem as there is with voom. After taking the logCPM, you just proceed as you would for microarray data, as detailed in the manual.

I suppose you might still need to use the trick of iterating back and forth until convergence if you want to use both arrayWeights and duplicateCorrelation in the same pipeline. (Edit: never mind, this would only be needed for voomWithQualityWeights instead of arrayWeights.)

ADD COMMENTlink modified 2.1 years ago • written 2.1 years ago by Ryan C. Thompson7.4k
1

+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.

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by Gordon Smyth38k

In line with the original question, which sequence of function calls would you use in that case for combining duplicateCorrelation with voom/voomWithArrayWeights?