基于高通量测序数据的结构变异检测方法的研究

张雨豪　王亚东

摘要：随着高通量测序数据技术的发展，人类全基因组的测序成本在不断降低，测序速度也有了较为显著地提升。运用生物信息学的手段处理这些海量基因组数据的需求也越来越迫切，而对于基因组结构变异的检测更是这个领域的核心内容。由高通量测序数据特征入手，介绍了当前主流的生物信息学结构变异检测方法，并阐述了有关基因组结构变异检测结果的评测指标和手段，最后，结合个人基因组的发展，对于该领域未来的发展提出了改进建议。

关键词：高通量测序； 结构变异检测； 生物信息学

中图分类号：TP391 文献标识码：A 文章编号：2095-2163（2013）05-0035-04

0引言

随着人类基因组计划的宣告完成，对于人类基因组海量数据的研究工作也逐步拉开了序幕，这给生物信息学的发展提供了很好的发展机遇，同时也带来了诸多挑战。之后的千人基因组计划更提供了大量第一手的人类基因组数据，这些数据既可以为生物学数据处理提供原始输入，又能为处理生物学数据所得的结果提供了良好的验证。

当利用高通量测序数据来检测结构变异时，主要有以下几种思路。第一种是单纯依靠覆盖率信息的方法，这种方法是最早提出检测结构变异的方法，现在已很少单独利用。第二种主要是依靠双末端测序数据中非一致序列并通过聚类来发现结构变异信息，这种方法很难发现具体的机构变异位点信息。第三种方法是利用split read来精确发现结构变异，这种方法可以精确发现结构变异信息，但是重复序列对其影响很大。现在大多数结构变异检测软件都会集成整合上述几种方法，取长补短，并会相应地构建一套独特的数据筛选处理流程，运用更快捷更有效的算法，由此而不断提高基因组结构变异检测的能力。

1高通量测序数据介绍

1.1高通量测序技术的介绍[HTSS]

对于人类基因组的全测序技术是解决基因组生物信息学的一个至关重要的前提。传统意义上最著名应用、最广泛的测序方法是Sanger测序法[1]，这种方法起源于上世纪70年代，已经过不断地改进而逐步趋于完善。而且，在2001年得到的第一条人类全基因组序列主要采用的就是这种方法，不过，这一过程是通过全球多个研究机构的共同努力，且耗费了数年时间花费巨资才完成。

随着对于更廉价、更快捷测序技术的需求激增，并经过该领域科学家的通力协作，高通量测序技术应运而生。高通量测序技术的出现极大地降低了全基因组的测序时间以及测序花费。

表1中显示了几种高通量测序技术的花费和优缺点，最后一列是第一代Sanger测序技术。从表中可以发现，虽然设备较贵，但是Illumina测序仪还是有相对便宜的价格和时间开销，并且由于Illumina测序仪可以使用户根据其需求生成不同的测序数据，因此，在结构变异检测中，原始数据大多数是通过Illumina测序仪得到的。

1.2双末端测序数据介绍

在Illumina测序仪的结果中主要会产生两种数据，一种是单末端数据（single end），一种是双末端测序数据（pair end）。这两种数据分别是根据不同的生物学手段得到的，其中双末端测序数据不仅有短序列（read）信息，而且还包含了插入距离信息，这对于同一组序列的位置关系提供了新的一种依靠和保证。在此重点介绍有关双末端测序数据的相关信息。

在双末端测序数据中，主要包含了相对基因组的上游序列信息、下游序列信息和插入距离信息，而且数据总是成对出现。由于在处理单末端数据时，主要通过短序列覆盖率信息和短序列自身信息来检测结构变异，在利用双末端测序时，不仅可以使用单末端数据中的信息，更能通过对于插入距离的信息来有效地检测结构变异，因此，在检测结构变异的时候大量采用了双末端测序数据。

2基因组结构变异类型介绍

随着人类基因组测序技术的进步，全基因组的数据每天都以海量的规模在增长。即使是两个不同人种的同性个体，其基因组之间的差别也是相当小的，虽然比例非常低，但是由于人类全基因组有30亿碱基序列，所以其数目仍是非常可观的，也正是这些差别导致了人类所有个体之间的万千差别。因此，开展这些差异的研究对于无论是疾病、或是医学等其他领域都有着至关重要的深远意义。

将参考基因组作为比对依据，由此得到的差异信息主要分为两类。第一类是SNP（单核苷酸多态性）；第二类是结构变异，在结构变异中较为常见的则是如图1所示的片段删除和片段插入。

一般来说，将某个体的基因组序列同参考序列进行比对，如果在一段序列区间内仅有一个位点不同，就将认定为SNP信息。如今的主要检测方法是基于贝叶斯估计进行分类，这种方法当1-5bp的结构变异时，就会产生一个基于统计学的较准确的结果，不过对于长序列问题的复杂度却会迅速增加，分析难度也会显著加大，此时该方法就不再可取。

3主流结构变异检测方法

相较于数量众多、性能优良的DNA序列比对工具，结构变异的检测工具一方面由于其发展起步较晚的影响，另一方面则由于结构变异事件的情况相对SNP和DNA比对更为复杂，因此，直到双末端测序数据的大量使用才出现了很多基于双末端测序数据的方法。

针对双末端测序数据，当前主流结构变异检测工具都已进行了较为充分的利用。SVseq2[2]是一款在低覆盖率情况下，主要通过对双末端测序数据中产生的split read和非一致的序列的分析来精确检测结构变异的工具。Pindel[3]是一款主要根据双末端测序数据来检测大长度的片段删除和中长度的片段插入的结构变异的软件。Delly[4]是一款使用短插入距离双末端测序数据、破碎序列片段数据相结合，既能检测基因组片段删除、片段倒置、片段连续重复，又能在平衡的基因重排数据中检测倒易移位事件（reciprocal translocations）的软件。PRISM[5]是一款既采用了双末端测序数据的聚类分析又使用了破碎片段序列的比对分析结果，并使用改良的Needleman-Wunsch算法而以精确到位点的方式来检测片段删除、片段插入为主的结构变异信息的软件。

