一种基于基因组环境确定遗传变异功能影响的方法与流程

本发明涉及生物信息学领域，更具体地涉及一种基于基因组环境确定遗传变异功能影响的方法。

背景技术：

随着以深度测序为代表的高通量遗传变异检测技术的快速发展，目前已可以快速鉴定个体基因组上的遗传变异。然而，如何准确确定这些遗传变异对生物分子功能的影响，从而为后续的个性化医疗、分子育种等应用提供线索、指导与支持，仍是目前该领域面临的重大挑战。

目前在变异注释领域常用的方法(例如，VEP^[1]、ANNOVAR^[2])通常是以变异为单位，基于参考基因模型(reference gene model)独立处理每个变异产生的影响。显然，这种假定每个变异独立工作产生影响的做法忽略了变异所在的基因组环境，是不符合生物学实际情况的。大规模人群基因型数据的分析结果显示，有大量的变异影响被这种独立处理每个变异的注释方法错误处理。

另外，在生物学通路水平，目前的注释方法(例如，DAVID^[3])仍然是以富集分析为主，其利用统计显著性检验方法，根据用户提交的大量受变异影响的基因列表中找出在某些通路中显著富集的结果。然而，统计检验分析旨在找出那些被变异显著影响的生物学通路，并不能准确直接地指出变异对生物学通路的具体影响。

因此，需要有克服上述缺陷的准确确定遗传变异对生物分子功能的影响的方法。

技术实现要素：

本发明针对上述缺陷，提供了一种基于基因组环境确定遗传变异影响的方法，所述方法以每个基因为单位注释该基因上所有变异共同的影响，其包括：1)将所有变异根据其坐标位置映射到给定基因模型的各个基因上；2)根据各基因上的所有变异重构出各基因的个体化序列；3)对所得个体化序列进行分析以得到变异对该基因的影响。

本发明还提供了一种基于基因组环境确定遗传变异对生物学通路影响的方法，所述方法包括如下步骤：1)将基因/蛋白相互作用通路抽象成有向无环图；2)删除功能缺失基因对应的图节点以及相应的边；3)找出因节点删除造成的最远的不连通路径。

本发明的方法充分考虑了变异所在的基因组环境，避免了大量的注释错误，提高了注释变异影响的准确性。

附图说明

图1示出了使用本发明方法对蛋白编码基因进行注释的流程图。

图2示出了使用本发明方法确定变异对基因CHD7影响的结果。

图3示出了使用本发明方法重新分析1000基因组的基因组数据的结果。

图4示出了使用本发明方法对转录因子结合位点进行注释的流程图。

图5示出了使用本发明方法确定变异对转录因子结合位点TFAP2结合位点影响的结果。

图6示出了利用1000基因组和GTEx项目的基因组数据对TFBS注释冲突进行具体分析的结果。

图7示出了使用本发明方法对microRNA进行注释的流程图。

图8示出了使用本发明方法所证明的SNP rs56301829与SNP rs2276448变异不会导致microRNA MIMAT0027683失去对基因ZNF716的调控的示意图。

图9示出了利用1000基因组数据分析来自TargetScan和miRanda的转录因子结合点上的互补突变的结果。

图10示出了使用本发明的方法确定变异对生物学通路的影响的一般性流程图。

具体实施方式

如上所述，本发明提供了一种基于基因组环境确定遗传变异功能影响的方法，所述方法以每个基因为单位注释该基因上所有变异共同的影响，其包括：1)将所有变异根据其坐标位置映射到给定基因模型的各个基因上；2)根据各基因上的所有变异重构出各基因的个体化序列；3)对所得个体化序列进行分析以得到变异对该基因的影响。

在本发明的方法中，术语“给定基因模型”是指用户指定的一套完整的基于参考基因组的基因结构描述，包括基因中可变剪切转录本外显子、内含子在参考基因组上的坐标。

在本发明的方法中，术语“个体化序列”是指根据个人基因型得到的个体基因组中的真实序列。

在一个实施方案中，所述基因为蛋白编码基因。在一个具体的实施方案中，所述蛋白编码基因为CHD7。在蛋白编码基因的情况下，本发明方法的步骤2)通过根据各个基因上的所有变异推断该基因的蛋白编码区，并将其翻译成蛋白序列来进行，并且步骤3)通过将所得蛋白序列与已知的参考蛋白序列比较来进行。

在一个实施方案中，所述基因为转录因子结合位点。在一个具体的实施方案中，所述转录因子结合位点为转录因子TFAP2结合位点。

