CN109493917A - 一种基因突变有害性预测值的害阶位计算方法 - Google Patents

Abstract

本发明是关于基因突变有害性的软件预测结果的标准化处理的算法。该算法基于分位数思想，利用基因突变的诸预测软件有害性预测值，计算基因突变各软件的害阶位。所谓害阶位，指某突变的某软件有害性预测值在样品测得所有突变的该软件有害性预测值中的排列位置，反映了比该突变更有害或更无害的突变的比例。该算法有利于解决各预测软件的取值范围、取值方向、致病与否的阈值不统一的问题，使各个软件的结果可以换算为统一的取值模式。在遗传病高通量测序分析中，需要判断大量突变的有害性。该算法可以准确地和充分地利用预测软件的预测，判断突变的有害性，实现自动化地提取致病突变。因此，害阶位概念和算法可成为遗传病高通量测序分析的有效工具。

Description

一种基因突变有害性预测值的害阶位计算方法

技术领域

本发明是解决基因突变有害性的软件预测结果的标准化处理问题，属于医学中遗传病基因突变检测领域。

技术背景

当前以全外显子测序为代表的高通量测序是人类遗传病检测和诊断的主要手段。高通量测序可以发现海量的变异，判断这些变异的有害性，找出和确认导致患者的致病基因和致病位点，是遗传病基因检测的主要任务。

在测序得到的变异中，有些变异在人群中有较高的分布频率。较高的分布频率意味着，如果该变异是有害的和不利于生存的，在显性遗传的情况下携带该突变的个体将趋于淘汰；在隐性遗传的情况下很容易与同一变异携带者婚配，从而生下纯合型携带者而趋于淘汰，自然地携带该突变的个体数量将趋于减少，因此高频率一般代表无害。在对基因变异进行注释时，根据变异在千人基因组数据库和其它同类数据库中的频率，就可以将大部分变异从致病变异候选者中排除出去。

在较低频率的变异中，有些变异属于移码突变，可导致其后所编码的肽链完全改变且在多数情况下很快中断；或者突变产生一个终止密码子引起所编码的肽链截断，除个别蛋白质外，这通常会导致蛋白质的失活，故该类型突变在大部分的已知机理下是有害的。因此在多数情况下，包括移码突变在内的截断突变是可快速判断有害性的突变类型。

有些突变发生在外显子前面或后面的内含子区域一两个碱基的位置。这些位置是与转录剪切有关的位置，该位置的突变将改变蛋白质序列，通常也是有害的，是较易判断的突变类型。

还有一大类分布于内含子或外显子上的碱基置换型的变异，或者插入或缺失但不引起移码突变的变异，这些变异中的一部分已被深入地研究过，对突变有害性给出判断。在引用这些已经研究过的突变时，既可以检索文献，也可以查阅或通过程序自动调取收录突变致病性的数据库，如 HGMD , ClinVar , dbSNP , Ensembl 数据库等。

然而还有大量非同义突变(即使同义突变，也不一定是无害的)，在人群中分布频率较低，又无文献支持，其致病性是难以判断的。在未经实验研究的情况下为了初步判断突变的致病性，一些基因突变致病性预测软件应运而生，从比较早期的SIFT，PolyPhen，MutationTaster，到最近几年发展起来的LRT,FATHMM,PROVEAN,VEST3,fathmm-MKL,MetaSVM,MetaLR等等，这个领域的软件越来越多了。

这些软件有些有自己的网站，根据要求的格式输入突变，经过一段时间的等待，可以给出预测的结果。有些软件没有网站，需要下载安装使用。在进行大规模突变的预测运算时，大多数情况都需要服务器，普通用户使用这些软件是有一定困难的。

因此，Xiaoming Liu等人开发了dbNSFP数据库(Liu X,Jian X,BoerwinkleE.dbNSFP:a lightweight database of human nonsynonymous SNPs and theirfunctional predictions.Hum Mutat.2011 Aug；32(8):894-9.)。他们整合了多种基因突变有害性预测软件，收录大量的人类变异，用这些预测软件计算所收录的突变并尽可能给出每个突变每种软件的预测值，进而构建 dbNSFP数据库。用户可以下载该数据库，在样本进行高通量测序并拿到VCF格式的数据后，使用dbNSFP数据库对突变进行有害性预测的注释。

然而每种软件的算法原理是不同的，所给出的取值范围是不同的，取值有害性的方向是不同的(有些软件值越大越有害，有些正相反)，是否致病的阈值设定是不同的；甚至有些软件的计算值不能用于判断致病性，而是另外给出了判断结果。这些都给用户的使用带来不便。

为了解决这个问题，dbNSFP数据库于2016年使用了评级值(rank score)算法(LiuX,Wu C,Li C,Boerwinkle E.dbNSFP v3.0:A One-Stop Database of FunctionalPredictions and Annotations for Human Nonsynonymous and Splice-Site SNVs.HumMutat.2016Mar；37(3):235-41.)。评级值是将dbNSFP收录的所有突变，按照某个软件的预测值排序；某个突变在整个数据库中的位置即为评级值。例如某突变的评级值是0.9，代表该突变在该软件的预测算法中比90％的突变更有害。由于每个软件的预测值几乎都可以通过这一方法排序并给出评级值，这就为将不同软件的预测值用统一的标准进行比较和计算提供了可能。

这个算法有几个局限，一是使用者必须使用dbNSFP数据库，二是突变必须是dbNSFP 数据库已经收录并计算了的，三是突变的预测软件必须是dbNSFP收录的。又由于dbNSFP 数据库还在不断收录新的突变，这些新收录的突变绝大部分是稀有突变，稀有突变的加入是否影响原有突变在某软件中的评级值，也是不得而知的。

鉴于此，我们开发了不依赖于大型数据库的基因突变有害性预测值的害阶位的计算方法，并希望在此基础上优化致病突变的筛选策略。

附图说明

图1软件害预值与害阶位关系散点图1

横坐标为害阶位，纵坐标为各软件的害预值，共包含了5款软件。

图2软件害预值与害阶位关系散点图2

横坐标为害阶位，纵坐标为各软件的害预值，共包含了6款软件。

图3分样本计算害预值害阶位曲线的重合情况

横坐标为害阶位，纵坐标为该软件的害预值。曲线中包含了20个单独样本和一个总体样本的结果曲线。

图4不同百分位分位值的标准差曲线图

横坐标为害阶位的序号，从1开始到100结束，纵坐标为害阶位的标准差的值。

图5害阶位计算的近似值求取流程

图6害阶位计算的精确值求取流程

图7fathmm-MKL-coding软件害阶位与评级值关系散点图

横坐标为评级值，纵坐标为害阶位。

图8PROVEN软件害阶位与评级值关系散点图

横坐标为评级值，纵坐标为害阶位。

图9有害突变集中度与害阶位倒数的关系

横坐标为对集中度的不同算法，纵坐标为集中度的值。

图10害阶位在遗传病高通量测序数据分析中的应用

第一行为excel表格中的位置编号，从AM-BJ，所对应的第二行的题头依次是：

预测的和，SIFT_score，SIFT_pred，Polyphen2_HDIV_score，Polyphen2_HDIV_pred， Polyphen2_HVAR_rankscore，Polyphen2_HVAR_pred，LRT_score，LRT_pred，MutationTaster_score，MutationTaster_pred，MutationAssessor_score_rankscore，MutationAssessor_pred，FATHMM_score，FATHMM_pred，PROVEAN_score，PROVEAN_pred，MetaSVM_score，MetaSVM_pred，MetaLR_score，MetaLR_pred，fathmm-MKL_coding_score，fathmm-MKL_pred

具体实施方式

具体实施方式，首先依赖于对分位数思想的理解。

分位数思想

在包含一定数量个体的样本中，当样本某个属性可以用不同的测量工具测量时，则不同个体的该属性的排序，不依测量工具的不同而不同。简单举例，一群人的身高可以用英尺测量，也可以用中国的尺寸测量，用两种测量体系分别准确测量后对身高数据排序得到两个序列，则任一个体在两个序列中将有着相同的身高排名。来源于同一群样本中的两个个体，如果使用了不同的测量方法，分别得到不同的测量数值，而这两种数值又不能或不方便直接换算，则为了比较这两个个体数值的大小，可以各自分别换算成它在样本中的排序再进行比较。由于排序的序号值是随着样本数量的变化而变动的，为了解决这一问题，引入了分位数思想。即将样本平均分成多少份，以便将某个个体的数值归类到相应的区间。

