Hello,

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).

https://drive.google.com/file/d/0B2X31Z4HBqmqN0UwNUhORGJWVEk/view?usp=sharing

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

DEXSeq: https://drive.google.com/file/d/0B2X31Z4HBqmqUERvZUJwTTNuWUE/view?usp=sharing

JunctionSeq/DEXSeq: https://drive.google.com/file/d/0B2X31Z4HBqmqWFBpTWs3RDYxYlE/view?usp=sharing

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.

Best,

 Lukas

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

locale:
[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      

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.