htseq-count error, AssertionError

[saskak]

Hi all,

I am having a bit of a problem with htseq-count. I have single-end stranded RNA-seq experiment, aligned with tophat-2.0.14. SAM file is names sorted and using HTSeq-0.6.1p1 on Mac 10.7.5 with python 2.7. I also tried a few other versions of HTSeq, but it comes out with the same result.

When trying to run htseq-count:
htseq-count --stranded=reverse -t exon --idattr=gene_id --mode=intersection-nonempty ./test.sam ./gtf.gtf > ./htseq_count.txt

I get this error message:
Error:
[Exception type: AssertionError, raised in __init__.py:599]

Here are a few lines of my sam file:
HWI-ST1146:240:C5HD2ACXX:1:1101:1040:33602 16 X 13759664 50 100M * 0 0 GGATGCCTCCAATGACCGACTGAAGCCGACCAAGTGGCTGACCGCCACCGAGCTGGAGAACGTGCCCTNCCTCAACGACATCACCTGGGAGCGTTTGGAG
C<C@BBCCCCECB>9BCCDCCC<<5BBC>CCCCCBCC?@>9CDA<HGGIFIIGHFF<HHIIIIGGB:0#IGIHGFHFEHEDHEIIIIHHH@HFFFFF@@@ AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:68C31 YT:Z:UU NH:i:1 XS:A:+
HWI-ST1146:240:C5HD2ACXX:1:1101:1041:2544 0 2L 450124 50 68M1I30M * 0 0 GATCCATTCGCCCACCGGCTTGGTTACGTACNACTTTCCGCCTTCATGGTTCAGCATAAAGTTTATGATTTTTTTTGACGATCTTGCGCTGGTCTCGCA
CCCFFFFFHFHHHJJJIIIJJHIBBFEHHHI#08?FHIIJJJGIJEHIJ=@AEA?AHFBE@7;@CCCCEEEDCB@@>?CBB@BDDCDDDBBBCDCDC<@ AS:i:-9 XN:i:0 XM:i:1 XO:i:1 XG:i:1 NM:i:2 MD:Z:31A66 YT:Z:UU NH:i:1 XS:A:-
HWI-ST1146:240:C5HD2ACXX:1:1101:1041:10782 0 2L 2191942 50 100M * 0 0 CGCTCATTCACATAGGTCACCTTGGTGGCCGNTTGACCCTGGCCCTGGCTCTTGTCCTTCTCCTTTTCGGATATCGTCCTCTCGATGGAGCAGGAGGCCG CC@FFFFFGDFHHJJJEGGJI@HIJHJJJJJ#1?DFHIGIJIIIGHIIGIJJJJIJIJJJJJIJEEHDE6BFCCE=CBDDDDDDCDDDCBDDDDD?BBB5 AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:31C68 YT:Z:UU NH:i:1 XS:A:-
HWI-ST1146:240:C5HD2ACXX:1:1101:1041:30995 16 3R 7900087 50 100M * 0 0 AAGGCAAGGATCCACCTTTCCTGCTAGAAGTCAGCGCACAAATCACTTTTGTTCTGTTCCTGCTGTTCNTGACCATCATTCTGATGAACCTGCTCGTGGG @:?CA>ACA><B@:9DDDBADDCA@CDDCEADFFHCHCIHG@@@7FHDHJIJIJIIHIGFIHEGB?11#IHGIIIIHFCHJJHGFEBHHHFFDDFFFCC@ AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:68G31 YT:Z:UU NH:i:1 XS:A:+
HWI-ST1146:240:C5HD2ACXX:1:1101:1041:48103 16 2L 14617161 50 100M * 0 0 AACCAGAACGGAGCCATCTGGAAACTGGACTTGGGCGCCCTGGAGGCCATCCAGTGGACCAAGCACTGNGACTCCGGCATCTAAGAAGTGATACTCCCAA
>BBBABBBBBBA4>BAA>;A@;5>@:BB?B?695(ECECEE@FAC@CCFEFFFDEFIEF??;F?:0#IIFGFCBAIHHEBEEIIIFFFA:=DDDD@@@ AS:i:-3 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:36A31G31 YT:Z:UU NH:i:1 XS:A:+

If I use another experiment with paired-end reads it works without a problem. As far as I can get it is an error related to "paired-end handling", although my reads are single-end.
It would appear that htseq-count recognises my reads as paired-end rather than single-end. It also does not matter if I have selected stranded= yes,no or reverse, same error.

Would really appreciate any help.

[saskak]

Forgot to say that tophat was supplied with --library fr-firststrand

[dpryan]

What happens if you do a "samtools view -f 1 -c test.sam"? If you have even a single alignment with an incorrect pairing flag then that'll break things.

[saskak]

What happens if you do a "samtools view -f 1 -c test.sam"? If you have even a single alignment with an incorrect pairing flag then that'll break things.

I get 0 with samtools view -f 1 -c test.sam

All of the flags are either 0 or 16. htseq-count breaks on the first read it encounters.

[saskak]

It is probably something to do with my htseq installation, as the same input works on fedora 20.
I still cannot understand, if it is installation problem, why would the same installation work with paired-end data and not with single-end.

[dpryan]

Very strange indeed. To tell you the truth, I'm not sure where exactly that assert is being called. The most likely __init__.py on that version of HTSeq that I have installed locally doesn't have an assert() on that line. I suspect that you're correct about this being a weird installation problem. Perhaps just deinstalling/reinstalling will fix things.

[GenoMax]

@saskak: Are those two installs using two different default versions of python?

[saskak]

@saskak: Are those two installs using two different default versions of python?

These are two different installs, one on MacOsX and another on fedora 20. Both are the default versions of python (2.7)