分位数，亦称分位值、分位点，是指将一个随机变量的概率分布范围分为几个等份的数值点，常用的有中位数(中值)、四分位数、百分位数等。其定义为:当随机变量X的概率密度函数为Pr(X)，若变量x满足

Pr(X<x)＝k/q，

则，x为X的第k个q分位数对应的值，记为x_k。对于有限离散样本空间，样本数量为N，将随机变量X按照从小到大的顺序依次标记为X₁，X₂，…，X_N，则x_k可以通过I_k＝N_k /q计算而得，即寻找X₁，X₂，…，X_N中对应的第I_k个值。通常，如果q＝4，称为四分位数；如果q＝100，称为百分位数(孙洪泉，吕娟，苏志诚,等.分位数法对多指标干旱等级划分一致性的作用[J].灾害学，2017，32(2):13-17，53.)。分位数思想在统计学领域，如测量身高(马军,吴双胜,宋逸,胡佩瑾,张兵，1985～2005年中国7～18岁学生身高、体重变化趋势分析，北京大学学报(医学版)，Vol.42 No.3 Jun.2010：318-322)，收入调查(李培林，朱迪：努力形成橄榄型分配格局——基于2006-2013年中国社会状况调查数据的分析，中国社会科学，2015年第一期：45-65)等领域有广泛用途。

在使用高通量测序进行遗传病检测时，所得到的大量突变经多种致病性预测软件计算可得到预测值，由于这些软件使用了不同的算法原理，不同的取值范围、不同的取值方向，使得不同软件对同一突变的预测的有害性很难比较。根据分位数的思想，我们希望将这些预测数值转化为可相互比较的分位数。其原理是，来自样本的一定数目的基因突变，经预测软件预测所得的突变有害性预测值(简称：害预值)，从有害性向可耐受性的方向排列编秩，某一突变害预值X在害预值编秩中的位置，可用该害预值的序号值与编秩中最大序号值的比值表示，取值范围在0-1之间；我们称该值为害预值X对应的害阶位Px；当害预值存在大量重复，使得某一害预值对应多个可能的害阶位时，则取害阶位最大值。因此，某一突变A的害预值X的害阶位Px的意义在于：较之突变A，有Px％的突变具有同样严重或更严重的有害性，有(1-Px)％的突变具有较轻的有害性。

同一害阶位对应多个害预值的可能性是极小的，但也有这种情况。某些预测软件的预测极大值与(或)极小值未设定，在有限样本获得害预值极值不足以代表理论极值时，可以对害阶位的极值(0或1)同时给予实际的害预值和理论的害预值(这两个值可以有很大的差别)，从而在以后的检测中，当害预值落到实际的害预值极值之外，理论的害预值极值之内时，仍可以计算得到正确的害阶位。在这种情况下，对害阶位的极值(0或1)是可以对应两个或多个害预值的。

根据以上原理，为了建立从基因突变的害预值计算害阶位，以及从害阶位计算害预值的方法，并考察这些方法的可靠性，我们循着以下步骤和逻辑展开工作：(1)获取可用于分析的数据； (2)展示不同软件判断致病性的差异；(3)获取害预值-害阶位表格(害阶尺)的算法；(4) 评估根据不同样本来源的数据计算得到的害阶位的稳定性；(5)根据害阶尺从害预值计算害阶位的算法；(6)根据害阶尺从害阶位计算各软件害预值的算法；(7)害阶位与评级值的比较；(8)集中度的计算；(9)害阶位与HGMD，ClinVar数据库注释的比较分析；(10)害阶位倒数与HGMD等数据库注释的比较分析；(11)“损低固高法”集中有害突变；(12)害阶位概念在二代测序分析中的应用。

1获取可用于分析的数据

目的：获取可用于分析的数据。

材料：我们用于分析的数据来自20个病人的全外显子二代测序的数据。这些病人或多或少存在遗传方面的疾病。每个病人经过全外显子二代测序后，对突变进行初步分析。凡是位于内含子区域，且在内含子-外显子交界处从外显子方向数起，10个碱基之外的突变，不予分析。这样每个样品的全部突变在23000个左右。

方法：使用dbNSFP数据库对这些突变进行有害性预测注释，有些突变在数据库中无记录，故没有注释，有些突变在数据库中有一个软件的注释，有些则有多个软件的注释。对这些突变进行分类，分别统计不同类型的突变的注释情况。应注意的是，这些突变同时进行了其它注释，如突变在HGMD数据库中是否有收录，HGMD标注的致病性情况等。

结果：我们发现每个样品至少有一个有害性预测注释的突变数量在9000左右，主要是错义突变。而完全没有有害性预测注释的突变，每个样品大约有14000个，主要是内含子突变和同义突变。详细结果统计如下：

表1，至少注释了一条害预值的突变情况一览

表2，没有任何害预值的突变情况一览

以上表格中的突变类型的中文名称-英文名称对应关系如下：

注释时名称	中文名称
		frameshift deletion	移码缺失
frameshift insertion	移码插入
		nonframeshift deletion	缺失
nonframeshift insertion	插入
		nonsynonymous SNV	错义突变
stopgain	终止
		stoploss	延长
synonymous SNV	同义
		unknown	未知
-	内含子

讨论：我们在进行害阶位分析时，只分析至少有一条害预值注释的突变，而忽略没有任何有害性预测注释的突变。

在我们的数据中，对突变注释的软件罗列如下：

SIFT

Polyphen2_HDIV

Polyphen2_HVAR

LRT

MutationTaster

MutationAssessor

FATHMM

PROVEAN

MetaSVM

MetaLR

fathmm-MKL coding

结论：我们获得了用于分析的数据。

2展示不同软件判断致病性的差异

目的：展示不同软件判断致病性的差异，以说明计算害阶位的必要性。

材料：已经获取的20个样本的全外测序结果(未去除不同样本的相同突变)。

方法：对不同软件的预测区域、突变类型、预测原理、取值范围，临界值等进行基于文献的比较和总结。

结果：不同的基因突变有害性预测软件，存在着极大的差异，这些差异包括：

(1)可预测区域的差异和可预测突变类型的差异。

从功能区域角度看，所有的预测软件所能给出预测结果的突变都主要集中在外显子区域，且主要为错义突变，大部分软件都或多或少地可以预测发生在内含子和剪切敏感区域的突变，其中MutationTaster，fathmm-MKL软件的功能略强。对存在于内含子和剪切敏感区域的所有突变而言，这些软件所分析的突变只占其中极少一部分。在20个样品的突变中，凡属含有至少一个软件注释的突变，我们统计了各功能区域含有预测结果的突变数目，结果见表“分软件分功能区域的突变统计”。

从突变类型角度看，所有软件都主要分析错义突变，但有些软件如MutationTaster， fathmm-MKL等软件也分析相当数量的剪切突变和同义突变。在20个样品的突变中，凡属含有至少一个软件注释的突变，我们统计了各突变类型的软件预测数目，结果见表“分软件分突变类型的突变统计”。从表中看，各软件所能分析的突变的总数也有一定差异。

表3，分软件分功能区域的突变统计

表4，分软件分突变类型的突变统计

(2)预测软件算法原理的差异。

不同的预测软件训练使用的数据来源不同，所用数据特征不同，预测模型也不相同，罗列如下：

表5，各预测软件算法原理一览

(文献：Liu X,Wu C,Li C,Boerwinkle E.dbNSFP v3.0:A One-Stop Database ofFunctional Predictions and Annotations for Human Nonsynonymous and Splice-Site SNVs.Hum Mutat.2016 Mar；37(3):235-41.)

