阅读说明：本技术 一种变异序列的注释方法 (Variant sequence annotation method ) 是由 文文 王红阳 朱赢 陈淑桢 何慧斯 高勇 汪德鹏 于 2020-05-25 设计创作，主要内容包括：本发明属于生物信息技术领域,具体涉及一种变异序列注释方法,所述方法包括：(1)确定变异序列信息：获得变异序列,整合参考序列信息,标准化变异信息；(2)变异注释,注释结果包括注释功能区域、变异类型、核酸序列、氨基酸序列。该方法不仅能实现行业金标准ANNOVAR的现有功能,而且克服了ANNOVAR中的缺点,在区别剪接位点和剪接区域变异、CDS边缘变异、注释frameshift和stoploss/stopgain等方面进行了完善,而且使用了规范的表示方式,还增加了基因编号Entrez ID,具有更好的应用价值。(The invention belongs to the technical field of biological information, and particularly relates to a variant sequence annotation method, which comprises the following steps: (1) determining variant sequence information: obtaining variant sequences, integrating reference sequence information and standardizing variant information; (2) and (4) variant annotation, wherein the annotation result comprises an annotated functional region, variant types, nucleic acid sequences and amino acid sequences. The method can not only realize the existing functions of the ANNOVAR of the trade golden standard, but also overcome the defects in the ANNOVAR, is perfect in distinguishing the variation of splicing sites and splicing regions, CDS edge variation, annotation of frameshift, stoppages/stopgain and the like, uses a standard representation mode, increases the gene number Entrez ID, and has better application value.)

一种变异序列的注释方法

技术领域

本发明属于生物信息技术领域，具体涉及一种变异序列的注释方法。

背景技术

随着测序技术的发展，测序通量不断上升、测序成本持续下降，有越来越多的物种已获得基因组、转录组信息。在细分领域，有越来越多的研究关注在同一物种的不同品种或群体，乃至差异化个体之间的变异，以寻求大的遗传背景下，个别遗传信息的变异所导致的表型差异。这就对变异序列的查找及注释提出了挑战。

以人类为例，ANNOVAR是现在主流的对变异进行注释的软件，被行业内认为是金标准，但是在实际使用时，发明人发现ANNOVAR未能解决以下几个问题：

在转录本形成过程中，mRNA前体中通过不同剪接方式，选择不同的剪接位点组合，产生不同的剪接异构体；其中，剪接位点即为对应元件的边缘。行业内共识认为，处于剪接位点±2bp以内的变异会对基因的剪接产生影响。但已有不少研究表明，在剪接位点邻近区域、处于剪接位点±2bp以外的变异也会对基因的剪接产生影响。也就是说，对紧邻剪接位点处的变异、剪接位点附近的变异进行区分注释，是更加科学合理的。但是，ANNOVAR仅是笼统的进行了剪接区域的注释，并未区分。

另外，研究表明：除了处于剪接区域内的变异对基因的剪接产生影响，处于CDS边缘的变异也会对基因的剪接产生影响。ANNVOAR未对此类位点进行特定的注释或标记。

对于某些InDel变异，变异类型同时出现frameshift和stoploss/stopgain时，ANNOVAR会丢失frameshift或stoploss/stopgain其中一个信息，造成注释信息缺失。

在基因的后续研究中，由于基因命名规则的特点，同一基因的基因名称(symbol)经常发生改变，这导致在不同版本的数据库的注释下同一变异所注释到的基因名称不同。目前业内很多权威数据库如NCBI、OMIM等已经开始带入gene的Entrez ID对基因名称加以标记，以保证注释结果的唯一性。

