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.

Best,

Dimitra

Thanks Dimitra for your report and for the fix! I have implemented it to the latest versions of DEXSeq.

Alejandro

Hi,

I am using the 3.2 version but I got the same problem:

python /home/julian/R/x86_64-pc-linux-gnu-library/3.2/DEXSeq/python_scripts/dexseq_count.py -f bam dexseq/annotation/RefGene.fixed.gff 2monthsCTX/WT/WT1/mm10.bam dexseq/annotation/dexseq/WT1fb.txt

Traceback (most recent call last):
  File "/home/julian/R/x86_64-pc-linux-gnu-library/3.2/DEXSeq/python_scripts/dexseq_count.py", line 187, in <module>
    if a.optional_field("NH") > 1:
  File "_HTSeq.pyx", line 1399, in HTSeq._HTSeq.SAM_Alignment.optional_field (src/_HTSeq.c:26483)
KeyError: 'SAM optional field tag NH not found

The alignment was done with the STAR aligner.

Any Ideas how it can be approached?

Thanks, Julian