CN111462823B

CN111462823B - 一种基于dna测序数据的同源重组缺陷判定方法

Info

Publication number: CN111462823B
Application number: CN202010270712.0A
Authority: CN
Inventors: 赵仲孟; 戴道成; 易鑫; 易玉婷; 管彦芳; 王嘉寅; 张选平
Original assignee: Geneplus-Beijing; Xian Jiaotong University
Current assignee: Geneplus-Beijing; Xian Jiaotong University
Priority date: 2020-04-08
Filing date: 2020-04-08
Publication date: 2022-07-12
Anticipated expiration: 2040-04-08
Also published as: CN111462823A

Abstract

本发明公开了一种基于DNA测序数据的同源重组缺陷判定方法，获取特征属性；提取有效数据；基于三重学习法框架，考虑到较好的泛化能力、较高的准确度和对多维特征属性的处理效率，选择三个不同的基分类器H₁、H₂、H₃；对H₁、H₂、H₃进行迭代训练得到扩充训练集，由此对模型进行更新，完成训练过程；使用所训练的模型对未标记样本集U进行标记，根据标记结果完成HRD状态的判定。本发明解决了使用单一或少量基因组不稳定性状态等局部特征来进行HRD状态判定的局限性，克服临床上已知HRD状态的样本数量极少的难点，实现已有样本数据下的多特征属性的学习，能够提高HRD判定方法的性能。

Description

一种基于DNA测序数据的同源重组缺陷判定方法

技术领域

本发明属于以精准医学为应用背景的数据科学技术领域，具体涉及一种基于DNA测序数据的同源重组缺陷判定方法。

背景技术

同源重组(英文名称：Homologous Recombination，英文缩写：HR)是指发生在非姐妹染色单体之间或同一染色体上含有同源序列的DNA分子之间或分子之内的重新组合，其意义在于修复基因损伤，如图1所示。若参与同源重组的基因发生突变，理论上都有可能影响DNA的同源重组过程，继而发生同源重组缺陷(英文名称：Homologous RecombinationDeficiency，英文缩写：HRD)，又称HRD阳性。近年来，一些研究报道了HRD不仅可以用于癌症诊断，而且可以作为用药决策和患者预后的重要临床指标。因此，近年来，HRD是否呈阳性是肿瘤精准治疗中的重要指标之一，特别适用于卵巢癌、乳腺癌等典型同源重组缺陷型肿瘤的诊断分型。鉴于其重要的临床应用价值，判定基因组同源重组缺陷非常重要。

DNA测序技术作为分子生物学相关研究中最常用的技术手段之一。随着第二代测序技术(英文名称：Next Generation Sequencing，英文缩写：NGS)日益普及，临床实践普遍认为，利用癌症患者组织或血液样本的NGS测序数据分析和判断HRD阳性是一种经济和实用的方法。已有的方法研研究指出，从测序数据中检测获得的突变和基因组不稳定状态是鉴定肿瘤患者的HRD状态的关键数据信号。请参阅图2，主要有如下两类算法：

1)记录并分析基因组不稳定性状态，包括基因组杂合性缺失(英文名称：Loss ofheterozygosis，英文缩写：LOH)、端粒等位基因不平衡(英文名称：Teromeric-allelicimbalance,英文缩写：TAI)和大规模状态转移(英文名称：Large-Scale StateTransitions，英文缩写：LST)三项基因组不稳定性状态的值，再综合乳腺癌易感基因(英文名称：Breast Cancer Susceptibility Gene，英文缩写：BRCA)状态联合判定HRD的状态。此类方法目前以FoundationFocus CDx BRCA LOH检测法和MyChoice HRD检测法最为常见。

2)分析有关HRD突变信号，即根据突变过程留在肿瘤基因组上的特征模式的突变特征(signature)来判定HRD状态，或用专门设计的测序探针识别HR基因的点突变并进行分析。

目前，第1类方法使用比较普遍。第2类方法受癌种的局限性、突变特征基线、点突变集的覆盖度等因素影响，尚不成熟，仍处在科研摸索阶段。

