samtools sorting issue or HTSeq-count problem?

[bbl]

Hello all,

This has puzzled me for days and I couldnt find any explanation on the internet.
I noticed HTSeq-count gives the reads counts of a gene as 0. But when viewed in IGV with the accepted_hits.bam from tophat2 alignment, there are hundreds to thousands reads aligned in that gene region (from different samples).

This is how I used HTSeq-count from HTSeq/0.6.1p1 and samtools/1.1
samtools sort -T tmp -o accepted_hits_nsort.bam -n accepted_hits.bam
htseq-count --format=bam --strand=no --order=name -a 5 --mode=union accepted_hits_nsort.bam Homo_sapiens.GRCh37.72.gtf > htseq_count_0.6.txt

However, when I break down the accepted_hits.bam and extract the chromosome where the gene is, HTSeq-count gives the right reads counts for this BAM file.

My accepted_hits.bam contains 60,845,216 alginemnt using samtools view -c accepted_hits.bam. I dont think for PE RNA-seq, file at size of 4.5G would be a problem for samtools or HTSeq-count to deal with?

My suspicions are:
1. samtools name sort problem
2. htseq-count bug
3. out of memory

Could any one give any hint here? Any suggestion would be appreciated. Thanks.

[dpryan]

Try using the -o option and see what happens to some of those reads that should be getting counted.

BTW, featureCounts is MUCH faster, which is especially convenient on larger files.

[kmcarr]

What gene are you looking at in IGV? The most likely answer is that you are looking at a duplicated gene and the 100's to 1000's of reads which map to it also map to the other copies elsewhere in the genome. When a read maps to multiple locations HTSeq-count will not count it for any of the genes it maps to, rather it will be be counted among the 'alignment_not_unique'. Your experiment of running HTSeq-count on just one chromosome shows that the additional mappings for these reads fall on other chromosomes. When HTSeq-count is given an incomplete alignment set with only one valid mapping for these reads it will naturally count them for that one alignment. This demonstrates why it is a bad idea to perform read counting on only partial alignment sets, unless you are prepared to deal with these types of situations.

If you want to confirm this identify one of the reads aligned to the gene you have identified in IGV and then search for that read name in your accepted_hits.bam to see if it is multiply mapped.

[bbl]

Try using the -o option and see what happens to some of those reads that should be getting counted.

BTW, featureCounts is MUCH faster, which is especially convenient on larger files.

Thanks very much dpryan. I shall try featureCounts too.

[bbl]

Thanks very much kmcarr. Indeed what you suspected is what happens!
1st, CEBPA is the gene I'm intereted. From samout by HTSeq-count, reads are mapped to two and three different genes incl. CEBPA. Therefore, they are counted as ambiguous. Actually these reads are also mapped to 1. 'retired' gene, 2. antisense in my .gff ref downloaded from Ensembl.
Since I'm looking at RNAseq samples treated using polyA+tail , I suppose it makes sense to 'manually' modify the .gff ref? For example, remove those pseudogenes, non-coding RNA etc...Would you suggest some way to modify it to be more suitable, if it is a right thing to do?

Many thanks again.