Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • HTseq-count : how to get unique reads?

    Hello,

    I used HTseq-count to get number of reads / gene from my RNA-seq data.

    The problem is that the sum of all counted reads (on every gene) is much bigger then the total sequencing depth I got from Illumina. After reading other posts, this could be due to the "duplicate" alignemnt... Anyway, my question is: how to get the true number of reads for some genes (in my case rRNA) in order to see what percentage of my total sequencing depth is actually "lost" on rRNA?

    Thanks in advance for any hints,

    TP

  • #2
    Perhaps, I addressed the question in the wrong way. Simply, i just wanna know how to get the Ribosomal RNA contamination in each of my libraries? I want to be able to say that 1E07 reads out of 1E09 are mapping to rRNA, so my samples are "contaminated" by 1% rRNA. See?

    Thanks

    Comment


    • #3
      This will require additional work but could get rid of the contamination in the process.
      http://seqanswers.com/forums/showthread.php?t=26960 Use the rRNA sequence for the species you are working with.

      Comment


      • #4
        Yeah I saw that thread... I was hoping for faster solution, but well
        Thanks

        Comment


        • #5
          htseq-count should identify the number of multimapped reads and include this number in its own category so that the total number should not exceed your starting number of reads....unless, whatever you mapped the reads with does not output a NH:i field in the SAM file.

          What did you use to map the reads?

          Are you specifically interested in only rRNA reads...not everything that maps multiple times will be an rRNA read.

          Comment


          • #6
            Originally posted by chadn737 View Post
            htseq-count should identify the number of multimapped reads and include this number in its own category so that the total number should not exceed your starting number of reads....unless, whatever you mapped the reads with does not output a NH:i field in the SAM file.

            What did you use to map the reads?

            Are you specifically interested in only rRNA reads...not everything that maps multiple times will be an rRNA read.
            I used bowtie for mapping. After that, I used Htseq-count to get number of reads/gene. If I sum all the reads for one replicate from Htseq table I get like double of what I'm supposed to have. I guess it's because of all multimapped reads. I will try to align my libraries only to rRNA, that should give me number of reads mapping there.

            TP

            Comment


            • #7
              Originally posted by ThePresident View Post
              I used bowtie for mapping. After that, I used Htseq-count to get number of reads/gene. If I sum all the reads for one replicate from Htseq table I get like double of what I'm supposed to have. I guess it's because of all multimapped reads. I will try to align my libraries only to rRNA, that should give me number of reads mapping there.

              TP
              Bowtie does not output a NH:i field in the SAM file and so htseq-count will not be able to identify multimapping reads.

              Furthermore, Bowtie is not ideally suited for RNA-seq data. You should use Tophat instead, Tophat puts out a NH:i field, and so then your htseq-count results will make sense. I would not use htseq-count on data aligned with bowtie unless multi-mapped reads were first removed or identified. Otherwise, your read counts for each gene will be incorrect. Furthermore, bowtie is not a spliced aligner and will preferentially align to the best match, even if its a pseudogene. With tophat, you can deal with splicing. Unless you have good reason to use bowtie, I would realign everything with tophat instead and problem solved.

              Also, as I said previously, not all multimapping reads will be due to rRNA, so if you are trying to identify/eliminate multimapping reads, aligning to rRNA does not solve your problem.

              Depending on whether you used bowtie or bowtie 2 (and if bowtie 2, what mode you ran it in) it is possible to remove all multimapping reads and get the unique reads only by filtering on quality scores.

              Comment


              • #8
                Well, I'm dealing with bacterial RNA-seq data, so a spliced aligner such as Tophat was irrelevant in this case. However, I didn't realize that bowtie is messing up with multimapped reads

                I think that overall my alignment was not so bad, 'cause many genes identified as transcriptionally altered in my RNA-seq were confirmed by qPCR. However, I might try tophat as well.

                About rRNA, it's just that I want to present some data in a paper I'm working on, and would like to give the number of reads that aligned on rRNA in my libraries. That way, it would give the audience some idea about abundance of "real" transcripts in my libraries.

                Comment


                • #9
                  Originally posted by ThePresident View Post
                  Well, I'm dealing with bacterial RNA-seq data, so a spliced aligner such as Tophat was irrelevant in this case. However, I didn't realize that bowtie is messing up with multimapped reads

                  I think that overall my alignment was not so bad, 'cause many genes identified as transcriptionally altered in my RNA-seq were confirmed by qPCR. However, I might try tophat as well.

                  About rRNA, it's just that I want to present some data in a paper I'm working on, and would like to give the number of reads that aligned on rRNA in my libraries. That way, it would give the audience some idea about abundance of "real" transcripts in my libraries.
                  There seems no reason to use tophat then.

                  What version of bowtie did you use and how did you use it. If you want to present data on rRNA, you can align specifically to these, however, to quickly identify uniquely mapped reads, you can filter based on quality score depending on what version of bowtie you used and how you ran it.

                  Comment


                  • #10
                    I used bowtie (not bowtie2) with following command line:

                    Code:
                    bowtie -q -a --best –m50 –p6 -t index 6_Index-1.WT_1_R1.fastq -S AlignWT1
                    So, I used -a for bowtie to report all alignments and --best to give only best hits.

                    Now, how you would filter for quality scores with bowtie?

                    Comment


                    • #11
                      Originally posted by ThePresident View Post
                      I used bowtie (not bowtie2) with following command line:

                      Code:
                      bowtie -q -a --best –m50 –p6 -t index 6_Index-1.WT_1_R1.fastq -S AlignWT1
                      So, I used -a for bowtie to report all alignments and --best to give only best hits.

                      Now, how you would filter for quality scores with bowtie?
                      By default, bowtie should assign a mapping quality of 255 to reads that map once. You can filter using samtools view:

                      Code:
                      samtools view -hb -q 255 input.bam > output.bam
                      Your output should have only uniquely aligned reads.

                      Comment


                      • #12
                        Thanks, I'll try that instead!

                        Comment

                        Latest Articles

                        Collapse

                        • seqadmin
                          Non-Coding RNA Research and Technologies
                          by seqadmin




                          Non-coding RNAs (ncRNAs) do not code for proteins but play important roles in numerous cellular processes including gene silencing, developmental pathways, and more. There are numerous types including microRNA (miRNA), long ncRNA (lncRNA), circular RNA (circRNA), and more. In this article, we discuss innovative ncRNA research and explore recent technological advancements that improve the study of ncRNAs.

                          Nobel Prize for MicroRNA Discovery
                          This week,...
                          10-07-2024, 08:07 AM
                        • seqadmin
                          Recent Developments in Metagenomics
                          by seqadmin





                          Metagenomics has improved the way researchers study microorganisms across diverse environments. Historically, studying microorganisms relied on culturing them in the lab, a method that limits the investigation of many species since most are unculturable1. Metagenomics overcomes these issues by allowing the study of microorganisms regardless of their ability to be cultured or the environments they inhabit. Over time, the field has evolved, especially with the advent...
                          09-23-2024, 06:35 AM

                        ad_right_rmr

                        Collapse

                        News

                        Collapse

                        Topics Statistics Last Post
                        Started by seqadmin, 10-02-2024, 04:51 AM
                        0 responses
                        103 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 10-01-2024, 07:10 AM
                        0 responses
                        112 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 09-30-2024, 08:33 AM
                        1 response
                        114 views
                        0 likes
                        Last Post EmiTom
                        by EmiTom
                         
                        Started by seqadmin, 09-26-2024, 12:57 PM
                        0 responses
                        21 views
                        0 likes
                        Last Post seqadmin  
                        Working...
                        X