Standard RNA-seq processing

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

Limma Analysis

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.

Dream Analysis

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

  1. voomWithDreamWeights() replaces voom() to estimate precision weights
  2. dream() replaces lmFit() to estimate regression coefficients.

Otherwise dream uses the same workflow as limma with topTable(), since any statistical differences are handled behind the scenes.

##           (Intercept) Disease1
## sample_01           1        0
## sample_02           1        0
## sample_03           1        0
##                                    logFC  AveExpr        t      P.Value    adj.P.Val     z.std
## ENST00000283033.5 gene=TXNDC11  1.556233 3.567624 68.62592 3.701410e-27 3.701410e-24 10.793327
## ENST00000257181.9 gene=PRPF38A  1.380549 4.398270 35.53035 6.311837e-21 3.155919e-18  9.384662
## ENST00000421974.2 gene=ATP6V0E2 1.386821 3.478030 22.35157 1.289327e-16 3.584329e-14  8.274558

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.

Advanced hypothesis testing

Using contrasts to compare coefficients

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.

Make sure to inspect your contrast matrix to confirm it is testing what you intend.

## [1] "L1"              "DiseaseSubtype0" "DiseaseSubtype1" "DiseaseSubtype2" "SexM"
##                                       logFC  AveExpr         t      P.Value  adj.P.Val     z.std
## ENST00000355624.3 gene=RAB11FIP2 -0.9493146 5.260280 -5.384718 2.857043e-05 0.02857043 -4.184572
## ENST00000593466.1 gene=DDA1      -1.7265709 3.901579 -3.516848 2.168824e-03 0.99955525 -3.066084
## ENST00000200676.3 gene=CETP       1.4777422 3.723438  3.732332 8.085140e-03 0.99955525  2.648494

Multiple contrasts can be evaluated at the same time, in order to save computation time:

##                                       logFC  AveExpr         t      P.Value  adj.P.Val     z.std
## ENST00000355624.3 gene=RAB11FIP2 -0.9493146 5.260280 -5.384718 2.857043e-05 0.02857043 -4.184572
## ENST00000593466.1 gene=DDA1      -1.7265709 3.901579 -3.516848 2.168824e-03 0.99955525 -3.066084
## ENST00000200676.3 gene=CETP       1.4777422 3.723438  3.732332 8.085140e-03 0.99955525  2.648494

Comparing multiple coefficients

So far contrasts have only involved the difference between two coefficients. But contrasts can also compare linear combinations of coefficients. A user can create contrasts ‘manually’, as long as the elements sum to 0. Here, consider comparing DiseaseSubtype0 to the mean of DiseaseSubtype1 and DiseaseSubtype2

##                                     logFC  AveExpr         t      P.Value    adj.P.Val     z.std
## ENST00000456159.1 gene=MET     -0.9830788 2.458926 -6.349967 3.385843e-06 3.385843e-05 -4.645908
## ENST00000418210.2 gene=TMEM64  -1.0343236 4.715367 -8.899159 1.603605e-05 5.403575e-05 -4.313954
## ENST00000555834.1 gene=RPS6KL1 -0.8954794 5.272063 -5.636526 1.621072e-05 5.403575e-05 -4.311559

Joint hypothesis test of multiple coefficients

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()

##                                DiseaseSubtype2 DiseaseSubtype1  AveExpr         F      P.Value
## ENST00000418210.2 gene=TMEM64         5.301001        5.211674 4.715367 1330.8471 6.020801e-12
## ENST00000555834.1 gene=RPS6KL1        5.662699        5.719196 5.272063  846.4267 4.670507e-11
## ENST00000589123.1 gene=NFIC           6.545195        6.181023 5.855335  531.9997 3.802396e-10
##                                   adj.P.Val    F.std
## ENST00000418210.2 gene=TMEM64  6.020801e-11 25.83580
## ENST00000555834.1 gene=RPS6KL1 2.335253e-10 23.78717
## ENST00000589123.1 gene=NFIC    1.267465e-09 21.69022

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={\text{number of coefficiets tested}}\) and \(df_2=\infty\). This can be used for downstream analysis.

Small-sample method

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

variancePartition plot

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.

Compare p-values from dream and duplicateCorrelation

In order to understand the empircal difference between dream and duplication correlation, we can plot the \(-\log_{10}\) p-values from both methods.

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.

Parallel processing

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.

There are multiple ways to use parallel processing depending on your needs

By default BPPARAM and the global setttings are set the results of bpparam(). But note that using SnowParam() can dramatically reduce the memory usage needed for parallel processing because it reduces memory redundancy between threads.

Session info

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.13-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.13-bioc/R/lib/libRlapack.so
## 
## 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       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] edgeR_3.34.0             pander_0.6.3             variancePartition_1.22.0
##  [4] Biobase_2.52.0           BiocGenerics_0.38.0      scales_1.1.1            
##  [7] BiocParallel_1.26.0      limma_3.48.0             ggplot2_3.3.3           
## [10] knitr_1.33              
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.0          tidyr_1.1.3         jsonlite_1.7.2      splines_4.1.0      
##  [5] foreach_1.5.1       gtools_3.8.2        bslib_0.2.5.1       assertthat_0.2.1   
##  [9] statmod_1.4.36      highr_0.9           yaml_2.2.1          progress_1.2.2     
## [13] numDeriv_2016.8-1.1 pillar_1.6.1        backports_1.2.1     lattice_0.20-44    
## [17] glue_1.4.2          digest_0.6.27       minqa_1.2.4         colorspace_2.0-1   
## [21] htmltools_0.5.1.1   Matrix_1.3-3        plyr_1.8.6          pkgconfig_2.0.3    
## [25] broom_0.7.6         purrr_0.3.4         snow_0.4-3          lme4_1.1-27        
## [29] tibble_3.1.2        farver_2.1.0        generics_0.1.0      ellipsis_0.3.2     
## [33] withr_2.4.2         pbkrtest_0.5.1      magrittr_2.0.1      crayon_1.4.1       
## [37] evaluate_0.14       fansi_0.4.2         doParallel_1.0.16   nlme_3.1-152       
## [41] MASS_7.3-54         gplots_3.1.1        tools_4.1.0         prettyunits_1.1.1  
## [45] hms_1.1.0           lifecycle_1.0.0     stringr_1.4.0       locfit_1.5-9.4     
## [49] munsell_0.5.0       colorRamps_2.3      compiler_4.1.0      jquerylib_0.1.4    
## [53] caTools_1.18.2      rlang_0.4.11        grid_4.1.0          nloptr_1.2.2.2     
## [57] iterators_1.0.13    labeling_0.4.2      bitops_1.0-7        rmarkdown_2.8      
## [61] boot_1.3-28         lmerTest_3.1-3      gtable_0.3.0        codetools_0.2-18   
## [65] DBI_1.1.1           reshape2_1.4.4      R6_2.5.0            dplyr_1.0.6        
## [69] utf8_1.2.1          KernSmooth_2.23-20  stringi_1.6.2       Rcpp_1.0.6         
## [73] vctrs_0.3.8         tidyselect_1.1.1    xfun_0.23

References

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.