Trình tự tương tự Python

*Lưu ý, nếu bạn muốn bỏ qua tính toán nền/căn chỉnh và đi thẳng đến nơi mã bắt đầu, chỉ cần nhấp vào đây

Lập trình động và DNA

Lập trình động có nhiều ứng dụng, bao gồm xác định sự giống nhau giữa hai chuỗi DNA hoặc RNA khác nhau, liên kết protein và trong nhiều ứng dụng khác trong tin sinh học [ngoài nhiều lĩnh vực khác]. Đối với bất kỳ ai ít quen thuộc hơn, lập trình động là một mô hình mã hóa giải quyết các vấn đề đệ quy bằng cách chia nhỏ chúng thành các vấn đề con bằng cách sử dụng một số loại cấu trúc dữ liệu để lưu trữ các kết quả của vấn đề con. Theo cách này, các bài toán đệ quy [chẳng hạn như dãy Fibonacci] có thể được lập trình hiệu quả hơn nhiều vì lập trình động cho phép bạn tránh các tính toán trùng lặp [và do đó, lãng phí] trong mã của bạn. Bấm vào đây để đọc thêm về lập trình động

Hãy đi sâu vào chủ đề chính của bài đăng này bằng cách triển khai thuật toán để đo lường sự tương đồng giữa hai chuỗi DNA

Động lực

Việc đo lường mức độ giống nhau giữa hai trình tự DNA khác nhau rất hữu ích vì nó có thể giúp chúng ta biết mức độ liên quan chặt chẽ [hoặc không] của các trình tự DNA đó và nguồn gốc của chúng [e. g. so sánh DNA của hai loài khác nhau, hoặc hai gen khác nhau]. Điều này giúp hiểu rõ hơn về cấu trúc và chức năng của một gen hoặc trình tự DNA mới được xác định

Tiểu sử

Đối với ví dụ về độ tương tự của chúng tôi, chúng tôi sẽ đo lường sự liên kết toàn cầu giữa hai chuỗi DNA và sau đó chỉ ra cách điều này có thể được thực hiện tổng quát hơn với lập trình động. “Toàn cầu” có nghĩa là chúng tôi sẽ sử dụng toàn bộ từng chuỗi DNA. “Căn chỉnh” có nghĩa là chúng tôi đang cố gắng sắp xếp các chuỗi DNA sao cho một số điểm được tối đa hóa dựa trên số lượng khớp/không khớp/khoảng trống mà chúng tôi có giữa các chuỗi. Bằng cách sắp xếp các chuỗi DNA theo cặp, chúng ta có thể tạo ra các khoảng trống giữa các nucleotide để tạo ra sự sắp xếp tốt hơn. Chúng ta sẽ hiểu thêm về điều đó trong giây lát. Có nhiều cách khác để đánh giá sự tương đồng của DNA và tôi có thể đề cập đến một số trong số chúng trong một bài viết trong tương lai, nhưng hiện tại chúng ta sẽ gắn bó với sự liên kết toàn cầu

DNA được tạo thành từ bốn loại nucleotide - adenine, guanine, cytosine và thymine. Chúng thường được viết tắt bằng các chữ cái A, G, C và T, tương ứng. Về mặt lập trình, DNA có thể được biểu diễn dưới dạng một chuỗi ký tự, trong đó mỗi ký tự phải là một trong A, G, C hoặc T. Sau đó, giả sử rằng chúng ta có hai chuỗi DNA như được thấy bên dưới

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Cách đo độ giống nhau giữa các chuỗi DNA

Chúng tôi cần một số liệu để sử dụng tính toán sự liên kết toàn cầu giữa các chuỗi DNA. Chúng tôi sẽ sử dụng phương pháp tính điểm [xem bên dưới] kết hợp với thuật toán Needlman-Wunsch để tìm ra sự liên kết toàn cầu tối ưu. Một ví dụ về việc thực hiện căn chỉnh toàn cầu không tối ưu trên các chuỗi của chúng tôi như sau

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Trong trường hợp trên, chúng tôi chỉ đơn thuần sắp xếp các chuỗi DNA theo cặp và làm nổi bật sự trùng khớp giữa các chuỗi. Bằng cách này, chúng tôi có hai kết quả phù hợp [được tô sáng] và bảy điểm không phù hợp

Hãy xác định điểm số mà chúng tôi đang cố gắng tối đa hóa như sau.

Cộng 1 cho mỗi điểm trùng khớp giữa các chuỗi
Trừ 1 cho mỗi điểm không khớp
Trừ 1 cho mỗi khoảng cách i. e. chèn/xóa [hiển thị bằng dấu “-“]

Số liệu tính điểm chính xác [thêm một điểm cho trận đấu, phạt một điểm vì không khớp / cách biệt] có thể khác nhau giữa các hệ thống tính điểm khác nhau, nhưng chúng tôi sẽ tuân thủ phương pháp tính điểm đã xác định của mình cho mục đích của bài đăng này

