阅读说明：本技术 一种考虑拷贝数变异因素的基因组结构变异分型方法 (Genome structure variation typing method considering copy number variation factors ) 是由 王嘉寅 郑田 朱晓燕 童瑶 刘佳琦 张选平 于 2020-05-06 设计创作，主要内容包括：本发明公开了一种考虑拷贝数变异因素的基因组结构变异分型方法,输入序列比对文件和突变识别文件并统计记录各变异位点的特征值；根据输入文件提取特征值,从突变识别文件VCF中提取基因型作为分类监督,通过Python提取VCF文件中第八列type后的基因型信息,一行对应一个变异依次将特征值以空格分隔,基因型以分隔符存储到txt文件中；确定核函数和核函数参数；将数据分为M-RVM模型的训练集和测试集；采用快速二类极大似然估计求解先验参数,采用最大期望估计算法求解核参数；输出分型结果、估计概率和总体精度。本方法全面理清了考虑拷贝数变异因素的基因组结构变异分型问题,利用多分类相关向量机设计了一种高准确率、高效率的解法。(The invention discloses a genome structure variation typing method considering copy number variation factors, which comprises the steps of inputting a sequence comparison file and a mutation identification file and counting and recording characteristic values of various variation sites; extracting characteristic values according to an input file, extracting genotypes from a mutation identification file VCF as classification supervision, extracting genotype information behind an eighth row of types in the VCF file through Python, sequentially separating the characteristic values by spaces in a row corresponding to one mutation, and storing the genotypes into a txt file by separators; determining a kernel function and kernel function parameters; dividing data into a training set and a testing set of an M-RVM model; solving prior parameters by adopting fast second-class maximum likelihood estimation, and solving kernel parameters by adopting a maximum expectation estimation algorithm; and outputting a typing result, an estimated probability and overall accuracy. The method comprehensively solves the problem of genome structure variation typing considering copy number variation factors, and designs a high-accuracy and high-efficiency solution by using a multi-classification relevance vector machine.)

一种考虑拷贝数变异因素的基因组结构变异分型方法

技术领域

本发明属于数据科学技术领域，具体涉及一种考虑拷贝数变异因素的基因组结构变异分型方法。

背景技术

拷贝数变异(英文全称：Copy Number Variant；英文缩写：CNV)指一种拷贝数变化的基因组结构变异，导致基因组增加(复制或插入)、缺失(删除)或复杂重排，涉及的DNA片段通常长于1千个碱基对(bps)，广泛存在于人类肿瘤数据中。结构变异基因分型是一种确定结构变异基因的类型的技术，二倍体生物人类的基因型有两种：同一基因位点两条链碱基类型相同为纯合型，不同为杂合型，区分结构变异纯合或杂合型的过程即为基因分型。正确的基因分型可广泛应用于下游分析，如填补空白基因、填补缺失的基因型、药物设计、计算基因组多样性与连锁不平衡、疾病诊断、指导治疗和管理药物反应等。

现有结构变异分型方法主要有基于双末端读对的BreakDancer、PEMer和Ulysses；基于拆分读对来识别断点信号，包括Pindel、Gustaf、SVseq2和Prism；基于读对深度信号，包括BIC-seq、CNVnator、CNVrd2、ERDS、ReadDepth、EWT、JointSLM、CNVeg和cn.MOPS；以及基于装配的方法，包括Velvet、EULER-USR、SOAPdenovo、ALLPATHS.LG和Magnolya。从算法的角度概括，已有方法可以分为三类，第一类以Pindel-C为代表，基于局部重叠点和断点估计基因型；第二类包括piCALL和MATE-CLEVER，利用基于人群的面板或先验系谱的贝叶斯框架；第三类以Gindel和CIGenotyper为代表，利用机器学习模型集成一系列特征来识别基因型，全面利用序列数据信息。

然而，现有方法存在以下问题：

1)现有方法信号利用不全面，大多数方法只能利用一种识别特征，没有全面利用识别信号；

2)不适用于存在拷贝数变异的基因分型，现有的方法中没有一种试图对CNV进行基因分型，大多数方法都集中在简单的某种结构变异上，例如GINDEL只识别插入、删除两种变异；在对所有结构变异进行基因分型时普遍呈现严重的精度损失，识别效果差，不能解决拷贝数变异引起的等位基因失衡问题。CNV会导致人类基因分型中杂合变异被误判为纯合，当CNV出现在杂合变异的变异侧，插入或游离于基因组时，测序数据会显示变异区出现比正常区更多的读对，导致各种工具将这种杂合情况误判为双链纯合变异。当CNV出现在杂合变异位点的正常侧时，测序结果会观察到更多正常区域的读对，淹没变异侧读对信息，进而导致误判为无变异纯合。拷贝数变异会使现有方法分型特征失效，出现假阳性和假阴性，这对基因分型在肿瘤数据中的应用准确性影响极大；

