CN105989246B - 一种基于基因组组装的变异检测方法和装置 - Google Patents
一种基于基因组组装的变异检测方法和装置 Download PDFInfo
- Publication number
- CN105989246B CN105989246B CN201510043893.2A CN201510043893A CN105989246B CN 105989246 B CN105989246 B CN 105989246B CN 201510043893 A CN201510043893 A CN 201510043893A CN 105989246 B CN105989246 B CN 105989246B
- Authority
- CN
- China
- Prior art keywords
- sequence
- sequencing
- preset value
- library
- sample
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 31
- 230000035772 mutation Effects 0.000 title claims abstract description 13
- 238000012163 sequencing technique Methods 0.000 claims abstract description 128
- 238000000034 method Methods 0.000 claims abstract description 37
- 238000002864 sequence alignment Methods 0.000 claims abstract description 18
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 12
- 238000003780 insertion Methods 0.000 claims description 11
- 230000037431 insertion Effects 0.000 claims description 11
- 238000003860 storage Methods 0.000 claims description 5
- 230000015556 catabolic process Effects 0.000 claims description 4
- 238000006731 degradation reaction Methods 0.000 claims description 4
- 230000004069 differentiation Effects 0.000 claims description 3
- 230000002159 abnormal effect Effects 0.000 claims description 2
- 230000011218 segmentation Effects 0.000 claims description 2
- 239000002585 base Substances 0.000 claims 8
- 239000003513 alkali Substances 0.000 claims 1
- 238000005516 engineering process Methods 0.000 description 17
- 239000012634 fragment Substances 0.000 description 13
- 238000002493 microarray Methods 0.000 description 11
- 238000012217 deletion Methods 0.000 description 7
- 230000037430 deletion Effects 0.000 description 7
- 238000011160 research Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 6
- 238000013459 approach Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 238000009826 distribution Methods 0.000 description 4
- 239000003814 drug Substances 0.000 description 4
- 238000001914 filtration Methods 0.000 description 4
- 238000001712 DNA sequencing Methods 0.000 description 3
- 239000008280 blood Substances 0.000 description 3
- 210000004369 blood Anatomy 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000003252 repetitive effect Effects 0.000 description 3
- 230000000717 retained effect Effects 0.000 description 3
- 238000005303 weighing Methods 0.000 description 3
- 108020004414 DNA Proteins 0.000 description 2
- 108091028043 Nucleic acid sequence Proteins 0.000 description 2
- 230000001186 cumulative effect Effects 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 230000007812 deficiency Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000003908 quality control method Methods 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 241000532838 Platypus Species 0.000 description 1
- 230000004308 accommodation Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000003766 bioinformatics method Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 230000035800 maturation Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000002773 nucleotide Substances 0.000 description 1
- 125000003729 nucleotide group Chemical group 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 239000000344 soap Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000011895 specific detection Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Abstract
本发明公开了一种基于基因组组装的变异检测方法和装置，所述方法包括：获取来源于梯度测序文库的测序读段序列；对测序读段序列进行过滤；将已经过滤的读段序列拼接成平均长度达到第四预设值的长序列；将拼接得到的长序列比对到参考基因组上；和对序列比对结果进行变异检测，获取序列变异。本发明的方法有效解决了变异检测“暗区”、长序列插入和复杂结构性变异的检测难题。
Description
技术领域
本发明涉及基因组学及生物信息学技术领域，具体涉及基于基因组组装的变异检测方法和装置。
背景技术
随着国际人类基因组计划的完成以及DNA测序成本的急剧下降，借助DNA测序技术实现临床个性化医疗的时代已经临近。在理想的情况下，期望的DNA测序技术应该是一种能够将细胞中的DNA一次性连续、完整地读取出来的技术。但事实上，即便是科技如此发达的今天也还是做不到，当前最先进的高通量测序平台，也只能将细胞中完整的DNA序列预先切割成一小段一小段，然后再在机器上进行测序，并且当前每次读取的序列长度最大也只能达到250个碱基对（英文简称：bp），但目前用得最广还是100bp和150bp这两种，而每个人的基因组总共约有30亿个碱基对。因此，在测序完成之后，实际上只是得到了一大堆非常零碎的小序列而已。
在这种情况下，为了对测得的人类基因组进行研究，往往需要先借助DNA序列比对软件（如BWA），将所有这些小序列片段都比对到人类基因组参考序列上，定出每一个小序列在基因组上的位置，之后才能进行下一步的分析和研究。不过，考虑到人与人之间绝大部分DNA序列是相同的，很多时候，比如针对不同人群的药物研发以及即将到来的个性化医疗时代，都不需要研究所有的30亿个碱基对，而只需要关注存在差异的部分——包括序列碱基的插入（insert，INS）、缺失（delete，DEL）、重复（repeat，REP）和倒位（inversion，INV）等。这种策略可以很好地抓住关键点，同时也减少了不必要的劳动输出，在科学研究中是非常有意义的。因此，如何能够有效、全面地检测基因组上的差异，发现人与人之间的不同点，就成了所有这些研究的一个重要基础。
目前，已有很多研究人员和科研机构提出了许多种不同的方案来解决这样的问题，并发布了相关软件，如GATK、Platypus、Pindel、Dindel、Breakdancer等基因组序列变异检测软件。但是它们都存在一些局限性，就是只能检测到一些长度很短的序列插入和缺失（一般在50bp以下）以及一些发生了1000bp以上的序列缺失的变异区域，而无法很有效地解决50bp~1000bp的序列变异、大长度的序列插入、序列倒位、位置变换（Trans）以及其他更为复杂的序列结构性变异。之所以存在这种局限，是由于：（1）测序的序列读长（read）比较短，一个变异要能够被检测出来，最有效的证据就是测得的序列能够完整地将一个变异覆盖过去，但限于当前的测序技术，测序读长很难满足这一需求；（2）人类基因组本身存在着约50%的重复序列，这些重复序列会影响比对软件将短序列定位到基因组上的准确性；（3）这些软件本身的算法存在局限性，它们主要是通过统计学的方法来推断这些序列变异的情况，在推断过程中要考虑测序文库的片段长度，如Illumina的测序仪，测序文库构建的序列片段长度一般是300bp~500bp，软件在发现序列比对的结果与测序文库的序列长度存在差异时才会判断可能存在序列变异，但统计学的差异性都是有一定适应范围的，这就导致很多变异会被漏掉，从而形成变异检测的“暗区”，比如长度在50bp-1000bp之间的变异是最容易在这类测序文库中被漏掉的，这些变异并非不存在了，而是检测的精度在这一长度范围上变得很差，因此它也就成了一个序列变异检测的“暗区”。而对于那些大长度的序列插入，由于处在变异中的序列根本就没能比对上参考基因组，这些软件就更加无能为力了。
因此，开发一种全新且准确性高，并能全面获得基因组上所有类型变异的方法和流程成为现今后基因组时代一个亟待解决的问题。
发明内容
本发明提供一种基于基因组组装的变异检测方法和装置，有效解决变异检测“暗区”、长序列插入和复杂结构性变异的检测难题。
根据本发明的第一方面，本发明提供一种基于基因组组装的变异检测方法，包括：获取来源于梯度测序文库的测序读段序列，该梯度测序文库为包括至少3个梯度的测序文库的集合，上述测序读段序列为去除样本接头序列并进行样本区分的读段序列；对上述测序读段序列进行过滤，除去测序质量低于第一预设值的碱基个数超过整条序列碱基个数的第二预设值的序列，以及序列中测序结果不确定的碱基个数超过整条序列碱基个数的第三预设值的序列；按照上述梯度测序文库从小到大的方式，逐步将各个文库中已经过滤的读段序列加进拼接过程拼接成平均长度达到第四预设值的长序列；将拼接得到的长序列比对到参考基因组上；对序列比对结果进行变异检测，获取序列变异。
根据本发明的第二方面，本发明提供一种基于基因组组装的变异检测装置，包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括可执行的程序；处理器，与数据输入单元、数据输出单元及存储单元数据连接，用于执行可执行的程序，程序的执行包括完成上述基于基因组组装的变异检测方法。
根据本发明的方法，通过将短小的读段序列拼接成长序列完整地将一个变异覆盖过去，因此有效解决变异检测“暗区”、长序列插入和复杂结构性变异的检测难题，高效、快速地实现序列变异的全面检测，为解码个人基因组和实现个性化医疗打下重要基础，并填补了现有生物信息学方法在完整检测基因组序列变异方面的不足。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施方式的描述中将变得明显和容易理解，其中：
图1为依据本发明的一种实施方式的基于基因组组装的变异检测方法流程图；
图2为依据本发明的一个实施例的梯度测序文库拼接序列的长度变化图；
图3为依据本发明的一个实施例检测到的基因组序列变异的数量分布图；
图4为依据本发明的一个实施例检测到的基因组序列变异的长度分布图。
具体实施方式
依据本发明的一种实施方式，提供一种基于基因组组装的变异检测方法，参考图1，包括如下步骤：
S1：获取来源于梯度测序文库的测序读段序列。
梯度测序文库通常是将来自目标个体的基因组样本经过打断成片段，并根据所选用的测序方法进行相应的文库（library）制备获得的，可选用的测序方法根据来自的测序平台包括但不限于CG（Complete Genomics）、Illumina/ Solexa、ABI/ SOLiD和Roche 454，依据所选测序平台进行单端或双端测序文库的制备。本发明一个实施例中，以Illumina/Solexa的新一代测序仪Genome Analyzer作为高通量测序方法实现梯度测序文库，该测序技术是一种基于边合成边测序（SBS，Sequencing-By-Synthesis）的新型测序方法，通过利用单分子阵列实现在小型芯片（Flow Cell）上进行桥式PCR反应。新的可阻断技术可实现每次只合并一个碱基，不需要标记荧光基团，再利用相应的激光激发荧光基团捕获激发光，从而读取碱基信息。应当理解，测序方法并不构成对本发明的限制。
本发明中梯度测序文库是指文库的平均片段长度呈梯度化分布的多个文库的集合，例如至少3个梯度的测序文库的集合，比如3至7个平均片段长度逐渐递增的文库，优选7个平均片段长度逐渐递增的文库，例如片段平均长度为160-200bp、400-600bp、700-1000bp、2kbp、5kbp、10kbp和20kbp共7个梯度的测序文库。一般来讲，至少应该有3个梯度的测序文库，比如160-200bp、400-600bp和2kbp这3个梯度，最好上述7个梯度的测序文库均有。梯度测序文库的数量和平均片段长度决定了拼接出的长片段的长度和质量情况，如果只有160-200bp、400-600bp和2kbp这3个梯度的小片段文库，一般N50值达到100kbp左右；如果上述7个梯度的测序文库均有，N50值可以达到10Mbp。其中N50定义为：用于衡量序列拼接结果的完整程度，具体是所有拼接出来的各个长序列片段，按照其长度从大到小排序，然后从最长的序列开始往下累加，当累加到某条序列时，累加出来的长度就达到或超过所有这些序列的总长的50%时，这条临界序列的长度就称为N50。上述7个梯度的测序文库，可以比较有效地保证最后序列拼接的结果。需要指出的是，本发明中N50与平均长度是不同的概念。平均长度定义为拼接出来的序列总长度除以序列的条数。
梯度测序文库构建过程中的DNA样本的提取和纯化以及接头连接等操作，可视具体测序平台而定，目前已经是成熟的技术。梯度测序文库的片段大小控制也有专门的仪器可以实现，例如Covaris的超声波打断仪，通过对具体工作参数的设置即可得到特定长度的打断片段。
梯度测序文库的测序深度，一般需要满足小片段（<2kbp）文库测得深一些，大片段（≥2kbp）文库测得浅一些，例如160-200bp的测序文库的测序深度可以是20×~30×，例如20×、22×、25×、27×或30×，优选20×或30×；400-600bp和700-1000bp的测序文库的测序深度在20×以下即可，例如18×、15×、12×、10×或8×，优选10×~20×；大片段测序文库的测序深度一般控制在10×以下，例如9×、8×、7×、6×、5×或3×，优选5×~10×。而且，总测序深度最好不要超过80×，例如70×~80×，不然测序错误率虽然低，但如果测序深度过高，其影响反而也会被放大。
依据测序平台的情况，下机的测序读段序列已经去除了测序接头序列，并且已经依据样本接头序列进行了样本区分，需要进行去除样本接头序列并进行样本区分，是由于样本接头序列是在高通量测序中引入的用来区别不同样本的序列，由于第二代高通量测序平台的测序通量巨大，一般都是将不同样本混合在一起上机测序。测序完成以后，需要依据样本接头序列将不同来源的样本区分开来，并需要将样本接头序列去除，以利于组装。来源于同一样本的读段序列具有相同的样本接头序列，而不同来源的读段序列的样本接头序列不同，这些不同的样本接头序列组成样本接头序列库。
本发明提供一种去除样本接头序列并进行样本区分的方法，包括：去除样本接头序列中测序质量低于第五预设值的碱基个数大于第六预设值的读段序列；顺序任意地a）将样本接头序列与样本接头序列库中的序列进行完全匹配，b）假设样本接头序列降解第七预设值个碱基再与样本接头序列库中序列对应部分进行完全匹配，c）允许样本接头序列有第八预设值个碱基的插入进行完全匹配，以及d）允许样本接头序列有第九预设值个碱基的缺失进行完全匹配；舍弃a）~d）中均无匹配结果、a）~d）中其中一个同时匹配到两个结果以及仅有且同时c）和d）匹配出结果的整条序列；将匹配到样本接头序列库中同一序列的读段序列作为同一样本来源的序列而实现样本区分，并去除读段序列中的样本接头序列。
去除样本接头序列中测序质量低于第五预设值的碱基个数大于第六预设值的读段序列这一步骤的实质是以样本接头序列的状况作为测序的质控指标。其中，此处的测序质量是衡量测序平台的测序结果的质量状况的指标，依据具体测序平台，如本发明一个实施例中的Illumina/ Solexa测序平台，可定义为：q=-10log 10 p error ，其中，q表示测序质量，p error 表示测序错误率，由测序平台根据测序过程中多个不同因素给出，p error 越小测序质量就越高。
依据具体测序技术及测序环境，对第五预设值和第六预设值进行设置。本发明经研究确定在下列数值范围内既能去除不合格的序列，同时又能保证尽可能多的合格序列得以保留，提高测序读段序列的利用率。具体地，第五预设值可以设置为4~6，第五预设值越大，意味着去除标准越严格，相应去除的读段序列越多，在本发明一个具体实施例中第五预设值设置为5；第六预设值可以设置为2~4，第六预设值越大，意味着去除标准越宽松，相应去除的读段序列越少，在本发明一个具体实施例中第六预设值设置为3。
上述完全匹配操作包括4个步骤，即a）、b）、c）和d），可以按任意顺序进行。一般步骤a）最先进行，因为这一步骤是最严格的匹配，称得上是真正意义上的“完全匹配”，而步骤b）、c）和d）一般在步骤a）之后进行。本发明的一个具体实施例中，优选地按a）、b）、c）和d）的顺序进行匹配操作。
之所以进行步骤b）的匹配，是因为考虑到一系列实验过程中，样本接头序列可能出现降解情况，在降解不严重的情况下，如果只是按照步骤a）最严格的匹配，可能会使得一部分合格的读段序列被排出，降低了读段序列的利用率。依据具体测序技术及测序环境，对步骤b）中的第七预设值进行设置。第七预设值越大，意味着匹配标准越宽松。本发明经研究确定第七预设值在1~2范围内，即假设样本接头序列降解1或2bp，再与样本接头序列库中序列对应部分进行完全匹配，既能去除不合格的序列，同时又能保证尽可能多的合格序列得以保留，提高测序读段序列的利用率。在本发明一个具体实施例中第七预设值设置为1。
之所以进行步骤c）的匹配，是因为考虑到一系列实验过程中，样本接头序列发生碱基插入，在碱基插入不严重的情况下，如果只是按照步骤a）最严格的匹配，可能会使得一部分合格的读段序列被排出，降低了读段序列的利用率。依据具体测序技术及测序环境，对步骤c）中的第八预设值进行设置。第八预设值越大，意味着匹配标准越宽松。本发明经研究确定第八预设值在1~2范围内，即允许样本接头序列有1或2bp碱基的插入进行完全匹配，既能去除不合格的序列，同时又能保证尽可能多的合格序列得以保留，提高测序读段序列的利用率。在本发明一个具体实施例中第八预设值设置为1，也就是允许样本接头序列仅有1个碱基的插入，从样本接头序列起始端进行完全匹配操作，当出现某碱基无法匹配时认为该碱基为插入碱基，跳过此碱基后继续严格的完全匹配操作。
之所以进行步骤d）的匹配，是因为考虑到一系列实验过程中，样本接头序列发生碱基缺失，在碱基缺失不严重的情况下，如果只是按照步骤a）最严格的匹配，可能会使得一部分合格的读段序列被排出，降低了读段序列的利用率。依据具体测序技术及测序环境，对步骤d）中的第九预设值进行设置。第九预设值越大，意味着匹配标准越宽松。本发明经研究确定第九预设值在1~2范围内，即允许样本接头序列有1或2bp碱基的缺失进行完全匹配，既能去除不合格的序列，同时又能保证尽可能多的合格序列得以保留，提高测序读段序列的利用率。在本发明一个具体实施例中第九预设值设置为1，也就是允许样本接头序列仅有1个碱基的缺失进行完全匹配操作。
完成上述a）、b）、c）和d）四步匹配操作后，可以优选地按照a）>b）>c）>d）的优先级顺序确定最终的样本接头序列的比对结果，而对于上述四步匹配操作中均无匹配结果，其中一个匹配步骤同时匹配到两个结果，或仅有步骤c）和d）的匹配结果，认为是无效信息，将相应的整条读段序列去除。将匹配到样本接头序列库中同一序列的读段序列作为同一样本来源的序列而实现样本区分，并去除读段序列中的样本接头序列。
S2：对测序读段序列进行过滤。
由于下机的测序读段序列可能包含一些不合格的序列，这些序列的存在会影响读段序列的组装质量，因此需要先进行过滤处理。
过滤处理主要是去除以下两类不合格的序列：测序质量低于第一预设值的碱基个数超过整条序列碱基个数的第二预设值的序列，以及序列中测序结果不确定的碱基个数（如测序结果中的N的个数）超过整条序列碱基个数的第三预设值的序列。
测序质量是衡量测序平台的测序结果的质量状况的指标，依据具体测序平台，如本发明一个实施例中的Illumina/ Solexa测序平台，可定义为：q=-10log 10 p error ，其中，q表示测序质量，p error 表示测序错误率，由测序平台根据测序过程中多个不同因素给出，p error 越小测序质量就越高。
依据具体测序技术及测序环境，对第一预设值、第二预设值和第三预设值进行设置。本发明经研究确定在下列数值范围内既能去除不合格的序列，同时又能保证尽可能多的合格序列得以保留，提高测序读段序列的利用率。具体地，第一预设值可以设置为4~6，第一预设值越大，意味着过滤标准越严格，相应去除的读段序列越多，在本发明一个具体实施例中第一预设值设置为5；第二预设值可以设置为40%~60%，第二预设值越大，意味着过滤标准越宽松，相应去除的读段序列越少，在本发明一个具体实施例中第二预设值设置为50%；第三预设值可以设置为5%~15%，第三预设值越大，意味着过滤标准越宽松，相应去除的读段序列越少，在本发明一个具体实施例中第三预设值设置为10%。
S3：将已经过滤的读段序列拼接成平均长度达到第四预设值的长序列。
本发明的变异检测方法，在完成以上质控（即步骤S2）之后，并不是与以往的变异检测软件一样，直接将这些测序的读段序列比对到参考基因组上，而是先组装成长片段，这是本发明最关键的构思之一。将短的读段序列拼接成长序列，最主要的目的是为了使尽可能多的变异序列能被完整地包含在拼接出来的长序列中，从而解决短序列在这一方面的不足，如背景技术中介绍的变异检测“暗区”、长序列插入和复杂结构性变异的检测难题。本发明一个优选实施例中具体是按照梯度测序文库从小到大逐步将各个文库加进拼接过程而实现长片段拼接。本发明一个实施例中，采用160-200bp、400-600bp、700-1000bp、2kbp、5kbp、10kbp和20kbp共7个梯度的测序文库进行拼接，使得拼接后的长序列平均长度在1百万（Mbp）个碱基对以上，是当前最先进的二代测序技术所测到的序列读长的4000倍！基本上已可以保证人类基因组上绝大部分的序列变异都能被这样长的序列所包括。用于读段序列拼接的技术目前有多种，本发明一个实施例中采用SOAPdenovo2（网址是http://soap.genomics.org.cn/soapdenovo.html；参考Ruibang Luo 等人, “SOAPdenovo2: AnEmpirically Improved Memory-Efficient Short-Read de Novo Assembler.,”GigaScience, 1 (2012), 18 <doi:10.1186/2047-217X-1-18>）进行读段序列拼接。
本发明中的第四预设值与梯度测序文库的数量以及大小有关，梯度测序文库的数量越多、大小越大，第四预设值能够达到的数值越大，例如本发明一个实施例中7个梯度测序文库进行拼接的结果是得到长序列平均长度为1百万个碱基对。
S4：将拼接得到的长序列比对到参考基因组上。
在本发明的变异检测方法中，经过序列拼接这一步骤之后，将拼接得到的长序列比对到参考基因组上，参考基因组序列可取于公共数据库（例如美国国家生物技术信息中心（NCBI，national center for biotechnology information）提供的HG19）。本发明中，使用拼接后的长序列还可以解决因参考基因组存在重复性序列而导致短序列出现比对错误的问题。由于短序列已变成很长的序列，这使得以前用于进行短序列比对的软件（比如BWA）不再适用于本发明的情景，而只能使用长序列比对软件，如LAST等。本发明一个实施例中使用LAST软件（网址是http://last.cbrc.jp/）将这些拼接后的长序列比对到参考基因组HG19上。
S5：对序列比对结果进行变异检测，获取序列变异。
在本发明的变异检测方法中，确定了长序列比对的结果之后，采用变异检测技术识别变异区域。主要检测原理是：扫描整个序列比对结果，获得每一段比对序列中所有与参考基因组存在差异的序列及其位置，检测出每一比对结果中的变异；对分段比对的结果按照拼接序列的物理位置顺序排序，挑选出顺序不正常的序列段和位置，得到结构性变异区域；对每个检测出来的变异进行局部重比对，重新确定其精确的序列位置。本发明中的该步骤可以采用目前已有的变异检测软件实现，本发明一个实施例中使用AsmVar（网址是http://www.stbioinf.com/AsmVar/），检测出所有可能的序列变异。
由于短序列已经拼接成了长序列，使用变异检测软件除了能够检测到小序列变异，还特别能够有效地检测出50bp~1000bp的序列变异、大长度的序列插入、序列倒位、位置变换以及其他更为复杂的序列结构性变异情况，并可以利用统计学的方法判定出这些变异在所研究的样本中的基因型。
需要说明的是，上述步骤S3~S5中使用的软件SOAPdenovo2、LAST和AsmVar仅仅是一个具体实施例中使用的软件，并不构成对本发明的限制，因为目前有其它替代技术可以实现同样的目的。
本领域普通技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。
依据本发明的另一方面还提供一种基于基因组组装的变异检测装置，包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括可执行的程序；处理器，与上述数据输入单元、数据输出单元及存储单元数据连接，用于执行上述可执行的程序，该程序的执行包括完成上述实施方式中各种方法的全部或部分步骤。
以下结合具体目标个体对依据本发明的具体检测方法的运行结果进行详细的描述。下述检测过程所使用的具体参数设置为：
1、实施例样本：10个家系共30个正常人的血液样本；
2、建立平均长度为160-200bp、400-600bp、700-1000bp、2kbp、5kbp、10kbp和20kbp共7个梯度测序文库，使用Illumina/ Solexa的新一代测序仪Genome Analyzer进行基于边合成边测序的新型测序方法的测序，其中各个文库的测序深度如下：160-200bp文库测30×，400-600bp文库测10×，700-1000bp文库测10×，2kbp~20kb这四个文库每个测5×，最终每个样本平均测70×。
3、设置第一预设值为5，第二预设值为50%，第三预设值为10%，第五预设值为5，第六预设值为3，第七预设值、第八预设值和第九预设值均为1，进行过滤、去除样本接头序列和样本区分，其中匹配操作按照上述a）、b）、c）和d）的顺序进行；
4、使用SOAPdenovo2序列拼接软件，将过滤之后的短测序读段序列拼接为精确的长序列，并将所有长度小于100bp的序列过滤掉；
5、采用长序列比对软件LAST将进行过拼接之后的长序列比对到参考基因组上，精确定位这些序列在基因组上的位置；
6、使用AsmVar软件对整个序列比对结果进行扫描，获得所有可能的变异列表，除此之外，AsmVar还对这些变异的基因型进行判定，并利用高斯混合模型计算每个变异的质量值，该值是变异位点可靠性的重要评断标准。
本实施例的主要结果如下。需要注意的是，这里为了简明起见给出的图示只是实施例的一部分结果，并非全部。
图2示出了在序列拼接过程中，随着梯度测序文库的逐步加入，拼接序列的N50随之不断增大。序列拼接过程中，严格按照梯度测序文库从小到大逐步将各个文库加进拼接的过程中。由图2所示的结果可以看到，N50有不断累积增长的趋势和过程，直到最后一个文库20kbp加进去之后，拼接出来的序列N50达到了10Mbp左右。图2中也可以看出，如果不是梯度文库，比如只有小片段文库（<2kbp），那么N50最大也达不到100kbp。另一方面，本发明中所用的这7个梯度文库，其间隔步长是比较合适的，不至于出现前一步与后一步差距过大，导致前后衔接不起来，反而不能起到明显增长N50的目的。图2中，附图符号“1006-01”等表示30个正常人的血液样本中其中一个的编号。
图3示出的是本实施例所检测到的不同序列变异的数量分布图。由图3所示的结果可以看出，使用本发明的方法检测到的序列变异类型更加全面，除了较容易检测到的序列缺失（即序列删除）和序列插入之外，还检测到序列倒位（INV）、位置变换（Trans）和序列片段替换（BSubstitution）——包括等长替换和不等长替换。图3中，附图符号“1006-01”等表示30个正常人的血液样本中其中一个的编号。
图4示出的是本实施例检测到的基因组序列变异的长度分布图，图形从中间左右划分，左边是序列缺失（即序列删除）的长度分布，右边是序列插入的长度分布。需要注意的是，为了作图的方便，本实施例将其他类型的变异，包括序列倒位、位置变换以及更为复杂的序列结构性变异都区分为序列缺失或者序列插入，具体方法是，如果变异的序列长度小于比对位置断点之间的长度则区分为序列缺失，反之则区分为序列插入。分布图中的已知转座子数据库、已知短序列插入缺失数据库、单碱基突变数据库、基因组结构性变异数据库和千人数据库表示的柱状部分是本方法和当前短序列结合统计推断的方法均能检测到的变异集合；而空白柱状部分（新发现变异）是本方法检测到的全新的变异，从中可以非常明显地看到，本发明的方法所能检测到的序列变异比起当前短序列结合统计推断的方法，长度范围更全面，数量也更多，有效地解决了当前序列变异检测方面的不足。
以上所述仅为本发明的较佳实施例，应当理解，这些实施例仅用以解释本发明，并不用于限定本发明。对于本领域的一般技术人员，依据本发明的思想，可以对上述具体实施方式进行变化。
Claims (14)
1.一种基于基因组组装的变异检测方法，其特征在于，包括：
获取来源于梯度测序文库的测序读段序列，所述梯度测序文库为包括至少3个梯度的测序文库的集合，所述测序读段序列为去除样本接头序列并进行样本区分的读段序列；
对所述测序读段序列进行过滤，除去测序质量低于第一预设值的碱基个数超过整条序列碱基个数的第二预设值的序列，以及序列中测序结果不确定的碱基个数超过整条序列碱基个数的第三预设值的序列；
按照所述梯度测序文库从小到大的方式，逐步将各个文库中已经过滤的读段序列加进拼接过程拼接成平均长度达到第四预设值的长序列；
将拼接得到的长序列比对到参考基因组上；
对序列比对结果进行变异检测，获取序列变异；
所述对序列比对结果进行变异检测包括：
扫描整个序列比对结果，获得每一段比对序列中所有与参考基因组存在差异的序列及其位置，检测出每一比对结果中的变异；
对分段比对的结果按照拼接序列的物理位置顺序排序，挑选出顺序不正常的序列段和位置，得到结构性变异区域；
对每个检测出来的变异进行局部重比对，重新确定其精确的序列位置。
2.根据权利要求1所述的方法，其特征在于，所述梯度测序文库为包括7个梯度的测序文库的集合。
3.根据权利要求1所述的方法，其特征在于，所述梯度测序文库为片段平均长度为160-200bp、400-600bp、700-1000bp、2kbp、5kbp、10kbp和20kbp共7个梯度的测序文库的集合。
4.根据权利要求3所述的方法，其特征在于，160-200bp文库的测序深度为20×～30×；400-600bp和700-1000bp文库的测序深度分别为20×以下；2kbp、5kbp、10kbp和20kbp文库的测序深度分别为10×以下；总测序深度为80×以下。
5.根据权利要求4所述的方法，其特征在于，400-600bp和700-1000bp文库的测序深度分别为10×～20×；2kbp、5kbp、10kbp和20kbp文库的测序深度分别为5×～10×；总测序深度为70×～80×。
6.根据权利要求1-3任一项所述的方法，其特征在于，所述测序质量通过如下公式计算：q＝-10log10perror，其中，q表示测序质量，perror表示测序错误率；所述第一预设值为4～6，所述第二预设值为40％～60％，所述第三预设值为5％～15％。
7.根据权利要求6所述的方法，其特征在于，所述第一预设值为5，所述第二预设值为50％，所述第三预设值为10％。
8.根据权利要求1-3任一项所述的方法，其特征在于，所述第四预设值为50万以上。
9.根据权利要求8所述的方法，其特征在于，所述第四预设值为100万以上。
10.根据权利要求1-3任一项所述的方法，其特征在于，所述去除样本接头序列并进行样本区分具体包括：
去除样本接头序列中测序质量低于第五预设值的碱基个数大于第六预设值的读段序列；
顺序任意地a)将样本接头序列与样本接头序列库中的序列进行完全匹配，b)假设样本接头序列降解第七预设值个碱基再与样本接头序列库中序列对应部分进行完全匹配，c)允许样本接头序列有第八预设值个碱基的插入进行完全匹配，以及d)允许样本接头序列有第九预设值个碱基的缺失进行完全匹配；
舍弃a)～d)中均无匹配结果、a)～d)中其中一个同时匹配到两个结果以及仅有且同时c)和d)匹配出结果的整条序列；
将匹配到样本接头序列库中同一序列的读段序列作为同一样本来源的序列而实现样本区分，并去除读段序列中的样本接头序列。
11.根据权利要求10所述的方法，其特征在于，按顺序a)、b)、c)和d)进行所述完全匹配。
12.根据权利要求11所述的方法，其特征在于，所述第五预设值为4～6，所述第六预设值为2～4，所述第七预设值、第八预设值和第九预设值均为1～2。
13.根据权利要求12所述的方法，其特征在于，所述第五预设值为5，所述第六预设值为3，所述第七预设值、第八预设值和第九预设值均为1。
14.一种基于基因组组装的变异检测装置，其特征在于，包括：
数据输入单元，用于输入数据；
数据输出单元，用于输出数据；
存储单元，用于存储数据，其中包括可执行的程序；
处理器，与所述数据输入单元、数据输出单元及存储单元数据连接，用于执行所述可执行的程序，所述程序的执行包括完成如权利要求1-13任意一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510043893.2A CN105989246B (zh) | 2015-01-28 | 2015-01-28 | 一种基于基因组组装的变异检测方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510043893.2A CN105989246B (zh) | 2015-01-28 | 2015-01-28 | 一种基于基因组组装的变异检测方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105989246A CN105989246A (zh) | 2016-10-05 |
| CN105989246B true CN105989246B (zh) | 2018-10-26 |
Family
ID=57034232
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510043893.2A Active CN105989246B (zh) | 2015-01-28 | 2015-01-28 | 一种基于基因组组装的变异检测方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105989246B (zh) |
Families Citing this family (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| SG10202104266VA (en) * | 2016-11-16 | 2021-05-28 | Illumina Inc | Methods of sequencing data read realignment |
| CN106611106B (zh) * | 2016-12-06 | 2019-05-03 | 北京荣之联科技股份有限公司 | 基因变异检测方法及装置 |
| CN106599616B (zh) * | 2017-01-03 | 2019-05-31 | 上海派森诺医学检验所有限公司 | 基于duplex-seq的超低频突变位点检测分析方法 |
| EP3571613A1 (en) * | 2017-01-17 | 2019-11-27 | Illumina, Inc. | Oncogenic splice variant determination |
| CN110462063B (zh) * | 2017-05-23 | 2023-06-23 | 深圳华大生命科学研究院 | 一种基于测序数据的变异检测方法、装置和存储介质 |
| CN107451428B (zh) * | 2017-08-02 | 2020-05-22 | 广东国盛医学科技有限公司 | 下一代测序中末端短串联序列的优化处理方法 |
| CN110021342B (zh) * | 2017-08-21 | 2020-12-15 | 北京哲源科技有限责任公司 | 用于加速变异位点的识别的方法及系统 |
| CN107784199A (zh) * | 2017-10-18 | 2018-03-09 | 中国科学院昆明植物研究所 | 一种基于总dna测序结果的细胞器基因组筛选方法 |
| CN107895104B (zh) * | 2017-11-13 | 2020-07-07 | 深圳华大基因科技服务有限公司 | 评估和校验三代测序的序列组装结果的方法与装置 |
| CN108171011B (zh) * | 2017-12-08 | 2020-09-29 | 志诺维思（北京）基因科技有限公司 | 一种dna复杂结构变异探测方法 |
| CN108073791B (zh) * | 2017-12-12 | 2019-02-05 | 元码基因科技(苏州)有限公司 | 基于二代测序数据检测目标基因结构变异的方法 |
| CN108830044B (zh) * | 2018-06-05 | 2020-06-26 | 序康医疗科技（苏州）有限公司 | 用于检测癌症样本基因融合的检测方法和装置 |
| CN109994155B (zh) * | 2019-03-29 | 2021-08-20 | 北京市商汤科技开发有限公司 | 一种基因变异识别方法、装置和存储介质 |
| CN110349629B (zh) * | 2019-06-20 | 2021-08-06 | 湖南赛哲医学检验所有限公司 | 一种利用宏基因组或宏转录组检测微生物的分析方法 |
| CN111243663B (zh) * | 2020-02-26 | 2022-06-07 | 西安交通大学 | 一种基于模式增长算法的基因变异检测方法 |
| CN113470742A (zh) * | 2020-03-31 | 2021-10-01 | 阿里巴巴集团控股有限公司 | 数据处理方法、装置、存储介质及计算机设备 |
| WO2021253346A1 (zh) * | 2020-06-18 | 2021-12-23 | 李雨澄 | 数据传输计算方法，装置及存储介质 |
| CN115831223B (zh) * | 2023-02-20 | 2023-06-13 | 吉林工商学院 | 一种挖掘近源物种间染色体结构变异的分析方法及系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102867134A (zh) * | 2012-08-16 | 2013-01-09 | 盛司潼 | 一种对基因序列片段进行拼接的系统和方法 |
| CN103014137A (zh) * | 2011-09-22 | 2013-04-03 | 深圳华大基因科技有限公司 | 一种分析基因表达定量的方法 |
-
2015
- 2015-01-28 CN CN201510043893.2A patent/CN105989246B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103014137A (zh) * | 2011-09-22 | 2013-04-03 | 深圳华大基因科技有限公司 | 一种分析基因表达定量的方法 |
| CN102867134A (zh) * | 2012-08-16 | 2013-01-09 | 盛司潼 | 一种对基因序列片段进行拼接的系统和方法 |
Non-Patent Citations (3)
| Title |
|---|
| "AGE:defining breakpoints of genomic structural variants at single-nucleotide resolution,through optimal alignments with gap excision";Alexej Abyzov et al;《BIOINFORMATICS》;20110113;第27卷(第5期);595-603 * |
| "Novel variation and de novo mutation rates in population-wide de novo assembled Danish trios";Soren Besenbacher et al;《NATURE COMMUNICATIONS》;20150119;1-9 * |
| "Pseudo-Sanger sequencing: massively parallel production of long and near error-free reads using NGS technology";Jue Ruan et al;《BMC Genomics》;20131017;1-9 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105989246A (zh) | 2016-10-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105989246B (zh) | 一种基于基因组组装的变异检测方法和装置 | |
| US11371074B2 (en) | Method and system for determining copy number variation | |
| De Coster et al. | Towards population-scale long-read sequencing | |
| AU2020202153B2 (en) | Single-molecule sequencing of plasma DNA | |
| Costa et al. | Uncovering the complexity of transcriptomes with RNA-Seq | |
| CN108138227A (zh) | 使用具有独特分子索引(umi)的冗余读段在测序dna片段中抑制误差 | |
| CA3050247A1 (en) | Methods and systems for generation and error-correction of unique molecular index sets with heterogeneous molecular lengths | |
| CN110168648A (zh) | 序列变异识别的验证方法和系统 | |
| EP4095262A2 (en) | Quality evaluation method, quality evaluation apparatus, program, storage medium, and quality control sample | |
| KR20190019219A (ko) | 모체 혈장으로부터의 비침습적 산전 분자 핵형분석 | |
| CN103902852A (zh) | 基因表达的定量方法及装置 | |
| CN108728522A (zh) | 药物基因组检测方法 | |
| CN113278706A (zh) | 一种用于区分体细胞突变和种系突变的方法 | |
| Goswami et al. | RNA-Seq for revealing the function of the transcriptome | |
| Ahmadian et al. | Massively parallel sequencing platforms using lab on a chip technologies | |
| CN104428423A (zh) | 确定外源基因在人类基因组中整合方式的方法和系统 | |
| WO2017009718A1 (en) | Automatic processing selection based on tagged genomic sequences | |
| CN109321646A (zh) | 基于ngs读段与参考序列比对的虚拟pcr方法 | |
| WO2022262569A1 (zh) | 一种用于区分体细胞突变和种系突变的方法 | |
| US20230144221A1 (en) | Methods and systems for detecting alternative splicing in sequencing data | |
| WO2016141516A1 (zh) | 获取子代特异性序列、检测子代新突变的方法和装置 | |
| Lin | Developing A Nanopore Sequencing Data Processing Pipeline for Structural Variation Identification | |
| Muzzey | Understanding the Basics of NGS in the Context of NIPT | |
| Hsu | Molecular Methods for Diagnosis of Genetic Diseases Involving the Immune System | |
| CN116403641A (zh) | 碱基测序错误的排除方法、低频率突变的鉴定方法及相关装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| REG | Reference to a national code |
Ref country code: HK Ref legal event code: DE Ref document number: 1228027 Country of ref document: HK |
|
| CB02 | Change of applicant information |
Address after: 518083 comprehensive building, Beishan Industrial Zone, Yantian District, Shenzhen, Guangdong Applicant after: BGI SHENZHEN Address before: 518083 comprehensive building, Beishan Industrial Zone, Yantian District, Shenzhen, Guangdong Applicant before: BGI SHENZHEN |
|
| CB02 | Change of applicant information | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20180530 Address after: 518083 the comprehensive building of Beishan industrial zone and 11 2 buildings in Yantian District, Shenzhen, Guangdong. Applicant after: MGI TECH Co.,Ltd. Address before: 518083 comprehensive building, Beishan Industrial Zone, Yantian District, Shenzhen, Guangdong Applicant before: BGI SHENZHEN |
|
| TA01 | Transfer of patent application right | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CP01 | Change in the name or title of a patent holder | ||
| CP01 | Change in the name or title of a patent holder |
Address after: 518083 comprehensive building of Beishan industrial zone and 11 Building 2, Yantian District, Guangdong, Shenzhen Patentee after: Shenzhen Huada Zhizao Technology Co.,Ltd. Address before: 518083 comprehensive building of Beishan industrial zone and 11 Building 2, Yantian District, Guangdong, Shenzhen Patentee before: MGI TECH Co.,Ltd. |
|
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20210607 Address after: 518000 2 / F, building 11, Beishan Industrial Zone, Yantian street, Yantian District, Shenzhen City, Guangdong Province Patentee after: Shenzhen Huada Zhizao Software Technology Co.,Ltd. Address before: 518083 the comprehensive building of Beishan industrial zone and 11 2 buildings in Yantian District, Shenzhen, Guangdong. Patentee before: Shenzhen Huada Zhizao Technology Co.,Ltd. |