一種Ｔｉｌｉｎｇ Ａｒｒａｙ信號識別流程的并行優(yōu)化

摘要：針對一種非常有效的信號識別算法——滑窗法（sliding window，SW），在其嵌入的串行信號識別流程基礎(chǔ)上提出并行處理方案，并基于MPI（message passing interface）平臺實現(xiàn)并行優(yōu)化版本。這在基因芯片研究領(lǐng)域是個嶄新的思路。最后，根據(jù)公共數(shù)據(jù)Affymetrix（2005）Tiling Array，展示SW并行流程的運行性能。實驗結(jié)果表明，信號識別流程的并行化對大規(guī)模的數(shù)據(jù)運算非常有效且十分必要。

關(guān)鍵詞：基因芯片; Tiling Array; 信號識別流程; 滑窗法; 消息傳遞接口

中圖分類號：TP312文獻(xiàn)標(biāo)志碼：A

文章編號：1001－3695(2008)06－1666－04

隨著人類和其他物種基因組測序的完成，基因組各種數(shù)據(jù)呈指數(shù)增長。如何從海量數(shù)據(jù)中提取重要信息，闡明千變?nèi)f化的生命現(xiàn)象，是生物信息學(xué)面臨的一個挑戰(zhàn)。其中，繪制人類基因組表達(dá)圖譜是在分子水平上理解遺傳信息以及單個細(xì)胞、組織特定功能的關(guān)鍵步驟。以此為目標(biāo)，基因組Tiling Array已經(jīng)成為一種識別表達(dá)區(qū)域的流行工具。在過去四年間，一系列基于Tiling Array的重大實驗和發(fā)現(xiàn)在《Science》上發(fā)表。雖然這些實驗細(xì)節(jié)不同，如芯片制作、RNA提取、雜交環(huán)境等，但它們有一個相同的設(shè)計理念，即構(gòu)造的Tiling Array探針覆蓋物種部分染色體的全部非重復(fù)序列。

設(shè)計Tiling Array芯片的最終目的在于對所得到的芯片亮度數(shù)據(jù)進(jìn)行挖掘和分析，得到新的生物信息學(xué)知識，而信號識別就是最先要做的工作。它把芯片亮度數(shù)據(jù)轉(zhuǎn)換成相應(yīng)的表達(dá)分值，作為判斷探針表達(dá)與否的指標(biāo)。已有的生物芯片信號識別算法更關(guān)注于精度和可信度，對難以容忍的程序運行時間沒有有效的研究。因此并行思想的引入在生物芯片研究領(lǐng)域是個全新的嘗試，具有很好的指導(dǎo)意義。同時，本文所受支持項目“特征信息發(fā)現(xiàn)的并行算法及實現(xiàn)研究”和“天文、生物信息和計算化學(xué)網(wǎng)格計算應(yīng)用系統(tǒng)建設(shè)”的宗旨之一是針對Ti－ling Array技術(shù)中非編碼遺傳信息識別和序列驗證等實際應(yīng)用需求，設(shè)計有效的并行算法，并集成在網(wǎng)絡(luò)計算應(yīng)用系統(tǒng)上為用戶提供服務(wù)。本文對Tiling Array信號識別流程的并行優(yōu)化正是在上述思想的指導(dǎo)下進(jìn)行的前期研究。

1Tiling Array構(gòu)建和多種信號識別方法

1．1Tiling Array實驗平臺的構(gòu)建方法

Tiling Array的構(gòu)建是整個研究過程的第一步。目前，有兩種典型的Tiling Array構(gòu)建方法。一種是寡核苷酸oligo Tiling Array，這種工藝生產(chǎn)出來的芯片探針長度為25~60 bp；另外一種Tiling Array構(gòu)建方法使用PCR產(chǎn)物或BAC(bacterial artificial chromosome)制作，PCR芯片序列長度大約在1 KB。PCR和BAC芯片的任一探針位置都包括目標(biāo)序列及目標(biāo)序列的互補(bǔ)序列，如果沒有大量的額外實驗，將很難確定表達(dá)序列是在雙鏈的哪個鏈上。例如2003年一篇關(guān)于人類22號染色體的Tiling Array文章，它的Tiling Array實驗利用了PCR模式，其中需要的PCR產(chǎn)物要超過2萬次的PCR反應(yīng)獲得。基于此，本文選取第一種oligo Tiling Array作為研究對象。