3)存在不同标签基因型的特征数据相等或差别不大的情况，无法区分拷贝数变异在杂合变异侧和纯合变异、正常无变异和拷贝数变异在杂合无变异侧等，加深了分型困难。

另外，现有方法需要大量的计算资源，很难应用于包含大量重复和片段重复的真核基因组，不能处理杂合基因型。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种考虑拷贝数变异因素的基因组结构变异分型方法，面向第二代基因测序数据，当基因组结构变异含有大量拷贝数变异因素的情况下，使用机器学习策略实现基因分型的问题。

本发明采用以下技术方案：

一种考虑拷贝数变异因素的基因组结构变异分型方法，包括以下步骤：

S1、输入序列比对文件和突变识别文件，序列比对文件符合SAM格式规范，和突变识别文件符合VCF格式规范；

S2、对于突变识别文件中的每个突变，分别从序列比对文件、突变识别文件中提取异常读对、正常读对、未完全匹配读对、完全匹配读对、拆分读对、单端读对、不能匹配读对、匹配质量、读对深度、加权读对深度、扩展加权读对深度、受影响读对特征、结构变异长度、读取变异区与5'端匹配的读对数量统计值和变异区与3'端匹配的读对数量统计值；

S3、确定核函数和核函数参数；

S4、将数据分为M-RVM模型的训练集和测试集；

S5、采用快速二类极大似然估计求解先验参数，采用最大期望估计算法求解核参数；

S6、输出分型结果、估计概率和总体精度。

具体的，步骤S2中，提取异常读对和正常读对具体为：

S20101、从VCF文件中提取第二列POS值和第八列INFO中End＝数值，每一行对应一个变异，分别表示变异的起始和终止匹配位置，将变异的起始和终止匹配位置[POS,End]区间定义为每个变异的变异区间；

S20102、从SAM文件中提取第四列POS和第九列ISIZE，对整个文件第九列全部ISIZE求和计算平均值记为μ，计算标准差记为σ；

S20103、对VCF中每个变异计算的变异区间[POS,End]，提取SAM文件第四列POS值中位于区间内的行，分别统计其中第九列ISIZE取值超出和位于μ±3σ区间的行的数量；

S20104、将超出区间的行数记为异常读对值，符合区间的数量记为正常读对值。

具体的，步骤S2中，提取未完全匹配读对、完全匹配读对、拆分读对、单端读对和不能匹配读对具体为：

S20201、从标准VCF文件中提取第二列POS值和第八列INFO中End＝数值，每一行对应一个变异，分别表示变异的起始和终止匹配位置，将变异的起始和终止匹配位置[POS,End]区间定义为每个变异的变异区间；

S20202、对VCF中每个变异对应的[POS,End]区间选择SAM文件中第四列POS值位于其内的行，并提取对应行中第六列CIGAR值，统计其中取值为100M的行数记为完全匹配读对特征值；统计CIGAR值不等于100M的行数记为未完全匹配读对；

S20203、SAM文件对VCF中每个变异对应的[POS,End]区间每两行连续读取SAM文件，统计第四列POS值一个位于变异区间、一个不在变异区间内的读对数量，记为拆分读对值；

S20204、对VCF中每个变异对应的[POS,End]区间，连续读取SAM文件的两行，统计一行的第六列CIGAR值等于100M，另一行的第六列CIGAR值不等于100M的行数，记为单端匹配特征值；

S20205、提取变异区间内的两个读段都不能匹配到参考序列上的读对，得到不能正常匹配到参考序列上的读对信息文件SAM文件二进制格式，在得到的SAM文件中，统计行中第四列POS值落在对应变异区间的行数，记为对应变异的不能匹配读对特征值。

具体的，步骤S2中，提取匹配质量具体为：

S20301、从VCF文件中提取第二列POS值和第八列INFO中End＝数值，每一行对应一个变异，分别表示变异的起始和终止匹配位置，将变异的起始和终止匹配位置[POS,End]区间定义为每个变异的变异区间；

S20302、对VCF中每个变异提取SAM文件中第四列POS值中位于[POS,End]区间的行相应的第五列MAPQ值，将行中对应的MAPQ值求和，记为变异对应的匹配质量特征值。

具体的，步骤S2中，提取读对深度、加权读对深度和扩展加权读对深度具体为：

