CN108304693B - 利用高通量测序数据分析基因融合的方法 - Google Patents

Abstract

本发明公开一种利用高通量测序数据分析基因融合的方法。本发明使用pair‑end数据作为输入，对soft‑clip功能不作要求，而且不必需要设置多个复杂参数，根据比对信息和后续过滤条件，检测融合并保留融合的序列和位置，对于单端数据量1G的双端测序数据来说，单对融合基因的分析仅需不到一小时，且每增加一对融合基因的分析，仅需增加5至10分钟。

Description

利用高通量测序数据分析基因融合的方法

技术领域

本发明通常涉及基因检测领域，特别地涉及利用高通量测序数据分析基因融合的方法。

背景技术

基因融合是指由于染色体易位、插入、缺失、倒置等现象引起两个基因互相融合成为新的基因的过程。基因融合对癌症的发生有重要的影响，而且和单核苷酸多态性(SNP，single-nucleotide polymorphism)一样，它也是诊断和治疗癌症的靶标之一。例如，在肺癌实体瘤中常见的融合基因有ALK-EML4、RET-CCDC6、ROS1-SLC34A2等。

传统的基因融合检测主基于聚合酶链式反应和荧光原位杂交(FISH)两种方法。这两种方法具有准确度高的优点，同时具有通量低、只能定性不能定量、不能检测出融合类型和位置等缺点。二代测序(NGS)具有通量高、价格低、检测范围广等优点，它的成熟在很大程度上加快了基因融合的研究。NGS不仅一次可检出多种融合类型，还可以得到具体的reads支持条数。比起传统检测方法只能得到“有融合”或者“没有融合”两个结论来说，NGS可以得到融合基因、融合位置、融合序列等信息，增加了结果的可靠性和可溯性。

二代测序是基于边合成边测序思想的一种测序方法，带有不同荧光基团的碱基ATCG，在进行PCR合成模板的互补链时，会发出不同颜色的荧光，从而确定基因序列。使用二代测序结果进行基因融合检测的一般步骤是，对原始数据，首先去除低质量的序列和接头，得到clean reads。再把它们比对到参考基因组上，找到比对的位置，然后对结果按照比对位置进行排序，标记由PCR扩增得到的重复序列并除去重复序列的影响，最后使用结构变异(SV)分析软件找到融合基因和断裂点。

目前分析结构变异的软件有三种思路：

read-depth方法，是假设序列映射的深度符合泊松分布，然后查看序列深度的分布符合情况。该思路主要用于检测缺失和重复，其他情况不适用。比如，和正常的二倍体相比，缺失区域的read depth会表现为显著降低，而重复区域的read depth则表现为显著提高。常见的使用read depth方法的软件如CNVnator。

read-pair方法，是利用pair-end的一对reads进行中间空缺长度和序列方向的检测，对跨度和方向不一致的reads，该方法会把它们聚类成为一个组。两个reads间的距离过长说明中间可能存在缺失，距离过短暗示存在仍然插入，方向不一致可能是因为存在插入和随机重复，如果只有一条read能比对到基因组上，说明可能存在新插入。该方法是应用最广泛的方法，常见的检测软件如BreakDancer、HYDRA和SPANNER。

split-read方法，是分割短序列的方法。Splitread是指在比对到参考基因组时包含一段空位的序列。要想应用此方法，需要所用的比对软件有soft-clip功能(如BWA)。当用pair-end模式使用BWA做比对时，如果一条read(记为R1)能比对到参考基因组上，但是另一条read(记为R2)却无法比上的时候，BWA会对R2在R1附近使用更为宽松的Smith-Waterman局部比对策略搜索可能的比对位置。如果这条read只有一部分能够比上，那么BWA会对其进行soft-clip，而这里也往往是包含结构性变异的断点之处。该方法主要用于检测缺失和小的插入，常见软件如pindel。

每种方法都有各自的优势，同时也有一定的局限性。比如Read-depth方法适用性很窄，只能检测拷贝数变异；read-pair方法对重复序列的处理有问题，且片段大小和分布决定了是否能找到断点；split-read方法在重复区域处理不准确，在唯一区域的结果更可靠。每个软件都有自己的侧重分析点，在各个参数中取平衡，目的是最大的分析出结构变异的数量和准确度。由于各个软件算法和参数的不同，不同软件产生的基因融合结果的交集很少，而且一般软件耗时多并消耗大量资源；软件得到的结果很冗杂、没有针对性，很少有我们关注的、重要的特定致病性的基因融合；从比对后的BAM文件到结构变异分析结果，一般软件用时从2小时到一天时间不等。目前应用split-read方法的融合检测软件均要求比对软件有soft-clip功能，对于已经确定的分析流程，需要很大改动，参数也需要重新测试，并且，有soft-clip功能的比对软件对于两个基因的断点附近的序列(一般长几bp)判断不准确，可能导致检测不出基因融合。所以，对于NGS基因融合检测的临床生产来说，目前的方法的适用性偏低。

