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