(3)致病程度分类、取值范围、方向及有害性临界值不同的预测软件不但算法模型不同，而且最终给出的预测结果的致病程度分类、取值范围、取值方向，区分致病性与否的临界值等也不尽相同。例如大部分软件的取值范围是0和1 之间，但另一些软件的取值范围则否，如Mutation Assessor的取值范围是[- 5.545,5.937]；还有一些软件未给出明确的取值范围界限，如PROVEAN，我们只能根据该软件大量的预测数值判断出该软件的大致取值范围。不同软件的取值与有害性的对应关系也是不同的。如SIFT值越小越有害，MutationAssessor值越大越有害，还有 MutationTaster2和LRT的取值与致病性分类之间并没有完全的对应关系。不同软件的致病性划分方法也是不同的，有的是二分法，基本分为有害和无害，如SIFT，FATHMM 等，有些则复杂一些，是三分法或四分法，如Mutation Assessor，Polyphen2_HDIV， Polyphen2_HVAR，LRT是三分法，FATHMM是四分法。自然的，不同软件在进行突变致病性分类时，如果直接根据评分大小划分，那么它们的用以分类的临界值是各自不同的。

表6，各预测软件致病性分类一览

(以上表格，主要来自http://annovar.openbioinformatics.org/en/latest/ user-guide/filter/)。

英文单词在翻译为中文词语时，由于英文词语和中文词语的含义会有差异，而使得阅读中文译文和阅读英文原文时，带来不同的理解。这种差异多数情况下是无关紧要的，但是在诸如致病程度分类方面，则相当重要。因此以上表格的致病性分类未做翻译。以下为意思不完全一致的中文翻译，读者对致病性的理解应以英文为准。

(4)有害性临界值及有害突变比例的差异

使用预测软件的目的，是给任意一突变以有害与否的判断，但不同的软件所给判断是很不一致的。对于同一批数据，如我们的20个病人的全外数据，用不同的软件，给出这批数据中有害突变的数量和比例，不同软件所给结果很不相同。这种差异不但在于不同软件有不同的分类方法，而且在于，即使是同样的二分法软件，或三分法软件，各个分类所包含的变异的数目和比例也有极大的差异。例如使用SIFT软件预测，有害突变比例为18.52％,而同样二分法的软件MetaSVM，有害突变只有0.99％了。这样，同一突变的有害性在不同的软件中可以得到不同的分类结果。我们假设某一款软件的预测是准确的，以此软件为标准去判断其它软件的准确准确性，那么其它软件都有很多假阳性或假阴性；反之亦然。因此，单纯地相信软件的判断，是不可取的。以下为各软件致病性分类的变异数目统计。

表7各软件致病性分类的变异数目统计

讨论：以上叙述已经明确，针对基因突变致病性的预测软件，在多个方面存在差异，这些差异包括：不同的预测软件存在预测区域和突变类型的差异，存在算法原理的差异，存在致病程度分类、取值范围、方向的差异，存在有害性临界值及有害突变比例的差异。这些差异给突变的致病性分析带来困难。

例如，有害突变的比例在不同软件中的不同，就会造成分析的混乱。在开篇介绍分位数思想时已提到，假设使用不同的测量系统都可以准确测量同一批样本中每个个体的某一特征，则某个体该特征数值在样本中排序位置，不因使用不同的测量系统而不同，但不同的测量系统可以使用不同的阈值将个体划分到不同的概念领域中。以身高为例，假设在美国人的眼里，身高排在全世界成人前10％的个体可以被称为高个子，其余人为非高个子；而在中国人眼里，身高排在全世界成人前5％的个体才有资格称为高个子，其余人等只能归入非高个子，那么，当持有以上观念的一个中国人和一个美国人讨论高个子时，对于处于世界排名7％位置的某个人，他们是可以分别用高个子和非高个子称呼的，但从排名上看可以是同一身高的人。这样，在不同的测量体系下人为地的划分高个子或非高个子，就可能成为交流的障碍。同样地，针对同一批基因突变，在进行有害性预测时，不同软件在给出的预测数值基础上进行人为的划分(每个软件在当初划分时可能都是有充分理由的)，分别归类为有害、无害或其它，而这种划分的界限又有很大的差异，表现为不同软件对同一批突变的预测的有害性比例很不一致，那么，也同样可以阻碍不同软件之间的比较。假设不同软件的预测是准确的，与实际情况都是相符的，那么不同软件的害预值经过排名后，直接比较排名，是更可靠的方法。由于某个个体的排名与样本的总数有关，为了消除这一影响，而使用分位数原理进行处理。

结论：针对同一批突变，不同软件所给预测结果有很大的差异，这为精确的分析带来困难。

因此，使用数学方法处理这些结果，获得可以比较的标准结果，是很有必要的。

3获取害预值-害阶位表格的算法

目的：使用某种算法获取害预值-害阶位表格

原理：分位数思想，即将害预值和害阶位对应起来。在我们研究的11个软件中，两个软件 (LRT，MutationTaster)给出的害预值与最终预测无关，因此我们只能直接使用最终预测。在进行数据处理时，对于LRT，我们分别用3,2,1替代D,U,N；对于Mutation Taster，我们分别用4.5,3.5,2.5,1.5代替A,D,P,N,以方便进行数据处理(这些数字设定不会影响个体的排名，只是数据处理的手段。不影响排名，也就不会影响计算所得到的害阶位)。

材料：已经获取的20个样本的全外测序结果，未去除不同样本的相同突变。

方法：在excel中，PERCENTILE函数可以根据设定的害阶位找到每个软件对应的害预值。使用20个病人的至少含有一个害预值的全外突变数据，我们计算得到了害阶位与害预值的对应表。在计算时，害阶位为0，代表最高的有害性，害阶位为1，代表最低的有害性。由于高有害性数据是我们的研究重点，我们以0.001的步长对前5％的害阶位进行计算，以0.01 的步长对后95％的害阶位进行计算。

结果：

我们得到了各软件害预值与害阶位对应表，我们将这个对应表简称害阶尺，它是用来测量某个软件的害预值以获取害阶位的工具。以20个样品的未去除重复的突变为基础获取的害阶尺，我们可称为普通害阶尺；以去掉重复的突变为基础获得的害阶尺，可称为去重害阶尺。

如下：

各软件害预值与害阶位对应表(普通害阶尺)

各软件的害预值与害阶位关系散点图见图1和图2。

讨论：从表中我们可以看到各软件的一个重要差异，即：对于取值范围同样在[0,1]之间，且越大代表越有害的四个软件Polyphen2_HDIV、Polyphen2_HVAR、MetaLR、Fathmm-MKL，取其中一个害阶位如0.03，对应的害预值分别是1.000、0.994、0.201、0.979，是有很大差异的。因此，不同软件的害预值，即使有同样的取值范围和取值方向，也是不宜拿来直接比较的。从散点图上看，这个关系更为清晰。

结论：

在excel中使用PERCENTILE函数可以得到害预值与害阶位对应表(普通害阶尺)。基于该表绘制的害预值与害阶位散点图，可以直观地反映各软件的害预值与害阶位的关系。

4评估根据不同样本来源数据计算得到的害阶位的稳定性

目的：评估根据不同样本来源数据计算得到的害阶位的稳定性各软件害预值与害阶位对应表是依赖20个病人的全外显子高通量测序的突变数据而计算得到的，具有根据此表就提供的上述任意一款预测软件的害预值计算相应的害阶位的潜在功能。倘若此表本身是不稳定的，在更换一批病人样本之后，会计算得到有较大差别的害预值 -害阶位对应数据，那么这个表的利用价值就降低了。因此我们需要研究使用不同样本的全外高通量测序数据得到的害预值-害阶位对应关系是否稳定，标准误差是否在可容忍的范围内。

材料：已经获取的20个样本的全外测序结果(未去除不同样本的相同的突变)。

方法：我们使用了FATHMM软件的数据分别计算了20个样本的害阶位-害预值对应表，分别于整体计算结果比较，并进行相互比较，以计算差异程度。。

结果：从数值上看，来自各样品的害阶尺具有很好的稳定性，各样本的害预值与根据整体数据计算的害预值是十分接近的，见图3。在每个害阶位，计算20个样品害预值所得到的标准差，可反映各样品在某一害阶位对应害预值的变异情况，见图4。

FATHMM软件的各样本及总体的害预值-害阶位表

(待续)

(续表)

(待续)

(续表)

讨论：对不同样品做害阶位-害预值图，其曲线几乎是完全重叠的。这显示就测试的FATHMM 软件而言，用不同的样本分别计算害预值-害阶位对应表，数值十分接近。

然而每个害阶位对应的不同样品的害预值是不可能完全一致的。对每个害阶位的不同样品的害预值取标准差，会发现标准差在接近0或者接近1的地方较高，而在中间比较低。这反映了不同害阶位的害预值的离散程度。我们知道，标准差与样本数量的平方根成反比(参见标准差的计算公式)，样本数量越多，标准差越小。当我们仅仅使用一个样品就可算出比较稳定的害阶尺时，那么取20个样品求害阶尺，其值是更稳定的，标准差是更小的。又由于预测软件的预测本身是有误差的，害阶尺因使用不同样本而带来的微小误差，与预测软件本身的预测偏差相比，就微不足道了。因此，使用有限样本计算所得害阶尺，是可靠的，是可作为尺子用于测量其它的害预值的。

但应该注意的是，当使用的用于生成害阶尺的原始材料发生改变时，害阶位-害预值对应关系必然也会改变。例如，原始材料中剔除高频突变，或者包含多个样本的原始材料去除重复突变(高频突变的重复率更高)，都会使整体的害预值对应的害阶位变大。类似于来自某群体的成人身高与排名的对应表中，当去除大量普通人群，在留下的个体中，大部分个体的相对身高排名反而降低。因为为了获得稳定的害阶尺，应固定测序Panel和突变入选规则。

结论：从全外测序Panel测到的突变，使用相同的突变筛选规则，以此为材料计算害阶位- 害预值对应关系，则即使只有一个样本，也可以获得比较准确的、可以与以20个样本为材料计算得到的结果相比拟的害阶尺。

5根据害阶尺从害预值计算害阶位的算法

目的：提出两种基于害阶尺从害预值到害阶位的算法

材料：各软件害预值与害阶位对应表(普通害阶尺)，已经获取的20个样本的全外测序结果 (未去除不同样本的相同突变)。

原理：从“各软件害预值与害阶位对应表(普通害阶尺)”构成的曲线可以看出，害预值和害阶位是单调函数，二者基本上具有一一对应的关系。因此，对于给定的任意一软件的害预值，推算出对应的害阶位，在理论上是可行的，是可以通过合适的算法予以解决和优化的。为解决这一问题，所输入的害预值应符合要求：(1)生成害预值的预测软件是表格中已有的预测软件；(2)害预值在每个预测软件的害预值取值范围之内；(3)对于如LRT和Mutation Taster 等只有几个孤立预测结果的软件，输入的数值也必须是这几个数值之一。

对于符合以上条件的给定害预值，主要有两种方法求取害阶位：(1)对于输入的某软件的害预值，在害阶尺的该软件列找到该数值的最佳近似值，再根据此近似值给出对应的害阶位数值。当某软件某个害预值在害阶尺中对应多个害阶位时，取所有对应害阶位中最大值。其流程图见图5。(2)对于输入的害预值，在害阶尺对应软件列中寻找是否有完全相同的数值，如有，则取对应的害阶位，如有多个，则取对应害阶位的最大值，如无，则找到上下相邻的两个害预值和对应的两个害阶位，害预值和害阶位分别充当平面直线坐标系中的X值和Y 值，上下相邻的两个害预值-害阶位构成坐标系中的两个点，两点决定一条直线，进而求得直线方程式。则待求害预值将落在两点之间，根据方程式求取害阶位。第二种的流程图见图 6。

以上两种方法，第一种不太精确，但计算方便，第二种比较精确，但计算复杂。这两种方法的自动计算我们都已经实现。

第一种是近似算法，其思路是：(1)根据20个样本的全外突变数据及相应的各软件害预值，使用Percentile函数生成分软件的害预值-害阶位对应表。(2)将某个需要计算害阶位的某软件的害预值输入害预值-害阶位对应表。(3)根据输入的软件的名称，从对应表找到该软件对应的一列数据。(4)用输入的害预值与该列数据比较，找到最接近的数值。(5)根据找到的数值，找出对应的害阶位。

第二种是精确算法，其思路是：(1)根据20个样本的全外突变数据及相应的各软件害预值，使用Percentile函数生成分软件的害预值-害阶位对应表。(2)将某个需要计算害阶位的某软件的害预值输入害预值-害阶位对应表。(3)根据输入的软件的名称，从对应表找到该软件对应的一列数据。(4)在找到的数据中查找输入害预值出现的次数。(5)如出现次数为1，则指向对应的害预值，再根据害预值提取害阶位。(6)如出现次数多于1次，则指向所有相同的害预值和对应的害阶位，取害阶位最大值作为输入害预值对应的害阶位。(7)如出现次数为0，则判断害预值排序为升序或降序。(8)使用Match函数找到害预值序列中与输入害预值相邻的两个数X1，X2。(9)使用Offset函数，找到X1，X2对应的害阶位Y1和Y2。(10)根据X1，Y1和X2，Y2，在平面直角坐标系中用Tread函数构建通过X1，Y1和X2，Y2两个点的直线方程式。(10)将输入的害预值输入该直线方程式，求取相应的害阶位值。

方法：通过excel公式编程，可分别给出近似算法和精确算法的编程公式，实现给定害预值的害阶位求取。

结果：利用excel编程得到的两种算法的公式，可用于计算害预值对应的害阶位，且结果与通过害预值-害阶位曲线验证的结果一致。

6根据害阶尺从害阶位计算各软件害预值的算法

根据害阶尺从害阶位计算各软件的害预值，是根据表格从害预值计算害阶位的逆运算。在很多情况下，害预值和害阶位是一一对应的，根据某害预值计算得到的害阶位，再用害阶位计算害预值，仍然可以得到与原来害预值相同的值。但在我们的系统中，这种情况并非必然发生的。这是因为在害预值-害阶位表格中，害预值与害阶位并非一一对应的，一个害预值可以对应多个害阶位，因此以某个害预值计算得到的多个害阶位中所选择的那个值，反过来求取害预值，就不一定是原来的害预值了。

由于前面已经给出从害预值计算害阶位的方法，而从害阶位计算害预值与前面思路相似，细节略有差异，且从害阶位计算害预值，其实际使用范围较小，本部分不再展开。

7害阶位与评级值的比较

目的：比较害阶位与评级值，评估二者的差异

材料：20个样本的全外测序结果中的已有评级值，根据害阶尺对所有全外测序结果计算所得到的害阶位。

原理：如前所述，在Xiaoming Liu等人开发的dbNSFP数据库中，为了解决不同预测软件所给的预测值难以直接比较的问题，使用了评级值(rank score)算法。评级值是将dbNSFP 收录的所有突变，分别按不同软件的预测值排序；某个突变的某个软件预测值在该软件所有突变预测值的位置即为评级值。例如某突变的评级值是0.9，代表该突变在该软件的预测算法中比90％的突变更有害。

与评级值相比，害阶位是一个取值范围与评级值相同，都是在0和1之间，但取值方向与评级值正好相反的概念。那么，是否可将害阶位看成评级值数值的简单的数学变形呢？我们将所有突变基于害阶尺计算得到的害阶位与数据库已经给的评级值相比较，通过散点图分析二者的关系。假设害阶位与评级值是同一算法原理的不同实现方式，那么对同一批突变，二者的散点图将构成一条直线或曲线。规则的曲线总可以被模拟出相应的方程式，根据方程式即可进行彼此的换算。如果二者只是相关关系，而不是严格的换算关系，那么二者的散点图将构成真正的散点关系。

方法：我们对我们使用的11个软件中的9个进行了散点图分析，而对害阶位来自有害性分类结果的两个软件，MutationTaster2和LRT，未进行同样分析。

结果：fathmm-MKL_coding软件的害阶位和评级值在平面直角坐标系中所绘制的散点图排列成一条位于第一象限的曲线，该曲线与Y轴相交于(0,1)点，与X轴相交于(1,0)点，散点构成单调平滑曲线，曲面向坐标原点凸出。见图7。根据害阶位和评级值的定义，二者都在(0,1)的区间之内，但大小是相反的，即害阶位越大，评级值越小，反之亦然。散点图的结果与之相符。这提示，对该软件而言，害阶位和评级值可能存在简单的换算关系。

但PROVEAN软件的害阶位与评级值所构成的散点图与fathmm-MKL_coding软件的散点图有较大差异。PROVEAN软件的散点图除了具有和fathmm-MKL_coding软件的散点图类似的曲线之外，另外还有大量的散点分布于曲线之上，见图8。这提示，对PROVEAN软件而言，有一部分突变的害阶位和评级值存在函数换算关系，另一部分突变的害阶位和评级值则不符合这样的换算关系。

讨论：从图中可以看出，fathmm-MKL_coding软件的突变的害阶位和评级值的散点构成平滑的单调曲线，而PROVEAN软件的害阶位与评级值的散点图除了有一条光滑曲线外，还有大量散点分布于曲线的上方。与PROVEAN软件的结果类似的还有FATHMM，SIFT，Polyphen2_HDIV，Polyphen2_HVAR四款软件；与fathmm-MKL_coding结果类似的还有MutationAssessor，metaSVM，MetaLR三款软件。如果说fathmm-MKL_coding软件的害阶位是评级值的另外的数学形式，那么PROVEAN等软件中，害阶位与评级值是虽有关但同时有较大差异的概念。对于PROVEAN等软件，假设dbNSFP数据库的确如其论文所宣称的那样，评级值是害预值的位置函数，是通过分位值思想从害预值计算而来，那么，我们使用类似方法从同一批害预值数据计算得到的害阶位，在散点图中应与评级值构成平滑曲线。然而大量的分布于平滑曲线之外的散点，提示评级值的计算可能并没有完全依赖dbNSFP数据库在论文中所公开的算法。这种实际使用的算法究竟是怎样的，我们不予猜测。

结论：害阶位在有些软件中，可视为评级值的另外的数学形式，而在另外的软件中，则是有相当差异的。因此从整体上看，害阶位不是评级值的简单的数学转换，而是有着自己独特算法的概念。

8集中度的计算

目的：提出一种计算DM-Pathogenic突变在害阶位编秩上分布集中度的算法

原理：使用害阶位算法处理突变在诸软件中的害预值的方法已经建立了。一批突变经由诸软件和害阶位算法给出了害预值和可以用统一的标准相互比较的害阶位后，如何利用害阶位检验各软件的预测准确性呢？其一般方法是，统计各个软件的每个害阶位区间(比如从0到 0.01区间，从0.5到0.51区间)包含的DM-Pathogenic突变的数量，以评估软件的预测能力。 DM-Pathogenic突变是指在HGMD数据库标记为DM或者在ClinVar数据库标记为Pathogenic 的突变，是业内大部分人都认可的致病突变。但两个数据库并不一致，具体见后面叙述。倘若DM-Pathogenic突变在某个软件的害阶位上主要分布在0值附近的区间，就说明该软件具有很好的预测能力；如果DM-Pathogenic突变是随机分布的，并没有集中分布在0值附近，说明该软件的预测在很大程度上与人为划分的突变的有害性不符，以DM-Pathogenic突变为标准，这款软件的预测能力是存疑的。但如何计算DM-Pathogenic突变在害阶位为0附近的集中度呢？

在为这个问题提供可行的算法之前，我们可将该问题抽象化为如下问题：在从A到B的方向划分连续的空格，每个空格分布不同数量的物品，如何从A点出发，走较少的空格获得较多的物品，就是物品在A向的集中度问题。所谓计算DM-Pathogenic突变在害阶位0值方向分布的集中度，就是计算A向的集中度。

方法：我们设计了一种经验算法以计算上述A向集中度问题。其布局如下表，位置参数是从大到小依次排列的整数，排列到最后是1，反映了每个空格的位置；各数据组是待计算集中度的数据。目的是计算各组数据在A方向的集中度。计算思路是，设数据组1共有n个空格，每个空格有一个正数或0，则这些数字自上而下依次是X1，X2，…Xn,这些数字的和是Sum。如求X1在数组中的比重，则简单地用X1/Sum即可。如求第一和第二个空格在数组中的比重，则用公式(X1+X2)/Sum.如果求A方向这两个空格数字的比重，也即越靠近A 方向越重要，则用公式(X1+X1+X2)/(Sum+Sum)。同样的，如计算整个空格在A方向的集中度，可以用公式[X1×n+X2×(n-1)+X3×(n-2)+…+Xn×1)]/(Sum×n)，这个公式是前一个公式的简单变形。此时N变为位置参数，如下表左侧所示。使用该经验公式，在位置参数固定的情况下，根据最后集中度的值的大小，可以比较几个数据组的A向集中度。如果设定为B 方向的数据越大，A方向的数据集中度越小，则可以在位置参数上使用负数，如下表右侧。这样更改后，改变的是集中度的数值，并不影响各个数字之间的相对大小。

