# Generating DNA sequences and looking for correlations

I've written a script to generate DNA sequences and then count the appearance of each step to see if there is any long range correlation.

My program runs really slow for a length 100000 sequence 100 times replicate. I already run it for more than 100 hours without completion.

#!/usr/bin/env python

import sys, random
import os
import math

length = 10000

initial_p = {'a':0.25,'c':0.25,'t':0.25,'g':0.25}

tran_matrix = {'a': {'a':0.495,'c':0.113,'g':0.129,'t':0.263},
'c': {'a':0.129,'c':0.063,'g':0.413,'t':0.395},
't': {'a':0.213,'c':0.495,'g':0.263,'t':0.029},
'g': {'a':0.263,'c':0.129,'g':0.295,'t':0.313}}

def fl():
def seq():
def choose(dist):
r = random.random()
sum = 0.0
keys = dist.keys()
for k in keys:
sum += dist[k]
if sum > r:
return k
return keys[-1]
c = choose(initial_p)
sequence = ''
for i in range(length):
sequence += c
c = choose(tran_matrix[c])
return sequence

sequence = seq()
# This program takes a DNA sequence calculate the DNA walk score.
#print sequence
#print len
u = 0
ls = []
for i in sequence:
if i == 'a' :
#print i
u = u + 1
if  i == 'g' :
#print i
u = u + 1
if  i== 'c' :
#print i
u = u - 1
if  i== 't' :
#print i
u = u - 1
#print u
ls.append(u)

#print ls
l = 1
f = []
for l in xrange(1,(length/2)+1):
lchange =1
sumdeltay = 0
sumsq = 0
for i in range(1,length/2):
deltay = ls[lchange + l ] - ls[lchange]
lchange = lchange + 1
sq = math.fabs(deltay*deltay)
sumsq = sumsq + sq
sumdeltay = sumdeltay + deltay
f.append(math.sqrt(math.fabs((sumsq/length/2) - math.fabs((sumdeltay/length/2)*(sumdeltay/length/2)))))
l = l + 1
return f

def performTrial(tries):
distLists = []
for i in range(0, tries):
fl()
distLists.append(fl())
return distLists

def main():
tries = 10
distLists = performTrial(tries)
#print distLists
#print distLists[0][0]
averageList = []
for i in range(0, length/2):
total = 0
for j in range(0, tries):
total += distLists[j][i]
#print distLists
average = total/tries
averageList.append(average)
# print total
return averageList

out_file = open('Markov1.result', 'w')
result = str(main())
out_file.write(result)
out_file.close()

• Could you explain what this is supposed to do in plain English? 1) At the moment, this is a bit long to read without explanation, 2) It may help you identify the problem yourself, and 3) Stop people taking a stab at stuff that looks like it may be wrong Dec 13, 2012 at 21:26
• Hi, the first function is seq(), it randomly generate DNA sequences based on the transition matrix follow markov chain, each step will be either a,c,t,or g. the length is ten thousands. after the dna sequence generated, the u will be calculated. u is the total score by DNA walk, for example, at each step, if there is a|g, u + 1, if there is a c|g, u -1. Therefore we will have u from step 1 to step 10 thousands. Then we calculate the fluctuation for u from l= 1 step to l =5000 step, to see if there is a long range correlation exist. The performTrial() is using to do replicate for fl() function.
– Frank
Dec 13, 2012 at 21:35
• Then the main() calculate the average score from the replications output from performTrial().
– Frank
Dec 13, 2012 at 21:38
• The f value return from fl() is the root mean fluctuation calculated from u on different length(step) scale.
– Frank
Dec 13, 2012 at 21:42
• @Frank Is there a reason your fl() function only operates over half the random sequence and does nothing with the other half?
– Iguananaut
Dec 13, 2012 at 22:52

One word: numpy.

A few more words:

numpy is a library that includes highly-tuned algorithms (often in Fortran, C, and C++) specifically for the purposes of accelerating numerical computations over matrices and other arrays. If it's at all applicable to what you're trying to do, it's almost certainly the best answer. Even when numpy itself isn't appropriate, look at the other components of scipy (the larger project it's part of).

At the moment, the numpy site seems to be down, and the scipy site seems to be responding but broken. The docs wiki is working, at least. (Please don't take this as an indication that numpy is some fly-by-night project; almost everyone who does numerical or scientific computing in Python relies on it. I'm sure it'll be fixed shortly.)

However, there are other possibilities for other use cases:

• Analyze your algorithms and make sure the problem isn't what you're doing, rather than how you're doing it. There's no point making everything 30% faster, or even 95%, when the real problem is that you're doing 1000000 loops instead of 1000.
• Profile your code, and don't worry about the parts that aren't slow or aren't called often.
• Whenever it's possible to use a built-in function, even if it's a bit more awkward, it's often faster. (For example, sticking an iterator in a collections.deque of size 0 is at least 30% faster than any pure-Python way of walking and disposing the iterator.)
• Cython gives you a way to write almost-Python code that gets compiled into C (and then into binary code).
• You can port (bits of) your code to C explicitly.
• Sometimes just running in PyPy (or maybe Jython or IronPython) instead of the normal CPython makes a big difference.
• Replace explicit lists with generators, and explicit list-building loops that can't be replaced by generators with list comprehensions.
• If you need that last 2%, there are all kinds of tricks, like rebinding members of objects to local variables, that aren't always more trouble than they're worth.

A couple more words on Cython:

The first step to accelerating some slow critical code with Cython is just to move it as-is into Cython. If that doesn't help enough, try adding explicit types to variables and turning function definitions from def into cdef. If that's still not enough, then you have to actually learn how to use Cython effectively—but you may get enough win that you don't even have to understand why it's helping. (Although obviously you should learn anyway.)

• I have only skimmed over his code, but some of the operations seem to be hard/impossible to vectorize, so for this particular case numpy may not be the answer. But still +1
– Jaime
Dec 13, 2012 at 21:36
• They should also take a look at PyMC for MCMC implementations.
– Iguananaut
Dec 13, 2012 at 21:39
• @Frank You might want to also look at Python for Data Analysis (note: I have no financial incentive in plugging this book). It's a good book for learning how to do the kinds of things you're trying to do with Python. In short, you're doing it all wrong. This should set you straight. See also if you can find some workshops in your area to learn using Python for scientific data analysis.
– Iguananaut
Dec 13, 2012 at 21:41
• @Iguananaut: You might want to write a separate answer if you've got specifics that are more appropriate to his use case than my general answer, because it may be more useful to a later reader (also, I can +1 you again that way).
– abarnert
Dec 13, 2012 at 21:44
• @Iguananaut (and OP) I think profiling the existing code would be the best first step -- and it could be very useful during the likely iterative process of speeding various parts of it up.
– martineau
Dec 13, 2012 at 22:53

Your algorithm is quadratic in length and linear in tries. The double loop in fl() dominates. O(length^2 * tries)

After Glenn's improvements (removing superfluous fl() halves the running time of main() and removing the function call from the inner loop further reduces the running time almost by another third. Your mileage may vary.) fl() runs about 6-7 secs. main totals ~70 secs.

increasing the length by ten (10,000 to 100,000), and increasing the tries by ten (10 to 100) (which I understand from "100 times replicate,") the running time should be about ~70,000 secs which is about 20hrs.

implementing the algorithm in C or another compiled language, Cython, (maybe PyPy, it can do JIT for pure python, as far as I know, which your inner loop is) should further half that time to about 8/10 hrs, which is acceptable.

Moreover,

lchange =1
sumdeltay = 0
sumsq = 0
for i in range(1,length/2):
deltay = ls[lchange + l ] - ls[lchange]
lchange = lchange + 1
sq = math.fabs(deltay*deltay)
sumsq = sumsq + sq
sumdeltay = sumdeltay + deltay
ff = math.sqrt(math.fabs((sumsq/length/2) - (sumdeltay/length/2)*(sumdeltay/length/2)))


portion of your code is read-only, and paralellizable. you can calculate ls beforehand and save to a text file. put the code in a separate python script. spawn new processes passing each l (or some interval of l values eg, 500-1000 depending on no of processes you want to spawn) as command-line parameters. and saving the results in text files, e.g. f_345.txt, if run with l=345. later combining the results. this can be done in python or bash or some other scripting language.

you can, possibly, further half the run time on a quad-core desktop. or split the work to multiple computers (in a work or school environment). also you will not lose partial results.

• I also think, since you mentioned correlation, that sumsq/length/2 should be sumsq/(length/2). / is left associative in python. But I guess you want to divide by (length/2) Dec 17, 2012 at 9:31
• another bad style or less probably a bug is the use of local variables in the loop lchange is always equal to i so why not just use i, if the name means something it can be called lchange then. l = l + 1 has no effect. l iterates over xrange(1,(length/2)+1) with or without it. if you are trying to iterate over odd numbers, which i do not think you do, you should use xrange(1,(length/2)+1, 2). Dec 18, 2012 at 7:58

Concerning just a little part of your code

u = 0
ls = []
for i in sequence:
u += (1 if i in 'ag' else (-1))
ls.append(u)

print ls


OR

def gen(L,u = 0):
for i in L:
u += (1 if i in 'ag' else (-1))
yield u

print list(gen(sequence))


At the risk of being frivolous, two suggestions with the posted code:

1. Did you really mean fl() on both occasions?

for i in range(0, tries):
fl()
distLists.append(fl())

2. There's no need to use math.fabs on values guaranteed to be positive!