CN108830044A - 用于检测癌症样本基因融合的检测方法和装置 - Google Patents

用于检测癌症样本基因融合的检测方法和装置 Download PDF

Info

Publication number
CN108830044A
CN108830044A CN201810570027.2A CN201810570027A CN108830044A CN 108830044 A CN108830044 A CN 108830044A CN 201810570027 A CN201810570027 A CN 201810570027A CN 108830044 A CN108830044 A CN 108830044A
Authority
CN
China
Prior art keywords
cluster
long
read
genome
breakpoint
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.)
Granted
Application number
CN201810570027.2A
Other languages
English (en)
Other versions
CN108830044B (zh
Inventor
黄炳顶
金林芳
戴珩
史耀舟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xukang medical technology (Suzhou) Co., Ltd
Original Assignee
Shanghai Whale Boat Gene Technology Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shanghai Whale Boat Gene Technology Co Ltd filed Critical Shanghai Whale Boat Gene Technology Co Ltd
Priority to CN201810570027.2A priority Critical patent/CN108830044B/zh
Publication of CN108830044A publication Critical patent/CN108830044A/zh
Application granted granted Critical
Publication of CN108830044B publication Critical patent/CN108830044B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了一种用于检测癌症样本基因融合的检测方法和装置。本发明方法包括步骤:(M1)将短序列与基因组进行比对,从而获得不正常比对的读长对;(M2)对不正常比对读长对(read pair)进行收集并聚类,从而获得不正常比对读长对簇;(M3)在不正常比对读长对簇所覆盖的基因组区域内,对断点对进行检测,从而获得唯一的融合事件断点位置对;(M4)对所述断点位置对所涉及的融合事件的进行注释和质控检测;和(M5)输出上一步骤的检测结果。本发明检测方法充分利用了双端测序数据的信息,只需一次比对,不需组装,不需重比对,确定的融合区域位置分界清晰,断点位置准确,结果稳定,灵敏度,特异性高。

Description

用于检测癌症样本基因融合的检测方法和装置
技术领域
本发明涉及生物信息领域,更具体地,涉及癌症基因组目标区域核酸片段融合的检测领域,尤其涉及用于检测癌症样本中基因融合的检测方法和装置。
背景技术
基因融合指的是生物基因组上两个在物理位置上分离的基因嵌合为一种新的基因的现象。
有文献报道证明,基因融合现象和癌症的发生有相关性,基因融合事件成为一些癌症的分子标记。比如在慢性髓细胞性白血病、前列腺癌、肺癌、恶性胶质瘤、肝癌、乳腺癌、皮肤癌、淋巴癌和肉瘤样本中发生的BCR-ABL1基因融合事件。这些融合事件是肿瘤诊断、预测、治疗过程中良好的靶标。一些典型的融合事件,比如BCR-ABL,PML-RAR和EML4-ALK都已经是靶向癌症治疗中的标准代表。
随着二代测序技术的发展,癌症细胞基因组上的部分或全部融合事件可以通过对DNA或者RNA测序数据集的挖掘来识别。数据集中的测序reads可以是单端或双端的。但双端reads的数据由于其反映了插入片段双端的信息,使其在鉴别中可以得到更高特异性的结果。
目前主要的高通量测序数据的基因融合的检测方法包括两类:
第一类是基于双末端(Pair End,PE)关系的检测方法:由于高通量测序文库构建时插入大小是确定的,那么如果根据PE测序所得到的序列(reads)的比对位置所判定的插入大小,显著偏离了测序文库构建时的插入大小的平均值(例如,一对reads分别比对到不同的染色体上),则有可能是发生了基因融合。此类方法主要利用这样的双末端关系来判断基因重排导致的异常双末端比对序列(reads),根据这些序列(reads)的比对位置、插入大小等信息来检测融合。
第二类是基于截断比对(split-mapping)的检测方法:主要利用非完全比对序列(soft-clipped reads—软截断序列)的序列信息进行融合断点识别,然后对断点上下游比对的reads做聚类分析以及拼接组装,最后重新对序列做定位分析,进而检测基因融合现象。
上述两种检测方法中,基于PE关系的检测方法,是根据异常双末端比对序列(reads)的信息来进行融合检测(例如,Break Dancer),只能大致给出融合位置,而不能确定准确的断点信息,并且此类方法的假阳性较高。
此外,目前的其他一些方法还存在一些其他明显的缺点,例如在不同数据集里测试得到的结果不稳定,会消耗大量的计算资源(计算时间+内存占用),特异性较低和灵敏度较低(例如,会失去一些真阳性事件,并出现许多假阳性事件)。
在融合检测算法中,融合断点的鉴别是核心任务。目前,鉴别融合断点的策略分为两大类:比对后处理和组装后处理。组装后处理耗时,耗内存,相比比对后处理存在劣势,比对后处理策略成为融合鉴定软件采用的主流策略。比对后,从比对结果文件中找到两端reads比对到较远位置的pair集合。从找到的pair集合中确定潜在融合区域范围,再在此范围中利用split reads确定融合位点的准确位置,是比较准确合理的策略。但在具体实施过程中,如何简化步骤,降低计算资源消耗,特别是低覆盖度样本种融合事件的准确检测成为本领域面临的一个重大难题。
现有的技术算法比较优秀的策略在发现不正常比对的pair集合后再确定候选融合区域集合的过程中,会存在区域元素距离较近,分界模糊的问题。在确定区域后,在确定断点准确位置的过程中,采用局部重比对方法,虽然合理,但消耗计算资源较大,需要进一步改进。
综上所述,现有的技术中检测基因融合的方法几乎都存在假阳性率高,流程耗时长,低融合频率情况下敏感性低的缺陷,不能满足临检生产的敏感性要求。
因此本领域迫切需要开发一种具有高稳定性、高灵敏度和高特异性,且计算资源消耗低和计算速度快的检测基因融合的检测方法和装置。
发明内容
本发明的目的就是提供一种具有高稳定性、高灵敏度和高特异性,且计算资源消耗低和计算速度快的检测基因融合的检测方法和装置。
本发明的第一方面,提供了一种检测基因组嵌合片段的检测方法,所述检测方法包括:
(M1)将测序获得的短序列与基因组进行比对,从而获得正常比对的读长对和不正常比对的读长对;
(M2)对不正常比对读长对(read pair)进行收集,并对所述不正常比对读长对进行聚类,从而获得不正常比对读长对簇;
(M3)在不正常比对读长对簇所覆盖的基因组区域内,对断点对进行检测,从而获得位于所述的被覆盖的基因组区域内的唯一的融合事件断点位置对;
(M4)对所述断点位置对所涉及的融合事件的进行注释和质控检测,其中如果满足质控标准则提示所述断点位置所涉及的融合事件是真的,而不满足质控标准则提示所述断点位置所涉及的融合事件是非真的;和
(M5)输出上一步骤的检测结果。
在另一优选例中,在步骤(M5)中,输出经质控检验的最终真融合事件相关信息。
在另一优选例中,在步骤(M1)中,在将所述短序列比对到基因组,从而获取测序读长比对信息文件。
在另一优选例中,所述的基因组包括动物、植物、微生物、或病毒的基因组。
在另一优选例中,所述的基因组为哺乳动物的基因组。
在另一优选例中,所述的哺乳动物包括人和非人哺乳动物(如啮齿动物、非人灵长动物)。
在另一优选例中,在步骤(M1)中,将所述测序读长使用BWA软件设置合理参数,比对到人类参考基因组hg19上。
在另一优选例中,在步骤(M1)中,将双端读长比对到相隔大于一定距离的信息储存到比对结果文件里的记录里,并标记分两部分比对到基因组上的读长,标记读长未比对到基因组的部分为软剪切,储存到比对记录里的CIGAR字符串中。
在另一优选例中,在所述的比对过程所得到的比对记录文件中,将所涉及的读长分为5类:
1)全部比对到基因组上的读长,视作A类读长。
2)一部分比对到基因组上,但另一部分未比对到基因组上的的读长,视作B类读长。
3)一部分比对到基因组上,另一部分也比对到基因组上的读长,视作C类读长;其中较长的部分称为主要比对部分,视作C1类,较短的部分称为次级比对部分,视作C2类。
4)读长对两端都比对到基因组上,读长对两端比对到相隔较远的位置上,或者读长对两端比对方向相同(正常情况读长对两端比对方向应相反)但读长对两端读长都不属于上述C类,视作D类,也叫不正常比对读长对(Discordant read pair)
5)不符合以上A,B,C,D类任何一类的读长视作E类。
在另一优选例中,在步骤(M2)中,即在不正常比对读长对的收集中,将比对短序列读长到基因组得到的所有远距离比对到基因组的D类读长对挑选出来。
在另一优选例中,所述的挑选原则如下:
1)两端读长均比对到基因组上;
2)两端读长质量都满足质控要求;
3)两端读长都不是PCR Duplicates;和
4)两端读长没有上述的C2类读长。
在另一优选例中,将按上述挑选原则所得到的不正常比对读长对收集起来,并将远距离比对读长对按照两端是否在同一染色体上,两端比对方向是否同向的情况将收集到的不正常比对读长对集合分为四类:
两端异向且比对到不同染色体上;
两端异向且比对到相同染色体上;
两端同向且比对到不同染色体上;
两端同向且比对到相同染色体上;
并在这四类产生的四个组的各自范围内,再依据组内不正常比对读长对的基因组上比对位置对它们进行无监督聚类。
在另一优选例中,在所述的聚类步骤之前,对不正常比对读长对的两端读长在基因组上坐标进行一维转换。
在另一优选例中,所述的一维转换指,将物种的序列(如各染色体序列)条目按照编号排序,并把这些序列依次连接,变为一维的线性长轴。
在另一优选例中,在所述的聚类步骤中,用选自下组的聚类方法进行聚类:排序聚类法和凝聚层次聚类法。
在另一优选例中,所述排序聚类法包括以下步骤:
1)将不正常比对读长对两端读长的比对位置二维坐标(x,y)集合内的元素按照x坐标的顺序从小到大排序;
2)定义上步排序后第一个元素的x坐标为previous position,并记录它为当前cluster的第一个元素;
3)从第二个元素开始遍历排序后的集合,比较当前元素的x坐标值是否与previous position相差设定阈值以内的数值,是的话则将此元素记录为当前cluster的一个新元素,并接着遍历下一个元素,对下一个元素重复此步骤,否的话则进入下一步;
4)判断当前cluster内元素个数是否大于等于已设定好的阈值(一般为2),如果是,则保留此cluster以及此cluster内的元素,并将cluster id赋给其内的每个元素,如果否,则丢弃此cluster以及此cluster内的元素;判断结束后,清空当前cluster,更新cluster id,将当前元素的X坐标设置为previous position;
5)遍历完一遍组内元素后,再将剩下的元素集合按照点的y坐标从小到大排列,重复3)和4),确定最后留下的参与成簇的不正常比对读长对元素集合以及它们所归属的簇id。
在另一优选例中,所述凝聚层次聚类法包括以下步骤:
1)将原始点集内的每个点初始化为单点簇;
2)遍历簇集,每次合并两个最近的簇为一个簇,这样离得相对较近的点都会逐渐的纳入到同一个簇里;其中簇元素之间的距离有三种:single-linkage,complete-linkage,average-linkage;其中,single-linkage为两个簇内对象之间的最小距离;complete-linkage为两个簇内对象之间的最大距离;complete-linkage为两个簇内对象之间的平均距离;
3)重复2)步骤,直到当前簇集内所有元素簇的簇间距离都大于设定的距离值,停止合并点簇步骤;
4)遍历当前合并停止后生成的簇集,过滤簇元素内点元素个数小于设定值的簇元素,将每个簇赋予一个簇元素ID;
5)遍历原始点集,记录当前点所属的簇元素ID;这样就得到一个符合需求的簇集合。
在另一优选例中,在所述远距离比对读长对簇覆盖基因组区域内断裂位置对的检测步骤中,此模块具体分为下面几个步骤:
a.遍历簇集,以每个簇的中心位置的x和y坐标分别向左右扩展一定距离,形成一个矩形范围;
b.遍历上步扩展的成对区域内所有的比对记录,找到这些比对记录中splitreads的比对记录,其中split reads是一部分比对到基因组上,另一部分也比对到基因组上的的读长;
c.分析上步获取的split reads比对记录,分析这些比对记录中读长两部分的比对模式(通过比较Cigar和SA tag的值),潜在断点就在两条记录中所描述的读长中部match和soft clipping的交界处,将读长的这个交界处对应的参考基因组序列上对应的坐标计算出来;
d.遍历成对区域中的所有的split reads记录,计算每一条记录所对应的断点位置,再按照读长名字将在两端都出现的读长挑出来组合在一起,组成一个断点对;其中,将断点对按照match部分的方向组合分为四类:(左,左),(左,右),(右,左)(右,右),同类内的断点对集合进行频次排序,挑选频次最高的那个断点对组合,最后在四类的频次最高的断点对组合中选取频次最高的断点对组合,作为最后这个不正常比对读长簇所覆盖区域中得到的唯一的融合事件断点位置对;
e)记录支持此断点对所代表的融合事件的支持不正常比对读长对的个数和splitreads个数,这部分信息会分别记录在结果文件的N_Discordant_Pairs字段和N_Split_Reads中;其中若支持的split reads个数不大于设定好的个数阈值,则不会被输出。
在另一优选例中,在所述远距离比对读长对簇覆盖基因组区域内断裂位置对的检测步骤后,还包括嵌合事件的鉴别步骤:
进一步分析断点左右的读长比对模式,从而鉴别此嵌合事件的嵌合方式。
在另一优选例中,所述的注释包含对选自下组的一个或多个或全部信息进行注释:
a)注释融合事件所涉及的两个断点所出的基因区域,外显子位置;
b)注释融合事件所涉及的两个断点附近是否为重复序列区域,如果是的话,则判定此融合事件为假阳性,在最终结果中过滤此事件,不报告出来;
c)如果测序样本物种为人,则会进而以当前融合事件的两个断点所处的基因融合对搜索热点人类融合基因对数据库,如果有匹配结果则加上hotspot的标签。
在另一优选例中,所述热点人类融合基因对数据库是装置内置的或在线的。
在另一优选例中,所述热点人类融合基因对数据库是定期更新的或实时更新的。
在另一优选例中,在所述的输出步骤中,将注释步骤中所得到的可能的符合质控标准的融合事件相关的选自下组的信息作为表格条目的方式写到结果文件中:
a)两个断点的坐标位置;
b)两个断点各自所处的基因及外显子位置注释(基因名以及外显子编号信息);
c)此融合事件的类型,包含颠换,易位和未知类型(inversion,translocation,unknown);
d)此融合事件的支持远距离比对读长对个数,断裂读长数目,此融合事件的等位基因频率;
e)融合事件所涉及的融合基因对与热点融合基因对数据库之间的比对结果,如果有match,则有hotspot标签。
在另一优选例中,所述的热点融合基因对数据库为人的热点融合基因对数据库。
在本发明的第二方面,提供了一种检测基因组嵌合片段的检测装置,所述检测装置包括:
(D1)序列比对模块,所述模块用于将测序的短序列比对到一预定的基因组;
(D2)不正常比对读长对的聚类模块,所述模块用于对不正常比对读长对(readpair)进行收集和进行聚类分析;
(D3)断点对的检测模块,所述模块用于在不正常比对读长对簇覆盖基因组区域内,对断点对进行检测,从而获得位于所述的被覆盖的基因组区域内的唯一的融合事件断点位置对;
(D4)注释-质控模块,所述模块对断点位置对所涉及的融合事件进行注释和进行质控处理,其中如果满足质控标准则提示所述断点位置所涉及的融合事件是真的,而不满足质控标准则提示所述断点位置所涉及的融合事件是非真的;和
(D5)输出模块,所述模块用于输出检测结果,所述检测结果包括最终真融合事件相关信息。
在另一优选例中,所述输出模块输出经质控检验的最终真融合事件相关信息。
在另一优选例中,序列比对模块将所述短序列比对到基因组,从而获取测序读长比对信息文件。
在另一优选例中,所述的基因组包括动物、植物、微生物、或病毒的基因组。
在另一优选例中,所述的基因组为哺乳动物的基因组。
在另一优选例中,所述的哺乳动物包括人和非人哺乳动物(如啮齿动物、非人灵长动物)。
在另一优选例中,在所述的比对过程所得到的比对记录文件中,将所涉及的读长分为如上定义5类,即A类读长、B类读长、C类读长、D类读长(即不正常比对读长对和E类读长。
在另一优选例中,所述的不正常比对读长对的聚类模块对将比对短序列读长到基因组得到的所有远距离比对到基因组的D类读长对挑选出来。
在另一优选例中,所述检测装置还包括:一个一维转换模块,所述模块用于对不正常比对读长对的两端读长在基因组上坐标进行一维转换。
在另一优选例中,所述聚类模块用选自下组的聚类方法进行聚类:排序聚类法和凝聚层次聚类法。
在另一优选例中,所述的断点对的检测模块执行以下步骤,从而获得断点对信息:
a.遍历簇集,以每个簇的中心位置的x和y坐标分别向左右扩展一定距离,形成一个矩形范围;
b.遍历上步扩展的成对区域内所有的比对记录,找到这些比对记录中splitreads的比对记录,其中split reads是一部分比对到基因组上,另一部分也比对到基因组上的的读长;
c.分析上步获取的split reads比对记录,分析这些比对记录中读长两部分的比对模式(通过比较Cigar和SA tag的值),潜在断点就在两条记录中所描述的读长中部match和soft clipping的交界处,将读长的这个交界处对应的参考基因组序列上对应的坐标计算出来;
d.遍历成对区域中的所有的split reads记录,计算每一条记录所对应的断点位置,再按照读长名字将在两端都出现的读长挑出来组合在一起,组成一个断点对;其中,将断点对按照match部分的方向组合分为四类:(左,左),(左,右),(右,左)(右,右),同类内的断点对集合进行频次排序,挑选频次最高的那个断点对组合,最后在四类的频次最高的断点对组合中选取频次最高的断点对组合,作为最后这个不正常比对读长簇所覆盖区域中得到的唯一的融合事件断点位置对;
e)记录支持此断点对所代表的融合事件的支持不正常比对读长对的个数和splitreads个数,这部分信息会分别记录在结果文件的N_Discordant_Pairs字段和N_Split_Reads中;其中若支持的split reads个数不大于设定好的个数阈值,则不会被输出。
在另一优选例中,所述的注释-质控模块对选自下组的一个或多个或全部信息进行注释:
a)注释融合事件所涉及的两个断点所出的基因区域,外显子位置;
b)注释融合事件所涉及的两个断点附近是否为重复序列区域,如果是的话,则判定此融合事件为假阳性,在最终结果中过滤此事件,不报告出来;
c)如果测序样本物种为人,则会进而以当前融合事件的两个断点所处的基因融合对搜索热点人类融合基因对数据库,如果有匹配结果则加上hotspot的标签。
在另一优选例中,所述的检测装置还包括一内置的储存器,所述储存器储存有热点融合基因对数据库(如人类融合基因对数据库)。
在另一优选例中,所述的输出模块将注释步骤中所得到的可能的符合质控标准的融合事件相关的选自下组的信息作为表格条目的方式写到结果文件中:
a)两个断点的坐标位置;
b)两个断点各自所处的基因及外显子位置注释(基因名以及外显子编号信息);
c)此融合事件的类型,包含颠换,易位和未知类型(inversion,translocation,unknown);
d)此融合事件的支持远距离比对读长对个数,断裂读长数目,此融合事件的等位基因频率;
e)融合事件所涉及的融合基因对与热点融合基因对数据库之间的比对结果,如果有match,则有hotspot标签。
在另一优选例中,所述的热点融合基因对数据库为人的热点融合基因对数据库。
应理解,在本发明范围内中,本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合,从而构成新的或优选的技术方案。限于篇幅,在此不再一一累述。
附图说明
图1显示了融合基因的形成以及“不正常比对读长对”(Discordant Reads Pair)和“断裂读长”(split reads,SR)。
图2显示了本发明用于检测癌症样本DNA融合位点事件的方法的整体流程示意图。
图3显示了本发明一个实施例中基因融合事件的检测方法流程图。
图4显示了本发明一个对比例中基因融合事件的检测方法流程图。
图5为基因组可视化结果的示意图。
具体实施方式
本发明人经过广泛而深入的研究,首次开发了一种高稳定性、高灵敏度和高特异性,且计算资源消耗低和计算速度快的检测基因融合的检测方法和装置。与已有的方法相比,本发明充分利用了双端测序数据的信息,只需一次比对,不需组装,不需重比对,确定的融合区域位置分界清晰,断点位置准确,结果稳定,灵敏度,特异性高。在此基础上完成了本发明。
术语
本发明中所述“参考基因组”是指所带侧样品相应的物种中已经公开发表的全基因组序列数据;
“读长”指的是测序平台测序过程一次测序的长度,也代指测序产生的序列内容;
“双端测序”指的是在构建待测DNA文库时在两端的接头上都加上测序引物结合位点,在第一轮测序完成后,去除第一轮测序的模板链,用对读测序模块(Paired-EndModule)引导互补链在原位置再生和扩增,以达到第二轮测序所用的模板量,进行第二轮互补链的合成测序,一个插入片段分子一个周期的测序会产生一对两端的测序读长(读长对);
“目标区域捕获测序”指的是针对感兴趣的基因组区域设计特异性探针,将其与基因组DNA进行杂交,富集目标基因组区域的DNA片段,而后利用高通量测序技术进行测序的检测方法;
“短序列比对”指的是将二代测序产生的高通量的短序列读长通过一定的算法匹配到参考基因组的某些区域上。本装置中限制配套使用bwa这个开源软件。
“不正常比对读长对”(Discordant Reads Pair)指的是双末端测序产生的成对读长在比对到基因组上后,算出的其所表示的文库插入片段分子大小明显大于实际建库插入片段长度分布的读长对或者两端读长比对的方向同向,如图1所示;
“断裂读长”(split reads,SR)指的是读长在比对到基因组序列上以后,一部分比对到基因组上的A区域,一部分比对到基因组的B区域(未与A相邻的区域)。
“cigar”是一种在比对结果文件里每条记录中表示读长在基因组上比对方式的字符串,通过它和记录中同时存在的读长比对起始位置信息的结合可以推断出读长上每个碱基在基因组上的匹配方式和匹配位置。Cigar字符串的格式一般为以“数字+匹配方式”为单位的一个或多个组合,匹配方式一般以约定好的特定字符集合里的任一一个字符表示,字符一般是“MSDIPH”当中的任意一个,bwa的结果文件中一般只会出现“MSID”这四种字符。其中,“M”代表“match or mismatch”,“S”代表“soft clipping”(长片段的未匹配到基因组上的部分),“D”代表在参考基因组序列上的删除,“I”代表参考基因组序列上的插入。举一个cigar字符串的例子,如图5的基因组可视化软件的截图,其中图中最底部的字母序列是人基因组DNA序列,字母序列上面堆叠的条形是测序产生的读长,举选中的那个红色的读长为例,此读长比对起始位点坐标为chr2:29447932,cigar为“83M67S”,代表有83个碱基与参考序列匹配match(即它左半部分没有字符标识的部分),67个碱基与它对应的参考序列不匹配(即它右半部分有字符标识的部分)。
“SA tag”即“secondary Alignment tag”的简称,次级比对记录标签。当一个读长分为两部分比对到参考基因组的两个相隔较远的区域时,那么这个读长两部分的比对会形成两个比对记录,较长的部分叫做“Primary Alignment”,较短的叫做“SecondaryAlignment”,此读长也叫作split reads,断裂读长。bwa会在比对记录的信息中加入SA tag来标识发生分裂比对的读长,两条记录中都在各自在SA tag中记录另一部分的比对方式,SA tag的内容格式为“SA:start_chr,start_pos,strand,cigar,mapping_quality,”,比如附图5中红色条带代表的补偿比对记录中就带有“SA:chr2,29448026,-,67M83S,60,0;”,代表此读长的另一部分(即上述右侧67S里未匹配的67个碱基部分)匹配到基因组的起始位点是2号染色体的29448026坐标,匹配方式是“67M83S”,匹配到负链(即此部分5'到3'的方向与是逆染色体坐标增长方向的),匹配质量是60。在本发明中,此标签内容是发现断点位置模块中的利用的关键信息。
检测方法和模块
为了解决这些缺陷,在本发明的一种典型的实施方式中,提供了一种基因融合的检测装置和方法。
典型地,参见图2,本发明提供一种检测基因组嵌合片段的检测方法,所述的方法主要包括步骤:将短序列比对到基因组,检测远距离比对读长对簇,根据得到的读长对簇找到断裂位置和注释断点。
本发明还提供相应的用于检测基因组嵌合片段的检测装置,它包括:短序列比对到基因组的模块,不正常比对读长对的收集和分组模块,不正常比对读长对比对位置的二维凝聚层次聚类模块,远距离比对读长对簇覆盖基因组区域内断裂位置对的检测模块,断裂位置对所涉及的融合事件的注释模块和最终真融合事件相关信息的输出模块。
不正常比对读长对的收集
在本发明中,需要对不正常比对读长对进行收集。
在该收集步骤中,将比对短序列读长到基因组得到的所有远距离比对到基因组的D类读长对挑选出来。
在本发明中,一种代表性的优选挑选原则包括如下条件:
1)两端读长均比对到基因组上;
2)两端读长质量都满足质控要求;
3)两端读长都不是PCR Duplicates;和
两端读长没有上述的C2类读长。
不正常比对读长对的分组
在本发明中,对于按照上述挑选原则所得到的不正常比对读长对收集起来,并将远距离比对读长对按照两端是否在同一染色体上,两端比对方向是否同向的情况将收集到的不正常比对读长对集合分为四类:
两端异向且比对到不同染色体上;
两端异向且比对到相同染色体上;
两端同向且比对到不同染色体上;
两端同向且比对到相同染色体上;
并在这四类产生的四个组的各自范围内,再依据组内不正常比对读长对的基因组上比对位置对它们进行无监督聚类。
在本发明中,通过上述分组处理,可以更精确地得到潜在断点可能存在的簇,不正常比对读长对的组内类型一致也为后来鉴定融合类型提供更可靠的依据。
一维转换
在本发明中,在聚类步骤之前,宜对不正常比对读长对的两端读长在基因组上坐标进行一维转换。
在本发明中,所述的“一维转换”指,将物种的序列(如各染色体序列)条目按照编号排序,并把这些序列依次连接,变为一维的线性长轴。
虽然读长对两端读长的比对位置都是按照“染色体编号:比对起始点在该染色体上所处的坐标数值”这种方式记录的。但是,这种格式不方便进行下一步的不正常比对读长对位置聚类步骤。于是,在本发明中,将所述物种的序列条目按照编号排序,并把这些序列依次连接,变为一维的线性长轴。
比如物种基因组存在chr1,chr2,chr3三条染色体序列,其中chr1长度为x1,chr2长度为x2,chr3长度为x3等,本发明会制作一个长度为x1+x2+x3的长轴,chr1上的坐标对应长轴的(1,x1)坐标范围,chr2上的坐标对应长轴的(x1+1,x1+x2)范围,chr3对应长轴的(x1+x2+1,x1+x2+x3)范围,这样基因组上的所有位置都可以用一个连续且唯一的数字表示。
聚类法
在本发明中,在所述远距离比对读长对位置的聚类步骤中,即在不正常比对读长对的局部聚类步骤中,有两种方法可独立或组合使用:排序聚类法和凝聚层次聚类法。
在本发明方法中,通过聚类处理可以充分利用不正常比对读长对的比对方式和比对位置的不同方面的信息,使得后续在不正常比对簇寻找断点位置变得更准确。
排序聚类法
在本发明中,所述排序聚类法是一种无监督的分类方法,在所述远距离比对读长对位置的聚类步骤中,将上述搜集的四组不正常比对读长对在组范围内对每一个组内的读长对两端比对位置构成的(x,y)二维坐标点集进行排序层次聚类,将点集里的读长对聚为一个读长对簇,簇内元素之间彼此可以认为是一簇支持一个潜在断点对,即同一个基因组嵌合事件的证据集合。
排序聚类法的计算速度较凝聚层次聚类方法快,算法复杂度为O(n),适用于检测基因组复杂度较高的样本,但没有凝聚性层次性聚类方法分类效果好。
聚层次聚类法
所述聚层次聚类法是一种无监督的分类方法,在所述远距离比对读长对位置的二维凝聚层次聚类步骤中,将搜集的四组不正常比对读长对在组范围内对每一个组内的读长对两端比对位置构成的(x,y)二维坐标点集进行凝聚层次聚类,将点集里的读长对聚为一个读长对簇,簇内元素之间彼此可以认为是一簇支持一个潜在断点对,即同一个基因组嵌合事件的证据集合
聚层次聚类法的聚类效果好,很少出现有交叉区域的情况。算法复杂度为O(n^2),在待测样本基因组复杂度水平不高的情况下会有很好的表现。
断裂位置对的检测
以图5为例,图5中的红色记录代表的潜在断点位置是它的匹配起始位点加上它左边匹配部分的长度,即2号染色体的42493957(42493874+83,匹配起始位点为42493874,匹配的部分长度为83,这个坐标恰好是67S未匹配的部分第一个位置的坐标)。它的另一部分的比对方式记录在它的SA tag里:chr2,29448026,-,67M83S,60,0,可以算出,配对断点位置应该在2号染色体上的29448093(29448026+67)位置上。
嵌合事件的鉴别
在本发明方法中,还可包括对嵌合事件的鉴别,该步骤通常包括:
结合簇内的不正常比对读长对的比对方向组合信息与所述的断点位置信息组合,推断出考察的当前嵌合事件的类型,其中断定标准如下:
a.远距离比对读长对两端读长不在一条染色体上,则鉴定融合事件的融合方式为易位(Translocation)。
b.不正常比对读长对两端在一条染色体上,其两端方向相同,且经过断点的splitreads两个部分的比对方向相反,则鉴定此融合事件的融合方式为颠换(Inversion)。
c.除去A和B中所述的两种类型之外,称为未知类型(Unknown)。
检测装置
本发明还提供了一种检测基因组嵌合片段的检测装置,该检测装置用于执行本发明第一方面中所述的检测方法。
典型地,所述检测装置包括包括以下主要模块:
(D1)序列比对模块,所述模块用于将测序的短序列比对到一预定的基因组;
(D2)不正常比对读长对的聚类模块,所述模块用于对不正常比对读长对(readpair)进行收集和进行聚类分析;
(D3)断点对的检测模块,所述模块用于在不正常比对读长对簇覆盖基因组区域内,对断点对进行检测,从而获得位于所述的被覆盖的基因组区域内的唯一的融合事件断点位置对;
(D4)注释-质控模块,所述模块对断点位置对所涉及的融合事件进行注释和进行质控处理,其中如果满足质控标准则提示所述断点位置所涉及的融合事件是真的,而不满足质控标准则提示所述断点位置所涉及的融合事件是非真的;和
(D5)输出模块,所述模块用于输出检测结果,所述检测结果包括最终真融合事件相关信息。
在本发明中,模块代码可采用常规的编程语言进行编写,其中包括(但并不限于):C/C++,R,Python。
本发明的主要优点包括:
(a)高稳定性;
(b)高灵敏度;
(c)高特异性;
(d)计算资源消耗低和计算速度快;
(e)充分利用了双端测序数据的信息只需一次比对,不需组装,不需重比对;
(f)确定的融合区域位置分界清晰,断点位置准确。
下面结合具体实施例,进一步阐述本发明。应理解,这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明具体条件的实验方法,通常按照常规条件,或按照制造厂商所建议的条件。除非另外说明,否则百分比和份数按重量计算。
实施例1
本实施例采用本发明所述的检测方法和装置,对来自细胞系的样本杰西验证。
在本实施例中,测序的读长来自表1和表2中所述的细胞系,这些细胞系含有一些基因融合。
表1
细胞系 基因融合 BreakPoint_Pair1 BreakPoint_Pair2
H2228 ALK-EML4 chr2:29448091 chr2:42493956
H2228 ALK-EML4 chr2:30003907 chr2:42508120
RT4 FGFR3-TACC3 chr4:1729558 chr4:1808765
HCC78 SLC34A2-ROS1 chr4:25666628 chr6:117658323
LC2AD RET-CCDC6 chr10:43609946 chr10:61638609
HD753 FGFR3-TACC3 chr4:1729558 chr4:1808765
HD753 SLC34A2-ROS1 chr4:25666628 chr6:117658323
表2
Sample Name FUSM1 FUSM2 FUSM3 FUSM4 FUSM5
NA18536 5% 10% 20% 30% 35%
H2228 10% 20% 30% 35% 5%
RT4 20% 30% 35% 5% 10%
HCC78 30% 35% 5% 10% 20%
LC2AD 35% 5% 10% 20% 30%
如表1和表2所示,本实施例采用H2228,RT4,HCC78,LC2AD,NA18536(表1)以及前面5种细胞系按照不同比例混合而形成的混合样本FUSM1,FUSM2,FUSM3,FUSM4,FUSM5,还有HD753细胞系,用这些样本的测序数据来检测装置的检测效果。
从表1中可以看出,H2228细胞系拥有已知的2个ALK-EML4融合事件,RT4拥有1个FGFR3-TACC3融合事件,HCC78拥有SLC34A2-ROS11个融合事件,LC2AD拥有1个RET-CCDC6融合事件,HD753具有FGFR3-TACC3和SLC34A2-ROS1两个融合事件,NA18536为阴性样本。
混合样本中FUSM2中包含5%的LC2AD样本,则可以知道LC2AD所携带的RET-CCDC6融合事件的频率在5%以下。以此类推,FUSM3含有频率在5%以下的SLC34A2-ROS1融合事件,FUSM4含有频率在5%以下的FGFR3-TACC3融合事件,FUSM5含有频率在5%以下的两个ALK-EML4融合事件,在5%敏感性容忍度下装置检测不到可以不计为假阴性并不参与敏感性计算。因此在最后计算本发明装置的敏感度和特异性数值是会计算出两种情况下的结果:真实情况和小于5%融合变异频率容忍度,如表3中所示“sensitivity”和“PPV”字段表示的是真实情况下的灵敏度和特异性,“sensitivity(5%)”和“PPV(5%)”字段表示的是小于5%融合变异频率容忍度情况下的灵敏度和特异性。
在上述短序列比对到基因组的模块中,本发明下面阐述的优选实施例中,所有样本的物种都为人,都采用目标捕获测序手段建库进行双端测序,使用的测序平台为NextSeq。得到的DNA测序数据读长均为150bp,质控结果符合要求的数据量为99.9%,平均覆盖度均为1700X,设计的panel完全覆盖已知的几个融合区域。
所有样本都是采用bwa 0.7.16版本,采用mem算法,开启“use soft clipping forsupplementary alignments”开关(-Y)和“mark shorter split hits as secondary”开关(-M,这样才可以在分两部分比对到基因组不同区域的读长的比对记录里出现携带另一部分比对位置和方式信息的SA标签),基因组版本为人类基因组hg19版本,其余bwa mem算法的参数为默认。这样产生的比对结果文件里的记录才会有足够下游装置要使用好的信息。
本发明的上述检测方法通过对现有的对现有的基因融合的检测方法进行改造。在上述不正常比对读长对的收集和分组模块,
在上述不正常比对读长对局部聚类模块之前,识别不正常比对读长对集合之后,本装置会按照不正常比对两端的读长比对的方向异同,是否在同一染色体上将不正常比对读长对分为四类:两端异向且比对到不同染色体上、两端异向且比对到相同染色体上、两端同向且比对到不同染色体上、两端同向且比对到相同染色体上。在这四类产生的四个组的各自范围内,再依据组内不正常比对读长对的基因组上比对位置对它们进行无监督聚类。这样就可以更精确的得到潜在断点可能存在的簇,不正常比对读长对的组内类型一致也为后来鉴定融合类型提供可靠依据。
本发明在识别不正常比对读长对集合且进行分组之后,将不正常比对读长对集合里的每端读长进行一维转化。在现有条件下,读长对两端读长的比对位置都是按照“染色体编号:比对起始点在该染色体上所处的坐标数值”这种方式记录的。这种格式不方便进行下一步的不正常比对读长对位置聚类步骤。于是,本发明中将所述物种的序列条目按照编号排序,并把这些序列依次连接,变为一维的线性长轴。比如物种基因组存在chr1,chr2,chr3三条染色体序列,其中chr1长度为x1,chr2长度为x2,chr3长度为x3等,本发明会制作一个长度为x1+x2+x3的长轴,chr1上的坐标对应长轴的(1,x1)坐标范围,chr2上的坐标对应长轴的(x1+1,x1+x2)范围,chr3对应长轴的(x1+x2+1,x1+x2+x3)范围,这样基因组上的所有位置都可以用一个连续且唯一的数字表示。
本发明在识别不正常比对读长对集合且进行分组之后改进了不正常比对读长对的局部聚类步骤,且聚类方法有两种可供用户选择和使用:排序聚类和凝聚层次聚类法。此聚类步骤可以充分利用不正常比对读长对的比对方式和比对位置的不同方面的信息,使得后续在不正常比对簇寻找断点位置变得更准确。在上述中提及的排序聚类法的步骤中,在不正常比对读长对收集步骤中产生的四组不正常比对读长对在组范围内对每一个组内的读长对两端比对位置构成的(x,y)二维坐标点集进行排序层次聚类,将点集里的读长对聚为一个读长对簇,簇内元素之间彼此可以认为是一簇支持一个潜在断点对,即同一个基因组嵌合事件的证据集合。具体步骤为:
1)将不正常比对读长对两端读长的比对位置二维坐标(x,y)集合内的元素按照x坐标的顺序从小到大排序;
2)定义上步排序后第一个元素的x坐标为previous position,并记录它为当前cluster的第一个元素;
3)从第二个元素开始遍历排序后的集合,比较当前元素的x坐标值是否与previous position相差设定阈值以内的数值,是的话则将此元素记录为当前cluster的一个新元素,并接着遍历下一个元素,对下一个元素重复此步骤,否的话则进入下一步;
4)判断当前cluster内元素个数是否大于等于已设定好的阈值(一般为2),如果是,则保留此cluster以及此cluster内的元素,并将cluster id赋给其内的每个元素,如果否,则丢弃此cluster以及此cluster内的元素。判断结束后,清空当前cluster,更新cluster id,将当前元素的X坐标设置为previous position;
5)遍历完一遍组内元素后,再将剩下的元素集合按照点的y坐标从小到大排列,重复3)和4),确定最后留下的参与成簇的不正常比对读长对元素集合以及它们所归属的簇id。
这种方法计算速度较凝聚层次聚类方法快,算法复杂度为O(n),适用于检测基因组复杂度较高的样本,但没有凝聚性层次性聚类方法分类效果好。
另外,在上述提及的凝聚型层次聚类算法中,搜集的四组不正常比对读长对在组范围内对每一个组内的读长对两端比对位置构成的(x,y)二维坐标点集进行凝聚层次聚类,将点集里的读长对聚为一个读长对簇,簇内元素之间彼此可以认为是一簇支持一个潜在断点对,即同一个基因组嵌合事件的证据集合。具体步骤为:
1)将原始点集内的每个点初始化为单点簇;
2)遍历簇集,每次合并两个最近的簇为一个簇,这样离得相对较近的点都会逐渐的纳入到同一个簇里。其中簇元素之间的距离有三种:single-linkage,complete-linkage,average-linkage:single-linkage为两个簇内对象之间的最小距离;complete-linkage两个簇内对象之间的最大距离;complete-linkage两个簇内对象之间的平均距离。
3)重复2)步骤,直到当前簇集内所有元素簇的簇间距离都大于设定的距离值,停止合并点簇步骤。
4)遍历当前合并停止后生成的簇集,过滤簇元素内点元素个数小于设定值的簇元素,将每个簇赋予一个簇元素ID;
5)遍历原始点集,记录当前点所属的簇元素ID;这样就得到一个符合需求的簇集合。
此种聚类方法聚类效果好,很少出现有交叉区域的情况。算法复杂度为O(n^2),在待测样本基因组复杂度水平不高的情况下会有很好的表现。
本次实施例中使用的是凝聚型层次聚类算法对不正常比对读长对集合进行聚类,在此步骤中,统一选用single-linkage距离作为簇元素之间距离的衡量距离,设定的簇间距离阈值统一为600bp,保留簇元素个数大于2的不正常比对读长对聚类簇,留作下一步寻找断点位置的候选簇。
随后进行的是远距离比对读长对簇覆盖基因组区域内断裂位置对的检测模块中,此模块具体分为下面几个步骤:
1)遍历簇集,以每个簇的中心位置的x和y坐标分别向左右扩展一定距离,形成一个矩形范围。
2)遍历上步扩展的成对区域内所有的比对记录,找到这些比对记录中splitreads的比对记录(即一部分比对到基因组上,另一部分也比对到基因组上的的读长)。
3)分析上步获取的split reads比对记录,分析这些比对记录中读长两部分的比对模式(通过比较Cigar和SA tag的值),潜在断点就在两条记录中所描述的读长中部match和soft clipping的交界处,将读长的这个交界处对应的参考基因组序列上对应的坐标计算出来,比如附图5中的红色记录它代表的潜在断点位置是它的匹配起始位点加上它左边匹配部分的长度,即2号染色体的42493957(42493874+83,匹配起始位点为42493874,匹配的部分长度为83,这个坐标恰好是67S未匹配的部分第一个位置的坐标)。它的另一部分的比对方式记录在它的SA tag里:chr2,29448026,-,67M83S,60,0,可以算的出来,配对断点位置应该在2号染色体上的29448093(29448026+67)位置上。
4)遍历成对区域中的所有的split reads记录,计算每一条记录所对应的断点位置,再按照读长名字将在两端都出现的读长挑出来组合在一起,组成一个断点对。将断点对按照match部分的方向组合分为四类:(左,左),(左,右),(右,左)(右,右),同类内的断点对集合进行频次排序,挑选频次最高的那个断点对组合,最后在四类的频次最高的断点对组合中选取频次最高的断点对组合,作为最后这个不正常比对读长簇所覆盖区域中得到的唯一的融合事件断点位置对。
5)记录支持此断点对所代表的融合事件的支持不正常比对读长对的个数和splitreads个数,这部分信息会分别记录在结果文件的N_Discordant_Pairs字段和N_Split_Reads中。若支持的split reads个数不大于设定好的个数阈值,则不会被输出。
在此模块中,本实施例设置的扩展距离为600bp,设定好的split reads个数阈值为5。在远距离比对读长对簇覆盖基因组区域内断裂位置对的检测步骤后进一步分析断点左右的读长比对模式,鉴别此嵌合事件的融合类型。其步骤包括:
结合簇内的不正常比对读长对的比对方向组合信息与断点位置信息组合,推断出考察的当前嵌合事件的类型,断定标准如下:
a.远距离比对读长对两端读长不在一条染色体上,则鉴定融合事件的融合方式为易位(Translocation)。
b.不正常比对读长对两端在一条染色体上,其两端方向相同,且经过断点的splitreads两个部分的比对方向相反,则鉴定此融合事件的融合方式为颠换(Inversion)。
c.除去A和B中所述的两种类型之外,称为未知类型(Unknown)。
本实施例中严格按照此规则进行操作。
随后是对断裂位置对所涉及的融合事件的注释模块,这一模块包含以下步骤:
a)注释融合事件所涉及的两个断点所出的基因区域,外显子位置
b)注释融合事件所涉及的两个断点附近是否为重复序列区域,如果是的话,则判定此融合事件为假阳性,在最终结果中过滤此事件,不报告出来。
c)如果测序样本物种为人,则会进而以当前融合事件的两个断点所处的基因融合对搜索该装置内置的热点人类融合基因对数据库(该数据库定期更新),如果有匹配结果则加上hotspot的标签。
最后是最终真融合事件相关信息的输出模块,此模块将注释模块中所得到的可能的符合质控标准的融合事件相关的下列信息作为表格条目的方式写到结果文件中:
a)两个断点的坐标位置。
b)两个断点各自所处的基因及外显子位置注释(基因名以及外显子编号信息)。
c)此融合事件的类型,包含颠换,易位和未知类型(inversion,translocation,unknown)。
d)此融合事件的支持远距离比对读长对个数,断裂读长数目,此融合事件的等位基因频率。
e)融合事件所涉及的融合基因对与软件内置的热点融合基因对数据库之间的比对结果,如果有match,则有hotspot标签。此部分信息仅在输入样本数据的物种为人类才有意义。
表3HCC78的结果文件内容
在本实施例中的HCC78样本full coverage的结果文件如表3。各字段的意义分别是:
a.Hotspot_Pair_Match:1代表此融合基因事件在人类热点融合基因对数据库里有匹配。(0代表无匹配)
b.BreakPoint_Pair1/BreakPoint_Pair2:融合断点1在基因组上的位置为chr4:25666630,融合断点2在基因组上的位置为chr6:117658326.
c.Behalf_Gene_Pair1/Behalf_Gene_Pair2:融合断点1在基因组上的位置所在的基因为SLC34A2,融合断点1在基因组上的位置所在的基因为ROS1.
d.BreakPoint_Info_Pair1/BreakPoint_Info_Pair2:这两个字段表示的是两个断点所处基因组位置上的外显子注释,其中融合断点1的注释内容为“+:NM_001177999.1:3-4”,即基因处在正链上,选中的转录本为NM_001177999.1,断点处在此转录本的3号和4号外显子之间的内含子;融合断点2的注释内容为“-:NM_002944.2:31-32”,即基因处在负链上,选中的转录本为NM_002944.2,断点处在此转录本的31号和32号外显子之间的内含子上。
e.N_Disconcordant_Pairs:支持此融合事件的不正常比对读长对个数为171。
f.N_Split_Reads:支持此融合事件的断裂读长个数为211.
g.P1_BreakPoint_Depth/P2_BreakPoint_Depth:融合断点位置上的测序深度。融合断点1为299,融合断点2为798.
h.P1_Alle_Freq/P2_Alle_Freq:两融合断点上分别算出来的融合等位基因频率。融合断点1中计算出来的融合事件频率为0.7057,融合断点2中计算出来的融合事件频率为0.264411。
i.Fusion_Type:此融合事件的融合类型为易位(Translocation)。
表4灵敏度和特异性结果
Depth TP FN TP(5%) FN(5%) FP Ignore sensitivity sensitivity(5%) PPV PPV(5%)
250 26 6 23 4 0 5 81.25% 85.19% 100.00% 100.00%
500 28 4 25 2 1 5 87.50% 92.59% 96.55% 96.15%
700 30 2 26 1 2 5 93.75% 96.30% 93.75% 92.86%
900 30 2 26 1 3 5 93.75% 96.30% 90.91% 89.66%
1000 31 1 27 0 4 5 96.88% 100.00% 88.57% 87.10%
1200 31 1 27 0 5 5 96.88% 100.00% 86.11% 84.38%
full 31 1 27 0 6 5 96.88% 100.00% 83.78% 81.82%
最后,按照图3所描述的方案操作,这10个样本经过核酸提取,捕获建库,产出的测序数据比对到基因组上,再对比对文件进行随机取样至不同深度:250X,500X,700X,900X,1000X,1200X.这样就可以得到上述各测序深度下本装置的检测灵敏度和特异性。得到表4的统计数字,可以看到在500X的测序深度下,本装置可达到87.50%的灵敏度,96.55%的特异性,在5%的稀释度容忍情况下,灵敏度和特异性可达到92.59%和96.15%;在1000X的测序深度下,本装置可达到96.88%的灵敏度,88.57%的特异性,在5%的稀释度容忍情况下,灵敏度和特异性可达到100.00%和87.10%。数据表明,本发明方法和装置的检测水平很高。
上述结果说明,本发明方法充分利用了双端测序数据的信息,只需一次比对,不需组装,不需重比对,确定的融合区域位置分界清晰,断点位置准确,结果稳定,灵敏度,特异性高。
对比例1
Factera软件是本领域常用的适用于检测基因组结构变异的装置。FACTERA使用BAM文件作为输入,产生bam文件的短序列比对软件必须可以对读长的未比对部分进行“soft clipping”标记。其他还需要输入你要考察的可能发生片段结构变异的潜在基因组区域的bed文件,2BIT格式的人类基因组序列文件。FACTERA的工作流程分为三个阶段:i.寻找DRP读长对簇;ii.单核苷酸分辨率的断点识别;iii.潜在融合事件的计算机模拟验证。
第一阶段,DRP识别以后,根据其两端读长距离最近的外显子待判定两端的读长分别属于哪个基因的范围内,从而将每个DRP归到一个特定的geneA-geneB基因对DRPs组。然后针对每个geneA-geneB DRPs组,定义该组两端的基因组考察区域Ri,使用的是对应两端的DRP读长或者读长最近外显子所处的位置的最3'端的位置和最5'端位置作为两端区域的最大值和最小值。每个Ri必须要有两个以上的DRP做为支持证据,否则进入不了下一个阶段。
第二阶段,FACTERA选出支持DRP数目排前n个(默认n=5)的Ri作为待考察基因组区域对。对于每个被选中的Ri,FACTERA将Ri的两个区域部分上的某些SR,筛选出来,筛选标准如下:a.SR读长全长的中间部位附近有一个断点;b.SR读长的被剪切掉的部分的长度大于15bp(默认),以尽可能去除非特异比对的序列。Ri两端的某两个SR,R1和R2,如果是来自于同一个融合片段,那么R1的未比对的部分应该和R2的比对部分在序列上匹配,反之亦然。判断这一点的部分FACTERA使用快速k-mer索引和比较算法。这个算法首先使用滑窗法将R1的比对上的序列部分分解为所有可能的长度为k(默认k=10)的短序列,又叫k-mer,然后将每个k-mer存储在一个哈希表中,并随之建立与此表相匹配的最小序列的索引。随后将R2的未比对的序列部分迭代的分解为k-mer组合,搜索前面R1的k-mer的哈希表,寻找匹配的序列。如果达到设置的匹配阈值,则判定R1,R2为支持一个融合事件的SR证据。两个SR的断点位置即为支持的融合事件断点对。
第三阶段,FACTERA模拟验证上步识别的融合事件。FACTERA将所有SR读长的未比对部分序列和未比对上的读长使用blastn与此区域(断点附近左右500bp的区域)的融合序列进行比对。比对结果中比对结果的一致性(identity)大于等于95%,且比对长度大于等于90%的读长长度的读长才会被保留。默认FACTERA会报告出支持的spanning-reads(经过断点的读长)数目大于5的融合事件。
为了比较Factera和本文中介绍的方法的效果差别,使用相同的人基因组目标区域捕获DNA二代测序数据测试Factera软件,按照图4所描述的方案操作,这10个样本经过核酸提取,捕获建库,产出的测序数据比对到基因组上,再对比对文件进行随机取样至不同深度:250X,500X,700X,900X,1000X,1200X.这样就可以得到上述各测序深度下Factera的检测灵敏度和特异性。得到表5的统计数字,可以看到在500X的测序深度下,Factera可达到68.75%的灵敏度,84.62%的特异性,远低于本装置同等条件下的87.50%的灵敏度,96.55%的特异性,在5%的稀释度容忍情况下,Factera灵敏度和特异性分别为74.07%和83.33%,远低于本装置同等条件下的92.59%的灵敏度和96.15%的特异性;在1000X的测序深度下,Factera的灵敏度为68.75%,特异性为73.33%,远低于本装置的96.88%的灵敏度,88.57%的特异性,在5%的稀释度容忍情况下,Factera的灵敏度也只能到70.37%,特异性只到70.37%,远低于本装置的100.00%和87.10%。
表5
从以上的描述中,可以看出,与对比例相比,本发明装置的实施例通过分析样本测序数据比对到基因组序列上的比对文件中带有SA标签的记录的比对模式以得到断点的具体断裂位置,提高了本发明装置的灵敏度(即降低了假阴性比率)和特异性(即降低了假阳性比例);并且较对比例装置来说,这一个策略使得本装置对外部装置环境的依赖性远远低于对比例中的Factera(此装置还依赖blastn软件,twoBitToFa软件),本装置不需要依赖外部任何软件;本发明在利用不正常比对读长对来确认候选短点所处区域的步骤前,对不正常比对读长对集合进行聚类,可以排除很多噪音,富集融合信号,这一步也是本装置高灵敏度和特异性的关键;并且在聚类这一步提供具有优秀聚类效果的凝聚层次聚类方法和快速准确的排序聚类法供用户根据样本基因组复杂程度的高低来选择,提高了装置的稳定性。本发明检测方法相比现有的融合检测方法具有灵敏度高,特异性高,检测周期短,运算资源消耗少等不可替代的优势。
显然,本领域的技术人员应该了解,本发明装置上述的一些模块或一些步骤可用通用的计算装置来实现,它们可以集中在单个的计算装置上,用计算装置可执行的程序代码来实现,从而可以将它们存储在存储装置中有机酸装置来执行。本发明不限制于任何特定的硬件和软件结合。
在本发明提及的所有文献都在本申请中引用作为参考,就如同每一篇文献被单独引用作为参考那样。此外应理解,在阅读了本发明的上述讲授内容之后,本领域技术人员可以对本发明作各种改动或修改,这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims (10)