结果：如下表所示，数据组2的非0数字整体向下移动一位，空出来的空格以0填充，相应的集中度也会变小。而数据组3虽然与前两个数据组包含相同的数字，但集中在B方向，因此计算得到的A向集中度最小。假设数据组1的与位置参数14对应的空格是某个正数，其余全部为0，则集中度的值为1。更改正整数的大小，集中度不变，仍为1。因此，集中度是与所计算对象的数值总和无关，只与数值在空格中的分布的相对大小有关的参数。在确定了空格数量，位置参数后，集中度数值大小可用来比较空格中数字分别的集中程度。

在excel表格中，假设以上表格的“14”位于B3单元格，则计算数据组1的集中度的公式是：SUMPRODUCT($B3:$B16,C3:C16)/(SUM(C3:C16)*14)；计算数据组2的集中度的公式是： SUMPRODUCT($B3:$B16,D3:D16)/(SUM(D3:D16)*14)；计算数据组3的集中度的公式是：SUMPRODUCT($B3:$B16,E3:E16)/(SUM(E3:E16)*14)

结论：以上这些结果，与人的主观判断基本一致，因此集中度算法作为经验算法，可以定量比较有害突变在软件害阶位排序中分布的集中程度，进而评估预测软件的预测准确性。

9害阶位与HGMD，ClinVar数据库注释的比较分析

