生物信息學(xué)序列拼接第1頁/共57頁序列拼接

序列拼接任務(wù)即將測序生成的reads短片段拼接起來，恢復(fù)出原始的序列。該問題是序列分析的最基本任務(wù)，是基因組研究成功與失敗的關(guān)鍵，拼接結(jié)果直接影響到序列標(biāo)注，基因預(yù)測、基因組比較等后續(xù)任務(wù)?；蚪M序列的拼接也是基因組研究必須解決的首要難題。其困難不僅來自它的海量數(shù)據(jù)(以人類基因組序列為例，從數(shù)量為10兆級的片斷恢復(fù)出長度為億級的原始序列)，而且源于它含有高度重復(fù)的序列。第2頁/共57頁第3頁/共57頁拼接問題的難點DNA測序數(shù)據(jù)有其固有的四個的特點，他們也正是解決實際的序列拼接問題的難點所在：

1．測序有誤差

2．不完全覆蓋性

3．序列所在鏈不確定

4．重復(fù)序列的干擾第4頁/共57頁1．測序有誤差

由于測序技術(shù)的局限，難免會出現(xiàn)測序錯誤，尤其是在序列的末端，一般錯誤率可控制在1％以下。所以對每個堿基一般有一個正確概率，以質(zhì)量打分的形式給出。因此每個ri都有個可信度。而read與read之間有不同程度的重疊，由此導(dǎo)致有的重疊可信度高，有的重疊可信度低。第5頁/共57頁2．不完全覆蓋性不是所有的堿基被測序的次數(shù)都等于平均測序覆蓋度。極端的情況，可能會出現(xiàn)源基因組序列上部分區(qū)域未被測序的情況(這段區(qū)域稱為gap)。即，測序的reads集合不是原始基因組序列一個完整覆蓋。此時需要借助于各種圖譜如：基因組指紋圖譜(genomefingerprintmap)，基因組級物理圖譜(genome-widephysicalmap)，細(xì)胞發(fā)生圖譜(cytogeneticmaps)等協(xié)助對reads進行定位．第6頁/共57頁3．序列所在鏈不確定由于測序過程中無法確定特定片斷屬于DNA雙鏈中的哪一條鏈上，所以我們在拼接過程中并不清楚使用的是read的正義鏈，還是其互補鏈。4．重復(fù)序列的干擾

DNA序列自身含有高度重復(fù)的子序列，它們一種表現(xiàn)為短序列的串級重復(fù)，比如：(GGAA)n?；駻mTn等。另一種表現(xiàn)為大量相似序列(其拷貝數(shù)可達(dá)幾十萬)散布在基因組的各個地方。Repeat的存在，將導(dǎo)致fragments間overlap的不真實性，進而產(chǎn)生錯拼的結(jié)果。因此在拼接過程中耍確定這些序列的形式及大小，才能保證以高概率恢復(fù)出其在原始真實序列中的位置．第7頁/共57頁拼接算法評價

以上拼接問題的四個難點不僅極大的增加了解決實際拼接問題的難度，而且從某種程度上說無法完整地恢復(fù)出原始DNA序列來。即實際上僅能構(gòu)建出若干個contig(重建的fragments的一種排列形式，它覆蓋基因組上一段連續(xù)區(qū)域)這些contig將指導(dǎo)測序項目finishing階段的實驗方法最終構(gòu)建DNA完整序列。第8頁/共57頁

目前，國際上對拼接軟件的公認(rèn)評價標(biāo)準(zhǔn)包括兩方面，即重建出的contig的數(shù)目和準(zhǔn)確度。我們發(fā)展的基因組序列拼接新算法的目標(biāo)是在確保準(zhǔn)確性的前提下，構(gòu)建盡量少的contig，以減少測序后期大量的人力和財力的投入。第9頁/共57頁基因組序列拼接算法研究現(xiàn)狀

現(xiàn)在最常用的拼接程序使用的拼接算法可分成兩類，一類是將拼接問題轉(zhuǎn)化為在圖中尋找的Hamilton路徑的問題；另一類是將拼接問題在某種特殊情況下轉(zhuǎn)化成尋求圖中的Euler路徑的問題。他們均有其成功的典型算法。第10頁/共57頁1.轉(zhuǎn)化為HamiltonPath問題

每個DNA片段(read)相當(dāng)于圖中一個結(jié)點，如果兩個片段之間存在著重疊(overlap)關(guān)系，則在兩個結(jié)點之間定義一條邊，而沿著DNA原始序列從頭到尾，則必然經(jīng)過每個結(jié)點一次且僅一次，即是一條Hamilton路徑。一條contig表示圖中一條簡單路，此類算法以Phrap，TIGRAssembler，CAP3，GigAssemble等為代表。第11頁/共57頁