虽然第1种方法已在临床诊断中应用，但在我国的临床实践中仍存在一些缺陷，主要问题是判定阈值难以确定和优化。Foundation Medicine公司的FoundationFocus CDxBRCA LOH检测法根据肿瘤基因组LOH值来判定HRD是否呈阳性，但是实验发现一些LOH低组的卵巢癌患者也是HRD阳性。Myraid公司的MyChoice HRD检测法根据高HRD评分和BRCA发生突变判定患者呈HRD阳性，低HRD评分和BRCA无突变判定患者呈HRD阴性，同时如果HRD评分失败且BRCA野生型，则HRD的状态无法确定。但是，关于HRD评分的特征中LOH和拷贝数的变化可能由非同源重组修复之间的染色体的重排引起的。综上所述，目前已有的两种方法不能完美解决利用基因组不稳定性状态等局部特征判定HRD状态的问题，存在较大的误差，效果亟待提高。

发明内容

本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于DNA测序数据的同源重组缺陷判定方法，将已有方法中的基因组不稳定性状态提取为一系列特征，利用已有数据库中少量的已明确获知HRD状态的样本，设计三重学习法框架；运用机器学习的结果优化判定阈值，即根据训练样本和待判定样本的共同特征优化判定阈值，解决固定阈值的局限性，实现更准确的判定HRD状态。

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

一种基于DNA测序数据的同源重组缺陷判定方法，包括以下步骤：

S1、对于测序数据下获得的样本序列比对文件，生成用于同源重组缺陷检测的特征及其属性文件，提取有效数据；

S2、基于三重训练法选择三个不同的监督学习算法构造对应的基分类器H₁、H₂、H₃，基于初始训练集训练基分类器H₁、H₂、H₃并生成对应的初始学习器；

S3、对步骤S2中的初始学习器进行测试训练，得到扩充训练集，使用扩充的训练集对学习器进行反复更新，直至满足迭代终止条件；

S4、对步骤S3完成训练的H₁、H₂、H₃中未标记的样本集U进行预测，采用软投票方式整合预测结果并对未标记样本赋予标记；根据每个个体全部样本的标记结果完成同源重组缺陷的判定：使用Python的机器学习库sklearn中的predict_proba()预测函数计算所预测类别的概率值，以三个基分类器的各自能力为权重，计算每个个体全部样本的权重和预测类别的概率值的乘积，获得其数学期望，选择数学期望值最高的结果作为个体的判定结果。

具体的，步骤S1具体为：

S101、读取每个个体的基因组测序数据，测序数据包含2个样本；将测序数据用生物信息学工具进行比对和变异检测，生成一组文件，每个文件包含1条染色体的单核苷酸突变的特征属性，用csv后缀表示，再生成另一组文件，每个文件包含1条染色体的等位基因特异性拷贝数的特征属性，用cncf.tsv后缀表示；

S102、采用数据归一化和数据关联分析的方法对步骤S101中生成的每个文件中的特征属性数据进行预处理。

进一步的，步骤S102中，使用统计软件对特征属性及其不同组合进行回归分析测试具体为：使用统计软件对cncf.tsv后缀文件中初步判断有关联性的13列特征属性进行步进回归分析，根据所得结果先后剔除chrom、seg、num.mark、segclust、cnlr.median.clust、cf.em共6类特征属性及其组合；再通过至少10次的预实验，最终选择nhet、cnlr.median、mafR、mafR.clust、start、end、cf.em、tcn.em、lcn.em共9类特征属性，用于构建三重学习法框架。

具体的，步骤S3中，更新训练集的方法具体为：

令初始标记样本集为L₀，样本数为|L₀|，按照三重训练法随机采样L₀中4/5的样本作为初始训练集L，初始训练集L的大小用|U|表示；设有三个基分类器H₁、H₂、H₃，L应满足H₁、H₂、H₃中任意1个所需的最小规模；首先，更新H₁训练集：令x为U中的任意一个元素，如果H₂、H₃对x的预测结果分别为H₂(x)、H₃(x)，且两者预测结果一致，将x的预测结果标记为H₂(x)并加入到H₁的训练集中；按此规则遍历U中的全部元素，形成H₁新的训练集L'₁；然后，更新H₂训练集：令x为U中的任意一个元素，如果H₁、H₃对x的预测结果分别为H₁(x)、H₃(x)，且两者预测结果一致，将x的预测结果标记为H₁(x)并加入到H₂的训练集中；按此规则遍历U中的全部元素，形成H₂新的训练集；更新H₃训练集；当H₁、H₂、H₃的训练集都被更新后，即认为完成1轮训练；接着，重新选取初始训练集L并重新训练和分别更新H₁、H₂、H₃的训练集；如此迭代直至满足迭代收敛条件：H₁、H₂、H₃的训练集与上一轮的训练集相同，则训练结束。