发明内容

针对现有技术中存在问题，本发明提供一种生物信息分析方法，可以只针对基因融合分析进行生物信息学流程的设计，避免了当前方法得到结果的复杂性，提高了分析的稳定性、实用性，缩短了时间，提高了分析效率。具体地本发明包括以下内容。

一种利用高通量测序数据分析基因融合的方法，其包括以下步骤：

(1)在由多个原始测序序列组成的集合中，针对各原始测序序列从5’端和3’端分别截取掉m个碱基，得到多个待比对序列，其中m为0-20之间的整数；

(2)以包含基因A的序列(例如，基因A的全长序列)为第一参考序列，以包含基因B的序列(例如，基因B的全长序列)为第二参考序列；

(3)将所述多个待比对序列分别与基因组序列进行比对，取出未完全比对的序列标记为R1，并将与其对应的互补序列标记为R2；

(4)取序列R1的前端x个碱基长度的序列标记为R1_x，取序列R1的前端y个碱基长度的序列标记为R1_y，并且取序列R2的前端x个碱基长度的序列标记为R2_x，取序列R2的前端y个碱基长度的序列标记为R2_y，其中x和y分别为20-80之间的整数，如果R1_x与第一参考序列完全匹配，且R2_y与第二参考序列完全匹配，或者R1_y与第二参考序列完全匹配，且R2_x与第一参考序列完全匹配，则将序列R1和序列R2合并得到融合候选序列；

(5)将所述融合候选序列比对到基因组，根据基因A和B在基因组上的位置过滤比对结果，过滤后剩下的融合候选序列作为融合序列，如果得到的融合序列的数量为1以上，则判定所述基因A和所述基因B发生融合，否则判定所述基因A和所述基因B未发生融合。

在某些实施方案中，所述高通量测序数据为双端测序数据，且所述高通量测序数据包括第一数据集和第二数据集。

在某些实施方案中，所述高通量测序数据包括全基因组测序数据、靶向测序数据和全外显子测序数据。

在某些实施方案中，所述各原始测序序列的长度分别为100-300bp。

在某些实施方案中，进一步包括从由多个原始测序序列组成的集合中去除重复的原始测序序列的步骤。

在某些实施方案中，通过将步骤(1)中得到的待比对序列比对到基因组来去除重复的原始测序序列。

在某些实施方案中，其中x与y相等，且分别为所述原始测序序列长度的1/3。

在某些实施方案中，所述基因A和所述基因B来源于同一物种基因组。

在某些实施方案中，R1和R2分别位于第一数据集和第二数据集中的不同数据集。

在某些实施方案中，所述基因A选自ALK、CD74和KIF5B，且所述基因B选自EML4、ROS1和RET。

本发明的方法使用pair-end数据作为输入，对soft-clip功能不作要求，而且不需要设置多个复杂参数，根据比对信息和后续过滤条件，检测融合并保留融合的序列和位置。在可接受的时间范围内，大大增加了敏感性。对于单端数据量1G的双端测序数据来说，单对融合基因的分析仅需不到一小时，且每增加一对融合基因的分析，仅需增加5至10分钟。

具体实施方式

现详细说明本发明的多种示例性实施方式，该详细说明不应认为是对本发明的限制，而应理解为是对本发明的某些方面、特性和实施方案的更详细的描述。

应理解本发明中所述的术语仅仅是为描述特别的实施方式，并非用于限制本发明。另外，对于本发明中的数值范围，应理解为具体公开了该范围的上限和下限以及它们之间的每个中间值。在任何陈述值或陈述范围内的中间值以及任何其他陈述值或在所述范围内的中间值之间的每个较小的范围也包括在本发明内。这些较小范围的上限和下限可独立地包括或排除在范围内。