目的：用HGMD，ClinVar数据库判断各软件对突变致病性预测的准确程度

材料：20个样本全外测序经去重处理的全部突变中，(1)具有HGMD数据库或ClinVar数据库注释，(2)具有软件预测注释；凡是符合这两个条件的所有突变，皆予以入选。

原理：将突变的预测软件害预值转化为各软件的害阶位数值，是否有利于判断突变的有害性呢？依害阶位数值判断的突变的有害性，是否可反映突变的实际有害性呢？解决这一问题的一般思路是将害阶位与突变数据库如HGMD，ClinVar进行比较，计算二者的一致性。倘若这些数据库注释为有害或可致病的突变，其害阶位数值较小，则提示软件预测与数据库记录具有较好的一致性；相反的，如数据库记录为有害，而软件害阶位较大，则二者一致性较差。这种方法是以数据库的记录为标准来判断软件预测结果的可靠性。如果数据库记录本身不准确，那么对软件预测结果的判断就有偏差。事实上，目前一般认为，HGMD数据库和ClinVar 数据库的判断大部分是准确的，少部分是不准确的。特别是HGMD数据库，由于起步较早，所收录的突变和突变致病性判断来自前期文献，此时用于判断突变有害性的国际通用规则如 ACMG尚未诞生，因此对突变致病性判断的标准较宽松。而ClinVar数据库最开始则是由几家诊断公司的内部数据集合而成，对突变的致病性判断总体上遵循了ACMG规则，故致病性判断的标准更严格一些。由于判断标准的不同，相当一部分在HGMD数据库中标记为“DM”的突变，在ClinVar数据库中标记为“Benign”。例如，我们收集的20个样本中的全部突变中，去除重复突变，HGMD数据库或ClinVar数据库有记载的突变共有3963个，其中DM-Pathogenic突变共有204个。这些突变整理如下：

204个突变的数据库标记

(HGMD数据库对致病性标记的解释：DM：文献支持的致病突变；DM？：文献支持的疑似致病突变；FP：体外实验支持的功能性SNP；DP：疾病相关SNP；DFP：体外实验支持的与疾病相关的功能性SNP。ACMG及ClinVar数据库对致病性标记的解释：Benign：良性； Likelybenign：疑似良性；Uncertain significance：不确定；Likely pathogenic：疑似致病；Pathogenic：致病；not provided：未提供；other：其它。)

