局部比对算法 —— 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)
序列比对算法
（1）Read Mapping
（2）Seed —— Kmer
（3）BLAST
在BLAST序列比对过程中，可以得到HSP（high-scoring segment pairs），同时在最后得到的结果中，我们只想保留“statistically significant HSPs”，此步骤的筛选，参考序列比对的得分。
Based on the scores of aligning 2 random seqs
在2条序列比对过程中，可能会出现多个HSP，如下：
参考文献：
http://www.sciencedirect.com/science/article/pii/S0022283605803602
https://academic.oup.com/nar/article/25/17/3389/1061651/
（4）Suffix Array
【标注】该方法用于构建基因组的索引（e.g. STAR）
1、Suffix Tree
软件：MUMmer
给定一条序列，并基于该序列，生成其对应的子序列。
e.g. BANANA
2、Suffix Array
软件：STAR
（5）Borrows-Wheeler transformation & LF mapping
软件：bwa，bowtie
1、Burrows-Wheeler Transform
BW算法概述：
- 给定一条序列T（acaacg$）；
- 将序列尾部的字符转移至序列头部，循环该过程直到得到原序列；
- 将序列进行排序（$最小，其他按字母表排序）
- 记录最后一列序列，舍弃其余信息
为什么使用BW算法？
将上述给定序列转化之后，可以达到对字符进行压缩（text compression）的目的（acaacg3ac）
2、LF mapping
LF mapping算法可以将BW算法记录的算法还原。
3、使用BW算法进行序列比对
e.g. Query Q = aac, DB = T acaacg3ac
给出的序列第一位是a，第两位是ac。需要注意的是，这边的查找是倒序的，即caa。
首先，基于最后一列数据找到第一列序列c的坐标范围（5 & 6），且由于c前2个字符为两个a，以此作为条件进行缩小搜索范围的前提。
以字符c作为开头的序列，对应字符a在第一列的坐标范围为（3 & 4）。此步骤找到了query序列中第二个a的范围，下一步基于此过程的比对结果，确定最后一个a的位置。
网友评论