1.一种检测基因组嵌合片段的检测方法,其特征在于,所述检测方法包括:
(M1)将测序获得的短序列与基因组进行比对,从而获得正常比对的读长对和不正常比对的读长对;
(M2)对不正常比对读长对(read pair)进行收集,并对所述不正常比对读长对进行聚类,从而获得不正常比对读长对簇;
(M3)在不正常比对读长对簇所覆盖的基因组区域内,对断点对进行检测,从而获得位于所述的被覆盖的基因组区域内的唯一的融合事件断点位置对;
(M4)对所述断点位置对所涉及的融合事件的进行注释和质控检测,其中如果满足质控标准则提示所述断点位置所涉及的融合事件是真的,而不满足质控标准则提示所述断点位置所涉及的融合事件是非真的;和
(M5)输出上一步骤的检测结果。
2.如权利要求1所述的检测方法,其特征在于,在步骤(M1)中,在将所述短序列比对到基因组,从而获取测序读长比对信息文件;
优选地,在所述的比对过程所得到的比对记录文件中,将所涉及的读长分为5类:
1)全部比对到基因组上的读长,视作A类读长;
2)一部分比对到基因组上,但另一部分未比对到基因组上的的读长,视作B类读长;
3)一部分比对到基因组上,另一部分也比对到基因组上的读长,视作C类读长;其中较长的部分称为主要比对部分,视作C1类,较短的部分称为次级比对部分,视作C2类;
4)读长对两端都比对到基因组上,读长对两端比对到相隔较远的位置上,或者读长对两端比对方向相同(正常情况读长对两端比对方向应相反)但读长对两端读长都不属于上述C类,视作D类,也叫不正常比对读长对(Discordant read pair);
5)不符合以上A,B,C,D类任何一类的读长视作E类。
3.如权利要求1所述的检测方法,其特征在于,在步骤(M2)中,即在不正常比对读长对的收集中,将比对短序列读长到基因组得到的所有远距离比对到基因组的D类读长对挑选出来。
4.如权利要求1所述的检测方法,其特征在于,在所述的聚类步骤之前,对不正常比对读长对的两端读长在基因组上坐标进行一维转换。
5.如权利要求1所述的检测方法,其特征在于,在所述的聚类步骤中,用选自下组的聚类方法进行聚类:排序聚类法和凝聚层次聚类法。
6.如权利要求5所述的检测方法,其特征在于,所述排序聚类法包括以下步骤:
1)将不正常比对读长对两端读长的比对位置二维坐标(x,y)集合内的元素按照x坐标的顺序从小到大排序;
2)定义上步排序后第一个元素的x坐标为previous position,并记录它为当前cluster的第一个元素;
3)从第二个元素开始遍历排序后的集合,比较当前元素的x坐标值是否与previousposition相差设定阈值以内的数值,是的话则将此元素记录为当前cluster的一个新元素,并接着遍历下一个元素,对下一个元素重复此步骤,否的话则进入下一步;
4)判断当前cluster内元素个数是否大于等于已设定好的阈值(一般为2),如果是,则保留此cluster以及此cluster内的元素,并将cluster id赋给其内的每个元素,如果否,则丢弃此cluster以及此cluster内的元素;判断结束后,清空当前cluster,更新cluster id,将当前元素的X坐标设置为previous position;
5)遍历完一遍组内元素后,再将剩下的元素集合按照点的y坐标从小到大排列,重复3)和4),确定最后留下的参与成簇的不正常比对读长对元素集合以及它们所归属的簇id。
7.如权利要求5所述的检测方法,其特征在于,所述凝聚层次聚类法包括以下步骤:
1)将原始点集内的每个点初始化为单点簇;
2)遍历簇集,每次合并两个最近的簇为一个簇,这样离得相对较近的点都会逐渐的纳入到同一个簇里;其中簇元素之间的距离有三种:single-linkage,complete-linkage,average-linkage;其中,single-linkage为两个簇内对象之间的最小距离;complete-linkage为两个簇内对象之间的最大距离;complete-linkage为两个簇内对象之间的平均距离;
3)重复2)步骤,直到当前簇集内所有元素簇的簇间距离都大于设定的距离值,停止合并点簇步骤;
4)遍历当前合并停止后生成的簇集,过滤簇元素内点元素个数小于设定值的簇元素,将每个簇赋予一个簇元素ID;
5)遍历原始点集,记录当前点所属的簇元素ID;这样就得到一个符合需求的簇集合。
8.如权利要求1所述的检测方法,其特征在于,在所述远距离比对读长对簇覆盖基因组区域内断裂位置对的检测步骤中,此模块具体分为下面几个步骤:
a.遍历簇集,以每个簇的中心位置的x和y坐标分别向左右扩展一定距离,形成一个矩形范围;
b.遍历上步扩展的成对区域内所有的比对记录,找到这些比对记录中split reads的比对记录,其中split reads是一部分比对到基因组上,另一部分也比对到基因组上的的读长;
c.分析上步获取的split reads比对记录,分析这些比对记录中读长两部分的比对模式(通过比较Cigar和SA tag的值),潜在断点就在两条记录中所描述的读长中部match和soft clipping的交界处,将读长的这个交界处对应的参考基因组序列上对应的坐标计算出来;
d.遍历成对区域中的所有的split reads记录,计算每一条记录所对应的断点位置,再按照读长名字将在两端都出现的读长挑出来组合在一起,组成一个断点对;其中,将断点对按照match部分的方向组合分为四类:(左,左),(左,右),(右,左)(右,右),同类内的断点对集合进行频次排序,挑选频次最高的那个断点对组合,最后在四类的频次最高的断点对组合中选取频次最高的断点对组合,作为最后这个不正常比对读长簇所覆盖区域中得到的唯一的融合事件断点位置对;
e)记录支持此断点对所代表的融合事件的支持不正常比对读长对的个数和splitreads个数,这部分信息会分别记录在结果文件的N_Discordant_Pairs字段和N_Split_Reads中;其中若支持的split reads个数不大于设定好的个数阈值,则不会被输出。
9.如权利要求1所述的检测方法,其特征在于,所述的注释包含对选自下组的一个或多个或全部信息进行注释:
a)注释融合事件所涉及的两个断点所出的基因区域,外显子位置;
b)注释融合事件所涉及的两个断点附近是否为重复序列区域,如果是的话,则判定此融合事件为假阳性,在最终结果中过滤此事件,不报告出来;
c)如果测序样本物种为人,则会进而以当前融合事件的两个断点所处的基因融合对搜索热点人类融合基因对数据库,如果有匹配结果则加上hotspot的标签。
10.一种检测基因组嵌合片段的检测装置,其特征在于,所述检测装置包括:
(D1)序列比对模块,所述模块用于将测序的短序列比对到一预定的基因组;
(D2)不正常比对读长对的聚类模块,所述模块用于对不正常比对读长对(read pair)进行收集和进行聚类分析;
(D3)断点对的检测模块,所述模块用于在不正常比对读长对簇覆盖基因组区域内,对断点对进行检测,从而获得位于所述的被覆盖的基因组区域内的唯一的融合事件断点位置对;
(D4)注释-质控模块,所述模块对断点位置对所涉及的融合事件进行注释和进行质控处理,其中如果满足质控标准则提示所述断点位置所涉及的融合事件是真的,而不满足质控标准则提示所述断点位置所涉及的融合事件是非真的;和
(D5)输出模块,所述模块用于输出检测结果,所述检测结果包括最终真融合事件相关信息。
CN201810570027.2A 2018-06-05 2018-06-05 用于检测癌症样本基因融合的检测方法和装置 Active CN108830044B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810570027.2A CN108830044B (zh) 2018-06-05 2018-06-05 用于检测癌症样本基因融合的检测方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810570027.2A CN108830044B (zh) 2018-06-05 2018-06-05 用于检测癌症样本基因融合的检测方法和装置

