text: user1770: I have a text file 'qnames.txt' with QNAMEs in the following format:
  EXAMPLE:QNAME1
  EXAMPLE:QNAME2
  EXAMPLE:QNAME3
  EXAMPLE:QNAME4
  EXAMPLE:QNAME5
  
  I would like to subset my BAM file.bam via all of these QNAMEs into a new SAM.
  Naturally, I can do this individually, e.g.
  samtools view file.bam | grep 'EXAMPLE:QNAME1' > subset.bam
  
  But I'm unsure how to do this for a list of QNAMES:
  
  How do I write a for loop that will do all of these queries, outputting the correct SAM needed?
  I could write a for loop that creates n SAM files, and then cats them…
  
  Is there a way to specifically grep by QNAME? The above may grep reads that may not be associated with the correct QNAME?
  
  How do I keep the BAM header?
  
  
  
  user425: You can use SAMsift:
  samsift -i file.bam -0 'q=open("qnames.txt").read().splitlines()' -f 'QNAME in q'
  
  -i file.bam specifies the input BAM file. -0 'q=open("qnames.txt").read().splitlines()' loads your list of QNAMEs. -f 'QNAME in q' specifies the filter for alignments.
  With SAMsift one can use any Python expression for filtering. The downside is that it's slower than the other tools.
  
  user6263: Most of the solutions posted thus far are either slow, or can yield spurious results in certain cases. (I haven't tested Picard, which I assume will work as intended. But personally I tend to recoil when I have to set up a JVM!) 
  Let's approach the problem from a different angle:
  Good hash functions have ~O(1) complexity for a lookup. If we can create a hash table of our QNAMEs, then searching a BAM file should be O(n) in the number of records. The Python developers are very proud of the hash function in cpython, so let's use that:
  #!/usr/bin/env python3
  
  import sys
  
  with open(sys.argv[1], 'r') as indexfile:
      ids = set(l.rstrip('\r\n') for l in indexfile)
  
  for line in sys.stdin:
      qname, _ = line.split('\t', 1)
      if qname in ids:
          sys.stdout.write(line)
  
  With a test set consisting of a ~3GB BAM file and a qnames.txt with ~65K entries, running samtools view input.bam | ./idfilter.py qnames.txt takes roughly 6.5 seconds on my test server. By comparison, piping the SAM to fgrep -f qnames.txt (using GNU grep) takes about 8.5 seconds (but may yield spurious results).
  (Why might grep yield spurious results? If you have a QNAME foo in your qnames.txt, but your BAM file also contains fooX and fooY, they will all be matched by fgrep. The solution is to anchor each of your search patterns between the beginning of the line and a tab, e.g. ^foo\t, but then you have to use standard grep which will have to construct an NFA for each pattern, and not fgrep which implements Aho-Corasick, and you'll be orders of magnitude slower.)
  I've rewritten my little Python script in Rust using std::collections::HashMap for the hash function, and using the rust_htslib crate to read and write straight from BAM, and this is more than twice as fast as the Python, even when reading and writing BAM. However, my Rust isn't really fit for public consumption… 
  
  user240: Performing a grep on $n$ alignments to $m$ qnames would give you $O(mn)$ operations. If you decide to sort your BAM and qnames first, you can reduce it to $O(\min\{m, n\})$ by popping reads and qnames off their stacks. You also don't read the entirety of the BAM at once this way, limiting memory consumption. I'm pretty sure Picard does something like this, as well.
  Here's an example in Python that I've written for this very purpose:
  import pysam
  
  def filter_qname(bamfile, idfile, outfile):
      '''
      Filter reads in a SAM/BAM file by their query names
  
      Parameters
      ----------
      bamfile : str
          Sorted BAM file path
      idfile : str
          Text file path containing qnames to keep
      outfile : str
          Output file to write to
      '''
      # read in name-sorted BAM file
      bam = pysam.AlignmentFile(bamfile, "rb")
      if outfile.endswith('.sam'):
          output = pysam.AlignmentFile(outputfile, "w", template=bam)
      elif outfile.endswith('.bam'):
          output = pysam.AlignmentFile(outputfile, "wb", template=bam)
      else:
          raise ValueError("Unknown output file format for `outfile`: {}".format(outfile))
      #  read in IDs to be removed (list(set(...)) to only keep unique IDs)
      ids = list(set([l.rstrip() for l in open(idfile, 'r').readlines()]))
      # sort IDs to match BAM for efficient processing (destructive function)
      ids.sort(key=str.lower, reverse=True)
      n_ids = len(ids)
      # variable for ensuring BAM is sorted by query name
      last_q = None
      reads = bam.fetch(until_eof=True)
      read = next(reads)
      while True:
          # if read name is greater than current top of stack
          if read.query_name > ids[-1]:
              # if this is the last ID
              if n_ids == 1:
                  # write remaining reads to new file and break out of while loop
                  output.write(read)
                  for r in bam:
                      output.write(read)
                  break
              # otherwise pop that ID, try again
              else:
                  ids.pop()
                  n_ids -= 1
          # skip reads that match the top of the stack
          elif read.query_name == ids[-1]:
              read = next(reads)
          # if read name is less that top of stack, write it and move on
          else:
              output.write(read)
              read = next(reads)
      output.close()
  
  I've implemented it in my own miscellaneous bioinformatics tools package here if you want to see it in action
  
  user964: (1)
  It won't be super fast but you can provide grep with a file of QNAMES.
  samtools view file.bam | grep -f 'qnames.txt > subset.sam
  
  where qnames.txt has 
  EXAMPLE:QNAME1
  EXAMPLE:QNAME2
  EXAMPLE:QNAME3
  EXAMPLE:QNAME4
  EXAMPLE:QNAME5
  
  (2)
  This would be a little more complicated but can you give an example where the grep might be have the correct QNAME?
  (3) 
  To keep the BAM header I would use samtools -H file.bam > header.txt to get the header and then cat the header file with the grep'ed file 
  
  user77: Below is a small bit of python code showing one naive but manageable way of doing this (N.B., I expect grep to be faster, though getting it to output the header will be annoying):
  import pysam
  qnames = set(...)  # read names go here
  bam = pysam.AlignmentFile("some input file.bam")
  obam = pysam.AlignmentFile("some output file.bam", "w", template=bam)
  
  for b in bam.fetch(until_eof=True):
      if b.query_name in qnames:
          obam.write(b)
  bam.close()
  obam.close()
  
  Now that'll work, but if you have a lot of names to look through then it'll end up being quite slow (as will grep -f ...). If you need to select a large list of reads by name, then a more efficient strategy is:
  
  Name sort the BAM file (samtools sort -n or picard), which can be done with multiple threads.
  Sort the list of read names to match (this will be more complicated than one would naively think, since samtools and picard will name sort differently, so be sure to put some thought into this).
  Perform the same sort of iteration as above, but only compare to the top read name in your list from step 2 (removing that element from the stack after a match).
  
  
  user2355: UPDATE 2021/06/28: since version 1.12, samtools now accepts option -N, which takes a file containing read names of interest.
  Using samtools 1.12 or greater:
  samtools view -N qnames_list.txt -o filtered_output.bam input.bam
  
  OLD ANSWER:
  When it comes to filter by a list, this is my favourite (much faster than grep):
  # given a bam file "aln.bam" and a list of read names "reads.txt":
  samtools view -h aln.bam \
  | awk 'FNR==NR {reads[$1];next} /^@/||($1 in reads)' reads.txt - \
  | samtools view -b - > filtered_aln.bam
  
  Basically, the awk command first parses the reads.txt file, creates an array of QNAMEs, then parses the alignments and print if the QNAME is present in the array (or if the line belongs to the header).
  
  user71: use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads
  
  READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null. 
  
  
original_text: 1770: <p>I have a text file 'qnames.txt' with QNAMEs in the following format:</p>
  <pre><code>EXAMPLE:QNAME1
  EXAMPLE:QNAME2
  EXAMPLE:QNAME3
  EXAMPLE:QNAME4
  EXAMPLE:QNAME5
  </code></pre>
  <p>I would like to subset my BAM <code>file.bam</code> via all of these QNAMEs into a new SAM.</p>
  <p>Naturally, I can do this individually, e.g.</p>
  <pre><code>samtools view file.bam | grep 'EXAMPLE:QNAME1' &gt; subset.bam
  </code></pre>
  <p>But I'm unsure how to do this for a list of QNAMES:</p>
  <ol>
  <li><p>How do I write a for loop that will do all of these queries, outputting the correct SAM needed?</p>
  <p>I could write a <code>for</code> loop that creates <em>n</em> SAM files, and then <code>cat</code>s them…</p>
  </li>
  <li><p>Is there a way to specifically <code>grep</code> by QNAME? The above may <code>grep</code> reads that may not be associated with the correct QNAME?</p>
  </li>
  <li><p>How do I keep the BAM header?</p>
  </li>
  </ol>
   425: <p>You can use <a href="https://github.com/karel-brinda/samsift" rel="nofollow noreferrer">SAMsift</a>:</p>
  
  <pre class="lang-bash prettyprint-override"><code>samsift -i file.bam -0 'q=open("qnames.txt").read().splitlines()' -f 'QNAME in q'
  </code></pre>
  
  <p><code>-i file.bam</code> specifies the input BAM file. <code>-0 'q=open("qnames.txt").read().splitlines()'</code> loads your list of QNAMEs. <code>-f 'QNAME in q'</code> specifies the filter for alignments.</p>
  
  <p>With SAMsift one can use any Python expression for filtering. The downside is that it's slower than the other tools.</p>
   6263: <p>Most of the solutions posted thus far are either slow, or can yield spurious results in certain cases. (I haven't tested Picard, which I assume will work as intended. But personally I tend to recoil when I have to set up a JVM!) </p>
  
  <p>Let's approach the problem from a different angle:</p>
  
  <p>Good hash functions have ~O(1) complexity for a lookup. If we can create a hash table of our QNAMEs, then searching a BAM file should be O(n) in the number of records. The Python developers are very proud of the hash function in cpython, so let's use that:</p>
  
  <pre><code>#!/usr/bin/env python3
  
  import sys
  
  with open(sys.argv[1], 'r') as indexfile:
      ids = set(l.rstrip('\r\n') for l in indexfile)
  
  for line in sys.stdin:
      qname, _ = line.split('\t', 1)
      if qname in ids:
          sys.stdout.write(line)
  </code></pre>
  
  <p>With a test set consisting of a ~3GB BAM file and a <code>qnames.txt</code> with ~65K entries, running <code>samtools view input.bam | ./idfilter.py qnames.txt</code> takes roughly <strong>6.5 seconds</strong> on my test server. By comparison, piping the SAM to <code>fgrep -f qnames.txt</code> (using GNU grep) takes about 8.5 seconds (but may yield spurious results).</p>
  
  <p>(Why might <code>grep</code> yield spurious results? If you have a QNAME <code>foo</code> in your <code>qnames.txt</code>, but your BAM file also contains <code>fooX</code> and <code>fooY</code>, they will all be matched by <code>fgrep</code>. The solution is to anchor each of your search patterns between the beginning of the line and a tab, e.g. <code>^foo\t</code>, but then you have to use standard <code>grep</code> which will have to construct an NFA for each pattern, and not <code>fgrep</code> which implements Aho-Corasick, and you'll be orders of magnitude slower.)</p>
  
  <p>I've rewritten my little Python script in Rust using <code>std::collections::HashMap</code> for the hash function, and using the rust_htslib crate to read and write straight from BAM, and this is more than twice as fast as the Python, even when reading and writing BAM. However, my Rust isn't really fit for public consumption… </p>
   240: <p>Performing a grep on <span class="math-container">$n$</span> alignments to <span class="math-container">$m$</span> qnames would give you <span class="math-container">$O(mn)$</span> operations. If you decide to sort your BAM and qnames first, you can reduce it to <span class="math-container">$O(\min\{m, n\})$</span> by popping reads and qnames off their stacks. You also don't read the entirety of the BAM at once this way, limiting memory consumption. I'm pretty sure Picard does something like this, as well.</p>
  
  <p>Here's an example in Python that I've written for this very purpose:</p>
  
  <pre><code>import pysam
  
  def filter_qname(bamfile, idfile, outfile):
      '''
      Filter reads in a SAM/BAM file by their query names
  
      Parameters
      ----------
      bamfile : str
          Sorted BAM file path
      idfile : str
          Text file path containing qnames to keep
      outfile : str
          Output file to write to
      '''
      # read in name-sorted BAM file
      bam = pysam.AlignmentFile(bamfile, "rb")
      if outfile.endswith('.sam'):
          output = pysam.AlignmentFile(outputfile, "w", template=bam)
      elif outfile.endswith('.bam'):
          output = pysam.AlignmentFile(outputfile, "wb", template=bam)
      else:
          raise ValueError("Unknown output file format for `outfile`: {}".format(outfile))
      #  read in IDs to be removed (list(set(...)) to only keep unique IDs)
      ids = list(set([l.rstrip() for l in open(idfile, 'r').readlines()]))
      # sort IDs to match BAM for efficient processing (destructive function)
      ids.sort(key=str.lower, reverse=True)
      n_ids = len(ids)
      # variable for ensuring BAM is sorted by query name
      last_q = None
      reads = bam.fetch(until_eof=True)
      read = next(reads)
      while True:
          # if read name is greater than current top of stack
          if read.query_name &gt; ids[-1]:
              # if this is the last ID
              if n_ids == 1:
                  # write remaining reads to new file and break out of while loop
                  output.write(read)
                  for r in bam:
                      output.write(read)
                  break
              # otherwise pop that ID, try again
              else:
                  ids.pop()
                  n_ids -= 1
          # skip reads that match the top of the stack
          elif read.query_name == ids[-1]:
              read = next(reads)
          # if read name is less that top of stack, write it and move on
          else:
              output.write(read)
              read = next(reads)
      output.close()
  </code></pre>
  
  <p>I've implemented it in my own miscellaneous bioinformatics tools package <a href="https://github.com/jrhawley/bio-jtools/blob/master/jtools/align/filter_qname.py#L23-L72" rel="nofollow noreferrer">here</a> if you want to see it in action</p>
   964: <p>(1)
  It won't be super fast but you can provide grep with a file of QNAMES.</p>
  
  <pre><code>samtools view file.bam | grep -f 'qnames.txt &gt; subset.sam
  </code></pre>
  
  <p>where qnames.txt has </p>
  
  <pre><code>EXAMPLE:QNAME1
  EXAMPLE:QNAME2
  EXAMPLE:QNAME3
  EXAMPLE:QNAME4
  EXAMPLE:QNAME5
  </code></pre>
  
  <p>(2)
  This would be a little more complicated but can you give an example where the grep might be have the correct QNAME?</p>
  
  <p>(3) 
  To keep the BAM header I would use <code>samtools -H file.bam &gt; header.txt</code> to get the header and then cat the header file with the grep'ed file </p>
   77: <p>Below is a small bit of python code showing one naive but manageable way of doing this (N.B., I expect <code>grep</code> to be faster, though getting it to output the header will be annoying):</p>
  
  <pre class="lang-py prettyprint-override"><code>import pysam
  qnames = set(...)  # read names go here
  bam = pysam.AlignmentFile("some input file.bam")
  obam = pysam.AlignmentFile("some output file.bam", "w", template=bam)
  
  for b in bam.fetch(until_eof=True):
      if b.query_name in qnames:
          obam.write(b)
  bam.close()
  obam.close()
  </code></pre>
  
  <p>Now that'll work, but if you have a lot of names to look through then it'll end up being quite slow (as will <code>grep -f ...</code>). If you need to select a large list of reads by name, then a more efficient strategy is:</p>
  
  <ol>
  <li>Name sort the BAM file (<code>samtools sort -n</code> or picard), which can be done with multiple threads.</li>
  <li>Sort the list of read names to match (this will be more complicated than one would naively think, since samtools and picard will name sort differently, so be sure to put some thought into this).</li>
  <li>Perform the same sort of iteration as above, but only compare to the top read name in your list from step 2 (removing that element from the stack after a match).</li>
  </ol>
   2355: <p><strong>UPDATE 2021/06/28:</strong> since version 1.12, samtools now accepts option <code>-N</code>, which takes a file containing read names of interest.</p>
  <p>Using samtools 1.12 or greater:</p>
  <pre><code>samtools view -N qnames_list.txt -o filtered_output.bam input.bam
  </code></pre>
  <p><strong>OLD ANSWER:</strong></p>
  <p>When it comes to filter by a list, this is my favourite (much faster than grep):</p>
  <pre><code># given a bam file &quot;aln.bam&quot; and a list of read names &quot;reads.txt&quot;:
  samtools view -h aln.bam \
  | awk 'FNR==NR {reads[$1];next} /^@/||($1 in reads)' reads.txt - \
  | samtools view -b - &gt; filtered_aln.bam
  </code></pre>
  <p>Basically, the awk command first parses the reads.txt file, creates an array of QNAMEs, then parses the alignments and print if the QNAME is present in the array (or if the line belongs to the header).</p>
   71: <p>use picard FilterSamReads <a href="http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads" rel="noreferrer">http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads</a></p>
  
  <blockquote>
    <p>READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null. </p>
  </blockquote>
  