Question: htseq-count for pair-end RNA-seq
0
gravatar for Ming Tang
5.5 years ago by
Ming Tang2.5k
Houston/MD Anderson Cancer Center
Ming Tang2.5k wrote:

Hi there,

I used htseq-count http://www-huber.embl.de/users/anders/HTSeq/doc/count.html to count the pair-end RNA-seq data.

at the end, there are warnings like 

Warning: 8612298 reads with missing mate encountered.

my question is how htseq-count deals with the pairs with only one end existing in the sam file?

does htseq-count simply discard those counts?

Thanks,

Ming

 

rna-seq • 10k views
ADD COMMENTlink modified 4.8 years ago by raphael.poujol10 • written 5.5 years ago by Ming Tang2.5k
are you sure that these reads really don't have a mate in the sam file? If you didn't sort the sam file by read name but you told htseq it is, htseq thinks there are lot of missing mates
ADD REPLYlink modified 5.5 years ago • written 5.5 years ago by Irsan7.0k

I sorted it with samtools srot -n   by name. any ideas? Thank you.

ADD REPLYlink modified 5.5 years ago • written 5.5 years ago by Ming Tang2.5k
1
gravatar for Ido Tamir
5.5 years ago by
Ido Tamir5.0k
Austria
Ido Tamir5.0k wrote:

HtSeq count counts them as 1

ADD COMMENTlink written 5.5 years ago by Ido Tamir5.0k

I know HTseq count one pair of reads as 1, but how about one pair missing the other end? Does HTseq discard it or still count it as 1? Thank you.

ADD REPLYlink written 5.5 years ago by Ming Tang2.5k
2

This is what I meant. My investigations made me believe that each read of a read pair where only one is mapped was counted as 1.

You could check by randomizing your bam file or simply taking the position sorted one. I bet, you will have 2x read count per gene. Please try it and tell me the result. - But I think the last version does not need a read name sorted one anymore, so try it with an old version.

ADD REPLYlink modified 5.5 years ago • written 5.5 years ago by Ido Tamir5.0k

what do you mean by randomizing the bam file? It is now sorted by name.

ADD REPLYlink written 5.5 years ago by Ming Tang2.5k
1

I mean, you can simply test the question how htseq-count counts paired end reads where only one was mapped by sorting this file by position or by randomly sorting it, then doing htseq-count and compare name sorted vs randomly sorted counts. randomly sorted counts should be 2x name sorted counts (in old versions of htseq-count). Or you filter the name sorted by read number (take only 1. read) and then count. Counts should be the same. This should work with all versions of htseq-count and is the better one. Or maybe not, because the "mate aligned flag" is set in those reads.

ADD REPLYlink modified 5.5 years ago • written 5.5 years ago by Ido Tamir5.0k
1
gravatar for raphael.poujol
4.8 years ago by
Canada
raphael.poujol10 wrote:

According to htseq documentation : http://www-huber.embl.de/users/anders/HTSeq/doc/counting.html

"The fact that the records describe the same fragment can be seen from the fact that they have the same read name"

 

So Tophat for instance is ouputing read paris with different names. (Adding 1 or 2 at the end of the name)

So simply do that:

 

(samtools view -H in.bam; in.bam  | awk '{print substr($1,1,length($1)-1),$0}' | sed 's/ [^   ]*//') | samtools view -bSh - > in.samereadname.bam

 

 

 

 

ADD COMMENTlink written 4.8 years ago by raphael.poujol10
1

just a small correction for people passing by, the line should be:
 

(samtools view -H accepted_hits.bam; samtools view accepted_hits.bam | awk '{print substr($1,1,length($1)-1),$0}' | sed 's/ [^   ]*//') | samtools view -bSh - > in.samereadname.bam
ADD REPLYlink written 3.9 years ago by Martombo2.6k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1691 users visited in the last hour