楊永潔,鐘 誠(chéng)

(廣西大學(xué) 計(jì)算機(jī)與電子信息學(xué)院 廣西高校并行分布式計(jì)算技術(shù)重點(diǎn)實(shí)驗(yàn)室,南寧 530004)

1 引 言

下一代測(cè)序技術(shù)的應(yīng)用[1]使得short reads(短序列)以巨大的數(shù)量產(chǎn)生.為了確定short reads來(lái)源于基因組的哪一部分,在進(jìn)行全基因組測(cè)序時(shí),它們必須被映射(比對(duì))到參考序列(reference sequences)中.固有的基因組變異性和測(cè)序錯(cuò)誤,尤其是短的連續(xù)的序列的插入和缺失(gaps,空位)使得short-read比對(duì)(alignment)性能下降.為了將海量的short reads映射到參考基因組,人們?cè)O(shè)計(jì)開(kāi)發(fā)了一批諸如Bowtie[2]、SOAP2[3]、REAL[4]、BWA[5]、SOAP3[6]、Bowtie2[7]等short-read比對(duì)算法.這些算法允許在short reads比對(duì)中存在誤配(mismatches);然而,允許在short reads比對(duì)過(guò)程中插入空位(gap)時(shí),其性能較差,有的甚至完全不允許在比對(duì)過(guò)程中插入gap.但序列中存在gap的位置可能蘊(yùn)含重要的生物學(xué)信息,精確識(shí)別gap可以獲得完整的基因組圖譜,有利于后續(xù)的基因組信息解讀[8].

空位(gap)序列是比對(duì)中連續(xù)插入或刪除字符(indels)的最大序列.許多生物自然過(guò)程導(dǎo)致DNA序列中產(chǎn)生gaps:通過(guò)單一突變事件可以復(fù)制和插入長(zhǎng)DNA片段,染色體間DNA的易位[9].

隨著查詢序列中g(shù)ap長(zhǎng)度的增加,gap的出現(xiàn)頻率大大減少.在包含大約25-300個(gè)堿基對(duì)(bp)的short reads中,出現(xiàn)大量gaps的概率很低.因此,運(yùn)用動(dòng)態(tài)規(guī)劃方法求解序列比對(duì)問(wèn)題時(shí),不能限制比對(duì)中插入的gap數(shù)量,會(huì)極大影響short-read比對(duì)映射置信度.

為解決查詢序列和參考序列均為short reads序列的比對(duì)問(wèn)題,Flouri和Iliopoulos等人基于動(dòng)態(tài)規(guī)劃方法,提出了允許序列中出現(xiàn)單個(gè)gap的半全局比對(duì)算法Gapmis[10],但該算法比對(duì)過(guò)程相當(dāng)耗時(shí).為此,Alachiotis和Berger等人提出了libgapmis算法[11],它實(shí)現(xiàn)GapMis算法的并行化,加速了short reads的比對(duì)過(guò)程.為了在比對(duì)過(guò)程中插入多個(gè)gaps,Flouri等人提出了允許k-誤配以及在比對(duì)中插入l個(gè)gaps的半全局比對(duì)算法GapsMis[12].相比于已有的同類short reads比對(duì)算法,GapsMis算法允許序列中存在誤配的同時(shí),也允許插入可變但有限數(shù)量的gaps.但是對(duì)于不同物種,其固有的gap發(fā)生頻率不同.對(duì)于不同長(zhǎng)度的short reads,其gap出現(xiàn)的概率也不盡相同.GapsMis算法不能根據(jù)不同物種以及不同的short reads長(zhǎng)度在比對(duì)過(guò)程中插入最佳gaps,其在比對(duì)精確度、所需存儲(chǔ)空間以及比對(duì)時(shí)間還有改進(jìn)空間.

對(duì)于查詢序列和參考序列均為short reads的序列比對(duì)問(wèn)題,本文從NCBI[注]http://www.ncbi.nlm.nih.gov/網(wǎng)站獲得免費(fèi)的真實(shí)物種的DNA序列,采取訓(xùn)練查詢序列樣本數(shù)據(jù)找到不同物種以及不同read長(zhǎng)度匹配的最優(yōu)插入空位(gap)數(shù)量的策略,基于動(dòng)態(tài)規(guī)劃方法對(duì)short reads進(jìn)行兩兩比對(duì),在比對(duì)過(guò)程中插入最優(yōu)數(shù)量的gaps,并采用向量存儲(chǔ)比對(duì)過(guò)程中產(chǎn)生的中間矩陣中對(duì)角線區(qū)域上的元素,以提升比對(duì)精確度又減少比對(duì)時(shí)間,且降低所需存儲(chǔ)空間容量.

