[BioC] limma: tissue specific genes using voom

[Gordon K Smyth]

Dear Jaaved,

limma and edgeR work very similarly.  If you use coef=2 in edgeR, then 
that's just what you use in limma also.

You only have 17 tissues, so obviously trying to refer to coefficient 
number 18 in limma will cause an error.  There is no tissue 18.

Best wishes
Gordon

> Date: Tue, 25 Jun 2013 13:18:31 -0400
> From: Jaaved Mohammed <jm889 at cornell.edu>
> To: <bioconductor at stat.math.ethz.ch>
> Subject: [BioC] limma: tissue specific genes using voom
>
> Hello everyone,
>
> I'm trying to compare limma vs edgeR in finding tissue specific genes
> from my RNAseq samples. For edgeR I've followed directions from a
> previous post and have managed to get good results. The thread is named
> "edgeR: finding tissue specific genes" and here is the link:
> http://permalink.gmane.org/gmane.science.biology.informatics.conductor/47950
>
> However, I'm having a difficult time creating the contrast for the last
> tissue in limma. For the last tissue, the contrast is a vector with 18
> elements which lmFit does not like.
>
> First of all, here is a snippet of code of what I did in edgeR
>
> contrasts(tiss_groups) <- contr.sum(tiss_groups)
> design <- model.matrix(~tiss_groups)
> y <- estimateGLMCommonDisp(y, design)
> y <- estimateGLMTrendedDisp(y, design)
> fit <- glmFit(y, design)
>
> ## Get enriched genes for 1st tissue:
> ql <- glmQLFTest(fit, coef=2)
> top1 <- topTags(ql)
>
> ## Getting tissue 2 - 17 is the same as above but with different coef
> values ###
>
> #Get enrichment for the last tissue
> cont <- rep(-1,18)
> cont[1] <- 0
> ql <- glmQLFTest(fit, contrast=cont)
> top18 <- topTags(de)
>
>
>
>
> Here is what I'm doing in limma:
>
> contrasts(tiss_groups) <- contr.sum(tiss_groups)
> design <- model.matrix(~tiss_groups)
>
> v <- voom(y, design, plot=T)
> fit <- lmFit(v,design)
> fit <- eBayes(fit)
>
> ##Get top genes in tissue j = 1:17
> top1 <- topTable(fit, coef=(j+1), adjust="BH", sort="p")
>
> ## How to get top genes the last tissue? Here is what I did:
> cont <- rep(-1, 18)
> cont[1] <- 0
> fit <- lmFit(v,design=cont)   ##<------- this fails with error
> "(subscript) logical subscript too long"
> fit <- eBayes(fit)
> top18 <- topTable(fit, n=nrow(enrich), sort="p")
>
>
>
> Thanks for your help,
> Jaaved Mohammed
> Graduate student of Computational Biology
> Cornell University


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