CN109949868B - 基于耐受性分析的基因等级排序方法和装置 - Google Patents

Abstract

一种基于耐受性分析的基因等级排序方法和装置，该方法包括：获取频率数据库和转录本数据库的数据；从转录本数据库获取各个基因的转录本作为编码序列数据；从频率数据库中提取目标基因序列；将编码序列数据映射到目标基因序列，检测得到目标基因序列上的功能性突变和非功能性突变；选取功能性突变和非功能性突变中等位基因频率在第三预设值以下的突变作为基因等级排序的分析对象；对分析对象建立基因等级排序表，得到多个不同基因的基因等级分数，以此表征该基因的耐受性。本发明以基因为单位，利用功能基因对疾病的耐受性，建立基因等级排序，在功能基因的预测上具有较高的可信度和准确性。

Description

基于耐受性分析的基因等级排序方法和装置

技术领域

本发明涉及生物技术领域，尤其涉及孟德尔遗传病的候选基因筛选方法，特别是基于耐受性分析的基因等级排序方法和装置。

背景技术

随着生物医疗水平的不断提升，DNA测序技术的地位不断提升，二代测序技术应运而生，包括全基因组测序(WGS)、全外显子测序(WES)及目标区域测序等。基因组测序是对未知基因组序列的物种进行个体的基因组测序，可发现基于DNA水平的点突变、微小片段的插入、缺失、拷贝数的变异及基因组结构变异，进而探寻基因突变所致的功能改变，与之对应孟德尔疾病或单基因疾病的致病基因发病机理。

WGS和WES检测技术已经非常成熟，变异检测流程在各个基因组检测中心都做得比较好，大多可以自动化形成分析流程。其困难和瓶颈在于后半部分医疗信息分析，需要从检测到的变异中筛选出可能的致病突变，这个过程非常依赖医学相关背景。现在缺乏的是能够将生物信息分析和医疗信息分析进行协调，并有能力从各个数据库遴选出有用的信息，在不同的分析阶段，进行不同程度的侧重注释。

功能基因的耐受性，是指对于一些有功能性突变的基因而言，在健康个体中含有较少的功能性突变的基因比含有较多的功能性突变的基因更能导致某种疾病的发生。例如，基因A有5个功能性突变，基因B有10个功能性突变，那么基因A的耐受性比较低，基因B的耐受性比较高，基因A的变异更容易导致疾病的发生，所以它的危害性高。功能性突变多数为有害突变，已经忍受了10个有害突变的基因B，耐受性更高，危害性更低。

目前，还没有基于耐受性分析对基因进行等级排序的技术用于涉及孟德尔遗传病的候选基因筛选中。

发明内容

本发明提供一种基于耐受性分析的基因等级排序方法和装置，以基因为单位，利用功能基因对疾病的耐受性，建立基因等级排序，在功能基因的预测上具有较高的可信度和准确性。

根据第一方面，一种实施例中提供一种基于耐受性分析的基因等级排序方法，包括：

获取频率数据库和转录本数据库的数据，其中上述频率数据库包括多个样本的测序数据，上述转录本数据库包括多个转录本的测序数据；

从上述转录本数据库获取各个基因的转录本作为编码序列数据，上述编码序列数据包括每个外显子的两端碱基；

从上述频率数据库中提取目标基因序列，上述目标基因序列是测序深度达到第一预设值的碱基位点的覆盖度达到第二预设值的基因序列；

将上述编码序列数据映射到上述目标基因序列得到上述两端碱基之间的序列，检测得到上述目标基因序列的上述两端碱基之间序列上的功能性突变和非功能性突变；

选取上述功能性突变和非功能性突变中等位基因频率在第三预设值以下的突变作为基因等级排序的分析对象；

对上述分析对象建立基因等级排序表，得到多个不同基因的基因等级分数，上述基因等级分数表征上述基因的耐受性。

在优选实施例中，上述频率数据库是ExAC数据库；上述转录本数据库CCDS数据库。

在优选实施例中，上述从上述转录本数据库获取转录本作为编码序列数据，包括：对每个外显子的两端进行两个碱基的延长得到上述每个外显子的两端碱基，上述两个碱基的延长用于提供剪接受体和供体位点。

在优选实施例中，上述基因包括多个上述转录本，上述从上述转录本数据库获取转录本作为编码序列数据，还包括：将该基因的所有转录本融合成为一个转录本作为上述编码序列数据。

