DEXSeq gene level counts

[Julien Roux]

Dear DEXSeq authors,
I have a question related to the following paragraph in the DEXSeq vignette:

Note that a read that mapped to several exons of a gene is counted for each of these exons by the dexseq_count.py script. The table returned geneCountTable will hence count the read several time for the gene, which may skew downstream analyses in subtle ways. Hence, we recommend to use geneCountTable with care and use more sophisticated counting schemes where appropriate.

I am wondering if this problem is also affecting the tests for differential expression of exons in DEXSeq? I have the feeling that the reads overlapping multiple exons are also counted several times when it comes to estimate gene expression levels (to control for overall differences in expression at the gene level when testing exons). Is this accounted for in the package functions?

Thanks for your help
Julien

[Simon Anders]

It's not really an issue because if a gene is expressed more strongle in one sample than in another, all exon counts, including those that result from double-counting due to overlaps, will be scaled up in a proportional manner, so that the expression ratio between the samples is still estimated correctly.

[Julien Roux]

Hi Simon,
Thanks for your answer. However I am not fully convinced that it is not an issue...

In my dataset at least 40% of the reads overlap an exon junction. If you count these read twice when estimating gene baseline levels, it can lead to major estimation errors for the expression ratio between conditions (especially when there is a large difference between conditions).

An easy fix would be to provide DEXSeq the correct gene count table. Is there a way to do this and use it in the function "testGeneForDEU"?
Thanks for your help
Julien

[Simon Anders]

As DEXSeq performs a hypothesis test, we only need to ensure that the "null" case is not affected. The null hypothesis in DEXSeq is that the relative proportions that the differnt isoforms contribute to overall gene expression does not depend on the condition, while the overall expression strength may change.

If now a gene is expressed twice as strongly in one condition than in the other, all per-exon counts go up by a factor of two. The number of reads that fall onto two exons goes up by two, as well, so that the ratio of each per-exon count to the sum of all per-exon counts stays as constant. This is all that we need. There is no need for the sum of all per-exon counts to reflect gene expression strength faithfully.