2 改進(jìn)的short-read比對(duì)算法

2.1 允許多gap的short reads比對(duì)問(wèn)題

Σ是一個(gè)有限的非空字符集,Σ上的字符串是Σ的有限元素序列,也可能是空串(空序列).空串用ε表示.

如果存在兩個(gè)字符串u和v,使得y=uxv,則字符串x是字符串y的因子.設(shè)字符串x,y,u和v,使得y=uxv.如果u=ε,那么稱x是y的前綴.如果v=ε,那么稱x是y的后綴.

設(shè)序列a=a0a1a2…ap-2ap-1,b=b0b1b2…bp-2bp-1.非空無(wú)gap的序列是一對(duì)對(duì)準(zhǔn)的長(zhǎng)度為p的有限非空序列:

(a0,b0),(a1,b1),…,(ap-1,bp-1),ai,bi∈Σ,0≤i

允許k-誤配和l-gaps(空位)的short-read半全局比對(duì)問(wèn)題定義如下:給定長(zhǎng)度為m的字符串x,長(zhǎng)度為n的字符串y(m≤n)以及整數(shù)k(0≤k

設(shè)序列x=GCGACGTCCGAA,序列y=GCAAGTAGA,k=4,l=2,并設(shè)x′=x[1,…,11]=GCGACGTCCGA,圖1給出了x′和y的對(duì)準(zhǔn)示意圖.

圖1 序列x′和y的對(duì)準(zhǔn)Fig.1 Alignment for x′ and y

圖1中“-”代表插入1個(gè)空位,其對(duì)準(zhǔn)實(shí)例的El(x,y)滿足最小且El(x,y)=4≤k,b=2≤l.對(duì)準(zhǔn)的序列z=GCAA-GT-AGA,其中z0=GCAA,g1=“-”,z1=GT,g2=“-”,z2=AGA.

對(duì)于兩個(gè)字符串x和y,xi為字符串x的第i個(gè)字符,yi為字符串y中的第i個(gè)字符,允許的編輯操作包括:

插入:向y中位置i插入一個(gè)字符yi,yi≠ε

刪除:刪除y中位置i的一個(gè)字符yi,yi≠ε

替換:用x中位置i的字符xi替換y的字符yi,xi,yi≠ε

為了記錄比對(duì)過(guò)程中的編輯操作次數(shù)以及將gaps插入序列的具體位置.構(gòu)造矩陣GS[0..n,0..m],其中GS[i,j]=ES(x[1..i],y[1..j])表示允許插入最多s個(gè)gaps時(shí),將x[1,…,i]轉(zhuǎn)化為y[1,…,j]所需的最小編輯操作次數(shù),s=1,2,…,l.同時(shí)構(gòu)造矩陣HS[0..n,0..m]記錄插入gaps的準(zhǔn)確位置.矩陣元素HS[i,j]計(jì)算如下[12]:

(1)

矩陣HS的計(jì)算表示插入gap的方向.將gap插入到y(tǒng)因子時(shí),HS[i,j]的值為正整數(shù).將gap插入到x因子時(shí),HS[i,j]的值為負(fù)整數(shù),s=1,2,…,l.

下面通過(guò)實(shí)例說(shuō)明short-read比對(duì)過(guò)程中,迭代計(jì)算具有最多l(xiāng)個(gè)gaps時(shí),所產(chǎn)生的矩陣G1,G2,…,GS,矩陣H1,H2,…,HS的具體形式,s=1,2,…,l.設(shè)序列x=ACATCGACG,序列y=CATTCGACG,l=2,圖2給出了矩陣G1、H1實(shí)例,圖3給出了矩陣G2、H2實(shí)例.

圖2 矩陣G1和矩陣H1實(shí)例Fig.2 Matrix G1 and matrix H1

從圖2和圖3可以看出,l=2時(shí),最多需要迭代兩次,其中矩陣G2和矩陣H2的值的計(jì)算依賴于矩陣G1和矩陣H1中的值.圖3中矩陣H2的粗體部分記錄了插入gaps的準(zhǔn)確位置.

圖3 矩陣G2和矩陣H2實(shí)例Fig.3 Matrix G2 and matrix H2

2.2 改進(jìn)的short-read比對(duì)算法與分析

我們首先分析GapsMis算法[12],然后給出本文改進(jìn)的算法思想、形式化描述及復(fù)雜度分析.

2.2.1 GapsMis算法分析

GapsMis算法[12]迭代計(jì)算具有最多l(xiāng)個(gè)gaps的最小成本的對(duì)準(zhǔn),插入至多l(xiāng)個(gè)gaps來(lái)計(jì)算矩陣G1,G2,…,Gl.首先計(jì)算插入第1個(gè)gap產(chǎn)生的矩陣G1和H1.接著依次計(jì)算插入第s-1個(gè)gap產(chǎn)生的矩陣Gs-1和Hs-1(1

矩陣GS、HS的計(jì)算只依賴于矩陣GS-1和HS-1,而算法迭代計(jì)算的次數(shù)取決于插入的gap數(shù)量l的大小.當(dāng)插入第s個(gè)gaps時(shí),將位置(i,j)的查詢序列與參考序列進(jìn)行比對(duì),若插入第s個(gè)gap的對(duì)準(zhǔn)結(jié)果GS[i,j]和插入第s-1個(gè)gaps的對(duì)準(zhǔn)結(jié)果GS-1[i,j]相同,則HS-1[i,j]=HS-1[i,j].設(shè)u為計(jì)算GS-1[0,j]到GS-1[i,j]的最小代價(jià)、v為GS-1[i,0]到GS-1[i,j]的最小代價(jià)、w為將GS-1[i-1,j-1]擴(kuò)展到GS-1[i,j]的最小代價(jià).HS[i,j]為{u,v,w,Gs-1[i,j]}中的最小值:

(2)

為了盡可能取得半全局最優(yōu)對(duì)準(zhǔn),GapsMis算法執(zhí)行允許k-誤配和l-gaps(空位)的成對(duì)reads的半全局比對(duì)的時(shí)間復(fù)雜度與產(chǎn)生多少個(gè)中間矩陣有關(guān),即與插入的gap數(shù)量l有關(guān).若比對(duì)過(guò)程中,固定允許插入的gap數(shù)量l,則只能得到局部最優(yōu),降低了比對(duì)精確度.但如果增加gap數(shù)量l進(jìn)行比對(duì),則又會(huì)大幅增加比對(duì)時(shí)間.

2.2.2 改進(jìn)算法的設(shè)計(jì)思想

對(duì)于不同的物種,其基因存在不同的gap發(fā)生頻率.對(duì)于不同長(zhǎng)度的short reads,gap的出現(xiàn)概率也不同.允許插入的gap數(shù)量l與算法迭代計(jì)算的次數(shù)正相關(guān).本文通過(guò)訓(xùn)練查詢序列樣本找到不同物種、不同長(zhǎng)度的short reads插入的最佳gaps數(shù)量lopt.通過(guò)利用最佳gaps數(shù)量lopt,使得比對(duì)算法在保證比對(duì)精確度的同時(shí),也減少了算法迭代次數(shù)、減少在比對(duì)過(guò)程中產(chǎn)生的中間矩陣,進(jìn)而減少中間矩陣計(jì)算量,從而使得算法在較少的時(shí)間內(nèi),獲得半全局比對(duì)最優(yōu)解.

同時(shí),為了減少比對(duì)算法所需的存儲(chǔ)空間,我們分析GapsMis算法的比對(duì)過(guò)程發(fā)現(xiàn),對(duì)于給定最大編輯操作距離k,只有矩陣G1,G2,…,GS中距離主對(duì)角線不超過(guò)k的那些路徑才能提供有效比對(duì),也就是說(shuō)存儲(chǔ)單元GS[(m-k),…,m,j]為一個(gè)有效解區(qū)域.例如,設(shè)x=ACATCGACG,y=CATTCGACG,k=2,l=2,圖4給出了矩陣G2和矩陣H2中寬度為“2k+1”的有效解區(qū)域.

給定允許誤配的數(shù)量k,GS矩陣以及HS矩陣元素的有效值為中心對(duì)角線2k+1的部分.由于矩陣GS、HS的計(jì)算只依賴于矩陣GS-1和矩陣HS-1,所以本文算法只需要分別使用(2k+1)個(gè)長(zhǎng)度為m的向量存儲(chǔ)中間矩陣GS-1、HS-1、GS-1、HS中的(2k+1)條對(duì)角線區(qū)域的矩陣元素?cái)?shù)據(jù),從而減少了算法所需的存儲(chǔ)空間容量.

