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

CC BY-SA 3.0
2018-01-24T22:05:01.120Z
2018-01-24T22:05:01.120Z
Pierre
<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>

CC BY-SA 4.0
2020-10-14T18:17:40.730Z
2021-10-22T15:06:13.123Z
2021-10-22T15:06:13.123Z
bepoli
<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>

CC BY-SA 3.0
2018-01-24T15:12:59.440Z
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`
Konrad Rudolph
CC BY-SA 3.0
2018-01-23T20:21:35.443Z
2018-01-23T20:21:35.443Z
Devon Ryan
<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, &quot;rb&quot;)
    if outfile.endswith('.sam'):
        output = pysam.AlignmentFile(outfile, &quot;w&quot;, template=bam)
    elif outfile.endswith('.bam'):
        output = pysam.AlignmentFile(outfile, &quot;wb&quot;, template=bam)
    else:
        raise ValueError(&quot;Unknown output file format for `outfile`: {}&quot;.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>

CC BY-SA 4.0
2019-09-27T13:32:04.073Z
2024-10-31T18:21:35.510Z
2024-10-31T18:21:35.510Z
James Hawley
<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>

CC BY-SA 4.0
2019-10-22T20:27:24.280Z
This is a great answer! It's quite comprehensive, and I agree with your distaste of setting up a JVM :)
EB2127
CC BY-SA 4.0
2019-09-25T08:28:47.760Z
2019-09-25T08:28:47.760Z
wjv
<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>

CC BY-SA 3.0
2018-02-08T21:30:12.453Z
2018-02-08T21:58:41.100Z
2018-02-08T21:58:41.100Z
Karel Břinda
<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>

CC BY-SA 4.0
2018-12-11T22:10:10.173Z
This is quite slow if the bam file or the list of reads to keep are of appreciable size.
conchoecia
CC BY-SA 4.0
2019-09-25T07:52:57.610Z
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
CC BY-SA 4.0
2022-10-24T09:28:01.423Z
As it was already said, this is very slow. You can easily do this much faster with samtools or picard (see mine or Pierre's answers).
bepoli
CC BY-SA 3.0
2018-01-23T18:49:22.583Z
2018-01-23T18:49:22.583Z
Bioathlete
<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>

CC BY-SA 3.0
2018-01-24T15:15:45.740Z
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.
Konrad Rudolph
CC BY-SA 4.0
2019-04-23T15:58:07.610Z
@KonradRudolph mind adding that solution over here?
Beeba
CC BY-SA 4.0
2019-04-23T19:58:25.403Z
@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.
Konrad Rudolph
CC BY-SA 4.0
2019-04-23T20:17:45.953Z
@KonradRudolph makes sense, no worries. I will look into samtools split. Thanks
Beeba
<sam><samtools>
How to subset a BAM by a list of QNAMEs?
# How to subset a BAM by a list of QNAMEs?

(Asked by: username_0 on 2018-01-23)

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?

**Tags:** sam, samtools
> **Comments:**
> **username_1**: 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.  
> **username_2**: @KonradRudolph mind adding that solution over here?  
> **username_1**: @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.  
> **username_2**: @KonradRudolph makes sense, no worries. I will look into samtools split. Thanks  

---

## 7 answers

### Answer by username_9 on 2018-01-23. Score: 4 | Accepted answer

(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

> **Comments:**
> **username_10**: This is quite slow if the bam file or the list of reads to keep are of appreciable size.  
> **username_7**: 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.  
> **username_4**: As it was already said, this is very slow. You can easily do this much faster with samtools or picard (see mine or Pierre's answers).  

---

### Answer by username_3 on 2018-01-24. Score: 9

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.


---

### Answer by username_4 on 2020-10-14. Score: 8

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


---

### Answer by username_5 on 2018-01-23. Score: 4

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

> **Comments:**
> **username_1**: 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`  

---

### Answer by username_6 on 2019-09-27. Score: 3

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(outfile, "w", template=bam)
    elif outfile.endswith('.bam'):
        output = pysam.AlignmentFile(outfile, "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


---

### Answer by username_7 on 2019-09-25. Score: 2

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…

> **Comments:**
> **username_0**: This is a great answer! It's quite comprehensive, and I agree with your distaste of setting up a JVM :)  

---

### Answer by username_8 on 2018-02-08. Score: 1

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.


---

