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