SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
HTSeq-count large proportion of mapped reads with no_feature edm1 Bioinformatics 4 05-26-2014 08:28 AM
"not_aligment" from the output of HTseq-count flyingoyster Bioinformatics 5 09-03-2013 12:21 AM
HTseq count - all reads are "alignment_not_unique" dmsoanes Bioinformatics 2 07-12-2012 05:07 AM
The position file formats ".clocs" and "_pos.txt"? Ist there any difference? elgor Illumina/Solexa 0 06-27-2011 07:55 AM
SEQanswers second "publication": "How to map billions of short reads onto genomes" ECO Literature Watch 0 06-29-2009 11:49 PM

Reply
 
Thread Tools
Old 01-14-2014, 11:03 AM   #1
zzhao2
Member
 
Location: Dallas

Join Date: Aug 2012
Posts: 21
Default htseq-count: many reads are "no_feature".

Hi,
I just started to use htseq-count for my paired-end RNA-Seq data. I've name-sorted uniquely mapped reads outputted by tophat, converted the sorted bam file to sam, and used the sam file with an ensembl gtf file (version 74 for human) as the input to htseq-count. Here's the command and output.
Command:
python -m HTSeq.scripts.count -o out.sam -i gene_name sorted.sam Ensembl_GRCh37.74.gtf > counts.txt

Output:
...
14098588 sam line pairs processed.
no_feature 11732479
ambiguous 52152
too_low_aQual 0
not_aligned 0
alignment_not_unique 0

There are 21150698 reads (some paired, some not) in my sorted sam file. I understand that paired reads are counted as one, and that's why there are only 14098588 processed by htseq-count. Here are my questions:
1) Why are there so many "no_feature" reads? Is this normal for a typical RNA-Seq experiment? Is 11732479 the actual number of reads (paired or not) or the number of single reads + pairs counted as one?
2) I also checked the output.sam but there are only 11585093 reads there. It looks like htseq-count didn't record all the reads that it processed. So how does htseq-count choose which reads to record?

Thank you,
Sylvia
zzhao2 is offline   Reply With Quote
Old 01-14-2014, 11:29 PM   #2
dariober
Senior Member
 
Location: Cambridge, UK

Join Date: May 2010
Posts: 310
Default

Quote:
Originally Posted by zzhao2 View Post
1) Why are there so many "no_feature" reads? Is this normal for a typical RNA-Seq experiment? Is 11732479 the actual number of reads (paired or not) or the number of single reads + pairs counted as one?
Hi- Do the chromosome names in the GTF file match those in the sam file? Ensembl files have chromosome names like 1, 2, 3 etc. while UCSC uses chr1, chr2, chr3 etc. Also, maybe check you are using the same version for both the reference fasta and the gtf (e.g. *both* are hg19).

I guess 11732479 refers to number of fragments so a single end read counts as 1 and the two ends of a paired-end read count 1 as well.

Good luck
Dario
dariober is offline   Reply With Quote
Old 01-15-2014, 03:09 AM   #3
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,133
Default

Quote:
Originally Posted by zzhao2 View Post
Command:
Code:
python -m HTSeq.scripts.count -o out.sam -i gene_name sorted.sam Ensembl_GRCh37.74.gtf > counts.txt
You did not specify a --stranded option for htseq-count which means it is using the default setting, --stranded=yes. There are three possible settings for --stranded (yes, no, reverse). This has to be set properly according to the protocol used to generate your RNA-Seq library. In fact --stranded=yes is probably the least likely. What kit/protocol was used to construct the library in this case?
kmcarr is offline   Reply With Quote
Old 01-16-2014, 07:57 AM   #4
zzhao2
Member
 
Location: Dallas

Join Date: Aug 2012
Posts: 21
Default

Hi Dario,
thanks for your reply. Yes I've noticed the chromosome name issue and changed them to chr1, chr2, etc. Without changing them I could only get a result with all genes' read counts being zero, and all reads as "no_feature".

I didn't use a reference fasta file. Should I use it, and where?

If I have 11732479 "no_feature" fragments, then it's even worse, because the total number of my fragments is 14098588, so only 17% of them are mapped to exons.

Quote:
Originally Posted by dariober View Post
Hi- Do the chromosome names in the GTF file match those in the sam file? Ensembl files have chromosome names like 1, 2, 3 etc. while UCSC uses chr1, chr2, chr3 etc. Also, maybe check you are using the same version for both the reference fasta and the gtf (e.g. *both* are hg19).

I guess 11732479 refers to number of fragments so a single end read counts as 1 and the two ends of a paired-end read count 1 as well.

Good luck
Dario
zzhao2 is offline   Reply With Quote
Old 01-16-2014, 08:05 AM   #5
zzhao2
Member
 
Location: Dallas

Join Date: Aug 2012
Posts: 21
Default

Hi kmcarr,
Thanks for pointing this out. I didn't specify the stranded option because my reads are stranded. I've actually tried all three settings: yes, no, and reverse, but none of them gave a number of "no_feature" reads that's significantly less than others.

I don't actually know the kit used, because what I have is only the data from our sequencing facility, but their report says that my reads are stranded.

Quote:
Originally Posted by kmcarr View Post
You did not specify a --stranded option for htseq-count which means it is using the default setting, --stranded=yes. There are three possible settings for --stranded (yes, no, reverse). This has to be set properly according to the protocol used to generate your RNA-Seq library. In fact --stranded=yes is probably the least likely. What kit/protocol was used to construct the library in this case?
zzhao2 is offline   Reply With Quote
Old 01-16-2014, 08:23 AM   #6
bob-loblaw
Member
 
Location: /home/bob

Join Date: Jun 2012
Posts: 59
Default

I wonder though in everyone's experience, what's an "acceptable" percent of no features? Surely no organism's annotation is perfect? Plus how often is it that the RNA that you get is 100% mRNA

In this nature paper they found about 86% mapped to known exons http://www.nature.com/nature/journal...ture08872.html
bob-loblaw is offline   Reply With Quote
Reply

Tags
htseq-count, rna-seq

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 06:52 AM.


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