Biopython для матрицы подобия - ищет лучшую производительность

Question

Я хочу рассчитать подобие между входной последовательностью и коротким фрагментом из последовательности. Результатом является матрица подобия с каждой позицией, являющейся оценкой выравнивания. Работает, но, к сожалению, медленно. Как я мог реализовать цикл более эффективно в python и numpy? Я также думаю использовать MPI, но многопоточность или даже лучшее внутреннее решение numpy было бы полезно. Ниже приведен код.

from Bio import pairwise2 
import numpy 

.... 

similarityMatrix = numpy.zeros(shape=(sequenceLength-fragmentLength,sequenceLength-fragmentLength)) 

for i in xrange(sequenceLength-fragmentLength): 
    currentFragment = sequence[i:i+fragmentLength] 

    for j in xrange(i,sequenceLength-fragmentLength): 
     aFragment = sequence[j:j+fragmentLength] 

     alns = pairwise2.align.globalds(aFragment, currentFragment, matrix, gap_open, gap_extend) 

     bestHit = alns[0] 
     score = bestHit[2] 

     similarityMatrix[i,j] = float(score) 
     similarityMatrix[j,i] = float(score) 

Answer 1

Сначала я попытался бы вычислить только половину диагонали, начиная внутреннего контура в i точке и избежать вычисления предыдущих рядов:

for i in range(full_size - frag_size): 
    curr_frag = seq[i:i + frag_size] 

    # ADD THIS ----vvvvv------------vvv 
    for j in range(i + 1, full_size + 1 - frag_size): 
     match_frag = seq[j:j + frag_size] 
     # Do the following calculation here 

Другое дело, может быть, чтобы сделать только вычисление оценки, но я считаю, что это слишком незначительное улучшение. Выполните профилирование:

score = pairwise2.align.globalds(aFragment, 
           currentFragment, 
           matrix, 
           gap_open, 
           gap_extend, 
           score_only=True) # <= ADD THIS 

---

Vectorize из numpy. http://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html

Определить векторизованную функцию, которая принимает последовательность вложенных объектов или Numpy массивов в качестве входных данных и возвращает массив Numpy в качестве вывода.

def chop(sequence, frag_size): 
    for i in range(full_size - frag_size): 
     yield sequence[i + 1:i + 1 + frag_size] 

def pairwise(seq1, seq2): 
    return pairwise2.align.globalds(
     seq1, seq2, MATRIX, -2, -1, score_only=True) 

query = numpy.array([x for x in chop(seq, frag_size)]) 
subject = numpy.array([x for x in chop(seq, frag_size)]) 

vfunc = numpy.vectorize(pairwise) 

results = [] 
for i in subject: 
    results.append(vfunc(i, query)) 

Но вы не получите не производительность. Как говорится:

Функция векторизации предоставляется в первую очередь для удобства, а не для производительности. Реализация по существу является циклом for.