text: user1770: 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? 
  
  user71: 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
  
  
  user29: 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!).
  
  user235: 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.
  
  user77: 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.
  
  user425: 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.
  
original_text: 1770: <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>
   71: <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>
   29: <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>
   235: <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>
   77: <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>
   425: <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>
  