基于支持向量機的基因組拼接分叉結構處理方法

朱曉　王亞東

摘 要：基因組拼接過程中通常會出現分叉，使拼接變得困難。測序堿基錯誤（sequencing error）是出現分叉的主要原因。針對分叉結構，研究分析分叉處的reads信息，建立SVM預測模型，提出基于支持向量機的分叉結構處理方法，取得了較好的效果。

關鍵詞：基因組拼接；分叉結構；支持向量機

中圖分類號： TP391.2

Research on Method for Dealing with Branches in Genome Assembly based on Support Vector Machine

ZHU Xiao， WANG Yadong

（School of Computer Science and Technology， Harbin Institute of Technology， Harbin 150001， China）

Abstract： There are usually some branches in genome assembly， which makes the genome assembly more difficult. Sequencing errors are the main reason. For branches， the paper analyzes the reads information at branches to build the SVM prediction models， and proposes the method for dealing with branches using SVM approach. Experimental results show that the SVM approach can obtain good performance.

Keywords： Genome Assembly；Branch； Support Vector Machine

0 引 言

基因組序列拼接（genome sequence assembly）是生物信息學領域的重要問題，測序產生的讀取片段（reads）經過序列拼接組裝，生成基因組的堿基序列。基因組拼接是研究新物種[1， 2]，以及研究已有物種的基因組結構變異[3]等的基礎。基因組拼接算法主要采用基于交疊的策略和基于De Bruijn圖的策略，去除圖中的分叉結構（branches），包括死路徑（dead ends）和氣泡結構（bubbles），如圖1所示。

圖1 基因組拼接過程中移除堿基錯誤[5]

Fig.1 Error removal scheme for genome assembly[5]

分叉結構是指圖中出度（outgoing）或入度（incoming）大于1的節點以及與之相鄰的邊所形成的結構。經過分叉節點的路徑（path）代表了不同的堿基序列。測序堿基錯誤（sequencing errors）是導致分叉結構的重要原因。由于測序數據中包含了大量的堿基錯誤（錯誤率約2%）[4]，而堿基錯誤將會導致拼接圖中的死路徑（dead ends）的出現，而大量的堿基錯誤將會使圖的結構復雜，拼接的最終結果往往很難產生唯一的基因組序列，而是一些組成基因組序列的子序列，可將其稱為contigs（或scaffolds）。基因組拼接的目標是構建出目標基因組的堿基序列，而分叉結構使基因組拼接變得更加復雜。基因組的拼接過程需要識別出由堿基錯誤導致的分叉，并將其移除，再選擇正確的路徑繼續contigs的擴展。本文分析分叉處的數據特征，提出基于支持向量機的分叉結構處理方法。

1 支持向量機

分叉對應于多條可選的路徑，每條路徑都有一定數量的支持該分叉的reads數據。reads數量少的路徑通常是由于測序錯誤導致的，而reads數量多的路徑最有可能是正確的，現有的拼接方法也都是基于這種思想，在遇到分叉的時候優先選擇reads數量多的路徑進行contigs的擴展。然而，在這些路徑都差別不大的情況下，這種擴展也可能是錯誤的，此種情況下即應停止擴展。對于一個分叉，其reads數量最多的路徑的擴展有兩種選擇，一種是繼續擴展，而另一種是停止擴展。這也各自對應于模式識別中的兩類分類的情況。支持向量機（Support Vector Machine， SVM）[6]建立在統計學習理論的結構風險最小原理基礎上的，根據有限的樣本信息在模型的復雜性（即對特定訓練樣本的學習精度）和學習能力（即無錯誤地識別任意樣本的能力）之間尋求最佳折中，以求獲得最好的推廣能力，在解決小樣本、非線性及高維模式識別中表現出許多特有的優勢，并能夠推廣應用到函數擬合等其他機器學習問題中。

支持向量機的基本思想是：首先，在線性可分情況下，基于原空間尋找兩類樣本的最優分類超平面。在線性不可分的情況下，加入了松弛變量進行分析，并通過非線性映射將低維輸入空間的樣本映射到高維屬性空間使其變為線性情況，使得在高維屬性空間采用線性算法對樣本的非線性進行分析成為可能，同時進一步在該特征空間中尋找最優分類超平面。其次，將通過使用結構風險最小化原理在屬性空間構建最優分類超平面，使得分類器得到全局最優，并在整個樣本空間的期望風險以某個概率滿足一定上界。