Publications (2)

Publication Number Publication Date
CN108830044A true CN108830044A (zh) 2018-11-16
CN108830044B CN108830044B (zh) 2020-06-26

Family

ID=64143927

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810570027.2A Active CN108830044B (zh) 2018-06-05 2018-06-05 用于检测癌症样本基因融合的检测方法和装置

Country Status (1)

Country Link
CN (1) CN108830044B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109712672A (zh) * 2018-12-29 2019-05-03 北京优迅医学检验实验室有限公司 检测基因重排的方法、装置、存储介质及处理器
CN109979528A (zh) * 2019-03-28 2019-07-05 广州基迪奥生物科技有限公司 一种单细胞免疫组库测序数据的分析方法
CN110010193A (zh) * 2019-05-06 2019-07-12 西安交通大学 一种基于混合策略的复杂结构变异检测方法
CN110322925A (zh) * 2019-07-18 2019-10-11 杭州纽安津生物科技有限公司 一种预测融合基因产生新生抗原的方法
CN111180013A (zh) * 2019-12-23 2020-05-19 北京橡鑫生物科技有限公司 检测血液病融合基因的装置
CN111292809A (zh) * 2020-01-20 2020-06-16 至本医疗科技(上海)有限公司 用于检测rna水平基因融合的方法、电子设备和计算机存储介质
CN112349346A (zh) * 2020-10-27 2021-02-09 广州燃石医学检验所有限公司 检测基因组区域中的结构变异的方法
CN113035273A (zh) * 2021-03-11 2021-06-25 南京先声医学检验有限公司 一种快速、超高灵敏度的dna融合基因检测方法
CN114464252A (zh) * 2022-01-26 2022-05-10 深圳吉因加医学检验实验室 一种检测结构变异的方法及装置

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101561878A (zh) * 2009-05-31 2009-10-21 河海大学 基于改进cure聚类算法的无监督异常检测方法和系统
CN104298892A (zh) * 2014-09-18 2015-01-21 天津诺禾致源生物信息科技有限公司 基因融合的检测装置和方法
CN104302781A (zh) * 2013-05-15 2015-01-21 深圳华大基因科技有限公司 一种检测染色体结构异常的方法及装置
CN104794371A (zh) * 2015-04-29 2015-07-22 深圳华大基因研究院 检测逆转座子插入多态性的方法和装置
CN105989246A (zh) * 2015-01-28 2016-10-05 深圳华大基因研究院 一种基于基因组组装的变异检测方法和装置
CN106909806A (zh) * 2015-12-22 2017-06-30 广州华大基因医学检验所有限公司 定点检测变异的方法和装置
CN106980763A (zh) * 2017-03-30 2017-07-25 大连理工大学 一种基于基因突变频率的癌症驱动基因的筛选方法
CN107267646A (zh) * 2017-08-02 2017-10-20 广东国盛医学科技有限公司 一种基于下一代测序的多基因融合检测方法
CN107992721A (zh) * 2017-11-10 2018-05-04 深圳裕策生物科技有限公司 用于检测目标区域基因融合的方法、装置和存储介质

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101561878A (zh) * 2009-05-31 2009-10-21 河海大学 基于改进cure聚类算法的无监督异常检测方法和系统
CN104302781A (zh) * 2013-05-15 2015-01-21 深圳华大基因科技有限公司 一种检测染色体结构异常的方法及装置
CN104298892A (zh) * 2014-09-18 2015-01-21 天津诺禾致源生物信息科技有限公司 基因融合的检测装置和方法
CN105989246A (zh) * 2015-01-28 2016-10-05 深圳华大基因研究院 一种基于基因组组装的变异检测方法和装置
CN104794371A (zh) * 2015-04-29 2015-07-22 深圳华大基因研究院 检测逆转座子插入多态性的方法和装置
CN106909806A (zh) * 2015-12-22 2017-06-30 广州华大基因医学检验所有限公司 定点检测变异的方法和装置
CN106980763A (zh) * 2017-03-30 2017-07-25 大连理工大学 一种基于基因突变频率的癌症驱动基因的筛选方法
CN107267646A (zh) * 2017-08-02 2017-10-20 广东国盛医学科技有限公司 一种基于下一代测序的多基因融合检测方法
CN107992721A (zh) * 2017-11-10 2018-05-04 深圳裕策生物科技有限公司 用于检测目标区域基因融合的方法、装置和存储介质

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109712672A (zh) * 2018-12-29 2019-05-03 北京优迅医学检验实验室有限公司 检测基因重排的方法、装置、存储介质及处理器
CN109979528A (zh) * 2019-03-28 2019-07-05 广州基迪奥生物科技有限公司 一种单细胞免疫组库测序数据的分析方法
CN109979528B (zh) * 2019-03-28 2021-06-25 广州基迪奥生物科技有限公司 一种单细胞免疫组库测序数据的分析方法
CN110010193A (zh) * 2019-05-06 2019-07-12 西安交通大学 一种基于混合策略的复杂结构变异检测方法
CN110010193B (zh) * 2019-05-06 2021-09-03 西安交通大学 一种基于混合策略的复杂结构变异检测方法
CN110322925A (zh) * 2019-07-18 2019-10-11 杭州纽安津生物科技有限公司 一种预测融合基因产生新生抗原的方法
CN111180013A (zh) * 2019-12-23 2020-05-19 北京橡鑫生物科技有限公司 检测血液病融合基因的装置
CN111180013B (zh) * 2019-12-23 2023-11-03 北京橡鑫生物科技有限公司 检测血液病融合基因的装置
CN111292809A (zh) * 2020-01-20 2020-06-16 至本医疗科技(上海)有限公司 用于检测rna水平基因融合的方法、电子设备和计算机存储介质
CN112349346A (zh) * 2020-10-27 2021-02-09 广州燃石医学检验所有限公司 检测基因组区域中的结构变异的方法
CN113035273A (zh) * 2021-03-11 2021-06-25 南京先声医学检验有限公司 一种快速、超高灵敏度的dna融合基因检测方法
CN114464252A (zh) * 2022-01-26 2022-05-10 深圳吉因加医学检验实验室 一种检测结构变异的方法及装置