在一个实施方案中，所述基因为microRNA。具体地，在确定变异对microRNA靶基因影响的情况下，本发明方法包括确定变异对microRNA生成和microRNA靶位点的影响。

本发明还提供了一种基于基因组环境确定遗传变异对生物学通路影响的方法，所述方法包括如下步骤：1)将基因/蛋白相互作用通路抽象成有向无环图；2)删除功能缺失基因对应的图节点以及相应的边；3)找出因节点删除造成的最远的不连通路径。

以下通过具体实施例来说明本发明的内容。应理解，所述具体实施例仅为说明目的，并不意味着本发明的内容仅限于具体实施例。

实施例1：对蛋白编码基因的注释

在本实施例中，使用本发明的方法对蛋白编码基因进行注释。

对于蛋白编码基因的注释，可通过以下步骤进行：1)将变异根据其坐标位置映射到给定基因模型的各个基因上；2)根据各个基因上的所有变异推断该基因的蛋白编码区；具体地，对于含有影响剪切的变异的转录本，在给定区间内(默认为+/-100bp)的范围内寻找隐藏的剪切位点，以及由其他变异造成的新的剪切位点。对于无法找到可替代的剪切位点的转录本则外显子跳过和内含子保留都会被考虑。3)根据所得的蛋白编码区序列将其翻译成蛋白序列；4)将所得蛋白序列与已知的参考蛋白序列比较以确定变异对该基因的影响。图1示出了使用本发明方法对蛋白编码基因进行注释的流程图。

以千人基因组项目中HG02861个体基因组中的变异为例进行具体说明。对于SNP rs549508773，VEP会将它注释为发生在基因CHD7上的终止密码子(gag→tag)获得突变(STOP-gained variant)。然而，本发明(COPE，基于基因组环境的变异注释工具)会根据HG02861个体所包含的所有变异，发现该SNP的旁边还存在另一个SNP rs567756521，两个SNP的综合结果是导致基因CHD7中的一个氨基酸替换(密码子gag→ttg)(如图2)。因而，我们的结果发现在考虑变异基因组环境后，相对于现有的VEP来说，可以更准确地注释变异影响。

为了进一步说明本发明方法的优势，我们使用本发明方法重新分析了千人基因组的基因组数据(下载于ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/)。我们选择了其中被VEP注释为功能缺失型的9728个终止密码子获得突变，2092个移码突变和5559个剪切突变重新用本发明方法进行分析，结果如图3所示。从图3中可以看出，有大约23.2％影响剪切的变异，6.5％移码突变以及2.1％终止密码子获得的变异是被以前算法错误注释的(图3a)。富集分析发现，这些被错误注释为功能缺失的基因与特定的生物功能相关(图3b)。另外，本发明还发现了38对由两个突变造成新的终止密码子获得的变异(被VEP忽略了的功能缺失变异)(图3c)。由此进一步说明了本发明方法的优势。

实施例2：对转录因子结合位点的注释

在本实施例中，使用本发明的方法对转录因子结合位点进行注释。

对转录因子结合位点的注释分为转录因子结合位点的丢失(TFBS loss)和转录因子结合位点获得(TFBS gain)。对于转录因子结合位点丢失的预测，通过对比变异前后序列对应的位点权重矩阵得分的大小判断是否存在转录因子结合位点丢失。另一方面，对于转录因子结合位点获得的预测，首先找出所有位于启动子内的变异，进而重构启动子序列。然后，根据位点权重矩阵预测转录因子结合位点获得。图4示出了使用本发明方法对转录因子结合位点进行注释的流程图。

以SNP rs9274541和rs9274542为例进行具体说明，这两个SNP都是位于转录因子TFAP2结合位点上的两个SNP。以往的注释方法，如SNP2TFBS[⁴]，会分别处理这两个变异，得到结果是rs9274541能降低TFAP2的结合强度(ΔPMW分数<0)，而rs9274542能增加它的结合强度(ΔPMW分数>0)(如图5a)，两者冲突。而本发明的方法综合考虑两者的影响，根据变异重构出个体化的转录因子结合位点，然后根据位点权重矩阵打分，比较参考基因组得分和个体化序列得分，最后判断为降低了TFAP2的结合强度(ΔPMW分数<0)(如图5b)。

接着，本发明人利用1000基因组和GTEx项目的基因组数据具体分析了上述注释出现冲突的情况，结果如图6所述，从图中可以看出大约有28％的含有多个变异的TFBSs会出现这类注释冲突，进一步说明了本发明方法的优势。

