# How to subset a BAM by a list of QNAMEs?

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:

1. 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…

2. Is there a way to specifically grep by QNAME? The above may grep reads that may not be associated with the correct QNAME?

3. How do I keep the BAM header?

• Are these QNAMEs completely unrelated? If, by chance, you simply want to select all reads that belong to the same read group (which is determined by the QNAME), then there’s a much more efficient and simpler solution. Jan 24 '18 at 15:15
• @KonradRudolph mind adding that solution over here? Apr 23 '19 at 15:58
• @Beeba Given that OP hasn’t confirmed whether splitting by read groups is indeed desired, I’m reluctant to post an unrelated answer here. That said, it’s simply what the command samtools split does. Apr 23 '19 at 19:58
• @KonradRudolph makes sense, no worries. I will look into samtools split. Thanks Apr 23 '19 at 20:17

(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

• This is quite slow if the bam file or the list of reads to keep are of appreciable size. Dec 11 '18 at 22:10
• Since the qnames are all fixed strings, using fgrep instead of grep will provide a significant speedup — several orders of magnitude. However, both approaches may lead to spurious results: If you have a qname foo in our list, and the BAM file also contains qnames fooX and fooY, the latter will also be returned in your results, erroneously.
– wjv
Sep 25 '19 at 7:52

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.

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:

1. Name sort the BAM file (samtools sort -n or picard), which can be done with multiple threads.
2. 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).
3. 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).
• Keeping the BAM header with the grep approach is simply a matter of echoing it first: (samtools -H input.bam; samtools input.bam | grep …) | samtools -b - -o output.bam Jan 24 '18 at 15:12

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


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).

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… 😬

• This is a great answer! It's quite comprehensive, and I agree with your distaste of setting up a JVM :) Oct 22 '19 at 20:27

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
while True:
# if read name is greater than current top of stack
# if this is the last ID
if n_ids == 1:
# write remaining reads to new file and break out of while loop
for r in bam:
break
# otherwise pop that ID, try again
else:
ids.pop()
n_ids -= 1
# skip reads that match the top of the stack
# if read name is less that top of stack, write it and move on
else:

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.