Also Published As

Publication number Publication date
CN108830044B (zh) 2020-06-26

Similar Documents

Publication Publication Date Title
CN108830044A (zh) 用于检测癌症样本基因融合的检测方法和装置
Lähnemann et al. Eleven grand challenges in single-cell data science
CN109022553B (zh) 用于肿瘤突变负荷检测的基因芯片及其制备方法和装置
CN103262086B (zh) 识别被测序基因组中的重排
CN109033749A (zh) 一种肿瘤突变负荷检测方法、装置和存储介质
CN106021984A (zh) 一种全外显子组测序数据分析系统
CN106971071A (zh) 一种临床决策支持系统及方法
CN109346130A (zh) 一种直接从全基因组重测序数据中得到微单体型及其分型的方法
CN106446597B (zh) 多物种特征选择及鉴定未知基因的方法
US20020064792A1 (en) Database for storage and analysis of full-length sequences
CN115083521B (zh) 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统
CN113362889A (zh) 基因组结构变异注释方法
CN106460045A (zh) 人类基因组常见拷贝数变异用于癌症易感风险评估
Liu et al. A comparison of topologically associating domain callers based on Hi-C data
CN106021986B (zh) 超低频突变分子一致性序列简并算法
CN109686414A (zh) 仅用于肝癌筛查的特异甲基化检测位点组合的选取方法
CN102831331B (zh) 基于酶切建库双末端测序的长度多态性标记的引物设计开发方法
Bueno-Sancho et al. Field pathogenomics: an advanced tool for wheat rust surveillance
Wu et al. InvBFM: finding genomic inversions from high-throughput sequence data based on feature mining
Valdes et al. Methods to detect transcribed pseudogenes: RNA-Seq discovery allows learning through features
Luebeck et al. AmpliconReconstructor: Integrated analysis of NGS and optical mapping resolves the complex structures of focal amplifications in cancer
CN115762641B (zh) 一种指纹图谱构建方法及系统
Olivieri et al. Iterative variable gene discovery from whole genome sequencing with a bootstrapped multiresolution algorithm
TW201920682A (zh) 多型之檢測方法
Li et al. Detecting genomic deletions from high-throughput sequence data with unsupervised learning

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20200417

Address after: Suzhou City, Jiangsu Province, Suzhou Industrial Park 215000 Xinghu Street No. 218 BioBAY building 201 unit B7

Applicant after: Xukang medical technology (Suzhou) Co., Ltd

Address before: 200120, 4-5 floor, building A, 19 lane, 3399 lane, Kang Xin Road, Pudong New Area, Shanghai.

Applicant before: SHANGHAI JINGZHOU GENE TECHNOLOGY Co.,Ltd.

GR01 Patent grant
GR01 Patent grant