Split FASTQ and matching BAM into matching chunks

Question

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.

Answer 1

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.