基因序列匹配算法發展及應用

仇瑛姿

摘 要：基因序列匹配是通過匹配算法，尋找兩個或者多個序列之間的相似性和同源性，常見于對蛋白質序列和DNA序列的匹配。基因序列匹配算法從始至今，在保證匹配結果精度的前提下，不斷演變出匹配時間較少的算法，以提供更優的匹配性能。文章即是介紹基因序列匹配算法的演變和實際的應用。

關鍵詞：基因匹配；算法演變；算法應用

中圖分類號：Q343.1 文獻標志碼：A 文章編號：2095-2945（2018）13-0063-02

Abstract： Gene sequence matching is to find the similarity and homology between two or more sequences by matching algorithm， which is usually used to match protein sequence and DNA sequence. From the beginning to now， gene sequence matching algorithms have evolved to provide better matching performance under the premise of ensuring the accuracy of matching results. This paper introduces the evolution and practical application of gene sequence matching algorithm.

Keywords： gene matching； algorithm evolution； algorithm application

緒論

隨著生物科學和計算機科學的迅猛發展，計算機和生物相結合形成一門新學科，利用計算機對數據的快速處理能力，挖掘大量而復雜的生物數據，為生物醫學的飛躍提供了有利條件。從1990年人類基因組計劃至今，已完成了人類基因組DNA30億個堿基對的測序，破譯了人類超過95%的遺傳信息，我國也在2017年12月啟動了“中國十萬人基因組計劃”，通過繪制中國人精細的基因組圖譜，來研究疾病健康和基因遺傳的關系。世界各國對生物基因組測序工作極快地展開，而且可以預測，今后DNA測序需求的增長將更為驚人，這些海量數據的積累和運算，挖掘與分析，離不開現代的計算機技術。

由于生物基因庫的不斷擴充和細化，提高基因序列匹配的性能迫在眉睫。基因序列匹配算法也由原來單一的動態規劃，全局匹配逐漸發展成啟發式的局部最優算法。同時為了獲得更準確的數據，在查詢匹配序列之前，添加了對基因序列的過濾條件，算法的不斷創新加之對前人算法的改良，使基因序列匹配算法逐步走向成熟。

1 全局匹配算法

全局匹配即整體考慮整條基因序列相似性，一條序列轉化成另一條序列需要的最小距離決定了該序列相似性的得分值，距離越小，則得分越高，相似度也越高。較為常見的是動態規劃算法，該算法通過拆解問題，之后處理問題之間不同狀態之間的關系，使問題可以通過遞歸的方式去分解[1]。

動態規劃算法適用于整體問題最優解包含子問題最優解和某一狀態不會被前一狀態影響，只與當前狀態有關的場景。在基因序列匹配中常見的動態規劃算法是：隱Markov模型（HMM）+Viterbi算法。

以下舉例利用動態規劃算法尋找序列相似性問題的方法（兩條DNA序列的比較）[2]：

Step1：將兩條DNA序列（Seq1， Seq2）的所有堿基對生成字符串矩陣，并得出兩個序列所有堿基對應的對應關系的距離值（Seq1轉換成Seq2的得分值）。計算兩個堿基距離一般有兩種方法：一種是漢明距離，即不考慮空位，移位的情況，一一對應；另一種是編輯距離，即將Seq1的一個位置上的堿基變換成Seq2的一個位置上刪除/替換/插入的最少基本操作數目。

Step2：通過矩陣計算出從初始位置（兩個序列的起點）走到矩陣下一個位置的距離（或者得分）。

Step3：遍歷整個矩陣的位置后得到所有位置的得分。

Step4：從終點位置向前找距離最小的位置（或者得分最高的位置）。

隱Markov模型是一種重要的統計方法，在生物信息學中應用于線性序列分析或基因發現等發main，利用它對普通動態規劃算法進行簡化，這是由于普通的動態規劃算法需要遍歷矩陣所有位置，對于大量的基因匹配，會耗費相當長的時間，而HMM模型，包含隱含在兩個狀態之間的轉移概率，可以減少遍歷節點，加速普通動態規劃的時間。

Viterbi算法應用到基因序列匹配[3][4]：

從HMM模型中可以推導出諸如核苷酸序列ACAC出現的概率：

P（ACAC）=P（A1）*P（2|1）*P（C2）*P（3|2）*P（A3）*P（4|3）*P（C4）

Viterbi算法要解決的就是HMM模型中從前t個堿基狀態到最終的堿基狀態k的觀察結果中，取最后可能對應狀態的序列，獲得Viterbi路徑。

2 局部匹配算法