Theo thước đo này, sự liên kết toàn cầu của chúng tôi ở trên cho chúng tôi số điểm là 2 điểm phù hợp trừ 7 điểm không khớp = -5. Chúng tôi có thể cải thiện điểm số này nếu chúng tôi tính đến việc mỗi chuỗi cụ thể có thể có các phần chèn hoặc xóa. e. chúng tôi có thể thay đổi các nucleotide [hoặc ký tự] trong từng trình tự cụ thể để tạo thành một sự liên kết tốt hơn, miễn là chúng tôi không thực sự sắp xếp lại một trong hai trình tự. Nói cách khác, chúng ta có thể dịch chuyển các nuclêôtit theo trình tự miễn là thứ tự các nuclêôtit không thay đổi. Nếu điều đó không rõ ràng, ví dụ dưới đây cho thấy cách chúng tôi có thể đạt được điểm căn chỉnh toàn cầu tối đa bằng cách căn chỉnh đã sửa đổi

Sequence 1 ==> G T C C A T A - C A -
Sequence 2 ==> - T C - A T A T C A G

Các chữ cái được đánh dấu [i. e. nucleotide] cho thấy nơi có sự khớp giữa các chuỗi DNA. Ở đây chúng tôi có bảy trận đấu và bốn trận đấu không phù hợp

Tỷ số của chúng ta là 7 – 4 = 3. Điều này tốt hơn nhiều so với điểm -5 trước đây của chúng tôi. Lưu ý thứ tự của các nucleotide [i. e. các ký tự trong “chuỗi DNA”] không thay đổi theo cả hai trình tự; . Lý do làm điều này có ý nghĩa về mặt sinh học là việc thêm hoặc xóa nucleotide là các dạng đột biến phổ biến. Ví dụ, hai chuỗi DNA có thể liên quan đến chức năng, nhưng đột biến khiến các nucleotide bổ sung được thêm vào [hoặc bị xóa] khỏi một hoặc cả hai chuỗi

Bây giờ, để đưa ra một cách tổng quát hơn để tìm giải pháp đạt điểm tối đa, hãy xây dựng một bảng như bên dưới sao cho một trong các chuỗi chạy ngang qua tiêu đề, trong khi chuỗi kia chạy dọc theo các hàng

Giá trị của mỗi ô cụ thể sẽ biểu thị số điểm tối đa đạt được bằng cách ghép nối từng chuỗi DNA cho đến khi có nhiều hàng và cột đó. Ví dụ: -3 ở hàng trên cùng [không tính đến hàng tiêu đề] là giá trị -3 vì đó là điểm đạt được bằng cách ghép một khoảng trống [“-“] với ba nucelotide đầu tiên của chuỗi DNA tiêu đề, GTC

Tương tự như vậy, chúng tôi nhận được từng giá trị -1, …, -9 bằng cách ghép từng chuỗi con liên tiếp của cột tiêu đề với khoảng trống, “-”

-, G ==> -1
-, GT ==> -2
-, GTC ==> -3, -4
-, GTCCA ==> -5
-, GTCCAT ==> -6
-, GTCCATA ==> -7
-, GTCCATAC ==> -8
-, GTCCATACA ==> -9

Trong mỗi dãy con, mỗi nucleotide được so sánh với “-“. Vì không có nucleotide nào trong bất kỳ chuỗi con nào khớp với “-“, điểm số của mỗi nucleotide khớp với “-” chỉ bằng -1 lần độ dài của chuỗi con

Điền vào một hàng và cột khác sẽ cho chúng tôi điều này

Các tính toán cho các mục đầu tiên của hàng mới có thể được nhìn thấy dưới đây

G ====> G ====> -1 = -1
T       T

G T ====> G T ====> -1 + 1 = 0
T         - T


G T C ====> G T C ====> -1 + 1 - 1 = -1
T           - T -


G T C C ====> G T C C ====> -1 + 1 - 1 - 1 = -2
T             - T - -


G T C C A ====> G T C C A ====> -1 + 1 - 1 - 1 - 1 = -3
T               - T - - -


G T C C A T ====> G T C C A T ====> -1 + 1 - 1 - 1 - 1 - 1 = -4
T                 - T - - - -


G T C C A T A ====> G T C C A T A ====> -1 + 1 - 1 - 1 - 1 - 1 - 1 = -5
T                   - T - - - - - 




G T C C A T A C ====> G T C C A T A C ====> -1 + 1 - 1 - 1 - 1 - 1 - 1 - 1 = -6
T                     - T - - - - - -



G T C C A T A C A ====> G T C C A T A C A ====> -1 + 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 = -7
T                       - T - - - - - - -

Hoàn thành phần còn lại của bảng

Bây giờ chúng tôi hoàn thành phần còn lại của bảng

Các mục đường chéo được tính như thế này

