How do I parallelize this code in python (GC calculator from FASTA)?

Question

System Monitor during process I am a novice when it comes to programming. I've worked through the book Practical Computing for Biologists and am playing around with some slightly more advanced concepts.

I've written a Python (2.7) script which reads in a .fasta file and calculates GC-content. The code is provided below.

The file I'm working with is cumbersome (~ 3.9 Gb), and I was wondering if there's a way to take advantage of multiple processors, or whether it would be worth-while. I have a four-core (hyperthreaded) Intel i-7 2600K processor.

I ran the code and looked at system resources (picture attached) to see what the load on my CPU is. Is this process CPU limited? Is it IO limited? These concepts are pretty new to me. I played around with the multiprocessing module and Pool(), to no avail (probably because my function returns a tuple).

Here's the code:

def GC_calc(InFile):
    Iteration = 0
    GC = 0
    Total = 0
    for Line in InFile:
        if Line[0] != ">":
            GC = GC + Line.count('G') + Line.count('C')
            Total = Total + len(Line)
            Iteration = Iteration + 1
            print Iteration
    GCC = 100 * GC / Total
    return (GC, Total, GCC)


InFileName = "WS_Genome_v1.fasta"
InFile = open(InFileName, 'r')
results = GC_calc(InFile)
print results

Answer 1

Currently, the major bottleneck of your code is print Iteration. Printing to stdout is really really slow. I would expect a major performance boost if you remove this line, or at least, if you absolutely need it, move it to another thread. However, thread management is an advanced topic and I would advice not to go into it right now.

Another possible bottleneck is the fact that you read data from a file. File IO can be slow, especially if you have a single HDD on your machine. With a single HDD you won't need to use multiprocessing at all, because you won't be able to provide enough data to processor cores. Performance-oriented RAIDs and SSDs can help here.

The final comment is to try to use grep and similar text-operating programs instead of python. They got through a decades of optimizing and have a good chance to work way more faster. There's a bunch of questions on SO where grep outperforms python. Or at least you can filter out FASTA headers before you pass data to the script:

$ grep "^[>]" WS_Genome_v1.fasta | python gc_calc.py

(Took grep "^[>]" from here.) In this case you shouldn't open files in the script, but rather read lines from sys.stdin, almost like you do it now.

Answer 2

Basically you're counting the number of C and G in every line and you're calculating the length of the line. Only at the end you calculate a total.

Such a process is easy to do in parallel, because the calculation for each line is independent of the others.

Assuming the calculations are done in CPython (the one from python.org), threading won't improve performance much because of the GIL.

These calculations could be done in parallel with multiprocessing.Pool. Processes don't share data like threads do. And we don't want to send parts of a 3.9 GB file to each worker process! So you want each worker process to open the file by itself. The operating system's cache should take care that pages from the same file aren't loaded into memory multiple times.

If you have N cores, I would create the worker function so as to process every N-th line, with an offset.

def worker(arguments):
    n = os.cpu_count() + 1
    infile, offset = arguments
    with open(infile) as f:
        cg = 0
        totlen = 0
        count = 1
        for line in f:
            if (count % n) - offset == 0:
                if not line.startswith('>'):
                    cg += line.count('C') +
                          line.count('G')
                    totlen += len(line)
            count += 1
    return (cg, totlen)

You could run the pool like this;

import multiprocessing as mp
from os import cpu_count

pool = mp.Pool()
results = pool.map(worker, [('infile', n) for n in range(1, cpu_count()+1)])

By default, a Pool creates as many workers as the CPU has cores.

The results would be a list of (cg, len) tuples, which you can easily sum.

Edit: updated to fix modulo zero error.

Answer 3

I now have a working code for parallelization (special thanks to @Roland Smith). I just had to make two small modifications to the code, and there's a caveat with respect to how .fasta files are structured. The final (working) code is below:

###ONLY WORKS WHEN THERE ARE NO BREAKS IN SEQUENCE LINES###

def GC_calc(arguments):
    n = mp.cpu_count()
    InFile, offset = arguments
    with open(InFile) as f:
        GC    = 0
        Total = 0
        count = 0
        for Line in f:
            if (count % n) - offset == 0:
                if Line[0] != ">":
                    Line = Line.strip('\n')
                    GC += Line.count('G') + Line.count('C')
                    Total += len(Line)
            count += 1
    return (GC, Total)

import time
import multiprocessing as mp

startTime = time.time()
pool = mp.Pool()
results = pool.map(GC_calc, [('WS_Genome_v2.fasta', n) for n in range(1, mp.cpu_count()+1)])
endTime = time.time()
workTime = endTime - startTime
#Takes the tuples, parses them out, adds them
GC_List = []
Tot_List = []
# x = GC count, y = total count: results = [(x,y), (x,y), (x,y),...(x,y)]
for x,y in results:
    GC_List.append(x)
    Tot_List.append(y)

GC_Final = sum(GC_List)
Tot_Final = sum(Tot_List)

GCC = 100*float(GC_Final)/float(Tot_Final)

print results
print
print "Number GC = ", GC_Final
print "Total bp = ", Tot_Final
print "GC Content = %.3f%%" % (GCC)
print
endTime = time.time()
workTime = endTime - startTime
print "The job took %.5f seconds to complete" % (workTime)

The caveat is that the .fasta files cannot have breaks within the sequences themselves. My original code didn't have an issue with it, but this code would not work properly when sequence was broken into multiple lines. That was simple enough to fix via the command line.

I also had to modify the code in two spots:

n = mp.cpu_count()

and

count = 0

Originally, count was set to 1 and n was set to mp.cpu_count()+1. This resulted in inaccurate counts, even after the file correction. The downside is that it also allowed all 8 cores (well, threads) to work. The new code only allows 4 to work at any given time.

But it DID speed up the process from about 23 seconds to about 13 seconds! So I'd say it was a success (except for the amount of time it took to correct the original .fasta file).