Question: JunctionSeq/DEXSeq: Dispersion plot looks strange (but is fine in DESeq2)
gravatar for l.frank
2.9 years ago by
l.frank0 wrote:


I have a single-end, stranded RNA-seq dataset consisting of two conditions (ctrl + treatment) and two biological replicates per condition and I would like to study how the treatment affects 

a) overall gene expression (--> finding differentially expressed genes)

b) splicing/exon usage (--> finding genes that show differential exon or exon-exon junction usage upon treatment)

To do so I am currently using DESeq2 for DE analysis. For analyzing changes in splicing I wanted to compare two different approaches: DEXSeq and JunctionSeq.

As far as I understood, beside using a slightly different normalisation and counting strategy, JunctionSeq is basically equally potent as DEXSeq with the addition of testing for differential splice junction usage.

Analysis with DESeq2 worked nicely.

Similarly, I overall managed to run DEXSeq and JunctionSeq on my dataset. Since DEXSeq can be run also "within" JunctionSeq, I made use of this possibility and run both analyses independently within JunctionSeq.

(Also I run DEXSeq "standalone" for comparison.)

The overall data looked ok. Also I confirmed the quality of the sequencing data by checking different quality metrics using the Java Tool "QoRTs". Quality was good for all samples, except for a slight (BUT UNIFORM) 3'bias in gene body coverage that I noticed.

Now to my actual problem. For estimating dispersion, DEXSeq uses the same methods/approach as DESeq2. Similarly, the vignette of the JunctionSeq package states that they use the approach from DEXSeq to estimate dispersion for both exons and splice junctions.

But what I observe when I generate a dispersion plot using plotDispEsts(Object) in R:

In DESeq2: The dispersion plot looks as "expected" - to say: as shown in the examples from various tutorials. A uniform cloud with few outliers and the fit line describes the distribution quite well. (See attached picture).

HOWEVER: When I do the same thing within JunctionSeq/DEXSeq

... I get a dispersion plot which I find difficult to interpret, since I have not much experience about the pitfalls/methods behind this. But somehow I find it looks a bit suspicious.

The plot looks exactly the same and equally strange for DEXSeq(standalone), DEXSeq (run within JunctionSeq) and JunctionSeq. But for DESEq2 it looked good. See attached pictures



Do you have any explanation for this? Or more precisely: can you comment on the dispersion plot - if it looks strange to you or if you can explain this behaviour.

Thank you a lot in advance.




R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X El Capitan 10.11.6

[1] C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] Biobase_2.34.0                httr_1.2.1                    AnnotationHub_2.6.4           splines_3.3.2                
 [5] Formula_1.2-1                 shiny_0.14.2                  assertthat_0.1                interactiveDisplayBase_1.12.0
 [9] stats4_3.3.2                  latticeExtra_0.6-28           RBGL_1.50.0                   BSgenome_1.42.0              
[13] Rsamtools_1.26.1              yaml_2.1.14                   RSQLite_1.1                   lattice_0.20-34              
[17] biovizBase_1.22.0             digest_0.6.10                 GenomicRanges_1.26.1          RColorBrewer_1.1-2           
[21] XVector_0.14.0                colorspace_1.3-1              ggbio_1.22.0                  htmltools_0.3.5              
[25] httpuv_1.3.3                  Matrix_1.2-7.1                plyr_1.8.4                    OrganismDbi_1.16.0           
[29] DESeq2_1.14.0                 XML_3.98-1.5                  biomaRt_2.30.0                genefilter_1.56.0            
[33] zlibbioc_1.20.0               xtable_1.8-2                  scales_0.4.1                  BiocParallel_1.8.1           
[37] htmlTable_1.7                 tibble_1.2                    annotate_1.52.0               IRanges_2.8.1                
[41] ggplot2_2.2.0                 SummarizedExperiment_1.4.0    GenomicFeatures_1.26.0        nnet_7.3-12                  
[45] BiocGenerics_0.20.0           lazyeval_0.2.0                survival_2.40-1               magrittr_1.5                 
[49] mime_0.5                      memoise_1.0.0                 GGally_1.3.0                  foreign_0.8-67               
[53] graph_1.52.0                  BiocInstaller_1.24.0          tools_3.3.2                   data.table_1.9.8             
[57] stringr_1.1.0                 S4Vectors_0.12.0              locfit_1.5-9.1                munsell_0.4.3                
[61] cluster_2.0.5                 AnnotationDbi_1.36.0          ensembldb_1.6.2               Biostrings_2.42.0            
[65] GenomeInfoDb_1.10.1           grid_3.3.2                    RCurl_1.95-4.8                dichromat_2.0-0              
[69] VariantAnnotation_1.20.2      bitops_1.0-6                  gtable_0.2.0                  DBI_0.5-1                    
[73] reshape_0.8.6                 reshape2_1.4.2                R6_2.2.0                      GenomicAlignments_1.10.0     
[77] gridExtra_2.2.1               knitr_1.15.1                  rtracklayer_1.34.1            Hmisc_4.0-0                  
[81] stringi_1.1.2                 parallel_3.3.2                Rcpp_0.12.8                   geneplotter_1.52.0           
[85] rpart_4.1-10                  acepack_1.4.1      


ADD COMMENTlink modified 2.9 years ago by Steve Lianoglou12k • written 2.9 years ago by l.frank0
Answer: JunctionSeq/DEXSeq: Dispersion plot looks strange (but is fine in DESeq2)
gravatar for Michael Love
2.9 years ago by
Michael Love26k
United States
Michael Love26k wrote:

The final dispersion estimates for JunctionSeq and DEXSeq look normal, clustered around the trend line. To be honest I'm not sure why there is the negative linear slope for the MLEs in the log dispersion over log mean scatterplot. But if I were you I would look at example plots for individual results (e.g. plotDEXSeq or the equivalent for JunctionSeq), and if those are as expected, I wouldn't think there is a problem here.

ADD COMMENTlink written 2.9 years ago by Michael Love26k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 391 users visited in the last hour