S20401、从VCF文件中提取第二列POS值和第八列INFO中End＝数值，每一行对应一个变异，分别表示变异的起始和终止匹配位置，将变异的起始和终止匹配位置[POS,End]区间定义为每个变异的变异区间；

S20402、对VCF中每个变异对应的[POS,End]区间统计SAM文件中第四列POS值位于变异区间的行数，记为n，计算n/|POS-End|并记为变异的读对深度特征值；

S20403、对VCF文件每一行代表的单个变异，存储匹配质量值特征记为匹配质量，计算SAM文件中第四列POS值位于[POS,End]区间的行的数量，记为匹配读对数量，提取SAM文件第六列CIGAR最大值记为完全匹配质量，计算|POS-End|记做变异长度；

S20404、将加权读对深度计算中的变异区域上下游各扩展100bps长度，对VCF文件每一行代表的单个变异，存储匹配质量值记为匹配质量特征值，将SAM文件中第四列POS值位于[[POS,End]-100，[POS,End]+100]区间的行的数量，记为扩展区域匹配读对数量，提取SAM文件第六列CIGAR最大值记为完全匹配质量，计算|POS-End|记做变异长度。

具体的，步骤S2中，提取受影响读对特征具体为：

S20501、从VCF文件中提取第二列POS值和第八列INFO中End＝数值，计算|POS-End|记做变异长度length，设定[POS-length/10,End+length/10]为变异影响区间；

S20502、每两行读取SAM文件，对VCF中每个变异提取SAM文件中第四列POS值位于变异影响区间[POS-length/10，End+length/10]的行，记为影响读对值。

具体的，步骤S2中，提取读取方向具体为：

S20701、从VCF文件中提取第二列POS值和第八列INFO中End＝数值，将此位置[POS,End]区间定义为每个变异的变异区间；

S20702、每两行一组，连续读取SAM文件，对VCF中每个变异提取SAM文件中第四列POS值位于[POS,End]区间的行，统计其中第二列FLAG为83和163的行数记为读取方向1；

S20703、每两行一组，连续读取SAM文件，对VCF中每个变异提取SAM文件中第四列POS值位于[POS,End]区间的行，统计第二列FLAG不为83和163的行数记为读取方向2。

具体的，步骤S3中，核函数具体为：

其中，γ为核函数参数。

具体的，步骤S4具体为：

S401、将基因类型向量化表示，将基因型状态分为五种类型：无变异的正常基因型、无CNV的纯合变异、无CNV的杂合变异、变异侧出现CNV的杂合变异、正常侧出现CNV的杂合变异，分别由向量[0,0,0,0,1]^T、[0,0,0,1,0]^T、[0,0,1,0,0]^T、[0,1,0,0,0]^T、[1,0,0,0,0]^T表示。

S402、利用k折交叉验证确定训练集和测试集，通过命令行自定义交叉验证的重数。

具体的，步骤S5中，通过求解先验参数，核参数为：

采用E步求得Y为：

其中，K指核函数集，K＝[k₁,…,k_n]^T,K∈R^N×N,为了确保模型的稀疏性，为权重向量引入零均值，方差为的标准正态先验分布，A由先验参数α_nc组成的矩阵，α_nc服从超参数为τ，υ的Gamma分布；Y为辅助回归目标，c为目标类别，n为训练样本数量，i为类别标记，为边缘似然函数。

与现有技术相比，本发明至少具有以下有益效果：

本发明一种考虑拷贝数变异因素的基因组结构变异分型方法，提出了将多分类相关向量机模型应用于基因分型的方法，并将分类类别从纯合型杂合型两类扩充到了五类，即不包含拷贝数变异的杂合型、不包含拷贝数变异的纯合型、包含拷贝数变异的杂合型、不包含拷贝数变异的纯合型和无变异纯合型，解决了对包含拷贝数变异在内的多种复杂结构变异的基因分型问题；融合了现有特征提取方法，并在综合利用现有特征信号的基础上，提出15种全新有效的特征信号，包括异常读对、正常读对、未完全匹配读对、完全匹配读对、拆分读对、单端读对、不能匹配读对、匹配质量、读对深度、加权读对深度、拓展的加权读对深度、受影响读对、变异长度和读对的方向。避免了使用组装方法，节省了大量计算资源，能将计算效率和结果精度平衡在高水平；利用期望最大化算法和快速Type-II最大似然来学习和求解，降低了相关向量依赖，优化计算效率。

进一步的，能够识别对读对insert size造成影响的变异，进而区分出变异和未变异，同时，拷贝数变异的存在会影响取值，根据数值大小，可以区分是否有拷贝数变异。