2.2.3 改進(jìn)算法的描述

圖4 矩陣G2和矩陣H2Fig.4 Matrix G2 and matrix H2

在描述本文改進(jìn)的算法之前,首先介紹涉及到的概念.對(duì)于每一個(gè)可能的核苷酸匹配或不匹配,使用評(píng)分矩陣(scoring matrix)NUC4.4[注]ftp://ftp.ncbi.nih.gov/blast/matrices/NUC.4.4為每個(gè)可能的核苷酸分配分?jǐn)?shù),以及對(duì)插入的gap進(jìn)行仿射間隙罰分,最終比對(duì)總分最大的是最優(yōu)比對(duì).設(shè)gappen為每個(gè)比對(duì)的gap罰分,當(dāng)gap序列長(zhǎng)度r>0時(shí),gappen計(jì)算如下:

gappen=gop+(r-1)×gep

(3)

其中g(shù)op為gap opening penalty(間隙基本懲分),gep為gap extending penalty(間隙擴(kuò)展懲分).最后調(diào)用函數(shù)GapPos[12]查找插入gaps的位置來(lái)確定該對(duì)準(zhǔn)結(jié)果.GapPos函數(shù)以矩陣HS作為輸入,輸出最終的gaps插入數(shù)量β(β≤s)以及以下3個(gè)數(shù)組:

數(shù)組gap_pos的大小為s,gap_pos[i]表示在第i個(gè)位置插入gap或0,0≤i≤s-1.

數(shù)組gap_len的大小為s.gap_len[i]表示第i個(gè)位置插入gap時(shí)的gap的長(zhǎng)度,0≤i≤s-1.

數(shù)組where的大小為s.如果第i個(gè)gap插入到y(tǒng)中,則where[i]=1;如果第i個(gè)gap插入到x中,則where[i]=2,0≤i≤s-1.

當(dāng)按照給定的誤配以及gap發(fā)生頻率構(gòu)造參考序列片段(目標(biāo)序列)時(shí),已知所有誤配和插入的gap的準(zhǔn)確位置以及目標(biāo)序列的長(zhǎng)度.定義插入gap的數(shù)量小于或等于目標(biāo)序列中插入的gaps數(shù)量的對(duì)準(zhǔn)是有效對(duì)準(zhǔn)(Valid).總長(zhǎng)度小于或等于目標(biāo)序列的長(zhǎng)度,并且誤配的數(shù)量小于或等于目標(biāo)序列誤配的數(shù)量,則稱那些有效的對(duì)準(zhǔn)是正確對(duì)準(zhǔn)(Correct).即對(duì)準(zhǔn)序列為有效對(duì)準(zhǔn)、正確對(duì)準(zhǔn)時(shí),滿足如下條件[12]:

Valid:l≤L

(4)

Correct:l≤L,z≤Z,Rlength≤Tlength

(5)

其中l(wèi)是對(duì)準(zhǔn)序列中插入的gaps數(shù)量,L是目標(biāo)序列的長(zhǎng)度,z是對(duì)準(zhǔn)序列中的誤配數(shù)量,Z是目標(biāo)序列中的誤配數(shù)量,Rlength是對(duì)準(zhǔn)序列長(zhǎng)度,Tlength是目標(biāo)序列長(zhǎng)度.

