序列匹配的演算法

❶ 有關匹配和排序的演算法，高手幫幫忙哈

一、插入排序(Insertion Sort)
1. 基本思想：
每次將一個待排序的數據元素，插入到前面已經排好序的數列中的適當位置，使數列依然有序；直到待排序數據元素全部插入完為止。
2. 排序過程：
【示例】：
[初始關鍵字] [49] 38 65 97 76 13 27 49
J=2(38) [38 49] 65 97 76 13 27 49
J=3(65) [38 49 65] 97 76 13 27 49
J=4(97) [38 49 65 97] 76 13 27 49
J=5(76) [38 49 65 76 97] 13 27 49
J=6(13) [13 38 49 65 76 97] 27 49
J=7(27) [13 27 38 49 65 76 97] 49
J=8(49) [13 27 38 49 49 65 76 97]

Procere InsertSort(Var R : FileType);
//對R[1..N]按遞增序進行插入排序, R[0]是監視哨//
Begin
for I := 2 To N Do //依次插入R[2],...,R[n]//
begin
R[0] := R[I]; J := I - 1;
While R[0] < R[J] Do //查找R[I]的插入位置//
begin
R[J+1] := R[J]; //將大於R[I]的元素後移//
J := J - 1
end
R[J + 1] := R[0] ; //插入R[I] //
end
End; //InsertSort //

二、選擇排序
1. 基本思想：
每一趟從待排序的數據元素中選出最小（或最大）的一個元素，順序放在已排好序的數列的最後，直到全部待排序的數據元素排完。
2. 排序過程：
【示例】：
初始關鍵字 [49 38 65 97 76 13 27 49]
第一趟排序後 13 〔38 65 97 76 49 27 49]
第二趟排序後 13 27 〔65 97 76 49 38 49]
第三趟排序後 13 27 38 [97 76 49 65 49]
第四趟排序後 13 27 38 49 [49 97 65 76]
第五趟排序後 13 27 38 49 49 [97 97 76]
第六趟排序後 13 27 38 49 49 76 [76 97]
第七趟排序後 13 27 38 49 49 76 76 [ 97]
最後排序結果 13 27 38 49 49 76 76 97

Procere SelectSort(Var R : FileType); //對R[1..N]進行直接選擇排序 //
Begin
for I := 1 To N - 1 Do //做N - 1趟選擇排序//
begin
K := I;
For J := I + 1 To N Do //在當前無序區R[I..N]中選最小的元素R[K]//
begin
If R[J] < R[K] Then K := J
end;
If K <>; I Then //交換R[I]和R[K] //
begin Temp := R[I]; R[I] := R[K]; R[K] := Temp; end;
end
End; //SelectSort //

三、冒泡排序(BubbleSort)
1. 基本思想：
兩兩比較待排序數據元素的大小，發現兩個數據元素的次序相反時即進行交換，直到沒有反序的數據元素為止。
2. 排序過程：
設想被排序的數組R〔1..N〕垂直豎立，將每個數據元素看作有重量的氣泡，根據輕氣泡不能在重氣泡之下的原則，從下往上掃描數組R，凡掃描到違反本原則的輕氣泡，就使其向上"漂浮"，如此反復進行，直至最後任何兩個氣泡都是輕者在上，重者在下為止。
【示例】：
49 13 13 13 13 13 13 13
38 49 27 27 27 27 27 27
65 38 49 38 38 38 38 38
97 65 38 49 49 49 49 49
76 97 65 49 49 49 49 49
13 76 97 65 65 65 65 65
27 27 76 97 76 76 76 76
49 49 49 76 97 97 97 97

