[의생명정보알고리즘] PA3 과제

---

✅ PA3.py 파일 생성 및 수정 과정 요약

1. 교수님 코드 복사 (두 가지 방법 중 하나 선택)

방법 1: PA2.py 파일 기반으로 복사→ 교수님이 제공한 PA2.py 파일을 현재 작업 디렉토리에 PA3.py라는 이름으로 복사

cp /mss/home_student/ABD2025/PA2.py PA3.py

 

방법 2: PA3.tmp.py 파일 기반으로 복사→ 교수님이 제공한 PA3.tmp.py 파일을 PA3.py라는 이름으로 복사

cp /mss/home_student/ABD2025/PA3.tmp.py PA3.py

1. PA3.py 파일 수정 시작
   • 복사된 PA3.py 파일을 열고 필요한 내용을 수정하면서 쌍서열 정렬(또는 Affine Gap Penalty 등)에 맞게 구현을 시작하면 됨.

---

교수님의 PA3.tmp.py 코드

#!/usr/bin/env python3

import sys
import random

match = int(sys.argv[1])
mismatch = int(sys.argv[2])
# gap = int(sys.argv[3])
# Variables for gap open and extension scores

# Change the index of the following two lines
seq1 = sys.argv[4]
seq2 = sys.argv[5]

#print(match, mismatch, gap, seq1, seq2)

# Creating DP matrix
# Change F to M
# Need to add two more matrices, IX and IY
M = {}
IX = ...
IY = ...
for i in range(len(seq1) + 1):
    M[i] = {}
    IX...
    IY...

# Initialization for M, IX, IY
M[0][0] = [0, '']
IX...
IY...

# Column initialization for M, IX, IY
for i in range(1, len(seq1) + 1):
    M[i][0] = [float('-inf'), '']
    IX...
    IY...
    IY...

# Row initialization for M, IX, IY
for j in range(1, len(seq2) + 1):
    M[0][j] = [float('-inf'), '']
    IX...
    IX...
    IY...

#print(F)
for i in range(1, len(seq1) + 1):
    nt1 = seq1[i-1]
    for j in range(1, len(seq2) + 1):
        nt2 = seq2[j-1]

        ## For M
        # Code for finding the max score for M
        # Direction for the max score

        M[i,j] = [..., ...]

        ## For IX
        # Code for finding the max score for IX
        # Direction for the max score

        IX[i,j] = [..., ...]

        ## For IY
        # Code for finding the max score for IY

        # Direction for the max score

        IY[i,j] = [..., ...]

# Finding the max score from M, IX, IY
maxscore = max(....)
print("Optimal alignment score =", maxscore)

# Track back
aseq1 = ''
aseq2 = ''

i = len(seq1)   # m
j = len(seq2)   # n

# Code for finding the starting matrix among M, IX, IY
# Need to perform random selection
cdirs = ''
if maxscore == M[...][...][0]:
    cdirs += 'm'
if maxscore == IX[...][...][0]:
    cdirs += 'x'
if maxscore == IY[...][...][0]:
    cdirs += 'y'

cmat = random.choice(cdirs)

while i != 0 or j != 0:
    # Need to find the current location (M, IX, or IY)
    # Need to get the directions from the current matrix

    dirs = ...
    next_dir = random.choice(dirs) # next_dir = 'u' or 'l' or 'd'

    # Code for handling directions: M(d, IX, IY), IX(M, u), IY(M, l)

print(aseq1)
print(aseq2)

---

 

🧬 전체 구조 요약

Step 1️⃣	사용자 입력 받기
Step 2️⃣	DP 테이블 초기화 (M, IX, IY)
Step 3️⃣	DP 테이블 채우기 (각 위치의 최적 점수 계산)
Step 4️⃣	최종 점수와 traceback 시작점 선택 (랜덤)
Step 5️⃣	traceback 실행 (정렬 결과 중 무작위 하나만 출력)

---

🔍 Step-by-Step 설명

---

✅ Step 1️⃣: 입력 받기

match = int(sys.argv[1])      # 일치 시 점수
mismatch = int(sys.argv[2])   # 불일치 시 점수
gap_open = int(sys.argv[3])   # 갭 시작 벌점
gap_extend = int(sys.argv[4]) # 갭 연장 벌점
seq1 = sys.argv[5]            # 첫 번째 서열
seq2 = sys.argv[6]            # 두 번째 서열

