Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • samtools sorting issue or HTSeq-count problem?

    Hello all,

    This has puzzled me for days and I couldnt find any explanation on the internet.
    I noticed HTSeq-count gives the reads counts of a gene as 0. But when viewed in IGV with the accepted_hits.bam from tophat2 alignment, there are hundreds to thousands reads aligned in that gene region (from different samples).

    This is how I used HTSeq-count from HTSeq/0.6.1p1 and samtools/1.1
    samtools sort -T tmp -o accepted_hits_nsort.bam -n accepted_hits.bam
    htseq-count --format=bam --strand=no --order=name -a 5 --mode=union accepted_hits_nsort.bam Homo_sapiens.GRCh37.72.gtf > htseq_count_0.6.txt

    However, when I break down the accepted_hits.bam and extract the chromosome where the gene is, HTSeq-count gives the right reads counts for this BAM file.

    My accepted_hits.bam contains 60,845,216 alginemnt using samtools view -c accepted_hits.bam. I dont think for PE RNA-seq, file at size of 4.5G would be a problem for samtools or HTSeq-count to deal with?

    My suspicions are:
    1. samtools name sort problem
    2. htseq-count bug
    3. out of memory

    Could any one give any hint here? Any suggestion would be appreciated. Thanks.

  • #2
    Try using the -o option and see what happens to some of those reads that should be getting counted.

    BTW, featureCounts is MUCH faster, which is especially convenient on larger files.

    Comment


    • #3
      What gene are you looking at in IGV? The most likely answer is that you are looking at a duplicated gene and the 100's to 1000's of reads which map to it also map to the other copies elsewhere in the genome. When a read maps to multiple locations HTSeq-count will not count it for any of the genes it maps to, rather it will be be counted among the 'alignment_not_unique'. Your experiment of running HTSeq-count on just one chromosome shows that the additional mappings for these reads fall on other chromosomes. When HTSeq-count is given an incomplete alignment set with only one valid mapping for these reads it will naturally count them for that one alignment. This demonstrates why it is a bad idea to perform read counting on only partial alignment sets, unless you are prepared to deal with these types of situations.

      If you want to confirm this identify one of the reads aligned to the gene you have identified in IGV and then search for that read name in your accepted_hits.bam to see if it is multiply mapped.

      Comment


      • #4
        Originally posted by dpryan View Post
        Try using the -o option and see what happens to some of those reads that should be getting counted.

        BTW, featureCounts is MUCH faster, which is especially convenient on larger files.
        Thanks very much dpryan. I shall try featureCounts too.

        Comment


        • #5
          Thanks very much kmcarr. Indeed what you suspected is what happens!
          1st, CEBPA is the gene I'm intereted. From samout by HTSeq-count, reads are mapped to two and three different genes incl. CEBPA. Therefore, they are counted as ambiguous. Actually these reads are also mapped to 1. 'retired' gene, 2. antisense in my .gff ref downloaded from Ensembl.
          Since I'm looking at RNAseq samples treated using polyA+tail , I suppose it makes sense to 'manually' modify the .gff ref? For example, remove those pseudogenes, non-coding RNA etc...Would you suggest some way to modify it to be more suitable, if it is a right thing to do?

          Many thanks again.

          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
          8 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, Yesterday, 06:07 PM
          0 responses
          8 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 03-22-2024, 10:03 AM
          0 responses
          49 views
          0 likes
          Last Post seqadmin  
          Started by seqadmin, 03-21-2024, 07:32 AM
          0 responses
          66 views
          0 likes
          Last Post seqadmin  
          Working...
          X