支持向量機中不同的內積核函數將形成不同的算法。目前常用的核函數主要有多項式核函數、徑向基核函數等。

多項式核函數：

（1）

徑向基核函數：

（2）

2 基于支持向量機的分叉結構處理方法

2.1 分叉結構及其特征選擇

由于測序數據中的堿基錯誤和基因組中的重復序列，當擴展contigs時，會導致拼接過程出現分叉。測序錯誤是導致分叉的主要原因，測序錯誤是隨機發生的，而測序深度通常很高（如>50x），對于基因組中的堿基，測序錯誤表現為出現的測序次數很低的隨機噪聲。分叉處的每條候選路徑分別對應一定數量的reads，分叉處reads數量少的路徑往往是由測序錯誤導致的，也就是reads數量多的路徑通常是正確的路徑，而reads數量少的通常是錯誤的路徑。處理分叉就是從這些候選的路徑中選擇正確的路徑，繼續進行contig的擴展，錯誤的路徑被移除，正確的路徑被合并到contig中。

研究將每個分叉處的最大reads數與次大reads數分別記為maxOcc和secOcc。通過實驗觀察，可以發現：如果maxOcc與secOcc相差很大（如maxOcc >> secOcc），maxOcc對應的路徑通常是正確的；相反，如果maxOcc與secOcc差別不大（如secOcc/maxOcc > 0.7），maxOcc對應的擴展有可能是不正確的，因此需要停止擴展并做進一步的檢查；而且，在重復區域contig的覆蓋度會升高，往往達到2倍以上，這種情況需要特別小心。因此，研究將以上的分叉信息抽取出來，建立SVM預測模型，用于引導contigs的擴展，如圖2所示。

圖2 Contig的末端的分叉信息

Fig.2 Branch information at contig end

研究將分叉信息記錄為特征（maxOcc， secOcc， covRatio， gapLen），其中，maxOcc和secOcc分別是支持候選路徑的最大reads數和次大reads數，covRatio是contig末端（兩個reads長度）的平均每堿基覆蓋的reads數與整條contig的平均每堿基的覆蓋的reads數的比值，gapLen是contig末端成功拼接到contig上的reads到分叉的最近距離。maxOcc與secOcc的差距越大，表明maxOcc對應的路徑正確的可能性就越高；gapLen越小，表明contig末端的reads數量越多，從而maxOcc對應的路徑就應該是正確的；covRatio越高，表明在contig的末端出現了重復序列，這時應停止contig的擴展。

2.2 支持向量機預測模型的生成

為生成SVM預測模型，本文使用廣泛被研究的模式生物大腸桿菌（Escherichia coli strain K-12 substrain MG1655）的基因組作為待測基因組，從NCBI （National Center for Biotechnology Information）網站（http：//www.ncbi.nlm.nih.gov/）下載其參考基因組（RefSeq： NC_000913.2），大小為4 639 675 bp （base pairs），包含一條染色體。研究中修改了模擬reads生成工具GemSim[7]，使其能夠輸出reads數據及其在基因組中的位置。使用修改后的該工具，生成40倍和100倍覆蓋深度的模擬配對reads數據，并用這些數據進行拼接。配對reads數據是指DNA分子的兩端分別進行測序得到的具有特定距離的reads數據。

在拼接的過程中記錄每個分叉的四個特征，將每個分叉視為四維（4D）空間的一個點。根據這些四維空間的數據點，通過SVM機器學習的方法繪制超平面，確定每個分叉是要繼續擴展還是應該停止。在拼接的過程中，根據拼接到contigs上的reads數據在基因組中的位置信息，通過contig與參考序列的比較，這些分叉可以被分為4類：正確擴展（correct extension），錯誤擴展（wrong extension），正確停止（correct stop），錯誤停止（wrong stop）。正確擴展和錯誤停止的分叉是應該被繼續擴展的，因此被標記為CONTINUE（正例）；相反，正確停止和錯誤擴展的分叉是應該被停止的，因此被標記為STOP（反例）。根據contigs上的reads的位置信息，可以確定每個分叉的正確擴展，用以標記這些樣本，正例標記為+1，反例標記為-1。

根據采集到的訓練樣本生成SVM預測模型，再基于SVM預測模型確定最大的堿基是否應該被擴展。研究使用Matlab （R2012b）訓練和測試SVM預測模型。支持向量機SVM的核函數使用Matlab自帶的3次多項式核函數 ，其中x， y是包含分叉信息的向量， 是向量x和y的內積（又稱點積），該預測模型用來決定分叉是否應該被擴展。

