前前后后接觸了一些基因組和轉(zhuǎn)錄組拼接的工作，而且后期還會持續(xù)進行,。期間遇到了各種各樣莫名其妙的坑,，也嘗試了一些不同的方法和軟件，簡單做一個階段性小結(jié),，本篇是原理部分,，下周的同一時間更新實戰(zhàn)部分。

6200 字,，約15分鐘,，文|思考問題的熊

拼接可以分為基因組和轉(zhuǎn)錄組拼接，基因組拼接對數(shù)據(jù)量和測序深度要求更高，而轉(zhuǎn)錄組用平時的 RNA-seq 數(shù)據(jù)就可以,。目前手里更多的是 RNA-seq 數(shù)據(jù),，所以做的也是和轉(zhuǎn)錄組拼接相關的內(nèi)容。 

無論拼基因組還是轉(zhuǎn)錄組,，歸根結(jié)底都是拼 DNA 序列,。拼接最簡單的邏輯就是一個由讓由長變短再由短到長的過程，基因打斷成一個個片段進行測序,，生成測序數(shù)據(jù),，然后再用reads拼成contigs再拼成scaffolds。

從 reads 到 contigs 的過程中,，需要進行多序列比對并將一致的 reads（consensus sequence）拼接起來來生成contigs，其中最大的問題是基因組上存在大量重復區(qū)域,，會對拼接帶來困擾,。Scaffolds數(shù)據(jù)則是通過pair end reads 信息來判斷contigs 的順序、方向和相鄰 contigs 之間的缺口（gap）大小來生成的,。從contigs到scaffolds是一個排序和定向的過程,。

關于拼接，Biostar 大牛 Istvan Albert在一個回答下面評論道：

目前常用的拼接算法都是基于數(shù)學中的圖論思想（Graph theory）產(chǎn)生,其中圖論中的兩個點表示兩個read,，而兩點之間的連線表示兩條read的重疊區(qū)域,。拼接要做的事情就是在所有的路線中找最優(yōu)解，類似于小時候玩過的一筆畫問題,。

如上圖左所示,，一個簡易的基因組產(chǎn)生了ABCD……若干read，在理想的情況下我們可以根據(jù)所有reads彼此的重疊區(qū)域重構(gòu)出圓圈所表示的基因組,。如果簡化成圖來表示則應該是上圖右的所示,，所有黑連接的是正確的基因組，但實際情況是基因組有很多區(qū)域比較相似,，以至于序列間會產(chǎn)生本沒有的聯(lián)系（如紅線所示）,。

在進行上述分析過程時，需要把所有的reads都進行比對,，以便找到重疊區(qū)域,，這個步驟非常耗費計算資源。Graph畫出來之后的問題就是如何從中得到最優(yōu)路徑,，即從有多種reads組合方式中找到從合適的一個形成contig,。

下面是尋找最優(yōu)路徑的常用算法。

貪婪算法（貪心圖）是早期提出的拼接算法,。首先選定初始read, 然后找和其重疊區(qū)域最高的read進行延伸,，直到拼接后的read兩端都不能再進行擴展為止。

每一次都是從最優(yōu)匹配開始,，然后次優(yōu)匹配,，到不能匹配時停止,。這樣一來，貪婪算法通常會得到局部最優(yōu)解,，而不是全局最優(yōu)解,。因此，這種算法在遇到重復序列時會出現(xiàn)比較大的問題,。

A greedy assembler compares all pairs of read fragments continually and replaces any pair with sufficient overlap between the edges with a combined sequence. The assembly completes, when the edges of the remaining sequences do not have any significant overlap.

OLC圖算法主要是用來針對長reads序列拼接,，如一代測序數(shù)據(jù)（三代測序數(shù)據(jù)），簡單理解就是把測序產(chǎn)生的長序列用彼此之間的overlap區(qū)域連接起來,。對于數(shù)據(jù)量很大的數(shù)據(jù)或者全基因組數(shù)據(jù)來說,，形成的olc圖非常復雜，會消耗大量內(nèi)存,。

