fitted expression in DEXSeq

[xguo]

Hi, there,
We have been using DEXSeq for splicing analysis. One thing we found puzzling is that the log2fold value derived from estimatelog2FoldChanges function is often not consistent with what can be inferred from fitted expression values shown by plotDEXSeq function. Often times, we can see clear upregulation of fitted expression value based on graph, while the log2fold value is negative. Does anyone know what's behind the log2fold and fitted expression value shown in the plotDEXSeq-derived graph?

thanks for any insight.

[Wolfgang Huber]

Dear Xguo,

thank you for the feedback. Can you please post here an example plot and the output of estimatelog2FoldChanges that you are refering to?

Best wishes
Wolfgang

[xguo]

In the following example, exon2 and exon11 have padjust less than 0.05. Exon 11 has log2fold of 31, which seems to be consistent with the graph. However, exon2 has decreased expression in tumor compared to normal (log2fold of -0.554), while this exon shows increased expression in the graph based on either fitted expression or normalized count.

GeneID exonID dispersion pvalue padjust meanBase log2fold(tumor/normal)
ENSG00000003987 E001 0.558 7.927e-01 1.00 2.108 0.480
ENSG00000003987 E002 0.071 9.037e-06 0.0067 75.078 -0.554
ENSG00000003987 E003 0.674 2.025e01 1.00 2.877 -1.379
ENSG00000003987 E004 0.199 7.090e-01 1.00 6.813 0.429
ENSG00000003987 E005 0.215 9.656e-01 1.00 6.096 0.380
ENSG00000003987 E006 0.295 1.591e-01 1.00 4.041 1.144
ENSG00000003987 E007 0.252 4.100e-02 1.00 4.940 2.242
ENSG00000003987 E008 0.196 1.652e-01 1.00 6.928 1.506
ENSG00000003987 E009 0.378 3.422e-02 1.00 2.998 2.127
ENSG00000003987 E010 0.235 3.124e-01 1.00 5.406 0.935
ENSG00000003987 E011 0.277 7.443e-05 0.033 4.382 31.106
ENSG00000003987 E012 0.555 3.291e-02 1.00 1.930 2.528
ENSG00000003987 E013 1.906 NA NA 0.519 -4.182
ENSG00000003987 E014 3.668 NA NA 0.265 -5.019
ENSG00000003987 E015 6.398 NA NA 0.151 -2.490
ENSG00000003987 E016 1.858 NA NA 0.533 -3.440
ENSG00000003987 E017 2.264 NA NA 0.435 -2.565
ENSG00000003987 E018 1.975 NA NA 0.500 -4.998
ENSG00000003987 E019 1.460 1.928e-02 1.00 0.861 -3.435

[Simon Anders]

Hi

the fold change for exon 11 may seem large, but this is to be expected as the expression of this exon in normal is essentially zero.

For exon 2, the normalized counts for tumor are larger than in normal, but this is the case for many exons. Hence, DEXSeq concludes that the gene has stronger overall expression in tumor than in normal and adjusts for this. If you use the plotDEXSeq function with option 'splicing=TRUE', you get the fitted values with the effect of overall expression averaged out. Then, the exon should be stronger in normal. (The log2 fold change in the result table corresponds to the difference depicted with splicing=TRUE.)

This gene seems to have a quite complex splicing pattern. In such cases, DEXSeq's approach helps to find the genes that have something interesting going on, but it does not do a good job in pinpointing which specific exons are affected. Here, we recommend to use the p values only as an indication that the gene as a whole is affected by differential exon usage, and rely on the splicing plot for manual interpretation.

[xguo]

Thanks a lot for the clarification. I'm clear about it now.

Another quick question about sorting for paired-end alignment. In DEXSeq manual, sort -k1,1 -k2,2n ... is used to sort paired-end alignment before generating exon count file. I'm wondering if samtools sort -n, or picard SortSam SortOrder=queryname can be used instead. Is the sorting on second field critical?

[Simon Anders]

No. Sorting with 'samtools sort -n' should be fine and is what I usually do. I admit that I don't remember why I wrote into the manual what you have quoted. I'm afraid I just copied and pasted it from some old documentation that I wrote for HTSeq long ago, when samtools sort did not yet support the '-n' option.