Oligo芯片的構(gòu)建包括生成代表目標(biāo)基因區(qū)域的核酸探針，以及將這些核酸探針固定在芯片上，如圖1所示。有三種方法在基因序列上選取探針，分別為重疊、首尾相接或者間隔為預(yù)先定義距離，如圖2所示。

本文所使用的數(shù)據(jù)是2005年發(fā)表于《Science》的Affymetrix Tiling Array實驗，下文簡稱Affy(2005) Tiling Array。此實驗芯片是一種高密度寡核苷酸芯片，采用原位光刻錄技術(shù)在芯片上直接合成PM、MM模式探針。PM（perfect match）探針表示完全最佳匹配，與之對應(yīng)的MM（mismatch）探針表示非完全匹配。PM與MM探針序列基本相同，只是MM探針序列中間位置的堿基與PM探針同處堿基互補(bǔ)，PM與MM組成了一對probe－pair。一個probe－set標(biāo)志基因組序列某個小塊區(qū)域，它包括很多的probe－pair，如圖3所示。設(shè)計MM探針的目的是為了減少非特異性交叉雜交（nonspecific cross－h(huán)ybridization）的干擾。芯片設(shè)計長度為25 mer和步長5 bp的探針覆蓋了人類10條染色體的非重復(fù)區(qū)域，并與八個細(xì)胞過程（cell lines）的表達(dá)進(jìn)行雜交實驗，得到芯片亮度數(shù)據(jù)。

1．2Tiling Array信號識別算法及比較

從2001年首個Tiling Array實驗設(shè)計完成開始，不同的信號識別算法相繼出現(xiàn)。這包括Kapranov[1]、Schadt[2]、Kampa[3]、Bertone[6]和Cheng[8]等人研究的算法。下面從基因組的表達(dá)識別出發(fā)，介紹從芯片亮度信號到表達(dá)識別判斷的兩個典型算法。

1） 滑窗算法[3]該算法本身并不對芯片的信號文件（.cell）作任何初始的背景修正和標(biāo)準(zhǔn)化處理。其優(yōu)勢在于利用查詢序列探針和相鄰序列探針來確定這個查詢序列探針是否被表達(dá)（positive）。SW算法可分為兩個部分描述：

a)探針的表達(dá)評估。這一步主要通過一個探針和其相鄰探針亮度確定此探針的表達(dá)分值。首先設(shè)定一個固定的窗口長度，其半徑為BW，則窗口長度為BW×2+1，對于任意探針Pj，設(shè)Zj=PMj-MMj。對于某個觀察探針Pi，設(shè)其中心在染色體上的位置為PTi，則包含在窗口區(qū)間[PTi-BW，PTi+BW]種的所有探針兩兩之間作平均值計算Am，n=(Zm+Zn)/2(m，n為窗口內(nèi)包含的探針編號)。最后考察探針的表達(dá)評估值Ei=Pseudo－median(Pi)=median(Am，n：m≤n)為窗口內(nèi)所有兩兩探針平均值的中值。其中確定合適的滑窗長度很重要，有兩個決定因素：探針間距離，即探針序列中心間距離；探針?biāo)谌旧w外顯子exon的平均長度（大約為150 bp）。以Affy(2005) Ti－ling Array為例，芯片探針長度為25 mer，步長為5 bp，所以若窗口長度設(shè)為150 bp，窗口內(nèi)大約包含31個探針。其探針表達(dá)評估的計算過程如圖4所示。

b)表達(dá)與非表達(dá)區(qū)域的劃分。利用在芯片上的細(xì)菌寡核苷酸測定一個threshold值，判斷查詢探針表達(dá)水平Ei是否高于這個threshold值。若高于則說明探針序列為表達(dá)。接下來，在滿足長度(minrun)≥90 bp，連續(xù)探針間隔(maxgap)≤40 bp的情況下，表達(dá)的探針可被認(rèn)定是組成表達(dá)片斷的一部分。