在204个突变中，HGMD数据库记录为DM而在ClinVar数据库中归类为Benign的突变有12个，占5.88％。在HGMD数据库中归类为DM而在ClinVar数据库中归类为Pathogenic 的致病类型，只有15.69％。这些差异使人有理由怀疑这些归类的准确性。但在没有其它更好的参考数据的情况下，只能以这些数据库的记录为标准检验预测软件对突变有害性预测的准确性。

方法：我们根据“各软件害预值与害阶位对应表(普通害阶尺)”，使用从害预值计算害阶位的第二种算法对20个样本不重复突变的害预值进行计算，得到对应的害阶位，然后按照百分位的方式，统计每个害阶位区间所包含的DM-Pathogenic突变数量，得到了每个害阶位中DM-Pathogenic突变数量的表格。再使用集中度计算公式计算DM-Pathogenic突变在有害方向的集中度，以比较各软件的突变有害性预测与人为划分的有害性判断之间的相对一致性。

结果：根据“各软件害预值与害阶位对应表(普通害阶尺)”计算的20个样品的DM-Pathogenic 突变的区间分布表

(注：以上右侧所谓“害阶位平均值”，是指将每个突变在各个软件计算得到的害阶位进行平均，通过平均值求取对应害阶位区间的DM-Pathogenic突变的数量)

讨论：通过以上计算，比较各个软件所得到的DM-Pathogenic突变分布和集中度，我们发现， MetaLR软件可得到最高的集中度。而LRT软件得到的集中度最低。这是因为LRT的数据分成D，N，U三种情况，每种情况的害阶位计算取最低害阶位为共同结果，从而使LRT软件以及类似的Mutation Taster软件出现集中度数值偏低的现象。

结论：通过以上计算，我们发现：(1)人为判断的突变有害性与软件判断的害阶位是相关的，有害的突变总是分布在害阶位0值附近。这显示以害阶位值反映突变有害性的程度是合适的。(2)以集中度为标准，不同软件对突变致病性的判断是有差异的。(3)使用各软件害阶位的平均值获取有害突变在害阶位上的分布，以集中度为标准，相比于其他软件并无优势。

10害阶位倒数与HGMD等数据库注释的比较分析

目的：使用害阶位倒数计算DM-pathognic突变的集中度

材料：各软件害预值与害阶位对应表(普通害阶尺)，已经获取的20个样本的全外测序结果 (对突变进行了去重处理)。

原理：

我们发现，一个DM-pathognic突变经诸软件预测并计算获得诸软件的害阶位后，不同软件的害阶位可有很大差异，这反映了不同预测软件对同一突变致病性判断的不一致性。如果某突变X有一两个害阶位值较大(反映这一两个软件判断该突变无害)，那么即使其它害阶位的值极小(反映其它软件判断该突变有害)，也会得到较大的平均值，故以害阶位平均值对突变排序，突变X会因为平均害阶位较大而排列在有害性较低的位置。我们认为，这种大值在平均值计算时对小值的掩盖作用，会降低突变有害性判断的准确性。为检验这一假设，我们尝试使用害阶位的倒数求平均值，这样某突变某软件的害阶位值越小，计算得到的倒数越大，进而得到较大的倒数平均数。而害阶位最大值为1,1的倒数为1，在动辄为十几或几十的害阶位倒数值中，为1的数字对平均值的影响是小的，这样害阶位小值对害阶位倒数平均值的影响就增加了。在计算时，由于某些害阶位的值是0，为了避免对0取倒数，我们将所有害阶位加上0.01再取倒数。经过这样的处理，大数(大的害阶位)变成小数(害阶位的倒数值较小)，小数变成大数，就某个软件的全部突变的数值而言，数列的方向因为从害阶位变成害阶位的倒数，按大小的排序发生了颠倒，但是每一个数据都不会因此而错位。而突变诸软件的害阶位的平均值，在更改为害阶位倒数的平均值后，不同的突变在排序中的位置将发生改变。例如假设突变A和B的两个软件的害阶位分别是A：0.1,0.1；B:0.01,0.19,则二者的平均值皆为0.1，取倒数后的平均值则分别为10和52.6，二者出现较大的差异，自然分别以诸软件害阶位平均值排序或者以诸软件害阶位倒数平均值排序，同一个突变在将有不同的排列位置，以百分位换算该位置，也会有不同的百分位数值。

害阶尺是基于20个病人的全外测序结果(含有178081个带害预值的未去重的突变)计算得到的害预值-害阶位对应表。从百分位定义可知害阶位的含义是，害阶位值在某个害阶位区间(如0.50-0.51)内的所有突变数目占总突变数目的比例，与区间的大小(0.50-0.51的区间大小为0.01)是一致的。使用诸软件害阶位倒数的平均值排序，突变的排序位置打乱，以害阶位各个百分位阈值的倒数为阈值，则害阶位倒数平均值落在两个阈值构成的区间内的突变数目，与这两个阈值换算前的理论数目，可能会有很大的偏差。因此，要计算某个突变的诸软件害阶位倒数的平均值在所有突变中的位置，需要重新计算害阶尺。

因此，我们在去掉重复突变的基础上重新计算新的害阶尺，可称为去重害阶尺，在此基础上分析我们的假设。

方法：(1)基于20个样本的全外测序结果中未去重的178081个带害预值的突变，计算害预值-害阶位对应表(普通害阶尺)。该步骤已经完成。(2)对20个样本的全外测序结果的178081 个带害预值的突变进行去重处理，不同样本的相同突变只保留一个，获得27615个突变。(3) 对去重的每个突变的诸软件害预值，重新计算害预值-害阶位对应表(去重害阶尺)。(4)根据去重害阶尺，计算每个突变诸软件害预值的害阶位。(5)根据计算的害阶位，换算为修正倒数。(6)对同一突变诸软件的修正倒数求平均值。(7)根据修正倒数和修正倒数平均值，求修正倒数-害阶位表，修正倒数平均值-害阶位表。(8)根据第七步的表格，求取各个害阶位区间含有的DM-Pathogenic突变的数目。(9)根据集中度计算方法，求各软件及修正倒数平均值的集中度。

结果：以下两个表格，是以去掉重复突变的害预值为材料，分别使用害阶位计算去重害阶尺和使用害阶位的倒数计算去重害阶尺，之后计算每个软件的DM-Pathogenic突变的集中度的结果。

使用害阶位统计每个害阶位区间的DM-Pathogenic突变的数目

使用害阶位倒数统计每个害阶位区间的DM-Pathogenic突变的数目

结果总结如下：(1)根据普通害阶尺和去重害阶尺计算得到的某软件的DM-Pathogenic突变的集中度是不同的。普通害阶尺计算得到的集中度偏大而去重害阶尺计算的集中度偏小。(2) 使用害阶位统计每个害阶位区间的DM-Pathogenic突变的数目，与使用害阶位倒数统计每个害阶位区间的DM-Pathogenic突变的数目，其集中度在各软件中都是完全一样的。(3)以各软件害阶位倒数的平均值求集中度，比以各软件害阶位平均值求集中度，其数值要高。

讨论：我们发现：(1)以20个标本的去掉重复的突变的害预值为材料重新计算害阶尺，称为去重害阶尺，以去重害阶尺重新计算每个突变的害预值的害阶位，进而计算得到不同软件的DM-Pathogenic突变的集中度，我们发现，相对于使用普通害阶尺计算的集中度，所有集中度数值都变小了，这是很容易理解的。在基因突变中，频率越高的突变，有害的可能性越小，这些突变在未去重的突变中，占据了大量的排序空间。去掉重复突变后，大量原先高频无害突变占据的排序空间大量释放，原先位于有害性上端的突变下移，那么分布于这些突变中的DM-Pathogenic突变自然其集中度也会降低。这是一种依赖于参照系不同而导致的变化，不是DM-Pathogenic突变分布本身的变化。(2)基于害阶位得到的集中度，与基于害阶位的倒数所得到的集中度，对于所有软件都是完全一样的。这表明害预值加上0.01后取倒数(修正倒数)，并不会打乱排序的相互位置，不会影响每个区间DM-Pathogenic突变的分布，自然也就不会影响集中度的值。(3)使用各软件害阶位倒数的平均值求集中度，集中度增加，由原来的0.756增加为0.789。我们是这样理解这一改变的：一个突变的诸软件害阶位，如果某个软件的害阶位很小(对应的倒数很大)，意味着这款软件有较高的把握认定该突变有害。此时即使其它软件都认为突变无害，仍然有很高的概率是有害的。而某个软件对突变有害性的较有把握的判断，是与人为判断的结果有更好的一致性的。但集中度增加的数值较小，也显示了这种效应并不是特别强。

