局部比對算法 —— Smith-Waterman Algorithm
Swimt-Waterman算法本質(zhì)上是一種Dynamic Programming（動態(tài)規(guī)劃算法）蘑险，和Needleman算法有許多相同之處滴肿。其分為3個步驟：Initialization —— Matrix Filling —— Trace Back。
Swith-Waterman算法相較于Needleman-Wunsch算法最大的區(qū)別就在于佃迄，在“matrix preparation階段”泼差，其用0表示mismatch的數(shù)值和比對過程中產(chǎn)生的gap數(shù)值。
SW算法呵俏，先構(gòu)建一個（i+1）×（j+1）的矩陣堆缘。從左上角開始，以0作為其初始值普碎，計算第一行和第一列的數(shù)值（此步驟計算為gap penalty的累加）吼肥。
【標(biāo)注】此過程主要就是填充矩陣的第一行和第一列
由于SW算法將所有的負(fù)值用0進(jìn)行替換，在進(jìn)行上述計算步驟過后麻车，其結(jié)果如下：
下一步進(jìn)行的操作缀皱，稱為“matrix filling”。將原始矩陣劃分為4個一組2×2矩陣绪氛，每一個矩陣內(nèi)右下角可以有3種相對位置唆鸡，分別為left、diagonal枣察、upper争占。
每一種相對位置的移動，代表著不同的序列比對方式序目，如下：
1臂痕、垂直方向上的移動 & 水平方向上的移動，其含義為Indel & gap
2猿涨、對角線移動握童，其含義為“substitution”，其可能為match 或 mismatch
之前給出的位置為A-A叛赚，因此為“match”的情況澡绩，且又因垂直和水平方向上的移動，其對應(yīng)的數(shù)值為0俺附，因此在該2×2矩陣內(nèi)右下角位置對應(yīng)的數(shù)值為1肥卡。
最終就可以將矩陣填充為如下結(jié)果：
SW比對的最后一步為“Trace Back”。整體過程可以描述為事镣，從矩陣中的最大值開始步鉴，依次往回推，得到序列比對的結(jié)果。
全局比對算法 —— Needleman-Wunsch算法
作者提了一嘴氛琢，為什么需要做序列比對
- Functionnal Relationship
- Structural Relationship
- Evolutionary Relationship
2種算法（ND & SW）之間的區(qū)別：
ND算法喊递，并不會像SW算法一樣，將比對過程中得到的gap penalty和mismatch負(fù)值轉(zhuǎn)化為0阳似，而是保留了對應(yīng)負(fù)值骚勘，最終得到的矩陣結(jié)果為：
在Trace Back過程中，如果矩陣對應(yīng)位置的堿基是匹配關(guān)系（相同）障般，則下一步繼續(xù)走對角線调鲸，如果矩陣對應(yīng)位置的堿基不是相同的，則下一步需要走到的位置是2×2矩陣中最大值的方向挽荡。
可能出現(xiàn)的情況如下：
Needleman-Wunsch算法的Python實現(xiàn)藐石，可以通過這個視頻來學(xué)習(xí)學(xué)習(xí)：
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）于微，同時在最后得到的結(jié)果中，我們只想保留“statistically significant HSPs”青自，此步驟的篩選株依，參考序列比對的得分。
Based on the scores of aligning 2 random seqs
在2條序列比對過程中延窜，可能會出現(xiàn)多個HSP恋腕，如下：
參考文獻(xiàn)：
http://www.sciencedirect.com/science/article/pii/S0022283605803602
https://academic.oup.com/nar/article/25/17/3389/1061651/
（4）Suffix Array
【標(biāo)注】該方法用于構(gòu)建基因組的索引（e.g. STAR）
1、Suffix Tree
軟件：MUMmer
給定一條序列逆瑞，并基于該序列荠藤，生成其對應(yīng)的子序列。
e.g. BANANA
2获高、Suffix Array
軟件：STAR
（5）Borrows-Wheeler transformation & LF mapping
軟件：bwa哈肖，bowtie
1、Burrows-Wheeler Transform
BW算法概述：
- 給定一條序列T（acaacg$）念秧；
- 將序列尾部的字符轉(zhuǎn)移至序列頭部淤井，循環(huán)該過程直到得到原序列；
- 將序列進(jìn)行排序（$最小摊趾，其他按字母表排序）
- 記錄最后一列序列币狠，舍棄其余信息
為什么使用BW算法？
將上述給定序列轉(zhuǎn)化之后砾层，可以達(dá)到對字符進(jìn)行壓縮（text compression）的目的（acaacg3ac）
2总寻、LF mapping
LF mapping算法可以將BW算法記錄的算法還原。
3梢为、使用BW算法進(jìn)行序列比對
e.g. Query Q = aac, DB = T acaacg3ac
給出的序列第一位是a，第兩位是ac。需要注意的是铸董，這邊的查找是倒序的祟印，即caa。
首先粟害，基于最后一列數(shù)據(jù)找到第一列序列c的坐標(biāo)范圍（5 & 6）蕴忆，且由于c前2個字符為兩個a，以此作為條件進(jìn)行縮小搜索范圍的前提悲幅。
以字符c作為開頭的序列套鹅，對應(yīng)字符a在第一列的坐標(biāo)范圍為（3 & 4）。此步驟找到了query序列中第二個a的范圍汰具，下一步基于此過程的比對結(jié)果卓鹿，確定最后一個a的位置。
