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.

Dear Brad, As noted in the help file for voom, voom can take any of the arguments that are passed onto lmFit, including a blocking variable and correlation. Using the limma and edgeR packages, an analysis on read counts using voom and a blocking variable could look something like this: > nf <- calcNormFactors(Counts) > design <- model.matrix(~Group) > y <- voom(Counts,design,lib.size=colSums(Counts)*nf) > corfit <- duplicateCorrelation(y,design,block=Subject) > y <- voom(Counts,design,plot=TRUE,lib.size=colSums(Counts)*nf,block= Subject,correlation=corfit$consensus) > fit <- lmFit(y,design,block=Subject,correlation=corfit$consensus) > fit <- eBayes(fit) Voom is run twice - the first time to obtain log-counts-per-million values and observational level weights to feed into duplicateCorrelation(); and a second time with the blocking variable and estimated correlation. I hope this helps! Regards, Charity > > Date: Tue, 28 May 2013 20:48:25 -0700 (PDT) > From: "Brad Rosenberg [guest]" <guest@bioconductor.org> > To: bioconductor@r-project.org, rosenbb@rockefeller.edu > Subject: [BioC] Block design with voom and limma for multi-level > experiment > > 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 ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:8}}