DEXSeq error: 'SAM optional field tag NH not found'

[Kwangwoo]

Hi all (and Simon & Alejandro, hopefully),

In my Bam files, optional NH tag was given to only mapped reads (NH:i:n; n>=1). There's no NH:i:0. Here is a part of the sam (PE; sorted by read name):

H13P7ADXX130822:1:1101:1122:27935 77 * 0 0 * * 0 0 CGGGNGGATGCCTTCTTTATCTTGGATCTTGGCNTTCACATTTTCGATGGTGTCACTGGGCTCCACCTCCAGAGTG CCCF#2ADHHHH
HJJJJJJJJJJJJIJJJJJJJ#1CFHIIJJJJJJJJJJGHHIJJHGEGGGIGHIIIEHHGDB7? PG:Z:MarkDuplicates.1E RG:Z:H13P7.1
H13P7ADXX130822:1:1101:1122:27935 141 * 0 0 * * 0 0 NNNCNCCATCGNAAATGTGAAGGCCAAGATCCAGGATAAAGAAGGCATCCCTCCCGACCAGCAGAGGCTCATCTTT ###2#2:?@@@#
4@@@@@@@@@@@@@@@@????8?????????????????1>?###################### PG:Z:MarkDuplicates.1E RG:Z:H13P7.1
H13P7ADXX130822:1:1101:1124:37458 1609 3 40500185 255 43M2670N33M = 40500185 0 AGTGNATGCAGCTCACTGATTTCATCCTCAAGTTTCCGCACAGTGCCCACCAGAAGTATGTCCGACAA
GCCTGGCA <<<@#2@((=?>7=<<6::?=9@=>>@<9><98)9@@6???8?33;>?>9;7::>6.6>?6)8=?8?8??;=0=;7 PG:Z:MarkDuplicates.1E RG:Z:H13P7.1 NH:i:1 NM:i:1 UQ:i:2 XS:A:+

I think DEXSeq-count.py doesn't like this format showing an error message below:

Traceback (most recent call last):
File "/home/kkim/softwares/R_packages/DEXSeq/python_scripts/dexseq_count.py", line 225, in <module>
rs = map_read_pair( af, ar )
File "/home/kkim/softwares/R_packages/DEXSeq/python_scripts/dexseq_count.py", line 134, in map_read_pair
if af != None and af.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'

Actually, this bam file has worked well with HTSeq-count without an error. The manual of HTSeq says ......If the aligner does not set this field, multiply aligned reads will be counted multiple times. ... I think HTSeq can handle my format.

Do you have any suggestion for avoiding this error message by dexseq-count? What if I grep the sam file by "NH" before running dexseq-count?? There was no error but I cannot figure out how does it affect to my result.

Thanks in advance for your reply.


Kwangwoo

[areyes]

Hi Kwangwoo,

Thanks for your post. This is a recent change introduced to the python script to make sure to only consider uniquely mapped reads.

Previously, the scripts were checking for a good alignment quality (e.g. more than 10) in the sam file, however we noticed that some aligners were giving very good alignment qualities to multiple mapping reads. If a read mapping multiple times should have a good quality score is debatable.

Therefore we decided to introduce this change and instead check for the NH flag, and a read will only be considered if its value is 1... but I did not know it was an optional flag? What aligner are you using?

A quick fix for you would be just to prefilter your reads for those with a unique alignment, e.g. NH:i:1. But I would first make sure that all the uniquely aligned reads have the NH flag!

Alejandro

[Jose Garcia]

NH not found

Same error for me. We are using Soapsplice and as Alejandro was saying, nothing happened with the old dexseq_count.py.
Jose

Hi Kwangwoo,

Thanks for your post. This is a recent change introduced to the python script to make sure to only consider uniquely mapped reads.

Previously, the scripts were checking for a good alignment quality (e.g. more than 10) in the sam file, however we noticed that some aligners were giving very good alignment qualities to multiple mapping reads. If a read mapping multiple times should have a good quality score is debatable.

Therefore we decided to introduce this change and instead check for the NH flag, and a read will only be considered if its value is 1... but I did not know it was an optional flag? What aligner are you using?

A quick fix for you would be just to prefilter your reads for those with a unique alignment, e.g. NH:i:1. But I would first make sure that all the uniquely aligned reads have the NH flag!

Alejandro

[dejavu2010]

where can i find the old python script?

[Richard Finney]

FYI: I run a custom program (getnh) to slap the NH tag on ...
samtools view original.bam | ./getnh | python dexseq_count.py -f sam --paired=yes hg19.gff - outputcounts.txt 2> errfile

This takes the bam file with missing nh tags and outputs as a sam file with new NH tags which gets piped into python dexseq_count.py

/*
how to compile : gcc -Wall -O2 -o getnh getnh.c
*/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main()
{
char *z;
char s[5000];

s[0] = (char)0;
while (gets(s))
{
z = strstr(s,"NH:i");
if (z)
{
printf("%s\n",s);
s[0] = (char)0;
continue;
}
printf("%s\tNH:i:1\n",s);
s[0] = (char)0;
}
return 0;
}

[areyes]

where can i find the old python script?

From old versions of DEXSeq, for instance:

http://www.bioconductor.org/packages/2.12/bioc/html/DEXSeq.html

But I would rather add the NH flag, rather than using old, not-mantained software!

[jrounds]

I aligned with bowtie2, and according to what I read it does not set this flag.

https://www.biostars.org/p/57313/