蔡 葵，楊進(jìn)才

(1.武漢理工大學(xué)華夏學(xué)院，湖北武漢 430223;2.華中師范大學(xué)計(jì)算機(jī)科學(xué)系，湖北武漢 430079)

針對 DNA 序列拼接，REPS［1］方法通過度量定長子串來確定重復(fù)序列。文獻(xiàn)［2-3］以此為基礎(chǔ)，提出了各自的定長為k的重復(fù)序列識別屏蔽方法。文獻(xiàn)［4］進(jìn)一步提出了預(yù)歸并的思想，但仍然局限于定長子串識別方法，對偏長的DNA重復(fù)序列處理能力不強(qiáng)，識別不夠精確。

筆者在PreMerged方法的基礎(chǔ)上，引入HUANG［5］等提出的 Superword array 概念，繼續(xù)預(yù)歸并思想，提出了利用變長子串來識別且屏蔽重復(fù)序列的VPreM方法。模擬實(shí)驗(yàn)證明，較之Pre-Merged方法，VPreM方法識別重復(fù)序列的精確度更高，并且可以更大規(guī)模地縮減shotgun集合規(guī)模。

1 基于變長子串的預(yù)歸并算法

1.1 基本思想

將shotgun集合中所有fragments片段末尾加上字符#，并串連成一大字符串F，為F中每個(gè)字符按順序編號。由字符集{A，G，C，T}上w個(gè)字符組成的字符串，就是一個(gè)長為w的word，長為w的word一共有4w個(gè)。將這4w個(gè)word按字母順序編號，形成一個(gè)word字查找表。表1為w值為2的word字查找表。

表1 w為2的word字查找表

大字符串 F上以p位置起始，長為 w的word，其在查找表中的編號Nu，就是p位置的編碼，記為Code(p)=Nu。若word中含#，則編號統(tǒng)一設(shè)為Nu=-1。圖1為一個(gè)大字符串F編碼示意圖。

通過編碼可以找出任意兩個(gè)fragments之間，連續(xù)h個(gè)長為w的word組成的重復(fù)部分repeats。假設(shè)shotgun集合中fragments總數(shù)目為n，可以推導(dǎo)出h及w的取值滿足h×w≥log4n。確定h和w的值之后，就可以為大字符串F編碼。然后根據(jù)不同fragments之間編碼的連續(xù)相同性(變長子串)，來進(jìn)行重復(fù)序列識別及預(yù)歸并屏蔽工作。

圖1 大字符串F編碼示意圖

1.2 構(gòu)建h×w變長子串統(tǒng)計(jì)表S

h×w變長子串統(tǒng)計(jì)表 S={Pos，F(xiàn)ragNum，Loc，CodeArray}，其中Pos為字符在F中的位置;FragNum為該字符所屬fragment的編號;Loc為該字符在所屬fragment中的位置;CodeArray={Code(Pos)，Code(Pos+w)，…，Code［Pos+(h-1) ×w］}，存儲連續(xù) h個(gè)長為 w的 word編碼。按CodeArray值排序后，對應(yīng)圖1的h×w變長子串統(tǒng)計(jì)表S1如表2所示。

表2 h×w變長子串統(tǒng)計(jì)表S1

先在表S1中去除所有CodeArray值中含-1的記錄，再去除僅部分與CodeArray值相同的記錄，增加FR項(xiàng)表示相同CodeArray值個(gè)數(shù)。得到優(yōu)化后的h×w變長子串統(tǒng)計(jì)表S2，如表3所示。

表3 優(yōu)化后的h×w變長子串統(tǒng)計(jì)表S2

設(shè)shotgun集合是DNA目標(biāo)序列的C條克隆隨機(jī)打散而成，其閾值T=C?？烧J(rèn)為表S2中FR值不超過T對應(yīng)的子串不屬于某個(gè)repeats;反之則認(rèn)為屬于某個(gè)repeats。

1.3 變長子串的直接預(yù)歸并

從表S2中取出FR值沒有超過閾值T且CodeArray項(xiàng)值相同的記錄，根據(jù)記錄中的Frag-Num和Loc值，直接在 shotgun集合 FRAG［n］中預(yù)歸并這些fragments為一個(gè)新的更大的fragment。其算法如下:

1.4 預(yù)歸并及識別重復(fù)序列

在表S2中去除所有FR值處于閾值T之內(nèi)的記錄，則剩余記錄均為與repeats相關(guān)的變長子串信息。根據(jù)每條記錄中的FragNum和Loc值，可以分別在shotgun集合中向左、向右預(yù)歸并對應(yīng)fragments，并借此識別出包含同一h×w變長子串的repeats的邊界。其算法如下:

i=1;

while(表S2未結(jié)束){

取表S2的第i條記錄，記為Si;

取Si的FR值，記為Time;

//以下為往左分析repeats信息

取第i至第i+Time-1條記錄中Loc值最大的記錄，記為St;

將St對應(yīng)的fragment編號記入NL數(shù)組;

取St的Pos值，記為Pt;

For(表S2的第i+1條至第i+Time-1條記錄){

取當(dāng)前記錄Sc的Pos值，記為Pc;

以 Code［Pc］與 Code［Pt］為起點(diǎn)，同步往左逐個(gè)比較Code值;

遇到第一個(gè)非負(fù)且不同的Code值，將Sc對應(yīng)的fragment編號記入NL數(shù)組;

截取Si與 Sc對應(yīng) fragment從起點(diǎn)到不同Code值位置片段，放入RepL集中;

截取Sc對應(yīng)fragment從不同Code值位置開始到Pc+h×w-1位置片段，記為REL;

}

//以下為往右分析repeats信息

取第i至第i+Time-1條記錄中對應(yīng)fragment長度減Loc值最大的記錄，記為St;

將St對應(yīng)的fragment編號記入NR數(shù)組;

取St的Pos值，記為Pt;

For(表S2的第i+1條至第i+Time-1條記錄){

取當(dāng)前記錄Sc的Pos值，記為Pc;

以 Code［Pc+h×w］與 Code［Pt+h×w］為起點(diǎn)，同步往右逐個(gè)比較Code值;

遇到第一個(gè)非負(fù)且不同的Code值，將Sc對應(yīng)的fragment編號記入NR數(shù)組;

截取Si與Sc對應(yīng)fragment從不同Code值位置到終點(diǎn)片段，放入RepR集中;

截取Sc對應(yīng)fragment從Pc+h×w開始到不同Code值位置片段，記為RER;

}

識別出第i至第i+Time-1條記錄中的Repeats=REL+RER;

順序提取 NL中的 fragment編號 u，更新FRAG［u］=RepL中對應(yīng)的左邊界;

順序提取 NR中的 fragment編號 v，更新FRAG［v］=RepR中對應(yīng)的右邊界;

i=i+Time;

}

1.5 方法的計(jì)算復(fù)雜性分析

h×w變長子串統(tǒng)計(jì)表S是通過掃描大字符串F得到的，其復(fù)雜性為Ο(l·n·h·4w)，其中l(wèi)為片段平均長度，n為shotgun集合中片段總數(shù)目，h為連續(xù)word的個(gè)數(shù)，w為字word的長度。原始表S經(jīng)優(yōu)化后，其中記錄數(shù)降為［T·4(h·w)］(T為閾值)，隨后的一系列預(yù)歸并及shotgun片段更新的復(fù)雜性，均在Ο(l·n·T)之內(nèi)，因此綜合得到的預(yù)歸并復(fù)雜性為Ο(l·n·T2·4(h·w))。

將其與文獻(xiàn)［2-3］方法的復(fù)雜性Ο(l·n·4k)相比較，由于h×w的值可視為與k相當(dāng)，因此筆者提出的方法復(fù)雜性比文獻(xiàn)［2-3］中的方法多了T2項(xiàng)，略大。再與文獻(xiàn)［4］方法的復(fù)雜性Ο(l·n·42k)比較，由于T2比4k小，因此筆者提出的方法復(fù)雜性比文獻(xiàn)［4］方法稍小。實(shí)質(zhì)上，l、T、h和 w均是有定值的，以上方法的復(fù)雜性最終都可以轉(zhuǎn)化為Ο(n)。筆者提出的方法是通過變長子串來識別repeats，識別精確性有所提高，對fragments的預(yù)歸并能力較之定長子串方法也有所增強(qiáng)。