局部匹配算法的核心是利用啟發式在局部獲得最優解，然后通過局部最優解得到全局最優解。這種方法精度沒有全局動態規劃的方法高，但比之前的方法快100倍左右，犧牲了少量精度換取了執行性能。要注意的是，動態規劃算法其實也可以通過分治法轉變成局部動態規劃匹配，但這種算法不在本文討論范圍之內。常見的局部匹配算法有FASTA算法和BLAST算法。

FASTA算法和BLAST算法均用到了啟發式，以簡化運算步驟。啟發式即利用先驗概率獲得預知能力，對不確定的狀態進行判斷，在某些極端情況下，啟發式會獲得很差的結果，但通常可以在合理的時間內獲得正確的答案。啟發式算法的結構是鄰域搜索，即從一組初始解中產生若干鄰域解，從鄰域解中獲得并更新當前的狀態，然后迭代并取滿足收斂條件的狀態。

FASTA算法[5]是由Lipman和Pearson推導實現的，具體實現步驟如下：

Step1：從兩個序列中尋找到長度為k的共同片段（k-tuple/ktup），諸如：

Seq1：GCATGCGCCCAC

Seq2：GAATGCGCCAAC

Seq1和Seq2的長度N1=N2=12，取k=2，則每一條序列都可以得到N-k+1個words，如下：

Seq1-words：GC（1）|CA（2）|AT（3）|TG（4）|GC（5）|CG（6）|GC（7）|CC（8）|CC（9）|CA（10）|AC（11）

Seq2-words：GA（1）|AA（2）|AT（3）|TG（4）|GC（5）|CG（6）|GC（7）|CC（8）|CA（10）|AA（11）|AC（12）

Step2：計算單個序列中相同單詞片段的數量：

Seq1中相同word的位置：GC（1，5，7） CA（2，10） AT（3） TG（4） CG（6） CC（8，9） AC（11）

Seq2中相同word的位置：GA（1） AA（2，11） AT（3） TG（4） GC（5，7） CC（8） CA（10） AC（12）

Step3：隨后將Seq1和Seq2中每一個單詞出現的數量兩兩相乘，得到一個wordlist（Seq1存在而Seq2不存在的單詞數量為0）：

Common words：GC=3*2=6；CA=2*1=2；AT=1*1=1；TG=1*1=1；CG=1*0=0；CC=2*1=2；AC=1*1=1

Step4：將所得結果求和，total=6+2+1+1+0+2+1=13

FASTA算法的數學公式為：Σw=|Lw（I）|*|Lw（J）|，其中I，J表示兩個序列，Lw表示序列每個單詞出現的次數，算法的復雜度為O（N1+N2），為了增加運行速度，可以將同一條對角線中的臨近單詞合并成一個。FASTA算法對DNA序列的搜索較于蛋白質更為靈敏，由于僅尋找最佳匹配，則有可能會遺漏掉其他有意義比對。

BLAST快速算法是另一種啟發式局部算法[6][7]，目前被廣泛應用并延伸出多種衍變算法，諸如：BLAST2，PSI-BLAST等等。基本思想是通過words的得分矩陣，取得質量更高的詞表，然后在該詞領域不斷延伸尋找高分片段，當得分低于設定值時，中止尋找流程。

算法的主要步驟是：

Step1：定長分割序列，一般DNA序列的單詞長度為3，蛋白質的單詞長度為11

Step2：利用得分矩陣取出分值>值域的單詞并組成單詞表

Step3：兩端延伸直至分值

Step4：將獲得的高分片段相近合并

利用BLAST算法可以獲得多條高分片段，能夠保留更多有意義的匹配。

參考文獻：

[1]Sniedovich， M. Dijkstra's algorithm revisited： the dynamic programming connexion[J]. Journal of Control and Cybernetics，2006，35（3）：599-620.

[2]Waterman， M. S. and M. Eggert . A new algorithm for best subsequence alignments with application to tRNA-rRNA comparisons[J]. J. Mol. Biol， 1987，197：723-728.

[3]Byung-Jun Yoon.Hidden Markov Models and their Applications in Biological Sequence Analysis[J]. PubMed Central， Curr Genomics，2009，10（6）：402-415.

[4]Krogh，A. An introduction to hidden Markov models for biological sequences[J]. Computational Methods in Molecular Biology . Elsevier， Amsterdam， The Netherlands， pp，1998：46-63.

[5]Lipman， DJ； Pearson， WR. Rapid and sensitive protein similarity searche[J]. Science， 1985，227（4693）：1435-41.

[6]Camacho， C.； Coulouris， G.； Avagyan， V.； Ma， N.； Papadopoulos， J.； Bealer， K.； Madden， T. L. BLAST+： Architecture and applications[J]. BMC Bioinformatics，2009，10：421.

[7]Kent， W. James. BLAT-The BLAST-Like Alignment Tool[J]. Genome Research，2002，12（4）：656-664.