结论：使用突变的害阶位的修正倒数的平均值，较之使用突变的害阶位的平均值，对于 DM-Pathogenic突变而言在排序后可以获得更高的集中度。

11“损低固高法”集中有害突变

目的：从以上分析中我们得到一个可能的规律：当几个软件同时对某突变进行预测时，如果一款或少数几款软件极有把握地判断该突变为有害突变，给出较小的害阶位，那么这个突变很可能是真有害的。这个判断的另一个方面是，如果多款软件都认为某个突变有害但都没有把握，给出较高的害阶位值，那么这个突变很可能没有害。使用多款软件的害阶位倒数求平均值，是一定程度上使个别软件的有把握的有害性判断不至于湮灭在平均数之中的一个方法。我们认为这可能是一个规律，为了验证这一假设，我们对数据做了进一步的处理。

原理：如前所述，我们发现，以20个样本的去重的27615个突变为材料，使用害阶位修正倒数的平均值，较之使用害阶位的平均值，DM-Pathogenic突变可以获得更高集中度。这是因为我们使用数学方法，放大了某突变诸软件预测值中，对高有害性比较有把握的预测值。这种数学方法是一种损低固高法，即越大的数值(害阶位的倒数)在经数学处理后保留的比例越大；越低的数值，同样处理后数值损失的比例越大。

方法：将前述的倒数平方后除以100(简称平除)。这样，数字越小，数字损失的比例越大；数字越大，在进行“损低固高法”处理后保留下来的数值越大，就越能在进行平均数处理中排列在较高的位置。对倒数进行平除处理后，对其结果再次进行平除处理，如此反复四次，每次处理都重新计算新值-害阶位表格，根据表格计算各软件的DM-Pathogenic突变的集中度。

结果：经过以上处理，各软件DM-Pathogenic突变的集中度如下表：

各软件DM-Pathogenic突变的集中度计算

在上表中，以害阶位加0.01后的倒数(修正倒数)为材料计算的DM-Pathogenic突变的集中度，在各软件中与以害阶位为材料计算的DM-Pathogenic突变的集中度并无差别，但各倒数的平均值达到0.789，比0.756有提升明显，且该数值高于所有软件的集中度。对害阶位的倒数进行一次、二次、三次、四次平除后，以各数值的平均值为材料计算害阶尺，在各自的害阶尺下计算所得DM-Pathogenic突变的集中度，是先上升而后下降的，而不是持续上升的。因此我们平除三次之后，不再计算。

对于以上数据，当我们仅仅关注Pathogenic突变而不是DM-Pathogenic突变的集中度时，其结果如下：

各软件Pathogenic突变的集中度计算

DM-Pathogenic突变和Pathogenic突变的集中度见图9。

讨论：可见无论是计算DM-Pathogenic突变的集中度，还是Pathogenic突变的集中度，都有共同的特点，即在取害阶位的倒数求平均值后，以平均值计算得到的集中度都明显增加，之后进行“损低固高法”进行处理，集中度不太明显地增加或降低。因此，使用“损低固高法”计算集中度，有一定价值，但不是无限的。

结论：在实际运算中，使用害阶位的修正倒数的平均值排序，足以获得比较好的排序形态，使有害突变集中地排列在序列的一端。

12害阶位概念在二代测序分析中的应用

目的：将害阶位原理和方法应用在二代测序的分析过程中

原理：普通人或病人经过全外显子测序，将所得数据与标准的人基因组参考版本(HG19， HG38等)比对，可以得到几万个变异位点。所谓遗传病高通量测序的分析，就是从几万个变异位点中，找到导致病人表型的变异位点。在基因型、表型、遗传型三个方面都符合条件的某个基因的突变位点，可判断为致病性突变位点。突变符合基因型的条件，是指突变的有害性得到证明，即突变后该基因编码的蛋白质结构发生改变，不能正常行使其应有的功能。这一点通常是通过对突变类型的分类及文献查阅进行判断的。所谓突变符合表型的条件，是指突变所在基因编码的蛋白质发生功能异常后，经过一系列传递，在人的相貌，代谢，智力、运动、语言等等一个或多个方面表现出的有规律的改变，这种有规律的改变记载于文献或数据库中，与待测病例的表型一致。突变符合遗传型的条件，是指突变符合遗传规律。例如父母表现正常而孩子有遗传病，那么对于隐性遗传的基因，孩子为纯合突变或复合杂合突变，父母为杂合突变或未携带突变就符合遗传型规律。再如父亲和孩子有同样的症状，母亲无症状，则对于显性遗传的基因，父亲，孩子携带而母亲不携带的突变，就符合遗传型规律。

所谓遗传病的高通量测序和分析，需要首先对海量的变异位点进行注释，注释包括：这个位点的位置是外显子还是内含子？是什么基因？突变的类型和名称是什么？是SNP吗，如果是，SNP的名称是什么？这个突变在人群数据库的频率是多少？这个突变在多个预测软件中的预测怎么样，预测是否致病？测序质量和测序深度怎么样？这个基因与什么疾病有关？是否有家系分析，家系分析怎么样？有些公司在这些注释的基础上，会给出更多的注释数据。再基于这些注释对变异进行分析时，其一般的分析流程是：(1)去除高频突变。凡是在人群中高频出现的突变，肯定是无害的。但对于不同遗传类型的基因，针对检测对象的种族，以及某些特殊基因，应对频率设置不同的cut-off值。去除了高频突变后，通常可去除90％的变异位点。(2)判断变异的类型、在各变异有害性预测软件的预测结果、在HGMD、ClinVar 等数据库的记录情况，并分别打分。变异的类型中，有些是极大可能致病的，如移码突变，与外显子临界的内含子上一两位的突变，则分值应该高一些。而同义突变、离外显子边界较远的内含子突变等，显然分值应该低一些。同样的，在HGMD数据库中记录为DM的突变，在ClinVar数据库记录为Pathogenic的突变，分值应高一些，而记录为良性突变的，分值应该低一些甚至负分。但是还有大量的非同义突变(占低频突变的80％左右)，人群中分布频率很低，ClinVar数据库等无记录，不能明确判断其有害性的，而只能通过突变有害性预测软件的预测，给出初步结果。(3)当使用突变类型、突变有害性预测软件预测、数据库记录情况进行判断并给出综合打分后，就可以排序，删除分值较低的变异，保留分值较高的变异。 (4)判断突变的遗传型是否符合范式。(5)将同一基因的不同突变，无论其分值高低，排列在一起。(6)使用表型数据，从高分值向低分值的变异逐步比对，研究这些突变基因育表型的契合度，初步给出可能导致表型的基因。如不存在这些突变，可扩大比对范围。如扩大范围仍找不到者，可初步判断为阴性结果。

在这个过程中，突变有害性预测软件的预测，虽然不能构成致病性的根本证据，但用做提示作用是足够的。因此将预测软件的预测定量化，而不是每个结果都通过人眼的判断，是实现二代测序数据分析的自动化或半自动化的的重要步骤。

方法：我们在将害阶位概念运用到二代测序分析时，数据分析步骤是：(1)根据害阶位的概念，将软件的害预值换算成为害阶位的修正倒数，同时保留各软件的判断；(2)将同一突变的各个软件的害阶位的修正倒数求平均值；(3)用平均值排序；(4)使用条件格式，用双色刻度标记不同的数值，用不同的颜色标记不同的判断结果(N,T,A等)。我们一般用红色表示高致病性预测，用绿色表示良性预测。在黑白图片下，则以数据条反映有害性的强弱。

结果：在进行以上处理后，其结果见图10。

