This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Best contrast strategy to identify condition-specific effects in limma

Hi, I’m analyzing an RNA-seq dataset with three conditions (C, D, and E) using a limma-voom pipeline, and I’m trying to understand whether condition C behaves differently from the other two.

I’m considering two possible approaches.

One is to run pairwise contrasts (C vs D and C vs E) and then take the intersection of differentially expressed genes, possibly requiring the same direction of logFC, to define a “C-specific signature”.

The other is to define a combined contrast such as C -(D + E)/2, which directly tests whether C differs from the average effect of the other two conditions.

From a statistical and biological interpretation standpoint, which approach is more appropriate for identifying C-specific effects?. Thanks!

rna differential expression drugs limma
• 493 views • link
written by fnepo • fnepo

1 answer

To identify genes that behave differently in condition C to the other two conditions, the second approach is more robust to threshold effects. If you take the intersection of C vs D and C vs E, you will exclude genes that had a p.adj of 0.051 in C vs D, but 0.049 in C vs E, while including genes that have a p.adj of 1e-10 in C vs D, but 0.049 in C vs E.

But the actaully way I would usually do this, is to do an Anova-like LRT test of ~condition vs ~1, and identify genes that vary in any way between conditions. I would then cluster those genes in to sets of genes that have similar behavoir, and look at those clusters on a heatmap to identify the cluster(s) that are different in C to D/E.

• 0 views • link
written by i.sudbery i.sudbery

That makes lot of sense! Thank you very much for the clear explanation :)

• 0 views • link
written by fnepo • fnepo

In limma,the F-test is done simply by testing both C vs D and C vs E contrasts at once, and that performs an actual ANOVA F-test instead of a anova-like LRT (which ignores the uncertainty of estimation of the dispersion parameter).

• 0 views • link
written by Gordon Smyth Gordon Smyth

Using limma with a design including condition and batch (~ batch + condition), I performed an F-test (topTableF) to identify genes varying across conditions while adjusting for batch effects. However, I obtained ~11,000 significant genes, which is essentially the total number of genes after filtering. I assume this may be due to high statistical power relative to variability (I have 4 samples per condition), leading to many small but statistically significant differences?

• 0 views • link
written by fnepo • fnepo

No, it is because you have conducted a test for the entire linear model, including the batch effect and the intercept, instead of testing for DE between conditions. Users are not supposed to use topTableF. If you used topTableF, you would have received the following message:

topTableF is obsolete and will be removed in a future version of limma. Please considering using topTable instead.

Just use topTable with the relevant coefficients. The help page tells you:

If coef takes two or more values, the moderated F-statistics for that set of coefficients or contrasts are used.

• 0 views • link
written by Gordon Smyth Gordon Smyth

To do an F-test in limma you could use design <- model.matrix(~condition+batch) and then later use topTable(fit, coef=2:3). That will give F-tests between the conditions.

However, the idea of clustering the F-test DE genes will give better results in situations where there are very large numbers of DE genes. I would not use such an approach in subtle situations with only moderate numbers of DE genes. I don't particularly recommend F-tests in situations where you want focus on one condition. When you want to focus on one condition, individual t-tests are generally more powerful and relevant.

• 0 views • link
written by Gordon Smyth Gordon Smyth

Understood! Thank you so much for your help! Keep up the good work :)

• 0 views • link
written by fnepo • fnepo

I performed an ANOVA-like likelihood ratio test by comparing a reduced model (~ batch) to a full model (~ batch + condition). In this case, I obtained around ~150 significant genes. Additionally, I corrected the expression values for batch effects using removeBatchEffect() for visualization and generated heatmaps based on the selected genes. However, I do not observe a clear cluster separating canagliflozin from the other two treatments.

• 0 views • link
written by fnepo • fnepo

Log in to answer this question.

More posts like this