[BioC] Block design with voom and limma for multi-level experiment

[Brad Rosenberg [guest]]

Hello all --

I'm attempting to implement the voom -> limma strategy for the analysis of a multilevel RNA-Seq experiment with a blocking setup, and I have a question about how voom fits in to the workflow. More specifically, for a paired experimental design that uses duplicateCorrelation() to handle blocking on the same subject, how is the appropriate experimental design (taking subject blocking into account) fed to voom?

An easily accessible example is in the limma users guide, section 8.7, "Multilevel experiments." Given the target frame provided:

   Subject Condition Tissue
1        1  Diseased      A
2        1  Diseased      B
3        2  Diseased      A
4        2  Diseased      B
5        3  Diseased      A
6        3  Diseased      B
7        4    Normal      A
8        4    Normal      B
9        5    Normal      A
10       5    Normal      B
11       6    Normal      A
12       6    Normal      B

The design referenced in the manual joins condition and tissue:
> Treat <-factor(paste(targets$Condition,targets$Tissue,sep="."))
> design <- model.matrix(~0+Treat)
> colnames(design) <- levels(Treat)

And then estimates subject correlation by setting a block by Subject:
> corfit <- duplicateCorrelation(eset,design,block=targets$Subject)
> corfit$consensus

The blocking is then input with the design into the fit:
> fit <- lmFit(eset,design,block=targets$Subject,correlation=corfit$consensus)

In order to translate this example to an RNA-Seq experiment, I would like to use voom prior to fitting. Given that voom takes "design" as an argument, but in this example the experimental design is divided amongst the design matrix and the blocking (by subject), what is the appropriate way to run voom such that it takes all of the necessary design information for its conversions? Is re-designing in a nested format the only solution? Or is there a way to maintain the blocking workflow and still use voom appropriately?

Thanks very much for your help.

Best,
Brad Rosenberg
The Rockefeller University

 -- output of sessionInfo(): 

> sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] annotate_1.36.0      AnnotationDbi_1.20.3 biomaRt_2.14.0      
 [4] edgeR_3.0.8          limma_3.14.4         DESeq_1.10.1        
 [7] lattice_0.20-10      locfit_1.5-9         Biobase_2.18.0      
[10] BiocGenerics_0.4.0  

loaded via a namespace (and not attached):
 [1] DBI_0.2-5          genefilter_1.40.0  geneplotter_1.36.0 grid_2.15.2       
 [5] IRanges_1.16.4     parallel_2.15.2    RColorBrewer_1.0-5 RCurl_1.95-3      
 [9] RSQLite_0.11.2     splines_2.15.2     stats4_2.15.2      survival_2.36-14  
[13] tools_2.15.2       XML_3.95-0.1       xtable_1.7-1     

--
Sent via the guest posting facility at bioconductor.org.