[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

[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.

[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]

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

[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?

[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

[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

[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

[Simon Anders]

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

[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]

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.

[zaki]

Hi members,

Reviving an old topic... I am bewildered about htseq-count output. The endpoint is I will be using DEseq for differential expression, but I am confused with the htseq-count results.

I compared the htseq-count results with what I view in IGV... and I see very inconsistent things which I can not explain.. I was hoping if anyone can give an insight into this matter..

Quick overview on what I did
Mapped the paired end .fastq file to a reference transcription using tophat2
[CODE]"tophat2 -G <transcript.gtf> -o <outputdir> --no-novel-juncs <bowtie2 indexed genome> <.fastq1> <.fastq2> [CODE]

Following this I converted the resulting .bam file to .sam in order to proceed with read counting using htseq-count
#Converting .bam to .sam

samtools view -h <bam> > <sam>

#Sorting the .sam (Would be great if anyone can tell me if this is the correct method or not, as sometimes it returns some error - I also asked here http://seqanswers.com/forums/showthr...t=12376&page=2)

sort -k1 <sam> > <sorted.sam>

#Running htseq-count

htseq-count -m intersection-strict -s no <sorted.sam> <transripts.gtf> > <result_count>

Above is the general code I used, I also tried running htseq-count with -m intersection-nonempty.

The problem I have now is when I look at the gene count and IGV it gives confusing statement

Location 1
Result from htseq-count gives - 14 (regardless which htseq-count mode)
Image from IGV


Location 2
Result from htseq-count gives - 7 (regardless which htseq-count mode)
Image from IGV


As you can see, at Location 2 has more reads, but how come the counts number only is 7??

Could it be to do with the .gtf file I am using as mentioned by

i'd guess it's due to the fact that htseq-count only reports one hit per aligned read assuming that read is aligned uniquely and does not overlap multiple features listed in your GTF. if an aligned read hits more than one feature in your GTF then it doesn't report that hit. bedtools gives you raw hits which includes every 1 hit for every intersection of every alignment with any features in the GTF no matter how many times it aligned or how many features it hit. you might think, "wow, htseq-count is dropping a lot of information". yes, it is! i've moved to using other tools to count hits to genes (RSEM/eXpress) since they disambiguate those ambiguous alignments and as a result you get counts for all of your aligned reads. in a genome with alternative splicing you lose too much data using htseq-count, in my opinion.

Or its because location 2 has tons of non-unique mapping and the read counts are very low??

Eventually if I would be doing the differential expression on DESeq, would this give me a false result at location 2 since there is only 7 count but there are much more reads??

Appreciate any feedback

Can anyone help explain why is there this discrepancy between IGV and count data??

[jingerlu]

usually , I go to ensemble(www.ensembl.org) => bioMart . CHOOSE DATABSE ( which is Ensemble gene 70 currently ) , CHOOSE DATASET (which is Homo sapiens for Human), once these two steps done , the left side of the webpage
will show Dataset / filter / Attributes , click Attributes (check Features) , expand "gene" , check "Ensembl gene id" "Associated gene name" "description" or other features you want . Click "Result " which is at the top left side of the webpage. you can "Export all result to " there .

How can I convert the ENSG00000261838 to real name of genes? Or is there any parameter do I need to use prior using htseq-count? Thank you

[wanfahmi]

How can I convert the ENSG00000261838 to real name of genes? Or is there any parameter do I need to use prior using htseq-count? Thank you

ENSG00000261838 0
ENSG00000261839 2
ENSG00000261840 3
ENSG00000261841 0
ENSG00000261842 0
no_feature 985727
ambiguous 435068
too_low_aQual 0
not_aligned 0
alignment_not_unique 778451

[sdriscoll]

That explains it! Thanks for pointing that out. It's unfortunate that the annotations are so lame. Cufflinks can actually help with this but people still will need to parse the gtf to build a translation table. You can use cuffcompare to compare the gtf to itself and that produces cuffcmp.combined.gtf. In that file gene ids are replaced with XLOC ids which are unique to loci. The gene names are set in the gene_name field, if they were in your gtf to start. You can also use the gffread program that comes packaged with Cufflinks. It has a --cluster-only option that adds a locus attribute to each gtf row grouping the features into unique loci. One could specify that attribute when running htseq-count and then translate the locus ids bak to gene names. I think ideally the counting script should do the clustering or at least provide the option. Then the annotation doesn't have to be perfect.