在优选实施例中，上述第一预设值是10倍深度，上述第二预设值是70％。

在优选实施例中，上述功能性突变包括错义突变、无义突变和剪切突变；上述非功能性突变包括同义突变。

在优选实施例中，上述第三预设值是1％以下的数值，优选0.01％以下的数值。

在优选实施例中，上述对上述分析对象建立基因等级排序表，包括：

以上述功能性突变的总数作为x轴，突变的总数作为y轴作图，其中上述突变的总数包括功能性突变和非功能性突变的总数；

上述x对上述y回归，将得到的标准残差作为基因等级分数来表征上述基因的耐受性。

根据第二方面，一种实施例中提供一种基于耐受性分析的基因等级排序装置，包括：

数据库数据获取单元，用于获取频率数据库和转录本数据库的数据，其中上述频率数据库包括多个样本的测序数据，上述转录本数据库包括多个转录本的测序数据；

编码序列数据获取单元，用于从上述转录本数据库获取各个基因的转录本作为编码序列数据，上述编码序列数据包括每个外显子的两端碱基；

目标基因序列提取单元，用于从上述频率数据库中提取目标基因序列，上述目标基因序列是测序深度达到第一预设值的碱基位点的覆盖度达到第二预设值的基因序列；

突变检测单元，用于将上述编码序列数据映射到上述目标基因序列得到上述两端碱基之间的序列，检测得到上述目标基因序列的上述两端碱基之间序列上的功能性突变和非功能性突变；

分析对象选取单元，用于选取上述功能性突变和非功能性突变中等位基因频率在第三预设值以下的突变作为基因等级排序的分析对象；

基因等级分数获取单元，用于对上述分析对象建立基因等级排序表，得到多个不同基因的基因等级分数，上述基因等级分数表征上述基因的耐受性。

根据第三方面，一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如第一方面的方法。

本发明的基于耐受性分析的基因等级排序方法，以基因为单位，利用功能基因对疾病的耐受性，建立基因等级排序，在功能基因的预测上具有较高的可信度和准确性，为孟德尔遗传病候选基因的筛选提供重要参考。

附图说明

图1为本发明实施例的基于耐受性分析的基因等级排序方法流程示意图；

图2为本发明实施例中一个具有3个CCDS转录本的基因ATP1A3根据坐标融合CCDS边界的结果图；

图3为本发明实施例中功能性突变x和突变的总数y的回归作图结果；

图4为本发明实施例中1至10号基因对应的标准残差结果图；

图5为本发明实施例的基于耐受性分析的基因等级排序装置结构框图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。在以下的实施方式中，很多细节描述是为了使得本发明能被更好的理解。然而，本领域技术人员可以毫不费力的认识到，其中部分特征在不同情况下是可以省略的，或者可以由其他材料、方法所替代。

另外，说明书中所描述的特点、操作或者特征可以以任意适当的方式结合形成各种实施方式。同时，方法描述中的各步骤或者动作也可以按照本领域技术人员所能显而易见的方式进行顺序调换或调整。因此，说明书和附图中的各种顺序只是为了清楚描述某一个实施例，并不意味着是必须的顺序，除非另有说明其中某个顺序是必须遵循的。

本文中为部件所编序号本身，例如“第一”、“第二”等，仅用于区分所描述的对象，不具有任何顺序或技术含义。

如图1所示，本发明实施例的基于耐受性分析的基因等级排序方法，包括：

S101：获取频率数据库和转录本数据库的数据，其中上述频率数据库包括多个样本的测序数据，上述转录本数据库包括多个转录本的测序数据。

