Subset smaller BAM to contain several thousand rows from multiple chromosomes

Question

There are many cases whereby I would like to subset a BAM to create a small file in order to work with (e.g. algorithmic testing, debugging, etc.)

Normally I do the following, which will subset the BAM file.bam and keep the header

samtools view -H file.bam > header.sam
samtools view file.bam | head -n 5000 | cat header.sam - | samtools view -Sb - > file.unique.bam

In this case, I would like 5000 rows in chromosome 1 and 5000 in chromosome 2.

I could first try grepping by individual chromosome, and then combining the two SAMs

e.g. here's complete BAM with grepped chr1 and (incorrect but complete) header

samtools view -H file.bam > header.sam
samtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - > file.unique.bam

but then I have two problems:

(1) I may not be grepping the alignments to chromsome 2---there may exist BAM rows which contain 'chr2' but are not alignments.

(2) I think one must manually edit the header. There's probably no way around this.

Is there an easy way, Bioinformatics SO?

Answer 1

If you're not too hung up on exact numbers like 5000 reads then you can do that with a single samtools command:

samtools view -bo subset.bam -s 123.4 alignments.bam chr1 chr2

That will select 40% (the .4 part) of the reads (123 is a seed, which is convenient for reproducibility). The convenient part of this is that it'll keep mates paired if you have paired-end reads. For 5000 reads per chromosome just change the .4 part to a sufficiently small number.

In general you don't really need to subset the header. Some tools will perform a bit better if you do, but you'll generally get the same results regardless.

Answer 2

You can use SAMsift:

samsift \
   -i file.bam \
   -0 'c={"chr1":5000,"chr2":5000}' \
   -f 'c[RNAME]>0' \
   -c 'c[RNAME]-=1' \
   -m nonstop-remove

Explanation:

• -i file.bam – input file
• -0 'c={"chr1":5000,"chr2":5000}' – initialization (create a count-down dictionary for the chromosomes of interest)
• -f 'c[RNAME]>0' – filtering criterion (is the counter for the current chromosome still >0?)
• -c 'c[RNAME]-=1' – code decrementing the counter of the current chromosome (5000 -> 4999 -> ... -> 1 -> 0 -> -1 -> ...)
• -m nonstop-remove – remove lines causing Python errors and don't stop (in this case, an error can be caused by accessing a non-existing counter for another chromosome, e.g., for chr3)

See the SAMsift readme for more info.

Answer 3

I'd generally recommend Devon Ryan's answer. However, if you do care about having the same number from each Chr, you could use the following python/pysam code (this will output approx 5000 from every Chr):

from pysam import AlignmentFile
from random import random

nreads = 5000

infile = AlignmentFile("mybam.bam")
outfile = AlignmentFile("outbam.bam", "wb", template=infile)

for chr in infile.reference_names:
     reads_in_chr = infile.count(chr)
     frac = float(nreads)/reads_in_chr
     count = 0
     #-- Replace this loop for the exactly 5000 first reads --#
     for read in infile.fetch(chr, multiple_iterators=True):
        if random() < frac:
            count += 1
            outfile.write(read)
     print("ouputted %i reads from %s" %(count, chr))
outfile.close()

This uses all chromosomes. If you only want to use, say chr1 and chr2 replace infile.reference_names with ["chr1","chr2"].

If you want exactly 5000 reads from each chr, but don't care if they are the first ones or not, then you could replace the inner for loop with:

for read in infile.fetch(chr, multiple_iterators=True):
    count += 1
    if count >= nreads:
        break
    outfile.write(read)

if you wanted the mates of those reads as well, you could add the following after outfile.write(read):

mate_read = infile.mate(read)
outfile.write(mate_read)

Note that this is slow.

Answer 4

using samjdk: http://lindenb.github.io/jvarkit/SamJdk.html

$ java -jar dist/samjdk.jar --body -e \
  'Map<String,Integer> c=new HashMap<>(); public Object apply(SAMRecord r) {int n=c.getOrDefault(r.getContig(),0);if(n>=5000) return false; c.put(r.getContig(),n+1); return r;}' \
 input.bam

Answer 5

If you just want to create a truncated BAM file with a header then you can significantly simplify your original code:

(samtools view -H input.bam; samtools view input.bam | head -5000) | samtools -bo output.bam

This command avoids the intermediate file and few gratuitous intermediate command invocations, at the cost of a subshell (invoked by (…)).

As above, but with formatting:

(
    samtools view -H input.bam
    samtools view input.bam | head -5000 # (*)
) \
| samtools -bo output.bam

… this can of course be extended to filter by multiple chromosomes by replacing the line marked with (*) above by one or multiple lines that subset by chromosome name (samtools view input.bam chrx, no need for grep if you have indexed the original BAM file!).