<p>You can use <a href="https://github.com/karel-brinda/samsift" rel="nofollow noreferrer">SAMsift</a>:</p>

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

<p>Explanation:</p>

<ul>
<li><code>-i file.bam</code> – input file</li>
<li><code>-0 'c={"chr1":5000,"chr2":5000}'</code> – initialization (create a count-down dictionary for the chromosomes of interest)</li>
<li><code>-f 'c[RNAME]&gt;0'</code> – filtering criterion (is the counter for the current chromosome still >0?)</li>
<li><code>-c 'c[RNAME]-=1'</code> – code decrementing the counter of the current chromosome (5000 -> 4999 -> ... -> 1 -> 0 -> -1 -> ...)</li>
<li><code>-m nonstop-remove</code> – 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)</li>
</ul>

<p>See the <a href="https://github.com/karel-brinda/samsift" rel="nofollow noreferrer">SAMsift readme</a> for more info.</p>

CC BY-SA 3.0
2018-02-26T23:08:03.597Z
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
llrs
CC BY-SA 3.0
2018-02-26T23:18:01.753Z
I've just added the explanation.
Karel Břinda
CC BY-SA 3.0
2018-02-26T22:30:14.907Z
2018-02-26T23:22:00.417Z
2018-02-26T23:22:00.417Z
Karel Břinda
<p>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):</p>



<pre class="lang-python prettyprint-override"><code>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() &lt; frac:
            count += 1
            outfile.write(read)
     print("ouputted %i reads from %s" %(count, chr))
outfile.close()
</code></pre>

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

<p>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:</p>

<pre class="lang-python prettyprint-override"><code>for read in infile.fetch(chr, multiple_iterators=True):
    count += 1
    if count &gt;= nreads:
        break
    outfile.write(read)
</code></pre>

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

<pre class="lang-python prettyprint-override"><code>mate_read = infile.mate(read)
outfile.write(mate_read)
</code></pre>

<p>Note that this is slow.</p>

CC BY-SA 3.0
2018-02-17T15:34:15.497Z
To speed things up, `infile.count(chr)` can be replaced with parsing `idxstats()` (unless one is using CRAM files).
Devon Ryan
CC BY-SA 3.0
2018-02-17T15:54:03.713Z
Thanks, I'd assumed that there would be an easy way to fetch this from the index.
Ian Sudbery
CC BY-SA 3.0
2018-02-19T21:00:28.350Z
"To speed things up, infile.count(chr) can be replaced with parsing idxstats() " I'm not quite following. What would the above look like?
EB2127
CC BY-SA 3.0
2018-02-17T15:20:09.120Z
2018-02-17T15:20:09.120Z
Ian Sudbery
<p>If you just want to create a truncated BAM file with a header then you can significantly simplify your original code:</p>

<pre class="lang-bash prettyprint-override"><code>(samtools view -H input.bam; samtools view input.bam | head -5000) | samtools -bo output.bam
</code></pre>

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

<p>As above, but with formatting:</p>

<pre class="lang-bash prettyprint-override"><code>(
    samtools view -H input.bam
    samtools view input.bam | head -5000 # (*)
) \
| samtools -bo output.bam
</code></pre>

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

CC BY-SA 3.0
2018-02-26T23:29:11.247Z
Is it `samtools -bo`?
EB2127
CC BY-SA 3.0
2018-02-20T15:08:42.420Z
2018-02-20T15:24:00.137Z
2018-02-20T15:24:00.137Z
Konrad Rudolph
<p>using <strong>samjdk</strong>: <a href="http://lindenb.github.io/jvarkit/SamJdk.html" rel="nofollow noreferrer">http://lindenb.github.io/jvarkit/SamJdk.html</a></p>

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

CC BY-SA 3.0
2018-02-19T14:28:40.233Z
It'd be nice if people down-voting answers would explain why they're doing so...
Devon Ryan
CC BY-SA 3.0
2018-02-26T23:02:54.243Z
@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
llrs
CC BY-SA 3.0
2018-02-19T13:49:23.270Z
2018-02-19T13:49:23.270Z
Pierre
<p>If you're not too hung up on exact numbers like 5000 reads then you can do that with a single samtools command:</p>

<pre><code>samtools view -bo subset.bam -s 123.4 alignments.bam chr1 chr2
</code></pre>

<p>That will select 40% (the <code>.4</code> part) of the reads (<code>123</code> 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 <code>.4</code> part to a sufficiently small number.</p>

<p>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.</p>

CC BY-SA 3.0
2018-02-17T12:29:12.107Z
2018-02-17T12:29:12.107Z
Devon Ryan
<p>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.)</p>

<p>Normally I do the following, which will subset the BAM <code>file.bam</code> and keep the header</p>

<pre><code>samtools view -H file.bam &gt; header.sam
samtools view file.bam | head -n 5000 | cat header.sam - | samtools view -Sb - &gt; file.unique.bam
</code></pre>

<p>In this case, I would like 5000 rows in chromosome 1 and 5000 in chromosome 2. </p>

<p>I could first try grepping by individual chromosome, and then combining the two SAMs</p>

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

<pre><code>samtools view -H file.bam &gt; header.sam
samtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - &gt; file.unique.bam
</code></pre>

<p>but then I have two problems: </p>

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

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

<p>Is there an easy way, Bioinformatics SO? </p>

<sam><samtools>
Subset smaller BAM to contain several thousand rows from multiple chromosomes
# Subset smaller BAM to contain several thousand rows from multiple chromosomes

(Asked by: username_0 on 2018-02-17)

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?

**Tags:** sam, samtools

---

## 5 answers

### Answer by username_4 on 2018-02-17. Score: 3 | Accepted answer

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 by username_1 on 2018-02-26. Score: 3

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.

> **Comments:**
> **username_2**: 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  
> **username_1**: I've just added the explanation.  

---

### Answer by username_3 on 2018-02-17. Score: 2

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.

> **Comments:**
> **username_4**: To speed things up, `infile.count(chr)` can be replaced with parsing `idxstats()` (unless one is using CRAM files).  
> **username_3**: Thanks, I'd assumed that there would be an easy way to fetch this from the index.  
> **username_0**: "To speed things up, infile.count(chr) can be replaced with parsing idxstats() " I'm not quite following. What would the above look like?  

---

### Answer by username_5 on 2018-02-20. Score: 1

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

> **Comments:**
> **username_0**: Is it `samtools -bo`?  

---

### Answer by username_6 on 2018-02-19. Score: 1

using samjdk: http://lindenb.github.io/jvarkit/SamJdk.html
$ java -jar dist/samjdk.jar --body -e \
  'Map 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

> **Comments:**
> **username_4**: It'd be nice if people down-voting answers would explain why they're doing so...  
> **username_2**: @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  

---

