SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
DEXSeq: dexseq_count.py produces lots of warnings (mate could not be found) enomis Bioinformatics 15 11-08-2013 01:48 PM
samtools XS:A optional field charlie_sequencing Bioinformatics 6 05-04-2013 11:19 PM
SAM Format - SEQ field '=' Bio.X2Y Bioinformatics 3 04-25-2012 04:26 AM
Bowtie optional fields in SAM question history_of_robots Bioinformatics 1 11-28-2011 07:56 AM
Extracting one field from a SAM file jdrum00 Bioinformatics 8 01-04-2010 08:40 PM

Reply
 
Thread Tools
Old 10-31-2013, 11:35 AM   #1
Kwangwoo
Junior Member
 
Location: Boston

Join Date: Oct 2013
Posts: 1
Smile DEXSeq error: 'SAM optional field tag NH not found'

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):

Quote:
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:


Quote:
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

Last edited by Kwangwoo; 10-31-2013 at 11:44 AM.
Kwangwoo is offline   Reply With Quote
Old 11-05-2013, 12:47 AM   #2
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

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
areyes is offline   Reply With Quote
Old 02-19-2014, 07:23 AM   #3
Jose Garcia
Junior Member
 
Location: Milan

Join Date: Jul 2013
Posts: 4
Default 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



Quote:
Originally Posted by areyes View Post
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 is offline   Reply With Quote
Old 08-25-2014, 01:08 PM   #4
dejavu2010
Member
 
Location: usa

Join Date: Jan 2012
Posts: 21
Default

where can i find the old python script?
dejavu2010 is offline   Reply With Quote
Old 08-25-2014, 01:32 PM   #5
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 700
Default

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;
}
Richard Finney is offline   Reply With Quote
Old 08-29-2014, 01:36 AM   #6
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

Quote:
where can i find the old python script?
From old versions of DEXSeq, for instance:

http://www.bioconductor.org/packages...ml/DEXSeq.html

But I would rather add the NH flag, rather than using old, not-mantained software!
areyes is offline   Reply With Quote
Reply

Tags
dexseq_count.py, nh tag

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:32 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO