1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Could you confirm what you mean by "chunks". Do you mean a genetic locus, e.g. X nucleotides is one "chunk". Alternatively, do you mean subsets of the read depth? $\endgroup$ – Michael G. Feb 17 at 21:53
  • 1
    $\begingroup$ What kind of parallel processing do you plan? Is is unusual to split by read names instead of loci for parallel processing. $\endgroup$ – finswimmer Feb 18 at 4:26
  • 1
    $\begingroup$ What have you tried? In particular, what python code using pysam have you written to do this? The steps are mostly (A) name sort and (B) iterate over and write a lot of files. $\endgroup$ – Devon Ryan Feb 18 at 9:23
  • 1
    $\begingroup$ I think it is a bit outside of the scope of this site to write a program for you. Perhaps if you show us what you have tried so far, it would be more possible for us to assist you. In the question's current state, you are asking us to write you a program which incorporates pysam. $\endgroup$ – d_kennetz Feb 18 at 15:23
  • 1
    $\begingroup$ Thank you all for your help. Did some search. I think that for fastq files or sam files, it should be easy based on the class here: blopig.com/blog/2016/08/… However, for bam files it's more challenging as I don't know how to use 'seek' in pysam. Or maybe just treat bam files as special zipped files and use general gzip reader for this purpose? Or is it possible to do so by using pysam? Thank you all. $\endgroup$ – jerry00 Feb 18 at 17:17
2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Thank you, this works great! What I can think of to optimize the code a little bit is 1) I don't need to check the name of each alignment. Only after reading a specified number of alignments, I check the read name, read again until I find a read with a different name. 2) It seems un-doable for pysam to pass the start and stop positions of chunks to multiprocessing.pool. $\endgroup$ – jerry00 Feb 19 at 16:12
  • $\begingroup$ Not sure how can I accept this answer. $\endgroup$ – jerry00 Feb 19 at 16:17
  • $\begingroup$ @jenny00 the comparison of read names will take next to no time. The slow bit is reading the read into python, once its in manipulations are fast. And as for 2: No, pysam is not compatible with python's simple multiprocessing API because you can only pass picklable objects through to the child processes and AlignedSegments (or reads) are binary objects - they are actually a thin adaptor layer on the C++ based htslib object. $\endgroup$ – Ian Sudbery Feb 20 at 18:16
  • $\begingroup$ Thanks, I just noticed that. So I have to put the necessary information from each alignment into a dictionary and pass dictionaries down for parallel processing. It does not save much time compared to just using one process. $\endgroup$ – jerry00 Feb 21 at 2:02
  • $\begingroup$ If you are using pythons multiprocessing.Pool infrastructure, you might find that the overhead associated is high - each submission to a Pool causes the creation of a new python process, which does its processing and then the process is closed down (i think). You might like to consider the more complex to implement asyncronous worker and processing queue type set up. Also try increasing the size of the chunks - the fewer chunks, the small the number of process spin ups and close downs. $\endgroup$ – Ian Sudbery Feb 21 at 14:47
2
$\begingroup$

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++;
}
$\endgroup$
  • $\begingroup$ Thank you so much! This works great as well, but I want to use pysam so the code can be integrated into a longer function. I don't even need to save the chunked alignments, but just pass them for parallel processing. It could be quite useful if I want to change my design to split and save a bam file into smaller bam files. I think that many other people may find this solution helpful as well. $\endgroup$ – jerry00 Feb 19 at 16:16
  • $\begingroup$ I was considering awk too. The problem with Perl/Python is the code will process the sequence data "[horizontal] line by [horizontal] line", when awk can do this "[vertical] line by [vertical] line " $\endgroup$ – Michael G. Feb 20 at 17:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.