📝 사용자로부터 match/mismatch/gap 점수와 두 개의 서열을 입력 받는다.

---

✅ Step 2️⃣: DP 테이블 초기화

M = {}   # 매치 or 미스매치 상태 점수
IX = {}  # seq1에서 갭 생긴 상태
IY = {}  # seq2에서 갭 생긴 상태

🔸 M, IX, IY 각각은 DP 테이블.
🔸 각 칸은 [점수, 이전 방향 리스트] 형태로 저장돼 ([score, ['m', 'x', 'y']]).

그리고 초기값은 아래처럼 설정:

M[0][0] = [0, []]
IX[0][0] = [gap_open, ['m']]
IY[0][0] = [gap_open, ['m']]

🧠 Affine gap이기 때문에 갭이 시작되면 gap_open, 이어질수록 gap_extend를 적용해!

---

✅ Step 3️⃣: DP 테이블 채우기

for i in range(1, len(seq1)+1):
    for j in range(1, len(seq2)+1):

각 위치 (i, j)에 대해:

• M[i][j]: 현재 문자가 매치/미스매치 될 경우
• IX[i][j]: seq1[i]에서 갭이 생길 경우 (위에서 내려옴)
• IY[i][j]: seq2[j]에서 갭이 생길 경우 (왼쪽에서 옴)

M[i][j] = max([
    (M[i-1][j-1][0] + s, 'm'),
    (IX[i-1][j-1][0] + s, 'x'),
    (IY[i-1][j-1][0] + s, 'y')
])

💡 이전 상태들로부터 가능한 점수를 모두 고려하고, 가장 큰 값을 저장해.

👉 IX, IY도 같은 방식으로 이전 상태를 참고하며 점수를 채워줘.

---

✅ Step 4️⃣: 최종 점수 & 시작 위치 선택 (랜덤)

i, j = len(seq1), len(seq2)
maxscore = max(M[i][j][0], IX[i][j][0], IY[i][j][0])

📌 마지막 위치에서 3개의 상태 중 최적 점수를 찾고,
그 중 랜덤으로 하나의 시작점(matrix) 을 고른다.

starts = [('m', i, j), ...]
cmat, i, j = random.choice(starts)

---

✅ Step 5️⃣: Traceback (정렬 경로 무작위로 1개만)

while i > 0 or j > 0:

각 상태(M, IX, IY)에서 이전 상태들 중 하나를 랜덤으로 선택하면서
서열을 거꾸로 따라가며 하나의 정렬을 만든다.

• M: 대각선으로 이동하며 문자 두 개 정렬
• IX: 위로 이동하며 seq1 문자와 - 정렬
• IY: 왼쪽으로 이동하며 -와 seq2 문자 정렬

마지막에 거꾸로 만든 문자열을 출력!

print(aseq1)
print(aseq2)

---

✅ 출력 예시

$ ./PA3.py 2 1 -2 -1 ATCTG TGCA
Optimal alignment score = 3
ATCTG
-TGCA

(또 실행하면 정렬 결과가 다를 수 있음)

---

🧬 최종 코드

#!/usr/bin/env python3

import sys
import random

# 1️⃣ 사용자 입력값 받기 (명령줄 인자)
match = int(sys.argv[1])         # 일치할 때 점수
mismatch = int(sys.argv[2])      # 불일치할 때 감점
gap_open = int(sys.argv[3])      # 갭 시작 시 감점
gap_extend = int(sys.argv[4])    # 갭 이어질 때 감점
seq1 = sys.argv[5]               # 첫 번째 염기서열
seq2 = sys.argv[6]               # 두 번째 염기서열

# 2️⃣ DP 테이블(M, IX, IY) 초기화
# M: 매치/미스매치 상태
# IX: seq1에서 gap 생긴 상태
# IY: seq2에서 gap 생긴 상태
M, IX, IY = {}, {}, {}

for i in range(len(seq1)+1):
    M[i], IX[i], IY[i] = {}, {}, {}
    for j in range(len(seq2)+1):
        # 각 칸은 [점수, 이전 방향 리스트] 형태로 저장됨
        M[i][j] = [-float('inf'), []]
        IX[i][j] = [-float('inf'), []]
        IY[i][j] = [-float('inf'), []]

