Can DEXSeq work for a single gene?

[Qiuyue]

Hi, folks!
Recently I'm using DEXSeq to perform differential exon usage.

My situation is a little different from the example of DEXSeq, that is, each gene of my design has a different two-class condition, so I first estimate sizeFactors for all samples across all genes. When testing for each gene, I give it the same sizeFactor as calculated above, because sizeFactor just correlated with sample. And then I estimate dispersions for each gene, and testForDEU gene by gene. Am I in the right way? Can estimateDispersions work for a single gene?

And by this way, some genes' I(1/means[good]) in the function fitDispersionFunction are negative, then the code will stop.

Another important point, the results for each exon of the same gene has a totally equal pvalue, and doesn't change after BH adjustment. How does this happen? My models just followed DEXSeq examples.

Appreciate any help from you guys! Thanks!

[Qiuyue]

Anybody help me?

[Simon Anders]

Different two-class condition for each gene? Sounds like you are trying to do some eQTL-like analysis, where you regress each exon or gene onto a nearby SNP, and because this changes from gene to gene, you need to keep changing the condition factor. Have I guessed right?

As you noticed, calling estimateDispersions for a single gene does not work. This is because we need information from all the other genes to fit a general dependence of dispersion on mean.

Your use case is not that uncommon, and it may be worth adding some functionality to DEXSeq to accommodate it. We actually have a number of preliminary ideas.

[Qiuyue]

Hi, Simon!
Thanks very much for your notice. You are totally right. We are trying to regress each exon or gene onto SNPs to map variants controlling differential exon usage. So for this topic, what's the role of NB distribution?
For another point, as DEXSeq work good for exon usage, but geneCountTable just sum up each exon's read count. Then junction reads would be counted many times. You also recommended to use geneCountTable with care, so do you know some good methods to get gene count table? Like htseq-count?
Anyway, I'm fond of DEXSeq and looking forward to seeing your good news.

[areyes]

HI Qiuyue - The double counting of reads in DEXSeq is true!, it is better to count reads on genes separately... for counting read numbers per gene, htseq-count is precisely what you need!