limma+voom flips the sign of the effect from CPM?

Question

[Jautis]

Hello, I'm having a hard time figuring out how limma can switch the sign of an effect compared to the raw cpm or log2(cpm) values. In the figure below, I'm showing 4 genes (2x2 matrix) in our data compared between two conditions (called infected and control for simplicity). On the top row, cpm and voom normalized data agree as I would expect. However, on the bottom, there appear to be pretty substantial differences.

Thanks for any insight you can add into this! I've already checked all the usual suspects like whether any of the labels are different, and I haven't found an explanation as to why an effect would switch sign (I get how an effect may become more/less evidence). Thanks!

[Gordon Smyth]

Cross-posted to Bioconductor https://support.bioconductor.org/p/9139613/

Answer 1

[Gordon Smyth]

It seems to me that the differences that you show are not possible, or at least cannot be due to voom. However it is hard to tell without any code. Rather than just occasional flipping of signs, your plots appear to show a global shift of all the infected cpms vs the control cpms for all the genes. The infected vs control fold change is increased vs the plots on the left for all four genes. If you believe this is a voom phenomenon then you are welcome to send me by email a count matrix and code illustrating the phenomenon.

Note that you cannot illustrate the effects of voom with plots like this. voom produces precision weights and plots like this do not take account of those weights.

It is possible to flip logFC signs when the counts are very small and the library sizes are systematically different between groups. However such flips are never statistically significant (at least I've never seen such a case) and I don't think that is what is happening in your example.

PS. If you want in-depth advice in the future with voom or limma you might consider posting to the Bioconductor support site. I have answered this question but I and the other limma authors can't guarantee to answer questions on Biostars.

[Jautis]

That's what I thought as well, but it's happening none the less which is why I'm confused. I don't think this is a consequence of small read depths because I have more 140-4k reads per sample for this gene and library sizes are relatively consistent.

I know that v$E here doesn't correspond directly to the lmFit results and normally I wouldn't plot it. That's how I discovered this was happening when my results from lmFit disagreed in sign with the log2(cpm) results that I was plotting for them.

Thanks for your offer to take a look! I'm e-mailing you a concise + clean copy of the code and the data needed to reproduce the effect.

[Gordon Smyth]

When you ask a questions on Biostars it is important to show the code you used because that allows us to figure out what the problem is.

As soon as you posted a few lines of code to Bioconductor https://support.bioconductor.org/p/9139613/ it became obvious that you have applied TMM normalization to the right hand side plots but not to the left hand side. The differences between the left and right hand sizes therefore have nothing to do with limma or voom.

If you look closely at your data, I think you will find that the TMM normalization factors for the infected samples are systematically smaller than those for the controls.

Let me add that TMM normalization usually is needed for RNA-seq data. Just plotting plain cpm values as in your left hand plots is too simplistic and will often give misleading results.