Seqanswers Leaderboard Ad

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts
  • jwhite
    Member
    • Jun 2012
    • 33

    flagstat

    Hi Folks,

    Our standard process for preparing BAM files involves the following:

    samtools sort
    samtools rmdup
    samtools sort
    samtools index

    20682963562 Mar 19 20:13 EAP120_R2.fastq.bam
    12717567745 Mar 19 21:50 EAP120_R2.fastq.srt.bam
    4529263043 Mar 19 22:11 EAP120_R2.fastq.nodup.bam
    4529263062 Mar 19 22:38 EAP120_R2.fastq.nodup.srt.bam

    So there will be 4 BAM files from this process. To evaluate them, I used samtools flagstat.
    What I found with one of our samples is that the number of mapped reads drops from 17.8 billion to 5.8 billion. (The flagstat data are below.)

    Why would that occur? Does that indicate poor quality data?
    It would seem to me to indicate that there are a LOT of duplicates. Is that correct? Should I normally see a 75 % reduction in the number of mapped reads after duplicate removal?

    So, I guess what I'm asking is how should I interpret flagstat results?

    Thanks for any light you can shed on this.

    Joe White

    -----------------------
    # initial BAM file
    samtools flagstat EAP120_R2.fastq.bam
    177833072 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
    172622167 + 0 mapped (97.07%:-nan%)
    177833072 + 0 paired in sequencing
    88916536 + 0 read1
    88916536 + 0 read2
    169943152 + 0 properly paired (95.56%:-nan%)
    171556513 + 0 with itself and mate mapped
    1065654 + 0 singletons (0.60%:-nan%)
    1376256 + 0 with mate mapped to a different chr
    1254219 + 0 with mate mapped to a different chr (mapQ>=5)

    # after the first sort
    samtools flagstat EAP120_R2.fastq.srt.bam
    177833072 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
    172622167 + 0 mapped (97.07%:-nan%)
    177833072 + 0 paired in sequencing
    88916536 + 0 read1
    88916536 + 0 read2
    169943152 + 0 properly paired (95.56%:-nan%)
    171556513 + 0 with itself and mate mapped
    1065654 + 0 singletons (0.60%:-nan%)
    1376256 + 0 with mate mapped to a different chr
    1254219 + 0 with mate mapped to a different chr (mapQ>=5)

    # after duplicate removal
    samtools flagstat EAP120_R2.fastq.nodup.bam &
    57792684 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
    52581779 + 0 mapped (90.98%:-nan%)
    57792684 + 0 paired in sequencing
    28898595 + 0 read1
    28894089 + 0 read2
    50015042 + 0 properly paired (86.54%:-nan%)
    51516125 + 0 with itself and mate mapped
    1065654 + 0 singletons (1.84%:-nan%)
    1376256 + 0 with mate mapped to a different chr
    1254219 + 0 with mate mapped to a different chr (mapQ>=5)

    # after the second sort
    samtools flagstat EAP120_R2.fastq.nodup.srt.bam &
    57792684 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
    52581779 + 0 mapped (90.98%:-nan%)
    57792684 + 0 paired in sequencing
    28898595 + 0 read1
    28894089 + 0 read2
    50015042 + 0 properly paired (86.54%:-nan%)
    51516125 + 0 with itself and mate mapped
    1065654 + 0 singletons (1.84%:-nan%)
    1376256 + 0 with mate mapped to a different chr
    1254219 + 0 with mate mapped to a different chr (mapQ>=5)
  • swbarnes2
    Senior Member
    • May 2008
    • 910

    #2
    You don't need to sort after the rmdup.

    But yes, sometimes a library really is that repetitive. Maybe you have way more reads than your small genome needs. Maybe your library prep got screwed up.

    Rule #1 Don't disbelieve your data just because you don't like what its telling you.

    If you like, use Picard MarkDuplicates instead, because it will only mark, and not remove the duplicates, and then you can eyeball and see the repetition. Or look at the first sort file in IGV, or something like that, and see if you can see the duplicates.

    Comment

    • mmmm
      Senior Member
      • Jul 2013
      • 131

      #3
      using flagstat- what does this mean?
      "76 reads have mapped to a different chromosome than their pair" (BWA was used for reads mapping to the reference bacterial genome)

      Comment

      • dpryan
        Devon Ryan
        • Jul 2011
        • 3478

        #4
        The origin of that is comparison of the RNAME and RNEXT columns. If they're not the same, then this is incremented. I would have to presume that the bacterial genome you aligned against contained not just the main genome but also some plasmids (or perhaps it was all just a bunch of contigs).

        Comment

        • mmmm
          Senior Member
          • Jul 2013
          • 131

          #5
          yes- there is a plasmid but it showed only 1 SNP

          I have got these stats for mapping to the chromosome and plasmid- does "76 + 0 with mate mapped to a different chr" refere to repitative sequences that are shared between plasmid and chromosome

          3753203 + 0 in total (QC-passed reads + QC-failed reads)
          0 + 0 duplicates
          3660519 + 0 mapped (97.53%:-nan%)
          3753203 + 0 paired in sequencing
          1876530 + 0 read1
          1876673 + 0 read2
          3644736 + 0 properly paired (97.11%:-nan%)
          3653854 + 0 with itself and mate mapped
          6665 + 0 singletons (0.18%:-nan%)
          76 + 0 with mate mapped to a different chr
          74 + 0 with mate mapped to a different chr (mapQ>=5)

          Comment

          • dpryan
            Devon Ryan
            • Jul 2011
            • 3478

            #6
            Given that 74/76 have MAPQ>=5 it would seem that these are unlikely to be that repetitive. However, they might be technical artefacts, since there are so few of them.

            Comment

            • mmmm
              Senior Member
              • Jul 2013
              • 131

              #7
              thanks- if I want to determine the (0.18%) of unmapped reads- what is the simplest way for this?- should I use velvet for denovo assembly (kmer 31) then annotate the output?

              Comment

              • dpryan
                Devon Ryan
                • Jul 2011
                • 3478

                #8
                ~2.5% and yeah that'd be one approach. You might have a look at the reads first to see if they just didn't align due to low quality or something like that.

                Comment

                • mmmm
                  Senior Member
                  • Jul 2013
                  • 131

                  #9
                  thank you- what do you mean by looking at the reads (do you mean viewing the bam file using IGV)

                  Comment

                  • dpryan
                    Devon Ryan
                    • Jul 2011
                    • 3478

                    #10
                    They're unmapped, so you would just "zcat unmapped.fq.gz | less" or even just run the unmapped reads through fastqc. The general idea is to first see if they're even worth assembling (if they're all repetitive regions or low quality then there's no reason to bother).

                    Comment

                    Latest Articles

                    Collapse

                    • seqadmin
                      Pathogen Surveillance with Advanced Genomic Tools
                      by seqadmin




                      The COVID-19 pandemic highlighted the need for proactive pathogen surveillance systems. As ongoing threats like avian influenza and newly emerging infections continue to pose risks, researchers are working to improve how quickly and accurately pathogens can be identified and tracked. In a recent SEQanswers webinar, two experts discussed how next-generation sequencing (NGS) and machine learning are shaping efforts to monitor viral variation and trace the origins of infectious...
                      03-24-2025, 11:48 AM
                    • seqadmin
                      New Genomics Tools and Methods Shared at AGBT 2025
                      by seqadmin


                      This year’s Advances in Genome Biology and Technology (AGBT) General Meeting commemorated the 25th anniversary of the event at its original venue on Marco Island, Florida. While this year’s event didn’t include high-profile musical performances, the industry announcements and cutting-edge research still drew the attention of leading scientists.

                      The Headliner
                      The biggest announcement was Roche stepping back into the sequencing platform market. In the years since...
                      03-03-2025, 01:39 PM

                    ad_right_rmr

                    Collapse

                    News

                    Collapse

                    Topics Statistics Last Post
                    Started by seqadmin, 03-20-2025, 05:03 AM
                    0 responses
                    49 views
                    0 reactions
                    Last Post seqadmin  
                    Started by seqadmin, 03-19-2025, 07:27 AM
                    0 responses
                    57 views
                    0 reactions
                    Last Post seqadmin  
                    Started by seqadmin, 03-18-2025, 12:50 PM
                    0 responses
                    50 views
                    0 reactions
                    Last Post seqadmin  
                    Started by seqadmin, 03-03-2025, 01:15 PM
                    0 responses
                    201 views
                    0 reactions
                    Last Post seqadmin  
                    Working...