本发明实施例中，频率数据库可以是全基因组测序(WGS)数据等，这样的数据包括多个样本的测序数据，例如外显子组整合数据库(the Exome Aggregation Consortium，ExAC)(http://exac.broadinstitute.org/)，其数据来源很丰富，包括6000多个个体的测序数据，人数多，各个碱基位点的质量值和测序深度都有提供，使用比较方便、准确性比较高。类似的数据库还有1000Genomes(http://www.internationalgenome.org/)等。转录本数据库可以是全外显子测序(WES)数据等，这样的数据包括多个转录本的测序数据，例如一致性编码序列数据库(Consensus coding sequence(CCDS)database)(https://www.ncbi.nlm.nih.gov/projects/CCDS/CcdsBrowse.cgi)，其是NCBI旗下的转录本数据库，CCDS数据都比较准确。类似的转录本数据库还包括hg19refGene(http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz),但是hg19refGene会有一些错误的剪切位点。

S102：从上述转录本数据库获取各个基因的转录本作为编码序列数据，上述编码序列数据包括每个外显子的两端碱基。

本发明实施例中，转录本数据库例如CCDS数据都比较准确，而且是外显子序列，正是感兴趣的突变所在的序列，因此从转录本数据库获取各个基因的转录本作为编码序列数据，该步骤针对数据库中每一个感兴趣的基因进行分析，每个基因可能有一个或多个转录本。

在基因具有一个转录本的情况下，对每个外显子的两端进行两个碱基的延长得到每个外显子的两端碱基，两个碱基的延长用于提供剪接受体和供体位点。这样经过两个碱基的延长得到的序列数据作为编码序列数据，编码序列数据包括每个外显子的两端碱基，即两个碱基的延长以后两个端点的碱基位置信息，该碱基位置信息作为转录本的边界，用于在后续步骤S104中对应到目标基因序列上，从而界定目标基因序列上有待分析突变情况的序列范围。之所以进行两个碱基的延长，原因如下：DNA在转录成RNA时，可选择性的剪切，将DNA“外显子”剪下来，“剪刀”识别的部分就是剪切受体和供体位点，剪切受体和供体位点的突变可能导致无法识别，产生突变的转录本，因此剪切受体和供体位点也是需要分析的位点，两个碱基的延长能够实现对剪切受体和供体位点的覆盖。

在基因具有多个转录本的情况下，还包括：将该基因的所有转录本融合成为一个转录本作为编码序列数据。在一些实施例中，对于Gene Nomenclature Committee(HGNC)上拥有多个CCDS转录本的基因，将所有转录本融合作为一个CCDS边界。HGNC提供官方认可的基因名(https://www.genenames.org/download/statistics-and-files/)，使用统一的基因名称，可以很好地去冗余。比如，某个标准的基因名称是WDR4，在数据库HGNC的记录是，HGNC:12756→WDR4→WD repeat domain 4→Approved→TRM82,TRMT82→21q22.3；在CCDS里面记录转录本的信息是，21→NC_000021.9→WDR4→10785→CCDS13691.1→Public→-→42850048→42879494→[42850048-42850241,42852254-42852323,42853568-42853751,42854561-42854625,42855681-42855779,42859661-42859721,42862281-42862393,42863439-42863595,42873550-42873690,42876701-42876766,42879406-42879494]→Identical。HGNC和CCDS可以互相用基因名字来索引提取。例如，图2示出了一个具有3个CCDS转录本的基因ATP1A3，根据坐标融合CCDS边界。

S103：从上述频率数据库中提取目标基因序列，上述目标基因序列是测序深度达到第一预设值的碱基位点的覆盖度达到第二预设值的基因序列。

目标基因序列需要满足至少两个条件，即测序深度和覆盖度，测序深度是指每个碱基位置被测到的次数，而覆盖度是指该基因序列上碱基位置达到设定的测序深度(第一预设值)的碱基数量占该基因序列总碱基数量的比例，如果覆盖度达到设定的数值(第二预设值)，那么该基因序列就是目标基因序列。在一个实施例中，第一预设值是10倍测序深度，即位点深度要求最低为10X；第二预设值是70％，即基因序列上至少70％以上的位点为10X以上的测序深度。通过该步骤，低质量数据被过滤掉，在一个实施例中，低质量数据即基因序列上达到10倍测序深度的位点不足70％的基因序列，这些基因序列数据被剔除掉。

需要说明的是，上述步骤S102和步骤S103的顺序没有特别限定，说明书和权利要求书中虽然步骤S102出现在步骤S103之前，但应当理解为先执行步骤S102再执行步骤S103，或先执行步骤S103再执行步骤S102，或步骤S102和步骤S103同时执行。

S104：将上述编码序列数据映射到上述目标基因序列得到上述两端碱基之间的序列，检测得到上述目标基因序列的上述两端碱基之间序列上的功能性突变和非功能性突变。

以在步骤S102中得到的编码序列数据的每个外显子的两端碱基的位置作为界限，能够找到目标基因序列上对应于上述两端碱基的位置之间的序列区域，作为感兴趣的突变分析区域，在该区域内分析功能性突变和非功能性突变。功能性突变包括错义突变、无义突变和剪切突变；上述非功能性突变包括同义突变。其中，错义突变是指由于DNA链上的碱基替换改变了信使RNA上特定的遗传密码，并引起合成的多肽链中的一个氨基酸被另一个氨基酸取代。无义突变是指由于单个碱基的替换引起出现了终止密码子，从而提前终止了多肽链的合成，产生的蛋白大都失去了活性或丧失了正常的功能。剪切突变是指上述提到的外显子两端两个碱基作为转录剪切的识别位点发生突变，不能识别，导致错误的剪切，产生不正常的多肽链。同义突变是指单个碱基的替换可能只改变了信使RNA上特定的密码子，但由于密码子具有简并性，因此并不影响氨基酸的正常编码。前三种突变都会影响氨基酸的正常编码，可能会导致这个基因丧失功能，称为功能性突变。最后一种突变不会影响氨基酸的正常编码，所以是非功能性突变。

S105：选取上述功能性突变和非功能性突变中等位基因频率在第三预设值以下的突变作为基因等级排序的分析对象。

由于功能性突变和非功能性突变中有大量突变的突变频率非常高，这样的突变一般都是无害的，不是本发明分析的对象，因此需要通过一定的方法和标准去除这样的突变。在一个实施例中，通过所谓“最小等位频率(MAF)的选取”，即选取等位基因频率在一定的预设值(第三预设值)以下的突变作为基因等级排序的分析对象。在一些实施例中，第三预设值是1％以下的数值，优选0.01％以下的数值，更优选0.01％。一般而言，第三预设值越小得到的分析对象越少，这部分突变更可能是有害的突变，而第三预设值越大得到的分析对象越多，得到的突变是有害的突变的准确性越小。但是，第三预设值过小可能导致分析对象过少，而第三预设值过大可能导致准确性降低。

例如，DMD基因在ExAC数据库中的突变如下链接中显示的信息所示：

http://exac.broadinstitute.org/gene/ENSG00000198947，其中最后一列就是等位基因频率，根据这一列的频率和设定的预设值(例如0.01％)选取分析对象，即低于预设值的突变。

S106：对上述分析对象建立基因等级排序表，得到多个不同基因的基因等级分数，上述基因等级分数表征上述基因的耐受性。

在本发明的一个实施例中，对分析对象建立基因等级排序表，通过如下方法实现：

以功能性突变的总数作为x轴，突变的总数作为y轴作图，其中突变的总数包括功能性突变和非功能性突变的总数；

x对y回归，将得到的标准残差作为基因等级分数来表征基因的耐受性，该基因等级分数可用于评估特定基因与单基因遗传特征的关系或影响大小。

基因等级分数小于0的为不耐受基因，且随着分数的降低，耐受性降低。基因等级分数大于0的为耐受基因，且随着分数的升高耐受性增高。

本发明实施例中，所谓“残差”是指观测值与预测值(拟合值)之间的差，即实际观察值与回归估计值的差，每个基因都有一个对应的残差和标准残差，以此评估该基因的耐受性。

例如，图3示出了一个实施例中功能性突变x和突变的总数y的回归作图结果，其中每一个圈表示一个基因的观测值，曲线表示预测值(拟合值)，每一个圈到曲线的距离表示残差。图4中示出了1至10号共10个基因对应的标准残差。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过硬件的方式实现，也可以通过计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘、光盘、硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。

对应于本发明实施例的基于耐受性分析的基因等级排序方法，本发明还提供一种基于耐受性分析的基因等级排序装置，如图5所示，包括：数据库数据获取单元501，用于获取频率数据库和转录本数据库的数据，其中上述频率数据库包括多个样本的测序数据，上述转录本数据库包括多个转录本的测序数据；编码序列数据获取单元502，用于从上述转录本数据库获取各个基因的转录本作为编码序列数据，上述编码序列数据包括每个外显子的两端碱基；目标基因序列提取单元503，用于从上述频率数据库中提取目标基因序列，上述目标基因序列是测序深度达到第一预设值的碱基位点的覆盖度达到第二预设值的基因序列；突变检测单元504，用于将上述编码序列数据映射到上述目标基因序列得到上述两端碱基之间的序列，检测得到上述目标基因序列的上述两端碱基之间序列上的功能性突变和非功能性突变；分析对象选取单元505，用于选取上述功能性突变和非功能性突变中等位基因频率在第三预设值以下的突变作为基因等级排序的分析对象；基因等级分数获取单元506，用于对上述分析对象建立基因等级排序表，得到多个不同基因的基因等级分数，上述基因等级分数表征上述基因的耐受性。

本发明的一种实施例中提供一种计算机可读存储介质，包括程序，该程序能够被处理器执行以实现如本发明实施例的基于耐受性分析的基因等级排序方法。

以上应用了具体个例对本发明进行阐述，只是用于帮助理解本发明，并不用以限制本发明。对于本发明所属技术领域的技术人员，依据本发明的思想，还可以做出若干简单推演、变形或替换。

Claims (11)

1.一种基于耐受性分析的基因等级排序方法，其特征在于，所述方法包括：

获取频率数据库和转录本数据库的数据，其中所述频率数据库包括多个样本的测序数据，所述转录本数据库包括多个转录本的测序数据；

从所述转录本数据库获取各个基因的转录本作为编码序列数据，所述编码序列数据包括每个外显子的两端碱基；

从所述频率数据库中提取目标基因序列，所述目标基因序列是测序深度达到第一预设值的碱基位点的覆盖度达到第二预设值的基因序列；

将所述编码序列数据映射到所述目标基因序列得到所述两端碱基之间的序列，检测得到所述目标基因序列的所述两端碱基之间序列上的功能性突变和非功能性突变；

选取所述功能性突变和非功能性突变中等位基因频率在第三预设值以下的突变作为基因等级排序的分析对象；

对所述分析对象建立基因等级排序表，得到多个不同基因的基因等级分数，所述基因等级分数表征所述基因的耐受性。

2.根据权利要求1所述的方法，其特征在于，所述频率数据库是ExAC数据库；所述转录本数据库是CCDS数据库。

3.根据权利要求1所述的方法，其特征在于，所述从所述转录本数据库获取转录本作为编码序列数据，包括：对每个外显子的两端进行两个碱基的延长得到所述每个外显子的两端碱基，所述两个碱基的延长用于提供剪接受体和供体位点。

4.根据权利要求3所述的方法，其特征在于，所述基因包括多个所述转录本，所述从所述转录本数据库获取转录本作为编码序列数据，还包括：将该基因的所有转录本融合成为一个转录本作为所述编码序列数据。

5.根据权利要求1所述的方法，其特征在于，所述第一预设值是10倍深度，所述第二预设值是70％。

6.根据权利要求1所述的方法，其特征在于，所述功能性突变包括错义突变、无义突变和剪切突变；所述非功能性突变包括同义突变。

7.根据权利要求1所述的方法，其特征在于，所述第三预设值是1％以下的数值。

8.根据权利要求7所述的方法，其特征在于，所述第三预设值是0.01％以下的数值。

9.根据权利要求1所述的方法，其特征在于，所述对所述分析对象建立基因等级排序表，包括：

以所述功能性突变的总数作为x轴，突变的总数作为y轴作图，其中所述突变的总数包括功能性突变和非功能性突变的总数；

所述x对所述y回归，将得到的标准残差作为基因等级分数来表征所述基因的耐受性。

10.一种基于耐受性分析的基因等级排序装置，其特征在于，所述装置包括：

数据库数据获取单元，用于获取频率数据库和转录本数据库的数据，其中所述频率数据库包括多个样本的测序数据，所述转录本数据库包括多个转录本的测序数据；

编码序列数据获取单元，用于从所述转录本数据库获取各个基因的转录本作为编码序列数据，所述编码序列数据包括每个外显子的两端碱基；

目标基因序列提取单元，用于从所述频率数据库中提取目标基因序列，所述目标基因序列是测序深度达到第一预设值的碱基位点的覆盖度达到第二预设值的基因序列；

突变检测单元，用于将所述编码序列数据映射到所述目标基因序列得到所述两端碱基之间的序列，检测得到所述目标基因序列的所述两端碱基之间序列上的功能性突变和非功能性突变；

分析对象选取单元，用于选取所述功能性突变和非功能性突变中等位基因频率在第三预设值以下的突变作为基因等级排序的分析对象；

基因等级分数获取单元，用于对所述分析对象建立基因等级排序表，得到多个不同基因的基因等级分数，所述基因等级分数表征所述基因的耐受性。

11.一种计算机可读存储介质，其特征在于，包括程序，所述程序能够被处理器执行以实现如权利要求1-9中任一项所述的方法。

Images