李大鵬　鞠穎　鄒權(quán)

1秦皇島市第四醫(yī)院腫瘤內(nèi)科，河北秦皇島066000

2廈門大學(xué)計(jì)算機(jī)科學(xué)系，福建廈門361005

3天津大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，天津300072

基于Map Reduce的多序列星比對方法在腫瘤研究中的應(yīng)用△

李大鵬1鞠穎2鄒權(quán)3#

1秦皇島市第四醫(yī)院腫瘤內(nèi)科，河北秦皇島066000

2廈門大學(xué)計(jì)算機(jī)科學(xué)系，福建廈門361005

3天津大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，天津300072

序列比對是生物信息學(xué)的基礎(chǔ)，通過多條序列比對可以挖掘出生物序列中的各種重要信息。大規(guī)模的基因序列比對方法對運(yùn)算能力要求較高，基于Map Reduce框架的多序列比對方法在多序列星比對算法的基礎(chǔ)上利用分布式并行計(jì)算來處理大規(guī)模數(shù)據(jù)。實(shí)驗(yàn)結(jié)果表明：相對于單機(jī)處理方法，基于Map Reduce的序列比對方法可以更快速地處理大規(guī)模數(shù)據(jù)，并且具有良好的硬件擴(kuò)展性。本文探討了多序列比對在腫瘤研究方面的應(yīng)用前景。

生物信息學(xué)；多序列比對；星比對算法；Map Reduce；并行計(jì)算；癌癥

序列比對是通過對基因序列的比較找出序列之間的相似性和同源性，以便通過已知序列的結(jié)構(gòu)和功能信息來預(yù)測未知基因序列所具有的生物學(xué)信息。雙序列比對通過添加空格或者刪除字母使得給定的兩條序列具有最大的相似性。多序列比對同時(shí)比對多條序列，使得這組序列之間具有最大的相似性。目前已有許多序列比對軟件被開發(fā)出來，如Clustal［1］、T-coffee［2］、MAFFT［3］和HMMER［4］等。但是隨著序列數(shù)目的急劇增加，海量的數(shù)據(jù)對這些軟件的計(jì)算成本造成了嚴(yán)峻的挑戰(zhàn)，因此開發(fā)出一種能處理大規(guī)模數(shù)據(jù)并具有良好擴(kuò)展性的軟件是一項(xiàng)重要的工作。

在多序列比對算法方面，目前已取得了一定的成果?；陉P(guān)鍵字樹的DNA多序列星比對算法［5］利用Aho-Corasick搜索算法（簡稱AC自動機(jī)）進(jìn)行模式匹配，以便于選取中心序列，進(jìn)一步用星比對算法進(jìn)行序列比對。還有研究人員提出利用啟發(fā)式算法解決多序列比對問題，包括基于A*算法的啟發(fā)式算法［6］，基于模擬退火的多序列比對算法［7］，基于遺傳算法的多序列比對算法［8］，基于粒子群算法的多序列比對算法［9］。也有研究人員利用概率模型進(jìn)行多序列比對，如基于隱馬爾可夫模型的比對算法［10］。然而以上算法都無法有效處理基因組規(guī)?；驇装贄l序列以上的大數(shù)據(jù)。

隨著科學(xué)技術(shù)的不斷發(fā)展，生物信息數(shù)據(jù)，氣象監(jiān)測數(shù)據(jù)，社交網(wǎng)絡(luò)數(shù)據(jù)呈爆炸性增長，快速有效的處理海量數(shù)據(jù)已成為急需解決的問題。Google公司提出了Map Reduce并行計(jì)算框架用于處理大規(guī)模的并行計(jì)算問題。Apache基金會基于Map Reduce框架開發(fā)出了java編程平臺Hadoop［11］，可以實(shí)現(xiàn)多機(jī)并行處理大規(guī)模數(shù)據(jù)。本文利用Map Reduce框架，結(jié)合基于k-band的多序列星比對算法［12-13］，引入并行計(jì)算來處理大規(guī)模數(shù)據(jù)，實(shí)驗(yàn)結(jié)果表明：相對于單機(jī)處理程序，本文提出的方法能顯著提高多序列比對的運(yùn)算效率。