他們都是遵循“overlap-layout-consensus”的框架。首先，為了構(gòu)建圖。計算任意兩個read間可能的比對情況。其次，通過去除歧義的或者不確信的邊得到較為準(zhǔn)確的圖，并在其上尋找非交叉的簡單路的集合，該集合對應(yīng)于contig的集合。最終，通過對包含在一個簡單路上的所有read進行多序列比對，為每一個contig構(gòu)建一個一致性序列(consensussequence)。第12頁/共57頁2.轉(zhuǎn)化為EulerPath問題EULER是這類算法的代表。與傳統(tǒng)方法沿著“Overlap—Layout—Consensus”路線不同，它不計算各個read之間的Overlap，即沒有Overlap步驟。它的大致想法如下：為了排除read中的錯誤，獲得Error-Free的read，將所有的read切割成小片n-mers。第13頁/共57頁

將每個read和Gk的近似進行比對，尋求read的最小改變能夠使得read的所有n-mers包含在Gk的近似集合中。從而構(gòu)建了高質(zhì)量序列，而對于Poorread，直接拋棄，對Chimericread(兩端在n-mers中但整體不在的reads)進行特殊處理。第14頁/共57頁

初始的想法是要實現(xiàn)去除reads中的測序錯誤的目的，如果知道原始序列G，那么直接使用測序獲得的read和G進行比較即可。但是實際上G并不可知，那么退而求其次，G的序列片斷Gk亦可，事實上Gk亦不可知。所以將所有的read切割成小片n-mers，所有Solid的n-mers形成的集合稱為Gk的近似。最后，構(gòu)造DeBruijn圖。第15頁/共57頁現(xiàn)有算法的主要問題

雖然已經(jīng)開發(fā)了以上的算法，基因組序列拼接問題尚未徹底解決，以上兩類算法都存在著各自的缺陷。第16頁/共57頁

對于第一類算法來說，實際上是在圖中尋找一條使得評價函數(shù)值最優(yōu)的Hamilton路徑，這是一個NP完全問題。一般都采用greedy-merging的算法近似求解。由于這種step-by-step的局部貪心算法，其明顯的局部特性忽略了reads間“長距離”或者整體性的聯(lián)系，從而導(dǎo)致了拼接錯誤，即拼接結(jié)果和真實的DNA原始序列不同。最近研究指出，在對已知序列的流行性感冒嗜血桿菌基因組的拼接過程中，無論是Phrap，TIGRAssembler，還是CAP3，都發(fā)生了拼接錯誤的現(xiàn)象。第17頁/共57頁

對于第二類算法來說，它只能在特殊的情況下，才能將問題簡化成尋求一條Euler路徑，最終的結(jié)果是從多條候選的Euler超路中選擇出來的。EULER算法依然存在拼接錯誤，且結(jié)果選擇的過程沒有理論依據(jù)。EULER軟件在實際數(shù)據(jù)集的運行速度上和第一類算法相當(dāng)。更重要的是，EULER采用的算法過于獨立，很難利用其他輔助生物信息，導(dǎo)致其實用性和流行性大打折扣。第18頁/共57頁局部搜索(LocalSearch)方法胡杰第19頁/共57頁

將局部搜索方法運用于一個具體的問題，需要對以下四項進行明確的定義。

1．將原同題表示成一個優(yōu)化問題，即定義一個可行域以及在該可行域上的一個目標(biāo)函數(shù)。

2．定義可行域中的鄰域結(jié)構(gòu)，即說明滿足什么條件的兩個點相鄰。第20頁/共57頁3．確定在鄰域中的搜索方式。4．局部極值點的處理。如果當(dāng)前解點鄰域中的所有點的目標(biāo)函數(shù)值都比當(dāng)前點大，那么這點就稱為局部極值點。在一些問題中，局部極值點就是全局最優(yōu)點。而對另一些問題而言，局部極值點已經(jīng)滿足求解實際問題的需求。第21頁/共57頁基于局部搜索的序列拼接算法框架

主要目標(biāo)是在layout階段盡可能準(zhǔn)確的前提下，獲得更長的contig。具體的，使用局部搜索算法求得數(shù)據(jù)集上近似全局的最短超串；然后根據(jù)求得的最短公共超串對應(yīng)的fragment的排列關(guān)系為基礎(chǔ)獲得“consensussegment”第22頁/共57頁1.Overlap定義

如果一個串的前綴是另一個串的后綴則認(rèn)為這兩個串之間存在overlap，并根據(jù)over-lap構(gòu)建超串。對給定的串f和g，存在多個可能的overlap關(guān)系．比如說，若f=ACTGGGAGCAGC，

g=AGCAGCTTTTACT，那么他們之間至少存在兩個overlap形式。第23頁/共57頁第24頁/共57頁