实施例3：对microRNA的注释

在本实施例中，使用本发明的方法对microRNA进行注释。

对miRNA的注释，同样根据给定的变异重构出个体化序列综合判断变异对miRNA生成和miRNA靶位点的影响。对于microRNA生成的注释，旨在预测位于pre-microRNA上的基因组变异对pre-microRNA二级结构最小自由能的影响。这里用于计算pre-microRNA最小自由能的工具是RNAfold。首先，找出所有位于同一个pre-microRNA上的变异。然后，重构出真实的pre-microRNA序列。最后，计算基因组变异发生前后，该pre-microRNA的最小自由能的变化，并以此作为对microRNA生成的影响。microRNA靶标结合注释是指预测变异对microRNA与3’UTR的相互作用的影响。为了避免预测假阳性过高，使用了两种常用的工具是TargetScan^[5]和miRanda^[6]。具体地，将靶标结合缺失定义为在参考序列下TargetScan和miRanda都可以预测到的靶标基因；而在变异发生后，TargetScan和miRanda都不能预测到的靶标基因。靶标结合获得则是指在参考序列下TargetScan和miRanda都不能预测到的靶标基因；而在变异发生后，TargetScan和miRanda都能预测到的靶标基因。具体地，首先同样根据变异重构出所有3’UTR和所有microRNA序列。然后，利用miRanda和TargetScan预测新的结合位点。最后，根据上面的定义找出所有的靶标丢失和靶标获得位点。图7示出了使用本发明方法对microRNA进行注释的流程图。

以SNP rs2276448和rs56301829为例进行具体说明，其分别是在microRNA MIMAT0027683和基因ZNF716上的两个变异。microRNA对靶位点的识别主要依靠microRNA种子区(seed region)内的几个碱基与靶位点碱基的互补配对。SNP rs2276448是一个发生在microRNA种子区的突变(可以影响microRNA靶位点)，在不考虑其他因素的情况下，基因ZNF716会被错误地认为不再受到该microRNA的调控。本发明的方法在预测microRNA上的变异对其靶位点的影响时会考虑靶位点上的序列环境，发现rs56301829恰好是一个发生在基因ZNF716上且与SNP rs2276448互补的突变。因此，这对互补的变异不会导致microRNA MIMAT0027683失去对基因ZNF716的调控(如图8)。利用1000基因组数据分析了来自TargetScan^[5]和miRanda^[6]的转录因子结合点上的互补突变，结果如图9所示，从图中可以看出超过半数的microRNA的靶位点会存在这种互补的变异，进一步说明了本发明方法的优势。

实施例4：确定变异对生物学通路的影响

在本实施例中，使用本发明的方法确定变异对生物学通路的影响，具体步骤如下：1)将基因/蛋白相互作用通路抽象成有向无环图；2)删除功能缺失基因对应的图节点以及相应的边；3)找出因节点删除造成的最远的不连通路径。图10示出了使用本发明的方法确定变异对生物学通路的影响的一般性流程图。

以基因GNAS发生功能缺失为例，本发明首先找到含有此基因的生物学通路，然后转化为有向无环图。最后我们根据删除基因GNAS对应的节点后，找出因该节点删除造成的最远不连通路径(如对胰腺分泌(Pancreatic secretion)通路造成的最远不连通路径是从SCTR(Secretin Receptor，分泌素受体)到ADCY1(Adenylate Cyclase 1，腺苷酸环化酶1))。本发明会直接指出变异影响的基因在发生功能缺失后对生物学通路造成的直接影响，相比于传统的富集分析来说，这样的注释更直观准确，有助于遗传诊断分析。

参考文献：

[1]Deriving the consequences of genomic variants with the Ensembl API and SNP Effect Predictor(McLaren et al.,2010).

[2]ANNOVAR:functional annotation of genetic variants from high-throughput sequencing data(Wang et al.,2010).

[3]DAVID Bioinformatics Resources:expanded annotation database and novel algorithms to better extract biology from large gene lists(Huang,D.W.et al.,2007).

[4]SNP2TFBS-a database of regulatory SNPs affecting predicted transcription factor binding site affinity(Kumar,S.et al.,2017).

[5]Human MicroRNA targets(John,B.et al.,2004).

[6]Conserved seed pairing,often flanked by adenosines,indicates that thousands of human genes are microRNA targets(Lewis,B.P.et al.,2004).