1　基于Map Reduce框架的多序列比對方法

1.1空隙罰分

由于序列比對是求解序列之間最大相似性的問題，所以需要一個(gè)判斷相似性的度量方法，于是引入了空隙罰分機(jī)制。對兩個(gè)或者多個(gè)字符串序列通過匹配相應(yīng)的字符或者插入空格從而得到序列之間的最大相似性。雖然加入空格匹配序列的方法能夠獲得較好的匹配效果，但是極可能會破壞其生物學(xué)意義?？障读P分就是用來防止無限制地引入空格，補(bǔ)償插入和缺失對序列相似性的影響。對于一個(gè)空格數(shù)為p的罰分計(jì)算公式如下：

其中a表示空位罰分（gap opening），即每插入一個(gè)空隙的罰分值，b表示空位擴(kuò)展罰分（gap extension），即每添加一個(gè)空格的罰分值，p為插入的空格數(shù)。

對于多序列比對假設(shè)有基因序列s=｛s1，s2，……sn｝，用score（x，y）表示兩條序列x和y之間的比對得分，那么這組基因序列的比對分值spscore（s）就表示為：

1.2多序列星比對算法

動態(tài)規(guī)劃算法是最常用于雙序列比對的算法，即Needleman和Wunsch提出的Neeedlemana-Wunsch算法［14］。由于動態(tài)規(guī)劃的計(jì)算時(shí)間復(fù)雜度為O（n^2），但是待比對的序列差異性通常不會太大，動態(tài)規(guī)劃的主要回溯路徑集中在二維表的主對角線上，為減小算法的時(shí)間復(fù)雜度以及空間復(fù)雜度，在計(jì)算動態(tài)規(guī)劃比對算法的時(shí)候可以以主對角線為中心，取寬度為k的一條窄帶上的相似性分值。對于兩條長度分度分別為n，m的不等長序列，假設(shè)m＞n，k的寬度初始時(shí)至少為m-n，此時(shí)最差的相似性分值為M（m-k-1）-2（k+1）（a+b），其中M為字符匹配的得分值，a為空隙罰分，b為空格罰分，如果計(jì)算出的相似性分值小于最差的相似性分值，就需要適當(dāng)?shù)財(cái)U(kuò)大k的數(shù)值，從而找到更好的序列比對方案。

星比對算法就是挑選一條序列作為中心序列，用這條中心序列和其他序列依次進(jìn)行比對，從而獲得比對后的更新序列以及新的中心序列。對于相似性較高的基因序列，先任意選取一條序列作為中心序列，用該中心序列依次與其他序列相互比對，然后匯總更新后的中心序列，獲得最終比對后的中心序列，再利用這條中心序列對第一次更新的序列進(jìn)行二次修正，獲得最終的更新序列組，具體計(jì)算方法如下。

星比對算法：

1.3Map Reduce的實(shí)現(xiàn)方法

由于Map Reduce每條記錄是以（key，value）形式輸入，其中key用序列名表示，value用基因序列所代表的堿基字符來表示。對于用戶輸入的文件可能會存在錯(cuò)誤，如文件格式不滿足要求，序列字符可能出現(xiàn)非法字符的情況要進(jìn)行辨別與修正。首先把文件轉(zhuǎn)化為適應(yīng)hadoop的鍵值對模型，在轉(zhuǎn)化文件結(jié)束后，提取序列組中的第一條序列并保存，作為中心序列。