G ====> G ====> -1 = -1
T       T


G T ====> G T - ====> -1 + 1 - 1 = -1
T C       - T -


G T C ====> G T C - ====> -1 + 1 + 1 - 1 = 0
T C A       - T C A


G T C C ====> G T C C - ====> -1 + 1 + 1 - 1 - 1 = -1
T C A T       - T C A T


G T C C A ====> G T C C - A ====> -1 + 1 + 1 - 1 - 1 + 1 = 0
T C A T A       - T C A T A


G T C C A T ====> G T C C - A T ====> -1 + 1 + 1 - 1 - 1 + 1 + 1 = 1
T C A T A T       - T C A T A T


G T C C A T A ====> G T C C A T A - - ====> -1 + 1 + 1 - 1 + 1 + 1 + 1 - 1 - 1 = 1
T C A T A T C       - T C - A T A T C


G T C C A T A C ====> G T C C A T A - C - ====> -1 + 1 + 1 - 1 + 1 + 1 + 1 - 1 + 1 - 1 = 2
T C A T A T C A       - T C - A T A T C A


G T C C A T A C A ====> G T C C A T A - C A - ====> -1 + 1 + 1 - 1 + 1 + 1 + 1 - 1 + 1 + 1 - 1 = 3
T C A T A T C A G       - T C - A T A T C A G

Triển khai Python

Bây giờ - làm thế nào để chúng tôi thực sự triển khai điều này trong Python? . Chúng tôi cũng sẽ xác định hàm GET_SCORE sẽ kiểm tra xem hai nucleotide đầu vào có bằng nhau hay không – trả về phần thưởng [mặc định là 1] nếu đúng và trả về hình phạt [mặc định -1] nếu không

import numpy as np

X = "GTCCATACA"
Y = "TCATATCAG"


def GET_SCORE[n1, n2, penalty = -1, reward = 1]:
    
    if n1 == n2:
        return reward
    else:
        return penalty


# initialize score matrix
score_matrix = np.ndarray[[len[X] + 1, len[Y] + 1]]
 
for i in range[len[X] + 1]:
	score_matrix[i, 0] = penalty * i

for j in range[len[Y] + 1]:
	score_matrix[0, j] = penalty * j

Tính giá trị của từng ô

Tiếp theo, chúng ta cần đặt giá trị của từng ô còn lại. Giá trị của mỗi ô được đệ quy dựa trên 1] các ô ở ngay bên trái, phía trên và theo đường chéo bên trái;

Cụ thể hơn, giá trị của ô ở vị trí thứ i,j của ma trận điểm [i. e. score_matrix[i, j]] được tính bằng cách lấy giá trị lớn nhất của các giá trị sau

score_matrix[i – 1, j – 1] + GET_SCORE[X[i – 1], Y[j – 1]]
score_matrix[i, j – 1] +
score_matrix[i – 1, j] + penalty

Các chỉ mục có thể hơi khó hiểu ở đây, vì vậy hãy tạm dừng một chút

G T C C
T C A T

G T C
T C A


G T C C
T C A


G T C
T C A T

def GET_SCORE[n1, n2, penalty = -1, reward = 1]:
    
    if n1 == n2:
        return reward
    else:
        return penalty


def global_alignment[X, Y, penalty = -1, reward = 1]:
    
    # initialize score matrix
    score_matrix = np.ndarray[[len[X] + 1, len[Y] + 1]]
     
    for i in range[len[X] + 1]:
        score_matrix[i, 0] = penalty * i
    
    for j in range[len[Y] + 1]:
        score_matrix[0, j] = penalty * j
        
    
    # define each cell in the matrix by as the max score possible in that stage
    for i in range[1, len[X] + 1]:
        for j in range[1, len[Y] + 1]:
            match = score_matrix[i - 1, j - 1] + GET_SCORE[X[i - 1], Y[j - 1], penalty, reward]
            delete = score_matrix[i -1, j] + penalty
            insert = score_matrix[i, j - 1] + penalty
            
            score_matrix[i, j] = max[[match, delete, insert]]
            
    
    i = len[X]
    j = len[Y]
    
    align_X = ""
    align_Y = ""
    
    while i > 0 or j > 0:
        
        current_score = score_matrix[i, j]
        left_score = score_matrix[i - 1, j]
        
        
        if i > 0 and j > 0 and X[i - 1] == Y[j - 1]:
            align_X = X[i - 1] + align_X
            align_Y = Y[j - 1] + align_Y
            i = i - 1
            j = j - 1
        
        elif i > 0 and current_score == left_score + penalty:
            align_X = X[i - 1] + align_X
            align_Y = "-" + align_Y
            i = i - 1
            
        else:
            align_X = "-" + align_X
            align_Y = Y[j - 1] + align_Y
            j = j - 1


    return align_X, align_Y

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G

Sequence 1 ==> G T C C A T A C A
Sequence 2 ==> T C A T A T C A G