进一步的，为了优化步骤S3中选择初始测试集的比例，选择L₀中L的补集作为初始测试集，记为L_t；使用L_t对每1轮迭代的模型进行评估：比较第t轮训练时的三重学习法框架在L_t上的分类精确度得分和第t-1轮训练时的三重学习法框架在L_t上的分类精确度得分，若第t轮的得分大于第t-1轮的得分，则对训练集进行扩充。

更进一步的，步骤S3中，对训练集进行扩充的具体方法为：

判断由H₂、H₃为H₁新标记的训练集L^t ₁是否应该加入到H₁的新训练集中：令

和

分别表示第t、t-1次迭代的错误率上限且

令L^t和L^t-1分别表示t、t-1轮迭代中用于训练的标签样本集合，其中样本个数分别为|L^t|和|L^t-1|；当|L^t|＞|L^t-1|时，有

将H₂和H₃从U中选择的样本纳入H₁扩充训练集；此时，令

和

分别表示在第t、t-1轮中用于H₁训练的标签样本集合，H₁的第t、t-1轮的训练集则分别为

和

对应的样本个数分别用

表示；在第t-1轮结束后，被标记的未标记的样本不放进原始的L中；第t轮结束后，

中所有的无标签的样本再次被放进U中。

更进一步的，根据下式判定由H₂、H₃为H₁新标记的训练集是否应该加入到H₁的新训练集中：

其中，

是H₂、H₃第t次迭代的错误率的上限。

更进一步的，H₁的第t、t-1轮的训练集的规模应满足如下约束条件：

其中，η_L是初始训练集L的噪声率，

的样本个数的限制s的计算方法如下：

其中，s满足不等式

当

时，

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

本发明一种基于DNA测序数据的同源重组缺陷(HRD)的判定方法，采用基因组测序获得的个体的测序数据，通过生物信息学软件检测数据中包含的突变和染色体上等位基因特异性拷贝数等特征属性，用以解决传统方法使用单一或少量基因组不稳定性状态的局部特征来进行HRD判定的局限性。使用三重学习法框架，能够用少量的已知HRD状态的个体来预测较多的未知个体，一定程度上能够克服生产中HRD个体的数据不足的难点。此外，采用软投票法进行样本预测结果的整合和判定输出：在以上训练的三个基分类器对未知样本进行类别预测时，调用Python的机器学习库sklearn中的predict_proba()预测函数输出所预测类别的概率值，同时以三个基分类器的各自能力为权重，通过计算权重和预测类别的概率值的乘积的平均值作为标准，选择平均值较高的预测结果作为判定结果，加权平均概率投票法给予高自信的投票更大的权重，因此较硬投票法具有更好的判定结果。

进一步的，针对从基因组测序数据难以直接获取用于HRD判定的特征的问题，通过snp-pileup工具生成每条染色体上单核苷酸突变(英文名称：Single NucleotideVariant，英文缩写：SNV)的特征属性，用csv后缀文件记录；通过facet工具生成每条染色体上等位基因特异性拷贝数的特征属性，用cncf.tsv后缀文件记录。

进一步的，特征数据的归一化和筛选是一个重要的质控步骤，有助于消除异常值、无关特征对数据分析精度的影响。

进一步的，用于HRD判定问题的特征属性较多，且特征属性的学习过程是对于多属性二分类的过程，考虑到较好的泛化能力、较高的准确度和对多维特征属性的处理效率，本发明基于预实验结果，选取支持向量机、决策树和随机森林三个监督学习算法作为基分类器。

