基因组高通量测序数据结构变异识别的方法分析

　　0 引 言
　　遗传变异是生命的基本特征，遗传变异与表型差异之间的关系，是现代生物学的一个基本问题。由基因决定生物体的遗传特征和主要个体差异的观念正在逐渐改变，过去几年的许多研究显示，基因组中大尺度的结构变异（structural variation， SV）与个体的表型差异和疾病等有一定的关联[1]。近年的研究更进一步，很多疾病和表型差异与基因组的结构变异有关系，人们已经逐渐认识到结构变异的普遍性和重要性[2-3]。
　　1 结构变异及其重要作用
　　基因组结构变异通常指基因组DNA分子上整片的片段差异，包括保持数量平衡的倒位（inversion）和易位（translocation）变化以及导致数量失衡的插入（insertion）、缺失（deletion）和重复（duplication）变化，其中数量失衡的变化也称为拷贝数变异（copy number variation， CNV）。一般认为结构变异的片段，长度约在一千到几百万个碱基之间，一致性在90%以上。2004年，全球范围内数个“人类基因组计划”研究基地意外地发现，表型正常的人群中，不同个体间在某些基因的拷贝数上存在差异，结构变异的发现及其对基因杂合性的重要作用，彻底改变了人们对基因型研究的认识[5]。
　　个体基因组之间的差异，除少数SNP和短序列的插入/缺失（indel）变化之外，大部分以结构变异为主[6]。相关研究也证实，在玉米自交系中存在异常高的结构变异[7]。虽然普遍认为当基因组规模相对较小时，发生结构变异的数量会较少，但对大豆的研究却特别发现，相对于大豆基因组的大小和其中SNP的数量来说，结构变异的发生率较高[8]。而且结构变异在大豆的不同染色体间的分布也不均匀，在基因分布较密的区域发生较多，并且与植物的抗性表现有一定的关联[9]。
　　2 相关研究
　　最近几年随着下一代测序技术中高通量测序平台在应用上的推广普及，基因组序列数据的获取速度越来越快。这些平台在测序中都会产生大量的相互交叠的对端读片段（paired-end read）。被测序的DNA分子将打断为很多碎片，其插入长度约2k～150kbp，设备仅测每个碎片两端长度为25～1kbp的对端读片段。使用这样的序列数据识别结构变异，在长度和精确度上均会高于芯片技术的使用效果[10]，因而适用于各种类型结构变异的识别以及新类型的发现。采用测序技术的结构变异识别可以使用重头组装和重测序两种方法。其中，重头组装方法最为理想，但由于目前测序技术的限制，完整的全基因组组装结果并非最优，不能很好地满足先组装、再识别的要求；而重测序方法则更加实际，可将得到的大量读片段定位到参考基因组并实施比较分析。读片段定位（read mapping）是序列比对的基本任务，对其研究已可堪称深入，但通用的读序列定位算法应用于结构变异检测时，却没有充分考虑断点（breakpoint）处的比对，从而会导致断点预测不准确甚至结构变异检测失败[11]。
　采用重测序技术的结构变异分析方法主要有三种，具体来说，分别是：
　　（1）覆盖度（read depth或depth of coverage， DoC）分析法，根据在参考基因组上读片段的覆盖度，判断是否有拷贝数变化的结构变异，如CNV-seq[12]等，这类方法无法识别拷贝数不变的倒位和易位变化；
　　（2）对端定位（paired-end mapping， PEM）法，根据对端读片段插入距离的一致性，判断是否发生结构变异，如VariationHunter等；
　　（3）读片段分割法（split read），利用读片段在参考基因组上定位时的断点判断结构变异，如Pindel[13]，但由于NGS测序的读长较短，不利于发现较大的结构变异。也有将这些信息结合的方法，如CNVer[14]等只关注了CNV一种结构变异；Break- Dancer[15]和GASV[16]使用了高置信度读片段，忽略了与重复区域有关的低置信度读片段；inGAP-sv和Gnome Strip[17]在精确性和特异性两方面都较好，但却都基于二倍体人类基因组设计和实现。因此，目前采用NGS数据进行基因组结构变异的检测与分析方法还不够成熟，缺少完整的分析、发现和评价体系。
　　3 结构变异识别方法
　　本文的基本的结构变异识别方法如图1所示，其原理可概括为：以目标物种基因组序列为参考，使用多个供体全基因组高通量末端配对测序读片段数据，采用PEM的方法识别结构变异，并利用对端信息进行读片段分隔方式的断点预测。本文设计基于单个参考基因组和多供体基因组重测序数据的结构变异识别方法。从识别方法的角度，本文可以分为三个主要方面：数据预处理，结构变异识别和精确的断点预测。下面将给出其具体实现。
　　3.1 数据预处理
　　数据预处理的主要任务是读片段定位，不一致和歧义读对（discordant and ambiguous pairs）的过滤，产生的中间结果为不同类型的候选读片段集，用于支持可能的多种结构变异的识别和断点精确预测。在此，对其过程解析可做如下展现。
　　（1）读片段定位
　　在参考基因组上定位读片段，对每个供体的读片段分别向参考基因组上定位，相应得到各自的定位结果。定位时除需考虑一般定位算法涉及的测序错误、soft clipping等问题之外，还需要考虑两方面问题。具体描述为：
　　① 多处匹配：若同一个读片段可能在多处得到较好的匹配，则可能因为有拷贝数变化，需采用基于种子及扩展的比对算法，保存并返回多个高分比对位置；
　　② 对读片段分割的惩罚：若某条读片段刚好在断点处，由于对空位的罚分，可能造成错误定位或无法定位，因此要调整比对算法的罚分策略，允许一处不罚分的长空格。
　　（2）不一致读片段对的过滤
　　由测序技术可知，同一碎片双端测序的读片段对，在参考基因组上的定位如若符合一定的方向和距离要求即称为一致的（concordant），否则称为不一致的（discordant），如图2所示。不一致的原因之一就是由于参考基因组和供体基因组之间的结构差异，主要分为三类，将其综合概述如下：
　　①方向异常的，可能因为倒位，利用定位的方向可筛选出；
　　② 距离异常，可能对应插入、缺失或易位，根据测序参数中平均插入长度L和标准差σ来处理，当距离超过L+3σ错误！未指定书签。和小于L-3σ错误！未指定书签。的都视为异常并筛选出；
　　③ 悬挂对（hanging pairs），即对端读片段中有一条未能成功定位的，可能对应较大的结构变异，由定位结果可直接 筛选出来。
　　（3）歧义读片段对过滤
　　若一条读片段在多处均有较好定位，又可将其称为歧义定位对。可能由于潜在的结构变异造成，同样按照类似不一致读对过滤的方法进行筛选。只是要依据在第（1）步中读片段定位算法所提供的结果，即将歧义读片段对的多处匹配结果统一、完整地实现返回。
　　（4）合并多供体的候选读片段
　　由于文库制备和测序实验的不同，每个供体测序的基本参数也是不同的，如插入长度、读片段长度等，所以读片段的定位和过滤即需分别设计运行。而读片段的定位结果参考同一基因组，因此可将过滤得到的候选读片段按类型合并，并保存供体来源信息，并在发生分歧时提供分类依据。而相应于悬挂对，则由于只有一条读片段能够成功定位，需要对未成功定位的读片段连同歧义读片段进行聚类，以发现较大的结构变异。
　　3.2 结构变异识别
　　（1）基本结构变异识别
　　由结构变异导致的不一致读对如图3所示。对于基本的插入、缺失和倒位等，需根据定位的座位进行分类并计算读对共同的比对得分。对每个座位，若有3组以上来自同一供体的同类读对，则认为是合法的结构变异；或者，另有多个供体且每个有2组以上的同类读对，也将认定为是合法的结构变异。
　　（2）读片段覆盖度分析
　　利用多供体数据，在提高等效的覆盖度方面具有显著优势。覆盖度分析的原理如图4所示，当覆盖度发生异常时，可能因为基因组的重复片段或存在较大的结构变异。分析定位到参考基因组上读片段的数量，就需要估计全部成功定位的读片段在参考基因组上的分布和特征，再根据总体及局部的分布以判断是否呈现有明显的变化，并在一定大小的窗口内计算读片段覆盖度，而当信号发生相对显著的增加和减少时，即可判断有结构变异发生。
　　覆盖度与对端读片段信息的结合分析，可发现较大结构变异的原理如图5所示。具体地，由于几个对端读片段的距离较大并且同时此处的覆盖度超过正常值。在覆盖度异常的区域，寻找方向不一致的对端读片段，并试图利用复制该区域的方法，将不一致的对端读片段一致排列，以发现结构变异。
　　（3） 复杂结构变异识别与重复片段的分析
　对于超过插入长度的结构变异，会导致悬挂对和歧义定位。因此对这些读片段，使用贪心算法将距离和方向类似的片段实现聚类后，即需进行基于de Bruijn图的拼接。在此图中对数据错误和重复序列进行处理，并根据路径与参考基因组实行比较，从而分析可能存在的结构变异。如图6所示，利用多供体信息，可以相互验证路径中的环路，是否是因为重复片段。若构成路径的片段相对供体来源是均匀的，可以认为是在供体中普遍存在的重复序列；否则即认为是由于结构变异，再与参考基因组相比较，并确认结构变异的类型和位置。
　　3.3 结构变异的断点预测
　　断点预测需要读片段分割分析和对端读片段分析相结合，来自多供体的对端读片段构成差异较大，更利于断点位置的精确预测。如果有多个读片段在断点处的一侧发生完美匹配，则可以通过该组读片段直接预测断点位置，误差在读片段间交叠的长度范围内；而如果在读片段覆盖范围之外，则需要使用对端读片段的插入长度进行分析。因为插入长度L有一定误差范围，假设在Lmin和Lmax之间且对端读片段的左右断点分别是（x， y），那么可以估计断点发生的位置（a， b）满足 Lmin ≤ （a-x） + （y-b） ≤ Lmax，即（a， b）落在x， y， Lmin， Lmax构成的梯形区域当中。若有多个对端读片段同时支持该区域的结构变异，则该断点所在位置需要在多个梯形的交叠区域中，如图7中分别为形成插入和缺失时的情况。此外，为了能够处理超过插入长度结构变异的断点，利用断点两翼序列的特点，根据断点处的微同源性，分析已知结构变异断点处的微同源特征和两翼序列的特征，再进行断点预测的。
　　4 实验结果及分析
　　本文实验利用模拟数据生成供体的测序片段，根据dbVar对已知结构变异大小和分布，在大豆1号染色体上构造了插入和缺失类型的结构变异，如表1所示。由于易位、倒位和重复事件现有数据中也很少，而且通常各种软件对这些事件识别的准确率都不高，所以没有进行模拟。本文即比较选择了主流的几个结构变异识别工具，分别是BreakDancer、Pindel和CREST。
　　对缺失事件的识别结果如表2所示。综合该结果可知，总体而言，Pindel、BreakDancer和本文方法的识别结果各有优势，不可概述；分别地，Pindel和本文方法在小型的缺失事件上结果较好，而在100bp以上的范围SVdetect、BreakDancer和本文方法的结果则更佳。其中100～1000bp的两组结果中，本文方法识别效果较好，而在100～500bp和 > 1 000b两组结果中当属本文方法的效能最优。主要原因是对于利用split-read方式的识别工具，其大小产生split-read的数量较少，而对于不一致读片段对的识别工具其长度却恰好落入插入长度的误差范围内。
　　综合分析实验结果，本文方法在不同的条件下均可以得出较好的识别结果。本文方法能够准确识别出大部分的结构变异事件，因而可以用于NGS数据中结构变异的识别。
　　参考文献：
　　. Nature， 2010， 464（7289）： 713–20.
　　[2] ABECASIS G R， ALTSHULER D， AUTON A， et al. A map of human genome variation from population-scale sequencing[J]. Nature， 2010， 467（7319）：1061–73.
　　[3] 周雪崖， 张学工. 基于拷贝数变异的遗传关联研究[J]. 科学通报， 2011， 56（6）：370-182.
　　[4] TREANGEN T J， SALZBERG S L. Repetitive DNA and next-generation sequencing： computational challenges and solutions[J]. Nature reviews. Genetics， 2012， 13（1）：36–46.
　　. Nature reviews. Genetics， 2006， 7（2）：85–97.
　　[6] Le ROUZIC A， BOUTIN T S， CAPY P. Long-term evolution of transposable elements[J]. P roceedings of the National Academy of Sciences of the United States of America， 2007， 104（49）：19375–80.
　　[7] SPRINGER N M， YING K， FU Y， et al. Maize inbreds exhibit high levels of copy number variation （CNV） and presence/absence variation （PAV） in genome content[J]. PLoS genetics， 2009， 5（11）： e1000734.
　　. Plant physiology， 2011， 155（2）：645–655.