Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Removing remaining read in pair after samtools mapq filtering

    Hi,

    I am filtering a bam file using samtools (and picard for duplicates) and am struggling to understand the results from flagstat (apologies in advance for the long post)

    This is my sorted bam file before any filtering:

    Code:
    samtools flagstat Sample.sorted.bam
    86744514 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 duplicates
    16772228 + 0 mapped (19.34%:nan%)
    86744514 + 0 paired in sequencing
    43372257 + 0 read1
    43372257 + 0 read2
    14664488 + 0 properly paired (16.91%:nan%)
    15934998 + 0 with itself and mate mapped
    837230 + 0 singletons (0.97%:nan%)
    833002 + 0 with mate mapped to a different chr
    317794 + 0 with mate mapped to a different chr (mapQ>=5)
    Now, if I filter with -f2 to keep only those reads mapping in proper pair, flagstsat gives this:

    Code:
    samtools flagstat Sample.sorted.mapped.bam
    14664488 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 duplicates
    14664488 + 0 mapped (100.00%:nan%)
    14664488 + 0 paired in sequencing
    7332244 + 0 read1
    7332244 + 0 read2
    14664488 + 0 properly paired (100.00%:nan%)
    14664488 + 0 with itself and mate mapped
    0 + 0 singletons (0.00%:nan%)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)
    Notice how I get the same number of read 1 and read 2, which makes sense to me.

    But if I then run -q10 to remove any reads that have mapq<10, I get this:

    Code:
    samtools flagstat Sample.sorted.mapped_qualfiltered.bam
    10803413 + 0 in total (QC-passed reads + QC-failed reads)
    0 + 0 duplicates
    10803413 + 0 mapped (100.00%:nan%)
    10803413 + 0 paired in sequencing
    5398606 + 0 read1
    5404807 + 0 read2
    10803413 + 0 properly paired (100.00%:nan%)
    10803413 + 0 with itself and mate mapped
    0 + 0 singletons (0.00%:nan%)
    0 + 0 with mate mapped to a different chr
    0 + 0 with mate mapped to a different chr (mapQ>=5)
    (btw, I get the exact same result if doing the flag and mapq filters in the same command)

    If I'm not doing something incredibly stupid here, my first question is this: I'm getting a different number of read1 and read2 so I'm assuming samtools is removing only the reads that fail the quality filter and not their mate? Is there a way to tell it to remove those? i can't use flags at this stage as they are still labelled as properly paired even though the mate has been removed.

  • #2
    What aligner are you using? My only guess is that it gives mates different MAPQ scores, so the filtering can then create orphans. I'm not aware of an aligner that does that, of course...

    Comment


    • #3
      I'm using bwa, so really don't understand what happened...

      Comment


      • #4
        bwa mem or bwa aln? The mem version can produce chimeric alignments, which might further complicate things. Aside from that, I don't know why that's happening, you might just have to peek at the data.

        Comment


        • #5
          bwa aln.
          Yep I'll extract those reads and take a look. In fact, is there a way using samtools to exclude reads that have eg: mapq>10 (in the same way that you can use -F to exclude reads with certain flags?), or will I have to script that?
          Thanks again for taking the time to help me with this.

          Comment


          • #6
            No prob, it's a nice break from my own stuff.
            I'm not aware of a samtools flag that would do that, but something like:

            Code:
            samtools view -h foo.bam | awk '{if(substr($1,0,0) == "@") { print $0; } else if($5< 10) { print $0 }}' | samtools view -Sbo filtered.bam -
            should do the job for you.

            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