进一步的，在使用少量的标记数据辅助较多的未标记数据进行机器学习时，为了提高分类器的性能，必须对模型进行反复迭代与更新，通过比较错误率，能够确保每轮扩充的数据集满足评分改进的条件，能够降低分类的错误率。

进一步的，在迭代训练模型的过程中，当基分类器训练尚不充分时，基分类器有可能将错误的分类结果引入其余分类器中，由此产生噪声样本；对此，本发明设置了样本数量约束条件；当满足约束条件时，引入的噪声样本能够逐步被正确标记的训练样本抵消，使得基分类器趋于收敛。

综上所述，本发明属于一种基于机器学习框架的判定模型，具体的，设计和使用了一种三重学习法的机器学习框架，实现了从基因组测序数据中提取特征属性，选取了三个高效的基分类器，实现了训练模型的迭代的控制，使用软投票法改进了判定结果。基于此，本发明突破了使用少量基因组不稳定性状态等局部特征来进行HRD判定的局限性，解决了已知HRD状态的样本数量极少的难点。

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

附图说明

图1为同源重组示意图；

图2为检测HRD方法中的基因组不稳定性状态的示意图；

图3为本发明流程示意图。

具体实施方式

本发明提供了一种基于DNA测序数据的同源重组缺陷判定方法，获取特征属性；提取有效数据；基于三重学习法，选取较好的泛化能力、较高的准确度和对多维特征属性的处理效率的三个基分类器H₁、H₂、H₃；逐个对H₁、H₂、H₃进行训练得到扩充训练集，对基分类器进行迭代更新；使用训练完成的三重学习框架对未标记样本集U进行标记，根据标记结果完成HRD状态的判定。本发明突破了使用少量基因组不稳定性状态等局部特征来进行HRD判定的局限性，解决了已知HRD状态的样本数量极少的难点，实现了有少量数据时的通过多特征属性的机器学习来提高HRD判断方法的性能。

请参阅图3，本发明一种基于DNA测序数据的同源重组缺陷判定方法，包括以下步骤：

S1、对用于HRD检测的特征属性进行选择和处理；

对于测序下机数据，使用生物信息学工具进行比对，获得bam文件，再通过生物信息学工具生成用于HRD检测的特征属性文件并进行有效数据的提取处理，具体包含以下步骤：

S101、读取每个个体的基因组测序数据，测序数据应包含2个样本，其中1个样本来自个体的肿瘤组织，另1个样本来自个体的癌旁对照组织；将测序数据用生物信息学BWA工具进行比对；通过snp-pileup工具完成变异检测，生成一组文件，每个文件包含1条染色体的单核苷酸突变的特征属性，用csv后缀表示；通过facet工具完成拷贝数变异检测，再生成一组文件，每个文件包含1条染色体的等位基因特异性拷贝数的特征属性，用cncf.tsv后缀表示；使用snp-pileup工具应配置的功能和参数如下：

snp-pileup-g-p-q20-Q20-P100-r25,20common_all_20180423.vcf.gzoutputfile normalbam tumorbam

其中，-g和-p为软件运行时的输出信息；考虑到捕获测序的深度为50×，按照每100个碱基位点(-P：100)统计SNV变异情况并进行输出；碱基质量值(Quality Score或Q-score)是碱基识别出错的概率的整数映射，其中Q＝-10*lgP，其中P为碱基识别出错的概率；当-q选用20，-Q选用20时，P＝1％，即选用99％的比对质量和测序质量进行质控，过滤掉质量低于1％的SNV；-r输出仅在足够深度处所捕获的SNV信息，对于50×的捕获测序深度，参照软件手册指定癌旁对照组织的最小测序深度的默认值为25×，考虑到肿瘤基因组中的特点，指定肿瘤组织的最小测序深度的默认值为20；common_all_20180423.vcf.gz是一个用于校准测序数据的vcf文件；outputfile是输出文件，以csv为文件后缀；normalbam、tumorbam分别表示来自个体的癌旁对照组织、肿瘤组织的bam文件。

调用facet工具对以上所得的csv文件进行进一步分析和处理，生成包含每条染色体上等位基因特异性拷贝数估计等特征属性的cncf.tsv后缀文件。使用facet工具应配置的功能和参数如下：