进一步的，在对五种基因型区分时，变异的存在会导致读对不能完美匹配到参考序列上，读对与变异的位置关系、变异的类型等因素会造成读对完全、完全不能、不能完全地与参考序列匹配，提取未完全匹配读对、完全匹配读对、拆分读对、单端读对和不能匹配读对出每种情况对应的读对的数量，能够提示变异是否发生以及变异的类型，区别于现有方法，将不同匹配情况的读对区分开能够扩大特征值的区分意义，实现高效精确的基因分型。

进一步的，提取匹配质量能够反映读对与参考基因的匹配情况，通过匹配质量值的大小和全部变异区域匹配质量值计算，不同基因型和拷贝数变异会造成匹配质量值有显著差异，这一特征值有利于基因分型。

进一步的，读对深度能反映基因正常、变异、拷贝数的差异，从而利于基因分型，通过匹配质量对读对深度信号进行加权，能够进一步扩大三者差异；同时，扩展变异区间，能够融合变异区间周围信息，且上下游各扩展一个读段的长度能够只融合可能被变异影响的所有读段，避免过多正常区域的影响，从而提高特征值的分型能力，有助于高效基因分型。

进一步的，提取受影响读对能够反映变异区间，反映变异类型，有助于基因分型。

进一步的，取读取方向特征有助于判断倒置等结构变异，对基因分型显著有效。

进一步的，确定核函数和核函数参数，能够在线性核函数、多项式核函数、高斯核函数、径向基核函数等中确定适合本方法应用问题的核函数和参数，并指导方法取得最好结果。

进一步的，将数据分为M-RVM模型的训练集和测试集，利用训练集拟合模型，通过设置分类器的参数，训练分类模型。通过训练集得出最优模型后，使用测试集进行模型预测，以用来衡量该最优模型的性能和分类能力，进行模型性能评价。

综上所述，本方法全面理清了考虑拷贝数变异因素的基因组结构变异分型问题，利用多分类相关向量机设计了一种高准确率、高效率的解法。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明流程图；

图2为正确基因型类别图；

图3为真实数据结果图；

图4为核参数对精度影响图；

图5为模型学习迭代过程中相关向量的变化图，其中，(a)为拷贝数为2，覆盖度为5时的相关向量变化图，(b)为拷贝数为3，覆盖度为10时的相关向量变化图，(c)为拷贝数为4，覆盖度为15的相关向量变化图，(d)为拷贝数为5，覆盖度为20的相关向量变化图。

具体实施方式

本发明提供了一种考虑拷贝数变异因素的基因组结构变异分型方法，在面对包含拷贝数变异在内的多种复杂结构变异的基因分型问题时，应区分五种可能的情况：

1)正常无变异型(又称野生型)；

2)不存在拷贝数变异的纯合变异型；

3)不存在拷贝数变异的杂合变异型；

4)存在拷贝数变异的纯合变异型；

5)存在拷贝数变异的杂合变异型。

判断拷贝数变异发生在变异侧或非变异侧，进而区分基因型，避免混为一谈的误判，提高分型正确率和精确度。

请参阅图1，本发明一种考虑拷贝数变异因素的基因组结构变异分型方法，包括以下步骤：

S1、输入序列比对文件(英文全称：sequence alignment/map format；英文缩写：SAM文件)和突变识别文件(英文全称：Variant call file；英文缩写：VCF文件)；

从SAM文件和VCF文件中统计并记录各变异位点的特征值，基于上述两个文件作为输入，本发明对每个变异提取了15个特征值。

S2、根据输入文件提取特征值；

S201、提取异常读对和正常读对；

特征异常读对和正常读对分别指插入大小异常的读对的统计值和插入距离服从正常分布的读对的统计值；SAM文件中ISIZE值服从正态分布，μ是平均插入大小，σ是标准偏差。

当发生结构变异时，变异区间与参考序列的不匹配会导致位于变异区间的读对匹配到参考序列的位置距离变化，反映到测序读对的ISIZE值变化。进一步地，拷贝数变异会导致变异或者非变异区间的数量扩大，会造成ISIZE值变化的读对的数量区别于普通变异情况，故而，统计这两种特征可以有效利于区分五种基因型。

S20101、从VCF文件中提取第二列POS值和第八列INFO中“End＝”数值，每一行对应一个变异，分别表示该变异的起始和终止匹配位置，将这一位置[POS,End]区间定义为每个变异的变异区间；

S20102、从SAM文件中提取第四列POS和第九列ISIZE，对整个文件第九列全部ISIZE求和计算平均值记为μ，计算标准差记为σ；

