基于支持向量机的基因组拼接分叉结构处理方法

摘 要：基因组拼接过程中通常会出现分叉，使拼接变得困难。测序碱基错误（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.