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

CC BY-SA 4.0
2019-02-19T16:12:35.207Z
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.
jerry00
CC BY-SA 4.0
2019-02-19T16:17:37.533Z
Not sure how can I accept this answer.
jerry00
CC BY-SA 4.0
2019-02-20T18:16:02.697Z
@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.
Ian Sudbery
CC BY-SA 4.0
2019-02-21T02:02:53.437Z
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.
jerry00
CC BY-SA 4.0
2019-02-21T14:47:52.703Z
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.
Ian Sudbery
CC BY-SA 4.0
2019-02-22T04:41:34.643Z
Thanks. Increasing the chunk size and using map_async help. But it is only twice faster using four cores, compared to using a single core.
jerry00
CC BY-SA 4.0
2019-02-19T11:55:18.650Z
2019-02-19T11:55:18.650Z
Ian Sudbery
<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>

CC BY-SA 4.0
2019-02-19T16:16:05.110Z
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.
jerry00
CC BY-SA 4.0
2019-02-20T17:35:03.250Z
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 "
M__
CC BY-SA 4.0
2019-02-19T07:02:12.363Z
2019-02-19T17:59:03.320Z
2019-02-19T17:59:03.320Z
finswimmer
<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>

CC BY-SA 4.0
2019-02-17T21:53:01.183Z
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?
M__
CC BY-SA 4.0
2019-02-17T23:30:48.900Z
I want to split the alignment into small segments, e.g., 10M alignments per segment. The expected input bam file is name sorted, so roughly 1-10M alignments in the first chunk,10M - 20M alignments in the second chunk and so on. The problem is that I don't want to put the alignments from the same read into different chunks. Thank you.
jerry00
CC BY-SA 4.0
2019-02-18T04:26:55.460Z
What kind of parallel processing do you plan? Is is unusual to split by read names instead of loci for parallel processing.
finswimmer
CC BY-SA 4.0
2019-02-18T09:23:26.147Z
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.
Devon Ryan
CC BY-SA 4.0
2019-02-18T15:23:01.620Z
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.
d_kennetz
CC BY-SA 4.0
2019-02-18T17:17:57.797Z
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: https://www.blopig.com/blog/2016/08/processing-large-files-using-python-part-duex/

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.
jerry00
CC BY-SA 4.0
2019-02-18T17:21:44.953Z
Basically, I want to know where all the alignments of a read mapped to, so I decide to process read by read. Here I distinguish reads from alignments.
jerry00
CC BY-SA 4.0
2019-02-18T18:03:50.007Z
Which aligner was used? `bwa` stored alternative mapping positions in the XA tag.
finswimmer
CC BY-SA 4.0
2019-02-18T18:20:00.043Z
These are RNA-seq data and I use STAR for alignment. I can see four tags, e.g., NH:i:6 HI:i:6 AS:i:49 nM:i:0
jerry00
CC BY-SA 4.0
2019-02-18T22:46:14.657Z
You can sort by read name instead of position. This is relevant for whatever you are doing. You should read [this thread](https://sourceforge.net/p/samtools/mailman/samtools-devel/thread/5A43E82B-1149-4F9C-B02E-790322FBB74E%40sanger.ac.uk/#msg30888139) to understand the bit flags for multiple primary alignments. It should be helpful.
d_kennetz
<python><sam><pysam><multithreading>
Chunk alignment in a name sorted bam for parallel processing
# Chunk alignment in a name sorted bam for parallel processing

(Asked by: username_0 on 2019-02-17)

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.

**Tags:** python, sam, pysam, multithreading
> **Comments:**
> **username_1**: 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?  
> **username_0**: I want to split the alignment into small segments, e.g., 10M alignments per segment. The expected input bam file is name sorted, so roughly 1-10M alignments in the first chunk,10M - 20M alignments in the second chunk and so on. The problem is that I don't want to put the alignments from the same read into different chunks. Thank you.  
> **username_2**: What kind of parallel processing do you plan? Is is unusual to split by read names instead of loci for parallel processing.  
> **username_3**: 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.  
> **username_4**: 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.  
> **username_0**: 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: https://www.blopig.com/blog/2016/08/processing-large-files-using-python-part-duex/  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.  
> **username_0**: Basically, I want to know where all the alignments of a read mapped to, so I decide to process read by read. Here I distinguish reads from alignments.  
> **username_2**: Which aligner was used? `bwa` stored alternative mapping positions in the XA tag.  
> **username_0**: These are RNA-seq data and I use STAR for alignment. I can see four tags, e.g., NH:i:6 HI:i:6 AS:i:49 nM:i:0  
> **username_4**: You can sort by read name instead of position. This is relevant for whatever you are doing. You should read [this thread](https://sourceforge.net/p/samtools/mailman/samtools-devel/thread/5A43E82B-1149-4F9C-B02E-790322FBB74E%40sanger.ac.uk/#msg30888139) to understand the bit flags for multiple primary alignments. It should be helpful.  

---

## 2 answers

### Answer by username_5 on 2019-02-19. Score: 2

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

> **Comments:**
> **username_0**: 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.  
> **username_0**: Not sure how can I accept this answer.  
> **username_5**: @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.  
> **username_0**: 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.  
> **username_5**: 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.  
> **username_0**: Thanks. Increasing the chunk size and using map_async help. But it is only twice faster using four cores, compared to using a single core.  

---

### Answer by username_2 on 2019-02-19. Score: 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++;
}

> **Comments:**
> **username_0**: 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.  
> **username_1**: 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 "  

---

