[yh253]

Hi Simon,

I was trying to use htseq-qa to assess the technique quality of my aligned sam file, but I've encountered the following errors. While, when I used the command on my solexa-fastq file, I got the quality plot successfully. My sam file was generated by bwa-0.5.7.

$htseq-qa -t sam q -r 30 s_8.sam
Traceback (most recent call last):
File "/Library/Frameworks/Python.framework/Versions/2.6/bin/htseq-qa", line 5, in <module>
pkg_resources.run_script('HTSeq==0.4.3-p4', 'htseq-qa')
File "/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/setuptools-0.6c11-py2.6.egg/pkg_resources.py", line 489, in run_script
File "/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/setuptools-0.6c11-py2.6.egg/pkg_resources.py", line 1207, in run_script
File "/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/HTSeq-0.4.3_p4-py2.6-macosx-10.3-fat.egg/EGG-INFO/scripts/htseq-qa", line 5, in <module>
HTSeq.scripts.qa.main()
File "/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/HTSeq-0.4.3_p4-py2.6-macosx-10.3-fat.egg/HTSeq/scripts/qa.py", line 124, in main
r.add_qual_to_count_array( qual_arr_A )
File "_HTSeq.pyx", line 715, in _HTSeq.SequenceWithQualities.add_qual_to_count_array (src/_HTSeq.c:12251)
File "_HTSeq.pyx", line 734, in _HTSeq.SequenceWithQualities.add_qual_to_count_array (src/_HTSeq.c:12169)
ValueError: Too large quality value encountered.

$ htseq-qa -t solexa-fastq -r 30 s_8_sequence.txt (This time, it works with fastq file).

I am not sure is this a problem with my BWA alignment or with htseq-qa. It would be very much appreciated if you could put some of your input here!

Yuan

[marcora]

Dear Simon,

thanx for the clear explanation. Is the lax part of GFF that makes going from GTF to GFF "difficult" sometimes, for example when a piece of software requires GFF with specific "comments" (tophat?).

Thanx again for your time and consideration,

Dado

[Simon Anders]

Hi Marcora

Yes, there is a very robust GTF->GFF converter available: Just don't do anything, because every GTF file is a GFF file as well.

GTF is a tightening of the GFF specification. This means: If your file has tab-separated fields with the contents <seqname> <source> <feature> <start> <end> <score> <strand> <frame> [attributes] [comments], it is a GFF file. The GFF specs are a bit lax about how certain columns are to be filled. Should the ID in the attributes field be called "ID" or "gene_ID" or "gene"? Which words should be used in the feature column? If you want to have a general format, it is hard to give clear rules, but once you have agreed that you want to describe not any kind of feature, but specifically gene models, you can be more explicit. This is what the GTF specification does: it explains how precisely a GFF file should look like if it is used to describe gene models, and if a GFF file follows these rules, it is called a GTF file.

Specifically for htseq-count: If you want to count reads in genes and have a GTF file, you can use it out of the box. If you want to count reads in some other kind of feature, and your GFF file hence cannot follow the GTF specs, you have to tell htseq-count which feature types it should use and how the field with the ID is named. (By default, it takes the lines with feature type "exon" and looks for the ID in the attribute field "gene_id", which is what makes sense for GTF files.)

I hope that clarifies it.

Simon

[marcora]

Thanx a lot Simon.

One more thing. It is unclear from your comments here and from the doc online whether HTseq handles both GTF and GFF interchangeably. I am new to this bioinformatics business, and already all these formats are giving me an headache, expecially when GTF files are easily available but no standard/robust GTF>GFF converter is readily available.

Cheers

[joro]

Great - thanks for the quick reply.

[Simon Anders]

Hi

This is the error I am getting:

Error: invalid literal for int() with base 10: '0.000000'
[Exception type: ValueError, raised in __init__.py:200]

I noticed this bug myself just yesterday and fixed it. Please try again with version 0.4.3-p4 and tell me whether this solves the issue.

Cheers
Simon

[marcora]

Dear Simon,

thanx for this package.

So far everything works except when I try to use htseq-count using tophat output sam file as input and a refseq gff file that has worked just fine with tophat.

This is the error I am getting:

Error: invalid literal for int() with base 10: '0.000000'
[Exception type: ValueError, raised in __init__.py:200]

[Simon Anders]

\I am using htseq-count with the -q option. However I'm still getting warnings telling me that htseq-count encountered a read, which has been aligned to a chromosome that did not appear in the GFF file.

I've just fixed this. It now runs properly quiet with the '-q' option.

Simon

[joro]

Hi Simon,

I am using htseq-count with the -q option. However I'm still getting warnings telling me that htseq-count encountered a read, which has been aligned to a chromosome that did not appear in the GFF file. Any ideas on how to resolve this?

The command I'm using is:
htseq-count -q <sam_file> <gff_file>

[Thomas Doktor]

Hi Simon

Thanks for the updated source. I did sort my sam file prior to analysis and although most read pairs seem to be in adjacent lines, some reads are lacking a mate. I suspect this is because TopHat does not discard unmated reads. The question is if I should remove these unmated reads or if the script considers them in the read count and merely displays a warning?

As it turns out, my GFF file is lacking the mitochondrial encoded genes.

On another note, the script seems to read the GFF file before checking if the sam file exists and as my GFF file is quite large it takes a couple of minutes for the script to exit when I - as happens - sometimes forget to supply an existing sam file. I think it would be nice to have the script check that the files exist as the first thing and then exit immediately upon error.

[Simon Anders]

Hi Thomas

I'm trying to use htseq-count version 0.4.2-p3 on a sam file produced by TopHat and a hg19 Ensembl GTF file. I'm analysing the reads in non-stranded mode and looking for exons in the gene_id features. The script runs for a while and outputs several warnings about reads incorrectly flagged as proper pairs, but then exits with the following error:
Is this an error in my sam file and if so how can I identify the read in question?

Line 100 in count.py, where the error occured, was supposed to tell you that htseq-count encountered a read, which has been aligned to a chromosome that did not appear in the GFF file, and that the read has hence been skipped. I've corrected the bug, in version 0.4.2-p4, you not get a proper warning, telling you the ID of the offending read and the name of the unknown chromosome. Please try again.

As for the warnings about improper pairs: Have you sorted your SAM file before calling htseq-count? This is necessary to make sure that the read pairs appear in adjacent lines (see man page).

Simon

[spenthil]

I'm trying to use htseq-count version 0.4.2-p3 on a sam file produced by TopHat and a hg19 Ensembl GTF file. I'm analysing the reads in non-stranded mode and looking for exons in the gene_id features. The script runs for a while and outputs several warnings about reads incorrectly flagged as proper pairs, but then exits with the following error:
Is this an error in my sam file and if so how can I identify the read in question?

Could you post the relevant code?

[Thomas Doktor]

I'm trying to use htseq-count version 0.4.2-p3 on a sam file produced by TopHat and a hg19 Ensembl GTF file. I'm analysing the reads in non-stranded mode and looking for exons in the gene_id features. The script runs for a while and outputs several warnings about reads incorrectly flagged as proper pairs, but then exits with the following error:

Error: 'tuple' object has no attribute 'read'
[Exception type: AttributeError, raised in count.py:100]

Is this an error in my sam file and if so how can I identify the read in question?

[fennan]

Nice work. I will start using it. I'll report my experience!

[Simon Anders]

Hi tcezard

thanks for the bug report and the nice code example. I've found and fixed the bug and uploaded a new release path, version now 0.4.2-p3.

If you find further bugs, just send me an e-mail.

Cheers
Simon (anders at embl dot de )