text: user4141: I have a bam file with 1 billion alignment reads of which there are 700 million unique reads. I want to split the alignments into chunks for parallel-processing. Multi-alignments of the same read should be in the same chunk. I want to use pysam to solve the problem and would be very grateful for suggestions.
  
  user1310: python will be to slow for this job. Here's a awk solution. One need to sort by read name and take track over the number of reads per chunk. If the number is reached the next chunk can only created of the read name is different.
  $ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
    BEGIN { part=0; line=n }       
    /^@/ {header = header$0"\n"; next;} 
    { if( line>=n && $1!=last_read ) {part++; line=1; printf header > part".sam"; print $0 >> part".sam" } 
      else { print $0 >> part".sam"; line++; } 
      last_read = $1;
    }'
  
  The output will be in sam format.
  You can controll the number of reads per chunk by setting the n value for the awk command.
  EDIT:
  A second version that will output in bam:
  $ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" ' 
  BEGIN { part=0; line=n }  
  /^@/ {header = header$0"\n"; next;}  
  { if( line>=n && $1!=last_read ) {print part,line; part++; line=1;}  
    print line==1 ? header""$0 : $0 | "samtools view -b -o "part".bam" 
    last_read = $1;
    line++;
  }
  
  
  user235: I don't know what the timings would be like, but the python code below will output in BAM rather than SAM, so you won't earn your PIs ire for using all that disk space, and I guess your processing code if going the be the slow bit.
  import pysam
  infile = pysam.AlignmentFile("input.bam")
  
  chunk_size = 10000000
  outfile_pattern = "output_segement%d.bam"
  
  chunk = 0
  reads_in_this_chunk = 0
  old_name = None
  outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)
  
  for read in infile.fetch(until_eof=True):
  
      if old_name != read.query_name and reads_in_this_chunk > chunk_size:
          reads_in_this_chunk = 0
          chunk += 1
          outfile.close()
          outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)
  
      outfile.write(read)
      old_name = read.query_name
      reads_in_this_chunk += 1
  
  outfile.close()
  
  
  
original_text: 4141: <p>I have a bam file with 1 billion alignment reads of which there are 700 million unique reads. I want to split the alignments into chunks for parallel-processing. Multi-alignments of the same read should be in the same chunk. I want to use pysam to solve the problem and would be very grateful for suggestions.</p>
   1310: <p>python will be to slow for this job. Here's a <code>awk</code> solution. One need to sort by read name and take track over the number of reads per chunk. If the number is reached the next chunk can only created of the read name is different.</p>
  
  <pre><code><span class="math-container">$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
    BEGIN { part=0; line=n }       
    /^@/ {header = header$0"\n"; next;} 
    { if( line&gt;=n &amp;&amp; $1!=last_read ) {part++; line=1; printf header &gt; part".sam"; print $0 &gt;&gt; part".sam" } 
      else { print $0 &gt;&gt; part".sam"; line++; } 
      last_read = $</span>1;
    }'
  </code></pre>
  
  <p>The output will be in sam format.</p>
  
  <p>You can controll the number of reads per chunk by setting the <code>n</code> value for the <code>awk</code> command.</p>
  
  <p>EDIT:</p>
  
  <p>A second version that will output in <code>bam</code>:</p>
  
  <pre><code><span class="math-container">$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" ' 
  BEGIN { part=0; line=n }  
  /^@/ {header = header$0"\n"; next;}  
  { if( line&gt;=n &amp;&amp; $1!=last_read ) {print part,line; part++; line=1;}  
    print line==1 ? header""$0 : $0 | "samtools view -b -o "part".bam" 
    last_read = $</span>1;
    line++;
  }
  </code></pre>
   235: <p>I don't know what the timings would be like, but the python code below will output in BAM rather than SAM, so you won't earn your PIs ire for using all that disk space, and I guess your processing code if going the be the slow bit.</p>
  
  <pre class="lang-py prettyprint-override"><code>import pysam
  infile = pysam.AlignmentFile("input.bam")
  
  chunk_size = 10000000
  outfile_pattern = "output_segement%d.bam"
  
  chunk = 0
  reads_in_this_chunk = 0
  old_name = None
  outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)
  
  for read in infile.fetch(until_eof=True):
  
      if old_name != read.query_name and reads_in_this_chunk &gt; chunk_size:
          reads_in_this_chunk = 0
          chunk += 1
          outfile.close()
          outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)
  
      outfile.write(read)
      old_name = read.query_name
      reads_in_this_chunk += 1
  
  outfile.close()
  
  </code></pre>
  