Hi all (and Alejandro!)
I have come across a problem with the DEXSeq count script (DEXSeq_1.14.2), namely getting the following error while wanting to produce the counts per exon bin by running dex_seq_count.py:
KeyError: 'SAM optional field tag NH not found'
I'm using gsnap for alignment and when working on paired end data the count script went through. It was only when working on single read that I got this type of error, which was corrected once I removed the unmapped reads from my bam file.
Looking at the script itself dexseq_count.py) lines 182-189, I think I know the reason:
if not is_PE: for a in reader( sam_file ): if a.optional_field("NH") > 1: continue if not a.aligned: counts[ '_notaligned' ] += 1 continue
So initially, the first if clause looks for the NH flag and then there is the control whether or not the read is aligned. Reversing the order of these two, so first checking whether the read has been aligned has now worked for me. I just wanted to let people with the same problem know how they could work around this, but I guess this should be as well changed in the next version.
Hi, the problem was that I used the wrong bam files. works now. Sorry for the useless post.