中 華 大 學 碩 士 論 文
考慮凹形間隔處罰函數之有效率生物序列 比對演算法
Efficient Algorithms of Sequence Alignments with Concave Gap Penalties
學 系 別：生物資訊學系碩士班 學號姓名：M10120001 蔡 維 謙 指導教授：吳 哲 賢 博 士
中 華 民 國 104 年 8 月
中文摘要
生物序列比對問題是生物資訊學中，非常基礎及重要的研究課題，如何讓序 列比對的時間減少，而達到有效率的生物序列比對，一直是一個被廣泛討論的問 題。
針對序列比對問題，假設 n 及 m 為比對兩序列長度，n ≧ m。首先是由 Needleman 及 Wunsch 於 1970 年，針對無間隔處罰函數，提出 O(nm)演算法。接 著 Waterman、Smith 及 Beyer 於 1976 年，提出考慮任意間隔處罰函數，時間複 雜度 O(nm2) 的演算法。Gotoh 於 1982 年針對線性間隔處罰，提出時間複雜度為 O(nm)的演算法。針對击形間隔處罰函數，Myers 及 Miller 於 1988 年、Galil 及 Giancarlo 於 1989 年，以及 Gusfield 於 1997 年，也分別提出了時間複雜度 O(nmlogm)的演算法。
然而到目前為止，尚未有學者針對凹形間隔處罰函數的序列比對問題，提出 相關研究。本篇論文中我們是針對击形間隔處罰函數序列比對問題，深入討論並 研究分析其特性後，進而推導出時間複雜度 O(nmlogm)有效率的演算法。
關鍵字：演算法、序列比對、凹形間隔處罰函數
ii
ABSTRACT
Biological sequence alignments problem is a very basic and important researches on Bioinformatics, how to reduce the excution time and obtain efficient sequence alignments has been widely discussed.
For sequence alignments problem, where n, m are the sequence lengths, n≧
m,Needleman and Wunsch first proposed the O(nm) algorithms in 1970. And then, Waterman, Smith and Beyer in 1976, proposed O(nm2) algorithms for affine gappenalties function. In 1982, Gotoh proposedO(nm) algorithmsfor linear gap penalties. For convex gappenalties,Myers and Miller in 1988, Galil and Giancarlo in 1989, and Gusfield in 1997, they also proposed O(nmlogm) algorithms.
Until now, there is no scholarsto propse researches about sequence alignments problem with concave gappenalties. In this paper, we focus our researches on sequence alignments problem with convcave gappenalties, and finally proposed O(nmlogm) efficient algorithms.
Keywords: Algortihms, Sequence alignments,Concavegap penalties
誌謝
我能夠順利的完成研究所的學業，發表論文，首先必頇感謝我的指導教授吳 哲賢教授，在我的研究所課業上的督促及提攜，碰到難題時也不厭其煩的為我講 解，甚至也會關心我在生活上的一些事情，讓我感覺教授不只是教授，像是我的 家人一樣，給了我心中滿滿的溫暖，最後完成了研究並順利將論文發表，真的非 常感謝教授。
再來我要感謝我的家人，家人對我在生活上的幫忙，平時間的噓寒問暖，雪 中送炭，都讓我感動不已，對我服完役還要回來念書這件事大力的支持與贊成，
讓我能夠全心全意把重心放在學業上，真的很感謝我的父母。
然後我還要謝謝我同期的研究生同學們，有的時候我沒有在學校的時候，如 果臨時有什麼事情，彼此間都會互相通知，讓我不會錯漏任何一絲消息或資訊，
平時相處起來也很融洽，還會相約打球放鬆身心，真的很謝謝他們。
iv
目錄
中文摘要 ... I ABSTRACT ... II 誌謝 ... III 目錄 ... IV 圖目錄 ... V
第一章導論... 1
1.1 序列比對 ... 2
1.2 間隔處罰函數 ... 10
1.3 任意間隔處罰函數 ... 12
1.4 研究目的與動機 ... 15
第二章相關研究... 17
2.1 击形間隔處罰函數研究 ... 18
2.2GUSFIELD的击形間隔處罰函數 ... 21
2.3 击函數演算法兩大步驟說明 ... 27
第三章我們的演算法... 28
3.1 凹形間隔處罰函數研究 ... 29
3.2 我們的演算法 ... 32
3.3 凹函數演算法兩大步驟說明 ... 38
第四章結論... 40
4.1 研究成果 ... 40
4.2 未來研究目標 ... 41
參考文獻 ... 43
圖目錄
圖 1.1：兩條序列比對的其中兩種得分可能。 ... 4
圖 1.2：V(I , J)的分數。 ... 6
圖 1.3：全域序列比對示意圖。 ... 8
圖 1.4：序列比對最後結果範例。 ... 9
圖 1.5：兩條序列插入間隔後下去比對的情況。 ... 10
圖 1.6：間隔處罰的三種函數。 ... 11
圖 1.7：目前為止的研究發展。 ... 16
圖 2.1：击形函數特性。 ... 18
圖 2.2：利用動態規劃處理击形間隔處罰函數的符號設定。 ... 21
圖 2.3：K、J、J'、J"在候選者表單裡的取捨關係。 ... 22
圖 2.4：B矩陣示意圖。 ... 23
圖 25：根據定理 2.1 可能遇到的三種情況。 ... 24
圖 2.6：B矩陣更新時的變化。 ... 25
圖 2.7：击函數的 E 表格。 ... 27
圖 2.8：击函數的Ē 表格。 ... 27
圖 2.9：箭頭指示的地區為使用BINARY SEARACH搜尋到B值的地方。 ... 28
圖 3.1：凹型函數特性。 ... 29
圖 3.2：在凹型間隔處罰函數裡K < J < J’< J”。 ... 32
圖 3.3：凹函數的B矩陣示意圖。 ... 34
圖 3.4：分區的時候會碰到的三種情形。 ... 35
圖 3.5：矩陣內內容變化的範例。 ... 36
圖 3.6：凹函數的 E 表格。 ... 36
圖 3.7：凹函數的Ē 表格。 ... 36
圖 3.8：箭頭指示的地區為使用BINARY SEARACH搜尋到B值的地方。 ... 36
圖 4.1：增函數示意圖。 ... 41
圖 4.2：未來的研究目標。 ... 42
1
第一章導論
生物序列比對（sequence alignments）問題，在生物資訊學中是一個很基礎 且很重要的研究課題，兩條長度有限但長度不同的生物序列，藉由插入間隔形成 兩條長度相同的序列，再根據評分矩陣以及間隔處罰（gap penalties）函數計算 分數，並找出分數最高的序列比對結果。
也因為插入間隔的關係，所以就有了間隔處罰函數的產生，間隔處罰函數主 要有分為任意（General）、線性（Affine）、以及击形（Convex）等方式。首先是 由 Needleman 及 Wunsch 在 1970 年的時候針對無間隔處罰函數之序列比對，利 用動態規劃(dynamic programming)的方法，先提出了 O(nm)的演算法，其中 n、
m 為序列長度，n≧m。接著 Waterman、Smith 及 Beyer 於 1976 年，則是提出了 考慮任意間隔處罰函數，時間複雜度為 O(nm2) 的演算法。Gotoh 於 1982 年針對 線性間隔處罰(affine gap penaltines) 函數，提出時間複雜度是 O(nm)的演算法。
Waterman 於 1984 年，提出击形間隔處罰(convex gap penaltines) 函數，並 預測存在時間複雜度為 O(nmlogm)。Myers 及 Miller 於 1988 年，針對击形間隔 處罰函數，首先提出時間複雜度為 O(nmlogm)的方法。爾後 Galil 及 Giancarlo 於 1989 年，以及 Gusfield 於 1997 年，也分別提出了時間複雜度是 O(nmlogm) 的改良演算法。
到目前為止，尚未有學者針對凹形間隔處罰(concave gap penaltines) 函數，
提出相關研究。本篇論文是針對討論考慮凹形間隔處罰函數之序列比對問題，利 用 Gusfield 針對击形間隔處罰設計的演算法，分析其特性並改良，進一步推導出 有效率的時間複雜度 O(nmlogm)的凹形演算法。
1.1 序列比對
基因序列，也就是去氧核醣核酸 DNA（Deoxyribonucleic acid），主要是存在 於染色體裡面，而在染色體中則包含有大量的基因序列，所謂的基因序列也就是 核酸序列，主要是由四種鹼基所組成，分別為，腺嘌呤（Adenine）、胸腺嘧啶
（Thymine）、鳥糞嘌呤（Guanine）、胞嘧啶（Cytosine），即是 A、T、C、G，由 四個以上的鹼基對所組合成的序列即為生物序列[1-4]。
序列比對是指將兩條或是多個序列排列在一起，兩條序列的長度不同的時候 藉由插入間隔（間隔通常用短橫線「-」表示）來讓兩條序列長度變成一致，來 對應蛋白質或胺基酸中的 A、T（U)、C、G，讓兩條甚至是多條序列的長短變的 一樣，再來觀察它們之間的相似程度，兩條序列在對齊的時候只會有三種情況，
字母碰字母、字母碰間隔、間隔碰字母，不能有間隔碰間隔的情況，因為這在比 對上並沒有任何的意義[5]。
所以，序列比對在我們現在的生物學裡面扮演了很重要的角色，因為藉由比 對以及分析，我們可以找出兩條序列之間的最大相似度，而相似度越高，則表示 這兩條序列的功能或是它們的表徵可能就越接近，甚至有可能代表這兩個生物的 祖先在久遠的遠古時期的時候是來自同一種生物所演化而成的現今的兩個不同 的物種。
現在存活在地球上的許多生物都是從以前開始經由無視次的演化而來的，而 演化裡最大的關鍵就在於基因序列裡產生了突變(mutation)，突變主要是生物的 序列在進行複製(duplication)以及轉譯(transcription)的時候發生的。會發生突變的 原因有很多種可能，例如細胞在進行細胞分裂的時候複製發生錯誤，或是受到化 學物質、輻射、病毒的影響而發生突變，突變主要的形式有三種：替代(substitution)、
插入(insert)、刪除(delete)；那也因為突變的關係所以才造就了現在地球上的生物 多樣性，而越相似的序列他們彼此之間也有可能具有相似的特徵或是功能。
3
比對分數
在這裡我們假設生物會以最容易發生突變的方式來進行突變，所以在序列比 對上我們將給予計分，計分方式是由我們研究人員設定的，一般來說，比對相似 的時候會給予正分，比對不相似的時候會給予負分，最後加總起來的就是總分，
總分越高則代表兩條序列在比對的結果上來說，相似的地方較多，間隔插入越多 以及比對不相同的地方，則扣的分數越多，相對的得到的總分會越低，基於符合 生物上的意義，在比對的時候，我們則是要找出扣的分數最少，也就是得分最高 的比對結果，即為此兩序列間最佳比對方式，藉此我們可以輕易的找到兩種生物 之間是否有同源關係，或是相似的構造或功能。在序列比對上會發生三種情況：
一、配對相同 二、配對不相同 三、插入間隔
在這裡我們將第三項考慮為跟第二項一樣是配對不相同的情況，配對不相同者我 們給予扣分(-1)，配對相同者給予正分(+5)，以先不考慮間隔處罰問題來進行比 對分數計算。
舉例說明：
圖 1.1：兩條序列比對的其中兩種得分可能。
從圖 1.1 裡我們可以看出，一樣的兩個序列進行比對，比對出來的結果有兩 種可能性，甚至還會有更多種的比對可能性存在，但是我們主要是找到比對得分 最高的比對結果，第一個 CASE 比對結果的得分較高，而第二個 CASE 的比對得 分則較低，很明顯的，第一個 CASE 的比對結果優於第二個比對結果，當然比對 的結果可能會有很多種，分數也會有很多種的可能性，所以我們要找出兩條序列 中最高的比對總分，因為那代表著兩條序列比對出來的最相似的結果[6-8]。
但是基於符合生物意義的形式來說，在生物演化上我們不僅要考慮單一一個 間隔，還要去考慮到間隔的長度問題，而且不同長度的間隔我們則給予不同的扣 分弧度，這樣才符合生物意義，所以間隔處罰函數也因此而產生，關於間隔處罰 函數，在稍後的 1.2 章節我們會再做討論。
5
動態規劃
若是我們以兩條序列長度同樣為 n 的情況下，來進行序列比對，比對結果的 時間複雜度為 n 的指數次方，任何一條胺基酸序列的鹼基排序皆由數萬甚至數十 萬以上個鹼基所組成，全部處理起來是個很龐大費時的工作，在以前沒有電腦的 年代，生物學家們若是要完成序列的比對，必頇將要處理的序列以人工的方式下 去對齊兩個甚至是多個序列，因為比對出來的結果會有很多種，所以這些作業處 理起來往往就是一年半載的，相當費時與耗功夫；雖然在現在科技進步的時代，
有電腦可以當作輔助，但是得到的效果有限，因為光是以人類的整個基因體來看，
光是總長度就將近有 30 億左右甚至以上，就算用當今最強的超級電腦可能都還 不一定跑的完。
為了縮短應付這龐大的數據處理的時間，所以在 1970 年的時候，Needleman 及 Wunsch 將動態規劃(Dynamic Programming)的概念應用到序列比對上，將處理 的時間縮短，讓時間複雜度降至只要 n2的時間即可處理完成，這個重大的突破，
讓後續有許多學者根據動態規劃的理念繼續持續著改良以及加速演算法的運作 [9-10]。
全局比對（Global alignments）
全局比對就是兩段長度分別為ｍ及ｎ的序列下去進行比對，我們假設m ≥ n，
找出兩條序列的最佳比對結果，1970 年時，Needleman 及 Wunsch 提出了利用動 態規劃(dynamic programming)及遞迴參數(recursive argument)來進行序列比對的 計分[11]。
我們假設有兩條序列 S 和 T，長度分別為 n 及 m ，給一個 V(i , j)的二維 矩陣，分別為 S 序列的第 i 個元素與 T 序列的第 j 個元素的比對結果分數，示意 圖如下：
圖 1.2：V(i , j)的分數是考慮來自 V(i-1,j-1)、V(i-1,j)、V(i,j-1)的分數並從裡面挑 出最大值來當 V(i , j)的值。
要計算出 V(i,j)的值首先必頇先考慮到來自三個方向的數值然後選擇裡面 的最大分數來當 V(i,j)的分數，V(i-1,j)+σ(S[i],-)表示 S 序列從第 1 至 i-1 位置的分 數與 T 序列第 1 至 j 位置的比對的分數加上 S 序列第 i 位置元素與空白對應的分 數後的總分；V(i,j-1)+σ(-,T[j])表示 S 序列第 1 至 i 位置的分數與 T 序列第 1 至 j-1
7
位 置 的 分 數 加 上 T 序 列 第 j 位 置 的 元 素 與 空 白 對 應 後 的 總 分 ； V(i-1,j-1)+σ(S[i],T[j])，則表示 S 序列從第 1 至 i-1 位置的分數與 T 序列第 1 至 j-1 位置比對之總分加上 S 序列第 i 位置的元素與 T 序列第 j 位置的比對分數；從三 個方向來的分數來取裡面分數最大者為 V(i,j)的現有分數，透過這樣的運算，我 們可以考慮到每個元素比對時的最佳狀況；因為是全局比對，所以在全部都完成 計分後，我們便可在 V(n,m)的表格上找到最後的分數結果，而且還可以用回溯 的方式找出比對後的序列結果，此演算法是由建立一個 n * m 的矩陣來計算，所 以時間複雜度為 O(mn)。
以下我們舉個例子，假設有兩條序列我們設它們為 S 和 T，而 S 序列為 ATGG，
T 序列為 AATCG，比對分數相同者我們給予+3 分，不同者我們給予-1 分，那這 邊我們總共要建 4 個表格，分別為 G 表格 E 表格 F 表格以及 V 表格，第一步驟 先從 G 表格開始，如圖 1.3。
圖 1.3：全域序列比對示意圖。
從圖 1.3 中我們可以得到最大比對分數為 7，我們以 V(4,5)的位置來看，要 找到 V(4,5)的分數首先我們要考慮來自三個地方的分數分別是 V(3,4) V(4,4) V(3,5)，而這三個位置的分數則是分別從 G、E、F 的表個裡面來找其最大的分數，
因為是全域序列比對，所以在比對的最後得到的分數我們可以在 V(n,m)上找到，
然後可以藉由回溯的動作來往前推下一個位置，直到 V(0,0)為止，藉由回溯的動 作做我們可以得到最後序列比對出來的最佳結果，從圖 1.3 回溯的動作我們可以 得到序列比對的結果為下圖所示
9
圖 1.4：序列比對最後結果範例。
序列長度分別為 n 以及 m，在計算 V(i,j)上的的分的時候，必頇考慮來自三 個方向而來的分數，分別為 V(i-1,j)、V(i,j-1)以及 V(i-1,j-1)，然後從這三個分數 裡面挑分數最高的得分來放在 V(i,j)，而最後的得分結果會在 V(n,m)上得到，然 後藉由回溯的動作我們可以從 V(n,m)追蹤回 V(0,0)，因為此方法是建立在 (n+1)*(m-1)的陣列上，所以它的時間複雜度是 O(nm)。
1.2 間隔處罰函數
前面我們提到經由生物序列比對來找到兩個序列之間的相似性，在演化的過 程中，序列可能只有一個元素會發生變異，也有可能會一次有多個元素同時發生 變異，也就是可能會發生多次的替代、插入、或是刪除的動作，所以為了要找出 最有可能的演化過程以及應付可能發生的多次變異的這種情況，所以才會有間隔 處罰函數的產生以及後續的演算法被定義出來[12]。
關於間隔處罰函數的研究與應用有很多，但是都是基於基本的三種函數之衍 伸應用的，此三種函數分別為：任意間隔處罰函數（General gap penalties function）、
線性間隔處罰函數（Affine gap penalties function）、以及击形間隔處罰函數
（Convex gap penalties function）[13-15]。
圖 1.5：兩條序列插入間隔後下去比對的情況。
11
圖 1.6：間隔處罰的三種函數。
1.3 任意間隔處罰函數
假設 n、m 為兩序列長度，n ≧ m。針對生物序列比對上，我們針對任意間 隔處罰函數的問題，來作討論，首先，我們會先定義一些變數：
S 分數= f(間隔長度)
G(i,j) = S[i]碰 T[j]的比對分數= V(i-1,j-1) + s(S[i],T[j]) F(i,j) = S[i]碰空白的分數= max( F(i-1,j)-s , V(i-1,j)-g-s ) E(i,j) =空白碰 T[j]的分數= max( E(i,j-1)-s , V(i,j-1)-g-s ) V(i,j) = max(G(i,j), F(i,j), E(i,j))
針對任意間隔處罰函數首先我們會先定義以下幾個符號：
符號定義：
V(i,j)：第一個序列的前 i 個元素和第二個序列的前 j 個元素的最大的序列比 對分數。
G(i,j)：在 V(i,j)裡的第一個序列的第 i 個元素和第二個序列的第 j 個元素的 比對分數。
E(i,j)：在 V(i,j)裡的第一個序列的間隔和第二個序列的第 j 個元素的比對分 數。
F(i,j)：在 V(i,j)裡的第一個序列的第 i 個元素和第二個序列的間隔的比對分 數。
S(i,j)：第一個序列的第 i 個元素和第二個序列的第 j 個元素的比對分數。
g(k)：間隔處罰函數，在線性間隔處罰的時候函數為 g(k)=A-Bk，但是在任 意間隔處罰以及击型間隔處罰函數的時候間隔處罰函數為
g(k)=-A-BFuunction(k)，A 代表的是間隔的扣分，而 B 代表的則是每一 個空白的扣分，Fuunction 則是代表我們使用的函數。
13
首先我們總共會有 4 個二微陣列，分別是 V(i,j)，G(i,j)，E(i,j)，F(i,j)，首先 我們從 G(i,j)開始，G(i,j)的第一個元素我們先從 V(i-1,j-1)的分數開始做加總，然 後將第 S[i]和 T[j]個位置去做比對，比對相同給予正分，比對不相同則給予負分，
計算出來的比對分數再加上前面一個位置，也就是 V(i-1,j-1)的位置的分數去加上 比對後的得分，做加總後的總得分即是現在的 G(i,j)位置的得分。
F(i,j)也就是 S[i]碰空白的比對得分總分，一個元素碰到空白的地方，也就是 比對不符合，我們會給予扣分，照字面上來看，F(i,j)是 S[i]碰空白的比對得分總 分好像會一直扣分到最後，但是實際上並不是如此，因為他除了扣分後算出來的 總分之外，他還要去和前面我們算出來的 G(i,j)位置的總分去比較，然後選擇兩 者中最大的分數來當 F(i,j)當下的分數，所以 F(i,j)的分數是從前面一個位置的分 數加上扣分函數後的總分也就是 F(i-1,j)-s 和 V(i,j-1)-g-s 的分數，兩個分數裡選 擇 最 大 的 分 數 ， 就 是 F(i,j) 位 置 的 總 分 ； E(i,j) 也 是 如 此 ， E(i,j)E(i,j-1)-s 和 V(i,j-1)-g-s 的分數，兩個分數裡選擇最大的分數，就是 E(i,j)位置的總分。
最後，三個二微陣列都好了之後，分別是 G(i,j)、F(i,j)、E(i,j)，再從裡面找 出每個位置的最大值，將他分別填入 V(i,j)的陣列裡，也就是說，V(i,j)的分數是 max(G(i,j), F(i,j), E(i,j))，V(i-1,j-1)的分數是 max(G(i-1,j-1), F(i-1,j-1), E(i-1,j-1))，
V(i+1,j+1)的分數是 max(G(i+1,j+1), F(i+1,j+1), E(i+1,j+1))，以此順序到最後一個 位置的分數。
稍微整理一下，我們可以得到以下的演算法：
V(i,j)：第一序列前 i 個字母和第二序列前 j 個字母，最高的序列比對分數 G(i,j)：在 V(i,j)中，第一序列第 i 個字母和第二序列第 j 個字母比對 E(i,j)：在 V(i,j)中，第一序列間隔和第二序列第 j 個字母比對
F(i,j)：在 V(i,j)中，第一序列第 i 個字母和第二序列間隔比對 S(i,j)：第一序列第 i 個字母和第二序列第 j 個字母比對分數 g(k)：間隔長度為 k 的比對分數
接著利用動態規劃演算法如下：
Begin
V(i,j) = max[E(i,j), F(i,j), G(i,j)], G(i,j) = V(i-1,j-1) + S(i,j),
E(i,j) = max[V(i,k) – g(j-k)], 0≦k≦j-1, F(i,j) = max[V(k,j) – g(i-k)], 0≦k≦i-1, V(i,0) = -g(i),
V(0,j) = -g(j), E(i,0) = -g(i), F(0,j) = -g(j).
End;
最後我們可以得到時間複雜度為 O(nm2) 的演算法。試著去考慮 E(i,j)矩陣中 的每一列總共需要 O(m2) 的計算時間，總共全部共有 n 個列，所以最後全部完 成總共需要 O(nm2) 的時間，而 F(i,j)也是一樣如此；而完成 G(i,j)只需要 O(nm) 的時間。所以最後求 V(m,n) 需要 O(nm2)的時間複雜度。
15
1.4 研究目的與動機
針 對 序 列 比 對 的 演 算 法 在 過 去 一 直 是 一 個 被 廣 泛 討 論 的 ， 從 一 開 始 Needleman 及 Wunsch 於 1970 年，首先針對無間隔處罰函數的序列比對，想到可 以利用動態規劃的方法來解決序列比對繁瑣的時間問題，並提出了 O(nm)的演算 法，其中 n 和 m 分別為序列的長度，且 n≧m。
接著於 1976 年的時候，Waterman、Smith 以及 Beyer，提出了考慮任意間隔 處罰函數，其時間複雜度為 O(nm2) 的演算法[16]。爾後 Gotoh 於 1982 年針對線 性函數的特性，提出時間複雜度為 O(nm)的線性間隔處罰函數的演算法。然後 Waterman 於 1984 年，提出了击形間隔處罰函數，並預測其可能有存在著時間複 雜度為 O(nmlogm)的特性[17]。隨後就在 1988 年，Myers 及 Miller，針對击形間 隔處罰函數，首先提出了時間複雜度為 O(nmlogm)的方法[18]。此方法算是一個 相當大的突破性，因為後期開始很多的演算法都是以他們的方法下去做改良，所 以像是後面在 1989 年的的時候 Galil 及 Giancarlo，以及後面還有 Gusfield，也 都針對击形間隔處罰函數分別提出了時間複雜度是 O(nmlogm)的改良演算法 [19]。
雖然如此但是我們並未看到有相關學者針對序列比對的最大比對分數下去 做關於凹形間隔處罰(concave gap penaltines) 函數的相關研究，因為基於符合生 物意義層面的問題，我們覺得生物在演化上也有可能會反向操作來做演化，所以 既然存在著最大比對分數的击形函數，想必可能也存在著最大比對分數的凹形函 數。
所以我們就利用最後 Gusfield 在 1997 年提出的時間複雜度是 O(nmlogm)的 击形間隔處罰函數下去做探討和研究，將此函數的特性分析後研究，試著找出凹 形間隔處罰函數且時間複雜度也不會增加，也就是函數的時間複雜度也控制在 O(nmlogm)的時間裡。讓整個在序列比對上會更符合貼切生物意義。
圖 1.7：目前為止的研究發展。
17
第二章相關研究
線性間隔處罰函數是目前生物學家最用來對準生物模型中的氨基酸序列的 方式，但是若基於符合生物意義上來說的話，击形間隔處罰函數會比線性間隔處 罰函數更來的是和處理序列比對的問題，因為線性間隔處罰函數若是在面對一大 段的長劍隔的時候，會做出過重的懲罰，例如 cDNA 在進行比對的時候因為序列 之間會有許多的長間隔，此時若是用線性函數下去處理的話，會發生懲罰過重的 情形，會影響到結果的正確性和準確率。
所以在 1984 年的時候，Waterman 提出了击形間隔函數：g(k+1)+g(k)≧
g(k)+g(k-1)，其中 k 為間隔的長度，這比線性間隔處罰函數來的好用許多，因為 它並不是只考慮單一一個間隔，而是會將碰到一連串長間隔的情況也考慮進去，
實用性比線性的來的高，且在當時 Waterman 還預測它的時間複雜度可以降至 O(nmlogm)甚至有可能存在著 O(nm)的可能性。
2.1 击形間隔處罰函數研究
Needleman 及 Wunsch 在 1970 年提出了利用動態規劃來處理序列比對的問 題，其中 n 和 m 為序列長度且 n≧m，因為兩條序列長度不一定相同，也因此後 面才發展出間隔處罰的這個概念，針對間隔處罰較完整的限制和定義，是由 Waterman 於 1984 年提出，並且還預測有存在著 O(nmlogm)的時間複雜度，Myers 及 Miller 於 1988 年，針對击形間隔處罰函數，首先提出時間複雜度為 O(nmlogm) 的方法。而後 Gusfield 於 1997 年針對了击形間隔處罰函數(Convex gap Penaltines ) 提出了改良演算法並將時間複雜度降至 O(nmlogm)[20]。
圖 2.1：击形函數特性。
從圖 2.1 我們可以發現，在 d 的條件被固定的時候，且 1≦d，此時 g(q'+d)-g(q')
≦g(q+d)-g(q)，也就是說間隔長度如果越長的話，處罰分數也會隨之的遞減。
今天我們假設兩條序列長度分別為 n、m，且 n≧m，然後利用以下的動態劃 來處理击形間隔處罰函數。
19
V(i,j)=max[E(i,j) , F(i,j) , G(i,j)] , G(i,j)=V(i-1,j-1)+S(i,j) ,
E(i,j)=max[V(i,k)-g(j-k)] , 0≦k≦j-1 F(i,j)=max[V(k,j)-g(i-k)] , 0≦k≦i-1 V(i,0) = -w(i),
V(0,j) = -w(j), E(i,0) = -w(i), F(0,j) = -w(j).
處理击形間隔處罰函數我們用的是跟任意間隔函數一樣的動態規劃，我們去 分析演算法的過程中可以發現，針對每一列來看，V(i,j)、G(i,j)的時間複雜度為 O(m)，假如使用的是任意間隔處罰函數的話，那處理 E(i,j)，F(i,j)的時間複雜度 為 O(m2)及 O(n2)。
針對此章節我們首先定義一些符號並在下一章節開始探討：
符號定義
V(i,j)：第一個序列的前 i 個元素和第二個序列的前 j 個元素的最大的序列比 對分數。
G(i,j)：在 V(i,j)裡的第一個序列的第 i 個元素和第二個序列的第 j 個元素的 比對分數。
E(i,j)：在 V(i,j)裡的第一個序列的間隔和第二個序列的第 j 個元素的比對分 數。
F(i,j)：在 V(i,j)裡的第一個序列的第 i 個元素和第二個序列的間隔的比對分 數。
S(i,j)：第一個序列的第 i 個元素和第二個序列的第 j 個元素的比對分數。
g(k)：間隔處罰函數，在線性間隔處罰的時候函數為 g(k)=A-Bk，但是在任 意間隔處罰以及击型間隔處罰函數的時候間隔處罰函數為
g(k)=-A-BFuunction(k)，A 代表的是間隔的扣分，而 B 代表的則是每一 個空白的扣分，Fuunction 則是代表我們使用的函數。
b(j)：儲存地 j 個位置屬於的那個區塊。
C(i,j)：候選人表格，同於 V(k)-g(j-k)。
E'(j)：固定 E 矩陣中的其中一列，儲存第 j 個位置暫時的最大值。
E(j)：固定 E 矩陣中的其中一列，儲存第 j 個位置的最大值。
End-of-block：儲存每個區塊的最右邊的位置。
Pointer list：等同於 b(j)值。
p：已知條件 C(i,p)>C(bs,p)，p 值為符合條件的最大值。
21
2.2 Gusfield 的击形間隔處罰函數
今天我們假設兩條序列長度分別為 n、m，且 n≧m，假設 S(i,j)位置的比對 結果得分，g(k)則是代表一個長度為 k 的許多間隔串在一起的空白條，然後我們 今天以击形間隔處罰函數具有的 g(k+1)-(k)≦g(k)-g(k-1)的特性下去看假設 V(i,j) 是在(i,j)位置上的最佳比對得分數，G(i,j)是在(i,j)位置上的比對分數，比對結果 若是相同則給予正分，布相同者則給予負分，而 F(i,j)則是在上述的過程中的第 一條序列插入間隔然後與第二條序列比對後所得到的分數，反之 E(i,j)則是在上 述的過程中的第二條序列插入間隔然後與第一條序列比對後所得到的分數，然後 利用以下的動態規劃來處理击形間隔處罰函數。
圖 2.2：利用動態規劃處理击形間隔處罰函數的符號設定。
先只看 E(i,j)的地方，假設 E(i,j)=maxC(k,j)，而 C(k,j)=V(i,k)-g(j-k)，C 代表 的意思是 candidate，也就是候選者，Gusfield 今天就是利用了候選者表單去修改
了 E(i,j)，因為修改 F(i,j)同於修改 E(i,j)的方式，所以今天只需要看如何從 E(i,j) 上去做修改。
我們今天先假設幾個變數，k、j、j'、j"，且 k<j<j'<j"，而這四個點皆在相同 的列裡的四個位置上，假如 C(j,j')≦C(k,j')，那也就表示 C(j,j")≦C(k,j")，如圖 2.3 所示。
圖 2.3：k、j、j'、j"四點皆在相同列裡的四個位置上，且在候選者表單裡的取捨 關係。
根據圖 2.3 的關係我們可以得到定理 2.1 如下：
[定理 2.1]：設 k<j<j'<j"且四點皆在相同列裡的四個位置，若 C(j,j')≦C(k,j')，
則 C(j,j")≦C(k,j")。
證明：當 C(k,j')＞C(j,j')，也就是 V(i,k)-g(j'-k)≧V(i,j)-g(j'-j)，整理一下算式 後 我 們 可 以 得 到 V(i,k)- V(i,j) ≧ g(j'-k)- g(j'-j) 。 (j'-k)=(j'-j)+(j-k) ， 所 以 (j"-k)=(j"-j)+(j-k)，所以(j'-k)＜(j"-k)。現在我們先設定三個變數分別為 q=(j'-j)、
q'=(j"-j)、和 d=(j-k)，j'＜j"以及 q＜q'。根據击形間隔處罰函數 g(q'+d)-g(q')≦
g(q+d)-g(q)，d≧1，因此我們得知 g(j'-k)- g(j'-j)≧g(j"-k)- g(j"-j)，所以根據前面提 到 的 證 明 我 們 可 以 得 知 V(i,k)-V(i,j) ≧ g(j"-k)-g(j"-j) 以 及 V(i,k)-g(j"-k) ≧ V(i,j)-g(j"-j)。
23
指標和列分區
根據來自動態規劃演算法的方式，他會召回在每一個小區中 j'的變量。這個 變量是指在 k＜j 裡一個指向最左邊的細胞，它會挑選一個小的區域內貢獻最大 的變數的最佳人選“利用這個方式將整個區塊合併成一個區塊最後達到加速的效 果。
舉例來說，當 C(j,j')≦C(k,j')，那 C(j,j")≦C(k,j")，所以,j"＞,j'。假設 h(代表 位於 k+1 至 m 的任何位置)符合 E'(h)＜C(k,h)，則將 k 存入 b 矩陣中 k+1 至 m 的 位置，可知在 k+1 至 m 的範圍內間隔長度為 h-k，可以得到 h 位置暫時的最大值，
根據這個特性可以得到一個 b 矩陣如圖 2.4 所示，圖中矩陣內的每一個數字都會 形成一個小方塊區域也就是一個分割(Partition)，同一個分割內的數字都是相同的，
假如 k=9 的情形下，那在第一個分割內的內容皆為 9，從這邊可以發現 b 矩陣具 有 b(j')≧b(j'+1)的特性。
圖 2.4：b 矩陣示意圖，內部文字為矩陣內容，矩陣內的內容會利用分割(Partition) 來分成一區一區的情形，同一個分割區域內的數字皆是相同的，從圖裡我們可以 發現 b 矩陣內具有 b(j')≧b(j'+1)的特性。
[定理 2.2]：考慮各點，但在此之前 j 會對於 j'從 j+1 至 m-1 發送轉發任何候選人
的值，其為 b(j')≧b(j'+1)。
證明：設 b(j')=k 和 b(j'+1)=k'，C(k',j)≦C(k,j')，若 k＜k'，根據定理 2.1 可以得到 C(k,j'+1)≧C(k',j'+1)，所以 b(j'+1)是在 k 的集合裡面而不是在 k'裡面，因此我們
可以得到 k≧k'。
再根據定理 2.2，我們總共會有以下的三種情況，如圖 2.5 所示
圖 2.5：根據定理 2.1 可能遇到的三種情況。
根據圖 2.5 我們可以發現，b 矩陣要更新的時候，會先找到符合 E'(ps)≧(j',ps) 裡的索引 s，由圖 2.6 裡我們可以看到，當 s=4 的時候，p3及 p4間會有新的分割(見 圖 2.5 裡的情況 2)，這邊我們設為 pi，從這邊我們可以看到分割後，pi的右邊區 塊保持不變，但是左邊的部分會與 pi結合成一個區塊分割。
25
圖 2.6：b 矩陣更新時的變化。
根據 Gusfield 在 1997 年的時候提出的 O(nmlogm)演算法，針對每一列利用 反前進動態規劃(Reversed Forward Dynamic Programming)的演算法，該演算法如 下：
Initialize the end-of-block list to contain the single number m.
Initialize the associated pointer list to contain the single number() For j = 1 to m do
begin
Set k to be the first pointer on the b-pointer list.
E( j ) = Cand( k , j );
V(j) = max[G( j ) , E( j ) , F( j )];
{As before we assume that the needed F and G values have been computed.) {Now see how j's candidates change the block-partition.}
Set j' equal to the first entry on the end-of-block list.
{look for the first index s in the end-of-block list where j loses)
If Cand(b( j’), j+1) < Cand( j , j+1) then {j's candidate wins one) begin
While
The end-of-block list is not empty and Cand(b( j’), j’) < Cand(j, j’) do
begin
remove the first entry on the end-of-block list, and remove the corresponding b-pointer If the end-of-block list is not empty then
set j’to the new first entry on the end-of-block list.
end;
end {while};
If the end-of-block list is empty then place m at the head of that list;
Else {when the end-of-block list is not empty) begin
Let ps denote the first end-of-block entry.
Using binary search over the cells in blocks , find the right-most point p in that block such that Cand( j , p )> Cand( bs , p).
Add p to the head of the end-of-block list;
end;
Add j to the head of the b pointer list.
end;
end.
上述的演算法裡 Gusfield 已證明 while loop 部分執行的時間，整個列只需 要 O(m)的時間複雜度。利用二元搜尋法尋找 p 值需要 O(logm)的時間，因為總共 有 m 個搜尋，所以需要 O(mlogm)。因為每一列總共是 O(m)+ O(mlogm)且總共 有 n 個列，所以時間複雜度是 O(nmlogm)。
27
2.3 击函數演算法兩大步驟說明
首先步驟一根據演算法的敘述，一開始的時候，C(k,j)=V(k)-w(k,j)。原本的 E 表格是考慮整個直排的時間要 n2，E(j)=Max Cand(k,j)，0≦k≦j-1，如圖 2.7。
圖 2.7：击函數的 E 表格。
今天先假設一個Ē(j)的暫存器，Ē(j)是以一整個列下去看，然後做下列判斷，
假如Ē(j)<cand(j,j’)，那 Ē(j)的值就放入 cand(j,j’)的值，如圖所示。
圖 2.8：击函數的 Ē 表格。
接著步驟二現在改成由一個列來看，然後根據圖 2.5 的三個情況來作 Ē(j)的 值的判定，如圖 2.9 所示，當 b points 值等於 2 的時候，因為在最左邊和最右邊 都找不到 2 的值，所以整條的分數都沒有任何變化，也就是圖 2.5 的第一個情況，
b points 等於 3 的時候，因為 3 的值比最左邊的值小，比最右邊的值大，此時就 用 binary searach 作搜尋動作，找到符合比左邊的值小，比右邊的值大的這個特 性的 b points 值的地方，如圖 2.9 裡箭頭指向的地方為 binary searach 找到的比左 小比右大的特性的地方，然後從該位置開始一直到最左邊的值全改成現在的 b points 值也就是 3，然後右邊的部分的值則保持不動。
圖 2.9：箭頭指示的地區為使用 binary searach 搜尋到 b 值的地方。
因為用 binary searach 作搜尋 p 值的動作，時間為 O(logm)，演算法有 m 迴 圈，所以共需 O(mlogm)時間。因此最後求出所有 n 個列的 E(i)及最後序列比對 分數 V(n, m)，共需 O(nmlogm)時間。
29
第三章我們的演算法
3.1 凹形間隔處罰函數研究
生物序列比對的問題在生物資訊學中是一個非常基礎，而且一直被廣泛討論 的問題，以最大比對分數來說，針對線性間隔處罰函數、任意間隔處罰函數、击 型間隔處罰函數，都有學者相繼提出相關演算法。
但是到目前為止，都尚未有學者針對凹形間隔處罰(concave gap penaltines) 函數，提出相關研究。本篇論文主要是討論考慮凹形間隔處罰函數之序列比對的 問題，利用上一章提到的 Gusfield 針對击形間隔處罰所設計的演算法，分析其特 性並加以改良，並進一步推導出有效率的時間複雜度 O(nmlogm)的凹形間隔處罰 演算法。
圖 3.1：凹型函數特性。
從圖 3.1 我們可以發現凹型函數的特性，在 d 的條件被固定的時候，且 1≦d，
此時 g(q'+d)-g(q')≧g(q+d)-g(q)，也就是說間隔長度如果越長的話，處罰分數也會
隨之的遞增。
今天我們假設兩條序列長度分別為 n、m，且 n≧m，然後利用以下的動態劃 來處理凹形間隔處罰函數。
V(i,j)=max[E(i,j) , F(i,j) , G(i,j)] , G(i,j)=V(i-1,j-1)+S(i,j) ,
E(i,j)=max[V(i,k)-g(j-k)] , 0≦k≦j-1 F(i,j)=max[V(k,j)-g(i-k)] , 0≦k≦i-1 H(j) = Cand(0, j);
V(i,0) = -w(i), V(0,j) = -w(j), E(i,0) = -w(i), F(0,j) = -w(j).
在處理凹形間隔處罰函數我們用的是跟任意間隔函數一樣的動態規劃，針對 凹形間隔處罰函數我們首先定義一些符號：
符號定義
V(i,j)：第一個序列的前 i 個元素和第二個序列的前 j 個元素的最大的序列比 對分數。
G(i,j)：在 V(i,j)裡的第一個序列的第 i 個元素和第二個序列的第 j 個元素的 比對分數。
E(i,j)：在 V(i,j)裡的第一個序列的間隔和第二個序列的第 j 個元素的比對分 數。
F(i,j)：在 V(i,j)裡的第一個序列的第 i 個元素和第二個序列的間隔的比對分 數。
S(i,j)：第一個序列的第 i 個元素和第二個序列的第 j 個元素的比對分數。
C(i,j)：候選人表格，同於 V(k)-g(j-k)。
31
H(j)：表示演算法目前各個 E(j)的值。
g(k)：間隔處罰函數，在線性間隔處罰的時候函數為 g(k)=A-Bk，但是在任 意間隔處罰以及击型間隔處罰函數的時候間隔處罰函數為
g(k)=-A-BFuunction(k)，A 代表的是間隔的扣分，而 B 代表的則是每一 個空白的扣分，Fuunction 則是代表我們使用的函數。
b(j)：儲存地 j 個位置屬於的那個區塊。
E'(j)：固定 E 矩陣中的其中一列，儲存第 j 個位置暫時的最大值。
E(j)：固定 E 矩陣中的其中一列，儲存第 j 個位置的最大值。
End-of-block：儲存每個區塊的最右邊的位置。
Pointer list：等同於 b(j)值。
p：已知條件 C(i,p)>C(bs,p)，p 值為符合條件的最大值。
3.2 我們的演算法
我們從 E(i,j)的地方來看，利用候選者表單，我們假設 E(i,j)=maxC(k,j)，而 C(k,j)=V(i,k)-g(j-k)，接著就看看如何從 E(i,j)上去做修改。
我們今天先假設幾個變數，k、j、j'、j"，且 k<j<j'<j"，而這四個點皆在相同 的列裡的四個位置上，假如 C(j,j')≧C(k,j')，那也就表示 C(j,j")≧C(k,j")，如圖 3.2 所示。
圖 3.2：在凹型間隔處罰函數裡 k < j < j’< j”。
根據圖 3.2 裡 k、j、j’、j”的關係我們得到定理 3.1 如下：
[定理 3.1]設 k<j<j'<j"且四點皆在相同列裡的四個位置，若 C(j,j')≧C(k,j')，則 C(j,j")
≧C(k,j")。
證明：
Let d = j – k q = j’– j
33
q’= j”- j so q + d = j’– k q’+ d = j”– k
Concave: w(q + d) - w(q) ≦ w(q’+ d) - w(q’)
w(j’– k) - w(j’– j) ≦ w(j”– k) - w(j”– j) 〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃(1) if Cand(j, j’) ≧ Cand(k, j’)
V(j) - w(j’– j) ≧ V(k) - w(j’– k)
V(k) - V(j) ≦ w(j’– k) - w(j’– j) 〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃(2) (1)+(2):V(k) - V(j) ≦ w(j”– k) - w(j”– j)
V(k) - w(j”– k) ≦ V(j) - w(j”– j) Then Cand (k , j”) ≦ Cand (j , j”)
我們可以看到當 C(j,j') ≧C(k,j')，那 C(j,j")≧C(k,j")，所以, j'＞, j"。假設 h(代 表位於 k+1 至 m 的任何位置)符合 E'(h)>C(k,h)，則將 k 存入 b 矩陣中 k+1 至 m 的位置，在 k+1 至 m 的範圍內間隔長度為 h-k，可以得到 h 位置暫時的最大值，
根據這個特性我們可以得到凹函數的 b 矩陣圖，從凹函數的 b 矩陣裡可以發現它 具有 b(j'+1)≧b(j')的特性。
圖 3.3：凹函數的 b 矩陣示意圖，內部文字為矩陣內容，矩陣內的內容會利用分 割分成一區一區的情形，同一個分割區域內的數字皆為相同的，從圖裡我們可以 發現 b 矩陣內具有 b(j'+1)≧b(j')的特性。
[定理 3.2]考慮各點，但在此之前 j 會對於 j'從 j+1 至 m-1 發送轉發任何候選人的
值，其為 b(j'+1)≧b(j')。
證明：設 b(j')=k 和 b(j'+1)=k'，C(k',j)≧C(k,j')，若 k＞k'，根據定理 3.1 可以得到 C(k,j'+1)≦C(k',j'+1)，所以 b(j'+1)是在 k'的集合裡面而不是在 k 裡面，因此我們 可以得到 k'≧k。
再根據定理 3.2 我們可以得到在分區的時候會碰到的三種可能情況如圖 3.4 所 示：
35
圖 3.4：分區的時候會碰到的三種情形。
圖 3.4 就是在說 b 矩陣內的內容是如何的變化，總共會有三個情況分別是(a) 保持不變，維持原本的內容，(b)k 到 1 的地方都改成新的值而 k 到 m 的地方的 職責保持不動，(c)m 以後的全部都改新的值，以下用圖 3.5 舉個簡單的例子簡單 說明(a)(b)(c)的三種情況：
圖 3.5：矩陣內內容變化的範例。
根據前面圖 2.5 所提到的候選者的挑選方法，在圖 3.5 這裡我們是向後搜尋整個 列的值，來搜尋Ē 的值。
針 對 每 一 列 利 用 倒 退 動 態 規 劃 演 算 法 (Reversed Backward Dynamic Programming)，該演算法如下：
Initialize the end-of-block list to contain the single number m.
Initialize the associated pointer list to contain the single number 0.
For j = 1 to m do begin
Set k to be the first pointer on the b-pointer list.
E(j) = Cand(k, j);
V(j) = max[G(j), E(j), F(j)];
{As before we assume that the needed F and G values have been computed.}
{Now see how j's candidates change the block-partition.}
Set j’ equal to the first entry on the end-of-block list.
{look for the first index s in the end-of-block list where j loses}
If Cand(b(j’), j+1) > Cand(j, j+1) then { j’s candidate wins one}
begin
While
The end-of-block list is not empty and Cand(b( j’), j’) > Cand(j, j’) do begin
remove the first entry on the end-of-block list, and remove the corresponding b-pointer. If the end-of-block list is not empty then set j’ to the new first entry on the end-of-block list.
end;
end {while};
37
If the end-of-block list is empty then place m at the head of that list;
Else {when the end-of-block list is not empty}
begin
Let ps denote the first end-of-block entry. Using binary search over the cells in block s, find the right-most point p in that block such that Cand(j, p) <
Cand(bs,p).Add p to the head of the end-of-block list;
end;
Add j to the head of the b pointer list.
end;
end.
上述演算法針對固定列，利用二元搜尋法找出符合 Cand(j, p) < Cand(bs, p) 條件 p 的位置，需要 O(logm)的計算時間。演算法有 m 迴圈，所以共需 O(mlogm) 時間。因此最後求出所有 n 個列的 E(i)及最後序列比對分數 V(n, m)，共需 O(nmlogm)時間。
3.3 凹函數演算法兩大步驟說明
首先步驟一根據演算法的敘述，原本的 E 表格是以考慮整個直排的方式去找 最大值，時間要 n2，E(j)=Max Cand(k,j)，0≦k≦j-1，如圖 3.6。
圖 3.6：凹函數的 E 表格。
先假設一個 Ē(j)的暫存器，Ē(j)是以一整個列下去看，然後做下列判斷，假 如Ē(j)<cand(j,j’)，那 Ē(j)的值就放入 cand(j,j’)的值，如圖 3.7 所示。
圖 3.7：凹函數的 Ē 表格。
39
接著步驟二現在改成由一個列來看，然後根據圖 3.4 的三個情況來作 Ē(j)的 值的判定，如圖 3.8 所示，當 b points 值等於 2 的時候，因為在最右邊和最左邊 都找不到 2 的值，所以整條的分數都沒有任何變化，也就是圖 3.4 的第一個情況，
b points 等於 3 的時候，因為 3 的值比最右邊的值小，比最左邊的值大，此時就 用 binary searach 作搜尋動作，找到符合比右邊的值小，比左邊的值大的這個特 性的 b points 值的位置，如圖 3.8 裡箭頭指向的地方為 binary searach 找到的比左 小比右大的特性的地方，然後從該位置開始一直到最右邊的值全改成現在的 b points 值也就是 3，然後左邊的部分的值則保持不動。
圖 3.8：箭頭指示的地區為使用 binary searach 搜尋到 b 值的地方。
這時候我們就會得到圖 3.3 的 b points 的分區特性，因為用 binary searach 作 搜尋 p 值的動作，時間為 O(logm)，演算法有 m 迴圈，所以共需 O(mlogm)時間。
因此最後求出所有 n 個列的 E(i)及最後序列比對分數 V(n, m)，共需 O(nmlogm) 時間。
第四章結論
4.1 研究成果
生物序列的比對問題，其中 n、m 為序列長度，n≧m，不管是在生物資訊學、
生物學或是計算機數學中，都是一個非常基礎及重要的開放未解的研究課題。首 先由 Needleman 及 Wunsch 於 1970 年，針對無間隔處罰函數提出 O(nm)演算法。
接著 Waterman、Smith 及 Beyer 於 1976 年，提出考慮任意間隔處罰函數，時間 複雜度 O(nm2) 的演算法。Gotoh 於 1982 年針對線性間隔處罰，提出時間複雜度 是 O(nm)的演算法。針對击形間隔處罰函數，Myers 及 Miller 於 1988 年、Galil 及 Giancarlo 於 1989 年，以及 Gusfield 於 1997 年的時候，也分別提出了時間複雜 度 O(nmlogm)的演算法。
然而到目前為止，生物序列比對問題，尚未有學者針對凹形間隔處罰函數，
提出相關研究。我們主要針對是針對考慮 Gusfield 於 1997 年提出的 O(nmlogm) 演算法研究分析其特性後並加以改良進一步推導出凹形間隔處罰函數之生物序 列比對的演算法，而且其時間複雜度一樣為 O(nmlogm)。
41
4.2 未來研究目標
我們首先針對 Gusfield 於 1997 年的時候所提出的 O(nmlogm)的击形間隔處 罰函數演算法，研究分析其特性之後並加以改良，進一步的推導出凹形間隔處罰 函數的演算法，其時間複雜度一樣是 O(nmlogm)的時間。
雖然如此，但是這樣還是不夠的，處理生物序列比對的問題是很繁瑣冗長的，
在很久以前的生物學家只能靠人力還進行比對工作，其耗費的時間耗人力動輒就 是好幾年以上的時間，所以後續就開始有學者提出了演算法來處理序列比對的問 題。
雖然後面陸續有學者針對如何降低演算法的時間複雜度陸續提出改良形演 算法，但是基於考慮到所謂的生物特性以及生物意義的層面來說，我們更希望的 是在未來的時候能更進一步的推出所謂的增函數(Increase Function)，變成函數不 只擁有击函數的特性，甚至還有包含了凹函數以及線性函數，三個函數交替而成 的一個增函數，並且將時間複雜度也控制在 O(nmlogm)甚至在更未來的時候能推 出一個函數且時間複雜度能控制在更快更短的時間，來讓整個評分變得可以更符 合以及貼切生物意義。
圖 4.1：增函數示意圖。
圖 4.2：未來的研究目標。
43
參考文獻
[1] J. Palmer and L. Herbon, “Plant mitochondrial DNA evolves rapidly in structure, but slowly in sequence”, J. Mol. Evolut, Vol. 27, pp. 87-97, 1988.
[2] W. Li, Z. Gu, H. Wang, and A. Nekrutenko, “Evolutionary analyses of the human genome”,Nature, Vol. 409, pp. 847-849, 2001.
[3] S. Altschul, W. Gish, W. Miller, E. Myers, and D. Lipman, “Basic Local Alignment search Tool”, J. Mol. Biol., 215, pp.403-410, 1990.
[4] G. Gonnet, M. Cohen, and S. Benner. "Exhaustive matching of the entire protein sequence database", Sci., Vol. 256, pp. 1443-1445, 1992.
[5] S. Altschul and B. Erickson, “Optimal sequence alignment using affinegap costs”, Bull. Math. Biol., Vol. 48, pp.603-616, 1986.
[6] J. Fickett, “Fast Optimal alignment”, Nucleic Acids Res., Vol. 12, pp.175-180,1984.
[7] Z. Galil and R. Giancarlo, “Speeding up dynamic programming with applications to molecular biology”, Theor. Comp. Sci., Vol. 64, pp. 107-118, 1989.
[8] W. Pearson and W. Miller. “Dynamic programming algorithms for biological sequence comparison”,Numerical Computer Methods, Vol. 210, pp. 575-601, 1992.
[9] S. Needleman and C. Wunsch,“A general method applicable to the search for similarities in the amino acid sequence of two proteins”,J. Mol. Biol, Vol. 48, pp.
443-453, 1970.
[10] S. Karlin and S. Altschul. “Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes”, Proceedings of the National Academy of Sciences of the U.S.A., Vol. 87, pp. 2264-2268, 1990.
[11] D. Gusfield, “Algorithms on strings, Trees and sequences”, Cambridge University Press, Cambridge, 1997.
[12] M. Waterman, T. Smith, and W. Beyer, “Some biological sequence metrices”, Advances in Mathematics, Vol. 20, pp. 367-387, 1976.
[13] O. Gotoh, “An improved algorithm for matching biological sequences”, J. Mol.
Boil., Vol. 162, pp. 705-708, 1982.
[14] M. Waterman, “Efficient sequence alignment algorithms”, J. Theor. Biol., Vol.
108, pp. 333-337, 1984.
[15] W. Miller and E. Myers, “Sequence comparisons with concave weighting functions”, Bull. Math. Biol., Vol. 50, pp. 97-120, 1988.
[16] D. Gusfield, K. Balasubramanian, and D. Naor.
“Parametricoptimizationofsequencealignment”, Journal of Algorithmica, Special Issue on String Algorithms, Vol. 12, pp. 312-326, 1994.
[17] M. Waterman and M. Eggert. “A new algorithm fo best subsequence alignments with application to tRNA-rRNA comparisons”, J. Mol. Boil., Vol. 197, pp.
723-725, 1987.
[18] M. Fredman, “Algorithms for computing evolutionary similarity measures with length independent gap penalities”, Bull. Math. Biol., Vol. 46, pp. 533-566, 1984.
[19] S. Altschul and B. Erickson, “Optimal sequence alignment using affine gap costs”, Bull. Math. Biol., Vol. 48, pp. 603-616, 1986.
[20] S. Altschul and W. Gish, “Local Alignment Statistics”, Methods Enzymol. Vol.
266, pp. 460-480, 1996.