按是否滿足式(4)和式(5)的條件,可以判斷對(duì)準(zhǔn)序列是有效對(duì)準(zhǔn)還是無(wú)效對(duì)準(zhǔn)、是正確對(duì)準(zhǔn)還是不正確對(duì)準(zhǔn).設(shè)變量V為有效對(duì)準(zhǔn)的序列數(shù)量,C為正確對(duì)準(zhǔn)的序列數(shù)量,G為參與比對(duì)的序列數(shù)目,則Short-read比對(duì)精確度acc計(jì)算如下[12]:

(6)

算法1形式描述了本文給出改進(jìn)的允許k-誤配和l-gaps(空位)的short-read的半全局比對(duì)算法(簡(jiǎn)記為GapsMis-adapt算法).

算法1.GapsMis-adapt

Input:物種種類sp,查詢序列集x0[0,…,m-1],x1[0,…,m-1],…,xT-1[0,…,m-1],參考序列集:y0[0,…,m-1],y1[0,…,m-1],…,yT-1[0,…,m-1],gop,gep,允許誤配數(shù)k

Output:比對(duì)精確度

Begin

1:訓(xùn)練查詢序列樣本數(shù)據(jù),獲取物種sp的查詢序列長(zhǎng)度m的最佳的插入gap的數(shù)目lopt;

2:Forp=0 tot-1 do

3: Fors=0 tolopt-1 do

4: 初始化矩陣GS和矩陣HS;

5: Forw=0 tos-2 do

6: Forj=0 tom-1 do

7: 計(jì)算矩陣元素GW[(m-k)…k,j]和HW[(m-k)…k,j],并分別用向量G0[0,…,m-1],G1[0,…,m-1],…,G2k[0,…,m-1]和H0[0,…,m-1],H1[0,…,m-1],…,H2k[0,…,m-1]存儲(chǔ)中間矩陣GW+1,GW+2,HW+1,HW+2中的(2k+1)條對(duì)角線區(qū)域的矩陣元素?cái)?shù)據(jù);

End

End

8: 以矩陣HS作為輸入,調(diào)用函數(shù)GapPos輸出gaps插入數(shù)量β(β≤s);

9: 按式(3)對(duì)插入的gap進(jìn)行仿射間隙罰分,最終總分最大的是最優(yōu)化比對(duì)為該查詢序列xp[0,…,m-1]的最優(yōu)對(duì)準(zhǔn)結(jié)果;

End

10: 按式(4)選出所有對(duì)準(zhǔn)序列中的有效對(duì)準(zhǔn)序列;

11: 在有效對(duì)準(zhǔn)序列的基礎(chǔ)上,按式(5)選出正確對(duì)準(zhǔn)序列;

12: 按式(6)計(jì)算比對(duì)精確度,并用數(shù)組A[0,…,lopt-1]保存比對(duì)精確度的值;

End

2.2.4 改進(jìn)算法的精確度、時(shí)間復(fù)雜度和空間復(fù)雜度分析

(7)

給定允許誤配數(shù)k,有效的對(duì)準(zhǔn)區(qū)域存在于GS矩陣、HS矩陣中心2k+1條對(duì)角線的區(qū)域,使用(2k+1)個(gè)向量保存該區(qū)域所占用的內(nèi)存空間為(2k+1)m.本文算法中矩陣GS、矩陣HS的計(jì)算只依賴于矩陣GS-1、矩陣HS-1,且僅需用(2k+1)個(gè)向量存儲(chǔ)矩陣中(2k+1)條對(duì)角線上的矩陣元素的值.因此,對(duì)于g條不同物種、不同長(zhǎng)度的short reads序列的比對(duì),改進(jìn)的GapsMis-adapt算法空間復(fù)雜度為O((2k+1)m)=O(km).這遠(yuǎn)遠(yuǎn)小于GapsMis算法的空間復(fù)雜度O(mn),m為查詢序列的長(zhǎng)度,n為參考序列的長(zhǎng)度.

3 實(shí) 驗(yàn)

實(shí)驗(yàn)使用的計(jì)算機(jī)為4核AMD Athlon(tm)ⅡX4 645 1900MHz處理器、內(nèi)存容量8GB.運(yùn)行操作系統(tǒng)Linux.采用C以及Python編程實(shí)現(xiàn)算法.為了模擬真實(shí)比對(duì)環(huán)境,通過(guò)網(wǎng)站NCBI下載了3種公開(kāi)的真實(shí)物種的DNA序列,分別為擬南芥(Arabidopsis thaliana,AT),甜菜(Beta vulgaris,BV)[13],智人(Homo sapiens,HS)[14].在實(shí)驗(yàn)將參考序列y與查詢序列x以FASTA格式作為輸入,然后運(yùn)行算法進(jìn)行比對(duì),最后生成一個(gè)對(duì)準(zhǔn)結(jié)果的文本文件.使用修飾符-l傳入允許插入gaps的最大數(shù)量l.