Procere BubbleSort(Var R : FileType) //從下往上掃描的起泡排序//
Begin
For I := 1 To N-1 Do //做N-1趟排序//
begin
NoSwap := True; //置未排序的標志//
For J := N - 1 DownTo 1 Do //從底部往上掃描//
begin
If R[J+1]< R[J] Then //交換元素//
begin
Temp := R[J+1]; R[J+1 := R[J]; R[J] := Temp;
NoSwap := False
end;
end;
If NoSwap Then Return//本趟排序中未發生交換，則終止演算法//
end
End; //BubbleSort//

四、快速排序（Quick Sort）
1. 基本思想：
在當前無序區R[1..H]中任取一個數據元素作為比較的"基準"(不妨記為X)，用此基準將當前無序區劃分為左右兩個較小的無序區：R[1..I-1]和R[I+1..H]，且左邊的無序子區中數據元素均小於等於基準元素，右邊的無序子區中數據元素均大於等於基準元素，而基準X則位於最終排序的位置上，即R[1..I-1]≤X.Key≤R[I+1..H](1≤I≤H)，當R[1..I-1]和R[I+1..H]均非空時，分別對它們進行上述的劃分過程，直至所有無序子區中的數據元素均已排序為止。
2. 排序過程：
【示例】：
初始關鍵字 [49 38 65 97 76 13 27 49〕
第一次交換後 〔27 38 65 97 76 13 49 49〕
第二次交換後 〔27 38 49 97 76 13 65 49〕
J向左掃描，位置不變，第三次交換後 〔27 38 13 97 76 49 65 49〕
I向右掃描，位置不變，第四次交換後 〔27 38 13 49 76 97 65 49〕
J向左掃描 〔27 38 13 49 76 97 65 49〕
（一次劃分過程）

初始關鍵字 〔49 38 65 97 76 13 27 49〕
一趟排序之後 〔27 38 13〕 49 〔76 97 65 49〕
二趟排序之後 〔13〕 27 〔38〕 49 〔49 65〕76 〔97〕
三趟排序之後 13 27 38 49 49 〔65〕76 97
最後的排序結果 13 27 38 49 49 65 76 97
各趟排序之後的狀態

Procere Parttion(Var R : FileType; L, H : Integer; Var I : Integer);
//對無序區R[1,H]做劃分，I給以出本次劃分後已被定位的基準元素的位置 //
Begin
I := 1; J := H; X := R[I] ;//初始化，X為基準//
Repeat
While (R[J] >;= X) And (I < J) Do
begin
J := J - 1 //從右向左掃描，查找第1個小於 X的元素//
If I < J Then //已找到R[J] 〈X//
begin
R[I] := R[J]; //相當於交換R[I]和R[J]//
I := I + 1
end;
While (R[I] <= X) And (I < J) Do
I := I + 1 //從左向右掃描，查找第1個大於 X的元素///
end;
If I < J Then //已找到R[I] >; X //
begin R[J] := R[I]; //相當於交換R[I]和R[J]//
J := J - 1
end
Until I = J;
R[I] := X //基準X已被最終定位//
End; //Parttion //

Procere QuickSort(Var R :FileType; S,T: Integer); //對R[S..T]快速排序//
Begin
If S < T Then //當R[S..T]為空或只有一個元素是無需排序//
begin
Partion(R, S, T, I); //對R[S..T]做劃分//
QuickSort(R, S, I-1);//遞歸處理左區間R[S,I-1]//
QuickSort(R, I+1,T);//遞歸處理右區間R[I+1..T] //
end;
End; //QuickSort//

五、堆排序(Heap Sort)
1. 基本思想：
堆排序是一樹形選擇排序，在排序過程中，將R[1..N]看成是一顆完全二叉樹的順序存儲結構，利用完全二叉樹中雙親結點和孩子結點之間的內在關系來選擇最小的元素。
2. 堆的定義: N個元素的序列K1,K2,K3,...,Kn.稱為堆，當且僅當該序列滿足特性：
Ki≤K2i Ki ≤K2i+1(1≤ I≤ [N/2])

堆實質上是滿足如下性質的完全二叉樹：樹中任一非葉子結點的關鍵字均大於等於其孩子結點的關鍵字。例如序列10,15,56,25,30,70就是一個堆，它對應的完全二叉樹如上圖所示。這種堆中根結點（稱為堆頂）的關鍵字最小，我們把它稱為小根堆。反之，若完全二叉樹中任一非葉子結點的關鍵字均大於等於其孩子的關鍵字，則稱之為大根堆。
3. 排序過程：
堆排序正是利用小根堆（或大根堆）來選取當前無序區中關鍵字小（或最大）的記錄實現排序的。我們不妨利用大根堆來排序。每一趟排序的基本操作是：將當前無序區調整為一個大根堆，選取關鍵字最大的堆頂記錄，將它和無序區中的最後一個記錄交換。這樣，正好和直接選擇排序相反，有序區是在原記錄區的尾部形成並逐步向前擴大到整個記錄區。
【示例】：對關鍵字序列42，13，91，23，24，16，05，88建堆

Procere Sift(Var R :FileType; I, M : Integer);
//在數組R[I..M]中調用R[I]，使得以它為完全二叉樹構成堆。事先已知其左、右子樹(2I+1 <=M時)均是堆//
Begin
X := R[I]; J := 2*I; //若J <=M, R[J]是R[I]的左孩子//
While J <= M Do //若當前被調整結點R[I]有左孩子R[J]//
begin
If (J < M) And R[J].Key < R[J+1].Key Then
J := J + 1 //令J指向關鍵字較大的右孩子//
//J指向R[I]的左、右孩子中關鍵字較大者//
If X.Key < R[J].Key Then //孩子結點關鍵字較大//
begin
R[I] := R[J]; //將R[J]換到雙親位置上//
I := J ; J := 2*I //繼續以R[J]為當前被調整結點往下層調整//
end;
Else
Exit//調整完畢，退出循環//
end
R[I] := X;//將最初被調整的結點放入正確位置//
End;//Sift//

Procere HeapSort(Var R : FileType); //對R[1..N]進行堆排序//
Begin
For I := N Div Downto 1 Do //建立初始堆//
Sift(R, I , N)
For I := N Downto 2 do //進行N-1趟排序//
begin
T := R[1]; R[1] := R[I]; R[I] := T;//將當前堆頂記錄和堆中最後一個記錄交換//
Sift(R, 1, I-1) //將R[1..I-1]重成堆//
end
End; //HeapSort//

六、幾種排序演算法的比較和選擇
1. 選取排序方法需要考慮的因素：
(1) 待排序的元素數目n；
(2) 元素本身信息量的大小；
(3) 關鍵字的結構及其分布情況；
(4) 語言工具的條件，輔助空間的大小等。
2. 小結：
(1) 若n較小(n <= 50)，則可以採用直接插入排序或直接選擇排序。由於直接插入排序所需的記錄移動操作較直接選擇排序多，因而當記錄本身信息量較大時，用直接選擇排序較好。
(2) 若文件的初始狀態已按關鍵字基本有序，則選用直接插入或冒泡排序為宜。
(3) 若n較大，則應採用時間復雜度為O(nlog2n)的排序方法：快速排序、堆排序或歸並排序。 快速排序是目前基於比較的內部排序法中被認為是最好的方法。
(4) 在基於比較排序方法中，每次比較兩個關鍵字的大小之後，僅僅出現兩種可能的轉移，因此可以用一棵二叉樹來描述比較判定過程，由此可以證明：當文件的n個關鍵字隨機分布時，任何藉助於"比較"的排序演算法，至少需要O(nlog2n)的時間。
(5) 當記錄本身信息量較大時，為避免耗費大量時間移動記錄，可以用鏈表作為存儲結構。

❷ Needleman-Wunsch 演算法及實現

Needleman-Wunsch (NW) 演算法常用於基因序列匹配以及單詞匹配，是一種實現全局最優匹配的動態規劃演算法。其最差情況下的時間復雜度是 , 空間復雜度是 ，其中 和 是兩個序列的長度。

不妨設有兩個序列：
GCATGCU
GATTACA

NW演算法的具體步驟如下：

運行結果：

❸ 基因組序列比對演算法介紹（一）

基因組重測序中序列比對介紹

重測序基因組數據比對，是指將測序儀下機fastq數據（NGS read序列，通常100-150bp），與人類參考基因組（reference）進行匹配，允許錯配（mismatch），插入缺失（indel），目的是在參考基因組找到序列最相似的位置，通常是基因組分析(包括 variation calling，ChIP-seq，RNA-seq，BS-seq)流程的第一步。

常用演算法

圖一

漢明距離（Hamming distance）表示兩個（相同長度）字對應位置不同的數量，我們以d（x,y）表示兩個字x,y之間的漢明距離。對兩個字元串進行異或運算，並統計結果為1的個數，那麼這個數就是漢明距離。圖中read1最佳位置的方法，就是通過查找最小漢明距離的實現的。

編輯距離（Edit distance）是針對二個字元串（例如英文字）的差異程度的量化量測，量測方式是看至少需要多少次的處理才能將一個字元串變成另一個字元串。圖中read3最佳位置，通過查找最我輯距離的方法實現。

圖二

全局比對（Global alignment）：全局比對是指將參與比對的兩條序列裡面的所有字元進行比對。全局比對在全局范圍內對兩條序列進行比對打分，找出最佳比對，主要被用來尋找關系密切的序列。其可以用來鑒別或證明新序列與已知序列家族的同源性，是進行分子進化分析的重要前提。其代表是Needleman-Wunsch演算法。圖一中，read3使用全部比對。

局部比對（Local alignment）：與全局比對不同，局部比對不必對兩個完整的序列進行比對，而是在每個序列中使用某些局部區域片段進行比對。其產生的需求在於、人們發現有的蛋白序列雖然在序列整體上表現出較大的差異性，但是在某些局部區域能獨立的發揮相同的功能，序列相當保守。這時候依靠全局比對明顯不能得到這些局部相似序列的。其次，在真核生物的基因中，內含子片段表現出了極大變異性，外顯子區域卻較為保守，這時候全局比對表現出了其局限性，無法找出這些局部相似性序列。其代表是Smith-Waterman局部比對演算法。圖一中，read2使用局部比對。

圖三

Smith-Waterman演算法介紹

Smith-Waterman是由Temple F. Smith和Michael S. Waterman於1981年提出的一種進行局部序列比對（相對於全局比對）的演算法，用於找出兩個核苷酸序列或蛋白質序列之間的相似區域。該演算法的目的不是進行全序列的比對，而是找出兩個序列中具有高相似度的片段。S-W演算法基於動態規劃，它接受任意長度、任意位置、任意序列的對齊，並確定是否能找到最優的比對。

簡單地說就是，動態規劃找到問題中較小部分的解，然後把它們放在一起，形成整個問題的一個完整的最優最終解。

它優於BLAST和FASTA演算法，因為它搜索了更大的可能性，具有更高的敏感性。

S-W演算法不是一次查看整個序列，而是對多個長度的片段進行比較，尋找能夠最大化得分的片段。演算法本身本質上是遞歸的：

圖四

演算法步驟如下：

基因組分析***** 微信 公眾號推出 《50篇文章深入理解NGS》系列文章， 第三篇文章 《基因組序列比對演算法介紹（一）》，爭取每周更新一篇高質量生信干貨帖子。

關注 "基因組分析" 微信公眾號，了解最新最全生信分析知識。