SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Evaluation of genomic high-throughput sequencing data generated on Illumina HiSeq and dan Literature Watch 1 11-09-2011 04:18 AM
ChIP-Seq: Systematic bias in high-throughput sequencing data and its correction by BE Newsbot! Literature Watch 0 06-08-2011 02:50 AM
ChIP-Seq: Data structures and compression algorithms for high-throughput sequencing t Newsbot! Literature Watch 0 10-16-2010 02:00 AM
ChIP-Seq: Savant: Genome Browser for High Throughput Sequencing Data. Newsbot! Literature Watch 0 06-22-2010 02:00 AM
Dave - interested in trial data sets and high-throughput sequencing tutorials Dave S. Introductions 4 02-02-2010 10:04 AM

Reply
 
Thread Tools
Old 05-16-2010, 03:24 PM   #21
marcora
Member
 
Location: Pasadena, CA USA

Join Date: Jan 2010
Posts: 52
Default

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:


Code:
Error: invalid literal for int() with base 10: '0.000000'
[Exception type: ValueError, raised in __init__.py:200]
marcora is offline   Reply With Quote
Old 05-16-2010, 11:14 PM   #22
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi

Quote:
Originally Posted by marcora View Post
This is the error I am getting:


Code:
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
Simon Anders is offline   Reply With Quote
Old 05-17-2010, 12:06 AM   #23
joro
Member
 
Location: UK

Join Date: Feb 2010
Posts: 28
Default

Great - thanks for the quick reply.
joro is offline   Reply With Quote
Old 05-17-2010, 03:18 AM   #24
marcora
Member
 
Location: Pasadena, CA USA

Join Date: Jan 2010
Posts: 52
Default

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
marcora is offline   Reply With Quote
Old 05-17-2010, 04:25 AM   #25
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

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
Simon Anders is offline   Reply With Quote
Old 05-17-2010, 05:24 AM   #26
marcora
Member
 
Location: Pasadena, CA USA

Join Date: Jan 2010
Posts: 52
Default

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
marcora is offline   Reply With Quote
Old 05-19-2010, 07:18 AM   #27
yh253
Member
 
Location: Ireland

Join Date: Jul 2009
Posts: 16
Default

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
yh253 is offline   Reply With Quote
Old 05-19-2010, 07:31 AM   #28
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi Yuan

Quote:
Originally Posted by yh253 View Post
ValueError: Too large quality value encountered.
Usually, SAM files don't contain quality values exceeding 40. When you aligned your _sequence.txt file, did you convert the quality strings from Solexa to Sanger scale? If not, BWA probably might not have found the optimal alignments, and htseq-qa gets confused, too.

If you did, and the large quality values are legitimate, I'd be interested to see your SAM file.

Simon
Simon Anders is offline   Reply With Quote
Old 05-19-2010, 07:41 AM   #29
yh253
Member
 
Location: Ireland

Join Date: Jul 2009
Posts: 16
Default

Hi Simon,

I am afraid not. Actually I was not aware of this issue until now.I will try to figure out how to convert the quality strings first then, and give a feedback on this later.Thanks again for your prompt reply!

Yuan
yh253 is offline   Reply With Quote
Old 05-25-2010, 07:16 AM   #30
gen2prot
Member
 
Location: Hyderabad, India

Join Date: Apr 2010
Posts: 64
Default

Hi Simon,

I am having difficulty in running the htseq-qa script. I think I have installed HTSeq correctly since I get no error message for "import HTSeq" command. Then on giving the "htseq-qa -t sam accepted.sam" command, I get a Syntax error. I have given the following export command in Unix

export PYTHONPATH=$PYTHONPATH:/Library/Python/2.6/

Is this wrong? On giving the command "whereis python", I get /usr/local/python. I am confused.

Thank you
Abhijit
gen2prot is offline   Reply With Quote
Old 05-27-2010, 12:56 PM   #31
gen2prot
Member
 
Location: Hyderabad, India

Join Date: Apr 2010
Posts: 64
Default

Hi Simon,

I have a mac 10.6, intel and python 2.6. Matplotlib is unavailable for 10.6, py 2.6.... as far as I know. I give up. The website (matplotlib) asks me to tinker around with System files which I don't want to tamper with. Thanks anyway for your help.

Abhijit
gen2prot is offline   Reply With Quote
Old 05-27-2010, 01:54 PM   #32
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi Abhijit

Quote:
Originally Posted by gen2prot View Post
I have a mac 10.6, intel and python 2.6. Matplotlib is unavailable for 10.6, py 2.6.... as far as I know.
Yours seems to be a pretty standard configuration, so it would be apity if matplotlib was not available. If you have Xcode installed, building from source might work. Just try:

Code:
wget http://sourceforge.net/projects/matplotlib/files/matplotlib/matplotlib-0.99.1/matplotlib-0.99.1.2.tar.gz/download
tar -xzvf matplotlib-0.99.1.2.tar.gz 
cd matplotlib-0.99.1.1/
python setup.py build
sudo python setup.py install
Simon
Simon Anders is offline   Reply With Quote
Old 05-27-2010, 02:25 PM   #33
gen2prot
Member
 
Location: Hyderabad, India

Join Date: Apr 2010
Posts: 64
Default

Got it to work... Coffee helps.
gen2prot is offline   Reply With Quote
Old 05-31-2010, 05:11 AM   #34
yh253
Member
 
