[BioC] limma: tissue specific genes using voom

[Jaaved Mohammed]

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