**D**ifferential expression for **re**pe**a**ted **m**easures (dream) uses a linear model model to increase power and decrease false positives for RNA-seq datasets with multiple measurements per individual. The analysis fits seamlessly into the widely used workflow of limma/voom (Law et al. 2014). Dream uses a linear model model to increase power and decrease false positives for RNA-seq datasets with repeated measurements. Dream achieves this by combining multiple statistical concepts into a single statistical model. The model includes:

- flexible modeling of repeated measures gene expression data
- precision weights to model measurement error in RNA-seq counts
- ability to model multiple random effects
- random effects estimated separately for each gene
- hypothesis testing for fixed effects in linear mixed models
- small sample size hypothesis test

Dream also includes multi-threaded analysis across thousands of genes on a multi-core machine.

This tutorial assumes that the reader is familiar with the limma/voom workflow for RNA-seq. Process raw count data using limma/voom.

```
library("variancePartition")
library("edgeR")
library("BiocParallel")
data(varPartDEdata)
# filter genes by number of counts
isexpr <- rowSums(cpm(countMatrix) > 0.1) >= 5
# Standard usage of limma/voom
dge <- DGEList(countMatrix[isexpr, ])
dge <- calcNormFactors(dge)
# make this vignette faster by analyzing a subset of genes
dge <- dge[1:1000, ]
```

Limma has a built-in approach for analyzing repeated measures data using `duplicateCorrelation()`

. The model can handle a single random effect, and forces the magnitude of the random effect to be the same across all genes.

```
# apply duplicateCorrelation is two rounds
design <- model.matrix(~Disease, metadata)
vobj_tmp <- voom(dge, 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(dge, 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)
```

The dream method replaces 4 core functions of limma with a linear mixed model.

`voomWithDreamWeights()`

replaces`voom()`

to estimate precision weights`dream()`

replaces`lmFit()`

to estimate regression coefficients.

`variancePartition::eBayes()`

replaces`limma::eBayes()`

to apply empircal Bayes shrinkage on linear mixed models.`variancePartition::topTable()`

replaces`limma::topTable()`

to give seamless access to results from`dream()`

.

For models with only fixed effects, `variancePartition::eBayes()`

, and `variancePartition::topTable()`

work seamlessly and give results equivalent to the `limma`

functions with the same name. From the user perspective, the `dream()`

workflow is the same as `limma`

since the statistical differences are handled behind the scenes.

```
# Specify parallel processing parameters
# this is used implicitly by dream() to run in parallel
param <- SnowParam(4, "SOCK", progressbar = TRUE)
# The variable to be tested must be a fixed effect
form <- ~ Disease + (1 | Individual)
# estimate weights using linear mixed model of dream
vobjDream <- voomWithDreamWeights(dge, form, metadata, BPPARAM = param)
# Fit the dream model on each gene
# For the hypothesis testing, by default,
# dream() uses the KR method for <= 20 samples,
# otherwise it uses the Satterthwaite approximation
fitmm <- dream(vobjDream, form, metadata)
fitmm <- eBayes(fitmm)
```

```
## (Intercept) Disease1
## sample_01 1 0
## sample_02 1 0
## sample_03 1 0
```

```
# Get results of hypothesis test on coefficients of interest
topTable(fitmm, coef = "Disease1", number = 3)
```

```
## logFC AveExpr t P.Value adj.P.Val B
## ENST00000283033.5 gene=TXNDC11 1.556169 3.567616 31.02010 3.101602e-21 3.101602e-18 38.15514
## ENST00000257181.9 gene=PRPF38A 1.380508 4.398265 21.59327 1.682434e-17 8.412170e-15 29.93328
## ENST00000525790.1 gene=TDRKH 1.508262 3.184918 19.51539 1.774506e-16 5.915021e-14 27.77815
## z.std
## ENST00000283033.5 gene=TXNDC11 9.459250
## ENST00000257181.9 gene=PRPF38A 8.513856
## ENST00000525790.1 gene=TDRKH 8.236413
```

Since dream uses an estimated degrees of freedom value for each hypothsis test, the degrees of freedom is different for each gene here. Therefore, the t-statistics are not directly comparable since they have different degrees of freedom. In order to be able to compare test statistics, we report `z.std`

which is the p-value transformed into a signed z-score. This can be used for downstream analysis.

Note that if a random effect is not specified, `dream()`

automatically uses `lmFit()`

, but the user must run `eBayes()`

afterward.

You can also perform a hypothesis test of the *difference* between two or more coefficients by using a contrast matrix. The contrasts are evaluated at the time of the model fit and the results can be extracted with `topTable()`

. This behaves like `makeContrasts()`

and `contrasts.fit()`

in limma.

Multiple contrasts can be evaluated at the same time, in order to save computation time. Make sure to inspect your contrast matrix to confirm it is testing what you intend.

