基因组断裂点识别方法及应用与流程

本发明涉及生物信息学分析
技术领域：
，特别是涉及一种基因组断裂点识别方法及应用。
背景技术：
：随着基因检测技术的成熟和普及，基因检测被越来越多的利用到致病基因研究、精准治疗中。通常来说，可通过将ngs测序数据比对到参考基因组上，找出个体与参考基因组的差异，这些差异除了常规的差异诸如点突变(snv)和微小插入缺失(indel)，还有一些比较复杂的结构差异，如染色体断裂、重排、倒位、易位、大片段插入缺失等，从而通过差异进行进一步的研究或诊断。对于结构性变异，与常规变异一样，结构性变异也分多态性和致病性，即每个人的基因组上都或多或少有一些多态性的结构变异，如根据文献报道，等位基因频率超过50％的结构变异就有1000个之多，而病人携带的致病性变异往往只有1个或2个，因此快速评估一个结构变异的致病性(排除非致病位点)为后续分析的前提。而由于基因组的复杂性，目前的结构变异识别软件的结果假阳性率较高，且由于多态性变异的存在，大量位点无法判断其致病性。并且，常规的结构变异识别软件主要用于全基因组测序，对于全外显子测序数据，由于测序数据在参考基因组上不连续，一般结构变异至少会有2个断裂点，如果发生断裂的位置不是都在外显子上，则通过外显子测序获得的数据是无法识别出具体变异的。但现实中的情况确是：外显子占全基因组比例较低，所有断裂点都发生在外显子的概率很低，所以通过常规方法处理外显子测序数据来识别结构变异很困难。且外显子捕获过程中产生的一些异常数据及偏向性会对目前软件的结果产生较大影响。技术实现要素：基于此，有必要针对上述问题，提供一种基因组断裂点识别方法，采用该方法，可迅速识别基因组中结构变异位点，且不受是否为全基因组测序或全外显子组测序等不同类型的限制，可适用于各种测序方案，如全基因组测序、全外显子测序、转录组测序等等。一种基因组断裂点识别方法，包括以下步骤：数据比对：获取样本测序下机数据，将样本reads与参考基因组进行对比，当一个read无法完整比对到参考基因组时，则分别按照read左、右两端最优匹配输出，则该read被分开成为左、右两端的分界即为read分割点；构建边缘坐标集合：汇总所有reads比对到参考基因组上的起始坐标和终止坐标，所述read分割点在参考基因组比对坐标轴上的位置即为read边缘坐标；识别断裂点：筛选reads边缘坐标集中的位点，当该位点覆盖超过阈值，且该处产生了分割点的reads占总reads比例大于预设比例时，则判断为断裂点；判定断裂点：将上述得到的断裂点在预设数据库中进行出现频次的查询，当出现频次小于预设频次时，则判断为高风险致病性结构变异。上述基因组断裂点识别方法，不直接识别突变类型，而是识别断裂点，并通过汇总足够样本的断裂点结果建立数据库，从而得到某一断裂点在人群中的频率。当分析一个新样本时，将识别到的断裂点到该数据库中检索，可以排除大量假阳性或者人群多态性的断裂点。这是因为假阳性位点经常会在多个样本中出现，其表现形式与多态性位点类似，所以两者都可以通过数据库频率排除。可以理解的，如样本中reads为正常read，其有左右各一个边缘坐标，而存在断裂点的read无法完整比对到参考基因组，将分别按照read左、右两端最优匹配分割后输出，被分割的read则会产生4个边缘坐标，即分割点会产生2个额外的边缘坐标。在其中一个实施例中，所述数据比对步骤中，所述样本测序下机数据选自：全基因组测序数据、全外显子测序数据、转录组测序数据中的任一种。常规技术中心，外显子上的结构性异常无法通过全外检测出，这是由于：一般结构变异至少会有2个断裂点，如果发生断裂的位置不是都在外显子上，则通过外显子测序获得的数据是无法识别出具体变异的。但常见情况确是：外显子占全基因组比例较低，所有断裂点都发生在外显子的概率很低，所以通过常规方法处理外显子测序数据来识别结构变异很困难。在其中一个实施例中，所述数据比对步骤中，所述参考基因组选自人类基因组序列。可以理解的，上述人类基因组序列为人类基因组计划完成的人类基因组序列，可由ucsc和ncbi获取。在其中一个实施例中，所述识别断裂点步骤中，当测序类型是全基因组测序时，所述预设比例为25％；当测序类型是全外显子测序时，所述预设比例为20％。人类有两条染色体，理论情况下当断裂点为纯合时，产生了分割点的reads(split-reads)占比为100％，当断裂点为杂合时，split-reads占比为50％。但由于断裂点处reads捕获率往往会更低，这会导致split-reads占比下降。因此对于不需捕获的全基因组测序，split-reads占比阈值可以提高至25-30％，而对于全外显子测序，split-reads占比阈值在20％较为适宜。在其中一个实施例中，所述识别断裂点步骤中，所述阈值为10个reads。在实践工作中，发明人发现，当碰到测序深度很低的地方，例如只有10x，则只需2个reads就占比超过20％，如此随机性太大，会造成大量假阳性。且鉴于目前全外测序的平均深度都在100x左右，10个reads占比只有10％，也就是正常情况下都能达到，上述设定阈值可用来排除低覆盖区域。即将阈值设定在为10个reads，既可以避免大量假阳性，又可排除低覆盖区域。在其中一个实施例中，所述预设数据库通过以下方法建立：获取超过1000例样本数据，按照上述数据比对、构建边缘坐标集合、识别断裂点的步骤，得到样本中所有断裂点，并统计相同断裂点出现的频次。优选地，获取超过5000例样本建立数据库具有更优的效果。可以理解的，上述预设数据库需满足一定的样本量，根据研究表明，致病的结构变异发生率低于千分之一，因此数据库样本量需大于1000例才能保证能计算得到低于千分之一的频率。样本量达到5000例则有较好的容错率，假阳和假阴性率都能控制在较低水平。在其中一个实施例中，所述预设频次为样本数的0.1％。将频次设定在该阈值，既能够降低假阳性比例，又可识别出真正具有高风险的变异位点。本发明还公开了一种基因组断裂点识别系统，包括：存储模块：用于获取样本测序下机数据，并进行数据存储；分析模块：按照上述的基因组断裂点识别方法，对数据进行分析；输出模块：用于将分析得到的高风险致病性结构变异输出。上述识别系统，可配合检测系统使用，也可仅作为后端数据分析工具使用。在其中一个实施例中，所述输出模块还输出与该样本匹配的患者临床信息，以及所述高风险致病性结构变异所在基因的功能信息。将上述判断为高风险致病性结构变异，与临床信息匹配对比后，可利于解释该变异的临床意义，特别是当变异是在外显子上发生的断裂，是比较好解释其临床意义的。若一个基因发生了断裂，那这个基因的功能大概率会受影响。所以发生在已知功能基因上的断裂点只需要看这个基因功能即可。例如：当fga基因发生断裂，其中，fibrinogenalphachain，即纤维蛋白原的alpha链，其功能缺失会导致先天性纤维蛋白原缺乏症，症状为血液不凝固，出血不止。如该fga基因发生断裂，与病人的临床信息可匹配，基本可以确定fga的断裂就是致病原因。本发明还公开了上述的基因组断裂点识别方法在人类基因组多态性研究中的应用。可以理解的，在基因组结构性变异中，缺乏临床信息支持或意义不明确的断裂点繁不胜数，这些断裂点也许并非都致病，但也可能具有某些生物学功能，是人类基因多态性的组成部分之一，具有研究价值。因此，上述基因组断裂点识别方法在人类基因组多态性研究中，具有重要的应用价值。与现有技术相比，本发明具有以下有益效果：本发明的一种基因组断裂点识别方法，不直接识别突变类型，而是识别断裂点，并通过汇总足够样本的断裂点结果建立数据库，从而得到某一断裂点在人群中的频率。当分析一个新样本时，将识别到的断裂点到该数据库中检索，可以排除大量假阳性或者人群多态性的断裂点。并且，以全外显子测序结果为例，在分析得到的千余个断裂点中，通过本发明的识别方法排除后，所剩高风险致病性结构变异的位点一般不超过20个，通常为10个左右，极大的降低了后续分析的难度。本发明的一种基因组断裂点识别系统，利用上述识别方法，可对测序数据进行分析，快速而准确的获得高风险致病性结构变异信息。附图说明图1为实施例1中可完整比对的全外显子测序后坐标集合的integrativegenomicsviewer图；图2为实施例1中无法完整比对的全外显子测序后坐标集合的integrativegenomicsviewer图。具体实施方式为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的较佳实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的
技术领域：
的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。本文所使用的术语“和/或”包括一个或多个相关的所列项目的任意的和所有的组合。以下实施例中所用测序数据，均为按照常规ngs测序方案得到，测序深度为100-200x之间，中位值在120x左右。临床样本来源和基本情况：样本共12385个，均为全外显子捕获ngs测序数据(illumina平台)，捕获范围为外显子加上侧翼50bp(约覆盖60mbp碱基)。实施例1一种基因组断裂点识别方法，包括以下步骤：一、数据比对。获取样本测序下机数据，使用比对软件(如bwa)将样本测序reads比对到参考基因组，覆盖断裂点位置的reads会被分割成两部分比对，具体为：当一个read无法完整比对到参考基因组时，则分别按照read左、右两端最优匹配输出，则该read被分开成为左、右两端的分界即为read分割点。二、构建边缘坐标集合。汇总所有reads比对到参考基因组上的起始坐标和终止坐标。如样本数据来源于全基因组测序，则因为全基因组是随机打断测序，该边缘坐标大致呈均匀分布。如果是全外显子捕获测序，该边缘坐标大致呈以捕获探针为中心对称的分布。如图1所示，图1为全外显子测序后坐标集合的integrativegenomicsviewer图，其中，上半部分表示覆盖深度，下半部分为一个个reads，浅灰色和深灰色分别表示正向和反向比对的reads。该图中显示，reads可完整比对到参考基因组，即read无边缘断裂点出现在坐标集合中。而当染色体上出现断裂点时，覆盖到断裂点上的reads无法完全比对到参考基因组上，只能从断裂点处分割成成左、右两部分分开比对，这会导致reads比对在参考基因组上的边缘坐标集中在断裂点处，如图2所示。图2为全外显子测序后坐标集合的integrativegenomicsviewer图，其中，上半部分表示覆盖深度，下半部分为一个个reads，浅灰色和深灰色分别表示正向和反向比对的reads。该图中显示，reads不可完整比对到参考基因组，即reads比对在参考基因组上的边缘坐标集中在断裂点处，即所述read分割点在参考基因组比对坐标轴上以read边缘出现。具体到图2所体现的断裂点，由于常染色体有2个拷贝，从该图2中可以看出有一半的reads是正常的，另外一半的reads被切断后才比对上。表示该染色体有一个拷贝的这个位置是正常的，另一个拷贝在这个位置断裂了。三、识别断裂点。筛选reads边缘坐标集中的位点，在本实施例中，当该位点覆盖超过10reads，且该处产生了分割点的reads(split-reads)占总reads比例超过20％时，则判断为断裂点。可以理解的，因人有两条染色体，所以理论情况下当断裂点为纯合时，split-reads占比为100％，当断裂点为杂合时，split-reads占比为50％。但由于断裂点处reads捕获率往往会更低，这会导致split-reads占比下降。因此对于不需捕获的全基因组测序，split-reads占比阈值可以提高至25-30％。四、判定断裂点。将上述得到的断裂点在预设数据库中进行出现频次的查询，当出现频次小于预设频次时，则判断为高风险致病性结构变异。1、预设数据库的建立。以全外显子测序方案为例，本发明以上述三个步骤分析了12385个样本，根据这些样本断裂点的坐标和方向为唯一标识进行汇总，并统计每个断裂点出现的频次、split-reads深度百分位、split-reads占比百分位。chr1中的部分统计数据示例如下表。表1.数据库的统计数据possidecntsrd-percentilesrr-percentile908199left6446[10,13,16,23,31,40,98][0.2,0.2,0.21,0.23,0.26,0.28,0.45]908644right1[64,64,64,64,64,64,64][0.35,0.35,0.35,0.35,0.35,0.35,0.35]908734right1[59,59,59,59,59,59,59][0.5,0.5,0.5,0.5,0.5,0.5,0.5]908857left1[56,56,56,56,56,56,56][0.52,0.52,0.52,0.52,0.52,0.52,0.52]908860left1[56,56,56,56,56,56,56][0.58,0.58,0.58,0.58,0.58,0.58,0.58]909264left3[12,16,23,34,39,42,45][0.21,0.22,0.24,0.26,0.26,0.27,0.27]909399left10[10,10,12,15,16,18,21][0.2,0.2,0.21,0.21,0.22,0.26,0.27]909418left21[10,11,12,16,21,23,26][0.2,0.21,0.23,0.25,0.29,0.33,0.4]911494left1[12,12,12,12,12,12,12][0.43,0.43,0.43,0.43,0.43,0.43,0.43]912061right1[84,84,84,84,84,84,84][0.42,0.42,0.42,0.42,0.42,0.42,0.42]912080right673[10,10,12,13,16,20,37][0.2,0.2,0.2,0.21,0.23,0.26,0.37]上表中，cnt为计数，srd-percentile即split-readsdepth百分位，srr-percentile即split-readsratio百分位。以第一行为例，该断裂点所在位置为1号染色体左链，908199位点，在6446个样本中有发现，6446个样本中该断裂点处srd的分布为[10,13,16,23,31,40,98]，即srd最大和最小值为98和10，上下十分位为40和13，上下四分位为31和16，中位值为23。srr同理。也即第一行中chr1上pos为908199，side为左的位点出现断裂点的频次为6446。2、查询、判断。将待分析样本中出现的断裂点，在上述数据库中进行检索比对，获得每个断裂点在数据库中的频次及srd/srr百分位信息。由于发生致病性突变的概率较低(通常认为低于千分之一)，可以排除频次超过12的断裂点。即但出现频次≤12次时，则可判断为高风险致病性结构变异。上述频次的计算方法为：由于本案样本数为12385，千分之一即为12。可以理解的，当样本数变化时，该频次也做相应调整，如样本数为20000，则频次相应调整为20。例如，将某一全外显子测序样本的断裂点在数据库中进行检索比对，初步分析得到1544个断裂点，排除频次超过12的断裂点后剩余14个位点，如下表所示。表2.数据库检索比对后得到的高风险致病性结构变异位点表3.高风险致病性结构变异位点情况上表中，omim为该基因在人类孟德尔遗传数据库所对应基因功能注释，hgvs为根据hgvs基因突变命名规则进行的命名。从上述结果可以看出，该样本在ppp1r8等13个基因中均出现了低频的高风险致病性断裂点。五、综合判断。将上述高风险致病性断裂点所在基因信息和srd/srr百分位信息，结合与该样本匹配的患者临床信息，以及断裂点所处基因功能注释综合判断。其中，该样本中fga基因发生断裂，此两个断裂点在数据库中没有检索到，也就是说这是个很罕见的断裂点，其功能缺失会导致先天性纤维蛋白原缺乏症，症状为血液不凝固，出血不止。并且与病人的临床症状可匹配，可以确定fga的断裂就是致病原因。如分析得到的断裂点是在数据库中能检索到的断裂点，则可以根据srd/srr的百分位信息检查该位点是否与数据库中其他样本出现的断裂点类似，以及与数据库中其他出现该断裂点的样本的临床信息做比较，如果srd/srr和临床信息都与数据库中类似，则该位点致病的风险较大，否则是良性位点的可能性更大。即srd/srr百分位信息可辅助进行判断。实施例2一种基因组断裂点识别系统，包括：存储模块：用于获取样本测序下机数据，并进行数据存储；分析模块：按照实施例1所述的基因组断裂点识别方法，对数据进行分析；输出模块：用于将分析得到的高风险致病性结构变异输出。以上述基因组断裂点识别软件系统，对实施例1中的12385例全外显子测序数据进行回顾。回顾结果发现了21例含断裂点的样本，且断裂位置所处基因与样本临床症状匹配，如下表所示。表4.含断裂点的样本情况其余缺乏临床信息以及断裂点所处基因意义不明确的样本多达千余例，这些断裂点也许并非都致病，但可能具有某些生物学功能，是人类基因多态性的组成部分之一，具有研究价值。上述样本中，样本3、7、8、11、12、17同时还进行了cnv及indel检测，其结论与本发明中发现的断裂点分析结果一致。样本13同时还进行了mlpa检测，结论为dmd基因30-43号外显子缺失，与本发明中发现的断裂点分析结果一致。上述结果表明，本发明的基因组断裂点识别方法和系统，可迅速识别基因组中结构变异位点，且可以排除大量假阳性或者人群多态性的断裂点，所得分析结果与临床表现相一致。以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。当前第1页12