Biopython for similarity matrix - 尋找更好的性能

我想計算輸入序列和序列中短片段之間的相似度。結果是一個相似度矩陣，每個位置都是對齊的分數。它可以工作，但不幸的是很慢。我如何更有效地在python和numpy中實現循環？我也在考慮使用MPI，但多線程或甚至更好的內部numpy解決方案將是有益的。以下是代碼。Biopython for similarity matrix - 尋找更好的性能

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)

來源

2014-09-30 El Dude

你真正在做的是本地對齊。爲此，您需要Smith-Waterman算法。 – wasserfeder 2014-10-01 03:51:19

我想幫忙，但我不明白你想達到什麼目的。無論如何，你的代碼可以使用pypy獲益很多。試試吧，你不需要改變任何東西。如果您需要更多幫助，請使用示例更新您的問題。 – tbrittoborges 2014-10-01 10:43:55

@wasserfeder。是的，這是本地調整，我對矩陣感興趣，而不是調整結果。但在Biopython中，我沒有返回矩陣的函數 - 因此我想自己生成它... – 2014-10-01 16:31:44

首先我想嘗試只計算對角線的一半，開始在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

另一件事可以做只計算得分，但我覺得這是太小的改進。從numpy

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

矢量化：做分析。 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循環。

來源

2014-10-01 20:22:31 xbello

是的，在發佈原始問題後，我將代碼更改爲nxn/2。我會看看globalds的功能，謝謝你的提示。 – 2014-10-02 18:19:30

我要添加你要求的'numpy'，但除了「美麗比醜陋好」之外沒什麼用處。同時保存一些for循環。 – xbello 2014-10-02 22:14:09

Biopython for similarity matrix - 尋找更好的性能

回答

相關問題