Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-26T23:08:03.597Z
Answers.Text: 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
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-26T23:18:01.753Z
Answers.Text: I've just added the explanation.
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-26T22:30:14.907Z
Answers.LastActivityDate: 2018-02-26T23:22:00.417Z
Answers.LastEditDate: 2018-02-26T23:22:00.417Z
Answers.UserDisplayName: Karel Břinda
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-17T15:34:15.497Z
Answers.Text: To speed things up, `infile.count(chr)` can be replaced with parsing `idxstats()` (unless one is using CRAM files).
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-17T15:54:03.713Z
Answers.Text: Thanks, I'd assumed that there would be an easy way to fetch this from the index.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-19T21:00:28.350Z
Answers.Text: "To speed things up, infile.count(chr) can be replaced with parsing idxstats() " I'm not quite following. What would the above look like?
Answers.UserDisplayName: EB2127
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-17T15:20:09.120Z
Answers.LastActivityDate: 2018-02-17T15:20:09.120Z
Answers.UserDisplayName: Ian Sudbery
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-26T23:29:11.247Z
Answers.Text: Is it `samtools -bo`?
Answers.UserDisplayName: EB2127
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-20T15:08:42.420Z
Answers.LastActivityDate: 2018-02-20T15:24:00.137Z
Answers.LastEditDate: 2018-02-20T15:24:00.137Z
Answers.UserDisplayName: Konrad Rudolph
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-19T14:28:40.233Z
Answers.Text: It'd be nice if people down-voting answers would explain why they're doing so...
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-26T23:02:54.243Z
Answers.Text: @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
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-19T13:49:23.270Z
Answers.LastActivityDate: 2018-02-19T13:49:23.270Z
Answers.UserDisplayName: Pierre
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-02-17T12:29:12.107Z
Answers.LastActivityDate: 2018-02-17T12:29:12.107Z
Answers.UserDisplayName: Devon Ryan
Body: <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>
  
Tags: <sam><samtools>
Title: Subset smaller BAM to contain several thousand rows from multiple chromosomes
ThreadText: # 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](https://github.com/username_1-brinda/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](https://github.com/username_1-brinda/samsift) 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 username_4 username_4'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):
  
  ```python
  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:
  
  ```python
  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)`:
  
  ```python
  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:
  
  ```bash
  (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:
  
  ```bash
  (
      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<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
  
  ```
  
  > **Comments:**
  > **username_4**: It'd be nice if people down-voting answers would explain why they're doing so...  
  > **username_2**: @username_4 (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  
  
  ---
  
  