htseq-count for pair-end RNA-seq

Question

[Ming Tommy Tang]

Hi there,

I used htseq-count 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

[Irsan]

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

[Ming Tommy Tang]

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

Answer 1

[Ido Tamir]

HtSeq count counts them as 1

[Ming Tommy Tang]

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.

[Ido Tamir]

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.

[Ming Tommy Tang]

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

[Ido Tamir]

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.

Answer 2

[raphael.poujol]

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

[Martombo]

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