我們將本文算法GapsMis-adapt與代表性的同類算法GapsMis、needle[15]進(jìn)行實(shí)驗(yàn)測(cè)試,并比較3種算法的有效對(duì)準(zhǔn)數(shù)量、精確度、運(yùn)行時(shí)間.

3.1 構(gòu)建參考序列和查詢序列

為了通過(guò)訓(xùn)練查詢序列樣本數(shù)據(jù),確定不同物種以及不同read長(zhǎng)度比對(duì)的最優(yōu)插入空位(gap)數(shù)量,并對(duì)本文算法GapsMis-adapt和算法GapsMis、needle的性能進(jìn)行評(píng)估.我們下載了5條AT的基因序列(AT1～AT5),18條BV基因序列(BV1～BV9,BV1_2320～BV9_2320),22條HS基因序列(HS1～HS22),共3.6GB.對(duì)于每一條基因序列構(gòu)造20000對(duì)read長(zhǎng)度分別為50bp、100bp、150bp、200bp、250bp、300bp的參考序列,并按照給定的誤配發(fā)生頻率和gap發(fā)生頻率構(gòu)造對(duì)應(yīng)的查詢序列.

表1給出了插入不同物種的誤配發(fā)生頻率以及gap發(fā)生頻率的信息.

表1 誤配及插入gap的頻率
Table 1 Frequency of inserted mismatches and gaps

物種名稱誤配發(fā)生頻率gap發(fā)生頻率AT1.6×10-32.4×10-5BV1.6×10-31.7×10-5HS1.6×10-35.7×10-6

3.2 算法比對(duì)精確度

許多自然過(guò)程導(dǎo)致DNA序列中產(chǎn)生gaps,精確識(shí)別gap是疾病研究的基礎(chǔ).我們通過(guò)模擬長(zhǎng)度分別為50bp、100bp、150bp、200bp、250bp、300bp的查詢序列,使用誤配發(fā)生頻率和3種不同的gap出現(xiàn)頻率構(gòu)造參考序列,并進(jìn)行實(shí)驗(yàn).算法needle(gop,gep)帶有參數(shù)gop和gep.算法GapsMis-adapt以及算法GapsMis和needle的比對(duì)精確度如表2所示.

表2 GapsMis、needle、GapsMis-adapt算法的比對(duì)精確度
Table 2 Accuracy of alignments for GapsMis,needle and GapsMis-adapt

speciesLengthGapsMis(l=1)GapsMis(l=2)GapsMis(l=3)GapsMis(l=4)GapsMis(l=6)Needle(10.0,0.5)GapsMis-adaptAT5099.900%99.900%99.900%99.900%99.900%99.895%99.900%AT10099.905%99.925%99.925%99.925%99.925%99.920%99.925%AT15099.845%99.905%99.905%99.905%99.905%99.895%99.905%AT20099.825%99.930%99.930%99.930%99.930%99.920%99.930%AT25099.730%99.885%99.885%99.885%99.885%99.865%99.885%AT30099.660%99.895%99.895%99.895%99.895%99.870%99.895%BV5099.950%99.955%99.955%99.955%99.955%99.950%99.955%BV10099.925%99.950%99.950%99.950%99.950%99.945%99.950%BV15099.905%99.940%99.940%99.940%99.940%99.935%99.940%BV20099.880%99.955%99.955%99.955%99.955%99.955%99.955%BV25099.825%99.935%99.935%99.935%99.935%99.930%99.935%BV30099.840%99.950%99.950%99.950%99.950%99.950%99.950%HS5099.795%99.825%99.825%99.825%99.825%99.810%99.825%HS10099.670%99.815%99.815%99.815%99.815%99.800%99.815%HS15099.505%99.805%99.820%99.820%99.820%99.805%99.820%HS20099.235%99.790%99.815%99.815%99.815%99.805%99.815%HS25098.920%99.775%99.820%99.820%99.820%99.805%99.820%HS30098.570%99.745%99.810%99.815%99.815%99.805%99.815%

