inst/doc/dream.R

In variancePartition: Quantify and interpret divers of variation in multilevel gene expression experiments

## ----setup, echo=FALSE, results="hide"------------------------------------------------------------
knitr::opts_chunk$set(tidy=FALSE, cache=TRUE,
                      dev="png",
                      package.startup.message = FALSE,
                      message=FALSE, error=FALSE, warning=TRUE)
options(width=100)

## ----kable, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'-----------------------------
library('pander')
tab = rbind(c('precision weights to model measurement error in RNA-seq counts', "limma/voom", "@Law2014"),
  c('ability to model multiple random effects', 'lme4', '@Bates2015'),
  c('random effects estimated separately for each gene', 'variancePartition', '@hoffman2016'),
  c('hypothesis testing for fixed effects in linear mixed models', 'lmerTest', 'Kuznetsova, et al. -@kuznetsova2017'),
  c('small sample size hypothesis test', 'pbkrtest' , 'Halekoh, et al. -@halekoh2014'), 
  # c('emprical Bayes moderated t-test', 'limma/eBayes', '@smyth2004'),
  c('','','')
  )
colnames(tab) = c("Model property", "Package", "Reference")

panderOptions('table.split.table', Inf)
panderOptions('table.alignment.default', 'left')

pander(tab, style = 'rmarkdown')

## ----preprocess, eval=TRUE, results='hide'--------------------------------------------------------
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
geneExpr = DGEList( countMatrix[isexpr,] )
geneExpr = calcNormFactors( geneExpr )

# make this vignette faster by analyzing a subset of genes
geneExpr = geneExpr[1:1000,]

## ----dupCor, eval=TRUE----------------------------------------------------------------------------
# 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 )

## ----lmm, eval=TRUE, message=FALSE, fig.height=2, results='hide'----------------------------------
# Specify parallel processing parameters
# this is used implicitly by dream() to run in parallel
param = SnowParam(4, "SOCK", progressbar=TRUE)
register(param)

# The variable to be tested must be a fixed effect
form <- ~ Disease + (1|Individual) 

# estimate weights using linear mixed model of dream
vobjDream = voomWithDreamWeights( geneExpr, form, metadata )

# Fit the dream model on each gene
# By default, uses the Satterthwaite approximation for the hypothesis test
fitmm = dream( vobjDream, form, metadata )

## ----lmm.print------------------------------------------------------------------------------------
# Examine design matrix
head(fitmm$design, 3)

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

## ----contrast, eval=TRUE, fig.width=7, fig.height=2-----------------------------------------------
form <- ~ 0 + DiseaseSubtype + Sex + (1|Individual)  
L = getContrast( vobjDream, form, metadata, c("DiseaseSubtype2", "DiseaseSubtype1"))

# Visualize contrast matrix
plotContrasts(L) 

## ----fit.contrast---------------------------------------------------------------------------------
# fit dream model with contrasts
fit = dream( vobjDream, form, metadata, L)

# get names of available coefficients and contrasts for testing
colnames(fit)

# extract results from first contrast
topTable( fit, coef="L1", number=3 )

## ----contrast.combine, eval=TRUE, fig.width=7, fig.height=3---------------------------------------
form <- ~ 0 + DiseaseSubtype + Sex + (1|Individual) 

# define and then cbind contrasts
L1 = getContrast( vobjDream, form, metadata, c("DiseaseSubtype2", "DiseaseSubtype1"))
L2 = getContrast( vobjDream, form, metadata, c("DiseaseSubtype1", "DiseaseSubtype0"))
L = cbind(L1, L2)     

# Visualize contrasts
plotContrasts(L) 

# fit both contrasts
fit = dream( vobjDream, form, metadata, L)

# extract results from first contrast
topTable( fit, coef="L1", number=3 )

## ----maual.contrasts, fig.width=8, fig.height=4---------------------------------------------------
# the tests is DiseaseSubtype0 - (DiseaseSubtype1/2 + DiseaseSubtype2/2)
# Note that the order of the coefficients must be the same as from getContrast()
L3 = c(1, -1/2, -1/2, 0)

# combine L1 and L2 contrasts with manually defind L3 contrast.
Lall = cbind(L, data.frame(L3 = L3))

# Note that each contrast must sum to 0
plotContrasts(Lall)

# fit dream model to evaluate contrasts
fit = dream( vobjDream[1:10,], form, metadata, L=Lall)

topTable(fit, coef="L3", number=3)

## ----joint.test, fig.height=3, message=FALSE------------------------------------------------------
# extract results from first contrast
topTable( fit, coef=c("DiseaseSubtype2", "DiseaseSubtype1"), number=3 )

## ----kr, eval=FALSE-------------------------------------------------------------------------------
#  fitmmKR = dream( vobjDream, form, metadata, ddf="Kenward-Roger")

## ----vp-------------------------------------------------------------------------------------------
# Note: this could be run with either vobj from voom()
# or vobjDream from voomWithDreamWeights()
# The resuylts are similar
form = ~ (1|Individual) + (1|Disease)
vp = fitExtractVarPartModel( vobj, form, metadata)

plotVarPart( sortCols(vp))

## ----define---------------------------------------------------------------------------------------
# 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)

## ----parallel, eval=FALSE-------------------------------------------------------------------------
#  # 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)
#  
#  # Or disable parallel processing and just do serial
#  param = SerialParam()
#  fitmm = dream( vobjDream, form, metadata, BPPARAM=param)

## ----parallel2, eval=FALSE------------------------------------------------------------------------
#  param = SnowParam(4, "SOCK", progressbar=TRUE)
#  register(param)
#  
#  # uses global parallel processing settings
#  fitmm = dream( vobjDream, form, metadata )

## ----sessionInfo----------------------------------------------------------------------------------
sessionInfo()