2） 亮度分布(signal distribution，SD)方法[6]它是Bertone等人用于識別人類基因組非注釋區(qū)域的表達(dá)序列片斷TARs（transcriptional active regions）的方法。TAR和transfrags雖然稱呼不同，但是意義相同。在Bertone等人的實驗中，探針的原始亮度達(dá)到整張芯片亮度的90%以上，才被認(rèn)為是表達(dá)探針。以當(dāng)時Bertone等人的實驗芯片為例，探針長度為36 mer，步長46 bp，則要求至少五個連續(xù)的表達(dá)探針可被稱為TARs。任何不表達(dá)探針的出現(xiàn)都將中止TARs的延展。

3） 兩種算法的比較SW適用于Affymetrix公司設(shè)計的芯片，其芯片中嵌入了可作為表達(dá)參照的物種序列探針；SD方法并不拘于芯片設(shè)計。另外，兩者在信號識別效果上也各有優(yōu)劣。SD方法在判斷探針表達(dá)上有些粗糙，容易受到交叉雜交的影響從而增加假陽性率，并且淘汰表達(dá)量不高的探針；而由于SW算法相應(yīng)的芯片探針設(shè)計方式，其本身對交叉雜交有一定制約作用，而且由于細(xì)菌基因組探針的輔助，低表達(dá)探針也能一定程度地被發(fā)現(xiàn)。一些方法測試的文獻(xiàn)也指出，SW算法比SD算法性能更突出。

2SW信號識別流程的串行實現(xiàn)

Tiling Array實驗產(chǎn)生大量的數(shù)據(jù)，對這些數(shù)據(jù)重新整合格式，提取有用的信息，并采用恰當(dāng)?shù)男盘栕R別算法加以處理，直至最后得到探針表達(dá)與否的衡量指標(biāo)，是Tiling Array實驗?zāi)康乃冢彩切盘栕R別的基本流程。

通過研究Affymetrix芯片設(shè)計特點和各種信號識別算法，本文在信號識別流程中嵌入較優(yōu)的SW算法并對整個流程加以實現(xiàn)，如圖5所示。

2．1從Tiling Array數(shù)據(jù)提取出滿足計算要求的數(shù)據(jù)格式

Affymetrix（2005）平臺共構(gòu)建98張芯片，覆蓋10條人類染色體，每張芯片上包含1 638 400個探針，所得到的結(jié)果數(shù)據(jù)字段包括探針堿基序列、探針基因坐標(biāo)位置、PM探針芯片坐標(biāo)位置、MM探針芯片坐標(biāo)位置以及與八個樣本雜交產(chǎn)生的亮度值。但是最初的數(shù)據(jù)并不是以方便程序處理的格式顯示的，有用的信息分散在不同的文件夾內(nèi)。以人類22號染色體為例，其相關(guān)數(shù)據(jù)分布在第79~83號芯片上，要從這些龐雜的數(shù)據(jù)文件中把需要處理的染色體相關(guān)信息提取出來，并設(shè)計合適的數(shù)據(jù)格式，把提取出來的信息加以存儲，這一系列過程本文借助于MySQL來實現(xiàn)。為方便說明問題，本文只提取22號染色體和其中一個細(xì)胞過程A375的細(xì)胞質(zhì)polyA表達(dá)樣本（A375_cytosolic_polyA+）進(jìn)行雜交得到的Tiling Array數(shù)據(jù)。

a）從第79~83號芯片上提取出22號染色體的坐標(biāo)和亮度相關(guān)信息。根據(jù)探針上堿基序列走向把這些信息分成“f”和“t”走向兩部分（在計算偽中值時作為兩個數(shù)據(jù)輸入文件）。

b）根據(jù)數(shù)據(jù)格式設(shè)計合理的數(shù)據(jù)庫表，并將坐標(biāo)數(shù)據(jù)和亮度信息分別導(dǎo)入。

c）依據(jù)PM探針坐標(biāo)值對坐標(biāo)數(shù)據(jù)表和亮度信息表進(jìn)行連接查詢。

d）依據(jù)MM探針坐標(biāo)值對坐標(biāo)數(shù)據(jù)表和亮度信息表進(jìn)行連接查詢。

