Answers.Body: <p>Hmm, it's hard to think of a super efficient way of doing this (assuming the files aren't ordered the same - if they are then this whole answer is basically redundant). And also assuming the read ids for both files aren't a perfect intersection.</p>
  
  <p>Off the top of my head you probably want to build a set of read ids for the fastq file and another for the bam. Some python code to get started on that:</p>
  
  <pre><code>import pysam
  import itertools
  
  def get_read_id_fastq(ref_path):
      """Extracts the read ids from a fastq file."""
      read_ids = set()
      with open(ref_path, 'r') as ref:
          for line in ref:
              if line.startswith('@'):  # i.e if line is header
                  # split the line on spaces, take the first element, remove @
                  read_id = line.split(' ')[0].replace('@', '')
                  read_ids.add(read_id)
  
      return read_ids
  
  def get_read_id_bam(ref_path):
      """Extract the read ids from a BAM file."""
  
      read_ids = set()
      with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
          for read in ref:
              # query_name is the query template name
              read_ids.add(read.query_name)
  
      return read_ids
  
  fastq_ids = get_read_id_fastq(fastq_path)
  bam_ids = get_read_id_bam(bam_path)
  </code></pre>
  
  <p>Then take the intersection of these two sets and that's your common read ids. </p>
  
  <pre><code>common_ids = fastq_ids &amp; bam_ids
  </code></pre>
  
  <p>The next part will be a bit more involved. You will have to iterate through each file one at a time. I would suggest that for the first one you iterate through, create a running dictionary with the read id that was written as the key and the 'chuck number' it was written to as the value. You could create a cycle for your chunk size to manage this effectively</p>
  
  <pre><code>chunk_cycle = itertools.chain(*zip(range(chunk_size)*len(common_ids)))
  write_idx = {}
  </code></pre>
  
  <p>The next part will probably require you to put in some hard-fought times getting it to work. I'll put in some rough pseudo-code to give you an idea.</p>
  
  <pre><code>for line in file:
      read_id = # get line read_id
      if read_id in common_ids:
          chunk_num = chunk_cycle.next()
          write_idx[read_id] = chunk_num
          file_to_write_to = 'out.{}.bam'.format(chunk_num)
          # open this file or write to it if already open
  </code></pre>
  
  <p>When you go to do the next file when you find a read id is in the common set, you look up it's value in <code>write_idx</code> and this will give you the chunk number to write to. </p>
  
  <p>The reads wont be in the exact same order between the two files, but you could sort later if you needed it (not sure it would matter?).</p>
  
  <p>Hope this helps. Sorry I couldn't give you more, but this should give you a head start hopefully.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-17T23:39:50.330Z
Answers.Text: Thanks! In fact, I wouldn't need the first part, since every read in the BAM came from the FASTQ. If there are FASTQ records that don't exist in the BAM, I don't mind.
Answers.UserDisplayName: Scott Gigante
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-20T03:08:54.737Z
Answers.Text: Update: anyone interested can grab my implementation of this solution at https://github.com/scottgigante/nanopore-scripts/blob/master/split_bam_and_fasta.py
Answers.UserDisplayName: Scott Gigante
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-17T07:33:10.510Z
Answers.LastActivityDate: 2017-08-17T13:00:35.147Z
Answers.LastEditDate: 2017-08-17T13:00:35.147Z
Answers.UserDisplayName: Michael Hall
Body: <p>I am running a slow downstream analysis on a large set of nanopore reads (approx 3 million), and would like to split them into smaller chunks, run the analysis in massively parallel, and then recombine. Originally I just split the FASTQ into chunks, re-aligned each chunk, and then merged the output, but here I would like to use an existing alignment so I can compare results with existing analyses (i.e. the alignments must be the same).</p>
  
  <p>How can I efficiently split a FASTQ file and a BAM file giving the alignment of the FASTA file into chunks, ensuring that all of the reads in FASTQ chunk 1 are in BAM chunk 1, vice versa and so on?</p>
  
  <p>My FASTQ is approximately 45GB and my BAM is 33GB, so I would prefer to avoid storing one of the two files in memory if possible.</p>
  
  <p>EDIT: Here's some pseudocode of exactly what I'm trying to do:</p>
  
  <pre><code># input: in.bam, in.fastq, chunk_size
  i &lt;- 0
  for fastq_read in in.fastq:
      bam_read &lt;- extract fastq_read.read_name from in.bam
      n &lt;- i modulo chunk_size
      write fastq_read to out.n.fastq
      write bam_read to out.n.bam
      i &lt;- i + 1
  </code></pre>
  
  <p>I could swap the above to iterate through the bam file and fetch from the fastq if that's easier.</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-17T07:02:32.077Z
