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
                        Strategies for Sequencing Challenging Samples
                        by seqadmin


                        Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                        03-22-2024, 06:39 AM
                      • seqadmin
                        Techniques and Challenges in Conservation Genomics
                        by seqadmin



                        The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                        Avian Conservation
                        Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                        03-08-2024, 10:41 AM

                      ad_right_rmr

                      Collapse

                      News

                      Collapse

                      Topics Statistics Last Post
                      Started by seqadmin, Yesterday, 06:37 PM
                      0 responses
                      10 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, Yesterday, 06:07 PM
                      0 responses
                      10 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 03-22-2024, 10:03 AM
                      0 responses
                      51 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 03-21-2024, 07:32 AM
                      0 responses
                      67 views
                      0 likes
                      Last Post seqadmin  
                      Working...
                      X