facet-s*.csv.gz-p prefix-o.

其中，-s选项为配对的SNV信息的csv格式的文件；-p为输出文件的前缀名称；-o为输出文件的目录。cncf.tsv文件包含13项特征属性，其物理意义如表1所示：

表1cncf.tsv文件包含的特征属性及其物理意义

S102、读取cncf.tsv文件中各个属性的数据，采用数据归一化和数据关联分析的方法对上述文件中的数据进行预处理。

由于无关特征属性或在多维特征属性的回归分析中占据极小比重的特征属性可能干扰机器学习过程，因此，考虑剔除这些特征属性。使用统计分析软件SPSS对表1中的特征属性及其不同组合进行回归分析测试。考虑多元线性回归模型，根据回归模型对自变量赋予的权重，设置权重阈值。随机抽取一部分训练样本，进行不少于10次预实验，低于阈值的自变量chrom、seg、num.mark、segclust、cnlr.median.clust、cf.em六类属性及组合属性予以排除。最终选择nhet、cnlr.median、mafR、mafR.clust、start、end、cf.em、tcn.em、lcn.em共9类特征属性，用于构建三重学习法框架。

在此基础上，采用线性归一化的方法进行数据归一化操作，即将相差较大的属性数据按比例缩放，使其落入一个较小的特定区间：

其中，x_min为本属性列的所有数据的最小值，x_max为本属性列的所有数据的最大值，x为任意一个属性数据。

通过上述操作，可以消除奇异样本数据(指相对于其他输入样本特别大或特别小的样本矢量)对机器学习的不良影响，可以提升算法的收敛速度和精度。

S2、基于三重训练法(英文名称：Tri-training)设计思想，令初始标记样本集为L₀，样本数为|L₀|，对L₀采用随机自助法(bootstrap)生成有差异的三份训练数据子集，然后使用不同的监督学习算法来构造三个不同的初始分类器H₁、H₂、H₃。考虑到较好的泛化能力、较高的准确度和对多维特征属性的处理效率，本发明基于预实验结果，选取支持向量机、决策树和随机森林三个监督学习算法作为基分类器。差异化的三个基分类器的生成有助于得到具有较高泛化能力的集成分类器。

S3、对步骤S2中的H₁、H₂、H₃进行训练得到扩充训练集，由此对模型进行更新，完成训练过程；

更新训练集的方法如下：

按照三重训练法随机采样L₀中4/5的样本作为初始训练集L，初始训练集L的大小用|U|表示；设有三个基分类器H₁、H₂、H₃，L应满足H₁、H₂、H₃中任意1个所需的最小规模；首先，更新的H₁训练集：令x为U中的任意一个元素(即样本)，如果H₂、H₃对x的预测结果分别为H₂(x)、H₃(x)，且两者预测结果一致，那么将x的预测结果标记为H₂(x)并加入到H₁的训练集中；按此规则遍历U中的全部元素，形成H₁新的训练集L'₁：

L'₁＝L∪{x|x∈U,H₂(x)＝H₃(x)}

然后，更新的H₂训练集：令x为U中的任意一个元素(即样本)，如果H₁、H₃对x的预测结果分别为H₁(x)、H₃(x)，且两者预测结果一致，那么将x的预测结果标记为H₁(x)并加入到H₂的训练集中；按此规则遍历U中的全部元素，形成H₂新的训练集；同理，更新的H₃训练集；当H₁、H₂、H₃的训练集都被更新后，即认为完成1轮训练；接着，重新选取初始训练集L并重新训练和分别更新H₁、H₂、H₃的训练集；如此迭代直至满足迭代收敛条件：H₁、H₂、H₃的训练集与上一轮的训练集相同，则训练结束。

选择L₀中L的补集作为初始测试集，记为L_t；使用L_t对每1轮迭代的模型进行评估：比较第t轮训练时的三重学习法框架在L_t上的分类精确度得分和第t-1轮训练时的三重学习法框架在L_t上的分类精确度得分，若第t轮的得分大于第t-1轮的得分，则可以对训练集进行扩充。

判断训练集是否可以扩充的步骤具体为：

