read counts per transcript, edgeR, tophat

[chrisbala]

Hi everyone,

I'm using Illumina RNAseq to test for differential expression between two conditions. I have successfully mapped reads with tophat/bowtie, and used cufflinks/cuffcompare to test for differential expression.

A number of methods (edgeR, DEGseq, DESseq), however, ask for raw read counts rather FPKM.

Is there an easy way to go from tophat output (say the bedgraph or .sam) to raw read counts for each transcript (say using a gtf that I would supply)?

Thanks!

[quinlana]

If I understand your question correctly, BEDTools has a utility called coverageBed that will compute the raw counts of reads in a BAM file among features in a BED or GTF file.

For example:

$ coverageBed -abam alignedReads.bam -b transcripts.gtf > transcript.cov.txt

There is also a histogram (-hist) mode that will give you a histogram of coverage for each feature in the -b file. I have not yet updated the manual , but the help (-h) for this tools describes the output format.

[Simon Anders]

Here is the script that I made for this purpose: http://www-huber.embl.de/users/ander...doc/count.html

Simon

[bioinfosm]

Yeah, I have used bedtools for similar analysis and it is quite intuitive

Originally Posted by quinlana

If I understand your question correctly, BEDTools has a utility called coverageBed that will compute the raw counts of reads in a BAM file among features in a BED or GTF file.

For example:

$ coverageBed -abam alignedReads.bam -b transcripts.gtf > transcript.cov.txt

There is also a histogram (-hist) mode that will give you a histogram of coverage for each feature in the -b file. I have not yet updated the manual , but the help (-h) for this tools describes the output format.

[Auction]

Originally Posted by Simon Anders

Here is the script that I made for this purpose: http://www-huber.embl.de/users/ander...doc/count.html

Simon

Simon, are there any specific reason to use Python2.5, because the cluster I can use only have Python2.4. And I trid to build the source codes using Python 2.4. When I ran "python setup.py build" I got the warning message which doesn't matter.

/usr/lib64/python2.4/distutils/dist.py:236: UserWarning: Unknown distribution option: 'requires'
warnings.warn(msg)

However, when I ran "python setup.py install --prefix=...", I got the error message

byte-compiling $HOME/librarys/lib64/python2.4/site-packages/HTSeq/StepVector.py to StepVector.pyc
File $HOME/librarys/lib64/python2.4/site-packages/HTSeq/StepVector.py", line 390
start = index.start if index.start is not None else self.start_index()
^
SyntaxError: invalid syntax

The reason seems to that "Conditional expressions" is not supported in Python2.4, how can I change the codes without upgrading to Python 2.5

Thanks

Ying

[Simon Anders]

When I started developing HTSeq, I had to make a decision on which Python version to require. I refrained from using any Python 2.6 features as many people still work with 2.5, but I didn't expect many 2.4 installations to be still around.

Python 2.5 brought a few nifty feature but I think I didn't use too many of them, and most are only syntactic sugar anyway. I guess that not too many changes would be needed to make it run on 2.4

But wouldn't it be much easier if you convinced your cluster admin to update the system? (Or install it yourself in your home directory: This takes only a couple of minute and Python typically builds from source out of the box.)

Python 2.5 was released 4 years ago, and by now, the release version is 2.6.5, and 2.7 is in alpha. You shouldn't need to waste your time to get software to run in such an out-of-date system.

[Siva]

Originally Posted by Simon Anders

Python 2.5 was released 4 years ago, and by now, the release version is 2.6.5, and 2.7 is in alpha. You shouldn't need to waste your time to get software to run in such an out-of-date system.

Simon:
You would be surprised to know that one of the clusters I use needs Python 2.3!! Hope fully they update it soon. Fortunately we have another cluster but some one else has to help me load HTseq on it as I do not have access to it.

Siva

[Auction]

Originally Posted by Simon Anders

When I started developing HTSeq, I had to make a decision on which Python version to require. I refrained from using any Python 2.6 features as many people still work with 2.5, but I didn't expect many 2.4 installations to be still around.

Python 2.5 brought a few nifty feature but I think I didn't use too many of them, and most are only syntactic sugar anyway. I guess that not too many changes would be needed to make it run on 2.4

But wouldn't it be much easier if you convinced your cluster admin to update the system? (Or install it yourself in your home directory: This takes only a couple of minute and Python typically builds from source out of the box.)

Python 2.5 was released 4 years ago, and by now, the release version is 2.6.5, and 2.7 is in alpha. You shouldn't need to waste your time to get software to run in such an out-of-date system.

The reason I stick on Python 2.4 is that several clusters I can use all use Rock 5.2+Python2.4. By the way, I upgraded my workstation to Python2.6 and installed HTSeq, but it doesn't indicate support for paired-ends. Are there specific options for paired-ends? Thanks

[Auction]

Originally Posted by quinlana

If I understand your question correctly, BEDTools has a utility called coverageBed that will compute the raw counts of reads in a BAM file among features in a BED or GTF file.

