[BioC] Gene expression analysis with edgeR with a large, nested design matrix

[Gordon K Smyth]

Dear Uli,

edgeR is probably the fastest of the glm negative binomial packages, as we 
have done a lot of work moving all the fitting code to the C++ level. 
Nevertheless, the ususal glm pipeline will become very slow when there are 
over 3000 samples and the design matrix has 1500 columns.

Here are two other possibilities.  First you could switch to voom instead. 
I ran voom on your simulated data example in a few minutes on my laptop 
computer:

   y <- DGEList(counts = reads)
   y <- calcNormFactors(y)
   v <- voom(y,expdesign,plot=TRUE)
   fit <- lmFit(v,expdesign)
   fit <- eBayes(fit)
   topTable(fit, coef=5)

and so on.  Alternatively, you would stick to edgeR and use the 
quasi-likelihood pipeline:

   y <- DGEList(counts = reads)
   y <- calcNormFactors(y)
   y <- estimateGLMCommonDisp(y, design=expdesign, method="Pearson")
   ql <- glmQLFTest(y,expdesign)
   topTags(ql,coef=5)

etc.  The glmQLFTest() function will take more time than voom but less 
time than estimateGLMCcommonDispersion on your data.  The common 
dispersion step is optional in this pipeline -- I include it because 
you've already done it.

Either of these approaches have good statistical motivations.  See the 
help pages for voom and glmQLFTest for references to published papers 
describing them.

Best wishes
Gordon

> Date: Fri, 23 May 2014 09:06:21 -0700 (PDT)
> From: "Uli Braunschweig [guest]" <guest at bioconductor.org>
> To: bioconductor at r-project.org, u.braunschweig at utoronto.ca
> Subject: [BioC] Gene expression analysis with edgeR with a large,
> 	nested design matrix
>
> Hi All,
>
> My problem is the following:
> I have expression counts for 50 genes in an RNAi screen with 1,536 treatments (which includes positive and negative controls, so really 1,416 unique treatments) in two replicates, done in 96-well format (2x16 plates). I know that plate effects and edge effects (whether a well was located on the edge of a plate) are significant, so the design should include treatment, plate, and location (edge or interior). Locations are identical between replicates. Each plate has two negative controls ("siNT"), as well as other controls.
>
> I am only interested in the contrasts of each of the treatments vs. the "siNT" control. I thought that the model of edgeR would be useful to score significant hits while at the same time dealing with the mentioned technical biases in a meaningful manner. However, I've had to kill the analysis because estimating the trended and tag-wise dispersions takes excessively long.
>
> My question is: Is it even feasible to try to adress a problem with a design that has so few genes and so many treatments using edgeR (or DESeq2)?
>
> library(edgeR)
>
> ## Data look similar to this:
> reads <- matrix(round(2000 * rexp(50 * 3072)), nrow=50)  # dense matrix of 50 genes x (1536 treatments in duplicate)
>
> ## Here is how I create my design factors:
> # (I've left out 'replicate' because it is a linear combination of plates 1-16 and 17-32)
> rows  <- rep(rep(1:8, each=12), 32)
> cols  <- rep(rep(1:12, 8), 32)
> plate <- rep(1:32, each=96)
>
> type.nt <-      rows == 1 & cols == 1 |   # the negative control to compare everything to; 2 per plate
>                rows == 4 & cols == 9
> type.posCtl <-  rows == 4 & cols == 3 |   # positive control; 2 per plate
>                rows == 5 & cols == 9
> type.mock <-    rows == 3 & cols == 3 |   # another control; 2 per plate
>                rows == 8 & cols == 12
> type.empty <-   plate %in% c(16, 32) & (  # another control; a bunch on only 2 plates
>                rows %in% c(2:3,6:8) & cols == 9 |
>                cols %in% c(10,11) |
>                rows %in% 1:7 & cols == 12
>                )
> type.edge <- rows %in% c(1,8) | cols %in% c(1,12)  # position on the plate
>
> treat <- rep(paste("T", 1:1536, sep=""), 2)  # treatments
> treat[type.nt]     <- "siNT"
> treat[type.mock]   <- "mock"
> treat[type.empty]  <- "empty"
> treat[type.posCtl] <- "siPosCtl"
> treatfac <- relevel(as.factor(treat), ref="siNT")
>
> edgefac  <- as.factor(ifelse(type.edge, yes="edge", no="interior"))
> platefac <- as.factor(paste("P", plate, sep=""))
>
> expfact <- data.frame(treatment = treatfac,
>                      platepos  = edgefac,
>                      plate     = platefac
>                      )
>
> expdesign <- model.matrix(formula(~ treatment + plate + platepos), data=expfact)
>
> ## Estimating the dispersions
> y <- DGEList(counts = reads)  # reads is the 50x3072 matrix
> y <- calcNormFactors(y)
> y <- estimateGLMCommonDisp(y, design=expdesign, method="Pearson", verbose=TRUE)  # faster than Cox-Reid and probably ok since there are many treatments
> y <- estimateGLMTrendedDisp(y, design=expdesign)
> y <- estimateGLMTagwiseDisp(y, design=expdesign)
> ## Neither of the last two steps finish running in a day; same for estimateGLMCommonDisp() if method="CoxReid"
>
> I was then hoping to extract the contrasts of each treatment against the "siNT" control.
> Would it make sense to combine the two technical factors, or subset the count and design matrices for each treatment in a way that reduces the number of treatments, and run them separately? Alternatively, I thought of doing the analysis separately for each treatment using the whole count matrix but amalgamating all other non-control treatments in an "other" group". This seems feasible when run in parallel, but it would be overkill...
> Any suggestions on how to proceed?
>
> Kind regards,
> Uli
>
> -- 
> Ulrich Braunschweig, PhD
>
> The Donnelly Centre
> University of Toronto
> 160 College Street, Room 1030
> Toronto, Ontario
> Canada M5S 3E1
>
> u.braunschweig at utoronto.ca
>
> -- output of sessionInfo():
>
> R version 3.1.0 (2014-04-10)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C
> [3] LC_TIME=en_IE.UTF-8        LC_COLLATE=en_CA.UTF-8
> [5] LC_MONETARY=en_IE.UTF-8    LC_MESSAGES=en_CA.UTF-8
> [7] LC_PAPER=en_IE.UTF-8       LC_NAME=C
> [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_IE.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] edgeR_3.6.2  limma_3.20.4

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}