The DEXSeq vignette describes several protocols for obtaining exon counts. One option is to use the
dexseq_count.py script (powered by htseq) that comes with the package. Another option is to use
summarizeOverlaps directly in R. I tried both methods with unstranded parameters and the results were identical (based on a few genes). However, stranded settings yield very different results. Some exons are off by a few reads, but some can be as low as half. I used the settings as specified in the vignette, which I assumed would be most appropriate. I also tried adjusting parameters, but it didn't help.
The command for dexseq_count.py:
python /path/DEXSeq/python_scripts/dexseq_count.py \ --format bam \ --order pos \ --paired yes \ --stranded reverse \ $ref_gff \ $bam \ $counts_txt
The command for GenomicAlignments:
txdb = makeTxDbFromGFF(file = genes_gtf, format = "gtf") exons = disjointExons(txdb, aggregateGenes = FALSE) SE = summarizeOverlaps(features = exons, reads = bam_list, singleEnd = FALSE, ignore.strand = FALSE, inter.feature = FALSE, fragments = TRUE)
Is there a mistake with one of the methods?