In the dream vignette it's described that for an RNA-seq repeated measure design the limma analysis needs two
duplicateCorrelation rounds, though this differs from the limma user's guide and manual. Below for convenience is the design:
> metadata Individual Disease Experiment Sex DiseaseSubtype sample_01 ID1 0 sample_01 F 0 sample_02 ID1 0 sample_02 F 0 sample_03 ID2 0 sample_03 M 0 sample_04 ID2 0 sample_04 M 0 sample_05 ID3 0 sample_05 F 0 sample_06 ID3 0 sample_06 F 0 sample_07 ID4 0 sample_07 M 0 sample_08 ID4 0 sample_08 M 0 sample_09 ID5 0 sample_09 F 0 sample_10 ID5 0 sample_10 F 0 sample_11 ID6 0 sample_11 M 0 sample_12 ID6 0 sample_12 M 0 sample_13 ID7 1 sample_13 F 1 sample_14 ID7 1 sample_14 F 1 sample_15 ID8 1 sample_15 M 1 sample_16 ID8 1 sample_16 M 1 sample_17 ID9 1 sample_17 F 1 sample_18 ID9 1 sample_18 F 1 sample_19 ID10 1 sample_19 M 2 sample_20 ID10 1 sample_20 M 2 sample_21 ID11 1 sample_21 F 2 sample_22 ID11 1 sample_22 F 2 sample_23 ID12 1 sample_23 M 2 sample_24 ID12 1 sample_24 M 2
and the limma two-round
duplicateCorrelation analysis code from the vignette, where the DE test is an ordinary two-group comparison
# apply duplicateCorrelation is two rounds design = model.matrix( ~ Disease, metadata) vobj_tmp = voom( geneExpr, design, plot=FALSE) dupcor <- duplicateCorrelation(vobj_tmp,design,block=metadata$Individual) # run voom considering the duplicateCorrelation results # in order to compute more accurate precision weights # Otherwise, use the results from the first voom run vobj = voom( geneExpr, design, plot=FALSE, block=metadata$Individual, correlation=dupcor$consensus) # Estimate linear mixed model with a single variance component # Fit the model for each gene, dupcor <- duplicateCorrelation(vobj, design, block=metadata$Individual) # But this step uses only the genome-wide average for the random effect fitDupCor <- lmFit(vobj, design, block=metadata$Individual, correlation=dupcor$consensus) # Fit Empirical Bayes for moderated t-statistics fitDupCor <- eBayes( fitDupCor )
This differs from the limma user's guide which, for similar designs and DE test questions, describes only one
duplicateCorrelation round. Also, the limma manual says that specifying parameters like
voom only passes them to
voom doesn't use it for precision weights calculation. This would also make sense for only needing one round. Which way is the correct way?
I compared the vignette example with two rounds and the limma guide with one round and the p-values are not identical, the one-round has slightly higher p-values and slightly lower log-odds.