【哈佛大学：计算生物学 & 生物信息学】学习记录（三）

作者: 陈有朴 | 来源:发表于2022-04-14 09:49 被阅读0次

【哈佛大学：计算生物学 & 生物信息学】学习记录（三）
【哈佛大学：计算生物学 & 生物信息学】学习记录（二）
【哈佛大学：计算生物学 & 生物信息学】学习记录（一）
【哈佛大学：计算生物学 & 生物信息学】学习记录（五）
生物信息入门基础
老司机辣评 | 研究中国人自己的基因数据库项目近日开启
领域内文章新消息（生信、表观，不定期更新）
生物信息学
Bioinformatics学习实战——Nature文章重复
比对算法总结(二)——基于BWT索引结构的比对算法-Bowite

局部比对算法 —— Smith-Waterman Algorithm

Swimt-Waterman算法本质上是一种Dynamic Programming（动态规划算法），和Needleman算法有许多相同之处。其分为3个步骤：Initialization —— Matrix Filling —— Trace Back。

Swith-Waterman算法相较于Needleman-Wunsch算法最大的区别就在于，在“matrix preparation阶段”，其用0表示mismatch的数值和比对过程中产生的gap数值。

SW算法，先构建一个（i+1）×（j+1）的矩阵。从左上角开始，以0作为其初始值，计算第一行和第一列的数值（此步骤计算为gap penalty的累加）。
【标注】此过程主要就是填充矩阵的第一行和第一列

由于SW算法将所有的负值用0进行替换，在进行上述计算步骤过后，其结果如下：

下一步进行的操作，称为“matrix filling”。将原始矩阵划分为4个一组2×2矩阵，每一个矩阵内右下角可以有3种相对位置，分别为left、diagonal、upper。

每一种相对位置的移动，代表着不同的序列比对方式，如下：

1、垂直方向上的移动 & 水平方向上的移动，其含义为Indel & gap

2、对角线移动，其含义为“substitution”，其可能为match 或 mismatch

之前给出的位置为A-A，因此为“match”的情况，且又因垂直和水平方向上的移动，其对应的数值为0，因此在该2×2矩阵内右下角位置对应的数值为1。

最终就可以将矩阵填充为如下结果：

SW比对的最后一步为“Trace Back”。整体过程可以描述为，从矩阵中的最大值开始，依次往回推，得到序列比对的结果。

全局比对算法 —— Needleman-Wunsch算法

作者提了一嘴，为什么需要做序列比对

Functionnal Relationship
Structural Relationship
Evolutionary Relationship

2种算法（ND & SW）之间的区别：

ND算法，并不会像SW算法一样，将比对过程中得到的gap penalty和mismatch负值转化为0，而是保留了对应负值，最终得到的矩阵结果为：

在Trace Back过程中，如果矩阵对应位置的碱基是匹配关系（相同），则下一步继续走对角线，如果矩阵对应位置的碱基不是相同的，则下一步需要走到的位置是2×2矩阵中最大值的方向。

可能出现的情况如下：

Needleman-Wunsch算法的Python实现，可以通过这个视频来学习学习：
https://www.youtube.com/watch?v=um8h3P216Fk

代码我也附带敲好了，如下：

#----------------Needleman-Wunsch-------------------#
# Importiing Modules
from threading import main_thread
import numpy as np
from pip import main

# Ask for sequences from the user
# sequence_1 = input("Enter or paste sequence 1:")
# sequence_2 = input("Enter or paste sequence 2:")


sequence_1 = "ATCGT"
sequence_2 = "ACGT"

# Create Matrices
main_matrix = np.zeros((len(sequence_1)+1, len(sequence_2)+1))
match_checker_matrix = np.zeros((len(sequence_1), len(sequence_2)))

# Providing the scores for match, mismatch and gap
match_reward = 1
mismatch_penalty = -1
gap_penalty = -2

# Fill the match checker matrix according to match or mismatch
for i in range(len(sequence_1)):
    for j in range(len(sequence_2)):
        if sequence_1[i] == sequence_2[j]:
            match_checker_matrix[i][j] = match_reward
        else:
            match_checker_matrix[i][j] = mismatch_penalty

# print(match_checker_matrix)

# Filling up the matrix using Needleman_Wunsch algorithm
# STEP 1: Initialization
for i in range(len(sequence_1)+1):
    main_matrix[i][0] = i*gap_penalty
for j in range(len(sequence_2)+1):
    main_matrix[0][j] = j*gap_penalty

# STEP 2: Matrix Filling
for i in range(1, len(sequence_1)+1):
    for j in range(1, len(sequence_2)+1):
        main_matrix[i][j] = max(main_matrix[i-1][j-1] + match_checker_matrix[i-1][j-1],
        main_matrix[i-1][j] + gap_penalty,
        main_matrix[i][j-1] + gap_penalty)

# print(main_matrix)

# STEP 3: TraceBack
aligned_1 = ""
aligned_2 = ""

ti = len(sequence_1)
tj = len(sequence_2)

while(ti > 0 and tj > 0):

    if (ti > 0 and tj > 0 and main_matrix[ti][tj] == main_matrix[ti-1][tj-1] + match_checker_matrix[ti-1][tj-1]):

        aligned_1 = sequence_1[ti - 1] + aligned_1
        aligned_2 = sequence_2[tj - 1] + aligned_2

        ti = ti - 1
        tj = tj - 1

    elif (ti > 0 and main_matrix[ti][tj] == main_matrix[ti-1][tj] + gap_penalty):
        aligned_1 = sequence_1[ti-1] + aligned_1
        aligned_2 = "-" + aligned_2
        
        ti = ti - 1
    else:
        aligned_1 = "-" + aligned_1
        aligned_2 = sequence_2[tj-1] + aligned_2

# test
print(aligned_1)
print(aligned_2)