Analyzing pair by pair DNA sequences with for loops in Python

Here is a code I would need to make far faster: it is made to analyse more than 3000 DNA sequences of more than 100 000 characters each.

The matter is I have to compare it by pair, so that would be around 3000*2999/2 = 4498500 calculation and each takes 1 to 2 seconds...

Could you think of another way of doing it ? I've seen that there are algorithms that are used to go faster (like Boyer-Moore), but could we apply it here ? Or replace the for loops with something faster ? Any idea would be more than welcome.

import numpy as np
import time
import random

def calculate_distance(genom1, genom2):
# in : two strings (sequences)
# out : score between those 2 strings

# example : 'AATTGCAT' compared to 'AAAACGTC'
# => 'AA' and 'AA' = 2
# => 'TT' and 'AA' = 0
# => 'GC' and 'CG' = 2
# => 'AT' and 'TC' = 1

score = 0
for i in range(0, len(genom1), 2):
if (genom1[i] in genom2[i:i+2]) or (genom1[i+1] in genom2[i:i+2]):
if sorted(genom1[i:i+2]) == sorted(genom2[i:i+2]):
score += 2
else :
score += 1
return score

def build_matrix(sequences, N):
# in : list of lists
# out : matrix of scores between each pair of sequences
matrix = np.zeros((N, N))
for i in range(N):
for j in range(i, N):
matrix[i][j] = calculate_distance(sequences[i], sequences[j])
return matrix

def test(len_seq, N):
sequences = []
for i in range(N):
sequences.append(''.join(random.choice(['0','A', 'T', 'C', 'G']) for x in range(len_seq)))

start = time.clock()
matrix = build_matrix(sequences, N)
elapsed = time.clock() - start
print('Run time for ' + str(N) + ' sequences of ' + str(len_seq) + ' characters : computed in ' + str(elapsed) + ' seconds')
return matrix

test(10**6, 2)


Returns :

Run time for 2 sequences of 1000000 characters : computed in 2.742817 seconds

• How/where do you keep the data. How are you loading it? Dec 6, 2016 at 15:38
• Maybe you should give it a try? en.wikipedia.org/wiki/Trie Dec 6, 2016 at 16:54
• @enedil: “Maybe you should give it a trie?” – FTFY! Dec 11, 2016 at 14:28

I notice something small, well you know the corner cases. You know that a DNA string will yield a score of that is the same as length. Obviously ^^. So don't check those cases.

You make 4 iterations in your code. While only one is interesting.

A big improvement in runtime could be to change:

def build_matrix(sequences, N):
# in : list of lists
# out : matrix of scores between each pair of sequences
matrix = np.zeros((N, N))
for i in range(N):
for j in range(i, N):
matrix[i][j] = calculate_distance(sequences[i], sequences[j])
return matrix


to this:

def build_matrix(sequences, N):
matrix = np.zeros((N, N))
for i in range(N):
for j in range(i+1, N):
matrix[i][j] = calculate_distance(sequences[i], sequences[j])
return matrix


def get_test_arr():
parts = ['0', 'A', 'T', 'C', 'G']
return ["".join(
sorted([random.choice(parts) for _ in range(2)]))
for _ in range(5*10**5)]


you could make them np arrays as:

a = np.array(get_test_arr(), dtype=object)
b = np.array(get_test_arr(), dtype=object)


and have your distance calculation as:

def calculate_distance1(genom1, genom2):
score = 0
for g1, g2 in zip(genom1, genom2):
if (g1[0] in g2) or (g1[1] in g2):
score += 2 if g1 == g2 else 1
return score


In the problem you've provided each induvidual pairs order does not seem to matter, so load them as such. If you do that, you don't have to sort the in the calculation function.

If you don't, you'll resort them later. Again and again. Them more samples you use, there more time is spent doing this.

Times:

0.575332 s for one loop in 100 iterations.
0.10583 s for one loop in 100 iterations.
• could this be reduced to a list comprehension? Something like [[calculate_distance(sequences[i], sequences[j]) for i in range(N)] for j in range(i+1, N)]? That's probably not correct, but I feel like it should be possible. Dec 6, 2016 at 16:04