软件对突变的有害性预测较低时，分值也低，excel用绿色标示；软件对突变的有害性预测较高时，分值也高，excel用红色标示。红色较多的突变，数据库认定属于有害的”DM”及”Pathogenic”的可能性也大，而绿色区域较多的突变，数据库认定属于无害的”DP”，“DFP”及”Benign”的可能性也大。用红色和绿色对应计算的数值，可以一目了然地反映突变在预测软件中致病的可能性。在本专利申请时，红色以较多的黑色数据条表示，绿色以无数据条的白色代替。

结论：在二代测序中，有大量的错义突变无数据库记录，快速地初步判断其致病性的方法是依赖预测软件的预测。害阶位计算提供了标准化软件预测结果的方法，为二代测序分析提供了有效工具。

Claims (8)

1.本发明公开了一种基于基因突变的多种软件有害性预测值（害预值）计算某突变的诸软件害阶位，以实现标准化处理基因突变的诸软件有害性预测值的算法；该算法的特征在于，由三个分算法构成：分算法1，根据分位值原理从大量基因突变的多种软件害预值注释中构建害预值－害阶位对应表（该表称为害阶尺）；分算法２，以害阶尺为基础，以害阶尺中含有的任意一款软件的任何真实的害预值为输入数据，计算该害预值的害阶位；分算法３，将前述分算法做适当整合和改造，联合集中度算法，以应用于基因突变的有害性预测之中；凡是用到上述分算法中的一种或几种的算法，视为被本权利要求书的权利要求所覆盖。

2.根据权利要求1所描述的害阶位，其特征是：基于某个设计的Ｐａｎｅｌ或全外显子、医学外显子、全基因组等所设定的基因测序的区域范围，某个样品或某批样品或几批样品的全部或经过某种规则过滤而来的突变中，为了判断突变的有害性，而使用一种或多种突变有害性预测软件对突变的有害性进行预测，在得到的有害性预测值的基础上，对于某突变的某软件有害性预测值在该软件本次所分析的突变的有害性预测值中的排列位置，即为该突变在本次分析的突变群体中在该软件中的害阶位，该值反映了比该突变更有害或更无害的突变的比例。

3.根据权利要求1，所描述的用于预测基因突变有害性预测值的软件，其特征是，包括但不限于以下软件：ＳＩＦＴ＿ｓｃｏｒｅ，Ｐｏｌｙｐｈｅｎ２＿ＨＤＩＶ，Ｐｏｌｙｐｈｅｎ２＿ＨＶＡＲ＿ｓｃｏｒｅ，Ｐｏｌｙｐｈｅｎ２＿ＨＶＡＲ＿ｐｒｅｄ，ＬＲＴ，ＭｕｔａｔｉｏｎＴａｓｔｅｒ，ＦＡＴＨＭＭ，ＰＲＯＶＥＡＮ，ＶＥＳＴ３，ＭｅｔａＳＶＭ，ＭｅｔａＬＲ，Ｍ－ＣＡＰ，ＣＡＤＤ＿ｒａｗ，ＤＡＮＮ，ｆａｔｈｍｍ－ＭＫＬ，ｆａｔｈｍｍ－ＭＫＬ＿ｃｏｄｉｎｇ，Ｅｉｇｅｎ＿ｃｏｄｉｎｇ＿ｏｒ＿ｎｏｎｃｏｄｉｎｇ，Ｅｉｇｅｎ－ｒａｗ，Ｅｉｇｅｎ－ＰＣ－ｒａｗ，ＧｅｎｏＣａｎｙｏｎ，ｉｎｔｅｇｒａｔｅｄ＿ｆｉｔＣｏｎｓ，ｉｎｔｅｇｒａｔｅｄ＿ｃｏｎｆｉｄｅｎｃｅ＿ｖａｌｕｅ，ＧＥＲＰ＋＋＿ＲＳ，ｐｈｙｌｏＰ1００ｗａｙ＿ｖｅｒｔｅｂｒａｔｅ，ｐｈｙｌｏＰ２０ｗａｙ＿ｍａｍｍａｌｉａｎ，ｐｈａｓｔＣｏｎｓ1００ｗａｙ＿ｖｅｒｔｅｂｒａｔｅ，ｐｈａｓｔＣｏｎｓ２０ｗａｙ＿ｍａｍｍａｌｉａｎ，ＳｉＰｈｙ＿２９ｗａｙ＿ｌｏｇＯｄｄｓ，Ｉｎｔｅｒｐｒｏ＿ｄｏｍａｉｎ，ＧＴＥｘ＿Ｖ６＿ｔｉｓｓｕｅ，ＲＥＶＥＬ。

4.根据权利要求1，分算法1的特征在于，为了构建害预值－害阶位对应表，在ｅｘｃｅｌ中可使用计算百分位的ＰＥＲＣＥＮＴＩＬＥ函数；一切软件、计算机语言或算法，针对基因突变诸软件的害预值计算害阶位，凡是可以达到与ｅｘｃｅｌ中ＰＥＲＣＥＮＴＩＬＥ函数类似效果的程序、软件、函数，即视为用到了分算法1；分算法1的特征还在于，害预值－害阶位对应表可以一次构建，多次使用，也可以在每次使用时，以所分析的对象为基础临时构建再行使用；分算法1的特征还在于，害预值－害阶位对应表可以以表格的形式存在，也可以以曲线的方式，或计算机字典、计算机数组的方式存在。

5.根据权利要求1，在分算法２中，以害阶尺为基础，根据输入的害预值计算害阶位的算法主要包括两种方式：一种是在害阶尺的该软件列中找到害预值的最接近的值，根据该值找出对应的害阶位，一种是在害阶尺的该软件列中找到与输入害预值相邻的两个害预值，通过建立局部曲线或直线的方式，计算出输入害预值的理论害阶位。

6.根据权利要求1，在分算法３中所谓集中度算法，是对人为判断的有害突变，在害阶位编秩上分布集中度的算法，该算法以人为判断的有害突变为标准，考察某软件或软件组合对突变有害性预测的准确性；所谓适当整合和改造，包括但不限于：通过对数据格式（色阶、数据条、图标集等）的改造以反映数据的大小；通过使用害阶位的倒数（或其它基于害阶位的数学变形）来联合多软件的分析结果；将所得数据做适当的数学变形，与其它的结果（突变类型、数据库数据、在人群中的分布频率等）一起计算和展示、分析等。

7.根据权利要求1，算法的特征是，可用任何计算机语言编写实现该算法，包括但不限于ＪＡＶＡ、ＶＢ、ｐｙｔｈｏｎ、ｑｂ、ｃ＋＋、ｖｃ＋＋、ｃ语言等；算法可运行的计算机操作系统包括但不限于：ｗｉｎｄｏｗｓ系列，ｄｏｓ，ｍａｃ ｏｓ系列，ｌｉｎｕｘ，ｕｎｉｘ、Ａｎｄｒｏｉｄ、ｉＯＳ等；算法的存在方式包括但不限于：单机版、网络在线版、内置于基因分析仪上、以模块方式或其它方式存在于其它软件之中。

8.根据权利要求1，所描述的基因突变的特征在于，在基因组水平，在ＲＮＡ水平，在ｃＤＮＡ水平，突变类型包括但不限于：碱基置换（Ｓｕｂｓｔｉｔｕｔｉｏｎｓ）；缺失（Ｄｅｌｅｔｉｏｎ）；重复（Ｄｕｐｌｉｃａｔｉｏｎ）；插入（Ｉｎｓｅｒｔｉｏｎ）；倒位（Ｉｎｖｅｒｓｉｏｎ）；缺失／插入（Ｉｎｄｅｌｓ）；突变位置包括但不限于：编码蛋白质的外显子区域，不编码蛋白质的外显子区域，内含子区域，５’－ＵＴＲ区，３’－ＵＴＲ区，以及其它尚未归类的ＤＮＡ序列；其变异来源包括但不限于：体细胞突变，癌症相关基因突变，遗传，以及其它未知因素造成的变异；用于发现变异的参考序列包括ＨＧ1８，ＨＧ1９，ＨＧ３８，ＨＧ３９等各个版本的人基因组序列，包括来源于ＮＣＢＩ数据库、ｅｎｓｅｍｂｌ数据库，ｕｃｓｃ数据库等各种数据库的基因序列，或其它来源的序列。