除非另有说明，否则本文使用的所有技术和科学术语具有本发明所述领域的常规技术人员通常理解的相同含义。虽然本发明仅描述了优选的方法和材料，但是在本发明的实施或测试中也可以使用与本文所述相似或等同的任何方法和材料。本说明书中提到的所有文献通过引用并入，用以公开和描述与所述文献相关的方法和/或材料。在与任何并入的文献冲突时，以本说明书的内容为准。

关于本文中所使用的“包含”、“包括”、“具有”、“含有”等等，均为开放性的用语，即意指包含但不限于。关于本文中所使用的“和/或”，包括所述事物的任一或全部组合。

针对现有分析基因融合方法所存在的问题，本发明致力于精准分析和缩短时间、提高分析效率，从而提出了一种新的利用高通量测序数据分析基因融合的方法，其包括以下5个步骤：

(1)在由多个原始测序序列组成的集合中，针对各原始测序序列从5’端和3’端分别截取掉m个碱基，得到多个待比对序列，其中m为0-20之间的整数。

(2)以包含基因A的序列为第一参考序列，以包含基因B的序列为第二参考序列。

(3)将所述多个待比对序列分别与基因组序列进行比对，取出未完全比对的序列标记为R1，并将与其配对的互补序列标记为R2。

(4)取序列R1的前端x个碱基长度的序列标记为R1_x，取序列R1的前端y个碱基长度的序列标记为R1_y，并且取序列R2的前端x个碱基长度的序列标记为R2_x，取序列R2的前端y个碱基长度的序列标记为R2_y，其中x和y分别为20-80之间的整数，如果R1_x与第一参考序列完全匹配，且R2_y与第二参考序列完全匹配，或者R1_y与第二参考序列完全匹配，且R2_x与第一参考序列完全匹配，则将序列R1和序列R2合并得到融合候选序列。

(5)将所述融合候选序列比对到基因组，根据基因A和B在基因组上的位置过滤比对结果，过滤后剩下的融合候选序列作为融合序列，如果得到的融合序列的数量为1以上，则判定所述基因A和所述基因B发生融合，否则判定所述基因A和所述基因B未发生融合。

下面详细说明各步骤。

第1步骤：

本发明的第1步骤包括在由多个原始测序序列组成的集合中，针对各原始测序序列从5’端和3’端分别截取掉m个碱基，得到多个待比对序列，其中m为0-20之间的整数。

本发明中，原始测序序列(有时也称作原始序列)为二代测序的序列，例如，高通量测序序列。各序列的长度一般在100-300bp之间，优选120-200bp之间，更优选150-200bp之间。

本发明中，所述高通量测序数据是指包括原始测序序列信息的数据。优选地，本发明中的原始序列为双端测序序列。更优选地，高通量测序数据包括第一数据集和第二数据集，且第一数据集包含双端测序中一端测序序列的数据，和第二数据集包含双端测序序列中另一端测序序列的数据。在某些实施方案中，第一数据集和第二数据集储存在两个不同的FASTQ文件中。优选地，所述高通量测序数据包括全基因组测序数据、靶向测序数据和全外显子测序数据。

本发明的方法中，需要对各原始测序序列进行处理，包括从5’端和3’端分别截取掉m个碱基，其中m为0-20之间的整数。在某些实施方案中，m可为1-15，更优选5-10，从而提高检测的灵敏度。在某些实施方案中，m可为10-20，从而去除测序时引入的标签序列等。在某些实施方案中，m可以为0，即可以使用直接获取的原始序列的5’端和3’端序列。

在某些实施方案中，本发明的方法进一步包括从由多个原始测序序列组成的集合中去除重复的原始测序序列，由此可提高检测的效率。可使用本领域内已知的任何手段去除重复的原始测序序列，例如，通过将第1步骤中得到的待比对序列比对到基因组来去除重复的原始测序序列。例如，当两个不同ID的序列均比对到基因组的相同位置，且序列完全相同时，则保留一对质量高的序列，去除重复的序列。

第2步骤：

本发明的第2步骤包括以包含基因A的序列为第一参考序列，以包含基因B的序列为第二参考序列。优选地，以基因A的全长序列为第一参考序列，以基因B的全长序列为第二参考序列。其中基因A和基因B是指待检测是否发生融合的两个基因。例如，当需要检测ALK基因和EML4基因是否发生融合时，ALK基因可以作为基因A，此时EML4基因作为基因B。另外，可选地，EML4基因可以作为基因A，此时ALK基因作为基因B。基因A和基因B的来源不特别限定，优选两个基因均来源于同一物种，即两个基因包含于同一基因组中。在某些实施方案中，基因A和基因B位于同一染色体上。在某些实施方案中，基因A和基因B位于不同染色体上。