Comments.Text: Does your BAM file already happen to have the same ordering as the fastq files (some aligners produce that by default)? If that's not the case you're pretty much left with (1) sorting files or (2) hashing read names.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-17T09:33:27.753Z
Comments.Text: Unfortunately they're not in the same order, and additionally some reads exist multiple times in the BAM (supplementary alignments) and a handful of reads do not exist in the BAM at all.
Comments.UserDisplayName: Scott Gigante
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-17T09:34:40.317Z
Comments.Text: Yikes, your only option is something like what Michael Hall posted then.
Comments.UserDisplayName: Devon Ryan
Tags: <sam><fasta><nanopore><fastq><subset>
Title: Split FASTQ and matching BAM into matching chunks
ThreadText: # Split FASTQ and matching BAM into matching chunks
  
  (Asked by: username_0 on 2017-08-17)
  
  I am running a slow downstream analysis on a large set of nanopore reads (approx 3 million), and would like to split them into smaller chunks, run the analysis in massively parallel, and then recombine. Originally I just split the FASTQ into chunks, re-aligned each chunk, and then merged the output, but here I would like to use an existing alignment so I can compare results with existing analyses (i.e. the alignments must be the same).
  
  How can I efficiently split a FASTQ file and a BAM file giving the alignment of the FASTA file into chunks, ensuring that all of the reads in FASTQ chunk 1 are in BAM chunk 1, vice versa and so on?
  
  My FASTQ is approximately 45GB and my BAM is 33GB, so I would prefer to avoid storing one of the two files in memory if possible.
  
  EDIT: Here's some pseudocode of exactly what I'm trying to do:
  
  ```
  # input: in.bam, in.fastq, chunk_size
  i <- 0
  for fastq_read in in.fastq:
      bam_read <- extract fastq_read.read_name from in.bam
      n <- i modulo chunk_size
      write fastq_read to out.n.fastq
      write bam_read to out.n.bam
      i <- i + 1
  
  ```
  
  I could swap the above to iterate through the bam file and fetch from the fastq if that's easier.
  
  **Tags:** sam, fasta, nanopore, fastq, subset
  
  > **Comments:**
  > **username_1**: Does your BAM file already happen to have the same ordering as the fastq files (some aligners produce that by default)? If that's not the case you're pretty much left with (1) sorting files or (2) hashing read names.  
  > **username_0**: Unfortunately they're not in the same order, and additionally some reads exist multiple times in the BAM (supplementary alignments) and a handful of reads do not exist in the BAM at all.  
  > **username_1**: Yikes, your only option is something like what username_2 username_2 posted then.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2017-08-17. Score: 5 | Accepted answer
  
  Hmm, it's hard to think of a super efficient way of doing this (assuming the files aren't ordered the same - if they are then this whole answer is basically redundant). And also assuming the read ids for both files aren't a perfect intersection.
  
  Off the top of my head you probably want to build a set of read ids for the fastq file and another for the bam. Some python code to get started on that:
  
  ```
  import pysam
  import itertools
  
  def get_read_id_fastq(ref_path):
      """Extracts the read ids from a fastq file."""
      read_ids = set()
      with open(ref_path, 'r') as ref:
          for line in ref:
              if line.startswith('@'):  # i.e if line is header
                  # split the line on spaces, take the first element, remove @
                  read_id = line.split(' ')[0].replace('@', '')
                  read_ids.add(read_id)
  
      return read_ids
  
  def get_read_id_bam(ref_path):
      """Extract the read ids from a BAM file."""
  
      read_ids = set()
      with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
          for read in ref:
              # query_name is the query template name
              read_ids.add(read.query_name)
  
      return read_ids
  
  fastq_ids = get_read_id_fastq(fastq_path)
  bam_ids = get_read_id_bam(bam_path)
  
  ```
  
  Then take the intersection of these two sets and that's your common read ids.
  
  ```
  common_ids = fastq_ids & bam_ids
  
  ```
  
  The next part will be a bit more involved. You will have to iterate through each file one at a time. I would suggest that for the first one you iterate through, create a running dictionary with the read id that was written as the key and the 'chuck number' it was written to as the value. You could create a cycle for your chunk size to manage this effectively
  
  ```
  chunk_cycle = itertools.chain(*zip(range(chunk_size)*len(common_ids)))
  write_idx = {}
  
  ```
  
  The next part will probably require you to put in some hard-fought times getting it to work. I'll put in some rough pseudo-code to give you an idea.
  
  ```
  for line in file:
      read_id = # get line read_id
      if read_id in common_ids:
          chunk_num = chunk_cycle.next()
          write_idx[read_id] = chunk_num
          file_to_write_to = 'out.{}.bam'.format(chunk_num)
          # open this file or write to it if already open
  
  ```
  
  When you go to do the next file when you find a read id is in the common set, you look up it's value in `write_idx` and this will give you the chunk number to write to.
  
  The reads wont be in the exact same order between the two files, but you could sort later if you needed it (not sure it would matter?).
  
  Hope this helps. Sorry I couldn't give you more, but this should give you a head start hopefully.
  
  > **Comments:**
  > **username_0**: Thanks! In fact, I wouldn't need the first part, since every read in the BAM came from the FASTQ. If there are FASTQ records that don't exist in the BAM, I don't mind.  
  > **username_0**: Update: anyone interested can grab my implementation of this solution at https://github.com/username_0/nanopore-scripts/blob/master/split_bam_and_fasta.py  
  
  ---
  
  