Counting exonic reads with DEXSeq and GenomicAlignments
Entering edit mode
igor ▴ 40
Last seen 8 months ago
United States

The DEXSeq vignette describes several protocols for obtaining exon counts. One option is to use the 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

python /path/DEXSeq/python_scripts/ \
--format bam \
--order pos \
--paired yes \
--stranded reverse \
$ref_gff \
$bam \

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?

dexseq genomicalignments exon usage • 949 views

Login before adding your answer.

Traffic: 1042 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6