在我們的算法中，僅考慮兩個串之間最大的overlap情況，并定義overlap(f，g)表示f和g之間存在的多個overlap關(guān)系中最長的一個overlap所包含的字符個數(shù)。在上面的例子中overlap(f，g)=6。如果f和g之間overlap區(qū)域長度小于M(M是一個足夠小的正整數(shù))，則overlap(f，g)=0。第25頁/共57頁2.優(yōu)化目標(biāo)定義

我們對reads集合S中的每個元素按照其左端在超串t中首次出現(xiàn)的位置進行排序，并沿超串t從左向右依次讀入，并將其對應(yīng)為序S={sl，s2…．，Sn)。我們用P(S)表示這個由字符串為元素構(gòu)成的序。在序P(S)中，對于每一個連續(xù)的字符串元素對si，si+1，都存在overlap(si，Si+1)。第26頁/共57頁

因此，字符串的一個排列等價于一個超串t及作用在其上的函數(shù)overlap，超串t的長度length(t)=│t│=∑s∈sIs|-overlap(si，Si+1),

為了求超串t具有最小的長度。因為在給定集合s中∑s∈sIs│是確定的值，我們可以進一步把問題轉(zhuǎn)化成尋找一個排列P使得最大化。第27頁/共57頁3.鄰域定義

在使用局部搜索方法前，我們首先定義排列的鄰域這一概念。我們稱一個排列為問題的一個解。設(shè)reads集合S是一個具有n個元素的字符串集合，P(S)是reads集合S所有解的集合。定義操作rshift(i，j)

，即對一個解P∈P(S)，將P第i位置的元素向右移至第j個位置(1≤i<j≤n)。第28頁/共57頁也就是說，如果我們應(yīng)用這一操作于給定的解P={sl，s2…．，sn)，則結(jié)果P’=rshift(i，j)(P)定義為；第29頁/共57頁

類似的，定義操作lshift(i,j)將P第i位置的元素向左移至第j個位置(1≤j<i≤n)，定義P’=lshift(i，j)(P)為；第30頁/共57頁

這樣，一個解P∈P(s)的鄰域可以定義為N(P)={rshift(i，j)(P)|1≤i<j≤n)｝∪{lshift(i，j)(P)|1≤j<i≤n｝．第31頁/共57頁4．局部搜索算法

通過計算從可行解P到其鄰居P’∈N(P)的過程中的?overlap，并采用“first-fit的搜索策略，我們在當(dāng)前解的鄰域范圍內(nèi)搜索第一個更好的解，并轉(zhuǎn)移到該解，因此可以通過迭代的方法改善解的質(zhì)量，并實現(xiàn)潛在的糾錯功能。第32頁/共57頁例如，在γshift(i,j)操作下將p移至p’，?

overlap=overlap(si-1,si+1)+overlap(sj,sj+1)+overlap(si,sj+1)-overlap(si-1,si)-overlap(si,si+1)-overlap(sj,sj+1)。當(dāng)?overlap>0時，p’比p更好，則由p轉(zhuǎn)移到p’第33頁/共57頁

此外，對給定的解P={s1’，s2’…．Sn’｝,內(nèi)部的元素均包含前趨和后繼。但是，頭元素僅有后繼而尾元素僅有前驅(qū)。在算法中，為了消除頭尾元素差異，我們將排列的頭尾元素連結(jié)起來并使用局部搜索方法尋找最優(yōu)的“Loop”超串。接著在overlap最小的地方切分該環(huán)狀超串，最終還原成一條線性超串。第34頁/共57頁5.局部極值點的處理

當(dāng)經(jīng)過以read為單元的搜索后，可以獲得一個當(dāng)前鄰域內(nèi)的局部最優(yōu)解{γ1…,γk1,γk1+1,…,γk2,γk2+1，…，γk3γkm｝。它對應(yīng)集合s上的一個superstring。該解對應(yīng)的reads的排列中，任意相鄰兩個reads間的overlop關(guān)系薄弱，即overlap(γi，γi+1)<M(其中，M是一個足夠小的正整數(shù)，1≤i≤n)。第35頁/共57頁例如：

如果overlap(rkl，rkI+1)<M，

overlap(rk2，rk2+1)<M,…，

overlap(rkm，rkm+1)<M，令fragment集合F=｛f1,f2,f3….,fm}，fi即為在序rki-1,…,rki下對應(yīng)的一個子超串。第36頁/共57頁

則在此薄弱處將superstring分割若干“sub-superstrings.它們構(gòu)成新的fragment集合F。此時，我們將F作為當(dāng)前最短超串求解的子串集合，可以獲得一個新的是最短公共超串(SCS)問題的實例。我們繼續(xù)應(yīng)用前述的Overlap定義和鄰域的定義。最終得到reads集合上的一個局部最優(yōu)解P’。第37頁/共57頁6．獲得“consensussegments”

在處理對“contiguoussegments”的識別過程中，有關(guān)鍵的兩個步驟：(1)在求得的超串中識別“weakjoins”．(2)在“weakjoins”位置，獲得contiguoussegments．第38頁/共57頁7.Repeat的識別策略

為了鑒別潛在的falsejoins，我們必須分析所有的包含不一致的前綴或者后綴的fragment。具體的fragmentγ2和γ3具有不一致的后綴，即至少存在兩個fragmentsr2和r3,且滿足overlap(r1，r2)>M1，overlap(r1，r3)>M1，但overlap(r2，r3)<M2和overlap(r3，r2)<M2。第39頁/共57頁

因為fragment是對基因組DNA序列的隨機采樣獲得的，包含顯著的大量fragment的區(qū)域很可能是repeat。因此，如果滿足以下兩個條件，我們就認(rèn)為該contiguoussegment是重復(fù)的：

1)其首部或者尾部分別是包含不一致的前綴和后綴的fragment；

2)其覆蓋度遠(yuǎn)遠(yuǎn)大于測序的平均覆蓋度。否則，我們就認(rèn)為其為unique區(qū)域．第40頁/共57頁