如果H₂和H₃对x的预测正确，则H₁得到一个新增的有效样本x，x被扩充进入训练集；否则H₁得到一个有噪音标签的样本。

对分类器的训练集进行扩充时，需满足至少有m个样本满足约束条件，样本数量m应满足的约束条件为：

其中，ε是假设的最坏情况下的分类错误率，η是分类噪声的上界，默认值是0.5，即允许某次分类预测错误的样本数量不超过本次实验待判定样本数量的一半，N是本次实验待判定样本的数量，δ是置信度。

判断于任意一个训练样本是否可以扩充进入训练集的步骤具体为：

在模型的训练过程中，判定由H₂、H₃为H₁新标记的训练集

是否应该加入到H₁的训练集中。

令

和

分别表示第t、t-1次迭代的错误率上限且

因此，根据下式判定由H₂、H₃为H₁新标记的训练样本集是否应该加入到H₁的训练集中：

L_n'中的样本被分类器预测的结果只有正确和错误两类：若预测正确，其对于被训练分类器而言是增加了一个正确的训练样本；若预测错误，其对于被训练分类器而言则是增加了一个噪声。为了提升分类器的分类能力，对于基分类器训练尚不充分时所产生的错误分类结果，即模型噪声，必须进行过滤处理。

在某轮训练中，H₂和H₃从U中选择的样本应纳入H₁扩充训练集；此时，令

和

和

对应的样本个数分别用

表示；在第t-1轮结束后，被标记的未标记的样本就不再放进原始的L中；第t轮结束后，

中所有的无标签的样本再次被放进U中。同理，可以判断由H₁、H₃为H₂新标记的训练集是否应该加入到H₂的新训练集中，以及由H₁、H₂为H₃新标记的训练集是否应该加入到H₃的新训练集中。

H₁的第t次和第t-1次的训练集的满足约束条件如下：

其中，η_L是初始训练集L的噪声率，其他变量如上所定义。

为了控制每轮标记的样本数；

样本个数的限制s的计算如下：

其中，s满足不等式

在此基础上，为了满足

则需满足条件：

基于机器学习模型的噪声理论，当

样本个数的限制s和

成立，噪声数据所带来的分类错误便能够逐步被正确标记的训练集所抵消，使得机器学习框架在迭代过程中保证了分类器的分类错误越来越少，从而趋于收敛。

按照以上步骤对三重学习法框架进行重复迭代，直至满足迭代收敛条件：H₁、H₂、H₃的训练集与上一轮的训练集相同，则训练结束。

S4、使用步骤S3中训练的模型对未标记样本集U进行标记，根据标记结果完成HRD状态的判定。

将上述迭代后的三个稳定的分类器集成后对未标记的样本集进行标记，根据较高性能的软投票的标记结果完成HRD状态的判定。

采用软投票方式整合预测结果并对未标记样本赋予标记；根据每个个体的全部样本的标记结果完成同源重组缺陷的判定：使用Python的机器学习库sklearn中的predict_proba()预测函数计算所预测类别的概率值，以三个基分类器的各自能力为权重，计算每个个体的全部样本的权重和预测类别的概率值的乘积，获得其数学期望，选择数学期望值最高的结果作为个体的判定结果。

对于标注集，输出样本中每条染色体的标记结果，并统计该算法对每条染色体的预测结果为阴性n₀和预测结果为阳性n₁的总数，当n₀≥n₁时，该样本的HRD判定结果Pred_Type归为阴性；当n₀<n₁时，该样本的HRD判定结果Pred_Type归为阳性。再对照原有标注集中HRD的状态，计算评价指标TP、TN、FP、FN、FPR、FNR、PPV、NPV、Accuracy、Sensitivity、Specificity的值。对于未标注集，输出样本中每条染色体的预测结果，其中n₀、n₁和Pred_Type的统计方法如上。

以上评价指标的定义分别是：

TP：true positive，原本为HRD阳性的样本经过本判定方法的检测结果也为阳性。

FP：false positive，原本为HRD阴性的样本经过本判定方法的检测结果为阳性。

TN：true negative，原本为HRD阴性的样本经过本判定方法的检测结果也为阴性。

FN：false negative，原本为HRD阳性的样本经过本判定方法的检测结果为阴性。