OLC算法共有三步：

De Bruijn 圖 是目前最常用的二代測序拼接算法。比較流行的拼接軟件如 Velvet,、Abyss 和 SOAP denovo 都使用該算法,。

與OLC不同之處在于，這個算法將已經(jīng)非常短的reads再分割成更多個kmer短序列（k 小于reads 序列的長度）,，相鄰的kmers序列通過（k-1）個堿基連接到一起（即每次只移動一個位置）,，進而降低算法計算重疊區(qū)域的復雜度，降低內(nèi)存消耗,。

這里的kmer首先不能太短,，比如2個堿基肯定拼不出來基因組，它的長度既需要能夠使其攜帶足夠的基因組的信息，也要短到可以進行后續(xù)的錯誤矯正,。除此之外,，一個read中小的片段被分割之后還不會丟失原來reads 的前后位置信息。

總體而言,，該算法將reads打斷成長度為K的核酸片段,，再用Kmer間的overlap關系構(gòu)建DBG，最后通過DBG得到基因組序列,。

拼接步驟通常包括：

• 構(gòu)建DBG圖,，將read分割為一系列連續(xù)kmer

如下圖，1和2兩個序列會拼接成一條長的序列,，在引入第三條序列后會出現(xiàn)一個圓環(huán),。

下圖是一個最簡單的DBG拓撲結(jié)構(gòu)，兩球一線（一進一出）代表相鄰的兩個kmer,，圓圈則代表有多種連接方式,。

下圖是三種常見的DBG結(jié)構(gòu)，分別代表了拼接過程中的不同情況,。

• 合并DBG圖

合并路徑中出度入度唯一（one incoming and one outcoming ）的節(jié)點，去除段末端,，低覆蓋度節(jié)點和泡狀結(jié)構(gòu),。

• 構(gòu)建contig

尋找最優(yōu)路徑（經(jīng)過每個節(jié)點且僅經(jīng)過一次），最優(yōu)路徑對應的堿基序列構(gòu)成一個contig

• 構(gòu)建scaffold

通過PE reads 位置信息確定contig之間的相對位置和方向,，組裝contig,，填充contig之間的gap，得到scaffold序列,。

在拼接相關的文章中，kmer是出現(xiàn)頻率非常高的一個詞,。而kmer在整個生物信息分析過程中的用處也是非常之多,。

k值越大可辨別更多的小重復序列，越容易把DBG轉(zhuǎn)換為唯一的序列,，但得到的拼接過程含有更多的gaps,；小的k值對應的DBG能夠得到較好的連通性，但是算法的復雜度會提高，repeats序列處理會更復雜,，增加了錯拼的可能性,。

在拼接數(shù)據(jù)預處理軟件khmer的文獻中有這樣一段關于kmer和內(nèi)存大小與處理結(jié)果關系的描述：

The interaction between these three parameters and the filtering process is complex and depends on the data set being processed, but higher coverage levels and longer k-mer sizes result in less data being removed. Lower memory allocation increases the rate at which reads are removed due to erroneous estimates of their abundance, but this process is very robust in practice

kmer越大需要的內(nèi)存就越多，所以計算機的內(nèi)存大小也會限制kmer的取值,。這里需要說明的是,，輸入數(shù)據(jù)的多少不會影響memory用量，但是輸入數(shù)據(jù)的錯誤越少,，占用的內(nèi)存也就越少,，假設所有測序數(shù)據(jù)都沒有任何錯誤，那么DBG的大小并不會因為測序深度的增加,，因為不需要將因為幾個堿基不一致的kmer存入到DBG中（下一部分會具體提到）,。至于需要多大的RAM則取決于DBG的大小和組裝基因組的大小。

另外,，在拼接的過程中盡量避免使用偶數(shù)kmer,，否則容易是kmer產(chǎn)生回文序列，特別是在鏈特意性的數(shù)據(jù)中,。

