Thanks for the updates. I have a feeling that the best answer you'll get is from Alejandro Reyes by posting to the bioconductor support site. One thing I should note is that I expect you can get weird results like this for genes that are also differentially expressed, as is the case in your example. The ~sample part of the model essentially corrects for group mean differences, so you then will get some weird results and fold-change directions if a gene is differentially expressed (e.g., that's probably why the "exon usage" plot looks like a shifted version of the normalized counts/expression...since the sample means are then more similar).
Now having said that, looking into the code for plotDEXSeq() is telling. The main difference is what happens if we specify splicing=T versus expression=T. If you specify expression=T, you end up invoking:
coeff <- as.matrix( t(getEffectsForPlotting(es, averageOutExpression=FALSE, groupingVar=fitExpToVar) )[featuresInGene,] )
whereas with splicing=T, you get:
coeff <- as.matrix( t( getEffectsForPlotting(es, averageOutExpression=TRUE, groupingVar=fitExpToVar) )[featuresInGene,] )
You can see that the only difference is the averageOutExpression parameter, which multiplies the fitted values by their means (well, it adds the mean, but since this is on the log2 scale this ends up being a scaling operation). Now if you're using normalized values to begin with, it would seem to me that doing that is just going to produce funky results.
Edit: I should add that this really emphasizes the importance of producing and looking at plots like this. Clearly the statistical results in this case are non-sense due to the differential expression of that gene...even though there's also legitimate differential exon usage going on. A better model might be to try to account for the sample nesting with groups, but that's a whole other can of worms and I suspect has a lot of issues of its own (as is often the case with mixed-effects models).
Edit2: Actually, the more I think about it, the "better model" I mentioned above wouldn't be better. Ignore that.
Edit3: To more directly answer your question: differential expression shows the fit normalized counts, whereas differential usage shows that after accounting for differences due to each fit coefficient (specifically, those interacting with exon). The former is closer to the raw data but the latter will better match the p-values and fold-changes. I would argue that it's best to compare the statistical results to the expression values (or just the normalized counts), for the reasons demonstrated in this example.
•
written
by
andrew.j.skelton73
could you put your code and your sampleTable please
Edited the post
Can you give an example of exons that seem to be differentially used but not differentially expressed? The point of differential use is to control for gene-level expression differences between samples that aren't then relevant to isoform usage/splicing.
I've added an example image of what I mean, if that helps!
N.B., I asked Alejandro to have a look at this question (finally found him on twitter thanks to Mike Love).