2 計(jì)算機(jī)模擬分析

筆者提出的算法命名為VPreM，使用VC語言編程，DNA目標(biāo)序列和shotgun集合均用SQL Server 2003存儲，所有模擬實(shí)現(xiàn)均在Intel Core2 E4300(1.8 GHz)，1 G內(nèi)存(667 MHz)的PC機(jī)上完成。選取與文獻(xiàn)［4］相同的5條DNA序列來做模擬分析。根據(jù)克隆條數(shù)C和分解片段數(shù)目n，確定各個(gè)關(guān)鍵值:T=C;k=［log4n］+2;w=3;h=［k/w］+1。

用VPreM、PreMerged以及當(dāng)前流行的RECON［6］方法分別處理以上5個(gè) shotgun集合，最終識別結(jié)果如表4所示。3種方法的各自運(yùn)行時(shí)間如圖2所示。

表4 VPreM、PreMerged和RECON的識別結(jié)果

從表4可以看出，筆者提出的VPreM方法，在識別重復(fù)段的能力上明顯強(qiáng)于RECON，比Pre-Merged略強(qiáng)。由圖2可知VPreM方法的計(jì)算速度稍稍優(yōu)于PreMerged，更優(yōu)于RECON。

將以上5個(gè)shotgun集合分別用VPreM、Pre-Merged、RECON 3種方法預(yù)處理后，再交給Phrap拼接，運(yùn)行所耗費(fèi)時(shí)間比較如圖3所示。

圖2 VPreM、PreMerged和RECON的運(yùn)行時(shí)間

圖3 分別用VPreM、PreMerged和RECON預(yù)處理后Phrap拼接時(shí)間比較

由圖3可以觀察到，經(jīng)過VPreM方法預(yù)歸并再用Phrap工具進(jìn)行拼接，較之采用PreMerged和RECON兩種方法預(yù)處理后的shotgun集合，具有更快的速度，這也進(jìn)一步驗(yàn)證了筆者所提出方法的正確性。

3 結(jié)論

VPreM方法充分利用h×w的靈活可變性，能精確識別出不同長度repeats的邊界，較之基于定長子串的識別方法，效率更高;與預(yù)歸并方法的結(jié)合，也使得識別過程中獲取的信息得到了更充分利用?？梢詰?yīng)用于文獻(xiàn)［7-10］等對應(yīng)的拼接算法的前期處理工作之中，對提高拼接效率，減少拼接錯(cuò)誤，有一定的積極意義。

［1］WANG J.REPS:a sequence assembler that masks exact repeats identified from the shotgun data［J］.Genome Research，2002(12):824-831.

［2］張博峰，王正華.DNA片段拼接中基于定長特征子串的重復(fù)序列信息屏蔽方法［J］.國防科技大學(xué)學(xué)報(bào)，2002，24(6):67-70.

［3］王磊，張祖平，陳建二.DNA片段拼接中重復(fù)序列算法研究［J］.計(jì)算機(jī)科學(xué)，2006，33(7):164-170.

［4］蔡葵，楊進(jìn)才.DNA片段拼接中的預(yù)歸并重復(fù)序列屏蔽方法［J］.計(jì)算機(jī)工程，2009，35(4):88-90.

［5］HUANG X Q，YANG S P.Application of a superword array in genome assembly［J］.Nucleic Acids Research，2006，34(1):201-205.

［6］BAO Z R.RECON Documentation［EB/OL］.［2011-07-08］.http://selab.janelia.org/recon.html.

［7］PHIL G.Phrap documentation［EB/OL］.［2011-07-08］.http://www.phrap.org/phredphrapconsed.html.

［8］GRANGER G S，OWEN W.TIGR assembler:a new tool for assembling large shotgun sequencing projects［J］.Genome Science ＆ Technology，1995，1(1):9-19.

［9］HUANG X Q，ANUP M.CAP3:A DNA sequence assembly program［J］.Genome Research，1999(9):868-877.

［10］PEVZNER P A，TANG H X，WATERMAN M S.An eulerian path approach to DNA fragment assembly［J］.Proceedings of National Academy of Sciences，2001，98(17):9487-9753.