e）設(shè)計數(shù)據(jù)庫表，將d）e）的查詢結(jié)果導(dǎo)入數(shù)據(jù)庫中，得到關(guān)于PM和MM探針全部信息的兩張表。

f）依據(jù)MM探針的Y基因坐標(biāo)值比PM探針的Y基因坐標(biāo)值大一個bp單位且兩者的X基因坐標(biāo)值相等，對PM信息表和MM信息表連接查詢。

g）取這幾個亮度值的平均值作為該探針亮度值，并且只保留其中一條記錄。

步驟g）是后處理階段。抽樣觀察上述步驟的提取結(jié)果，發(fā)現(xiàn)擁有相同基因坐標(biāo)位置的相同探針序列被放置在芯片的不同坐標(biāo)位置上。這種情況下，經(jīng)SW算法處理后，該探針是否顯示就無法確定了。

h）把經(jīng)過如上處理步驟得到的數(shù)據(jù)庫表導(dǎo)出為結(jié)果數(shù)據(jù)文件。

2．2應(yīng)用SW算法處理提取出的數(shù)據(jù)

提取后的格式化數(shù)據(jù)滿足SW算法對輸入數(shù)據(jù)的要求。當(dāng)SW算法需要輸入數(shù)據(jù)時，只需按順序依次讀入數(shù)據(jù)文件，每31條記錄對應(yīng)計算出一個偽中值，直至到達(dá)文件末尾。設(shè)計串行SW程序流程如圖6所示。

3SW信號識別流程的并行優(yōu)化

3．1并行可行化分析

從SW流程描述可以看出，大部分運行時間耗費在對龐大數(shù)據(jù)文件的操作上面，導(dǎo)致有效運算時間所占比重較低，程序運行效率很低。

具體來講，對于2．1節(jié)數(shù)據(jù)提取階段，占用時間最多的是對兩個數(shù)據(jù)庫表的連接查詢。

若一個SQL查詢同時涉及兩個以上的表，則稱之為連接查詢。連接查詢中用來連接兩個表的條件稱為連接條件。從概念上講，若對兩張表進(jìn)行連接查詢，數(shù)據(jù)庫管理系統(tǒng)（DBMS）執(zhí)行連接操作的過程如下：

a）在第一張表中找到第一條記錄。

b）從頭開始掃描第二張，逐一查找滿足連接條件的記錄。

c）找到后把第一張表中的第一個記錄和第二張表的該記錄拼接起來，形成結(jié)果表中的一個記錄。

d）重復(fù)b）和c），直到第二張表中的全部記錄查找完畢。

e）再找第一張表中第二條記錄，跳到b）。

f）重復(fù)上述操作，直到第一張表中的全部記錄都處理完畢為止。

假定第一張表中記錄條數(shù)為N1，第二張表中記錄條數(shù)為N2，那么對這兩個表進(jìn)行連接查詢的時間復(fù)雜度為O（N1×N2）。如果N1和N2都很大，那么它們的乘積N1×N2相當(dāng)可觀。這會造成SW算法在數(shù)據(jù)提取階段花費大量時間。

對2．2節(jié)的數(shù)據(jù)處理階段，探針分值計算復(fù)雜度并不高。但是由于提取出來的Tiling Array格式化文件非常大，而每一個探針分值的計算都要從文件中讀取一次數(shù)據(jù)，頻繁的磁盤I/O導(dǎo)致程序運行效率較低。

Affy(2005 )Tiling Array原始數(shù)據(jù)不存在相關(guān)性，所以程序運行效率的決定性因素就是數(shù)據(jù)量。假設(shè)有n條數(shù)據(jù)記錄，如果只有一個處理器對它們進(jìn)行處理，那么工作負(fù)載為S＝n；如果有np個處理器相互協(xié)作處理這n條記錄，那么每個處理器的工作負(fù)載就降至S/np，則程序的執(zhí)行時間也隨之減少。這種做法也正和數(shù)據(jù)庫提高查詢效率的優(yōu)化手段之一——大表數(shù)據(jù)分割不謀而合。因此，可以利用數(shù)據(jù)分割的思路實現(xiàn)SW識別流程的并行化。