在平時分析中,，一般會設置一個kmer的梯度（21，23,，25,，27,2931），來解決DBG算法loss of read coherence的問題,。然后從中選擇最好的結(jié)果,。另外,，還有一種說法是在進行拼接過程時，kmer應該選擇read長度的1/2到2/3大小，否則可能拼接出過多的Contig,。這一點,，也可能是我們平常使用trinity拼接時拼出Contig 過多的原因,，trinity的默認拼接大小是25,。上限是32？（有待確定）,。如果kmer有上線,，是否也可以考慮在預處理的時候，處理的力度大一點,，把序列截短一些,？

在實際情況中，拼接往往是在覆蓋度不均勻且含噪聲的數(shù)據(jù)中進行,，這為拼接帶來了三個方面的困難：

1. 增加了大量假kmer從而提高了對緩存的要求

2. 錯誤的reads通過增加tips,，bubbles和corss-links等改變了DBG的結(jié)構(gòu)

3. 不統(tǒng)一的read覆蓋度使得拼接參數(shù)對拼接結(jié)果有非常大的影響

在一次測序得到的數(shù)據(jù)中,，kmer matches 的數(shù)量和測序深度以及read長度相關。假設在完全沒有測序錯誤的數(shù)據(jù)中,，read長度是100,，測序深度是50X，選取kmer值為21,，那么一個只匹配基因組一次的kmer 出現(xiàn)的次數(shù)應該是$(100-21+1)*50/100=40$（基因組的每個位置被測了50次,，100bp的read有80個21bpkmer）,如果匹配基因組兩次的kmer應該出現(xiàn)80次。因此,，峰值應該在40,，80的位置有一個小峰。

但是,，當存在測序錯誤的時候,，會在匹配1次的位置出現(xiàn)大量的kmer,就是由于測序的誤差導致的。為什么說是測序錯誤呢,，因為在50x的測序中只出現(xiàn)了一次,，如果一個read中有一個堿基錯了，那么這個read就會產(chǎn)生21個錯誤的21kmer,。更大的問題是,，隨著測序深度的增加，這樣錯的kmer數(shù)量也會增加,，

錯誤的reads 會為DBG引入三種類型的錯誤

tips

所謂tips指一個小分支,，下圖所示，我們有5條10bp的reads,，其中第5條有一個堿基測序錯誤。使用7kmer會產(chǎn)生11個節(jié)點,。如果用前4個read來拼是正常的,，一旦引入第五條就會因為一個錯誤的堿基出現(xiàn)一條錯誤的有三個節(jié)點的分支。

在大量拼接的過程中,，大量的read會匹配到正確的位置,，一小部分會比配到錯誤的位置，因此可以對錯誤的tip進行清除,。

bubbles

當read和kmer相比足夠長時,，錯誤可能出現(xiàn)read中間。此時會出現(xiàn)bubble的情況

需要注意的是,，除了測序錯誤以外,，可變剪切和snp以及插入缺失等也會導致tips的出現(xiàn)，因此增加了拼接的難度,。因此,，在進行拼接之間有必要對原始數(shù)據(jù)進行預處理,。

覆蓋度不均一

在實際拼接過程中，會去除一些低頻率的kmer,，這一操作在刪除了大量錯誤kmer的同時,，我們也不可避免的刪除了很多是因為覆蓋度低所以出現(xiàn)次數(shù)低的kmer。

綜上,，如果調(diào)低cutoff,，高覆蓋區(qū)域出錯的可能就高，但是低覆蓋區(qū)域的質(zhì)量提升,。kmer長度增加會使低覆蓋區(qū)域更加分散,，因為kmer的覆蓋度會因為低于設置的cutoff被刪除。

隨著kmer的增加,，分支會逐漸減少,，DBG會越來越趨向于線性。對于一個很大的kmer,可能就完全線性,，同時按照染色體分開,。拼接質(zhì)量通常會隨著k值的增加先變好然后再變壞，因為這個過程中存在兩種競爭性的過程,。一方面,，kmer的增加可以更好的處理重復區(qū)域，但是另一方面,，由于覆蓋度的原因,，kmer的增加在一些區(qū)域使得他們出現(xiàn)的次數(shù)越來越低直到低于篩選的閾值。所以,，kmer 的選擇非常重要,。