S20103、对VCF中每个变异(每一行)计算的变异区间[POS,End]，提取SAM文件第四列POS值中位于区间内的行，分别统计其中第九列ISIZE取值超出和位于(μ±3σ)区间的行的数量；

S20104、将超出区间的行数记为异常读对值，符合区间的数量记为正常读对值。

这一对特征，能够识别对读对的插入大小(英文全称：inset size)造成影响的变异，进而区分出变异和未变异，同时，拷贝数变异的存在会影响取值，根据数值大小，可以进一步区分是否有拷贝数变异，进而完成基因分型。

S202、提取未完全匹配读对、完全匹配读对、拆分读对、单端读对和不能匹配读对；

特征未完全匹配读对指不能完全匹配的读对的统计值；完全匹配读对指匹配质量为百分之百匹配的读对的统计值；拆分读对指不能完全匹配，拆分匹配的读对数量；单端读对指一对读对中只有一个能完全匹配的读对数量；不能匹配读对指不能完全匹配的读对统计值。

变异与参考序列的不匹配会导致变异区间读对不能完全匹配到参考基因，读对与变异区间的不同位置关系会导致读对呈现完全匹配、完全不能匹配、不能完全匹配等比对结果，不同比对结果的读对的数量能够揭示不同基因型。

S20201、从VCF文件中提取第二列POS值和第八列INFO中“End＝”数值，每一行对应一个变异，分别表示该变异的起始和终止匹配位置，将这一位置[POS,End]区间定义为每个变异的变异区间；

S20202、对VCF中每个变异(每一行)对应的[POS,End]区间选择SAM文件中第四列POS值位于其内的行，并提取对应行中第六列CIGAR值，统计其中取值为“100M”的行数记为完全匹配读对特征值；统计CIGAR值不等于“100M”的行数记为未完全匹配读对；

S20203、基于SAM文件格式和双末端测序原理，对VCF中每个变异(每一行)对应的[POS,End]区间每两行(一对读对)连续读取SAM文件，统计第四列POS值一个位于变异区间、一个不在变异区间内的读对数量，记为拆分读对值。

S20204、类似地，对VCF中每个变异(每一行)对应的[POS,End]区间，连续读取SAM文件的两行，统计一行的第六列CIGAR值等于“100M”，另一行的第六列CIGAR值不等于“100M”的行数，记为单端匹配特征值。

S20205、利用SAMtools工具提取变异区间内的两个读段都不能匹配到参考序列上的读对，得到不能正常匹配到参考序列上的读对信息文件SAM文件二进制格式(bam文件)，在得到的SAM文件中，统计行中第四列POS值落在对应变异区间的行数，记为对应变异的不能匹配读对特征值。

在对五种基因型区分时，变异的存在会导致读对不能完美匹配到参考序列上，读对与变异的位置关系、变异的类型等因素会造成读对完全、完全不能、不能完全地与参考序列匹配，一一提取出每种情况对应的读对的数量能够提示变异是否发生以及变异的类型，区别于现有方法，将不同匹配情况的读对区分能够进一步地扩大特征值的区分意义，实现高效精确的基因分型。

S203、提取匹配质量；

特征匹配质量的定义是变异区间内读对的匹配质量之和。SAM文件中MAPQ表示读对匹配质量，对VCF中每个变异，统计变异区间中读对的匹配质量之和记为匹配质量特征。由于结构变异的存在，不同区域中读对的匹配质量存在显著差异。不同区域读对的匹配质量的和可以用作识别基因型的特征。

S20301、从VCF文件中提取第二列POS值和第八列INFO中“End＝”数值，每一行对应一个变异，分别表示该变异的起始和终止匹配位置，将这一位置[POS,End]区间定义为每个变异的变异区间；

S20302、对VCF中每个变异(每一行)提取SAM文件中第四列POS值中位于[POS,End]区间的行相应的第五列MAPQ值，将这些行中对应的MAPQ值求和，记为该变异对应的匹配质量特征值；

这一特征能够反映读对与参考基因的匹配情况，通过匹配质量值的大小和全部变异区间匹配质量值计算，不同基因型和拷贝数变异会造成匹配质量值有显著差异，因此，这一特征值有利于基因分型。

S204、提取读对深度、加权读对深度和扩展加权读对深度；

