Biopython: How to split a very large .fasta file into roughly equally sized files?

[dmenning]

I do a lot of primer design and checking so I have written scripts to parse through .fasta files. In some cases these are fairly large 10+gb. I would like to break these up into smaller more manageable files. I was looking for a method to break up large files by size and came along this:
http://stackoverflow.com/questions/291740/how-do-i-split-a-huge-text-file-in-python

    import logging.handlers
    log = logging.getLogger()
    fh = logging.handlers.RotatingFileHandler("verylarge.fasta", 
        maxBytes=2**20*100, backupCount=100) 
    log.addHandler(fh)
    log.setLevel(logging.INFO)
    f = open("D://biglog.txt")
    while True:
        log.info(f.readline().strip())

My question is, can I use this to break up a large .fasta file into smaller, equally sized files? I think it would just read the text and possibly cut off the DNA sequence before finishing if it got close to the max size and am not sure how to handle this using Biopython. Is there a way to get a running output file size, current size of data to be written, then if there is enough room in the file to write the data, and move on?

Thanks

Damian

[peterjc]

It sounds like you want to divide the large FASTA file into chunks, splitting at a record boundary (i.e. never breaking up any of the individual sequence entries).

I would script this setting in advance either the number of records per file (e.g. 1000 sequences per chunk), or a maximum file size on disk per chunk file. There a doubtless many existing scripts out there doing this already - do you have a particular reason to write your own?

An example (not using Biopython) for the fixed number of records per file approach: https://github.com/peterjc/pico_galaxy/blob/e8a3b88c1d33c75b99edc17c27b9fc4d1260696e/tools/protein_analysis/seq_analysis_utils.py#L102

[peterjc]

Also if you are dealing with 10+ gb files, I would avoid parsing the records at all and simply use a line-by-line state machine approach (is it a header line line.startswith(">") or sequence line).

[dmenning]

Thanks for the info Peter. The script in the link is much more elegant that what I was using.

One of the reasons I was using Biopython was the ability to search a DNA sequence for specific primers using regular expressions and if found to write the description information and sequence data to another file using Biopython. It seemed easier since Biopython makes access to that info simple and keeps all the data together.

[dmenning]

Ok, so this is what I wrote to split a large .fasta file in order to keep the data intact and not split by lines. The original file "Bony Fish 18S 160914.fasta" is a little over 5GB so it should end up with 10 files that are approximately 500 MB each. The problem I am running into now is that it is getting hung up on the last file. It is splitting everything properly so I see through .fasta.9 but the last fasta file should be "Bony Fish 18S 160914 split.fasta" is not finishing. I have tried changing the output file extension to .txt but am getting the same problem. Any ideas?

import logging.handlers
from Bio import SeqIO
log = logging.getLogger()
fh = logging.handlers.RotatingFileHandler("Bony Fish 18S 160914 split.fasta", maxBytes=2**20*500, backupCount=100) 
# 500 MB each, up to 100 files
log.addHandler(fh)
log.setLevel(logging.INFO)
f = open("Bony Fish 18S 160914.fasta")
while True:
    for seq_record in SeqIO.parse(f, "fasta"):
        desc = str(seq_record.description)
        sequence = str(seq_record.seq)
        to_write = desc + "\n" + sequence
        log.info(to_write)

[peterjc]

I've not used RotatingFileHandler but my guess is you need to flush the data to disk, either via .flush() or implicitly via .close() to trigger the final batch of data gets written to disk immediately. That ought to happen anyway when the script ends (I'm guessing you might be showing only part of your script here).

[dmenning]

I've redone my script without RotatingFileHandler and it is working fine now. Thanks for the help.

import os, sys
from os import listdir
from stat import *
from Bio import SeqIO

#checks to make sure input file is in the folder
try:
    name = raw_input("\nEnter file name with fasta sequences you want to split: ")
    handle = open(name, 'r')
except:
    print "File does not exist in folder! Check file name and extension."
    quit()

#output file(s)
n = 1
outmid = os.path.splitext(name)[0]+" split "+str(n)+".fasta"
output_handle = open(outmid, "a")
print "\n"

#main script
for seq_record in SeqIO.parse(handle, "fasta"):
    output_handle.write(">"+str(seq_record.description) + "\n" + str(seq_record.seq) + "\n")
    size = os.stat(outmid).st_size
    if size > 500000000:
        #user feedback
        print str(outmid) + " file size = \t" + str(size)
        n += 1
        output_handle.close()
        outmid = os.path.splitext(name)[0]+" split "+str(n)+".fasta"
        output_handle = open(outmid, "a")

print str(outmid) + " file size = \t" + str(size)        
output_handle.close()

[dmenning]

Ran into a memory error with a 13GB fasta file. Rewrote the main potion below without SeqIO.

#main script
for line in handle:
    if line.startswith(">"):
        size = os.stat(outmid).st_size
        if size > 500000000:
            #user feedback
            print str(outmid) + "\tfile size = " + str(size)
            n += 1
            output_handle.close()
            outmid = os.path.splitext(name)[0]+" split "+str(n)+".fasta"
            output_handle = open(outmid, "a")
    output_handle.write(str(line))

[peterjc]

I was going to say since you only use strings and don't care about line wrapping in the output, use:

from Bio.SeqIO.FastaIO import SimpleFastaParser
for title, seq in SimpleFastaParser(handle):
    # etc

Anyway, glad you've got something working.