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