Location: Ireland

Join Date: Jul 2009
Posts: 16
Default

Hi Simon,

I used BWA to do the alignment. If a read mapped to the chromosomal junction (end of one chromosome and beginning of another chromosome), BWA will produce a Flag = 4, but you can still see the other tags: "chr", "CIGAR","MAPQ", etc. in the sam file. This causes a problem for htseq-qa in that it can't process an alignment with flag=4, but mapping position does "NOT" equal "*". One option for me is to pre-process my sam file by excluding all such alignments before giving to htseq-qa, but I am wondering is it possible to turn off this requirement in htseq-qa so that I don't have to change my sam file in advance?Thank you very much for your advice in advance!

Yuan
yh253 is offline   Reply With Quote
Old 06-02-2010, 06:12 PM   #35
gen2prot
Member
 
Location: Hyderabad, India

Join Date: Apr 2010
Posts: 64
Default

Hi Simon,

I made a Drosophila Gene Index file and ran all the solexa reads against this. I have a 5.4 GB sam file. Can I use the Htseq-qa program on this sam file? or will this crash since the subject to which the alignment was done are genes and not chromosomes. Secondly how can I get the read count in my case. The GTF file that I have has the start and end coordinates wrt to the genome. Therefore this won't work with the sam output. Any suggestions?

thanks
Abhijit
gen2prot is offline   Reply With Quote
Old 06-05-2010, 05:49 AM   #36
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi Yuan

Quote:
Originally Posted by yh253 View Post
I used BWA to do the alignment. If a read mapped to the chromosomal junction (end of one chromosome and beginning of another chromosome), BWA will produce a Flag = 4, but you can still see the other tags: "chr", "CIGAR","MAPQ", etc. in the sam file. This causes a problem for htseq-qa in that it can't process an alignment with flag=4, but mapping position does "NOT" equal "*". One option for me is to pre-process my sam file by excluding all such alignments before giving to htseq-qa, but I am wondering is it possible to turn off this requirement in htseq-qa so that I don't have to change my sam file in advance?Thank you very much for your advice in advance!
I've changed the SAM parser such that it now only writes a warning rather than stopping with an error if this case is encountered. Please try again with version 0.4.4p2.

However, I don't quite understand how such SAM lines come about. What is a chromosomal junction? How could a read get mapped partly to one and partly to another chromosome?

Cheers
Simon
Simon Anders is offline   Reply With Quote
Old 06-05-2010, 05:53 AM   #37
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi Abhijit

Quote:
Originally Posted by gen2prot View Post
I made a Drosophila Gene Index file and ran all the solexa reads against this. I have a 5.4 GB sam file. Can I use the Htseq-qa program on this sam file? or will this crash since the subject to which the alignment was done are genes and not chromosomes.
Just try it out. But it should work. htseq-count does not care what you have aligned against.

Quote:
Secondly how can I get the read count in my case. The GTF file that I have has the start and end coordinates wrt to the genome. Therefore this won't work with the sam output. Any suggestions?
So, basically, wherever you usually have a chromosome name (i.e., in the third column of your SAM file) you now have a gene name, and you want to count how many reads fall onto each gene? Why don't you just cut out this third column and count how often each gene appears there?

Of course, your approach has other dangers. For example, how does your aligner handle multiple transcripts? If the same exon appears in several transcripts, the aligner might think it is a repeat.

Cheers
Simon
Simon Anders is offline   Reply With Quote
Old 06-05-2010, 06:11 AM   #38
yh253
Member
 
Location: Ireland

Join Date: Jul 2009
Posts: 16
Default

Hi Simon,

Thanks very much for your reply! Quoted from BWA:"Internally BWA concatenates all reference sequences into one long sequence. A read may be mapped to the junction of two adjacent reference sequences. In this case, BWA will flag the read as unmapped, but you will see position, CIGAR and all the tags". This has happened to me that some reads mapped to the end of chrY and the beginning of chrM by using BWA. However, I didn't see these mappings by using Bowtie. This might arise from mapping strategy adopted by BWA, which converts "N" in the reference sequence into random bases. During mapping, BWA must have converted "N"s at the end of ChrY into random bases which happened to match the beginning bases of some reads.

Cheers,
Yuan

Last edited by yh253; 06-05-2010 at 06:17 AM.
yh253 is offline   Reply With Quote
Old 07-13-2010, 12:21 PM   #39
rkusko
Junior Member
 
Location: USA

Join Date: Jul 2010
Posts: 4
Default

Hi Simon,
I've used HTSeq-counts to extract raw counts from Cufflinks output, and it worked correctly for 40/42 of my samples. For some reason, two samples raised the error:

Error: 'generator' object has no attribute 'get_line_number_string'
[Exception type: AttributeError, raised in count.py:126]


I haven't been able to find anything wrong with the data in those two samples. Any ideas? Thank you!
rkusko is offline   Reply With Quote
Old 07-14-2010, 08:43 AM   #40
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi

Quote:
Originally Posted by rkusko View Post
I've used HTSeq-counts to extract raw counts from Cufflinks output, and it worked correctly for 40/42 of my samples. For some reason, two samples raised the error
Could you maybe send me (be e-mail to anders(at)embl(dot)de) some excerpts from the data that produce the error? Then I could investigate.

Cheers
Simon
Simon Anders is offline   Reply With Quote
Reply

Tags
htseq, python

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 09:11 AM.


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