對(duì)于不同物種的基因,不同長(zhǎng)度的read其生物自然過(guò)程導(dǎo)致DNA序列中產(chǎn)生gaps的概率不同.故設(shè)定每組實(shí)驗(yàn)在序列比對(duì)過(guò)程中允許插入1個(gè)gap、2個(gè)gaps、3個(gè) gaps、4個(gè)gaps以及6個(gè)gaps.允許插入6個(gè)gaps進(jìn)行實(shí)驗(yàn),是為確保算法GapsMis與GapsMis-adapt進(jìn)行比對(duì)時(shí),GapsMis算法能盡可能選取半全局比對(duì)最優(yōu)解.

表2的結(jié)果表明算法GapsMis隨著short reads中允許插入gap數(shù)量l的增加,比對(duì)精確度逐步提高,但是當(dāng)l≥4時(shí),比對(duì)精確度不再提高.此外,當(dāng)read長(zhǎng)度大于50bp時(shí),GapsMis(l=1)在任何情況下都具有最低的比對(duì)精確度.

另一方面,從表2的結(jié)果可知,本文算法GapsMis-adapt的比對(duì)精確度只有1個(gè)和needle算法的精確度相同,其他的全部高于needle算法;當(dāng)l>2時(shí),GapsMis算法的比對(duì)精確度只有2個(gè)和needle算法的精確度相同,其他的全部高于needle算法.當(dāng)l<4時(shí),算法GapsMis的比對(duì)精確度低于或等于本文算法GapsMis-adapt;當(dāng)l≥4時(shí),算法GapsMis的精確度與本文算法GapsMis-adapt的比對(duì)精確度相同.

正如文獻(xiàn)[13]指出,算法needle為了將short reads序列與參考序列正確對(duì)準(zhǔn),在比對(duì)過(guò)程中,會(huì)插入過(guò)多數(shù)量的gaps,導(dǎo)致無(wú)法得到滿足有效對(duì)準(zhǔn)(Valid)定義的結(jié)果.在默認(rèn)情況下,needle算法的gop為10.0,gep為0.5.按式(3),可以增加gop的值,限制needle算法在比對(duì)過(guò)程中插入gaps的數(shù)量.下面我們通過(guò)實(shí)驗(yàn)來(lái)測(cè)試3種不同gop值對(duì)needle算法性能的影響,同時(shí)與算法GapsMis、 GapsMis-adapt進(jìn)行有效對(duì)準(zhǔn)、比對(duì)精確度的比較.

表3 算法needle、GapsMis、GapsMis-adapt的有效對(duì)準(zhǔn)數(shù)和比對(duì)精確度
Table 3 Values ofVandaccof alignment for algorithms needle,GapsMis,and GapsMis-adapt

specieslengthneedle(10.0,0.5)(V/acc)needle(15.0,0.5)(V/ acc)needle(20.0,0.5)(V/ acc)GapsMis(l=4)(V/ acc)GapsMis-adapt(V/ acc)AT5019998/99.895%20000/99.890%20000/99.835%20000/99.900%20000/99.900%AT10019999/99.920%20000/99.915%20000/99.865%20000/99.925%20000/99.925%AT15019999/99.895%20000/99.885%20000/99.825%20000/99.905%20000/99.905%AT20019999/99.920%20000/99.910%20000/99.850%20000/99.930%20000/99.930%AT25019998/99.865%20000/99.865%20000/99.815%19999/99.885%19999/99.885%AT30019997/99.870%19999/99.880%20000/99.825%19999/99.895%19999/99.895%BV5020000/99.950%20000/99.940%20000/99.910%20000/99.955%20000/99.955%BV10019999/99.945%19999/99.935%20000/99.895%20000/99.950%20000/99.950%BV15020000/99.935%20000/99.930%20000/99.890%20000/99.940%20000/99.940%BV20020000/99.955%20000/99.950%20000/99.915%20000/99.955%20000/99.955%BV25020000/99.930%20000/99.915%20000/99.875%20000/99.935%20000/99.935%BV30020000/99.950%20000/99.940%20000/99.910%20000/99.950%20000/99.950%HS5019997/99.810%19998/99.780%19999/99.645%20000/99.825%20000/99.825%HS10019996/99.800%19997/99.775%19999/99.645%19999/99.815%19999/99.815%HS15019997/99.805%19998/99.780%20000/99.640%19999/99.820%19999/99.820%HS20019998/99.805%19998/99.775%20000/99.640%20000/99.815%20000/99.815%HS25019997/99.805%19998/99.775%20000/99.625%19999/99.820%19999/99.820%HS30019997/99.805%19998/99.775%20000/99.635%20000/99.815%20000/99.815%

