使用一组字符而不是一个字符的序列对齐算法

Summary:

我从一些有关对齐算法的细节开始，最后我提出了我的问题。如果您了解对齐算法，请从头开始。

考虑我们有两个字符串，例如：

ACCGAATCGA
ACCGGTATTAAC

有一些算法，例如：史密斯-沃特曼 https://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm Or 尼德曼-温施 https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm，对齐这两个序列并创建一个矩阵。查看以下部分的结果：

Smith-Waterman Matrix
§   §   A   C   C   G   A   A   T   C   G   A   
§   0   0   0   0   0   0   0   0   0   0   0   
A   0   4   0   0   0   4   4   0   0   0   4   
C   0   0   13  9   4   0   4   3   9   4   0   
C   0   0   9   22  17  12  7   3   12  7   4   
G   0   0   4   17  28  23  18  13  8   18  13  
G   0   0   0   12  23  28  23  18  13  14  18  
T   0   0   0   7   18  23  28  28  23  18  14  
A   0   4   0   2   13  22  27  28  28  23  22  
T   0   0   3   0   8   17  22  32  27  26  23  
T   0   0   0   2   3   12  17  27  31  26  26  
A   0   4   0   0   2   7   16  22  27  31  30  
A   0   4   4   0   0   6   11  17  22  27  35  
C   0   0   13  13  8   3   6   12  26  22  30  

Optimal Alignments
A   C   C   G   A   -   A   T   C   G   A   
A   C   C   G   G   A   A   T   T   A   A   

问题：

我的问题很简单，但也许答案并不像看起来那么容易。我想将一组字符用作单个字符，例如：[A0][C0][A1][B1]。但在这些算法中，我们必须使用单个字符。我们怎样才能做到这一点？

附：考虑我们有这个序列：#read #write #add #write。然后我将其转换为类似的内容：#read to A .... #write to B.... #add to C。然后我的序列变为：ABCB。但我有很多不同的词开头#。而且 ASCII 表不足以全部转换。然后我需要更多的角色。唯一的方法是使用类似的东西[A0] ... [Z9]对于每个单词。或者使用数字。

P.S：Smith-Waterman 的一些示例代码存在于此link http://www.vivid-abstractions.net/logical/programming/smith-waterman-algorithm-program-and-source-code/

PS：有另一个帖子 https://stackoverflow.com/questions/29847642/sequence-alignment-with-one-to-many-characters他们想要类似的东西，但我想要的是不同的。在这个问题中，我们有一组以 a 开头的字符[并以]。并且不需要使用像这样的语义ee等于i.

我适应了这个Python实现 https://github.com/alevchuk/pairwise-alignment-in-pythonSmith-Waterman 和 Needleman-Wunsch 算法（GPL 版本 3 许可）支持具有多个字符组的序列：

#This software is a free software. Thus, it is licensed under GNU General Public License.
#Python implementation to Smith-Waterman Algorithm for Homework 1 of Bioinformatics class.
#Forrest Bao, Sept. 26 <http://fsbao.net> <forrest.bao aT gmail.com>

# zeros() was origianlly from NumPy.
# This version is implemented by alevchuk 2011-04-10
def zeros(shape):
    retval = []
    for x in range(shape[0]):
        retval.append([])
        for y in range(shape[1]):
            retval[-1].append(0)
    return retval

match_award      = 10
mismatch_penalty = -5
gap_penalty      = -5 # both for opening and extanding
gap = '----' # should be as long as your group of characters
space = '    ' # should be as long as your group of characters

def match_score(alpha, beta):
    if alpha == beta:
        return match_award
    elif alpha == gap or beta == gap:
        return gap_penalty
    else:
        return mismatch_penalty

def finalize(align1, align2):
    align1 = align1[::-1]    #reverse sequence 1
    align2 = align2[::-1]    #reverse sequence 2

    i,j = 0,0

    #calcuate identity, score and aligned sequeces
    symbol = []
    found = 0
    score = 0
    identity = 0
    for i in range(0,len(align1)):
        # if two AAs are the same, then output the letter
        if align1[i] == align2[i]:                
            symbol.append(align1[i])
            identity = identity + 1
            score += match_score(align1[i], align2[i])

        # if they are not identical and none of them is gap
        elif align1[i] != align2[i] and align1[i] != gap and align2[i] != gap:
            score += match_score(align1[i], align2[i])
            symbol.append(space)
            found = 0

        #if one of them is a gap, output a space
        elif align1[i] == gap or align2[i] == gap:
            symbol.append(space)
            score += gap_penalty

    identity = float(identity) / len(align1) * 100

    print 'Identity =', "%3.3f" % identity, 'percent'
    print 'Score =', score
    print ''.join(align1)
    # print ''.join(symbol)
    print ''.join(align2)


