## ----setup, echo=FALSE, results="hide"------------------------------------------------------------
knitr::opts_chunk$set(tidy=FALSE, cache=TRUE,
                      dev="png",
                      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------------------------------------------------------------------------
library('variancePartition')
library('edgeR')
library('doParallel')
data(varPartDEdata)

isexpr = rowSums(cpm(countMatrix)>0.1) >= 3

# Standard usage of limma/voom
genes = DGEList( countMatrix )
genes = calcNormFactors( genes )
design = model.matrix( ~ Disease, metadata)
vobj_tmp = voom( genes, design, plot=FALSE)

# apply duplicateCorrelation 
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( genes, design, plot=FALSE, block=metadata$Individual, correlation=dupcor$consensus)

## ----dupCor, eval=TRUE----------------------------------------------------------------------------

design = model.matrix( ~ Disease, metadata)

# 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----------------------------------------------------------------
cl <- makeCluster(4)
registerDoParallel(cl)

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

# Get the contrast matrix for the hypothesis test
L = getContrast( vobj, form, metadata, "Disease1")

# Fit the dream model on each gene
# Apply the contrast matrix L for the hypothesis test  
# By default, uses the Satterthwaite approximation for the hypothesis test
fitmm = dream( vobj, form, metadata, L)

# Fit Empirical Bayes for moderated t-statistics
fitmm = eBayes( fitmm )

# get results
topTable( fitmm )

## ----kr, eval=FALSE-------------------------------------------------------------------------------
#  fitmmKR = dream( vobj[1:10,], form, metadata, L, ddf="Kenward-Roger")

## ----contrast, eval=TRUE--------------------------------------------------------------------------
form <- ~ 0 + Disease + (1|Individual) 
L = getContrast( vobj, form, metadata, c("Disease1", "Disease0"))
L

## ----contrast.combine-----------------------------------------------------------------------------
form <- ~ 0 + Disease + (1|Individual) 

# define and then cbind contrasts
L1 = getContrast( vobj, form, metadata, "Disease0")
L2 = getContrast( vobj, form, metadata, "Disease1")
L = cbind(L1, L2)

# fit both contrasts
fit = dream( vobj[1:10,], form, metadata, L)

# empirical Bayes step
fiteb = eBayes( fit )

# extract results from first contrast
topTable( fiteb, coef="L1" )

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

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