不同的拼接內(nèi)容需要不同的拼接策略，其原因如上圖所示,，即不同的數(shù)據(jù)產(chǎn)生的DBG結(jié)構(gòu)和覆蓋度不同,。

對于轉(zhuǎn)錄組拼接來說，如果假設一個轉(zhuǎn)錄組只有兩個基因而且這兩個基因沒有重復區(qū)域,，那只需要構(gòu)建兩個沒有關聯(lián)的DBG就可以了,。

對于實際數(shù)據(jù)來說，一個轉(zhuǎn)錄組包含了成千上萬個基因,。且大多數(shù)基因沒有overlap,，我們構(gòu)建的DBG實際上是眾多不相關grphs的集合，每一個圖都代表一個基因或者一個基因家族,。下面的示意圖顯示構(gòu)建了五個基因的圖并且也展示了其覆蓋度的差別,，其中基因2和3可能是兩個可變剪切或者一個家族非常類似的兩個基因。

通過上圖也可以發(fā)現(xiàn),，轉(zhuǎn)錄組拼接最大的問題在于kmer數(shù)目的高低很大程度上取決于某個基因的表達量高低,。如果把kmer cutoff 設置的高一些,，那么低表達的基因就可能拼不起來。

在基因組拼接過程中,，kmer coverage 通常是單峰或者是雙峰,，但是轉(zhuǎn)錄組拼接不同。轉(zhuǎn)錄組拼接中,，每一個基因的kmer分布可能是一個單峰,，peak的位置取決于基因表達量的多少。而整體的分布則是所有這些單基因分布的總和,。

轉(zhuǎn)錄組拼接時應該注意這樣量個問題,。首先，轉(zhuǎn)錄組拼接已經(jīng)不滿足DBG算法中覆蓋度均一的假設,，表達量非常高的基因可能主宰拼接的結(jié)果,，因此需要能調(diào)整這種覆蓋度的差異。另外,，轉(zhuǎn)錄組拼接不用太擔心低復雜度區(qū)域,，不會有很多重復區(qū)域出現(xiàn)。

轉(zhuǎn)錄組拼接常用軟件：trinity,；SOAPdenovo-trans(華大),；Trans-ABySS；bridger

基因組層面的拼接質(zhì)量,，一般會比較看重長度相關的指標,。比如最大長度、平均長度,、拼接后的總長度和 contig N50長度,。

Contig N50 指 reads 拼接后獲得一些不同長度的 contigs，將所有的 contigs 長度相加得到總長度,。然后將所有的 contigs 按照長度從長到短進行排序,，再將 contigs 按照這個順序相加，當長度等于 contigs 總長度的一半時,，最后一個加上的 contig 的長度稱為 contig N50。對于總長度不同的兩個拼接數(shù)據(jù),，直接對比N50 的數(shù)值沒有什么意義,。

對于轉(zhuǎn)錄組拼接而言，并不是越長越好,，我們更在意的是拼接的質(zhì)量,，方向和回帖率等等信息。如果我們在轉(zhuǎn)錄組拼接過程中使用了kmer=25這個參數(shù),，在拼接好后應該用拼接的fastq文件mapping回拼接好的轉(zhuǎn)錄組,，測試mapping效率,，這里推薦使用salmon軟件，需要注意salmon中的kmer應該和拼接時采用的kmer保持一致,。

另外,，transrate是一個專門用來評價拼接質(zhì)量的軟件。在后續(xù)的實際應用部分會有介紹,。

你可能還想看 思考問題的熊 其它好文

1. 貫穿你我一生的基因之謎
   

2. 轉(zhuǎn)錄組學習建議及入門指北
   

3. 四步把實驗室服務器搬回家
   

4. 鏈特異建庫那點事
   

5. ChIP-seq陰陽-正負對照
   

6. 來了,，總要留下些什么
   

7. 可能是最全的JBrowse基因瀏覽器介紹
   

8. 師妹 你要的都在這里了