Question: Inconsistency between htseq-count output sorted by position or by name
1
gravatar for lolec
2.8 years ago by
lolec10
France
lolec10 wrote:

Hello everyone!

I am using htseq-count (version 0.6.1p1) on paired-end RNA seq data. I did the analysis sorting by name and sorting by position, however I obtained different results.

By name:

samtools sort -n accepted_hits.bam -o accepted_nsorted.bam
python -m HTSeq.scripts.count -f bam -r name -s reverse accepted_nsorted.bam GRCh37.75.gtf > htseq_accepted_nsorted.txt

tail htseq_accepted_nsorted.txt
ENSG00000273487 2
ENSG00000273488 2 
ENSG00000273489 3
ENSG00000273492 5
ENSG00000273493 1
__no_feature    1372372
__ambiguous     820117
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  1836341

By position:

samtools sort accepted_hits.bam -o accepted_psorted.bam
python -m HTSeq.scripts.count -f bam -r pos -s reverse accepted_psorted.bam GRCh37.75.gtf > htseq_accepted_psorted.txt

Message:

Warning: Mate pairing was ambiguous for 23096 records; mate key for first such record: ('HWI-1KL149:106:C8CUGACXX:5:2102:5222:87333', 'first', '1', 146410, '1', 146469, 9416).

29201104 SAM alignment pairs processed.

tail htseq_accepted_psorted.txt
ENSG00000273487 2
ENSG00000273488 2
ENSG00000273489 3
ENSG00000273492 5
ENSG00000273493 1
__no_feature    1467721
__ambiguous     866870
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  2028312

What is wrong in this analysis ? How can we explain the different counts obtained ? Is there a solution to this problem ?

rna-seq htseq-count • 1.4k views
ADD COMMENTlink written 2.8 years ago by lolec10

Looks weird. Probably related to the warning you received.

See htseq-count manual:

-r <order>, --order=<order>¶ For paired-end data, the alignment have to be sorted either by read name or by alignment position. If your data is not sorted, use the samtools sort function of samtools to sort it. Use this option, with name or pos for <order> to indicate how the input data has been sorted. The default is name.

If name is indicated, htseq-count expects all the alignments for the reads of a given read pair to appear in adjacent records in the input data. For pos, this is not expected; rather, read alignments whose mate alignment have not yet been seen are kept in a buffer in memory until the mate is found. While, strictly speaking, the latter will also work with unsorted data, sorting ensures that most alignment mates appear close to each other in the data and hence the buffer is much less likely to overflow.

Makes me guess that the names of your read (pairs) are not unique? Just something I would check. How did you obtain the bam?

ADD REPLYlink written 2.8 years ago by WouterDeCoster39k

There is no 'warning' when sorted by name. Then I don't thing the problem is coming from reads' name. The bam was obtained with TOPHAT2.

ADD REPLYlink written 2.8 years ago by lolec10

I wouldn't expect a warning when sorting by name, since in that case htseq-count expects the both pairs to be adjacent. If sorted by position, htseq-count will search for the other read of the pair (to perform the mate pairing). And this gave a warning as you can see. Maybe it's the best to contact the author of the tool.

ADD REPLYlink written 2.8 years ago by WouterDeCoster39k

Okay, I will do that. Thank you for your help.

ADD REPLYlink written 2.8 years ago by lolec10

Did you manage to find any solution? (thanks in advance)

ADD REPLYlink written 2.4 years ago by al-ash100

Unfortunately no. It might be inherent to how HTSEQ-count works.

ADD REPLYlink written 2.4 years ago by lolec10
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: 844 users visited in the last hour