由于在拼接的起始階段主要使用單端reads數據擴展contig，直到contig足夠長時才使用配對reads數據。單端數據和配對數據的使用具有明顯的區別，如在contig的起始拼接階段，還未有reads能夠成功拼接到contig上，此時的分叉節點的gapLen將會很大，甚至達到100 bp （read長度），而配對數據的gapLen將難以接受如此大的gapLen。因此研究生成了兩個SVM預測模型，分別對應配對數據和單端數據，即在有配對數據時使用配對數據的SVM預測模型，無配對數據時使用單端數據的SVM預測模型。這兩個預測模型的生成方法相同，如圖3所示。

圖3 SVM預測模型的生成方法

Fig.3 SVM prediction model generation

研究首先生成單端數據的SVM預測模型。使用40倍的模擬數據，并將其看作單端數據。在拼接的過程中，一共記錄了96 897個分叉，其中94 973個正例樣本（應該被擴展的分叉樣本），1 924個反例樣本（應該被停止的分叉樣本）。因為SVM預測模型的訓練較為耗時，為降低計算的時間，即選擇隨機選取了其中的8 845個正例樣本和1 769個反例樣本，以此來訓練SVM預測模型。給定新的分叉樣本，SVM模型將被用于確定該分叉是否應該繼續擴展。進一步地，再將訓練之后的SVM預測模型用于全部的分叉樣本，得到正確預測的樣本96 667個，正確率99.76%。

然后，用同樣的方法生成了配對數據的SVM預測模型。使用100倍的配對數據，在拼接的過程中，記錄了136 637個分叉，其中136 539個正例樣本（應該被擴展的分叉樣本），98個反例樣本（應該被停止的分叉樣本）。研究過程隨機選取了其中的2 800個正例樣本和56個反例樣本，以此來訓練SVM預測模型。而且將訓練之后的配對數據的SVM預測模型用于全部的分叉樣本，得到正確預測的樣本136 539個，正確率99.93%。

3 實驗結果與分析

3.1 實驗數據

為評價SVM預測模型的性能，并將其與已有的分類方法進行比較，包括隨機森林（random forest）、k-近鄰（k-nearest neighbor，KNN），以及樸素貝葉斯（Na?ve Bayes）分類方法。研究使用模式生物大腸桿菌（Escherichia coli strain K-12 substrain MG1655）的基因組作為目標基因組（RefSeq： NC_000913.2），使用GemSim[7]生成模擬數據，覆蓋深度分別為40×和100×，分別使用單端數據和配對數據拼接，并在拼接的過程中記錄出現的分叉信息，如表1所示。

表1 大腸桿菌40×和100×模擬數據的分叉信息

Tab.1 Branches for 40× and 100× E.coli simulated reads

訓練集 測試集

正例數 反例數 正例數 反例數

SE_40× 8845 1769 94，973 1924

PE_40× 2935 59 56，439 71

SE_100× 1910 382 231，900 434

PE_100× 2800 56 136，539 98

按照同樣的方法，使用GemSim[7]生成模擬數據，覆蓋深度分別為50×、60×和100×，分別標記為D1，D2和D3，測試SVM預測模型在不同覆蓋深度的數據集上的性能表現。

3.2 實驗結果與分析

3.2.1 不同分類方法的比較

試驗中使用40倍和100倍測序深度的4組分叉數據集，比較SVM預測模型和其他分類方法的性能，這些分類方法包括隨機森林（random forest）、k-近鄰（k-nearest neighbor，KNN），以及樸素貝葉斯（Na?ve Bayes）分類方法，并且比較了SVM預測模型在選擇不同核函數時的性能表現。SVM核函數分別使用3次多項式核函數（polynomial kernel function） 和高斯徑向基核函數（radial basis kernel function） ，其中x， y是包含分叉信息的向量， 是向量x和y的內積（又稱點積）。

在此基礎上，又進一步比較了不同分類方法的分類錯誤率（mis-classification rate），即錯誤識別的樣本數量占總樣本的比例，結果如表2所示，多項式核函數用“SVM_poly”表示，高斯徑向基核函數用“SVM_rbf”表示，最優結果加粗表示。從表2可以看出，多項式核函數的SVM預測模型總體上具有最好的性能表現，其次是隨機森林和k-近鄰方法，而樸素貝葉斯方法的錯誤率最高，性能表現最差。徑向基核函數的SVM預測模型的分類錯誤率較高，是多項式核函數SVM預測模型的2倍以上。

