Counting Optimal Alignments

https://rosalind.info/problems/ctea/

ROSALIND | Counting Optimal Alignments

Problem

Recall from “Edit Distance Alignment” that if s′s′ and t′t′ are the augmented strings corresponding to an alignment of strings ss and tt, then the edit alignment score of s′s′ and t′t′ was given by the Hamming distance dH(s′,t′)dH(s′,t′) (because s′s′ and t′t′ have the same length and already include gap symbols to denote insertions/deletions).

As a result, we obtain 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. Strings s′s′ and t′t′ achieving this minimum correspond to an optimal alignment with respect to edit alignment score.

Given: Two protein strings ss and tt in FASTA format, each of length at most 1000 aa.

Return: The total number of optimal alignments of ss and tt with respect to edit alignment score, modulo 134,217,727 (227-1).

Sample Dataset

>Rosalind_78
PLEASANTLY
>Rosalind_33
MEANLY

Sample Output

4

 

 이 문제는 edit distance가 가장 작은 optimal alignment의 개수를 구하는 문제이다. 이 문제는 optimal alignment를 구하는 코드에서 조금만 더 변형해서 푼다면 쉽게 풀 수 있다. 

 우리가 optimal alignment를 구할 때에는 배열에서 [i,j]의 값을 구하기 위해서 [i-1,j],[i,j-1],[i-1,j-1]의 값을 비교하였다. 지금 비교하는 문자열이 같다면 i-1,j-1의 값을 그대로 아니라면 저 3개 값 중 최솟값을 취하는 식으로 구했다. 하지만 제거나 치환을 해도 그대로 두는 경우랑 edit distance가 같을 수 있는 경우가 있기에, 같다면 i-1,j-1 값 그대로 아니라면 +1 하는식으로 한다음 그 중 최솟값을 구했다. 그리고 optimal alignment의 개수는 일종의 여기까지 올 수 있는 배열에서의 경로를 이야기한다. 예를 들어 2x2의 정사각형에서 제2사분면에서 제4사분면으로 향하는 경우의 수는 3개이다. 이는 2->1->4,           2->3->4, 2->4  같은 식이다. 이는 맞닿아 있는 3개의 배열에서의 경로의 수를 합한 것이다. 

 그러니 우리는 최솟값을 가지고 있는 배열의 경로의 수를 계속해서 더해가는 식으로 optimal alignment의 개수를 구할 것이다.

from Bio import SeqIO

def edit_distance(s1,s2):
    array=[[0 for _ in range(len(s2))] for _ in range(len(s1))]#[s1][s2]로 넣기
    dp=[[0 for _ in range(len(s2))] for _ in range(len(s1))]
    modulus = 2**27 - 1

    for i in range(len(s1)):
        for j in range(len(s2)):
            if i==0:#i가 0이면
                array[i][j] = j
                dp[i][j] = 1
            elif j==0:#j가 0이라
                array[i][j] = i
                dp[i][j] =1
            else:#그외
                temp=[array[i-1][j]+1,array[i][j-1]+1,array[i-1][j-1]+(s1[i]!=s2[j])]#맞닿이있는 3개의 배열     
                array[i][j] = min(temp)#edit distance = 최솟값

                dp[i][j] += [0,dp[i-1][j]][temp[0]==array[i][j]]#다르면0, 같다면 경로의 수 더하기
                dp[i][j] += [0,dp[i][j-1]][temp[1]==array[i][j]]
                dp[i][j] += [0,dp[i-1][j-1]][temp[2]==array[i][j]]

                dp[i][j] = dp[i][j] % modulus#지속적으로 나머지로

    return dp[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])

print(edit_distance(S1,S2))