def needle(seq1, seq2):
    m, n = len(seq1), len(seq2)  # length of two sequences

    # Generate DP table and traceback path pointer matrix
    score = zeros((m+1, n+1))      # the DP table

    # Calculate DP table
    for i in range(0, m + 1):
        score[i][0] = gap_penalty * i
    for j in range(0, n + 1):
        score[0][j] = gap_penalty * j
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            match = score[i - 1][j - 1] + match_score(seq1[i-1], seq2[j-1])
            delete = score[i - 1][j] + gap_penalty
            insert = score[i][j - 1] + gap_penalty
            score[i][j] = max(match, delete, insert)

    # Traceback and compute the alignment 
    align1, align2 = [], []
    i,j = m,n # start from the bottom right cell
    while i > 0 and j > 0: # end toching the top or the left edge
        score_current = score[i][j]
        score_diagonal = score[i-1][j-1]
        score_up = score[i][j-1]
        score_left = score[i-1][j]

        if score_current == score_diagonal + match_score(seq1[i-1], seq2[j-1]):
            align1.append(seq1[i-1])
            align2.append(seq2[j-1])
            i -= 1
            j -= 1
        elif score_current == score_left + gap_penalty:
            align1.append(seq1[i-1])
            align2.append(gap)
            i -= 1
        elif score_current == score_up + gap_penalty:
            align1.append(gap)
            align2.append(seq2[j-1])
            j -= 1

    # Finish tracing up to the top left cell
    while i > 0:
        align1.append(seq1[i-1])
        align2.append(gap)
        i -= 1
    while j > 0:
        align1.append(gap)
        align2.append(seq2[j-1])
        j -= 1

    finalize(align1, align2)

def water(seq1, seq2):
    m, n = len(seq1), len(seq2)  # length of two sequences

    # Generate DP table and traceback path pointer matrix
    score = zeros((m+1, n+1))      # the DP table
    pointer = zeros((m+1, n+1))    # to store the traceback path

    max_score = 0        # initial maximum score in DP table
    # Calculate DP table and mark pointers
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            score_diagonal = score[i-1][j-1] + match_score(seq1[i-1], seq2[j-1])
            score_up = score[i][j-1] + gap_penalty
            score_left = score[i-1][j] + gap_penalty
            score[i][j] = max(0,score_left, score_up, score_diagonal)
            if score[i][j] == 0:
                pointer[i][j] = 0 # 0 means end of the path
            if score[i][j] == score_left:
                pointer[i][j] = 1 # 1 means trace up
            if score[i][j] == score_up:
                pointer[i][j] = 2 # 2 means trace left
            if score[i][j] == score_diagonal:
                pointer[i][j] = 3 # 3 means trace diagonal
            if score[i][j] >= max_score:
                max_i = i
                max_j = j
                max_score = score[i][j];

    align1, align2 = [], []    # initial sequences

    i,j = max_i,max_j    # indices of path starting point

    #traceback, follow pointers
    while pointer[i][j] != 0:
        if pointer[i][j] == 3:
            align1.append(seq1[i-1])
            align2.append(seq2[j-1])
            i -= 1
            j -= 1
        elif pointer[i][j] == 2:
            align1.append(gap)
            align2.append(seq2[j-1])
            j -= 1
        elif pointer[i][j] == 1:
            align1.append(seq1[i-1])
            align2.append(gap)
            i -= 1

    finalize(align1, align2)

如果我们使用以下输入运行它：

seq1 = ['[A0]', '[C0]', '[A1]', '[B1]']
seq2 = ['[A0]', '[A1]', '[B1]', '[C1]']

print "Needleman-Wunsch"
needle(seq1, seq2)
print
print "Smith-Waterman"
water(seq1, seq2)

我们得到这个输出：

Needleman-Wunsch
Identity = 60.000 percent
Score = 20
[A0][C0][A1][B1]----
[A0]----[A1][B1][C1]

Smith-Waterman
Identity = 75.000 percent
Score = 25
[A0][C0][A1][B1]
[A0]----[A1][B1]

我所做的具体修改参见：这个 GitHub 存储库 https://github.com/BioGeek/pairwise-alignment-in-python/commits/master/alignment.py.