表2 不同分類模型的錯誤率比較

Tab.2 Mis-classification rates for different classification models

Dataset SVM_poly SVM_rbf Random forest KNN Na?ve Bayes

SE_40× 0.24 % 0.76 % 0.29 % 0.39 % 4.58 %

PE_40× 0.25 % 0.73 % 0.18 % 0.20 % 2.68 %

SE_100× 0.22 % 0.57 % 0.24 % 0.23 % 1.43 %

PE_100× 0.07 % 0.34 % 0.12 % 0.10 % 1.61 %

通過將不同的模型應用于基因組拼接，發現40×單端數據（SE_40×）和100×配對數據（PE_100×）的SVM預測模型能夠生成最好的拼接結果，因此研究選擇這兩個模型作為最終要使用的SVM預測模型。

研究選擇SVM學習方法作為拼接過程中的分叉結構處理方法的原因主要有：（1） 多項式核函數的SVM預測模型能夠給出更加準確的預測結果；（2） SVM具有清晰的數學表達式，便于將已有的預測模型嵌入到所需要的拼接程序中。

3.2.2 SVM預測模型的性能表現

分別在拼接過程中記錄D1～D3的分叉信息，并采用之前生成的SVM預測模型進行預測，結果如表3所示。針對每個數據集，統計如下四類樣本：

（1） 真陽性（True Positives，TP）：正確擴展的分叉數量（Correct extensions）；

（2） 假陽性（False Positives，FP）：錯誤擴展的分叉數量（Incorrect extensions）；

（3） 真陰性（True Negatives，TN）：正確停止的分叉數量（Correct stops）；

（4） 假陰性（False Negatives，FN）：錯誤停止的分叉數量（Incorrect stops）。

由表3，就可以得到在這3組不同測序深度的數據集上，SVM預測模型的正確擴展的比例達到99.7%以上，同時SVM預測模型也一樣能夠確定停止的分叉情況，并且具有很小的錯誤率，表現性能良好。

表3 SVM預測模型的統計結果

Tab.3 Statistical results for SVM prediction model

Dataset Correct extensions （TP） Incorrect extensions （FP） Correct stops （TN） Incorrect stops （FN）

D1 （50×） 70 299 （99.70%） 60 （0.09%） 123 （0.17%） 26 （0.04%）

D2 （60×） 84 829 （99.74%） 46 （0.05%） 148 （0.18%） 25 （0.03%）

D3 （100×） 136 309 （99.82%） 48 （0.04%） 169 （0.12%） 27 （0.02%）

4 結束語

本文介紹了分叉結構及其特征，并提取該特征，提出了基于支持向量機的分叉結構處理方法，并利用支持向量機生成SVM預測模型，對給定的分叉結構是應該進行擴展還是應該停止進行預測，取得了很好的結果，準確率可以達到99%以上，對于處理由于測序錯誤導致的分叉具有很好的作用。

參考文獻：

[1] LI R Q， FAN W， TIAN G， et al. The sequence and De Novo assembly of the Giant Panda Genome[J]. Nature， 2010， 463（7279）：311-317.

[2] SWART E C， BRACHT J R， MAGRINI V， et al. The Oxytricha Trifallax Macronuclear Genome： A complex Eukaryotic Genome with 16，000 tiny chromosomes[J]. PLoS Biol， 2013， 11（1）：e1001473.

[3] LI Y R， ZHENG H C， LUO R B， et al. Structural variation in two human genomes mapped at single-nucleotide resolution by whole Genome De Novo Assembly[J]. Nat Biotechnol， 2011， 29（8）：723-730.

[4] SHENDURE J， JI H. Next-generation DNA sequencing[J]. Nat Biotechnol，2008， 26（10）：1135-1145.

[5] FLICEK P， BIRNEY E. Sense from Sequence Reads： Methods for alignment and assembly[J]. Nat Methods，2009， 6（11 Suppl）：S6-S12.

[6] V Vapnik. The Nature of Statistical Learning Theory[M]. 2nd ed. Berlin： Springer； 1995.

[7] MCELROY K E， LUCIANI F， THOMAS T. Gemsim： general， error-model based simulator of next-generation sequencing data[J]. BMC Genomics， 2012， 13：74.