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.
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.
#!/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()
fl()
function only operates over half the random sequence and does nothing with the other half? \$\endgroup\$