表3給出了算法needle、GapsMis、GapsMis-adapt的有效對(duì)準(zhǔn)數(shù)V和比對(duì)精確度acc,其中needle(10.0,0.5)、 needle(15.0,0.5)、needle(20.0,0.5)對(duì)應(yīng)于gop分別為10.0、15.0、20.0和gep為0.5的情形,算法GapsMis、GapsMis-adapt中的gop為10.0、gep為0.5.

表3的結(jié)果表明,通過(guò)增加gop值,減少了needle算法在比對(duì)過(guò)程中插入gaps的數(shù)量,增加了有效(Valid)對(duì)準(zhǔn),但降低了算法needle的精確度(acc),這會(huì)導(dǎo)致大量的錯(cuò)配被低估,對(duì)映射質(zhì)量產(chǎn)生潛在的嚴(yán)重影響.

另一方面,表3的結(jié)果與needle算法相比,算法GapsMis-adapt和GapsMis都能得到更高精確度的 short-read對(duì)準(zhǔn).

為了評(píng)估本文算法GapsMis-adapt和算法needle、 GapsMis的運(yùn)行時(shí)間,我們采用AT、BV、HS物種的基因數(shù)據(jù)模擬生成5,400,000對(duì)長(zhǎng)度分別為50bp、100bp、 150bp、200bp、250bp、300bp的查詢序列.3個(gè)算法完成short-read比對(duì)所需運(yùn)行時(shí)間如圖5所示,圖5中A1表示needle(20.0,0.5),A2表示needle(15.0,0.5),A3表示needle(10.0,0.5),A4表示GapsMis(l=6),A5表示GapsMis(l=4),A6表示GapsMis(l=3),A7表示GapsMis(l=2),A8表示GapsMis(l=1),A9表示GapsMis-adapt.

圖5以及表3的結(jié)果表面,雖然算法GapsMis(l=1)的運(yùn)行時(shí)間最少,但其比對(duì)精確度也最低;而在確保較高比對(duì)精確度的前提下,算法GapsMis-adapt比對(duì)算法 needle、GapsMis更快地完成short-read比對(duì).這表明在確保short-read比對(duì)精確度的前提下,本文算法GapsMis-adapt效率更高.本文算法之所以能確保高的精確度、高的比對(duì)效率,主要是它根據(jù)不同的物種以及不同的reads長(zhǎng)度,在比對(duì)過(guò)程中采用最佳插入的gaps數(shù)量作為算法的迭代次數(shù),從而減少了算法所需的中間矩陣計(jì)算量.

圖5 GapsMis、needle、GapsMis-adapt算法的運(yùn)行時(shí)間Fig.5 Required time for algorithms GapsMis,needle and GapsMis-adapt

4 結(jié)束語(yǔ)

對(duì)于查詢序列和參考序列均為short reads的序列比對(duì)問(wèn)題,本文給出一種改進(jìn)的short-read 比對(duì)算法,它的特點(diǎn)是通過(guò)訓(xùn)練查詢序列樣本數(shù)據(jù)尋找不同物種以及不同read長(zhǎng)度匹配的最優(yōu)插入空位數(shù)量,以減少算法比對(duì)過(guò)程的迭代次數(shù),從而降低算法所需的中間矩陣計(jì)算量,在確保比對(duì)精確度的前提下,減少了比對(duì)計(jì)算時(shí)間;此外,采用向量存儲(chǔ)中間矩陣中對(duì)角線上的矩陣元素值,節(jié)省了所需存儲(chǔ)空間,使得算法空間復(fù)雜度下降到O(km),k