人类基因组变异协会HGVS(Human Genome Variation Society)制定了目前学术界所公认的突变命名规则(http://varnomen.hgvs.org/)，但ANNOVAR默认未使用HGVS规范命名。同时，在蛋白命名的规则上，HGVS建议使用氨基酸三字母简写如p.Arg727Ser，而ANNOVAR使用了不符合规范建议的氨基酸一字母简写如p.R727S。

发明内容

有鉴于此，本发明提供了一种变异序列注释的方法，该方法不仅能实现ANNOVAR现有功能，还实现了区别剪接位点和剪接区域变异，新增CDS边缘变异注释，完整包含frameshift和stoploss/stopgain信息等功能，而且具有更规范的输出形式。

为了实现上述目的，本发明的技术方案具体如下：

一种序列变异注释方法，包括以下步骤：

(1)确定变异序列信息

(1.1)获得变异序列：

使用变异分析软件，将待分析序列与参考基因组进行比对，获得变异信息；

(1.2)整合参考序列信息：

获取参考基因组序列以及参考基因组注释文件；根据注释文件，从参考基因组序列中提取全部的参考基因组转录本和CDS序列；

根据基因的描述信息，获取参考基因组中基因对应的Entrez ID信息；

将参考基因组转录本和CDS序列、参考基因组注释文件、Entrez ID信息进行整合，获得整合参考基因组信息；

(1.3)标准化变异信息

对步骤(1.1)中获取的每一个变异信息，提取每个变异的染色体信息、参考基因组物理位置、参考基因组序列、变异序列信息，进行标准化处理，获得标准化变异信息；

所述标准化变异信息包括：染色体信息、起始位置、终止位置、标准化参考基因组序列、标准化变异序列；

(2)变异注释

(2.1)注释功能区域

根据标准化变异信息判断变异与元件相对位置，具体分为：变异位于元件边缘、变异位于元件深处；所述位于元件边缘是指所述起始位置或终止位置距离元件邻近边缘长度小于等于xbp，所述位于元件深处是所述起始位置或终止位置距离元件邻近边缘长度大于xbp；需要说明的是，由于每个元件都会有两个边缘，相当于该元件的开始位置和终止位置，在进行比较时，所述边缘指的是两个边缘中相对更近的那个边缘。

当起始位置或终止位置位于元件边缘时，进一步区分边缘位点和边缘区域；所述边缘位点是指所述起始位置或终止位置在元件邻近边缘±ybp，所述边缘区域是指所述起始位置或终止位置在元件邻近边缘-ybp至-xbp或+ybp至+xbp的区域，所述y小于x；

所述元件包括UTR、CDS和Intron；

(2.2)注释变异类型

若一个变异的起始、终止位置均位于非CDS区域，注释为空；

若一个变异的起始位置和/或终止位置位于CDS区域内，则将参考cDNA序列翻译为参考氨基酸序列，将参考cDNA中的碱基替换为变异碱基得到变异cDNA序列，并翻译为变异氨基酸序列；然后通过比对参考cDNA序列和变异cDNA序列、参考氨基酸序列和变异氨基酸序列，按照单碱基变异、***变异和删除变异对变异类型进行分类注释；

(2.3)注释核酸序列变异

对比参考cDNA序列与变异cDNA序列，根据HGVS规则注释变异cDNA序列的核酸变异信息；

(2.4)注释氨基酸序列变异

对比参考氨基酸序列与变异氨基酸序列，根据HGVS规则注释变异氨基酸序列的氨基酸变异信息，其中氨基酸使用三字母缩写表示。

在上述技术方案中，步骤(1.2)所述提取参考基因组转录本和CDS序列的方法为：根据参考基因组注释文件中每个转录本的物理位置信息，以染色体为单位，从参考基因组序列中提取全部的参考基因组转录本和CDS序列；或一次性读取全部的参考基因组序列，然后根据参考基因组注释文件中每个转录本的物理位置，提取参考基因组转录本和CDS序列；对比两种方案，第一种提取方法消耗的内存资源少，速度更快。

在上述技术方案中，对于步骤(1.2)所述的整合参考基因组信息，建立信息索引，具体方法为：以染色体为单位，以一定步长，将参考基因组切割为若干窗口，根据参考基因组注释文件中的物理位置信息，获取每一个窗口所包含的转录本信息；更进一步地，所述步长为300kb。建立索引更便于快速调取信息，步长的大小直接影响索引数目、计算机运算速度和内存。

在上述技术方案中，步骤(1.3)所述标准化处理方法如下：

当参考基因组序列与变异序列的长度同时等于1时，起始位置＝终止位置＝参考基因组物理位置；

当参考基因组序列与变异序列的长度不同或相同但不等于1时，删除两者中相同的碱基，所删除的参考基因组序列的左侧碱基长度记为LEN，起始位置和终止位置的确定方式如下：

当标准化参考基因组序列长度为0时，起始位置＝参考基因组物理位置+LEN-1；当标准化参考基因组序列长度大于0时，起始位置＝参考基因组物理位置+LEN；

当标准化参考基因组序列长度小于等于1时，终止位置＝起始位置；当标准化参考基因组序列长度大于1时，终止位置＝起始位置+标准化参考基因组序列长度-1。

在上述技术方案中，注释功能区域的目的在于确定变异序列是位于基因的哪个功能区，步骤(2.1)所述注释功能区域的方法具体为：

a.变异位于上游或下游，注释为UpStream或DownStream；

b.变异位于UTR元件，

变异位于UTR深处，注释为UTR3或UTR5；

变异位于UTR边缘，且与该边缘相邻的元件为非Intron区域，则注释为UTR3或UTR5；

变异位于UTR边缘，且该与边缘相邻的元件为Intron：若距离边缘长度小于或等于y，注释为UTR3_splicing_site或UTR5_splicing_site；若距离边缘长度大于y小于x,注释为UTR3_splicing_region或UTR5_splicing_region；

c.变异位于CDS元件，

变异位于CDS深处，注释为exonic；

变异位于CDS边缘，且与该边缘相邻的元件为非Intron区域，则注释为exonic；

变异位于CDS边缘，且该与边缘相邻的元件为Intron：若距离边缘长度小于或等于y，注释为CDS_splicing_site；若距离边缘长度大于y小于x，注释为CDS_splicing_region；

d.变异位于Intron元件，

变异位于Intron深处，注释为intronic；

变异位于Intron边缘：若距离边缘长度小于或等于y，注释为splicing_site；若距离边缘长度大于y小于x，注释为splicing_region；

变异跨过Intron与相邻元件的连接点，注释为splicing_site；

所述变异为标准化变异信息中变异的起始位置或终止位置。

在上述技术方案中，所述x＝10，y＝2。

在上述技术方案中，步骤(2.2)所述注释变异类型的方法具体为：

a.对于单碱基变异

若参考氨基酸序列与变异氨基酸序列相同，注释为synonymous_snv

若参考氨基酸序列与变异氨基酸序列不同，注释为nonsynonymous_snv

b.对于***变异

当变异cDNA序列与参考cDNA序列长度之差为3的倍数时，比较变异cDNA序列与参考cDNA序列终止密码子的位置：若变异cDNA序列中提前出现终止密码子，注释为ins_nonframeshift_stopgain；若变异cDNA序列中终止密码子消失，注释为ins_nonframeshift_stoploss；若变异cDNA序列的终止密码正常出现在末端，设置变异类型为ins_nonframeshift；

当变异cDNA序列与参考cDNA序列长度之差不为3的倍数时，比较变异cDNA序列与参考cDNA序列终止密码子的位置：若变异cDNA序列提前出现终止密码子，注释为ins_frameshift_stopgain；若变异cDNA序列中终止密码子消失，注释为ins_frameshift_stoploss；若变异cDNA序列的终止密码正常出现在末端，注释为ins_frameshift；

c.对于删除变异

当变异cDNA序列与参考cDNA序列长度之差为3的倍数时，比较变异cDNA序列与参考cDNA序列终止密码子的位置：若变异cDNA序列中提前出现终止密码子，注释为del_nonframeshift_stopgain；若变异cDNA序列中终止密码子消失，注释为del_nonframeshift_stoploss；若变异cDNA序列的终止密码正常出现在末端，注释为del_nonframeshift；

当变异cDNA序列与参考cDNA序列长度之差不为3的倍数时，比较变异cDNA序列与参考cDNA序列终止密码子的位置：若变异cDNA序列提前出现终止密码子，注释为del_frameshift_stopgain；若变异cDNA序列中终止密码子消失，注释为del_frameshift_stoploss；若变异cDNA序列的终止密码正常出现在末端，注释为del_frameshift。

本发明的有益效果为：本发明提供的方法与行业金标准ANNOVAR相比，不仅注释数目、大的分类一致，而且克服了ANNOVAR中的缺点，在区别剪接位点和剪接区域变异、CDS边缘变异、frameshift和stoploss/stopgain缺失等方面进行了科学、细致的分类，而且使用了规范的表示方式，还增加了基因编号Entrez ID。

具体实施方式

为了更好的理解本发明，下面结合实施例对本发明做进一步的详细说明。

实施例

(1)确定变异序列信息(call variants)

(1.1)获得变异序列

使用探针捕获技术，对人全外显子进行二代测序，获得待分析序列。使用变异序列分析软件(GATK，https://gatk.broadinstitute.org/hc/en-us)，将待分析序列与参考基因组进行比对，获取变异信息(call variants)。

(1.2)整合参考序列信息

获取hg19参考基因组序列和hg19参考基因组注释文件，该注释文件包括基因名称、转录本名称、物理位置、正负链、每个元件信息(元件包括UTR、Intron、CDS)等。其中，hg19参考基因组序列的下载地址为ftp://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz；hg19参考基因组注释文件的获取地址为：ftp://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz。

根据参考基因组注释文件中每个转录本的物理位置信息，以染色体为单位，从参考基因组序列中提取全部的参考基因组转录本和CDS序列。该提取方法消耗的内存资源少，速度快。

根据基因的描述信息获取参考基因组中每个基因对应的Entrez ID信息，其目的在于：使用统一的编号信息，避免不规范命名带来的混淆。

上述基因描述信息的下载地址为：ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz。

将参考基因组转录本和CDS序列、参考基因组注释文件、Entrez ID信息进行整合，获得整合参考基因组信息。

为了达到快速调取信息的目的，建立了参考基因组转录本信息索引，具体方法为：以染色体为单位，以300kb为步长将参考基因组的序列切割为若干窗口，根据参考基因组注释文件中的物理位置信息，获取每一个窗口所包含的转录本信息，作为索引。测试发现，步长越小，索引越多，计算速度越慢；步长越大，索引越少，计算所需内存越大；本实施例中采用300kb为步长，是平衡计算速度、计算内存后的最佳值。

(1.3)标准化变异信息

对步骤(1.1)中获取的每一个变异信息，提取每个变异的染色体信息(CHROM)、参考基因组物理位置(POS)、参考基因组序列(REF)、变异序列(ALT)，进行标准化处理，获得标准化变异信息。

标准化后处理后的每个变异的信息包括：染色体信息(CHROM)、起始位置(START)、终止位置(END)、标准化的参考基因组序列(REF)、标准化的变异序列(ALT)；

标准化处理前，当REF与ALT的长度同时等于1时，即此该变异为单碱基变异(SNP)，则标准化信息中的START和END与标准化前的POS数值相同；例如下表所示(“/”表示对应项无信息)：

标准化处理前，当REF与ALT的长度不同、或相同但不等于1时，即此变异为***或删除变异(InDel)，删除两者中相同的碱基，删除的REF左侧碱基长度记为LEN；标准化处理后，标准化的REF的长度记为LEN_REF，标准化的ALT的长度记为LEN_ALT，此时，START和END的确定方式如下：

START确定方式：若LEN_REF＝0，则START＝POS+LEN-1；若LEN_REF>0，则START＝POS+LEN；

END确定方式：若LEN_REF≤1，则END＝START；若LEN_REF>1，则END＝START+LEN_REF-1。

若LEN_REF为0，则标准化的REF为“-”；若LEN_ALT为0，则标准化的ALT为“-”；例如下表所示(“/”表示对应项无信息)：

(2)变异注释

注释结果包括注释功能区域、变异类型、核酸序列和氨基酸序列四种，具体步骤如下：

(2.1)注释功能区域

首先，需要明确的是：变异落在元件边缘是指该变异的起始位置或终止位置距离元件边缘长度小于或等于10bp，变异落在元件深处是指该变异的起始位置或终止位置距离元件边缘长度大于10bp；当位于元件边缘时，进一步区分边缘位点和边缘区域；边缘位点(splicing_site)是指所述的起始位置或终止位置在元件边缘±2bp的区域，即剪接位点±2bp的区域；边缘区域(splicing_region)是指所述起始位置或终止位置在元件边缘-ybp至-xbp或+ybp至+xbp的区域剪接区域即剪接位点-2bp至-10bp、+2bp至+10bp的区域；由于每个元件都有两个边缘，本步骤中所述的边缘是个两个边缘中相对更近的那个边缘。

在上述规定的基础上，具体分为以下情形(本步骤中所述的变异是指每个变异的起始位置或终止位置)：

a.变异位于上游或下游，注释为UpStream或DownStream。这里所述的上游或下游具体指非元件区，即不属于UTR、CDS和Intron的都归在上下游内。

b.变异位于UTR元件：

变异位于UTR深处，注释为UTR3或UTR5；

变异位于UTR边缘，且与该边缘相邻的元件为非Intron区域，则注释为UTR3或UTR5；

变异位于UTR边缘，且该与边缘相邻的元件为Intron：若距离边缘长度小于或等于2，注释为UTR3_splicing_site或UTR5_splicing_site；若距离边缘长度大于2小于10,注释为UTR3_splicing_region或UTR5_splicing_region。

c.变异位于CDS元件：

变异位于CDS深处，注释为exonic；

变异位于CDS边缘，且与该边缘相邻的元件为非Intron区域，则注释为exonic；

变异位于CDS边缘，且该与边缘相邻的元件为Intron：若距离边缘长度小于或等于2，注释为CDS_splicing_site；若距离边缘长度大于2小于10，注释为CDS_splicing_region。

d.变异位于Intron元件：

变异位于Intron深处，注释为intronic；

变异位于Intron边缘：若距离边缘长度小于或等于2，注释为splicing_site；若距离边缘长度大于2小于10，注释为splicing_region；

变异跨过Intron与相邻元件的连接点，注释为splicing_site。

(2.2)注释变异类型

当一个变异的起始、终止位置均位于非CDS区域，注释为空；

当一个变异的起始位置和/或终止位置位于CDS区域内，则将参考cDNA序列(NS1)翻译为参考氨基酸序列(AS1)，将参考cDNA中的碱基替换为变异后的碱基得到变异cDNA序列(NS2)，并翻译为变异氨基酸序列(AS2)；

a.对于单碱基变异

若AS1＝AS2，注释为synonymous_snv

若AS1≠AS2，注释为nonsynonymous_snv

b.对于***变异

当NS2与NS1长度之差为3的倍数时，比较NS2与NS1终止密码子的位置：若NS2中提前出现终止密码子，注释为ins_nonframeshift_stopgain；若NS2中终止密码子消失，注释为ins_nonframeshift_stoploss；若NS2的终止密码正常出现在末端，设置变异类型为ins_nonframeshift；

当NS2与NS1长度之差不为3的倍数时，比较NS2与NS1终止密码子的位置：若NS2提前出现终止密码子，注释为ins_frameshift_stopgain；若NS2中终止密码子消失，注释为ins_frameshift_stoploss；若NS2的终止密码正常出现在末端，注释为ins_frameshift。

c.对于删除变异

当NS2与NS1长度之差为3的倍数时，比较NS2与NS1终止密码子的位置：若NS2中提前出现终止密码子，注释为del_nonframeshift_stopgain；若NS2中终止密码子消失，注释为del_nonframeshift_stoploss；若NS2的终止密码正常出现在末端，注释为del_nonframeshift；

当NS2与NS1长度之差不为3的倍数时，比较NS2与NS1终止密码子的位置：若NS2提前出现终止密码子，注释为del_frameshift_stopgain；若NS2中终止密码子消失，注释为del_frameshift_stoploss；若NS2的终止密码正常出现在末端，注释为del_frameshift。

(2.3)注释核酸序列变异

对比NS1与NS2，根据HGVS规则设置NS2的核酸变异信息。

(2.4)注释氨基酸序列变异

对比AS1与AS2，根据HGVS规则设置AS2的氨基酸变异信息，其中氨基酸使用三字母缩写形式。

对比例

采用实施例步骤(1.1)中获得的变异序列数据，使用ANNOVAR软件进行变异位点注释，作为对比例。

结果分析：

实施例和对比例获得的注释结果均为119,586条，且在变异位点、变异类型上的判断全部一致，说明本发明的注释方法是准确可信的。但是，实施例实现了更为细致的注释，且采取更为规范的表现方式，克服了ANNOVAR软件中的诸多问题。具体见一下示例：

示例1

变异位点：1号染色体第69511位A突变为G(步骤1.3标准化变异信息，1:69511:69511:A:G)

氨基酸序列变异注释结果：

对比例：OR4F5:NM_001005484:exon1:c.421A>G:p.T141A

实施例：OR4F5:NM_001005484:exon1:c.421A>G:p.Thr141Ala

注释结果一致，但对比例中ANNOVAR使用氨基酸一字母简写不符合规范。

核酸序列变异注释结果：

对比例：Symbol：OR4F5

实施例：Symbol：OR4F5，EntrezID：79501

注释结果一致，但对比例中ANNOVAR无EntrezID信息。

功能区域注释：均为exonic；结果一致。

变异类型注释：均为nonsynonymous SNV；结果一致。

示例2

变异位点：第9号染色体70176769位G缺失(步骤1.3标准化变异信息，9:70176769:70176769:G:-)

核酸及氨基酸序列变异注释结果：

对比例：FOXD4L5:NM_001126334:exon1:c.1215delC:p.W406Gfs*21

实施例：FOXD4L5:NM_001126334:exon1:c.1215_1215del:p.Trp406Glyfs

注释结果一致，但实施例中显示缺失起始及终止位点，氨基酸使用三字母。

功能区域注释：均为exonic；结果一致。

变异类型注释：

对比例：frameshift deletion

实施例：del_frameshift_stoploss

实施例中成功注释stoploss信息。

基因信息：

对比例：Symbol：FOXD4L5

实施例：Symbol：FOXD4L5，EntrezID：653427

对比例中ANNOVAR无EntrezID信息。

示例3

变异位点：第1号染色体883625位A突变为G(步骤1.3标准化变异信息，1:883625:883625:A:G)

核酸序列注释结果：

对比例：NM_015658:exon14:c.1558-13T>C

实施例：NOC2L:NM_015658:exon14:c.1558-13T>C

实施例中能提供基因名，对比例无此功能。

功能区域注释：

对比例：splicing

实施例：splicing_region

实施例成功注释剪切区域变异信息。

变异类型注释：均为空；一致。

基因信息：

对比例：Symbol：NOC2L

实施例：Symbol：NOC2L，EntrezID：26155

示例4

变异位点：第4号染色体190874281位G突变为T(步骤1.3标准化变异信息，4:190874281:190874281:G:T)

核酸序列注释结果：

对比例：NM_004477:exon4:c.317+1G>T

实施例：FRG1:NM_004477:exon4:c.317+1G>T

功能区域注释：

对比例：splicing

实施例：splicing_site

实施例中成功注释剪切位点变异信息。

基因信息：

对比例：Symbol：FRG1

实施例：Symbol：FRG1，EntrezID：2483。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。