第3步骤：

本发明的第3步骤包括将多个待比对序列分别与基因组序列进行比对，取出未完全比对的序列标记为R1，并将与其对应的互补序列标记为R2。

本发明中，R1和R2为分别从两端测序时所得结果对应的序列。即，R1和R2具有相同的ID，R1和R2是否属于相同的ID的判断是本领域内已知的。例如，选取未比对上基因组的序列ID，然后通过该ID从第一数据集和第二数据集中挑出R1和R2。因此，当R1包含于第一数据集时，R2包含于第二数据集。当R2包含于第一数据集时，R1包含于第二数据集。需要说明的中，R1和R2所表示的序列并不一定完全互补。

第4步骤：

本发明的第4步骤包括取序列R1的前端x个碱基长度的序列标记为R1_x，取序列R1的前端y个碱基长度的序列标记为R1_y，并且取序列R2的前端x个碱基长度的序列标记为R2_x，取序列R2的前端y个碱基长度的序列标记为R2_y。其中所述“前端”优选是指序列的5’端。例如，取序列R1的前端x个碱基长度的序列是指从序列R1的5’端起截取x个bp的碱基序列。

在某些实施方案中，x和y分别为20-80之间的整数，优选25-70，更优选30-60，例如50。如果长度过长则可能增加假阴性结果。如果过短则检测灵敏度会降低。在某些实施方案中，x和y的数值为原始测序序列长度的1/3，此时可同时保证检测灵敏度及特异性。x和y的数值可以相同，也可以不同。

本发明的第4步骤还包括如果R1_x与第一参考序列完全匹配，且R2_y与第二参考序列完全匹配，或者R1_y与第二参考序列完全匹配，且R2_x与第一参考序列完全匹配，则将相同ID的序列R1和序列R2合并得到融合候选序列。

第5步骤：

本发明的第5步骤包括将融合候选序列比对到基因组，根据基因A和B在基因组上的位置过滤比对结果，过滤后剩下的融合候选序列作为融合序列，如果得到的融合序列的数量为1以上，则判定所述基因A和所述基因B发生融合，否则判定所述基因A和所述基因B未发生融合。

如上所述虽然本发明仅详细描述了上述5个步骤，但是本发明的方法并不仅局限于所述5个步骤，还可包括其他步骤。另外，5个步骤之间的顺序并不特别限定，可以依次进行第(1)至第(5)步骤，也可以任意的顺序进行所述5个步骤。只要能够实现本发明的目的，则对于步骤之间的顺序不特别限定。

本发明的利用高通量测序数据分析基因融合的方法可用于检测各种基因之间是否发生融合。对于待检测的基因A和B的类型不特别限定。例如，所述基因A和基因B可选自ALK、CD74、KIF5B、EML4、ROS1和RET。本发明的方法可同时检测是否存在多个融合。多个融合是指在同一测序数据集合内，或同一样本得到的数据中存在多对基因的融合。

实施例

一、样本信息

选取实际肺癌检测的探针捕获的双端(Paired-end)测序且读长为151bp的数据Sample作为分析举例，运用本发明的方法及较常用的Cutadapt去接头算法、BWA(aln)比对算法，对该测序数据进行分析检测肺癌中最常见的ALK-EML4融合。

二、操作步骤

1)测序后，使用cutadapt去除原始reads中的接头、质量低的碱基得到cleanreads，再使用bwa aln算法使用双端模式(pair-end)进行比对到hg19人类基因组参考序列上，进行去重后得到比对结果，记为alignment。

2)从hg19参考序列上，截取ALK和EML4的序列，记为ALK.fa和EML4.fa。

3)从alignment中选择未比对上的序列R1和R2，分别截取前50bp，记为R1_50，R2_50。

4)R1_50比对到ALK.fa上，R2_50比对到EML4.fa上，得到结果严格过滤，错配设为0，且是唯一比对。

5)同理R1_50比对到EML4.fa上，R2_50比对到ALK.fa上，严格过滤，即，错配设为0，且是唯一比对。

6)得到两种情况都比对上的reads ID信息，然后根据reads ID从clean reads中选出双端序列合并序列，记为candidate_reads。