FPR＝FP/(FP+TN)，FP rate，原本为HRD阴性的样本经过本判定方法的检测结果为阳性的比例。

FNR＝FN/(FN+TP)，FN rate，原本为HRD阳性的样本经过本判定方法的检测结果为阴性的比例。

PPV＝TP/(TP+FP)，positive predict value，HRD样本经过本判定方法的检测结果为阳性的正确比例。

NPV＝TN/(TN+FN)，negative predict value，HRD样本经过本判定方法的检测结果为阴性的正确比例。

Accuracy＝(TP+TN)/(TP+TN+FN+FP)，HRD样本经过本判定方法的检测结果与其原本状态一致的比例。通常来说，此比例越高，准确率越高，即本判定方法的检测效果越好。

Sensitivity＝TP/(TP+FN)，原本为HRD阳性的样本经过本判定方法的检测结果为正确的比例。通常来说，此比例越高，本判定方法对HRD阳性样本的检测效果越好。

Specificity＝TN/(TN+FP)，原本为HRD阴性的样本经过本判定方法的检测结果为正确的比例。通常来说，此比例越高，本判定方法对HRD阴性样本的检测效果越好。

根据标注集的结果中precision、recall、f1-score值和其微均值、宏均值及加权均值的最优值选择相应最精确的预测结果。

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

收集了两组实验数据，一组是公开的临床数据库中经过人工筛选的已知HRD状态的56例肿瘤样本，其中卵巢癌7例、乳腺癌21例、肺腺癌28例，用于检验本发明所提方法的性能。另外一组是公开的临床数据库中已有的未知HRD状态的188例肿瘤样本，其中卵巢癌48例、乳腺癌112例、肺腺癌28例，用于实测判定结果。取标注集和为未标注集的比例为1:4，以及二分类问题中两类别的比例1:1。选取已知HRD状态的56例肿瘤样本中2例卵巢癌样本和4例乳腺癌样本作为标注集的HRD阳性样本，选取6例肺腺癌样本作为标注集的HRD阴性样本，同时选取剩余的5例卵巢癌样本、17例乳腺癌样本和22例肺腺癌样本作为未标注集的HRD状态待判定样本。

按照本发明方法的步骤，首先根据56例样本的bam文件，通过snp-pileup工具、facet工具生成cncf.tsv文件。选择cncf.tsv文件中的所有属性列作为HRD的特征属性时，应用本发明方法对未标注集进行检测实验，10次实验结果如下：

表2本发明方法应用13列特征属性对44例未标注集检测的10次实验结果

使用SPSS工具对特征属性进行回归分析，剔除chrom、seg、num.mark、segclust、cnlr.median.clust特征属性后，使用完全相同的步骤进行检测实验，10次实验结果如下：

表3本发明方法应用9列特征属性对44例未标注集检测的10次实验结果

结合表2和表3可见，去除了chrom、seg、num.mark、segclust、cnlr.median.clust特征属性后，44例未标注集检测的平均精度有所上升，说明本发明筛选特征属性，选择nhet、cnlr.median、mafR、mafR.clust、start、end、cf.em、tcn.em、lcn.em特征属性是有益的。

选取12例标注集的HRD样本和44例未标注集的HRD样本，使用本发明方法对44例未标注集的HRD样本进行检测实验，15次的实验结果如下：

表4本发明对44例未标注集的检测结果

从表4可见，本发明方法对44例未标注集的检测精确度平均达到97.73％，平均检测精确度达到99.39％。

选取目前已有的48例卵巢癌样本，112例乳腺癌样本，28例肺腺癌样本通过本发明方法进行HRD状态判定，得到188例样本的HRD的3次判定结果。按照本发明方法中标注集和未标注集的比例，本次预测实验的标注集选用以上检测实验中的部分标注集数据。

表5本发明方法对188例样本HRD的3次预测结果

3次重复实验下，57例卵巢癌样本中HRD的阳性比例为59.65％，131例乳腺癌样本中HRD的阳性比例为71.76％。这与癌症和肿瘤基因图谱(英文名称：The Cancer GenomeAtlas，英文缩写：TCGA)项目研究发现超过50％高级别浆液性卵巢癌存在HRD，以及与GeparSixto研究结果表明三阴性乳腺癌患者中同源重组缺陷率为70.5％的结论基本一致。