3．2并行優(yōu)化

數(shù)據(jù)分割是并行編程中的一個基本思想，對需要處理大數(shù)據(jù)量的算法來講是行之有效的并行方法。首先，把龐大的數(shù)據(jù)文件劃分成若干部分，并盡可能平均地分配到各個處理器上，每個處理器計算出自己應(yīng)該處理的記錄條數(shù)，并進(jìn)行處理。數(shù)據(jù)分割要盡量保證每個處理器所要處理的記錄條數(shù)和別的處理器相差最小，以此來保證各處理器運行時間的一致性，提高運行效率。

并行編程模式有三種類型：主從模式、對等模式和多程序模式。本文采用主從模式和對等模式相結(jié)合的方法對SW進(jìn)行并行優(yōu)化。在程序的開始和數(shù)據(jù)提取部分的結(jié)束，由主進(jìn)程負(fù)責(zé)程序初始參數(shù)的廣播和數(shù)據(jù)提取結(jié)果的收集工作。在數(shù)據(jù)處理部分，因為每個處理器的處理時間不盡相同，并沒有讓主進(jìn)程進(jìn)行結(jié)果收集，而是讓各個進(jìn)程分別把自己的計算結(jié)果寫入輸出文件。這樣可以避免在通信和相互等待上花費不必要的時間。SW整個流程并行程序流程如圖7所示。

每個處理器都要完成的任務(wù)是：

a）多次的部分?jǐn)?shù)據(jù)庫表和另一數(shù)據(jù)庫表的連接查詢。

b）調(diào)用SW算法對部分提取出來的數(shù)據(jù)進(jìn)行處理，并把結(jié)果寫入輸出文件。

主進(jìn)程除了完成以上的任務(wù)之外，還要進(jìn)行并行程序初始化時參數(shù)的廣播，以及部分聯(lián)合查詢結(jié)果收集和整合。并行SW整個流程中涉及到的進(jìn)程間通信只有主進(jìn)程對各個進(jìn)程的參數(shù)廣播以及數(shù)據(jù)庫查詢結(jié)果整合，其他各個進(jìn)程間并不涉及通信。并行優(yōu)化性能與主進(jìn)程執(zhí)行收集時所進(jìn)行通信的時間開銷有直接關(guān)系。若進(jìn)行數(shù)據(jù)庫查詢結(jié)果收集的時間表示為T(DB， pid)。其中：pid為進(jìn)程序號。各個進(jìn)程執(zhí)行對等模式的任務(wù)時間表示為PT(sametask， pid)。設(shè)有np個處理器，則只有當(dāng)∑np-1pid=0T(DB，pid)<<∑np-1pid=1PT(sametask，pid)，才有明顯并行優(yōu)化性能體現(xiàn)。具體實驗環(huán)境的系統(tǒng)互聯(lián)網(wǎng)絡(luò)對進(jìn)程間通信效率起著決定性作用。

3．3SW流程并行優(yōu)化性能測試

測試環(huán)境：Dell SC1435服務(wù)器

a)pc－cluster 2 nodes， 2 processors per node；

b)內(nèi)存：4 GB DDR2 667

測試數(shù)據(jù)：人類22號染色體第79號芯片Affy(2005) Ti－ling Array數(shù)據(jù)

測試結(jié)果（串行程序運行時間：7 084．79 s）如表1和圖8所示。

當(dāng)處理器個數(shù)NP＝1時，執(zhí)行并行程序所用的時間要比串行程序長。這主要是因為系統(tǒng)在處理MPI初始化以及處理器間數(shù)據(jù)收集發(fā)送過程中有一定的時間開銷。隨著處理器個數(shù)增多，分配給每個處理器處理的記錄條數(shù)減少，每個處理器的執(zhí)行時間也相應(yīng)減少，且都遠(yuǎn)大于系統(tǒng)通信時間，因此總的并行運行時間呈減少趨勢，加速比也相應(yīng)地不斷增大。當(dāng)處理器個數(shù)增加時，初始分發(fā)數(shù)據(jù)的讀取開銷增大，以及多處理器與主處理器通信成本增高，使得并行效率有所降低。當(dāng)處理器個數(shù)增加到一定數(shù)量時，則并行處理的優(yōu)勢將與各種代價相抵。但是依然可以依靠并行劃分?jǐn)?shù)據(jù)庫的方法獲得較高的加速比。Tiling Array原始數(shù)據(jù)庫越大，則加速比越明顯。若計算機(jī)單CPU內(nèi)存足夠大，可以考慮將提取后的數(shù)據(jù)文件在內(nèi)存中備份，這樣在信號識別階段可以直接到內(nèi)存中讀取，而不需要I/O操作，可以節(jié)省很多讀取開銷。

