Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • samtools filter PE reads on mapping quality of only one of the reads

    Hi,

    I'm trying to filter a BAM file using samtools for mapping quality, but I want to keep both reads in a pair if at least one of them has MAPQ>=10.

    My understanding is that using
    Code:
    -q 10
    as an option would remove all reads with MAPQ<10 regardless of whether the other read in the pair is of high quality. Is there a way to get around that?

    Thanks

  • #2
    There's no way to do that directly with samtools, you'll need to script something. It'd be pretty straight forward with pysam, if you know a bit of python (I'm in the middle of a bunch of C code at the moment or I'd through together a quick script for you).

    Comment


    • #3
      ahh, I was afraid that would be the case...Thanks for your reply anyway!
      I don't know python so I'm going to have to do it painfully in R.

      Comment


      • #4
        Too bad, I recall that the R method of dealing with SAM/BAM files is to read the whole thing into memory, which would really suck. Post the script you use when you have it working in case someone else ends up needing this in the future.

        Comment


        • #5
          Yep, that's exactly why it's going to be painful...Sure, I'll post it if/when it works!

          Comment


          • #6
            I ended up getting side-tracked by a dissection anyway, so here's a quick method using python and the pysam package. It assumes that you have everything name-sorted and that there aren't any random single-end reads in the file:

            Code:
            #!/usr/bin/env python
            import pysam
            import argparse
            import sys
            
            parser = argparse.ArgumentParser(description='Filter a paired-end BAM file such that both reads have MAPQs >= a threshold.')
            parser.add_argument('threshold', type=int, help="MAPQ threshold")
            parser.add_argument('input_file', help="input file name")
            parser.add_argument('output_file', help="output file name")
            args = parser.parse_args()
            if(args.threshold == None or args.input_file == None or args.output_file == None) :
                parser.print_help()
                sys.exit()
            
            #One could look at the file name too see if we're dealing with a SAM or BAM file, but this is simpler
            ifile = pysam.Samfile(args.input_file, "rb")
            ofile = pysam.Samfile(args.output_file, "wb", template=ifile)
            
            #iterate
            while(1) :
                try :
                    r1 = ifile.next()
                    r2 = ifile.next()
                except :
                    break
            
                if(r1.mapq >= int(args.threshold) and r2.mapq >= int(args.threshold)) :
                    ofile.write(r1)
                    ofile.write(r2)
            
            ifile.close()
            ofile.close()

            Comment


            • #7
              Thanks a lot for doing this.
              Quick question: what I'd like is to filter the file so that AT LEAST ONE of the reads has MAPQs>= threshold (as opposed to BOTH);

              For that should I just change this section of your code:
              Code:
              if(r1.mapq >= int(args.threshold) and r2.mapq >= int(args.threshold)) :
              to this:?
              Code:
              if(r1.mapq >= int(args.threshold) or r2.mapq >= int(args.threshold)) :
              or would that not keep pairs where both are above the threshold?
              thanks again!

              Comment


              • #8
                That should work. "or" isn't an "exclusive or" in python.

                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
                67 views
                0 likes
                Last Post seqadmin  
                Working...
                X