然后進(jìn)入Map Reduce框架進(jìn)行序列比對計(jì)算，在第一個(gè)Map函數(shù)處理階段，數(shù)據(jù)文件會被自動分割成一個(gè)個(gè)大小為64 M的split文件，每個(gè)Map函數(shù)在不同的Data Node節(jié)點(diǎn)機(jī)上執(zhí)行，中心序列同時(shí)和每個(gè)Map塊上的序列進(jìn)行比對，每次比對會生成兩條更新序列，一條是原中心序列的更新（New_center_sequence［i］），一條是序列i的更新（New_sequence［i］），然后將比對的結(jié)果保存為（sequence_name，value）的格式，其中value的值為New_center_sequence［i］與New_sequence［i］，這兩條序列之間用tab分隔符隔開，Map_1函數(shù)輸入與輸出結(jié)果，詳見圖1。

圖1　第一階段Map函數(shù)輸入與輸出

進(jìn)入reduce階段后，reduce函數(shù)并不對數(shù)據(jù)進(jìn)行處理，直接把結(jié)果輸出到HDFS文件系統(tǒng)下。隨后從HDFS文件系統(tǒng)中將結(jié)果復(fù)制回本地文件系統(tǒng)，并從文件中分離出更新中心序列（New_center_sequence［i］）。隨后對更新中心序列進(jìn)行匯總，匯總的方法即對于n條更新序列，統(tǒng)計(jì)各條序列各個(gè)空隙所具有的空格數(shù)，求出相應(yīng)空隙所具有的最大空格數(shù)，即為最終中心序列（Final Center Sequence）的空隙插入方式，接著進(jìn)入第二階段Map Reduce，具體計(jì)算方法如下。

中心序列匯總：

第二階段Map Reduce的過程基本和第一階段相同，采用匯總后的中心序列與第一階段輸出的更新序列（New_sequence［i］）再次進(jìn)行比對。不同之處在于，Map函數(shù)并不保存比對后的中心序列，只保存二次比對序列i的最終更新結(jié)果（Final_new_sequence［i］），Map_2函數(shù)輸入與輸出，詳見圖2。

圖2　第二階段Map函數(shù)輸入與輸出

2　實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)基于Hadoop1.0版本完成，集群配置分別采用1個(gè)Name Node，2個(gè)Data Node以及1個(gè)Name Node，4個(gè)Data Node兩組對比實(shí)驗(yàn)，用于對比硬件擴(kuò)展對于運(yùn)算效率的影響。同時(shí)添加了一組單機(jī)版本的實(shí)驗(yàn)，以觀測并行計(jì)算相對于單處理在運(yùn)算效率上是否有所提高。實(shí)驗(yàn)運(yùn)行環(huán)境配置如下：單機(jī)以及集群的Name Node節(jié)點(diǎn)均采用64位Ubuntu操作系統(tǒng)，CPU為Intel（R）Core（TM）i7-3820@3.60 GHz，內(nèi)存8 G*8。集群的Data Node節(jié)點(diǎn)采用64位Ubuntu操作系統(tǒng)，CPU為Intel（R）Xeon（R）E31230@3.2 GHz，內(nèi)存4 G*8。

為便于比對不同數(shù)據(jù)規(guī)模對于運(yùn)算性能的提升，實(shí)驗(yàn)采用大小分別為1 G、2 G、4 G、8 G的人類Mitochondrial DNA數(shù)據(jù)（每條序列的長度約為16 kbp）進(jìn)行分組觀察，統(tǒng)計(jì)程序5次運(yùn)算時(shí)間并選取其中位數(shù)（圖3），從圖中可以觀測得到，在相同數(shù)據(jù)規(guī)模的情況下，基于Map Reduce框架的多序列比對程序相對于單機(jī)程序所需運(yùn)算時(shí)間降低，同時(shí)隨著數(shù)據(jù)規(guī)模的增大，程序的運(yùn)算效率進(jìn)一步提高，這表明基于Map Reduce框架的多序列比對程序能很好地處理大規(guī)模數(shù)據(jù)。