最后可進(jìn)行一個推算。Affymetrix（2005）Tiling Array79號芯片上有1 638 400個探針，單CPU運行信號識別流程需要將近2 h。人類22號染色體總共被覆蓋在五張芯片上，且它是人類第二短染色體，所以想要完成人類整個基因組Tiling Array的信號識別流程，至少需要花費數(shù)天時間，這與生物信息學(xué)高時效性的要求相悖。為了能夠提高信號識別流程的效率，借助于計算機(jī)強(qiáng)大的計算能力和高效并行算法，把幾天縮為幾小時，為后續(xù)的全基因組表達(dá)信息識別等工作贏得先機(jī)。

4結(jié)束語

當(dāng)前，基因芯片已經(jīng)是生物信息學(xué)領(lǐng)域非常具有生命力的研究方向，Tiling Array僅是其中的一個分支，它可以高效識別大規(guī)模表達(dá)區(qū)域。但是Tiling Array數(shù)據(jù)量龐大，需要高性能計算機(jī)的支持和優(yōu)秀并行算法的協(xié)助。本文介紹并比較了SW和SD兩種Tiling Array信號識別算法，選擇較優(yōu)的一個加以實現(xiàn)，并用數(shù)據(jù)并行分割的方法進(jìn)行并行優(yōu)化，對大數(shù)據(jù)量的信號識別有效縮減了運行時間，對其他信號識別算法也提供了有意義的方法支持，從一個嶄新的角度對基因芯片數(shù)據(jù)分析和處理進(jìn)行了啟發(fā)式的研究。

參考文獻(xiàn)：

［1］KAPRANOV P，CAWLEY S E，DRENKOW J，et al．Large－scale transcriptional activity in chromosome 21 and 22[J].Science，2002，296:916－919.

[2]SCHADT E E， EDWARDS S W，GUHA－THAKURTA D， et al. A comprehensive transcript index of the human genome generated using microarray and computational approaches[J].Genome Biol，2004，5(10):73.

[3]KAMPA D， CHENG J，KAPRANOV P， et al． Novel RNAs identified from in－depth analysis of the transcriptome of human chromosomes 21 and 22[J].Genome Res，2004，14:331－342.

[4]HUBER W，TOEDLING J，STEINMETZ L M， et al.Transcription mapping with high－density oligonucleotide Tiling Arrays[J].Gene Expression，2006，22(16):1963－1970.

[5]EMANUELSSON O，NAGALAKSHMI U， ZHENG De－you， et al. Assessing the performance of different high－density Tiling Microarray strategies for mapping transcribed regions of the human genome[J].Genome Res，2007，17:886－897.

[6]BERTONE P， STOLC V， ROYCE T E， et al.Global identification of novel transcribed sequences in human using high－resolution genomic Tiling Arrays[J].Science，2004，306:2242－2246.

[7]ROYCE T E， ROZOWSKY J S， BERTONE P， et al. Issues in the analysis of oligonucleotide Tiling Microarrays for transcript mapping[J].Trends Genet，2005，21(8):466－475.

[8]CHENG J， KAPRANOV P， DRENKOW J， et al. Transcriptional maps of 10 human chromosomes at 5－nucleotide resolution[J].Science，2005，308:1149－1154.

[9]QUINN M J．Parallel programming in C with MPI and OpenMP[M]. Beijing:Tsinghua University Press，2005.

[10]陳國良.并行計算——結(jié)構(gòu)·算法·編程[M].修訂版．北京:高等教育出版社，2004.

注：本文中所涉及到的圖表、注解、公式等內(nèi)容請以PDF格式閱讀原文