技术摘要：
本发明公开了一种考虑拷贝数变异因素的基因组结构变异分型方法，输入序列比对文件和突变识别文件并统计记录各变异位点的特征值；根据输入文件提取特征值，从突变识别文件VCF中提取基因型作为分类监督，通过Python提取VCF文件中第八列type后的基因型信息，一行对应一个 全部
背景技术：
拷贝数变异(英文全称：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)存在不同标签基因型的特征数据相等或差别不大的情况，无法区分拷贝数变异 在杂合变异侧和纯合变异、正常无变异和拷贝数变异在杂合无变异侧等，加深了分型困难。 另外，现有方法需要大量的计算资源，很难应用于包含大量重复和片段重复的真 核基因组，不能处理杂合基因型。 5 CN 111583998 A 说 明 书 2/15 页
技术实现要素：
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种考虑拷贝 数变异因素的基因组结构变异分型方法，面向第二代基因测序数据，当基因组结构变异含 有大量拷贝数变异因素的情况下，使用机器学习策略实现基因分型的问题。 本发明采用以下技术方案： 一种考虑拷贝数变异因素的基因组结构变异分型方法，包括以下步骤： 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的行数，记为单端匹 6 CN 111583998 A 说 明 书 3/15 页 配特征值； 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具体为： 7 CN 111583998 A 说 明 书 4/15 页 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＝[k1,…,k Tn] ,K∈RN×N,为了确保模型的稀疏性，为权重向量 引入零均值，方差为 的标准正态先验分布，A由先验参数αnc组成的矩阵，αnc服从超参数 为τ，υ的Gamma分布；Y为辅助回归目标，c为目标类别，n为训练样本数量，i为类别标记， 为边缘似然函数。 与现有技术相比，本发明至少具有以下有益效果： 本发明一种考虑拷贝数变异因素的基因组结构变异分型方法，提出了将多分类相 关向量机模型应用于基因分型的方法，并将分类类别从纯合型杂合型两类扩充到了五类， 即不包含拷贝数变异的杂合型、不包含拷贝数变异的纯合型、包含拷贝数变异的杂合型、不 包含拷贝数变异的纯合型和无变异纯合型，解决了对包含拷贝数变异在内的多种复杂结构 变异的基因分型问题；融合了现有特征提取方法，并在综合利用现有特征信号的基础上，提 出15种全新有效的特征信号，包括异常读对、正常读对、未完全匹配读对、完全匹配读对、拆 分读对、单端读对、不能匹配读对、匹配质量、读对深度、加权读对深度、拓展的加权读对深 度、受影响读对、变异长度和读对的方向。避免了使用组装方法，节省了大量计算资源，能将 计算效率和结果精度平衡在高水平；利用期望最大化算法和快速Type-II最大似然来学习 和求解，降低了相关向量依赖，优化计算效率。 进一步的，能够识别对读对insert size造成影响的变异，进而区分出变异和未变 异，同时，拷贝数变异的存在会影响取值，根据数值大小，可以区分是否有拷贝数变异。 进一步的，在对五种基因型区分时，变异的存在会导致读对不能完美匹配到参考 序列上，读对与变异的位置关系、变异的类型等因素会造成读对完全、完全不能、不能完全 地与参考序列匹配，提取未完全匹配读对、完全匹配读对、拆分读对、单端读对和不能匹配 读对出每种情况对应的读对的数量，能够提示变异是否发生以及变异的类型，区别于现有 方法，将不同匹配情况的读对区分开能够扩大特征值的区分意义，实现高效精确的基因分 型。 进一步的，提取匹配质量能够反映读对与参考基因的匹配情况，通过匹配质量值 8 CN 111583998 A 说 明 书 5/15 页 的大小和全部变异区域匹配质量值计算，不同基因型和拷贝数变异会造成匹配质量值有显 著差异，这一特征值有利于基因分型。 进一步的，读对深度能反映基因正常、变异、拷贝数的差异，从而利于基因分型，通 过匹配质量对读对深度信号进行加权，能够进一步扩大三者差异；同时，扩展变异区间，能 够融合变异区间周围信息，且上下游各扩展一个读段的长度能够只融合可能被变异影响的 所有读段，避免过多正常区域的影响，从而提高特征值的分型能力，有助于高效基因分型。 进一步的，提取受影响读对能够反映变异区间，反映变异类型，有助于基因分型。 进一步的，取读取方向特征有助于判断倒置等结构变异，对基因分型显著有效。 进一步的，确定核函数和核函数参数，能够在线性核函数、多项式核函数、高斯核 函数、径向基核函数等中确定适合本方法应用问题的核函数和参数，并指导方法取得最好 结果。 进一步的，将数据分为M-RVM模型的训练集和测试集，利用训练集拟合模型，通过 设置分类器的参数，训练分类模型。通过训练集得出最优模型后，使用测试集进行模型预 测，以用来衡量该最优模型的性能和分类能力，进行模型性能评价。 综上所述，本方法全面理清了考虑拷贝数变异因素的基因组结构变异分型问题， 利用多分类相关向量机设计了一种高准确率、高效率的解法。 下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。 附图说明 图1为本发明流程图； 图2为正确基因型类别图； 图3为真实数据结果图； 图4为核参数对精度影响图； 图5为模型学习迭代过程中相关向量的变化图，其中，(a)为拷贝数为2，覆盖度为5 时的相关向量变化图，(b)为拷贝数为3，覆盖度为10时的相关向量变化图，(c)为拷贝数为 4，覆盖度为15的相关向量变化图，(d)为拷贝数为5，覆盖度为20的相关向量变化图。
本发明公开了一种考虑拷贝数变异因素的基因组结构变异分型方法，输入序列比对文件和突变识别文件并统计记录各变异位点的特征值；根据输入文件提取特征值，从突变识别文件VCF中提取基因型作为分类监督，通过Python提取VCF文件中第八列type后的基因型信息，一行对应一个 全部
背景技术：
拷贝数变异(英文全称：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)存在不同标签基因型的特征数据相等或差别不大的情况，无法区分拷贝数变异 在杂合变异侧和纯合变异、正常无变异和拷贝数变异在杂合无变异侧等，加深了分型困难。 另外，现有方法需要大量的计算资源，很难应用于包含大量重复和片段重复的真 核基因组，不能处理杂合基因型。 5 CN 111583998 A 说 明 书 2/15 页
技术实现要素：
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种考虑拷贝 数变异因素的基因组结构变异分型方法，面向第二代基因测序数据，当基因组结构变异含 有大量拷贝数变异因素的情况下，使用机器学习策略实现基因分型的问题。 本发明采用以下技术方案： 一种考虑拷贝数变异因素的基因组结构变异分型方法，包括以下步骤： 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的行数，记为单端匹 6 CN 111583998 A 说 明 书 3/15 页 配特征值； 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具体为： 7 CN 111583998 A 说 明 书 4/15 页 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＝[k1,…,k Tn] ,K∈RN×N,为了确保模型的稀疏性，为权重向量 引入零均值，方差为 的标准正态先验分布，A由先验参数αnc组成的矩阵，αnc服从超参数 为τ，υ的Gamma分布；Y为辅助回归目标，c为目标类别，n为训练样本数量，i为类别标记， 为边缘似然函数。 与现有技术相比，本发明至少具有以下有益效果： 本发明一种考虑拷贝数变异因素的基因组结构变异分型方法，提出了将多分类相 关向量机模型应用于基因分型的方法，并将分类类别从纯合型杂合型两类扩充到了五类， 即不包含拷贝数变异的杂合型、不包含拷贝数变异的纯合型、包含拷贝数变异的杂合型、不 包含拷贝数变异的纯合型和无变异纯合型，解决了对包含拷贝数变异在内的多种复杂结构 变异的基因分型问题；融合了现有特征提取方法，并在综合利用现有特征信号的基础上，提 出15种全新有效的特征信号，包括异常读对、正常读对、未完全匹配读对、完全匹配读对、拆 分读对、单端读对、不能匹配读对、匹配质量、读对深度、加权读对深度、拓展的加权读对深 度、受影响读对、变异长度和读对的方向。避免了使用组装方法，节省了大量计算资源，能将 计算效率和结果精度平衡在高水平；利用期望最大化算法和快速Type-II最大似然来学习 和求解，降低了相关向量依赖，优化计算效率。 进一步的，能够识别对读对insert size造成影响的变异，进而区分出变异和未变 异，同时，拷贝数变异的存在会影响取值，根据数值大小，可以区分是否有拷贝数变异。 进一步的，在对五种基因型区分时，变异的存在会导致读对不能完美匹配到参考 序列上，读对与变异的位置关系、变异的类型等因素会造成读对完全、完全不能、不能完全 地与参考序列匹配，提取未完全匹配读对、完全匹配读对、拆分读对、单端读对和不能匹配 读对出每种情况对应的读对的数量，能够提示变异是否发生以及变异的类型，区别于现有 方法，将不同匹配情况的读对区分开能够扩大特征值的区分意义，实现高效精确的基因分 型。 进一步的，提取匹配质量能够反映读对与参考基因的匹配情况，通过匹配质量值 8 CN 111583998 A 说 明 书 5/15 页 的大小和全部变异区域匹配质量值计算，不同基因型和拷贝数变异会造成匹配质量值有显 著差异，这一特征值有利于基因分型。 进一步的，读对深度能反映基因正常、变异、拷贝数的差异，从而利于基因分型，通 过匹配质量对读对深度信号进行加权，能够进一步扩大三者差异；同时，扩展变异区间，能 够融合变异区间周围信息，且上下游各扩展一个读段的长度能够只融合可能被变异影响的 所有读段，避免过多正常区域的影响，从而提高特征值的分型能力，有助于高效基因分型。 进一步的，提取受影响读对能够反映变异区间，反映变异类型，有助于基因分型。 进一步的，取读取方向特征有助于判断倒置等结构变异，对基因分型显著有效。 进一步的，确定核函数和核函数参数，能够在线性核函数、多项式核函数、高斯核 函数、径向基核函数等中确定适合本方法应用问题的核函数和参数，并指导方法取得最好 结果。 进一步的，将数据分为M-RVM模型的训练集和测试集，利用训练集拟合模型，通过 设置分类器的参数，训练分类模型。通过训练集得出最优模型后，使用测试集进行模型预 测，以用来衡量该最优模型的性能和分类能力，进行模型性能评价。 综上所述，本方法全面理清了考虑拷贝数变异因素的基因组结构变异分型问题， 利用多分类相关向量机设计了一种高准确率、高效率的解法。 下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。 附图说明 图1为本发明流程图； 图2为正确基因型类别图； 图3为真实数据结果图； 图4为核参数对精度影响图； 图5为模型学习迭代过程中相关向量的变化图，其中，(a)为拷贝数为2，覆盖度为5 时的相关向量变化图，(b)为拷贝数为3，覆盖度为10时的相关向量变化图，(c)为拷贝数为 4，覆盖度为15的相关向量变化图，(d)为拷贝数为5，覆盖度为20的相关向量变化图。