```
form <- ~ 0 + DiseaseSubtype + Sex + (1 | Individual)
L <- makeContrastsDream(form, metadata,
contrasts = c(
compare2_1 = "DiseaseSubtype2 - DiseaseSubtype1",
compare1_0 = "DiseaseSubtype1 - DiseaseSubtype0"
)
)
# Visualize contrast matrix
plotContrasts(L)
```

```
# fit dream model with contrasts
fit <- dream(vobjDream, form, metadata, L)
fit <- eBayes(fit)
# get names of available coefficients and contrasts for testing
colnames(fit)
```

```
## [1] "compare2_1" "compare1_0" "DiseaseSubtype0" "DiseaseSubtype1" "DiseaseSubtype2"
## [6] "SexM"
```

```
## logFC AveExpr t P.Value adj.P.Val B
## ENST00000355624.3 gene=RAB11FIP2 -0.9493298 5.260267 -5.114693 3.697208e-05 0.03697208 0.1471393
## ENST00000593466.1 gene=DDA1 -1.7268844 3.901505 -3.700715 1.206260e-03 0.60313007 -1.2754458
## ENST00000200676.3 gene=CETP 1.4778546 3.723402 3.910328 3.396490e-03 0.99951203 -1.4380854
## z.std
## ENST00000355624.3 gene=RAB11FIP2 -4.125630
## ENST00000593466.1 gene=DDA1 -3.237396
## ENST00000200676.3 gene=CETP 2.929371
```

So far contrasts have only involved the difference between two coefficients. But contrasts can also compare any linear combination of coefficients. Here, consider comparing `DiseaseSubtype0`

to the mean of `DiseaseSubtype1`

and `DiseaseSubtype2`

. Note you can also customize the name of the contrast.

```
L2 <- makeContrastsDream(form, metadata,
contrasts =
c(Test1 = "DiseaseSubtype0 - (DiseaseSubtype1 + DiseaseSubtype2)/2")
)
plotContrasts(L2)
```

```
# fit dream model to evaluate contrasts
fit <- dream(vobjDream[1:10, ], form, metadata, L = L2)
fit <- eBayes(fit)
topTable(fit, coef = "Test1", number = 3)
```

```
## logFC AveExpr t P.Value adj.P.Val B
## ENST00000456159.1 gene=MET -0.9828664 2.458905 -6.974652 3.534137e-08 3.534137e-07 8.661088
## ENST00000570099.1 gene=YPEL3 -0.9510245 2.063287 -6.648276 6.068289e-07 2.176720e-06 7.229240
## ENST00000418210.2 gene=TMEM64 -1.0342889 4.715361 -6.654887 6.530161e-07 2.176720e-06 7.691971
## z.std
## ENST00000456159.1 gene=MET -5.512678
## ENST00000570099.1 gene=YPEL3 -4.989031
## ENST00000418210.2 gene=TMEM64 -4.974840
```

Joint hypothesis testing of multiple coefficients at the same time can be performed by using an F-test. Just like in limma, the results can be extracted using `topTable()`

```
# extract results from first contrast
topTable(fit, coef = c("DiseaseSubtype2", "DiseaseSubtype1"), number = 3)
```

```
## DiseaseSubtype2 DiseaseSubtype1 AveExpr F P.Value
## ENST00000418210.2 gene=TMEM64 5.300955 5.211628 4.715361 744.2802 5.226122e-23
## ENST00000555834.1 gene=RPS6KL1 5.662671 5.719164 5.272059 657.9501 2.376578e-22
## ENST00000589123.1 gene=NFIC 6.545179 6.180997 5.855330 519.2397 4.312572e-21
## adj.P.Val F.std
## ENST00000418210.2 gene=TMEM64 5.226122e-22 51.30579
## ENST00000555834.1 gene=RPS6KL1 1.188289e-21 49.79121
## ENST00000589123.1 gene=NFIC 1.437524e-20 46.89275
```

Since dream uses an estimated degrees of freedom value for each hypothsis test, the degrees of freedom is different for each gene here. Therefore, the F-statistics are not directly comparable since they have different degrees of freedom. In order to be able to compare test statistics, we report `F.std`

which is the p-value transformed into an F-statistic with \(df_1\) the number of coefficients tested and \(df_2=\infty\). This can be used for downstream analysis.

For small datasets, the Kenward-Roger method can be more powerful. But it is **substantially** more computationally intensive.

Dream and variancePartition share the same underlying linear mixed model framework. A variancePartition analysis can indicate important variables that should be included as fixed or random effects in the dream analysis.

Here we compare p-values from from `dream()`

and `duplicateCorrelation`

. In order to understand the empircal difference between them, we can plot the \(-\log_{10}\) p-values from each method.

```
# Compare p-values and make plot
p1 <- topTable(fitDupCor, coef = "Disease1", number = Inf, sort.by = "none")$P.Value
p2 <- topTable(fitmm, number = Inf, sort.by = "none")$P.Value
plotCompareP(p1, p2, vp$Individual, dupcor$consensus)
```