For example:

$ coverageBed -abam alignedReads.bam -b transcripts.gtf > transcript.cov.txt

There is also a histogram (-hist) mode that will give you a histogram of coverage for each feature in the -b file. I have not yet updated the manual , but the help (-h) for this tools describes the output format.

I tried to run coverageBed for my BAM file from pair-end RNA-seq. When I compared the result in a certain gene to my calculation using Bio:B::Sam. It seems coverageBed don't consider the mate properly paired flag in the BAM file. Are there any options control for it? Thanks

[quinlana]

Originally Posted by Auction

I tried to run coverageBed for my BAM file from pair-end RNA-seq. When I compared the result in a certain gene to my calculation using Bio:B::Sam. It seems coverageBed don't consider the mate properly paired flag in the BAM file. Are there any options control for it? Thanks

Hi,
You are correct, coverageBed processes every BAM entry in the BAM file provided. The rationale here is that one can use samtools to filter your data before passing to coverageBed. For example, to only count proper pairs, do the following (passing the output of samtools to coverageBed via a pipe, hence the "stdin" for the -abam parameter:

$ samtools -f 0x2 <BAM> | coverageBed -abam stdin -b transcripts.bed > transcripts.proper.cov

Does this help?

[Auction]

Originally Posted by quinlana

Hi,
You are correct, coverageBed processes every BAM entry in the BAM file provided. The rationale here is that one can use samtools to filter your data before passing to coverageBed. For example, to only count proper pairs, do the following (passing the output of samtools to coverageBed via a pipe, hence the "stdin" for the -abam parameter:

$ samtools -f 0x2 <BAM> | coverageBed -abam stdin -b transcripts.bed > transcripts.proper.cov

Does this help?

I see, thanks.

[chrisbala]

Hi Simon,

Sorry to trouble you with with what is surely a silly question. I've just installed HT-Seq from source on my Mac. I've just tried to run a couple of basic commands to see if things were working:

python -m HTSeq.scripts.qa s_3_sequence.txt

but I'm getting an error message:

Traceback (most recent call last):
  File "/System/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/runpy.py", line 85, in run_module
    loader = get_loader(mod_name)
  File "/System/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/pkgutil.py", line 456, in get_loader
    return find_loader(fullname)
  File "/System/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/pkgutil.py", line 466, in find_loader
    for importer in iter_importers(fullname):
  File "/System/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/pkgutil.py", line 422, in iter_importers
    __import__(pkg)
ImportError: No module named HTSeq.scripts

I'm running the script from within the HT-Seq directory. I placed the HT-Seq folder within my Applications directory.

Any idea what I've done wrong?

[Simon Anders]

Originally Posted by chrisbala

I've just installed HT-Seq from source on my Mac. I've just tried to run a couple of basic commands to see if things were working:

python -m HTSeq.scripts.qa s_3_sequence.txt

but I'm getting an error message:

It seems that Python does not know were you put the HTSeq modules.

Once you have unpacked the HTSeq tarball, you should run the commands

   python setup.py build
   python setup.py install

The second command puts it somewhere where Python should find it.

If this does not work, use only the first command. It makes a subdirectory 'build', and in there something like 'lib.macos' (or called similarly), and in there, 'HTSeq'. You can just tell Python to search for HTSeq there by writing

export PYTHONPATH=$PYTHONPATH:/<the_path_to_your_directory/build/lib.macos/

and Python will find it there.

Do not run the script in the directory that you have just unpacked, go somewhere else, as Python otherwise tends to get confused.

[chrisbala]

Thanks Simon,

THIngs seem to be working...

Will be try DESeq as soon as I get these counts finalized.

One question: for the counting script, you ask for a gff, and elsewhere you use gtf. I've noticed in tophat as well that is some cases gffs are requested (or output) and sometimes gtf.

Is there are particular aspect of gff that makes it more suitable for this counting algorithm?

just curious...

[Simon Anders]

Hi Chris,

GTF is just a tightening of the GFF specification, i.e., every GTF file is also a GFF file. The tightening is mainly with respect to the "feature type" column: in GTF, it has to be one of "exon", "CDS", "stop_codon", etc., while in GFF, it is not clearly said which words should be used.

In the htseq-count script, I give the user the option to tell the script how the features are named, and hence, it works with GFF as well.

See here:
GFF: http://www.sanger.ac.uk/resources/so.../gff/spec.html
GTF: http://mblab.wustl.edu/GTF22.html

Simon

[wenhuang]

Hi Simon,

Can I ask you a question that if I have paired end reads who are overlapping with each other, for a overlapped region, does HTSeq count it as single coverages if reads are from the same fragment (but different ends) or count it as covered twice?

Thanks!

Originally Posted by Simon Anders

Here is the script that I made for this purpose: http://www-huber.embl.de/users/ander...doc/count.html

Simon

[Simon Anders]

Originally Posted by wenhuang

Can I ask you a question that if I have paired end reads who are overlapping with each other, for a overlapped region, does HTSeq count it as single coverages if reads are from the same fragment (but different ends) or count it as covered twice?

If you use the full Python framework (as opposed to the the simple htseq-count script) you can adjust such stuff as you need them.

At the moment, however, HTSeq does not support paired-end reads at all. I'll add this functionality very soon, hopefully next week.

Your data would be very useful as test data. Could you maybe send me your SAM file (or parts of it)?

By the way, last time I tried, Bowtie was unable to deal with overlapping ends from the same fragment and did not align these at all. Have you used another software, or did they improve that?

Cheers
Simon

[wenhuang]

Hi Simon,

Thanks for your reply. I ended up using tophat and it has no problem aligning overlapping reads by specifying a negative -r as suggested by Cole Trapnell.

I attached an example SAM file below. It is from a highly expressed gene on bovine chromosome one (ATPO). The Ensembl annotation is as below.

I am very eager to see what you get.

Thanks a lot.


chr1 protein_coding exon 720699 720764 . + . gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "1"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding CDS 720729 720764 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "1"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN"; protein_id "ENSBTAP00000024326";
chr1 protein_coding start_codon 720729 720731 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "1"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding exon 722101 722151 . + . gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "2"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding CDS 722101 722151 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "2"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN"; protein_id "ENSBTAP00000024326";
chr1 protein_coding exon 723098 723208 . + . gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "3"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding CDS 723098 723208 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "3"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN"; protein_id "ENSBTAP00000024326";
chr1 protein_coding exon 725362 725491 . + . gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "4"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding CDS 725362 725491 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "4"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN"; protein_id "ENSBTAP00000024326";
chr1 protein_coding exon 726174 726286 . + . gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "5"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding CDS 726174 726286 . + 2 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "5"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN"; protein_id "ENSBTAP00000024326";
chr1 protein_coding exon 727512 727598 . + . gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "6"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding CDS 727512 727598 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "6"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN"; protein_id "ENSBTAP00000024326";
chr1 protein_coding exon 727876 728057 . + . gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "7"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";
chr1 protein_coding CDS 727876 727986 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "7"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN"; protein_id "ENSBTAP00000024326";
chr1 protein_coding stop_codon 727987 727989 . + 0 gene_id "ENSBTAG00000018278"; transcript_id "ENSBTAT00000024326"; exon_number "7"; gene_name "ATPO_BOVIN"; transcript_name "ATPO_BOVIN";

Originally Posted by Simon Anders

If you use the full Python framework (as opposed to the the simple htseq-count script) you can adjust such stuff as you need them.

At the moment, however, HTSeq does not support paired-end reads at all. I'll add this functionality very soon, hopefully next week.

Your data would be very useful as test data. Could you maybe send me your SAM file (or parts of it)?

By the way, last time I tried, Bowtie was unable to deal with overlapping ends from the same fragment and did not align these at all. Have you used another software, or did they improve that?

Cheers
Simon

[Simon Anders]

Hi wenhuang

I've now added paired-end support for the SAM fomat to HTSeq and htseq-count. It would be great if you could download the new version 0.4.1 of HTSeq, try to do your counting job with it and tell me whether it worked correctly.

Cheers
Simon

[wenhuang]

Hi Simon,

Thank you for your work. Unfortunately, although I am able to get the yeast data provided on your website to work with htseq-count, I am not able to get my own data to work. The SAM file I used is two posts ago and I attach the GTF file to this post.

The error message I got was as below.

Again, thank you for your help.

Traceback (most recent call last):
File "/usr/local/bin/htseq-count", line 5, in <module>
pkg_resources.run_script('HTSeq==0.4.1', 'htseq-count')
File "/System/Library/Frameworks/Python.framework/Versions/2.5/Extras/lib/python/pkg_resources.py", line 489, in run_script

File "/System/Library/Frameworks/Python.framework/Versions/2.5/Extras/lib/python/pkg_resources.py", line 1207, in run_script
path = self.module_path
File "/Library/Python/2.5/site-packages/HTSeq-0.4.1-py2.5-macosx-10.5-i386.egg/EGG-INFO/scripts/htseq-count", line 5, in <module>
HTSeq.scripts.count.main()
File "/Library/Python/2.5/site-packages/HTSeq-0.4.1-py2.5-macosx-10.5-i386.egg/HTSeq/scripts/count.py", line 156, in main
opts.mode, opts.featuretype, opts.idattr, opts.quiet )
File "/Library/Python/2.5/site-packages/HTSeq-0.4.1-py2.5-macosx-10.5-i386.egg/HTSeq/scripts/count.py", line 68, in count_reads_in_features
if iv.chrom not in features.step_vectors:
AttributeError: 'NoneType' object has no attribute 'chrom'

Originally Posted by Simon Anders

Hi wenhuang

I've now added paired-end support for the SAM fomat to HTSeq and htseq-count. It would be great if you could download the new version 0.4.1 of HTSeq, try to do your counting job with it and tell me whether it worked correctly.

Cheers
Simon