# 시작 위치 (0,0) 초기화
M[0][0] = [0, []]
IX[0][0] = [gap_open, ['m']]   # M에서 시작해 IX로 간 경우
IY[0][0] = [gap_open, ['m']]   # M에서 시작해 IY로 간 경우

# 3️⃣ 첫 번째 열 초기화 (seq2가 비었을 때)
for i in range(1, len(seq1)+1):
    IX[i][0] = [gap_open + (i-1)*gap_extend, ['x']]
    IY[i][0] = [-float('inf'), []]  # 왼쪽으로 못 옴

# 첫 번째 행 초기화 (seq1이 비었을 때)
for j in range(1, len(seq2)+1):
    IY[0][j] = [gap_open + (j-1)*gap_extend, ['y']]
    IX[0][j] = [-float('inf'), []]  # 위로 못 옴

# 4️⃣ DP 테이블 채우기 (각 칸마다 최적 점수 계산)
for i in range(1, len(seq1)+1):
    nt1 = seq1[i-1]  # 현재 seq1 문자
    for j in range(1, len(seq2)+1):
        nt2 = seq2[j-1]  # 현재 seq2 문자
        s = match if nt1 == nt2 else mismatch

        # M[i][j]: 대각선에서 매치/미스매치 점수 계산
        opts = [
            (M[i-1][j-1][0] + s, 'm'),
            (IX[i-1][j-1][0] + s, 'x'),
            (IY[i-1][j-1][0] + s, 'y')
        ]
        best = max(x[0] for x in opts)
        M[i][j] = [best, [d for sc, d in opts if sc == best]]

        # IX[i][j]: 위쪽에서 gap 열기/연장
        opts = [
            (M[i-1][j][0] + gap_open, 'm'),
            (IX[i-1][j][0] + gap_extend, 'x')
        ]
        best = max(x[0] for x in opts)
        IX[i][j] = [best, [d for sc, d in opts if sc == best]]

        # IY[i][j]: 왼쪽에서 gap 열기/연장
        opts = [
            (M[i][j-1][0] + gap_open, 'm'),
            (IY[i][j-1][0] + gap_extend, 'y')
        ]
        best = max(x[0] for x in opts)
        IY[i][j] = [best, [d for sc, d in opts if sc == best]]

# 5️⃣ 마지막 칸에서 최적 점수 및 시작 matrix 선택
i, j = len(seq1), len(seq2)
maxscore = max(M[i][j][0], IX[i][j][0], IY[i][j][0])

# 최종 점수 출력
print("Optimal alignment score =", maxscore)

# 최적 점수인 matrix 후보들 중 하나를 무작위로 선택
starts = []
if M[i][j][0] == maxscore:
    starts.append(('m', i, j))
if IX[i][j][0] == maxscore:
    starts.append(('x', i, j))
if IY[i][j][0] == maxscore:
    starts.append(('y', i, j))

cmat, i, j = random.choice(starts)  # (matrix, i, j)

# 6️⃣ Traceback: 하나의 경로만 무작위로 따라가며 정렬 복원
aseq1 = ''
aseq2 = ''

while i > 0 or j > 0:
    if cmat == 'm':
        prev_dirs = M[i][j][1]
        cmat = random.choice(prev_dirs)
        i -= 1
        j -= 1
        aseq1 = seq1[i] + aseq1
        aseq2 = seq2[j] + aseq2
    elif cmat == 'x':
        prev_dirs = IX[i][j][1]
        cmat = random.choice(prev_dirs)
        i -= 1
        aseq1 = seq1[i] + aseq1
        aseq2 = '-' + aseq2
    elif cmat == 'y':
        prev_dirs = IY[i][j][1]
        cmat = random.choice(prev_dirs)
        j -= 1
        aseq1 = '-' + aseq1
        aseq2 = seq2[j] + aseq2

# 7️⃣ 결과 출력
print(aseq1)
print(aseq2)

---

🎯 예시 실행

$ ./PA3.py 2 1 -2 -1 ATCTG TGCA
Optimal alignment score = 3
ATCTG
-TGCA

---

🧠 핵심 요약

• Affine Gap: 갭을 "시작할 때"와 "계속될 때" 점수를 다르게 줌
• M / IX / IY: 세 가지 상태별 점수를 따로 관리
• Traceback: 여러 최적 경로 중 하나를 무작위로 따라가 정렬 생성
• 정렬 1개만 출력: 성능, 간결성 모두 확보

---

=> 이렇게 3가지 나오면 과제 끝!