6
$\begingroup$

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?

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ If you could add some explanation about what are those options and why they are necessary. For instance I don't understand why the argument to f is c[RNAME] and not other things $\endgroup$ – llrs Feb 26 '18 at 23:08
  • 1
    $\begingroup$ I've just added the explanation. $\endgroup$ – Karel Brinda Feb 26 '18 at 23:18
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ To speed things up, infile.count(chr) can be replaced with parsing idxstats() (unless one is using CRAM files). $\endgroup$ – Devon Ryan Feb 17 '18 at 15:34
  • $\begingroup$ Thanks, I'd assumed that there would be an easy way to fetch this from the index. $\endgroup$ – Ian Sudbery Feb 17 '18 at 15:54
  • $\begingroup$ "To speed things up, infile.count(chr) can be replaced with parsing idxstats() " I'm not quite following. What would the above look like? $\endgroup$ – EB2127 Feb 19 '18 at 21:00
1
$\begingroup$

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
$\endgroup$
  • 3
    $\begingroup$ It'd be nice if people down-voting answers would explain why they're doing so... $\endgroup$ – Devon Ryan Feb 19 '18 at 14:28
  • 1
    $\begingroup$ @Devon (I have seen the comment until now) I downvoted because I find that a long line of code without explaining what it does internally is usually harder to understand, the code as is would improve much (IMHO) broken in several lines and with comments. It is hard for me that I don't know java to understand what is happening $\endgroup$ – llrs Feb 26 '18 at 23:02
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Is it samtools -bo? $\endgroup$ – EB2127 Feb 26 '18 at 23:29

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.