Unsuccessful DE analysis using limma - voom
1
0
Entering edit mode
loainaom • 0
@d6b8183e
Last seen 3 days ago
Israel

This might be a bit long, please bare with me. I'm conducting a differential expression analysis using limma - voom. My comparison is regarding response vs non-response to a cancer drug. However, I'm not getting any DE genes, absolute zeros. Someone here once recommended not to use contrast matrix for such a simple task, so I did that but then I get almost all 18 thousand genes as DE, so that's not good either. I don't know what is the problem here. The counts matrix and metadata matrix have the same samples with the exact same order. I filtered my counts matrix from noise genes earlier. I did everything according to this guide: https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html

<h5>Note: This problem is not only in this dataset, this is just one of many. So the problem is probably in my code.</h5>

Let me show you the code and the results I get:

## Code:

d0 = DGEList(counts)
group = metadata$Response d0$samples$group = group geneid = rownames(d0) d0$genes = geneid
d0 =  calcNormFactors(d0, method = 'TMM')
dim(d0)

## Design
design = model.matrix(~ 0 + Response, metadata) %>% as.data.frame()
colnames(design)[c(1,2)] = c('NoResponse','Response')

## Contrast
contrast = makeContrasts(NoResp_VS_Resp = NoResponse - Response, levels = colnames(design))

Voom <- voom(d0, design, plot = TRUE)
vfit <- lmFit(Voom, design = design)
vfit  <- contrasts.fit(vfit , contrasts = contrast)
efit <- eBayes(vfit)
plotSA(efit, main = 'final model: Mean-Variance trend')

summary(decideTests(efit))

deg <- topTable(efit,
coef = 'NoResp_VS_Resp',
p.value = 0.05,
number = Inf)


### And this is without contrast matrix: ( Almost 17 thousand significant genes, out of 18 thousand)

Can anyone help me figure out what is the problem ? If any additional info is needed I'll add it!

limma • 194 views
0
Entering edit mode
@gordon-smyth
Last seen 46 minutes ago
WEHI, Melbourne, Australia

Failing to get DE genes is not an error. Indeed one of the key purposes of the limma software is to avoid giving DE genes when there isn't worthwhile evidence of differences. It is very common for drug response studies with human patients to show either no differences or weak differences. That isn't the fault of the software, just the reality of biological variation between human subjects and drug responses.

You asked the same question a few weeks ago and I replied:

"The limma documentation provides several ways to plot and trouble-shoot your data. At very least you should be making the plots shown in the limma workflow paper, especially the MDS plot."

It makes no difference to the results whether you use a contrast or not, provided you do it correctly. Steve Lianoglou showed you complete code last time to get the same results with and without a contrast. If you don't want to use a contrast then you have to reform the design matrix as Steve showed you.

0
Entering edit mode

I see.. one final question, regarding gene filtering.

Could you tell me why this isn't working?

keep.exprs <- filterByExpr(d0, group = group)
d0 <- d0[keep.exprs,, keep.lib.sizes=FALSE]


Gives this error : Error in object[[a]][i, , drop = FALSE] : incorrect number of dimensions

1
Entering edit mode

BTW, with human data like yours, I would almost always use voomWithQualityWeights instead of voom in order to allow for poor quality samples or outliers. That will generally increase power to detect DE genes.

0
Entering edit mode

Thanks a ton, your help is realy appreciated!

0
Entering edit mode

That would be because you set d0$genes to be a vector instead of a data.frame, so that the genes component can't be subsetted by rows. In other words, your d0 isn't a valid DGEList object. See help("DGEList-class"). To fix this, just delete the lines geneid = rownames(d0) d0$genes = geneid


They're unnecessary and they make d0 unsubsettable.

Later on, you convert design into a data.frame. That won't cause any errors, but it is supposed to remain a matrix.