Longitudinal differential expression experiment: several subjects, measured once before treatment and twice after treatment at set time points. Goal is the treatment effect. DESeq2 with design ~ treatment + subject underestimates p-values because the post-treatment samples from the same subject are not independent. duplicateCorrelation was suggested, but there is a concern that it inflates degrees of freedom and gives liberal p-values. An alternative under consideration: run voom, then collapse to an inverse-variance weighted mean per subject x treatment cell and fit with lmFit using the summed weights as precision. Question: is the weighted collapse valid, or is duplicateCorrelation the right route, and how should subject be handled (fixed effect vs block)?
Modelling replicate measurements in a longitudinal RNA-seq experiment with voom
nice writeup, and the weighted-collapse idea is reasonable. a couple of thoughts:
your approach is basically the manual version of what limma already does for technical replicates (avearrays / arrayWeights). one caveat: lmFit will treat the collapsed values as independent with known weights, which is a little optimistic since the voom weights were estimated. in practice the effect is usually small, but it means your p-values lean slightly liberal rather than being exact.
on duplicateCorrelation: it doesnt inflate degrees of freedom. it estimates one consensus correlation across all genes and plugs it into the precision weights, treating that estimate as known. thats the source of the "slightly liberal" comment, not df inflation. for repeated measures on the same subjects its the recommended route in the limma user guide, so i wouldnt rule it out:
design <- model.matrix(~ treatment, sample_info)
v <- voom(dge, design)
corfit <- duplicateCorrelation(v, design, block = sample_info$subject)
fit <- lmFit(v, design, block = sample_info$subject, correlation = corfit$consensus)
fit <- eBayes(fit)one detail: pick one way to handle subject, not both. either subject as a fixed effect in the design, or block = subject (random effect). mixing them double-counts. the block version is the cleaner match for "what is the effect of treatment".
on power: with a few subjects theres no free lunch. collapsing to one value per subject x treatment (your approach, or plain avearrays) is the conservative option and behaves like a paired comparison. duplicateCorrelation keeps the individual samples and will usually give you a bit more power, at the cost of that slight liberality. id run duplicateCorrelation as the main analysis and the collapsed version as a sensitivity check, and see whether the conclusions agree.
Source: https://support.bioconductor.org/p/9163807/