[zaki]

Please take a few of these many reads in location 2, hover your mouse over them in IGV to get the read name, then grep for the read name to find the corresponding lines in the SAM file. Maybe they are all non-unique mapping reads.

The fact that htseq-count does not count these reads is deliberate, of course. See earlier threads for an explanation.

Many thanks Simon! I think I need to get my basic understanding of mapping solid.

Forgive my naivety, If multi-mapping becomes an issue when we are translating mapping reads to count data (expression level)...what than would be an acceptable standard to extract count data?? (is it using htseq-count ,bedtools, cuffdiff)

Also....during the mapping stage..what bias would be introduced, if we discard multimap reads?


Edit - you are correct in saying that the reads at location 2 is multimapped, many thanks for the insight!

[Simon Anders]

See my post #4 in this thread: http://seqanswers.com/forums/showthread.php?t=9129

[lpachter]

Dear Zaki,

You ask an excellent question. This paper answers it in detail:

http://www.nature.com/nbt/journal/v31/n1/abs/nbt.2450.html

Lior

[zaki]

Thanks Simon and Lior for the suggested readings!

Taking a step back and reading the fundamentals behind each approach is really essential. I think its something I overlooked at first.

Since the protocol and vignette of each approach is very well described and easy to follow, I guess I jumped the gun and started running commands without understanding much behind it (I think I still have a lot yet to understand... statistics...mapping..biology...etc...)

Thanks again for your comments and direction! It helped a lot.



Zaki

[Jose Garcia]

Please take a few of these many reads in location 2, hover your mouse over them in IGV to get the read name, then grep for the read name to find the corresponding lines in the SAM file. Maybe they are all non-unique mapping reads.

The fact that htseq-count does not count these reads is deliberate, of course. See earlier threads for an explanation.

Indeed. I was having the same problem. I was getting crazy because I had already checked what Simon pointed about the possibility of having non-unique mapping reads. I think I have found out why.
If you use htseq-count with -i gene_id, it will construct the set where to add up reads based on the ENSEMBL GENE ID, which will sum reads that could come from different ENSEMBL transcript ID. You may lose the information about transcripts at this level (maybe you can use next DEXSeq approach to look inside...) but at least, you do not lose reads mapping to different identifiers(transcript id), as was my case because I used the UCSC gtf file where gene_id=trascritpt_id.

However, I still encountered IGV regions in my BAM files with many reads that htseq count anyhow flagged as ambigous, but were not multi aligned. I think I have found out why by looking at the locus with the ensembl genome browser and looking for the "gene_name", that is, the SYMBOL. I realized that the Ensembl annotation (present in the GTF) contained for the same locus, with the same gene_name (SYMBOL), two different ENSEMBLGENE Ids for what it called: "protein coding" and "non-sense-mediated decay", as gene biotype. Both lying in the same locus with high overlap. I guessed then that htseq-count had found reads aligning to both ENSEMBLGENE identifiers and discarded them as ambiguous.

While discarding reads mapping to different locus or different genes which overlap (hence also different gene Symbols) may have sense, I do not see the sense in discarding reads merely because we have a maybe too exhaustive overlapping annotation. We might also have other features like this, antisense ncRNAs, pseudogenes, for instance, which may have a big overlapping with its parental, and hence, all counts in regions with the same sequence will be discarded.
For the pseudogenes, unfortunately, if they have a huge overlap, all those reads will be discarded by htseq-count -i gene_name either, but if we are not interested in looking at pseudogenes we could choose to use a more limited annotation like RefSeq.

So, what do you think of using as a global strategy with htseq-count when interested in having also all non coding, pseudogenes et. al a command line like this:

htseq-count -m union -t exon -i gene_name - $GTF > $out

It is always possible to try and deconvolve counts thereafter,
What do you think about it?
Thanks
Best
Jose

[vyellapa]

In order to get gene level expression which is a better tool to use? HTSeq or Bedtools. Or has anyone tried any other tool for this purpose?

[dpryan]

Don't use bedtools, the results will be wrong for what you want to do. Use htseq-count or featureCounts (from subread).

[vyellapa]

Thank you for your reply.
For gene expression measurement, htseq values I believe could be conservative as only unique mapping events are used along with non-overalapping reads. Feature counts is fast and through there is substantial concordance of count values with htseq count. About 20 features per sample have count differences of greater than 10000.

Which one does one believe?

[dpryan]

As a general rule, htseq-count. You might use the debug options ("-o" for htseq-count, I don't recall what it is for featureCounts) to see why there's such a big difference between the two. Given the same thresholds (i.e., MAPQ filtering) and intersection model, the counts should be the same.

[shi]

Thank you for your reply.
For gene expression measurement, htseq values I believe could be conservative as only unique mapping events are used along with non-overalapping reads. Feature counts is fast and through there is substantial concordance of count values with htseq count. About 20 features per sample have count differences of greater than 10000.

Which one does one believe?

The following paper explains in details the difference between featureCounts and htseq-count:

http://www.ncbi.nlm.nih.gov/pubmed/24227677