The duplicateCorrelation method estimates a single variance term genome-wide even though the donor contribution of a particular gene can vary substantially from the genome-wide trend. Using a single value genome-wide for the within-donor variance can reduce power and increase the false positive rate in a particular, reproducible way. Let \(\tau^2_g\) be the value of the donor component for gene \(g\) and \(\bar{\tau}^2\) be the genome-wide mean. For genes where \(\tau^2_g>\bar{\tau}^2\), using \(\bar{\tau}^2\) under-corrects for the donor component so that it increases the false positive rate compared to using \(\tau^2_g\). Conversely, for genes where \(\tau^2_g<\bar{\tau}^2\), using \(\bar{\tau}^2\) over-corrects for the donor component so that it decreases power. Increasing sample size does not overcome this issue. The dream method overcomes this issue by using a \(\tau^2_g\).

Here, the \(-\log_{10}\) p-values from both methods are plotted and colored by the donor contribution estiamted by variancePartition. The green value indicates \(\bar{\tau}^2\), while red and blue indicate higher and lower values, respectively. When only one variance component is used and the contrast matrix is simple, the effect of using dream versus duplicateCorrelation is determined by the comparison of \(\tau^2_g\) to \(\bar{\tau}^2\):

dream can increase the \(-\log_{10}\) p-value for genes with a lower donor component (i.e. \(\tau^2_g<\bar{\tau}^2\)) and decrease \(-\log_{10}\) p-value for genes with a higher donor component (i.e. \(\tau^2_g>\bar{\tau}^2\))

Note that using more variance components or a more complicated contrast matrix can make the relationship more complicated.

variancePartition functions including `dream()`

, `fitExtractVarPartModel()`

and `fitVarPartModel()`

can take advange of multicore machines to speed up analysis. It uses the BiocParallel package to manage the parallelization.

- Specify parameters with the BPPARAM argument.

```
# Request 4 cores, and enable the progress bar
# This is the ideal for Linux, OS X and Windows
param <- SnowParam(4, "SOCK", progressbar = TRUE)
fitmm <- dream(vobjDream, form, metadata, BPPARAM = param)
```

By default `BPPARAM = SerialParam()`

.

```
## R Under development (unstable) (2023-11-11 r85510)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_GB
## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] edgeR_4.1.1 pander_0.6.5 variancePartition_1.33.2
## [4] BiocParallel_1.37.0 limma_3.59.1 ggplot2_3.4.4
## [7] knitr_1.45
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 farver_2.1.1 dplyr_1.1.3 bitops_1.0-7
## [5] fastmap_1.1.1 digest_0.6.33 lifecycle_1.0.4 statmod_1.5.0
## [9] magrittr_2.0.3 compiler_4.4.0 rlang_1.1.2 sass_0.4.7
## [13] tools_4.4.0 utf8_1.2.4 yaml_2.3.7 labeling_0.4.3
## [17] plyr_1.8.9 KernSmooth_2.23-22 withr_2.5.2 purrr_1.0.2
## [21] numDeriv_2016.8-1.1 BiocGenerics_0.49.1 grid_4.4.0 aod_1.3.2
## [25] fansi_1.0.5 caTools_1.18.2 colorspace_2.1-0 scales_1.2.1
## [29] gtools_3.9.4 iterators_1.0.14 MASS_7.3-60.1 cli_3.6.1
## [33] mvtnorm_1.2-3 rmarkdown_2.25 generics_0.1.3 reshape2_1.4.4
## [37] minqa_1.2.6 cachem_1.0.8 stringr_1.5.0 splines_4.4.0
## [41] parallel_4.4.0 matrixStats_1.1.0 vctrs_0.6.4 boot_1.3-28.1
## [45] Matrix_1.6-3 jsonlite_1.8.7 pbkrtest_0.5.2 locfit_1.5-9.8
## [49] jquerylib_0.1.4 tidyr_1.3.0 snow_0.4-4 glue_1.6.2
## [53] nloptr_2.0.3 codetools_0.2-19 stringi_1.8.1 gtable_0.3.4
## [57] EnvStats_2.8.1 lme4_1.1-35.1 lmerTest_3.1-3 munsell_0.5.0
## [61] tibble_3.2.1 remaCor_0.0.16 pillar_1.9.0 htmltools_0.5.7
## [65] gplots_3.1.3 R6_2.5.1 Rdpack_2.6 evaluate_0.23
## [69] lattice_0.22-5 Biobase_2.63.0 highr_0.10 rbibutils_2.2.16
## [73] backports_1.4.1 RhpcBLASctl_0.23-42 broom_1.0.5 fANCOVA_0.6-1
## [77] corpcor_1.6.10 bslib_0.5.1 Rcpp_1.0.11 nlme_3.1-163
## [81] xfun_0.41 pkgconfig_2.0.3
```

<>

Law, C. W., Y. Chen, W. Shi, and G. K. Smyth. 2014. “Voom: precision weights unlock linear model analysis tools for RNA-seq read counts.” *Genome Biology* 15 (2): R29. https://doi.org/10.1186/gb-2014-15-2-r29.