# Split FASTQ and matching BAM into matching chunks

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
n <- i modulo chunk_size
i <- i + 1


I could swap the above to iterate through the bam file and fetch from the fastq if that's easier.

• 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. Aug 17 '17 at 7:02
• 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. Aug 17 '17 at 9:33
• Yikes, your only option is something like what Michael Hall posted then. Aug 17 '17 at 9:34

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

"""Extracts the read ids from a fastq file."""
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 @

"""Extract the read ids from a BAM file."""

with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
# query_name is the query template name



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:
chunk_num = chunk_cycle.next()

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.