序列匹配的算法

❶ 有关匹配和排序的算法，高手帮帮忙哈

一、插入排序(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》系列文章， 第三篇文章 《基因组序列比对算法介绍（一）》，争取每周更新一篇高质量生信干货帖子。

关注 "基因组分析" 微信公众号，了解最新最全生信分析知识。