SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
htseq-count only not-unique counts moser Bioinformatics 7 08-09-2013 08:28 AM
total number of counts with HTseq oliviera Bioinformatics 17 07-26-2013 07:33 AM
htseq-count : 0 counts per miRNA JPerales RNA Sequencing 4 06-04-2013 06:00 AM
Compare RNA counts: HTSeq vs Partek schaffer RNA Sequencing 2 12-02-2011 09:01 AM
understanding HTSeq counts nimmi Bioinformatics 3 11-27-2010 07:24 PM

Reply
 
Thread Tools
Old 01-19-2014, 06:17 PM   #1
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default HTseq: Very few counts recognised

Hi!
Ive seen a lot of threads on this, but I can't figure it out. I got 16-60 millions single end reads in each library. Ive used Tophat 2 with UCSC GTF file for hg19.

This is my code:

samtools view accepted_hits.bam | \
htseq-count -m intersection-nonempty -s no -a 10 \
- UCSC/hg19/genes.gtf \
> Out.txt

Here is a typical result, its propotional to the library size:

no_feature 7013689
ambiguous 269370
too_low_aQual 0
not_aligned 0
alignment_not_unique 6645341

How come i get on average 25 - 50% reads that is "no_feature",
"ambiguous" or "alignment_not_unique".

This is RNAseq, and if I must visually inspect, how to precede?
sindrle is offline   Reply With Quote
Old 01-20-2014, 12:38 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Perhaps you have a lot of immature mRNAs or a lot of expressed repeat regions. The general idea is to look at some of the alignments in IGV and see if they really don't match anything. Also ensure that the chromosome names in the BAM file and GTF file match (that probably causes this sort of thing half the time).
dpryan is offline   Reply With Quote
Old 01-20-2014, 12:38 PM   #3
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Hi again!

I have now tested HTSeq with all modes, also upgraded to Python 2.7.6 and inspected using IGV.

Screen Shot 2014-01-20 at 22.43.27.png

Here is in total 4 reads, one with mapping quality 50 and three with 3.
I used HTSeq option -a 0, so they should been picked up..

All three modes only counts 1 read.. How can this be?
sindrle is offline   Reply With Quote
Old 01-20-2014, 12:42 PM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

HTSeq-count also looks at the NH auxiliary tag. With a MAPQ of 3, it's likely that three of those are multimappers (this will be the case if you used tophat2) and would be (properly) ignored.
dpryan is offline   Reply With Quote
Old 01-20-2014, 01:26 PM   #5
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

Oh yeah, that make sense.

How come you know so much about everything? Where have you learned?

But, its pretty sure something is wrong here right, so I should keep looking? I have checked my GTF, the chromosome names are the same.
sindrle is offline   Reply With Quote
Old 01-20-2014, 05:00 PM   #6
sindrle
Senior Member
 
Location: Norway

Join Date: Aug 2013
Posts: 266
Default

I have made a SAM file with -samout option and checked around a bit..

From Tophat.log I get 21.8 mill. total kept reads
In my SAM I get a total of: 20.90 mill., wonder where the rest are?

Of the 20.9 mill. I have:

17.7 mill. NH:i:1 of which 158.000 ambiguous
I also get 3.4 mill. alignment_not_unique &
1.6 mill. no_feature

Looking at the HTSeq output file I get:

no_feature: 3.8 mill.
ambiguous: 158.000
alignment_not_unique: 3.4 mill.

So the SAM has 1 million reads less than the BAM. Also "no_feature" is different in the SAM and HTSeq output..

I tried to watch specific reads in IGV, but selecting reads by read name (right click the BAM track and choose "select by name", does not change the view....Annoying).

But anyone have something to add on this?
sindrle is offline   Reply With Quote
Reply

Tags
htseq-count

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 08:26 PM.


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