Chunk alignment in a name sorted bam for parallel processing

Question

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.

Answer 1

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()

Answer 2

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++;
}