dexseq_count.py vs htseq-count

Question

[anduril50]

Dear All,

I have some RNA paired-end files which I have used tophat to align

~/tophat-2.0.7.Linux_x86_64/./tophat -p10 --segment-length 18 \
    --no-coverage-search --segment-mismatches 1 -G

then samtools

samtools sort -n /media//PRE_6_8_2/accepted_hits.bam \
    /media/PRE_6_8_2/PRE_6_8_2.sort.hg19.bam

and then

samtools view /media/PRE_6_8_2/PRE_6_8_2.sort.hg19.bam.bam \
    /media/PRE_6_8_2/PRE_6_8_2.sort.hg19.sam

then I am trying to use DEXSeq and using the

python2.7 /home/3.0/DEXSeq/python_scripts/dexseq_count.py -p yes -r name -s no \
    /media/PRE_6_8_2/hg19/hg19_IlluminaAnnotation_genes.gtf \
    /media/PRE_6_8_2/PRE_6_8_2.sort.hg19.sam \
    /media/PRE_6_8_2//PRE_6_8_2/untreated1fb.txt

and gives me

_ambiguous 0
_empty 39845364
_lowaqual 0
_notaligned 0

without any reads for any exons. 0 for everything

whereas by using the htseq-count -r name -s no everything looks fine (reads in genes)...

could you please advice me how to solve it? I really need to run the DEXseq as I want reads per exon to detect splice differences.

Thank you in advance for your help

Answer 1

[komal.rathi]

I am assuming you did not run the python script dexseq_prepare_annotation.py on your actual GTF file. Run the command below and use the GFF file instead of GTF file as input in dexseq_count.py:

# prepare annotation
python2.7 /home/3.0/DEXSeq/python_scripts/dexseq_prepare_annotation.py /media/PRE_6_8_2/hg19/hg19_IlluminaAnnotation_genes.gtf /media/PRE_6_8_2/hg19/hg19_IlluminaAnnotation_genes.DEXSeq.gff
# run dexseq_count.py
python2.7 /home/3.0/DEXSeq/python_scripts/dexseq_count.py -p yes -r name -s no /media/PRE_6_8_2/hg19/hg19_IlluminaAnnotation_genes.DEXSeq.gff /media/PRE_6_8_2/PRE_6_8_2.sort.hg19.sam /media/PRE_6_8_2/PRE_6_8_2/untreated1fb.txt

[anduril50]

That's the problem...even though I did it, it still doesn't work.....

Could you please suggest me something else?

[komal.rathi]

Must you use Illumina Annotation? Can you try using any other annotation like Ensembl (Recommended) or Gencode? Section 2.4 of DEXSeq manual says that they have tested the script on Ensembl annotations and it may fail when used on other annotations.