7)使用BLAST将candidate_reads比对到hg19参考序列上通过ALK和EML4在hg19上的位置信息进行过滤，从而得到真正的融合序列(记为fusion_seq)和位置信息。

8)结果判定：fusion_seq中序列数量大于或等于1，则有目标基因的融合发生；fusion_seq中序列数量为0，则无目标基因的融合发生。

三、结果总结

测试样本均选用实际肺癌探针捕获的双端(Paired-end)测序且读长为151bp的测序数据，共50例，测序数据量从40M到1G大小不等，其中10例用ARMS法验证为ALK融合。验证使用了艾德的试剂盒，可检测ALK-EML4，ALK-KIF5B，ALK-TFG和ALK-KLC1四种关于ALK的融合类型。

为了测试本发明的实用性，将本发明的方法与CREST+ANNOVAR和FACTERA两个方法在敏感性、特异性和运行时间上做了比对，总结如下表1。由于发明中步骤一用时相等，故表格中所示的时间均为去除步骤一后所用时间。对每一个样本的总结请见附表。

表1本发明的方法与现有方法的比较数据

表2本发明的方法与现有方法的比较数据汇总结果

总结	本发明的方法	CREST+ANNOVAR	FACTERA
				敏感性	90.00％	30.00％	10.00％
特异性	92.00％	100.00％	100.00％
				平均用时(mins)	12.43	14.87	0.76

如表格所见，CREST+ANNOVAR在用时多于本发明的情况下，敏感性远低于本发明。FACTERA虽然用时极短，但是它的敏感性仅有10％，基本遗漏掉了所有有融合的阳性样本，结果不可靠。本发明在可接受的时间内，取得的敏感性远远超过其他两个软件，虽然在特异性上仍有些不足，但是可以通过做验证来确定融合的发生。

在不背离本发明的范围或精神的情况下，可对本发明说明书的具体实施方式做多种改进和变化，这对本领域技术人员而言是显而易见的。由本发明的说明书得到的其他实施方式对技术人员而言是显而易见得的。本申请说明书和实施例仅是示例性的。

Claims (7)

1.一种利用全基因组高通量测序数据分析基因融合的方法，其包括以下步骤：

(1)在由多个原始测序序列组成的集合中，针对各原始测序序列从5’端和3’端分别截取掉m个碱基，得到多个待比对序列，其中m为0-20之间的整数；

(2)以包含基因A的序列为第一参考序列，以包含基因B的序列为第二参考序列；

(3)将所述多个待比对序列分别与基因组序列进行比对，取出未完全比对的序列标记为R1，并将与其对应的互补序列标记为R2；

(4)取序列R1的前端x个碱基长度的序列标记为R1_x，取序列R1的前端y个碱基长度的序列标记为R1_y，并且取序列R2的前端x个碱基长度的序列标记为R2_x，取序列R2的前端y个碱基长度的序列标记为R2_y，其中x和y分别为20-80之间的整数，如果R1_x与第一参考序列完全匹配，且R2_y与第二参考序列完全匹配，或者R1_y与第二参考序列完全匹配，且R2_x与第一参考序列完全匹配，则将序列R1和序列R2合并得到融合候选序列，其中所述前端是指序列的5’端；

(5)将所述融合候选序列比对到基因组，根据基因A和基因B在基因组上的位置过滤比对结果，过滤后剩下的融合候选序列作为融合序列，如果得到的融合序列的数量为1以上，则判定所述基因A和所述基因B发生融合，否则判定所述基因A和所述基因B未发生融合；

其中所述高通量测序数据为双端测序数据，且所述高通量测序数据包括第一数据集和第二数据集，R1位于第一数据集，且R2位于第二数据集；或者R1位于第二数据集，且R2位于第一数据集。

2.根据权利要求1所述的方法，其中所述各原始测序序列的长度分别为100-300bp。

3.根据权利要求1所述的方法，进一步包括从由多个原始测序序列组成的集合中去除重复的原始测序序列的步骤。

4.根据权利要求3所述的方法，其中通过将步骤(1)中得到的待比对序列比对到基因组来去除重复的原始测序序列。

5.根据权利要求1所述的方法，其中x与y相等，且分别为所述原始测序序列长度的1/3。

6.根据权利要求1所述的方法，其中所述基因A和所述基因B来源于同一物种基因组。

7.根据权利要求1所述的方法，其中所述基因A选自ALK、CD74和KIF5B，且所述基因B选自EML4、ROS1和RET。