Edit Distance Alignment

https://rosalind.info/problems/edta/

ROSALIND | Edit Distance Alignment

Problem

An alignment of two strings ss and tt is defined by two strings s′s′ and t′t′ satisfying the following three conditions: 1. s′s′ and t′t′ must be formed from adding gap symbols "-" to each of ss and tt, respectively; as a result, ss and tt will form subsequences of s′s′ and t′t′. 2. s′s′ and t′t′ must have the same length. 3. Two gap symbols may not be aligned; that is, if s′[j]s′[j] is a gap symbol, then t′[j]t′[j] cannot be a gap symbol, and vice-versa.

We say that s′s′ and t′t′ augment ss and tt. Writing s′s′ directly over t′t′ so that symbols are aligned provides us with a scenario for transforming ss into tt. Mismatched symbols from ss and tt correspond to symbol substitutions; a gap symbol s′[j]s′[j] aligned with a non-gap symbol t′[j]t′[j] implies the insertion of this symbol into tt; a gap symbol t′[j]t′[j] aligned with a non-gap symbol s′[j]s′[j] implies the deletion of this symbol from ss.

Thus, an alignment represents a transformation of ss into tt via edit operations. We define the corresponding edit alignment score of s′s′ and t′t′ as dH(s′,t′)dH(s′,t′) (Hamming distance is used because the gap symbol has been introduced for insertions and deletions). It follows that dE(s,t)=mins′,t′dH(s′,t′)dE(s,t)=mins′,t′dH(s′,t′), where the minimum is taken over all alignments of ss and tt. We call such a minimum score alignment an optimal alignment (with respect to edit distance).

Given: Two protein strings ss and tt in FASTA format (with each string having length at most 1000 aa).

Return: The edit distance dE(s,t)dE(s,t) followed by two augmented strings s′s′ and t′t′ representing an optimal alignment of ss and tt.

Sample Dataset

>Rosalind_43
PRETTY
>Rosalind_97
PRTTEIN

Sample Output

4
PRETTY--
PR-TTEIN

 

 이 문제는 주어진 두 문자열의 edit distance와 optimal alignment를 출력하는 문제이다. edit distnace문제의 코드를 응용해서 풀면 쉽게 풀 수 있다. optimal alignment는 edit distance가 가장 적게 하는 alignment이다. alignment는 말 그대로 두 문자열을 정렬하는 것으로 두 문자열에 gap sequence '-'를 추가하여 같은 길이로 만들어 정렬하는 것이다. 

 이 문제는 edit distance처럼 2차원 배열에서 다이나믹 프로그래밍 기법을 사용하여 하나하나 훑으면서 푼다. 말로는 어려우니 코드의 주석을 참고하도록 해라.

from Bio import SeqIO

def edit_distance(s1,s2):
    array=[[[0,'',''] for _ in range(len(s2))] for _ in range(len(s1))]

    for i in range(len(s1)):
        for j in range(len(s2)):
            if i==0:
                array[i][j][0] = j
            elif j==0:
                array[i][j][0] = i
            elif s1[i] == s2[j]:
                array[i][j][0] = array[i-1][j-1][0]
                array[i][j][1] = array[i-1][j-1][1] + s1[i]
                array[i][j][2] = array[i-1][j-1][2] + s2[j]
            else:
                temp=min([array[i-1][j][0],array[i][j-1][0],array[i-1][j-1][0]])#이 중 최소
                array[i][j][0] = 1+temp
                if array[i-1][j][0] == temp:#1에 add
                    array[i][j][1] = array[i-1][j][1] + s1[i]
                    array[i][j][2] = array[i-1][j][2] + '-'
                elif array[i][j-1][0] == temp:#2에 add
                    array[i][j][1] = array[i][j-1][1] + '-'
                    array[i][j][2] = array[i][j-1][2] + s2[j]
                elif array[i-1][j-1][0] == temp:#치환이므로 각각
                    array[i][j][1] = array[i-1][j-1][1] + s1[i]
                    array[i][j][2] = array[i-1][j-1][2] + s2[j]



    return array[len(s1)-1][len(s2)-1]

if __name__ == '__main__':
    with open(r"파일경로",'r') as fa:
        seqs=[]
        for s in SeqIO.parse(fa,'fasta'):
            seqs.append(s.seq)
        S1=' '+str(seqs[0])
        S2=' '+str(seqs[1])

f=open(r"파일경로",'w')
for i in edit_distance(S1,S2):
    print(i,file=f)

이 중 유의해야 할 점은 치환하는 경우를 가장 마지막에 두어, gap을 넣는 것을 우선으로 두어 최대한 두 서열에 같은 부분이 겹치도록 해야한다.