Dear all,

I am assigned to analyse RNA-Seq data from a study with a repeated measures design with three time points. From what I understand (after going through Q&As in this Bioconductor forum and reading several manuals), the limma package is able to incorporate the within-subjects correlation via the blocking argument and the duplicateCorrelation command. However, a statistician in my university advised me to use lme4 for loops instead to fit every gene one by one (after voom transformation & using sample-specific weights obtained from the command).

The latter approach took a LOT of time to run (few hours), while limma only took a few seconds. In addition, working within the limma workflow enables me to use the included commands such as camera and goana seamlessly.

What do you usually do when you analyse repeated measures RNA-Seq data? How many of you use custom-scripts with lme4/nlme to fit mixed model for each genes? How common is the practice anyway?

Thank you very much for your time.

Best regards,
Mikhael

It is very rare to use lme4 for a microarray differential expression analysis and even rarer for RNA-seq. Some of the reasons are explained in the previous thread linked to by Aaron Lun.

I'm the limma author, so you might consider me biased, but I would not seriously consider using lme4 for an RNA-seq differential expression analysis. Apart from speed and convenience, the loss of power and the increase in false discovery rate would be unacceptable to me.

I have used lme4 myself for one particular microarray study:

    https://doi.org/10.1371/journal.pone.0019556

but the purpose of that study was to estimate the variance components, not to do differential expression analysis.

There are trade-offs to both approaches. If you fit a separate model for every gene, you lose the core feature of limma: empirical Bayes information sharing when estimating the variance for each gene. The other main difference is that limma's mixed model functionality estimates a single correlation value for the random effect that is shared across all genes when fitting the model. If you have hundreds of samples, it might make sense to fit a separate model for each gene using lme4 or some other package. However, most RNA-seq experiments don't have sufficient sample sizes to robustly estimate parameters for each gene independently, so you generally want to stick to limma or other methods (like edgeR and DESeq2) that share information between genes to make up for the small sample sizes.