由于單機(jī)處理所需運(yùn)算時(shí)間的基數(shù)過大，圖3所示的不同配置集群在運(yùn)算效率提升上的差異性不是特別明顯，由此引入加速比來衡量并行系統(tǒng)或程序并行化性能。加速比是同一個(gè)任務(wù)在單處理器系統(tǒng)和并行處理器系統(tǒng)中運(yùn)行消耗時(shí)間的比率，計(jì)算公式如下：

圖3　不同集群配置在不同數(shù)據(jù)規(guī)模上的運(yùn)算時(shí)間

其中Sp為系統(tǒng)加速比，T1是單處理器下的運(yùn)行時(shí)間，Tp是并行系統(tǒng)中的運(yùn)行時(shí)間。統(tǒng)計(jì)不同配置集群的加速比結(jié)果（圖4），可以直觀的觀測到隨著硬件配置的擴(kuò)展，系統(tǒng)運(yùn)算效率的提升表現(xiàn)更加明顯。當(dāng)數(shù)據(jù)規(guī)模進(jìn)一步增大時(shí)，適當(dāng)?shù)臄U(kuò)充現(xiàn)有的硬件條件，可以獲取更高的運(yùn)算效率。

本文提出了一種基于Map Reduce框架，結(jié)合星比對算法的多序列比對方法，通過引入并行計(jì)算來提升程序的運(yùn)行效率。實(shí)驗(yàn)證明：隨著問題規(guī)模的擴(kuò)大，與單機(jī)單線程的多序列比對軟件相比［1，5］，基于Map Reduce的多序列比對方法能提升序列比對的運(yùn)算效率，同時(shí)程序具有良好的硬件擴(kuò)展能力。為了方便研究使用，作者還開發(fā)了可以在Hadoop平臺上使用的多序列比對軟件MSA Hadoop1.0（http://datamining.xmu.edu.cn/software/MSA Hadoop），研究者可以在上述網(wǎng)址中下載到本文實(shí)驗(yàn)中提及的軟件以及操作指南。

3　多序列比對在腫瘤研究中的應(yīng)用

多序列比對的結(jié)果有助于發(fā)現(xiàn)單核苷酸多態(tài)性（SNP），拷貝數(shù)差異（CNV）和序列缺失（indel）。而已有大量研究表明，腫瘤的產(chǎn)生機(jī)制，遺傳效應(yīng)和變異程度均與上述現(xiàn)象相關(guān)。2004年，日本研究者通過對人類的線粒體DNA進(jìn)行測序，然后利用多序列比對的方法發(fā)現(xiàn)，阿茲海默癥和糖尿病均與線粒體DNA上的若干敏感位點(diǎn)相關(guān)［15］。此外，通過多序列比對和將測序序列向基因組上回貼（mapping），也可以找出腫瘤患者與正常人在部分基因上的表達(dá)差異，相關(guān)研究發(fā)現(xiàn)MicroRNA的表達(dá)差異與腫瘤有密切關(guān)系［16-17］。還有研究發(fā)現(xiàn)，MicroRNA不但在表達(dá)量上產(chǎn)生差異，而且成熟體的序列上也有細(xì)微差異，這種差異被稱為isomiR（miRNA的多樣性表達(dá)）［18］，這些isomiR的比較分析依然需要多序列比對的方法。目前在腫瘤遺傳學(xué)研究中，盡管對SNP難以定義，但全基因組關(guān)聯(lián)分析（GWAS）依然是主要的研究方法。從腫瘤細(xì)胞到癌旁細(xì)胞，從未轉(zhuǎn)移的腫瘤細(xì)胞到轉(zhuǎn)移的腫瘤細(xì)胞的對比分析依然需要依靠GWAS，而相關(guān)的測序數(shù)據(jù)不但要映射回基因組找出差異，還需要在立體的層面對所有的位點(diǎn)進(jìn)行對齊和分析，隨著GWAS研究的深入，已不僅是對SNP位點(diǎn)進(jìn)行關(guān)聯(lián)分析［19］，還需要更多的情況考慮序列的插入和缺失，以及各種基因的表達(dá)量等，這對多序列比對提出了更高的挑戰(zhàn)。因此，關(guān)于多序列比對的理論與內(nèi)容還要面對更大規(guī)模的數(shù)據(jù)研究，以及融合入更多信息的復(fù)雜模型。

