diff --git a/README.md b/README.md index 82a361d..937068a 100644 --- a/README.md +++ b/README.md @@ -31,30 +31,8 @@ Terdapat dua titik capai dalam tugas ini : [3] Anda dapat membuat 2 buah algoritma untuk perbanding antar sekuens dan perbandingan antar profil. Coba anda ubah agar anda melakukan keduanya secara langsung, sehingga mengurangi beban kerja anda.
[4] Anda bisa saja menggunakan Python untuk mengerjakan tugas, namun ingat bahwa **kinerja C++ dan C jauh lebih cepat**. Sebagai pengalaman, asisten menjalankan skoring pensejajaran global 2 DNA dengan panjang ~29000 nukleotida. Algoritma berjalan 100 menit untuk Python, dan algoritma berjalan hanya 3 menit untuk bahasa C++ dengan flag -O3 (optimization) ketika kompilasi (**33 x speedup !**). Sebagai saran (bila anda keukeuh menggunakan Python) , anda bisa menggunakan Python optimizer (misal Numba ataupun Cython) untuk mempercepat eksekusi algoritma anda. -## Pengumpulan -### Pengerjaan -Silahkan lakukan *fork* dari *repository* ini. - -### Deliverables -1. File yang berisi hasil pensejajaran sekuens global. Misal , hasil pensejajaran antara sekuens pertama (file1.fasta) dan sekuens kedua (file2.fasta) ditulis dalam folder /result/file1_file2/. Lalu dalam folder file1_file2, tuliksan hasil pensejajaran masing-masing file1.fasta dan file2.fasta sebagai file1.txt dan file2.txt. Tuliskan score dalam sebuah file score.txt. -2. Kode sumber. Tuliskan cara kompilasi bila menggunakan *compiled language*, namun lebih baik dilengkapi dengan Makefile. -3. Ubah Readme ini. Tuliskan pendekatan pensejajaran yang anda lakukan dan cara menjalankan program. - -### Teknis Pengumpulan -      Kumpulkan dengan membuat *merge request* pada *repository* ini. Batas pengumpulan dan demo adalah 29 Juli 2020. - -### Demo -      Setelah selesai, jadwalkan demo dengan asisten. Kontak dapat dilihat pada Readme ini. Demo berlangsung 15-30 menit. Demo akan berisi tanya jawab, namun belum tentu akan diisi oleh pengujian, tergantung *runtime* dari algoritma anda. - -## Penilaian -      Uji cobakan algoritma anda dengan data yang telah disediakan di *repository* ini. Percobaan minimal menghasilkan 3 pensejajaran global untuk kasus 2 sekuens, 2 pensejajaran global kasus 3 sekuens, 1 pensejajaran global untuk kasus 4 sekuens.
-      Nilai maksimal 1350 untuk *milestone* pertama, dan nilai maksimal 1800 untuk milestone kedua. Nilai maksimal demo adalah 850, sehingga nilai maksimal total adalah 4000. Algoritma anda **wajib optimum** untuk pensejajaran global 2 sekuens. Akan tetapi, algoritma anda tidak harus optimum untuk MSA. Implementasi MSA lebih kepada *proof-of-concept*. - -## Kontak -      Silahkan hubungi asisten lewat line @alamhasabiebaru atau lewat email 13517096@std.stei.itb.ac.id dengan subjek \[SELEKSI IRK - SEQUENCE ALIGNMENT\] . *Note : waktu menjawab bervariasi, namun email biasanya akan dibalas kurang dari sehari. Line mungkin tidak dibalas dalam waktu satu-dua hari. Mohon bersabar :)*. Pertanyaan juga dipersilahkan. Jawaban akan diposting dalam bagian QnA README ini. - -## QnA -null. +## Run +python main.py ## Referensi Silahkan gunakan referensi berikut sebagai awal pengerjaan tugas:
diff --git a/main.py b/main.py new file mode 100644 index 0000000..8beff07 --- /dev/null +++ b/main.py @@ -0,0 +1,110 @@ +def read_fasta(fp): + name, seq = None, [] + for line in fp: + line = line.rstrip() + if line.startswith(">"): + if name: yield (name, ''.join(seq)) + name, seq = line, [] + else: + seq.append(line) + if name: yield (name, ''.join(seq)) + +def given_matrices_inserter(filename): + try: + data = open(filename, 'r') + dimensions = data.readline() + try: + n = int(dimensions) + except ValueError: + print("Wrong file") + exit() + letters = data.readline() + letters = letters.replace('\n', '') + letters_arr = letters.split(' ') + score_matrix = np.zeros((n,n)) + for i in range(0,n): + arr = data.readline().split(" ") + for j in range(0,n): + score_matrix[i][j] = float(arr[j]) + return letters_arr, score_matrix + except FileNotFoundError: + print("There is no matrix file.") + exit() + +import numpy as np + +def write_to_file(filename, value): + with open(filename, "w+") as f: + f.write(value) + +def match_score(c1, c2, m, mm, letters, score_matrix): + mapped_values = {} + for i, v in enumerate(letters): + mapped_values[v] = i + a = mapped_values[c1] + b = mapped_values[c2] + return score_matrix[a][b] + +def global_alignment(first, second, letters, matrix, protein): + match, mismatch, gap, ind = 1,-1,-1,-1 + if ind == 0 and protein == 0: + letters, matrix = matrix_chooser(letters, match, mismatch) + n = len(first) + m = len(second) + backtrack = [[(-1,1) for j in range(m+1)] for i in range(n+1)] + s = [[0 for j in range(m+1)] for i in range(n+1)] + + for i in range(1, n+1): + s[i][0] = s[i-1][0] + gap + backtrack[i][0] = (i-1, 0) + for j in range(1, m+1): + s[0][j] = s[0][j-1] + gap + backtrack[0][j] = (0, j-1) + + for i in range(1, n+1): + for j in range(1, m+1): + s[i][j] = max(s[i-1][j] + gap, s[i][j-1] + gap, s[i-1][j-1] + match_score(first[i-1], second[j-1], match, mismatch, letters, matrix)) + if s[i][j] == s[i-1][j] + gap: + backtrack[i][j] = (i-1, j) + elif s[i][j] == s[i][j-1] + gap: + backtrack[i][j] = (i, j-1) + else: + backtrack[i][j] = (i-1, j-1) + + first_p = "" + second_p = "" + i = n + j = m + while (i,j) != (0,0): + if backtrack[i][j] == (i-1, j-1): + first_p = first[i-1] + first_p + second_p = second[j-1] + second_p + elif backtrack[i][j] == (i-1, j): + first_p = first[i-1] + first_p + second_p = '-' + second_p + else: + first_p = '-' + first_p + second_p = second[j-1] + second_p + (i, j) = backtrack[i][j] + + return first_p, second_p, s[n][m] + +file1 = 'AGP04929.1.fasta' +file2 = 'AYN64561.1.fasta' + +with open(f'./data/protein/{file1}') as fp: + A = '' + for name, seq in read_fasta(fp): + A += seq + +with open(f'./data/protein/{file2}') as fp: + B = '' + for name, seq in read_fasta(fp): + B += seq + +protein = 1 +letters, matrix = given_matrices_inserter('pam250.txt') +resA, resB, score = global_alignment(A, B, letters, matrix, protein) +write_to_file(f'./outputs/{file1}.txt', resA) +write_to_file(f'./outputs/{file2}.txt', resB) +print(f'Score: {score}') diff --git a/pam250.txt b/pam250.txt new file mode 100644 index 0000000..1a4e1cb --- /dev/null +++ b/pam250.txt @@ -0,0 +1,26 @@ +23 +G A V L I P S T D E N Q K R H F Y W M C B Z X +5 1 -1 -4 -3 0 1 0 1 0 0 -1 -2 -3 -2 -5 -5 -7 -3 -3 0 0 -1 -8 +1 2 0 -2 -1 1 1 1 0 0 0 0 -1 -2 -1 -3 -3 -6 -1 -2 0 0 0 -8 +-1 0 4 2 4 -1 -1 0 -2 -2 -2 -2 -2 -2 -2 -1 -2 -6 2 -2 -2 -2 -1 -8 +-4 -2 2 6 2 -3 -3 -2 -4 -3 -3 -2 -3 -3 -2 2 -1 -2 4 -6 -3 -3 -1 -8 +-3 -1 4 2 5 -2 -1 0 -2 -2 -2 -2 -2 -2 -2 1 -1 -5 2 -2 -2 -2 -1 -8 +0 1 -1 -3 -2 6 1 0 -1 -1 0 0 -1 0 0 -5 -5 -6 -2 -3 -1 0 -1 -8 +1 1 -1 -3 -1 1 2 1 0 0 1 -1 0 0 -1 -3 -3 -2 -2 0 0 0 0 -8 +0 1 0 -2 0 0 1 3 0 0 0 -1 0 -1 -1 -3 -3 -5 -1 -2 0 -1 0 -8 +1 0 -2 -4 -2 -1 0 0 4 3 2 2 0 -1 1 -6 -4 -7 -3 -5 3 3 -1 -8 +0 0 -2 -3 -2 -1 0 0 3 4 1 2 0 -1 1 -5 -4 -7 -2 -5 3 3 -1 -8 +0 0 -2 -3 -2 0 1 0 2 1 2 1 1 0 2 -3 -2 -4 -2 -4 2 1 0 -8 +-1 0 -2 -2 -2 0 -1 -1 2 2 1 4 1 1 3 -5 -4 -5 -1 -5 1 3 -1 -8 +-2 -1 -2 -3 -2 -1 0 0 0 0 1 1 5 3 0 -5 -4 -3 0 -5 1 0 -1 -8 +-3 -2 -2 -3 -2 0 0 -1 -1 -1 0 1 3 6 2 -4 -4 -2 0 -4 -1 0 -1 -8 +-2 -1 -2 -2 -2 0 -1 -1 1 1 2 3 0 2 6 -2 0 -3 -2 -3 1 2 -1 -8 +-5 -3 -1 2 1 -5 -3 -3 -6 -5 -3 -5 -5 -4 -2 0 7 0 0 -4 -4 -5 -2 -8 +-5 -3 -2 -1 -1 -5 -3 -3 -4 -4 -2 -4 -4 -4 0 7 10 0 -2 0 -3 -4 -2 -8 +-7 -6 -6 -2 -5 -6 -2 -5 -7 -7 -4 -5 -3 -2 -3 0 0 17 -4 8 -5 -6 -4 -8 +-3 -1 2 4 2 -2 -2 -1 -3 -2 -2 -1 0 0 -2 0 -2 -4 6 -5 -2 -2 -1 -8 +-3 -2 -2 -6 -2 -3 0 -2 -5 -5 -4 -5 -5 -4 -3 -4 0 8 -5 12 -4 -5 -3 -8 +0 0 -2 -3 -2 -1 0 0 3 3 2 1 1 -1 1 -4 -3 -5 -2 -4 3 2 -1 -8 +0 0 -2 -3 -2 0 0 -1 3 3 1 3 0 0 2 -5 -4 -6 -2 -5 2 3 -1 -8 +-1 0 -1 -1 -1 -1 0 0 -1 -1 0 -1 -1 -1 -1 -2 -2 -4 -1 -3 -1 -1 -1 -8 +-8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 -8 1 \ No newline at end of file