读对深度(英文全称：Read Depth；英文缩写：RD)特征指匹配到特定位点或基因区域的读对数量；加权读对深度特征指通过映射质量加权的变量区域的读对深度；拓展的加权读对深度指将计算加权读对深度的区间上下拓展100个碱基长度得到的统计值。假设测序过程是一致的，读对深度遵循随机分布(泊松分布或修正泊松分布)，且变异区间的读对匹配数量预计与该区域在测序样本中出现的次数成比例，与该区域的拷贝数之间存在相关性。因此，与非变异区间相比，变异区间的读对深度将减少或为零，拷贝数变异区间的读对深度增加，利用变异区间内所有读对数量和与变异区间长度的比值计算。

此外，本发明提出加权读对深度作为一个新特性，其公式为：加权RD＝(∑匹配质量×匹配读对数量)÷(完全匹配质量×变量长度)，各因子定义及计算方法见提取计算步骤。

S20401、从VCF文件中提取第二列POS值和第八列INFO中“End＝”数值，每一行对应一个变异，分别表示该变异的起始和终止匹配位置，将这一位置[POS,End]区间定义为每个变异的变异区间；

S20402、对VCF中每个变异对应的[POS,End]区间统计SAM文件中第四列POS值位于变异区间的行数，记为n，计算n/|POS-End|并记为该变异的读对深度特征值；

S20403、对VCF文件每一行代表的单个变异，存储匹配质量值特征记为匹配质量，计算SAM文件中第四列POS值位于[POS,End]区间的行的数量，记为匹配读对数量，提取SAM文件第六列CIGAR最大值记为完全匹配质量(一般为100M，记做100；60M记做60)，计算|POS-End|记做变异长度。即加权RD＝(∑匹配质量×匹配读对数量)÷(完全匹配质量×变异长度)

计算求得加权读段深度特征值。

S20404、考虑到变异区间和非变异区间的差异，将加权读对深度计算中的变异区域上下游各扩展100bps长度，该扩展区域的加权读取对深度可以更好地反映多个变异的基因型特征。与step3类似，对VCF文件每一行代表的单个变异，存储匹配质量值记为匹配质量特征值，将SAM文件中第四列POS值位于[[POS,End]-100,[POS,End]+100]区间的行的数量，记为扩展区域匹配读对数量，提取SAM文件第六列CIGAR最大值记为完全匹配质量，计算|POS-End|记做变异长度。计算扩展加权读段深度特征值依据：扩展加权读段深度＝(∑匹配质量×扩展区域匹配读对数量)/(完全匹配质量×变异长度)

读对深度能反映基因正常、变异、拷贝数的差异，从而利于基因分型，通过匹配质量对读对深度信号进行加权，能够进一步扩大三者差异；同时，扩展变异区间，能够融合变异区间周围信息，且上下游各扩展一个读段的长度能够只融合可能被变异影响的所有读段，避免过多正常区域的影响，从而提高特征值的分型能力，有助于高效基因分型。

S205、提取受影响读对特征；

受影响读对特征指受变异影响的读对数量统计值。结构变异发生时，结构变异附近的碱基变异概率较高且不同的变异可能导致不同的碱基变异概率。因此，本发明将受到结构变异影响的读对匹配数作为一个特征，可以很好地反映结构变异的存在，有利于基因型的识别。统计变异区间内读对的CIGAR值记录的不匹配的数量，记为该变异的受影响读对特征值。

S20501、从VCF文件中提取第二列POS值和第八列INFO中“End＝”数值，计算|POS-End|记做变异长度length。因为二代测序读对长度固定为100bps，本发明设定[POS-length/10,End+length/10]为变异影响区间。

S20502、每两行读取SAM文件，对VCF中每个变异提取SAM文件中第四列POS值位于变异影响区间[POS-length/10,End+length/10]的行，记为影响读对值。本发明将处于这一区域的读对定义为受影响读对，理论上，在这一区域的读对必定会由于变异的影响而不能正常匹配，不能正确显示POS值，甚至不会出现在SAM文件中，因此，提取这一特征能够反映变异区间和变异类型，有助于基因分型。

S206、提取结构变异长度；

变异长度特征即指结构变异的长度。结构变异长度可以反映变异本身的信息，从VCF文件中END与POS之间的差值绝对值提取每个变异的长度信息。

从VCF文件中提取第二列POS值和第八列INFO中“End＝”数值，对每一行计算|POS-End|值，记为变异长度。已有研究证实，不同类型的结构变异的变异长度有区别，拷贝数变异一般在三千碱基对以上，提取变异长度能够有助于基因分型。

S207、提取读取方向；

本发明提取方向1特征，为变异区与5'端匹配的读对数量统计值；方向2特征为变异区与3'端匹配的读对数量统计值。根据双端测序原理，一对读对读取的方向相反，当包含倒置变异时，由于反向互补性，读对方向会出现例外，因此定义读对方向作为特征，统计变异区间内与参考基因方向一致的从5'端开始的读对数量作为第十四个特征，和初始映射时从3'端的读对数量作为第十五个特征。