［1］ThomposonJD，GibsonTJ，HigginsD.CLUSTALW:improving the sensitivity of progressive multiple sequence alignment through sequence weighting position-specific gap penalties and weight matrix choice［J］.Nucleic Acids Res，1994，22（22）:4673-4680.

［2］Noredama C，Holm L，Higgins DG.COFFEE:an objective function for multiple sequence alignments［J］.Bioinformatics，1998，14（5）:407-422.

［3］Katoh K，Toh H.Parallelization of the MAFFT multiple sequence alignment program［J］.Bioinformatics，2010，26（15）:1899-1900.

［4］Hirosawa M，Hoshida M，Ishikawa M，et al.MASCOT: multiple alignment system for protein sequences based on three-way dynamic programming［J］.Bioinformatics，1993，9（2）:161-167.

［5］鄒權(quán)，郭茂祖，王曉凱，等.基于關(guān)鍵字樹的DNA多序列星比對算法［J］.電子學(xué)報(bào)，2009，37（8）:1746-1750.

［6］袁激光，金人超，李紅濤.基于A*算法的啟發(fā)式算法求解多序列比對問題［J］.華中科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2003，31（9）:50-52.

［7］Kim J，Pramanik S，Chung MJ.Multiple sequence alignment using simulated annealing［J］.Bioinformatics，1994，10（4）:419-426.

［8］Notredame C，O'Brien EA，Higgins DG.RAGA:RNA sequence alignment by genetic algorithm［J］.Nucleic Acids Research，1997，25（22）:4570-4580.

［9］張鵬帥，霍紅衛(wèi).多序列比對問題的粒子群優(yōu)化算法求解［J］.計(jì)算機(jī)工程與應(yīng)用，2005，18:84-87.

［10］Krogh A，Brown M，Mian IS，et al.Hidden Markov models in computational biology:Applications to protein modeling［J］.Journal of Molecular Biology，1994，235（5）: 1501-1531.

［11］Zou Q，Li XB，Jiang WR，et al.Survey of Map Reduce frame operation in bioinformatics［J］.Brief Bioinfor，2014，15（4）:637-647.

［12］Zou Q，Shan X，Jiang Y.A novel center star multiple sequence alignment algorithm based on affine gap penalty and K-band［J］.Physics Procedia，2012，33:322-327.

［13］Zou Q，Hu Q，Guo M，et al.HAlign:Fast multiple similar DNA/RNA sequence alignment based on the centre star strategy［J］.Bioinformatics，2015，31（15）:2475-2481.

［14］Needleman SB，Wunsch CD.A general method applicable to the search for similarities in the amino acid sequence of two proteins［J］.J Mol Biol，1970，48（3）:443-453.

［15］Tanaka M，Cabrera VM，González AM，et al.Mitochondrial genome variation in eastern Asia and the peopling of Japan［J］.Genome Res，2004，14（10a）:1832-1850.

［16］Zou Q，Mao Y，Hu L，et al.miRClassify:An advanced web server for miRNA family classification and annotation［J］.Comput Biol Med，2014，45:157-160.

［17］Wang Q，Wei L，Guan X，et al.Briefing in family characteristics of microRNAs and their applications in cancer research［J］.BBA-Proteins Proteom，2014，1844（1）:191-197.

［18］Guo L，Chen F.A challenge for miRNA:multiple isomiRs in miRNAomics［J］.Gene，2014，544（1）:1-7.

［19］Li P，Guo M，Wang C，et al.An overview of SNP interactions in genome-wide association studies［J］.Brief Funct Genomics，2015，14（2）:143-155.

R730

A

10.11877/j.issn.1672-1535.2016.14.06.03

2015-10-15）

國家自然科學(xué)基金（61370010）

（corresponding author），郵箱：zouquan@nclab.net