HtSeq counts 0 features for mapped reads in sam fliles

[rags_brown]

Hello,
My input :
50 bp RNA Seq single end illumina reads mapped to GRCh 38 from ensemble.
Mapping was performed both by STAR and bowtie-2.For star almost 70% reads were uniquely mapped and I used HtSeq to count the features:
python -m HTSeq.scripts.count -t exon -m union -s no 1.sam $ref > 1.sam_1.counts

and for the gene "ENSG00000100219" I got a feature count of 0. I tried to see if there was any mapping for this gene

samtools view 1.bam 22:28794555-28800597 | wc -l

and found that there were 2899 reads mapped.

I tried counting just the transcript and the gene to no avail.

I am quite not sure why this would be the case? Please let me know if my understanding is right. Is there an issue with the alignment or htseq flags?

[alexdobin]

Hello,
My input :
50 bp RNA Seq single end illumina reads mapped to GRCh 38 from ensemble.
Mapping was performed both by STAR and bowtie-2.For star almost 70% reads were uniquely mapped and I used HtSeq to count the features:
python -m HTSeq.scripts.count -t exon -m union -s no 1.sam $ref > 1.sam_1.counts

and for the gene "ENSG00000100219" I got a feature count of 0. I tried to see if there was any mapping for this gene

samtools view 1.bam 22:28794555-28800597 | wc -l

and found that there were 2899 reads mapped.

I tried counting just the transcript and the gene to no avail.

I am quite not sure why this would be the case? Please let me know if my understanding is right. Is there an issue with the alignment or htseq flags?

Hi,

I can think of a few of reasons HTseq might not count these reads: (i) they are multimappers, (ii) they are not overlapping exons of the gene (i.e. they are all intronic) (iii) they overlap exons of more than one gene.

If you make a new BAM file with just the selected reads, using
samtools view -b 1.bam 22:28794555-28800597 > 2.bam
and then run HTseq with 2.bam, it will tell you how it classified the unassigned reads (after all the gene counts):
__no_feature --> (ii)
__ambiguous --> (iii)
__too_low_aQual
__not_aligned
__alignment_not_unique -->(i)

Also, you can check for (i) by looking at the NH:i:<Nmult> tag of the reads selected over the gene locus.
To check for (ii) or (iii), you can load your alignments into a browser and see if any of them overlap exons of the gene.

Cheers
Alex

[Michael.Ante]

Hi,

next to Alex' suggestions, you might also check the chromosomes' name in your gtf-file. Sometimes they come as e.g. 'chr22' while your bam-file has a '22' instead.

Cheers,
Michael