本发明属于生物信息领域，具体涉及一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法。
背景技术：
1、结构变异(structural variation,sv)是基因组变异的一种基本形式，指基因组中大于50bp大小的dna重排事件，包括插入(insertion)、缺失(deletion)、倒位(inversion)、重复(duplication)和易位(translocation)等多种类型。越来越多的证据表明，sv与人类的许多疾病(如认知神经疾病、肥胖症、癌症等)存在紧密关联，并且对生物体的许多其他表型特征(如生长发育、适应性特征等)也有显著影响。因此，全面表征多种形式的sv对于充分了解它们对遗传多样性、物种差异和一些表型的贡献至关重要。
2、深入了解sv的特征、定位和功能，可以更好地理解其在人类健康和疾病中的作用，为疾病诊断、治疗和预防提供重要依据。随着太平洋生物科学(pacbio)和牛津纳米孔技术(oxford nanopore technologies，ont)平台等长读测序技术的快速发展，长读测序数据为更高分辨率和更全面的sv检测提供了机会。目前，主流方法通常根据测序平台的特性和不同sv事件的类型，从比对数据中手动设计启发式规则推断sv。然而，由于sv事件的复杂性，手动设计的启发式规则容易将测序或者比对产生的数据噪声误认为是真实的sv信号。因此，在噪声较大的长读测序数据中准确地检测sv仍然是一个重要挑战。
3、综上所述，现有方法对长读测序数据中sv检测的准确率仍然较低，提出一种新的检测方法以解决上述问题是十分必要的。
技术实现思路
1、本发明的目的是为解决现有方法对长读测序数据中sv检测的准确率低的问题，而提出了一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法。
2、本发明为解决上述技术问题所采取的技术方案是：
3、一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，所述方法具体包括以下步骤：
4、步骤一、将待检测的长读测序数据与参考基因组进行比对，得到cigar字符串以及缺口信息，再根据cigar字符串和缺口信息提取sv信号；
5、并分别记录根据cigar字符串提取的sv信号的特征以及根据缺口信息提取的sv信号的特征；
6、步骤二、利用根据缺口信息提取的sv信号的特征构建原始特征矩阵，对构建的原始特征矩阵进行最小最大值归一化处理，得到处理后的特征矩阵；将处理后的特征矩阵输入随机森林网络，根据随机森林网络的输出来最终确定根据缺口信息提取的各个sv信号是否为真实sv信号以及真实sv信号的类型；
7、步骤三、对步骤一中根据cigar字符串提取的sv信号和步骤二中确定出的真实sv信号进行聚类，得到包含sv信号的各个簇；
8、步骤四、对于步骤三中获得的任意一个簇，若簇内包括步骤二中确定出的真实sv信号，则不对该簇内的sv信号进行处理，直接计算该簇内全部sv信号的起点坐标的中位数以及该簇内全部sv信号的终点坐标的中位数，将计算出的中位数作为候选sv信号的最终断点坐标；
9、否则，利用该簇内的sv信号的特征构建特征矩阵，将构建的特征矩阵输入到训练好的卷积神经网络中，通过卷积神经网络的输出来最终确定该簇是否为真实的sv簇；
10、若该簇为真实的sv簇，则计算出该簇内全部sv信号的起点坐标中位数和终点坐标中位数，将计算出的中位数作为候选sv信号的最终断点坐标；
11、若该簇不为真实的sv簇，则不需要处理；
12、同理，对步骤三中获得的每个簇分别进行处理后，获得全部候选sv信号的最终断点坐标，完成基因组结构变异的检测。
13、进一步地，所述cigar字符串的获得方式为：
14、将待检测长读测序数据直接与参考基因组进行比对时，得到的比对记录即为cigar字符串；
15、所述根据cigar字符串提取sv信号的具体过程为：
16、cigar字符串是由数字和字母组成的字符串，根据cigar字符串中的字母定位出缺失和插入两种类型的sv信号，并根据cigar字符串中的数字确定缺失和插入两种类型的sv信号的长度。
17、进一步地，所述缺口信息的获得方式为：
18、将待检测长读测序数据拆分为多个片段后，每个片段中两端的碱基即为片段的断点，断点的坐标即为缺口信息；
19、所述根据缺口信息提取sv信号包括识别sv信号和提取sv信号；其中：
20、识别sv信号的具体过程为：
21、(1)对于任意相邻的两个片段来说，均将其中一个片段作为主比对片段，将另一个片段作为补充比对片段，将主比对片段与补充比对片段在待检测长读测序数据中的距离记为dis_read，将主比对片段与补充比对片段在参考基因组中的距离记为dis_ref；
22、比较dis_read与dis_ref的大小：
23、若dis_read<dis_ref，则主比对片段与补充比对片段之间存在缺失sv信号；
24、若dis_read>dis_ref，则主比对片段与补充比对片段之间存在插入sv信号；
25、若dis_read＝dis_ref，则主比对片段与补充比对片段之间不存在sv信号；
26、(2)若主比对片段在参考基因组上的比对位置与补充比对片段在参考基因组上的比对位置存在重叠，则主比对片段与补充比对片段之间存在重复sv信号；
27、(3)若主比对片段相对于参考基因组的方向和补充比对片段相对于参考基因组的方向不同，则主比对片段与补充比对片段之间存在倒位sv信号；
28、(4)若主比对片段在参考基因组中对应的染色体号与补充比对片段在参考基因组中对应的染色体号不相同，则主比对片段和补充比对片段之间存在易位sv信号；
29、提取sv信号的具体过程为：将识别出的sv信号提取出来。
30、进一步地，所述主比对片段与补充比对片段在待检测长读测序数据中的距离dis_read为：
31、若在待检测长读测序数据中，主比对片段位于补充比对片段之前，则距离dis_read为补充比对片段的起点碱基坐标与主比对片段的终点碱基坐标的差值；
32、若在待检测长读测序数据中，主比对片段位于补充比对片段之后，则距离dis_read为主比对片段的起点碱基坐标与补充比对片段的终点碱基坐标的差值。
33、进一步地，所述根据cigar字符串提取的sv信号的特征包括全局比对特征、局部cigar特征和碱基分布特征；其中：
34、全局比对特征包括sv信号的比对质量、软剪切占比、软剪切长度、flag值、待检测长读测序数据的长度、变异数量、比对的分区号、sv信号在待检测长读测序数据中的起始位置、sv信号在待检测长读测序数据中的结束位置、待检测长读测序数据在参考基因组的染色体号、待检测长读测序数据在参考基因组的起始位置、待检测长读测序数据在参考基因组的结束位置、比对的编辑距离以及不匹配碱基长度；
35、局部cigar特征包括待检测长读测序数据中sv信号在参考基因组的起点、sv信号在参考基因组终点、sv信号长度、sv信号类型、sv信号在cigar中的起始位置以及sv信号在cigar中的结束位置；
36、碱基分布特征包括sv信号中四个碱基分别出现的频率以及重复碱基的出现频率。
37、进一步地，所述根据缺口信息提取的sv信号的特征包括全局比对特征，拆分比对特征和碱基分布特征；其中：
38、全局比对特征与根据cigar字符串提取的sv信号的全局比对特征相同；
39、碱基分布特征与根据cigar字符串提取的sv信号的碱基分布特征相同；
40、拆分比对特征包括主比对与补充比对在参考基因组中对应的染色体号一致性、主比对与补充比对在参考基因组中方向一致性、主比对与补充比对在待检测长读测序数据中的距离d1、主比对与补充比对在参考基因组中的距离d2以及距离d1与d2之差。
41、进一步地，所述随机森林网络的训练过程为：
42、步骤1、将已知的长读测序数据与参考基因组进行比对，得到缺口信息，再根据缺口信息提取sv信号；并记录根据缺口信息提取的sv信号的特征；
43、步骤2、按照7：3的比例，将步骤1中提取的1号染色体至11号染色体上的全部sv信号划分为两部分，将其中一部分作为训练集，另一部分作为验证集，并根据基准集对训练集和验证集中的sv信号进行标注；将步骤1中提取的12号染色体至22号染色体、性染色体上的全部sv信号作为测试集；
44、sv信号的标注方式为：
45、缺失、插入、重复以及倒位的sv信号的标注方法相同，以插入的sv信号为例：
46、将提取到的任一插入sv信号的起点坐标记为s1，并将该插入sv信号的长度记为size1，判断基准集中是否存在与该插入sv信号在同一条染色体，且满足下式条件的插入sv信号：
47、
48、其中，s2是基准集中的插入sv信号的起点坐标，|·|代表取绝对值，size2是基准集中的插入sv信号的长度；
49、若基准集中存在与该插入sv信号在同一条染色体且满足条件的插入sv信号，则提取到的该插入sv信号为真实插入sv信号，将提取到的该插入sv信号标注为插入sv信号；否则，基准集中不存在与该插入sv信号在同一条染色体且满足条件的插入sv信号，则将提取到的该插入sv信号标注为误检；
50、易位的sv信号的标注方法为：
51、将提取到的任一易位sv信号的起点坐标记为s′1，判断基准集中是否存在与该易位sv信号在同一条染色体，且满足下式条件的易位sv信号：
52、|s′1-s′2|≤1000
53、其中，s′2是基准集中的易位sv信号的起点坐标；
54、若基准集中存在与该易位sv信号在同一条染色体且满足条件的易位sv信号，则提取到的该易位sv信号为真实易位sv信号，将提取到的该易位sv信号标注为易位sv信号；否则，基准集中不存在与该易位sv信号在同一条染色体且满足条件的易位sv信号，则将提取到的该易位sv信号标注为误检；
55、步骤3、利用步骤1中根据缺口信息提取的sv信号特征构造特征矩阵，再对构造的特征矩阵进行最小最大值归一化处理，得到处理后的特征矩阵；
56、利用训练集中标注后的sv信号和处理后的特征矩阵对随机森林网络进行训练。
57、进一步地，所述步骤三中的聚类过程是对各个类型的sv信号分别进行的；以插入sv信号为例：
58、步骤三一、将步骤一中根据cigar字符串提取和步骤二中确定出的每个插入sv信号分别表示为四元组(t，c，s，e)的形式，其中，t代表sv信号类型，c代表sv信号所在的染色体号，s代表sv信号在参考基因组中的起点坐标，e代表sv信号在参考基因组中的终点坐标；
59、对各个插入sv信号进行排序：即对于一个染色体上的全部插入sv信号，按照插入sv信号在参考基因组中的起点坐标由小到大的顺序对各个插入sv信号依次进行排序，若排序后相邻两个插入sv信号的空间距离小于1000bp，则相邻两个插入sv信号被划分到同一分区，否则相邻两个插入sv信号被划分到不同的分区；
60、同理，对各染色体上的插入sv信号分别进行排序和分区；
61、步骤三二、分别获取每个分区的sv局部深度以及全部分区的sv全局深度：
62、
63、其中，a是sv全局深度，bi是第i个分区的sv局部深度，n是分区的总个数；
64、步骤三三、对于第i个分区内的插入sv信号进行聚类，i＝1,2,…,n；聚类的具体过程为：
65、步骤三三一、对于第i个分区内的第k个插入sv信号以及第l个插入sv信号，计算第k个插入sv信号与第l个插入sv信号的相似度s(k,l)：
66、
67、其中，pos_dis是第k个插入sv信号与第l个插入sv信号的位点距离，span_dis是第k个插入sv信号与第l个插入sv信号的跨度距离，sk和ek分别为第k个插入sv信号在参考基因组中的起点坐标和终点坐标，sl和el分别为第l个插入sv信号在参考基因组中的起点坐标和终点坐标，λ是修正参数；
68、
69、span_dis＝|(ek-sk)-(el-sl)|
70、
71、步骤三三二、根据计算的相似度对第i个分区内的插入sv信号进行层次聚类，获得各个簇；
72、步骤三四，对步骤一中根据cigar字符串提取和步骤二中确定出的缺失sv信号、重复sv信号、倒位sv信号和易位sv信号分别执行步骤三一至步骤三三的过程。
73、进一步地，所述排序后，相邻两个插入sv信号的空间距离为sk+1-sk，其中，sk是第k个插入sv信号的起点坐标，sk+1是第k+1个插入sv信号的起点坐标；
74、相邻两个易位sv信号的空间距离计算方法与相邻两个插入sv信号的空间距离计算方法相同；
75、相邻两个缺失sv信号的空间距离为s′k+1-e′k，其中，e′k是第k个缺失sv信号的终点坐标，s′k+1是第k+1个缺失sv信号的起点坐标；
76、相邻两个倒位sv信号的空间距离、相邻两个重复sv信号的空间距离的计算方法均与相邻两个缺失sv信号的空间距离计算方法相同。
77、更进一步地，所述卷积神经网络的训练方法为：
78、步骤s1、将已知的长读测序数据与参考基因组进行比对，得到cigar字符串，再根据cigar字符串提取sv信号；并记录根据cigar字符串提取的sv信号的特征；
79、步骤s2、根据步骤s1中已知长读测序数据的缺口信息提取sv信号特征，利用根据缺口信息提取的sv信号的特征构造特征矩阵，对构造的特征矩阵进行最小最大值归一化处理后，将处理后的特征矩阵输入随机森林网络，通过随机森林网络预测出提取的sv信号中的真实sv信号；
80、步骤s3、对步骤s1中根据cigar字符串提取的sv信号和步骤s2确定的真实sv信号一起进行聚类，得到包含sv信号的各个簇；
81、步骤s4、确定出比对内部提取的sv信号的簇，根据每个簇的特征矩阵和标签对卷积神经网络进行训练；具体为：
82、步骤s41、对于确定出的任意一个簇，判断该簇内的sv信号个数是否处于[1,20]范围内；
83、若不处于[1,20]范围内，则执行步骤s42；
84、若处于[1,20]范围内，则执行步骤s45；
85、步骤s42、继续判断该簇内的sv信号个数是否小于等于100；
86、若该簇内的sv信号个数小于等于100，则执行步骤s43；
87、若该簇内的sv信号个数大于100，则执行步骤s44；
88、步骤s43、直接利用簇内的sv信号特征构建维度为(100,25)的特征矩阵，若簇内sv信号的个数小于100，需要通过填零操作获得特征矩阵，并给特征矩阵添加标签；
89、并对构建的特征矩阵进行最小-最大值归一化处理，得到一个处理后的特征矩阵；
90、步骤s44、从该簇内随机选取出100个sv信号来执行步骤s43，并判断该簇内剩余sv信号的个数是否处于[1,20]范围内；
91、若不处于[1,20]范围内，则执行步骤s42；
92、若处于[1,20]范围内，则执行步骤s45；
93、步骤s45、将簇内的sv信号的个数记为w，利用小于w的窗口大小对簇内的各个sv信号分别进行一次数据截断，再将每次截断获得的窗口内的数据分别作为一个新的sv信号；并利用簇内截取之前的sv信号和截取之后得到的新sv信号的特征一起构建特征矩阵，并给特征矩阵添加标签；
94、再对构建的特征矩阵进行最小-最大值归一化处理，得到一个处理后的特征矩阵；
95、步骤s5、采用步骤s4的方法对确定出的每个簇分别进行处理；
96、步骤s6、利用获得的处理后特征矩阵和添加的标签对卷积神经网络进行训练。
97、步骤s6、利用获得的处理后特征矩阵和添加的标签对卷积神经网络进行训练。
98、本发明的有益效果是：
99、本发明以个体长读测序数据作为输入，通过与参考基因组的比对初步获得sv信号，再将随机森林网络和卷积神经网络作为假阳性过滤模型，并重新设计sv信号的聚类方法，从初步获得的sv信号中过滤掉假阳性sv信号，可以有效降低假阳性事件的调用可能性，具备快速、准确地检测来自多种长读测序平台和各类sv信号的能力。
1.一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述方法具体包括以下步骤：
2.根据权利要求1所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述cigar字符串的获得方式为：
3.根据权利要求2所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述缺口信息的获得方式为：
4.根据权利要求3所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述主比对片段与补充比对片段在待检测长读测序数据中的距离dis_read为：
5.根据权利要求4所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述根据cigar字符串提取的sv信号的特征包括全局比对特征、局部cigar特征和碱基分布特征；其中：
6.根据权利要求5所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述根据缺口信息提取的sv信号的特征包括全局比对特征，拆分比对特征和碱基分布特征；其中：
7.根据权利要求6所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述随机森林网络的训练过程为：
8.根据权利要求7所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述步骤三中的聚类过程是对各个类型的sv信号分别进行的；以插入sv信号为例：
9.根据权利要求8所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述排序后，相邻两个插入sv信号的空间距离为sk+1-sk，其中，sk是第k个插入sv信号的起点坐标，sk+1是第k+1个插入sv信号的起点坐标；
10.根据权利要求9所述的一种基于长读测序数据和假阳性过滤模型的基因组结构变异检测方法，其特征在于，所述卷积神经网络的训练方法为：