原型系統(tǒng)試驗

趙芳方第41頁/共57頁試驗環(huán)境

1)原始基因組數(shù)據(jù)是從GenBank中隨機抽取并下載的；2)read序列是原始基因組序列中的隨機子串，長度為500到700個堿基對；3)read序列隨機的分布于原始基因組數(shù)據(jù)互補的兩鏈中的任意一條．將BasicLSA用于不同的數(shù)據(jù)中，得到的數(shù)據(jù)的三個特點：第42頁/共57頁試驗結(jié)果

我們在沒有雙端測信息的約束下，從生成的CONTIG的數(shù)量和準(zhǔn)確性兩個角度比較了BasicLSA.PHRAR.CAP3的性能。第43頁/共57頁原序列名序列長度覆蓋度系統(tǒng)Contigs總長Contigs數(shù)目長contigs數(shù)目(>50k)Contigs平均長度錯誤contigs數(shù)目錯誤contigs總長(bp)AE0078725428697.5CAP35414652132876300PHRAP5405881932845200B-LSA5410821442864800AE00065715513357.5CAP31539089777199886144539PHRAP1528509717215003104164B-LSA15459535011309192101102INRUENSS18301388CAP317957526812264085227941PHRAP17986456314287484188705B-LSA181116944124116215219AE00078221784007.5CAP3216969310311215096279380PHRAP21506509411223814212060B-LSA2171203751328949339270AE00786928415817.5CAP328269911261422436136274PHRAP28155571161524272229589B-LSA283368499203148500AL00912642148147.5CAP3416548417110243586138549PHRAP413779915519256955168824B-LSA419066712321314849122422CAP3、PHRAP、B-LAS的性能比較第44頁/共57頁結(jié)果一contigs的個數(shù)contig平均長度B-LSA40833.661PHRAP53325.483CAP356623.035----結(jié)果：產(chǎn)生更少的contig，而且contig更長有效的減少了測序項目后續(xù)工作的工作量第45頁/共57頁結(jié)果二總contigs數(shù)錯誤數(shù)目錯拼率B-LSA408153.67CAP3745526.97PHRAP566244.24----結(jié)果：包含更少的錯誤拼寫現(xiàn)象第46頁/共57頁

基于原型算法框架實現(xiàn)的BasicLSA系統(tǒng)在針對實際序列模擬的shotgunreads數(shù)據(jù)集上的拼接結(jié)果表明：在不引入其他輔助信息的前提下，BasicLSA比PHRAP(1999)及CAP3都有明顯的改進．這也證明了原型算法的有效性和可行性．第47頁/共57頁基于原型算法的優(yōu)化

原型算法的運行速度過慢，以至成為整個系統(tǒng)運行的瓶頸，從而導(dǎo)致該系統(tǒng)無法滿足實驗和可能的進一步應(yīng)用的需要．這也成了將原型系統(tǒng)實用化必須解決的問題.提速優(yōu)化的目標(biāo)是在不改變運行結(jié)果的前提下，提高系統(tǒng)運行速度．第48頁/共57頁“鄰域剪枝”

----(NeighborhoodPruning)

盡可能通過排除那些目標(biāo)函數(shù)值低于當(dāng)前值的鄰接點來縮小搜索空間．這里的關(guān)鍵在于找到了一個關(guān)于能夠增加目標(biāo)函數(shù)值的變換的必要條件．第49頁/共57頁Aftertransposition“領(lǐng)域剪枝”（NeighborhoodPruning）方法Befortransposition第50頁/共57頁

在搜索過程中，每一個transposition將導(dǎo)致目標(biāo)函數(shù)改變△overlap，Aoverlap=Overlap4+Overlap5+