Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-19T16:12:35.207Z
Answers.Text: 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.
Answers.UserDisplayName: jerry00
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-19T16:17:37.533Z
Answers.Text: Not sure how can I accept this answer.
Answers.UserDisplayName: jerry00
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-20T18:16:02.697Z
Answers.Text: @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.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-21T02:02:53.437Z
Answers.Text: 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.
Answers.UserDisplayName: jerry00
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-21T14:47:52.703Z
Answers.Text: 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.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-22T04:41:34.643Z
Answers.Text: 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.
Answers.UserDisplayName: jerry00
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-19T11:55:18.650Z
Answers.LastActivityDate: 2019-02-19T11:55:18.650Z
Answers.UserDisplayName: Ian Sudbery
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-19T16:16:05.110Z
Answers.Text: 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.
Answers.UserDisplayName: jerry00
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-20T17:35:03.250Z
Answers.Text: 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 "
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-02-19T07:02:12.363Z
Answers.LastActivityDate: 2019-02-19T17:59:03.320Z
Answers.LastEditDate: 2019-02-19T17:59:03.320Z
Answers.UserDisplayName: finswimmer
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-17T21:53:01.183Z
Comments.Text: 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?
Comments.UserDisplayName: M__
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-17T23:30:48.900Z
Comments.Text: 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.
Comments.UserDisplayName: jerry00
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T04:26:55.460Z
Comments.Text: What kind of parallel processing do you plan? Is is unusual to split by read names instead of loci for parallel processing.
Comments.UserDisplayName: finswimmer
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T09:23:26.147Z
Comments.Text: 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.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T15:23:01.620Z
Comments.Text: 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.
Comments.UserDisplayName: d_kennetz
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T17:17:57.797Z
Comments.Text: 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.
Comments.UserDisplayName: jerry00
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T17:21:44.953Z
Comments.Text: 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.
Comments.UserDisplayName: jerry00
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T18:03:50.007Z
Comments.Text: Which aligner was used? `bwa` stored alternative mapping positions in the XA tag.
Comments.UserDisplayName: finswimmer
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T18:20:00.043Z
Comments.Text: 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
Comments.UserDisplayName: jerry00
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-02-18T22:46:14.657Z
Comments.Text: 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.
Comments.UserDisplayName: d_kennetz
Tags: <python><sam><pysam><multithreading>
Title: Chunk alignment in a name sorted bam for parallel processing
ThreadText: # 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.
  
  ```py
  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 "  
  
  ---
  
  