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