其中，可利用的双末端测序数据都是经过BWA等[6]软件比对之后的SAM格式文件。文中将可利用的数据主要分为两类：非一致短序列对（discordant pair）和单映射双末端测序数据（hanging pair）。如果这两个序列片段的映射距离被认为是在插入距离的可接受范围内，而且两个片段的朝向都没有发生改变，即可认为这种序列对为一致的序列对（concordant pair），该种序列在绝大多数情况下均不会被认为覆盖了一个结构变异。除此之外，其他的双末端测序数据，无论是序列朝向问题、插入距离问题或者CIGAR值异常等问题发生时，均可认为产生的是非一致的序列对（discordant pair）。除此之外，一种特殊情况，就是双末端测序数据中仅有一个序列片段比对到参考序列上，而另一个却未能比对到参考序列上，由此将没有CIGAR值，这类特殊的序列可称为单映射双末端测序数据对[7]。

在此，以最近发布的一篇检测结构变异信息的软件PRISM为例，来讲述检测结果变异的大致流程：

（1）数据筛选。DNA全序列比对是检测结构变异的前提条件，因此输入是标准SAM格式文件。而其中的大多数序列不会覆盖结构变异信息，因此对于discordant pair和hanging pair的筛选工作将是后续的研究基础。[JP2]这里主要还是根据对SAM格式中的CIGAR值、插入距离、序列对正反朝向这三个方面的判断，其中之一有异常的状况，就将其归类并分别输出。

（2）数据聚类。虽然各个软件在聚类时采用的方法不尽相同，比如：PRISM采用了CNVer[8]作为聚类工具，并采用贪心策略将相似映射距离和朝向的序列对进行划分聚类，但却基本都采用了一种聚类手段以实现对数据进行更好地划分。

（3）Split read比对。通过上述两步得到了Split read数据，不同软件采用各自的算法将split read重新比对到参考序列上，再结合上一步的聚类结果，共同判断是否可以支持一个结构变异信息，并分别将其完整记录。

（4）发现结构变异信息。大多数软件会根据序列的质量数，并结合已经设定的可支持序列数的具体值来判定一个结构变异事件。如PRISM一般采用5作为支持一个结构变异信息的短序列最小值，由此一个结构变异事件只有被5个或者更多短序列支持的时候，才能判断是发生了结构变异。但Pinel却因其设计本就是为低覆盖率数据而实现准备的，因而通常选用2作为支持结构变异的短序列最小值。这些参数既跟数据的选用有密切关系，又与具体采用的聚类方法以及split read比对算法有关，因此往往需要不断调整以获得最理想数据。

4结构变异评测手段

在评测中，主要关注在一定覆盖率情况下该软件发现结构变异的精确率和召回率。首先，可认为包含实际存在的结构变异信息的集合为答案集，在模拟数据中，答案集可以

通过人为手段实现验证，并通过对于植入的结构变异信息的监控和追踪而对结构变异事件实时记录。同样，也可认为由结构变异检测工具报告的结构变异信息集合为结果集。共同集则表示既认定其存在、而又被结构变异检测工具报告的共同结果，三者之间的关系如图2所示。

其中，召回率=共同集/答案集；精确率=共同集/结果集。

文中以40x覆盖下chr1染色体上的模拟数据为基准通过使用PRISM、Pindel两种软件主要针对片段删除和片段插入进行检测。首先，人为地将Venter的结构变异信息植入hg18版本的1号染色体中，通过生成的序列与参考基因组相比较生成有关结构变异信息的答案集，再通过结果变异检测工具生成结果集，其后根据序列名称间的对应关系，由此而形成评估结果。其具体流程如图3所示。

这里分别使用15x和40x覆盖率的chr1染色体作为输入，对Pindel和PRISM两款结构变异检测软件的精确率和召回率进行比较。

结合图4和图5可以发现，当覆盖率较高的时候，主流的结构变异软件都具有良好的召回率和精确率。并且随着覆盖率的降低，精确率有小幅度的增加，召回率却有较明显的减少，因为其中包含了许多小型结构变异信息。这些信息采用split read比对的方法进行实验时会有一些噪音数据出现，所以，当剔除小型结构变异信息时，这些软件所能检测的数目会有所降低，但是结果却会有更好的精确率。

5基因组结构变异发展方向

虽然当今大多数软件在一定输入下都能够取得较满意的结果，但是，采用split read比对结合聚类的方法仍有一定的局限性。例如：对于被测对象序列中的大于插入距离的片段插入，现有的各类方法都是很难检测出来的。这种现象主要是因为，没有一个双末端测序序列能够覆盖这个片段插入事件。此外，由于聚类方法的采用，将很难精确发现一个结构变异事件的具体精确位点。而且，现在的方法主要是针对片段删除和片段插入的情况，对于其他结构变异信息却还未得到很优异的结果。

6结束语

随着高通量测序技术的迅猛发展，对于人类基因组的研究正在逐步迈入后基因组时代，对于测序数据的处理工作正变得尤为重要，利用生物信息学手段检测人类基因组结构变异这一核心领域的研究也将越来越深入和细致。本文从测序数据的介绍入手，对主流结构变异检测软件的流程和方法进行了分析，并且根据PRISM和Pindel在不同覆盖率情况下的输出结果对检测工具的性能评估进行了叙述。