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.

  • 1
    $\begingroup$ 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. $\endgroup$ – Devon Ryan Aug 17 '17 at 7:02
  • $\begingroup$ 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. $\endgroup$ – Scott Gigante Aug 17 '17 at 9:33
  • $\begingroup$ Yikes, your only option is something like what Michael Hall posted then. $\endgroup$ – Devon Ryan 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

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('@', '')

    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

    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.

  • 1
    $\begingroup$ 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. $\endgroup$ – Scott Gigante Aug 17 '17 at 23:39
  • 1
    $\begingroup$ Update: anyone interested can grab my implementation of this solution at github.com/scottgigante/nanopore-scripts/blob/master/… $\endgroup$ – Scott Gigante Aug 20 '17 at 3:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.