综上所述，本发明一种基于DNA测序数据的同源重组缺陷判定方法在各项指标上优于目前国内外的HRD判定方法。

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

Claims

1.一种基于DNA测序数据的同源重组缺陷判定方法，其特征在于，包括以下步骤：

S1、对于测序数据下获得的样本序列比对文件，生成用于同源重组缺陷检测的特征及其属性文件，提取有效数据，具体为：

S102、采用数据归一化和数据关联分析的方法对步骤S101中生成的每个文件中的特征属性数据进行预处理；

S3、对步骤S2中的初始学习器进行测试训练，得到扩充训练集，使用扩充的训练集对学习器进行反复更新，直至满足迭代终止条件，更新训练集的方法具体为：

令初始标记样本集为L₀，样本数为|L₀|，按照三重训练法随机采样L₀中4/5的样本作为初始训练集L，初始训练集L的大小用|U|表示；设有三个基分类器H₁、H₂、H₃，L应满足H₁、H₂、H₃中任意1个所需的最小规模；首先，更新H₁训练集：令x为U中的任意一个元素，如果H₂、H₃对x的预测结果分别为H₂(x)、H₃(x)，且两者预测结果一致，将x的预测结果标记为H₂(x)并加入到H₁的训练集中；按此规则遍历U中的全部元素，形成H₁新的训练集L'₁；然后，更新H₂训练集：令x为U中的任意一个元素，如果H₁、H₃对x的预测结果分别为H₁(x)、H₃(x)，且两者预测结果一致，将x的预测结果标记为H₁(x)并加入到H₂的训练集中；按此规则遍历U中的全部元素，形成H₂新的训练集；更新H₃训练集；当H₁、H₂、H₃的训练集都被更新后，即认为完成1轮训练；接着，重新选取初始训练集L并重新训练和分别更新H₁、H₂、H₃的训练集；如此迭代直至满足迭代收敛条件：H₁、H₂、H₃的训练集与上一轮的训练集相同，则训练结束；

为了优化步骤S3中选择初始测试集的比例，选择L₀中L的补集作为初始测试集，记为L_t；使用L_t对每1轮迭代的模型进行评估：比较第t轮训练时的三重学习法框架在L_t上的分类精确度得分和第t-1轮训练时的三重学习法框架在L_t上的分类精确度得分，若第t轮的得分大于第t-1轮的得分，则对训练集进行扩充，对训练集进行扩充的具体方法为：

判断由H₂、H₃为H₁新标记的训练集

是否应该加入到H₁的新训练集中：令

和

分别表示第t、t-1次迭代的错误率上限且

将H₂和H₃从U中选择的样本纳入H₁扩充训练集；此时，令

和

和

对应的样本个数分别用

中所有的无标签的样本再次被放进U中；

2.根据权利要求1所述的基于DNA测序数据的同源重组缺陷判定方法，其特征在于，步骤S102中，使用统计软件对特征属性及其不同组合进行回归分析测试具体为：使用统计软件对cncf.tsv后缀文件中初步判断有关联性的13列特征属性进行步进回归分析，根据所得结果先后剔除chrom、seg、num.mark、segclust、cnlr.median.clust、cf.em共6类特征属性及其组合；再通过至少10次的预实验，最终选择nhet、cnlr.median、mafR、mafR.clust、start、end、cf.em、tcn.em、lcn.em共9类特征属性，用于构建三重学习法框架。

3.根据权利要求1所述的基于DNA测序数据的同源重组缺陷判定方法，其特征在于，步骤S3中，根据下式判定由H₂、H₃为H₁新标记的训练集是否应该加入到H₁的新训练集中：

其中，

是H₂、H₃第t次迭代的错误率的上限。

4.根据权利要求3所述的基于DNA测序数据的同源重组缺陷判定方法，其特征在于，H₁的第t、t-1轮的训练集的规模应满足如下约束条件：

其中，η_L是初始训练集L的噪声率，

的样本个数的限制s的计算方法如下：

其中，s满足不等式

当

时，