5
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Konrad Rudolph Jan 24 '18 at 15:15
  • $\begingroup$ @KonradRudolph mind adding that solution over here? $\endgroup$ – Beeba Apr 23 at 15:58
  • $\begingroup$ @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. $\endgroup$ – Konrad Rudolph Apr 23 at 19:58
  • 1
    $\begingroup$ @KonradRudolph makes sense, no worries. I will look into samtools split. Thanks $\endgroup$ – Beeba Apr 23 at 20:17
3
$\begingroup$

(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

$\endgroup$
  • 1
    $\begingroup$ This is quite slow if the bam file or the list of reads to keep are of appreciable size. $\endgroup$ – conchoecia Dec 11 '18 at 22:10
  • 1
    $\begingroup$ 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. $\endgroup$ – wjv Sep 25 at 7:52
9
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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).
$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – Konrad Rudolph Jan 24 '18 at 15:12
2
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.