本发明涉及一种基于多重捕获的融合基因生物信息分析方法，属于生物信息。
背景技术：
1、基因融合或染色体重排是癌症中一类重要的体细胞改变，在肿瘤发生的初始阶段起着重要作用[1]。融合基因的类型与肿瘤表型之间的密切联系，使其可作为癌症诊断依据之一[2]。第一例癌症相关的基因重组是在1960年发现的9号和22号染色体重排，该重排会产生一个异常小的染色体，被称为费城染色体[3]。费城染色体会导致9号染色体的原癌基因abl1和22号染色体上的bcr基因组合成新的融合基因。在95％的cml患者中都发现了这个融合基因，它通过bcr启动子，使abl1酪氨酸激酶持续激活。这种不受控的酪氨酸激酶活性是诱发肿瘤的主要原因。第一个针对融合基因设计的靶药——伊马替尼，也在2001被fda批准，用于bcr-abl1阳性的cml患者。很多肿瘤相关的融合基因在其他血液肿瘤中也有发现，例如runx1-runx1t1、pml-rarα等。此外，在实体肿瘤包括肉瘤、上皮肿瘤和中枢神经系统肿瘤中也有癌症相关的融合基因。每种肿瘤类型都有肿瘤特异性的融合，例如尤文肉瘤的ewsr1-fli融合，也有几种癌症常见的融合，如fgfr3-tacc3在脑癌、肺癌和膀胱癌中均有发现[4]。
2、传统细胞遗传学方法染色体显带技术和分子细胞遗传性方法包括荧光原位杂交(fish)、基因融合微列阵、pcr/rt-pcr都被应用于临床检测融合基因。但是这些技术具有较低的分辨率，也无法检测新的基因融合事件。因此，可一次性测百万甚至上十亿序列的二代测序技术(ngs)，开始被用于融合基因的研究。很多研究运用全基因组测序、全外显子测序或者靶向测序等方法，结合不同的生物信息分析流程，发现了大量新的融合基因[5]。当这些技术运用到临床时，考虑到效率和成本问题，针对有临床意义的基因专门设计探针或者引物的靶向测序更具有性价比。靶向测序因为只局限于关注的部分基因组序列，测序量大大降低，且在有限数据量的情况下保持关注区域的高深度，可提高检测敏感性。尽管dna和rna都可用于检测染色体重排事件，但在dna层面的染色体改变并不一定能生成有效的基因融合。除了一些特殊的融合例如igh-myc，一般分析其临床意义时，通常要看是否有融合转录组的表达。因此靶向rna检测方法，在融合基因检测的临床应用上，更具有效率。本发明开发的生物信息方法，针对的就是基于多重捕获法的靶向rna检测方法。多重捕获法需根据目标区域序列设计高特异性的引物，再通过两步pcr扩增来完成目标区域序列的捕获。
3、从大量的ngs序列中找到融合基因需要一系列的生物信息分析处理，包括原始数据清洗、比对、比对后处理、融合基因分析、融合过滤及可视化。融合基因分析相关算法或者软件的选择，是影响结果准确性的关键步骤。很多的融合基因分析工具如starfusion、fusioncacher、arriba等均被广泛用于融合转录本的检测[6]。然而这些软件均是为通用的rna测序技术设计(全转录组、靶向捕获rna测序等)，并没有针对多重捕获方法专门设计算法。直接使用这些检测软件，一方面运行的效率较低，另一方面，因为多重捕获法的局限性，不适用于这些软件流程的过滤方法，它们会将一些真阳性融合认为是低质量融合而过滤。有些rna融合检测流程为降低漏检的概率，经常同时使用多个融合分析软件[7,8]，无疑又增加分析的时间和运算的压力。在此背景下，本发明旨在结合多重捕获方法的特性，开发出新的速度快、敏感性高的融合基因生物信息分析方法。
4、以上内容涉及的参考文献如下：
5、1.taniue k,akimitsu n:fusion genes and rnas in cancerdevelopment.noncoding rna 2021,7(1).
6、2.mertens f,johansson b,fioretos t,mitelman f:the emerging complexityof gene fusions in cancer.nat rev cancer 2015,15(6):371-381.
7、3.nowell pc,hungerford da:chromosome studies on normal and leukemichuman leukocytes.j natl cancer inst 1960,25:85-109.
8、4.parker bc,zhang w:fusion genes in solid tumors:an emerging targetfor cancer diagnosis and treatment.chin j cancer 2013,32(11):594-603.
9、5.j,kumar a,wong sq:overview of fusion detection strategiesusing next-generation sequencing.methods mol biol 2019,1908:125-138.
10、6.haas bj,dobin a,li b,stransky n,pochet n,regev a:accuracyassessment of fusion transcript detection via read-mapping and de novo fusiontranscript assembly-based methods.genome biology 2019,20(1):213.
11、7.arindrarto w,borràs dm,de groen ral,van den berg rr,locher ij,vandiessen s,van der holst r,van der meijden ed,honders mw,de leeuw rh et al:comprehensive diagnostics of acute myeloid leukemia by whole transcriptomerna sequencing.leukemia 2021,35(1):47-61.
12、8.vicente-garcés c,maynou j,fernández g,esperanza-cebollada e,torrebadell m,catalàa,rives s,camós m,vega-garcía n:fusion inpipe,anintegrative pipeline for gene fusion detection from rna-seq data in acutepediatric leukemia.front mol biosci 2023,10:1141310.
技术实现思路
1、本发明的主要目的是：克服现有技术存在的问题，提供一种基于多重捕获的融合基因生物信息分析方法，效率高，速度快，避免过度过滤，敏感性高。
2、本发明解决其技术问题的技术方案如下：
3、一种基于多重捕获的融合基因生物信息分析方法，包括以下步骤：
4、第一步、采用以多重捕获法进行靶向rna测序获得的ngs数据；对这些ngs数据进行预处理，以实现数据清洗；预处理所得ngs数据中含有双端测序数据、以及引物序列。
5、第二步、将预处理所得ngs数据中的双端测序数据与引物序列进行匹配，根据匹配结果，将匹配的双端测序数据与引物序列列入匹配组，将未找到匹配引物的双端测序数据列入unknown组。
6、第三步、将第二步所得的匹配组、unknown组分别与目标基因的参考转录本进行比对；当某对测序数据中的两条序列均与目标基因的参考转录本序列完全匹配时，仅计数但不保存该序列及其匹配结果；当某对测序数据的某条序列与目标基因的参考转录本序列部分匹配时，记录匹配结果，保存不匹配序列；当某对测序数据的两条序列与目标基因的参考转录本序列的比对结果不一致时，保存匹配结果和不匹配序列；根据匹配结果，将不匹配序列作为融合序列并推定其断点坐标，将可能来自同一个融合基因的序列进行组装，并记录支持序列数量。
7、第四步、将所有组装后的融合序列与参考基因组进行比对，分析比对结果，注释融合序列的断点位置，找到另一端融合伙伴基因，计算并输出结果信息；结果信息包括融合基因名称、断点位置、融合序列在两边基因的匹配长度即anchor长度、引物来源、支持序列数量、融合序列占比。
8、第五步、根据第四步得出的结果信息，通过设置过滤条件对所有融合基因进行过滤以移除假阳性融合基因，经排序后输出结果。
9、本发明进一步完善的技术方案如下：
10、优选地，第一步中，预处理的具体过程为：采用ngs数据专用数据质控软件进行处理；所述ngs数据专用数据质控软件为cutadat、fastp、或fastqc。
11、采用以上优选方案，可进一步优化第一步的具体技术特征。
12、优选地，第二步的具体过程包括：
13、s1、利用预处理所得ngs数据中的引物序列，以所有引物序列构建引物序列集合p；将各引物序列按预设步长和预设长度分别拆分为系列短片段序列，记作kmer序列，以所有kmer序列构建序列kmer序列集合k；以kmer序列为键，记录所有来源引物以及将kmer在引物中的起始位置作为值，构建引物索引集合index；以每个引物对应基因、长度、方向以及引物参考转录本中的位置构建信息集合pinfo。
14、s2、利用预处理所得ngs数据中的双端测序数据，双端测序数据为成对的测序数据，以所有成对的测序数据构建集合r；根据引物的理论位置，对每一对测序数据进行引物匹配；如果匹配结果为匹配到符合预设条件的引物，则将该对测序数据与该引物作为一条记录列入匹配组，也即匹配组的各条记录分别由一对测序数据与相应引物构成；如果匹配结果为未找到匹配的引物，则将该对测序数据作为一条记录列入unknown组，也即unknown组的各条记录分别由一对测序数据构成。
15、更优选地，s1中，引物序列的数量为i，引物序列集合p＝{p1,p2,p3...pi}；预设步长记作step，预设长度记作lk且不超过16；kmer序列集合k＝{ki1,ki2,ki3…kij}；来源引物记作kp，起始位置记作kploc，引物索引集合index＝{k1:(kp1，kploc1),k2:(kp2，kploc2),…,kn:(kpn，kplocn)}，其中，(kpn,kplocn)＝{(kp1，kploc1),…,(kpm，kplocm)}、且m≥1；对应基因记作gp，长度记作lp，方向记作d，位置记作ploc，信息集合pinfo＝{kp1:(gp1,lp1,d1,ploc1),…,kpi:(gpi,lpi,di,ploci)}。
16、或者，s2中，集合r＝{r1，r2}，成对的测序数据记作(r1,r2)，且(r1,r2)∈r；引物匹配的具体过程为：将成对的测序数据中与引物对应的序列打断成长度为lk的kmer，记录kmer在该序列中的位置krloc；在index中搜索该kmer，当找到该kmer所对应的引物集合(kpn，kplocn)，对每个可能来源的引物，以其位置为参考，计算引物在该序列中的起始终止位置(ps，pe)：
17、ps＝max(0,krloc-kploc)，
18、pe＝krloc-kploc+lp，
19、取该序列中的[ps:pe]部分与原始引物序列pn进行比对，若错配序列小于预设值，则认为(r1,r2)来源于该引物扩增序列，记录(r1,r2)与引物kp匹配，并列入匹配组；否则，判定(r1,r2)未找到引物来源，列入unknown组。
20、采用以上优选方案，可进一步优化第二步的具体技术特征。
21、优选地，第三步中，具体的比对过程包括：
22、t1、对于匹配组中的每条记录，利用参考smith-waterman动态规划算法的预设算法，将该条记录中的一对测序数据与相应引物所对应基因gp的参考转录本进行比对。
23、t2、对于unknown组的每条记录，将该条记录中的一对测序数据与所有引物所对应的所有基因gp的参考转录本按t1的预设算法分别进行比对；在每一次比对过程中，当一对测序数据中的两条序列均与某条参考转录本序列完全匹配时，跳过后续所有步骤；比对后，若该条记录中的一对测序数据与多个参考转录本部分匹配时，选取最佳结果，并记录疑似断点坐标。
24、t3、按照预设的阈值判断各断点坐标是否相近，将断点坐标相近的测序数据的融合序列进行比较，合并重复序列并计数，然后根据重叠区域，将这些序列组装成重叠群contig，并输出所有组装后的融合序列及其原始来源序列数量。
25、更优选地，t1中，所述预设算法的具体比对过程为：
26、将一对测序数据的两条序列分别作为测序序列a，其长度分别为n；参考转录本序列作为参考序列b，其长度为m；将各测序序列a分别与参考序列b进行比对，同时按下式构建并计算得分矩阵m：
27、
28、
29、i＝(0,1,2,…,n)；j＝(0,1,2,…,m)
30、其中，α、β和δ分别为匹配、不匹配、插入缺失得分；
31、从得分矩阵m中找出最佳得分mopt：
32、mopt＝max(mi,j)
33、记录与此对应的测序序列a和参考序列b的最佳匹配终止位置(iopt，jopt)。
34、上述比对过程可按以下要点进行优化：从引物在参考转录本序列中的起始位置ploc开始向扩增方向搜寻并进行比对；各基因转录本选择mane转录本；当一对测序数据中的两条序列均与参考转录本序列完全匹配即iopt＝n时，跳过后续所有步骤。
35、或者，t2中，疑似断点坐标b1＝(gp，ploc+jopt)。
36、或者，t3中，在比较时，还对融合序列进行过滤，去除低质量序列，低质量序列包括大量重复序列。
37、采用以上优选方案，可进一步优化第三步的具体技术特征。
38、优选地，第四步的具体过程如下：
39、利用star软件将所有组装后的融合序列与参考基因组进行比对后，根据chimeric.out.junction文件提供的断点信息，将所有融合序列的断点位置注释到具体基因和对应外显子，找到另一端融合伙伴基因，并计算融合序列在两边基因的匹配长度即anchor长度；合并重复比对结果和上一步原始序列数量结果，记录并输出以下结果信息：融合基因名称、断点位置、anchor长度、引物来源、支持序列数量、融合序列占比。
40、更优选地，第四步中，注释到具体基因和对应的外显子时利用从现有数据库中下载对应基因组版本的gtf文件；融合序列占比＝融合序列数量÷统一引物来源的序列总数；在注释时，检查断点是否在已知外显子连接处，判断是否为in-frame融合并记录入结果信息。
41、采用以上优选方案，可进一步优化第四步的具体技术特征。
42、优选地，第五步中，具体的过滤类型包括：去除位于非目标染色体的融合基因；去除同源序列间融合的融合基因；若一条融合序列在参考基因组的多个位置均有匹配，则去除按过滤条件判断为比对质量低的融合基因；去除按过滤条件判断为融合序列占比过低的融合基因；移除在多个样本重复出现的融合基因；在过滤条件中设定最小anchor长度，若融合基因的anchor长度低于最小anchor长度，则判断是否属于例外情况，若属于则不予移除，若不属于则将该融合基因移除，其中，例外情况为：融合基因的anchor长度小于或等于引物长度，或者，融合基因的anchor长度小于或等于测序读长减去引物与断点的距离。
43、更优选地，第五步中，在过滤后，先对所得融合基因按预设规则进行优先级排序，再输出结果；其中，预设规则为：优先输出支持序列数量多、融合序列占比高、属于已知与癌症相关的融合基因、属于in-frame融合的融合基因。
44、采用以上优选方案，可进一步优化第五步的具体技术特征。
45、与现有的检测手段相比，本发明方法具有速度快、敏感性高的特点，避免了所有ngs序列在基因组所有序列上的比对，通过提取部分疑似融合序列，并结合去重和组装的手段，极大的提高了比对效率，缩短分析的时间。同时，根据多重捕获的实验特点，重新设计了过滤的条件，避免过度过滤导致去除真阳性融合结果。
1.一种基于多重捕获的融合基因生物信息分析方法，其特征是，包括以下步骤：
2.根据权利要求1所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，第一步中，预处理的具体过程为：采用ngs数据专用数据质控软件进行处理；所述ngs数据专用数据质控软件为cutadat、fastp、或fastqc。
3.根据权利要求1所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，第二步的具体过程包括：
4.根据权利要求3所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，s1中，引物序列的数量为i，引物序列集合p＝{p1,p2,p3...pi}；预设步长记作step，预设长度记作lk且不超过16；kmer序列集合k＝{ki1,ki2,ki3…kij}；来源引物记作kp，起始位置记作kploc，引物索引集合index＝{k1:(kp1，kploc1),k2:(kp2，kploc2),…,kn:(kpn，kplocn)}，其中，(kpn,kplocn)＝{(kp1，kploc1),…,(kpm，kplocm)}、且m≥1；对应基因记作gp，长度记作lp，方向记作d，位置记作ploc，信息集合pinfo＝{kp1:(gp1,lp1,d1,ploc1),…,kpi:(gpi,lpi,di,ploci)}；
5.根据权利要求1所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，第三步中，具体的比对过程包括：
6.根据权利要求5所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，t1中，所述预设算法的具体比对过程为：
7.根据权利要求1所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，第四步的具体过程如下：
8.根据权利要求7所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，第四步中，注释到具体基因和对应的外显子时利用从现有数据库中下载对应基因组版本的gtf文件；融合序列占比＝融合序列数量÷统一引物来源的序列总数；在注释时，检查断点是否在已知外显子连接处，判断是否为in-frame融合并记录入结果信息。
9.根据权利要求1所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，第五步中，具体的过滤类型包括：去除位于非目标染色体的融合基因；去除同源序列间融合的融合基因；若一条融合序列在参考基因组的多个位置均有匹配，则去除按过滤条件判断为比对质量低的融合基因；去除按过滤条件判断为融合序列占比过低的融合基因；移除在多个样本重复出现的融合基因；在过滤条件中设定最小anchor长度，若融合基因的anchor长度低于最小anchor长度，则判断是否属于例外情况，若属于则不予移除，若不属于则将该融合基因移除，其中，例外情况为：融合基因的anchor长度小于或等于引物长度，或者，融合基因的anchor长度小于或等于测序读长减去引物与断点的距离。
10.根据权利要求9所述的一种基于多重捕获的融合基因生物信息分析方法，其特征是，第五步中，在过滤后，先对所得融合基因按预设规则进行优先级排序，再输出结果；其中，预设规则为：优先输出支持序列数量多、融合序列占比高、属于已知与癌症相关的融合基因、属于in-frame融合的融合基因。