S20701、从VCF文件中提取第二列POS值和第八列INFO中“End＝”数值，将这一位置[POS,End]区间定义为每个变异的变异区间；

S20702、每两行一组，连续读取SAM文件，对VCF中每个变异提取SAM文件中第四列POS值位于[POS,End]区间的行，统计其中第二列FLAG为“83”和“163”的行数记为读取方向1；

S20703、每两行一组，连续读取SAM文件，对VCF中每个变异(每一行)提取SAM文件中第四列POS值位于[POS,End]区间的行，统计其中第二列FLAG不为“83”和“163”的行数记为读取方向2；提取变异方向不同的读对数目有助于判断倒置等结构变异，对基因分型显著有效。

S208、提取特征值的同时，本发明从VCF文件中提取基因型作为分类监督，通过Python提取VCF文件中第八列“type”之后的基因型信息；

S209、一行对应一个变异依次将特征值以空格分隔，基因型以tab分隔存储到txt文件中。

S3、确定核函数和核函数参数；

本方法选取高斯核函数，具体为：

其中，γ为核函数参数。

S4、将数据分为M-RVM模型的训练集和测试集；

多分类支持向量机(M-RVM)是Damoulasy于2008年提出的一种有监督机器学习的模式识别新方法，截至申请时，还未应用于基因分型问题，本发明将应用结果与现有方法比较，能够解决基于拷贝数变异的复杂结构变异的分型问题，并且分型精确度和计算速度都有显著提升。本部分输入是本发明步骤S1中得到的特征值txt文件作为输入；输出为基因分型结果，计算时间，类别概率，方法精度和多重交叉验证结果。

S401、将基因类型向量化表示，将基因型状态分为五种类型：无变异的正常基因型(N)、无CNV的纯合变异(G1)、无CNV的杂合变异(G2)、变异侧出现CNV的杂合变异(G3)、正常侧出现CNV的杂合变异(G4)，分别由向量[0,0,0,0,1]^T、[0,0,0,1,0]^T、[0,0,1,0,0]^T、[0,1,0,0,0]^T、[1,0,0,0,0]^T表示。

S402、利用k折交叉验证(英文全称：k-fold cross validation)确定训练集和测试集，通过命令行自定义交叉验证的重数。

S5、采用快速二类极大似然估计求解先验参数，采用最大期望估计算法求解核参数；

引入辅助回归目标Y∈R^C×N和权重参数W∈R^C×N，由先验参数α_nc组成的矩阵记为A∈R^N×C,α_nc服从超参数为τ，υ的Gamma分布，边缘似然函数为

P(Y|K,A)＝log∫P(Y|K,W)P(W|A)dw

通过对数方法取得，其中，

求解即可得到先验参数。

最大期望估计算法采用M步，求得

采用E步求得Y：

S6、输出分型结果、估计概率和总体精度。

在得到分型结果之后，准确的SVs基因型就可以指导下游分析，如基因型的输入、基因组多样性的估计、连锁不平衡的计算以及疾病诊断、治疗管理和药物设计等临床实践。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面通过仿真数据和真实数据实验及对比结果证明本发明的实用性、准确性和高效性：

1)仿真数据实验结果

为了评估所提出方法的性能，本发明依照通用实验配置，从人类参考基因组19号染色体上随机抽取了1Mbps区域，然后为每个数据集随机植入300个结构变异。变异的类型包括插入，缺失，倒位，复杂插入/缺失和CNV。根据CNV的分布概率，本发明创建了60个N型变异(无变异的正常纯合基因型)，80个G1类型变异(无CNV的纯合变异)，80个G2类型变异(无CNV的杂合变异)，50个G3类型变异(变异侧出现CNV的杂合变异)和30个G4型变异(正常侧出现CNV的杂合变异)。将变异长度设置为0.5～5kbps之间，CNV的长度设置为1～5kbps之间。

对于每个变异，本发明设置一个比其自身长度长1000bps的可能突变区域，在其中种植一些相关的SNV并将区域突变率设置为0.01，背景突变率设置为0.0001。大约四分之一的复杂插入缺失的插入片段来自附近区域。将测序读取长度设置为100bps，读对插入大小的分布遵循均值为500bps，标准偏差为15bps的正态分布，测序错误率设为0.005。

为了全面评估本发明的性能，本发明选择了三种现有方法进行比较：

1)通过Facets获取CNV区域，然后通过20％-80％规则估算基因型，Facets所需的snp.VCF.gz文件从NCBI公开数据库下载；

