基因序列匹配算法發(fā)展及應(yīng)用

仇瑛姿

摘 要：基因序列匹配是通過(guò)匹配算法，尋找兩個(gè)或者多個(gè)序列之間的相似性和同源性，常見(jiàn)于對(duì)蛋白質(zhì)序列和DNA序列的匹配?；蛐蛄衅ヅ渌惴◤氖贾两?，在保證匹配結(jié)果精度的前提下，不斷演變出匹配時(shí)間較少的算法，以提供更優(yōu)的匹配性能。文章即是介紹基因序列匹配算法的演變和實(shí)際的應(yīng)用。

關(guān)鍵詞：基因匹配；算法演變；算法應(yīng)用

中圖分類(lèi)號(hào)：Q343.1 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：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

緒論

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

由于生物基因庫(kù)的不斷擴(kuò)充和細(xì)化，提高基因序列匹配的性能迫在眉睫?；蛐蛄衅ヅ渌惴ㄒ灿稍瓉?lái)單一的動(dòng)態(tài)規(guī)劃，全局匹配逐漸發(fā)展成啟發(fā)式的局部最優(yōu)算法。同時(shí)為了獲得更準(zhǔn)確的數(shù)據(jù)，在查詢(xún)匹配序列之前，添加了對(duì)基因序列的過(guò)濾條件，算法的不斷創(chuàng)新加之對(duì)前人算法的改良，使基因序列匹配算法逐步走向成熟。

1 全局匹配算法

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

動(dòng)態(tài)規(guī)劃算法適用于整體問(wèn)題最優(yōu)解包含子問(wèn)題最優(yōu)解和某一狀態(tài)不會(huì)被前一狀態(tài)影響，只與當(dāng)前狀態(tài)有關(guān)的場(chǎng)景。在基因序列匹配中常見(jiàn)的動(dòng)態(tài)規(guī)劃算法是：隱Markov模型（HMM）+Viterbi算法。

以下舉例利用動(dòng)態(tài)規(guī)劃算法尋找序列相似性問(wèn)題的方法（兩條DNA序列的比較）[2]：

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

Step2：通過(guò)矩陣計(jì)算出從初始位置（兩個(gè)序列的起點(diǎn)）走到矩陣下一個(gè)位置的距離（或者得分）。

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

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

隱Markov模型是一種重要的統(tǒng)計(jì)方法，在生物信息學(xué)中應(yīng)用于線(xiàn)性序列分析或基因發(fā)現(xiàn)等發(fā)main，利用它對(duì)普通動(dòng)態(tài)規(guī)劃算法進(jìn)行簡(jiǎn)化，這是由于普通的動(dòng)態(tài)規(guī)劃算法需要遍歷矩陣所有位置，對(duì)于大量的基因匹配，會(huì)耗費(fèi)相當(dāng)長(zhǎng)的時(shí)間，而HMM模型，包含隱含在兩個(gè)狀態(tài)之間的轉(zhuǎn)移概率，可以減少遍歷節(jié)點(diǎn)，加速普通動(dòng)態(tài)規(guī)劃的時(shí)間。

Viterbi算法應(yīng)用到基因序列匹配[3][4]：

從HMM模型中可以推導(dǎo)出諸如核苷酸序列ACAC出現(xiàn)的概率：

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

Viterbi算法要解決的就是HMM模型中從前t個(gè)堿基狀態(tài)到最終的堿基狀態(tài)k的觀(guān)察結(jié)果中，取最后可能對(duì)應(yīng)狀態(tài)的序列，獲得Viterbi路徑。

2 局部匹配算法

局部匹配算法的核心是利用啟發(fā)式在局部獲得最優(yōu)解，然后通過(guò)局部最優(yōu)解得到全局最優(yōu)解。這種方法精度沒(méi)有全局動(dòng)態(tài)規(guī)劃的方法高，但比之前的方法快100倍左右，犧牲了少量精度換取了執(zhí)行性能。要注意的是，動(dòng)態(tài)規(guī)劃算法其實(shí)也可以通過(guò)分治法轉(zhuǎn)變成局部動(dòng)態(tài)規(guī)劃匹配，但這種算法不在本文討論范圍之內(nèi)。常見(jiàn)的局部匹配算法有FASTA算法和BLAST算法。

FASTA算法和BLAST算法均用到了啟發(fā)式，以簡(jiǎn)化運(yùn)算步驟。啟發(fā)式即利用先驗(yàn)概率獲得預(yù)知能力，對(duì)不確定的狀態(tài)進(jìn)行判斷，在某些極端情況下，啟發(fā)式會(huì)獲得很差的結(jié)果，但通常可以在合理的時(shí)間內(nèi)獲得正確的答案。啟發(fā)式算法的結(jié)構(gòu)是鄰域搜索，即從一組初始解中產(chǎn)生若干鄰域解，從鄰域解中獲得并更新當(dāng)前的狀態(tài)，然后迭代并取滿(mǎn)足收斂條件的狀態(tài)。

FASTA算法[5]是由Lipman和Pearson推導(dǎo)實(shí)現(xiàn)的，具體實(shí)現(xiàn)步驟如下：

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

Seq1：GCATGCGCCCAC

Seq2：GAATGCGCCAAC

Seq1和Seq2的長(zhǎng)度N1=N2=12，取k=2，則每一條序列都可以得到N-k+1個(gè)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：計(jì)算單個(gè)序列中相同單詞片段的數(shù)量：

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中每一個(gè)單詞出現(xiàn)的數(shù)量?jī)蓛上喑耍玫揭粋€(gè)wordlist（Seq1存在而Seq2不存在的單詞數(shù)量為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：將所得結(jié)果求和，total=6+2+1+1+0+2+1=13

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

BLAST快速算法是另一種啟發(fā)式局部算法[6][7]，目前被廣泛應(yīng)用并延伸出多種衍變算法，諸如：BLAST2，PSI-BLAST等等。基本思想是通過(guò)words的得分矩陣，取得質(zhì)量更高的詞表，然后在該詞領(lǐng)域不斷延伸尋找高分片段，當(dāng)?shù)梅值陀谠O(shè)定值時(shí)，中止尋找流程。

算法的主要步驟是：

Step1：定長(zhǎng)分割序列，一般DNA序列的單詞長(zhǎng)度為3，蛋白質(zhì)的單詞長(zhǎng)度為11

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

Step3：兩端延伸直至分值

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

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

參考文獻(xiàn)：

[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.