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

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

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

中圖分類號：TP312文獻標志碼：A

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

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

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

1Tiling Array構建和多種信號識別方法

1．1Tiling Array實驗平臺的構建方法

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

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

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

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

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

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

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

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

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

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

2SW信號識別流程的串行實現

Tiling Array實驗產生大量的數據，對這些數據重新整合格式，提取有用的信息，并采用恰當的信號識別算法加以處理，直至最后得到探針表達與否的衡量指標，是Tiling Array實驗目的所在，也是信號識別的基本流程。

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

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

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

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

b）根據數據格式設計合理的數據庫表，并將坐標數據和亮度信息分別導入。

c）依據PM探針坐標值對坐標數據表和亮度信息表進行連接查詢。

d）依據MM探針坐標值對坐標數據表和亮度信息表進行連接查詢。

e）設計數據庫表，將d）e）的查詢結果導入數據庫中，得到關于PM和MM探針全部信息的兩張表。

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

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

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

h）把經過如上處理步驟得到的數據庫表導出為結果數據文件。

2．2應用SW算法處理提取出的數據

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

3SW信號識別流程的并行優化

3．1并行可行化分析

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

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

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

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

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

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

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

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

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

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

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

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

3．2并行優化

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

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

每個處理器都要完成的任務是：

a）多次的部分數據庫表和另一數據庫表的連接查詢。

b）調用SW算法對部分提取出來的數據進行處理，并把結果寫入輸出文件。

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

3．3SW流程并行優化性能測試

測試環境：Dell SC1435服務器

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

b)內存：4 GB DDR2 667

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

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

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

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

4結束語

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

參考文獻：

［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]陳國良.并行計算——結構·算法·編程[M].修訂版．北京:高等教育出版社，2004.

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