2)GATK，默认设置为HaplotypeCaller；

3)Gindel。

本发明比较了不同拷贝数的结果，对于每种拷贝数，本发明将覆盖度从5×提升到20×。实验结果如表1所示，每个结果均是五次重复实验的平均值。

表1不同覆盖率和拷贝数下的性能测试及对比实验

结果表明，本发明的准确性稳定在83％(±2％)以上。M-SVM的平均精度为45％，极差约为30％；Facets的平均值为52.79％，极差为10.70％；GATK的精度为63.8％，极差为57.95％。具体而言，随着覆盖率的增加，这些方法的准确性呈上升趋势，而随着拷贝数的增加，这些方法略有下降，这与理论原理是一致的。当覆盖率和拷贝数变化时，本发明显示出稳定的适应性，强大的鲁棒性，精度和稳定性都强于其他方法。

2)真实数据实验结果

本发明从公共数据库中获得了9组真实测序数据。原始数据已经在公共数据库上进行了脱敏处理，所有能够呈现和追踪个人信息的特征均全部删除，BWA-0.7.5a与参考基因匹配，GATK识别变异，CNVkit检测真实的结构变异信息。本发明通过处理离线数据获得了算法的输入SAM文件和VCF文件，并将所提出的方法与流行的机器学习框架M-SVM，Bayes，BP神经网络和Lanrange-SVM进行了比较，如表2所示。

表2真实数据实验结果

结果表明，本发明能够很好地适应真实数据，并且由于样本量大，真实数据覆盖率高，实验结果甚至优于模拟结果。与流行算法相比，本发明准确性的平均值为88.20％(±15％)，而M-SVM的平均值为79.17％(±20％)，朴素贝叶斯的精度均值为75.85％(±40％)，BP神经网络的平均精度为83.67％(±16％)，Lanrange-SVM的平均精度为68.43％(±28％)。这表明本发明保持了较高的准确性和稳定性，具有准确性和计算优势，可以很好地应用于临床。

请参阅图3，真实数据结果，其中X轴表示数据集序号，Y轴表示准确率(英文全称：Accuracy)，其中图例Our method代指本发明；M-SVM代指多分类支持向量机(英文全称：Multi-class Support Vector Machine)；NB代指朴素贝叶斯(英文全称：Bayes)；G-Features+NB代指朴素贝叶斯应用Gindel方法提取的简单特征；BP代指BP神经网络算法(英文全称Back-ProPagation Network)；G-Features+BP代指BP神经网络算法应用Gindel方法提取的简单特征；SVM with OVO代指一对一法支持向量机(英文全称Support VectorMachine with one-versus-one)算法；G-Features+SVM with OVO代指SVM with OVO方法应用Gindel提取的简单特征。

请参阅图4和图5，为模型学习迭代过程中相关向量的变化，其中X轴表示迭代计算次数，Y轴表示相关向量个数(Relevant Vectors)，其中，(a)为拷贝数为2，覆盖度为5时的相关向量变化图；(b)为拷贝数为3，覆盖度为10时的相关向量变化图；(c)为拷贝数为4，覆盖度为15的相关向量变化图；(d)为拷贝数为5，覆盖度为20的相关向量变化图。结果表明，该模型最终只需要少量相关向量，且具有较高的稀疏性。当模型确定后，新样本输入数据诊断的计算复杂度较低。

综上所述，本发明提出了将多分类相关向量机(英文全称：Multiclass RelevantVector Machine；英文缩写：M-RVM)模型应用于基因分型的方法，并将分类类别从纯合型杂合型两类扩充到了五类，即不包含拷贝数变异的杂合型、不包含拷贝数变异的纯合型、包含拷贝数变异的杂合型、不包含拷贝数变异的纯合型和无变异纯合型，解决了对包含拷贝数变异在内的多种复杂结构变异的基因分型问题。

本发明融合了现有特征提取方法，并在综合利用现有特征信号的基础上，提出15种全新有效的特征信号，包括异常读对、正常读对、未完全匹配读对、完全匹配读对、拆分读对、单端读对、不能匹配读对、匹配质量、读对深度、加权读对深度、拓展的加权读对深度、受影响读对、变异长度和读对的方向。避免了使用组装方法，节省了大量计算资源，能将计算效率和结果精度平衡在高水平。

本发明利用期望最大化算法(英文全称：Expectation-Maximum算法；英文缩写：EM算法)和快速Type-II最大似然(英文全称：Fast Type-II Maximum Likelihood；英文缩写：Fast Type-II ML)来学习和求解，降低了相关向量依赖，优化计算效率。

以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。