Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • htseq-count

    Hi,
    I was wondering about the output of htseq-count to count reads associated with transcripts.

    I run htseq-count to estimate reads associated with genes with the command below:

    htseq-count --stranded=yes --mode=intersection-nonempty file.sam file.gtf >> gene_id.txt

    example I found one gene ENSGxxx with associated 700 reads count.

    To check how many reads are associated with the different transcripts of ENSGxxx gene i run the command below:

    htseq-count --stranded=yes --idattr=transcript_id --mode=intersection-nonempty file.sam file.gtf >> transcript id.txt

    I would expect that the sum of the reads associated with each transcript is almost similar to the reads associated with the gene, but it seems not the case..

    example:
    gene 700 counts
    transcript1 0 counts
    transcript2 4 counts
    transcript3 72 counts


    any suggestions?
    Ciao!
    paolo

  • #2
    I can't say for certain, but this is what I bet is going on.

    Even with the interesection_nonempty option, reads that map equally well to more than will attribute will not be counted and instead will be counted in the ambiguous category.

    The different transcripts are likely going to share the majority of the same exons. So by using transcript_id as the attribute, most of those reads are going to be counted as ambiguous since they could be counted equally for any one of the three transcripts. The counts that you do get for each transcript are those that map to the exons unique to each of those transcripts.

    Look at the end of your count list at the category "ambiguous" I almost guarantee you that the number of counts in this category is vastly more when you use transcript_id versus gene_id.

    Comment


    • #3
      exactly
      Ambiguos reads in transcript id is much more bigger than gene id,

      I tried also to use --mode=union but nothing changed,

      so I think is useless to counts reads associated with transcript with htseq-count,

      any other program could I use?

      Paolo
      Last edited by paolo.kunder; 03-01-2012, 08:16 AM.

      Comment


      • #4
        What do you want to do with it after you have your counts? Is there a specific program you want to use?

        There is RSEM

        Comment


        • #5
          htseq-count is designed to count for genes, not for transcripts. Post #2 correctly explained why the counts are so low.

          Why would you want to count reads per transcript? Given that the vast majority of reads will overlap with more than one transcript, such count values have little to do with transcript expression strengths. (If it were as easy, nobody would bother with developing algorithms for isoform deconvolution.)

          Comment


          • #6
            Thanks for the replies, I'll have a look to RSEM and continue my analysis with DEXSeq and so on...




            Paolo

            Comment


            • #7
              Hello Folks,
              I had a question about HTSeq output which is kind of in tune with this blog....I have sam files as output as well as the count data in text format. So I wanted to look at the no_feature tags that HTSeq puts out to then match with GFF file to see what those no_feature are ( if they were transposons or intronic regions etc).
              My problem is that according to HTSeq count text it generates , i have a value of 33190 for no_feature. However when I either grep for 'no_feature' or reverse grep for 'ENSBTAG' from my sam file to get counts, I get 33055. Why this variation? I would love to get the sequences associated with this no_feature out of sam file to further process it . Any suggestions?
              Thanks
              Geneart.

              Comment


              • #8
                Coincidentally, I spotted a minor bug today when preparing the new release of HTSeq: If a read maps to a chromosome or contig for which no exons at all are listed in the GTF file, the read will be correctly counted as "no_feature" but will not be written to the "samout" file. This is fixed in version 0.6.0, available since an hour or so. Maybe this solves your issue.

                Comment


                • #9
                  Thanks very much!! Will rerun the process with new version.

                  Comment


                  • #10
                    Originally posted by Simon Anders View Post
                    Coincidentally, I spotted a minor bug today when preparing the new release of HTSeq: If a read maps to a chromosome or contig for which no exons at all are listed in the GTF file, the read will be correctly counted as "no_feature" but will not be written to the "samout" file. This is fixed in version 0.6.0, available since an hour or so. Maybe this solves your issue.
                    You are so right that it worked after i changed into exon in .gff. Thanks!
                    wwq

                    Comment


                    • #11
                      Dear All,

                      I am trying to run htseq-count but it keep failing.

                      Script:
                      htseq-count --stranded=no -o samout_test accepted_hits.sam Arabidopsis_thaliana.TAIR10.23.gtf

                      100000 GFF lines processed.
                      200000 GFF lines processed.
                      300000 GFF lines processed.
                      400000 GFF lines processed.
                      485101 GFF lines processed.
                      Error occured when processing SAM input (line 20 of file accepted_hits.sam):
                      ('need more than 10 values to unpack', 'line 20 of file accepted_hits.sam')
                      [Exception type: ValueError, raised in _HTSeq.pyx:1277]

                      I converted the tophat2 file accepted_hits.bam to accepted_hits.sam and used this as input. I am really struggling to understand what this means and how to solve the problem. Any ideas?

                      thanks in advance
                      G

                      Comment

                      Latest Articles

                      Collapse

                      • seqadmin
                        Latest Developments in Precision Medicine
                        by seqadmin



                        Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

                        Somatic Genomics
                        “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
                        Yesterday, 01:16 PM
                      • seqadmin
                        Recent Advances in Sequencing Analysis Tools
                        by seqadmin


                        The sequencing world is rapidly changing due to declining costs, enhanced accuracies, and the advent of newer, cutting-edge instruments. Equally important to these developments are improvements in sequencing analysis, a process that converts vast amounts of raw data into a comprehensible and meaningful form. This complex task requires expertise and the right analysis tools. In this article, we highlight the progress and innovation in sequencing analysis by reviewing several of the...
                        05-06-2024, 07:48 AM

                      ad_right_rmr

                      Collapse

                      News

                      Collapse

                      Topics Statistics Last Post
                      Started by seqadmin, Yesterday, 07:15 AM
                      0 responses
                      13 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 05-23-2024, 10:28 AM
                      0 responses
                      17 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 05-23-2024, 07:35 AM
                      0 responses
                      20 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 05-22-2024, 02:06 PM
                      0 responses
                      10 views
                      0 likes
                      Last Post seqadmin  
                      Working...
                      X