CN104598775A - 一种rna编辑事件识别机制及其高效实现方案 - Google Patents

Abstract

本发明涉及分子遗传学和信息技术领域，使用计算机信息技术解决分子遗传学中的RNA编辑事件的识别问题。本发明提出了一种系统整合的RNA编辑事件识别机制，以及其高效率实现的解决方案。本发明提出的机制包括规则型滤除算法和统计型滤除算法两大模块，其中规则型模块考虑测序的实际特点和可能变异位点在基因中的分布，并对可能变异位点设定规则进行约束；统计型模块则在规则型模块处理后，对符合规则约束的位点进行统计学计算，并根据结果判断是否落入拒绝域，确定可能变异位点是否被滤除。为了实现上述机制，本发明还提出了一种高效率的识别机制实现方案。即使用高级编程语言、数据库和统计学工具相结合的方式实现相应算法模块。

Description

一种RNA编辑事件识别机制及其高效实现方案

技术领域

本发明涉及分子遗传学和信息技术领域，使用计算机信息技术解决分子遗传学中的RNA编辑事件的识别问题。

相关背景

长期以来，在分子遗传学领域，从庞杂的生物遗传信息中获取目标结果是一项艰巨的挑战。由于基因数据测序昂贵，可获取的数据样本匮乏，很难发现一些发生概率较小但是具有显著意义的生物学事件。随着低成本、大数据量测序，即第二代高通量测序的发展，分子遗传学研究工作者能够以更低的成本获取更多的测序源数据。然而，高通量测序获取的海量原始数据已经远远超过了人力实验研究的能力范围，这使得相关研究人员必须借助计算机信息技术来对源数据进行分析处理，这项工作已经和通过传统的实验得到实验室、临床实践的观测结果同样重要。

高通量测序的普及和信息技术的应用，使得许多分支研究工作能够开展，如本发明聚焦的RNA编辑事件研究。RNA编辑事件代表了一种全新的生物性状调控表达方式，被公认为当今分子遗传学界的热点问题。

RNA编辑事件的产生在转录阶段，通过碱基置换改变前体RNA，从而产生与原始DNA模板携带信息存在差异的RNA，即DNA层面上没有发生变异，而RNA上产生了变异。RNA编辑可发生于多种生物体内，从低等的病毒、细菌、真菌，到高等植物、动物，包括范围较广。以高等动物最常见的A-I类型的RNA编辑事件为例，RNA腺苷脱氨酶(Adenosine Deaminase Acting on RNA，ADAR)起催化作用，特异性识别DNA双链底物上的腺嘌呤(A)，经脱氨酶的催化作用，A被转换为I(次黄嘌呤)，次黄嘌呤在翻译和剪接过程被识别为G，从而可能导致性状表达的偏差。

RNA编辑事件可能有多种生理作用，包括RNA编辑事件能改变氨基酸序列、翻译起始子、终止密码子，或者破坏或新建剪接信号、影响miRNA前体的加工及成熟体的靶向功能。异常的RNA编辑事件与多种疾病密切相关。

RNA编辑事件的重要性，使得从高通量测序数据中寻找RNA编辑位点十分重要。相关研究者提出了RNA编辑事件部分识别的规律和方法，并研发出了一些基本处理工具。但是，由于研究尚处于起步阶段，现有的识别规律和方法存在功能不完整，对编辑事件特点考虑不周全等问题。由Picardi E,Pesole G等人实现的信息处理工具REDItools则存在处理中间过程不清晰，执行效率较低等问题。相关领域研究人员迫切需要一个系统整合的识别机制，并具有高效实现的解决方案。

发明内容

本发明提出了一种系统整合RNA编辑事件的识别机制，及其高效率的实现方案。值得指出的是，虽然本发明应用领域为分子遗传学，但是识别机制和实现方案均不依赖于特定的核苷酸序列，即适用于任何符合本发明要求格式的测序数据。

本发明提出的RNA编辑事件识别机制流程图见说明书附图1。该机制可分为规则型和统计型滤除算法两大模块，其中规则型模块考虑测序的实际特点和可能变异位点在基因中的分布，并对可能变异位点设定规则进行约束；统计型模块则在规则型模块处理后，对符合规则约束的位点进行统计学计算，并根据结果判断是否落入拒绝域，确定候选位点是否被滤除。

根据研究数据的特点，本发明机制可支持两种模式，即De novo模式(源数据仅包括RNA测序数据)和DNA-RNA模式(源数据包括RNA和DNA测序数据)。

规则型和统计型两大模块均可细分为公共和专有两种子模块，其中公共子模块适用于两种模式，而专有子模块仅适用于DNA-RNA模式，是公共子模块的补充。每次对不同源数据进行RNA编辑事件的识别，产生的结果编辑位点分布和数量均可能不同。

本发明提出的机制中，提及的待分析源数据格式，一律为VCF(variant callingformat)格式。

本发明提出的机制基本原理如下：

1规则型滤除算法公共子模块

本发明提出的机制中，规则型滤除算法的公共子模块是包含内容最多的部分。该子模块算法流程图见说明书附图2，其原理如下：

1.1指定型滤除算法

RNA编辑事件包括12种类型。根据研究内容的不同，研究工作者往往不需 要同时分析全部类型，而是仅选取最重要、最相关的类型，比如高等动物常见的A-I(等效A-G)类型编辑。

根据VCF格式数据特点，源数据文件为方便对比亲族等分类遗传信息，往往将同一病人的不同组织合并为同一个文件，每个位点的突变基因型是所有组织的突变类型并集。如果某位点存在突变型，只能确定所有亲族的全体中至少一人在该位点发生了突变。如果要识别具体是否发生突变，需要读取位点基因型，如果为纯合(如AA)则没有发生RNA编辑事件。

1.2质量控制滤除算法

质量控制，主要针对测序数据的质量、测序的可靠水平进行约束控制。有质量保证的高通量测序，是后续数据处理、分析以及实验验证的基础。低质量的测序数据将会产生大量假阳性数据，干扰结果的判断；而同时也会导致假阴性数据的产生，遗漏有价值的信息。

虽然高质量的测序具有显著的意义，但是受测序技术的制约，测序的每个位点上的碱基都可能出现错误。根据专业测序统计规律，对不同测序位点、结合当时的测序条件，给出了每个位点的测序精度参考值，即碱基质量因数。碱基质量因数(Q)同测序点错误概率(P)的计算关系见公式(1)：

Q＝-lgp    (1)

除了测序质量外，还有一个关键的参数是测序覆盖深度。在测序过程中，每个点并不是均等被测到，被测到次数取决于增殖转录本数量。每个位点被测的次数越多，可提供的信息量越大，所产生的结果就越有说服力。同时，有足够的覆盖深度，能同时提供充足的参考位点和变异位点，使得RNA编辑水平的结果更为科学可信。

1.3重复区域滤除算法

人类的基因组分布结构，可以认为相当于基因的一组DNA序列和多组反复重现的DNA序列组成。多组重复重现的序列被称为重复序列，组成的区域即重复区域。这些序列组织形式有两种：串联重复和分散重复。前一种，成簇存在于DNA的特定区域，后一种分散于DNA的各个位点上。

这些重复序列在DNA上大量存在的，但是并没有转录活性，即无法通过转录生成转录本，进而指导蛋白质的合成。因此由于现阶段研究认为此种重复区域不影响性状的表达，故存在于该区域的编辑位点对后续的研究工作意义不大，所以予以滤除。

但是，有一种被命名为Alu的重复序列，被认为具有一定的生理功能。Alu重复序列，是哺乳动物基因组中SINE族的一部分。平均每4000到6000个碱基中，就有一组Alu重复序列。其表现特征为，序列组成相对固，典型的人鳄梨基 因组上该序列长282个碱基，由两个同源但有差别的亚基构成。

Alu序列被认为可能存在多种功能。(1)部分序列同某些基因内含子剪切接头较为相似，因此，可能参与了信使RNA的加工与成熟。(2)Alu序列在人类基因中，极为不寻常地大量存在，可能同遗传重组及染色体不稳定等现象有关。(3)Alu序列可能具有调节作用。由于Alu序列存在如此巨大的遗传学意义可能性，所以该序列予以保留。

1.4测序连接区域滤除算法

第二代高通量测序，产生大量数据方便识别RNA序列的变异。但是在测序过程中，会遇到转录本的拼接问题，即将大量的零散转录本拼接测序并映射到原序列上。正确的映射非常关键，这是避免当RNA映射到参考基因组的时候，出现不匹配而被错误地解读为解释为单核苷酸多态性。影响正确映射的原因主要有：(1)基因组中高度相似区域进行的干扰，迷惑正确的映射匹配；(2)多数据流，通过特定算法映射时，在拼接处导致拼接连接点的失真，这也是最大的障碍和挑战。

因此，当读取到跨越拼接点的转录本时，可信区段应当比测序样本略短，避免因映射时的拼接错误导致参考基因组的结果出现错误。

1.5核苷酸多态性滤除算法

从定义角度考虑，RNA编辑事件属于DNA层面上没有变异，而RNA上发生了突变的类型。如果之前的阶段，识别出的位点是人类基因组数据库中的DNA已知变异位点，那么就排除该位点是RNA编辑位点的可能性。

2规则型滤除算法专有子模块

当输入数据同时包括病人的DNA和RNA测序数据时，该子模块生效。

以人类最典型的A-I(等效A-G)类型RNA编辑事件为例。转录后，只有当A碱基受ADAR脱氨酶作用变为肌苷I，并在修饰加工过程中以G的形式参与RNA合成时，才能被称为RNA编辑事件。如果病人本身，在转录前就已经发生了变异，就不能被认为是RNA编辑事件。由于DNA的双链特性，如果是杂合的A型位点，则仍有另一条链是等位基因，从而破坏了RNA编辑事件成立的条件。所以，某位点的DNA测序为纯合，是该位点上发生RNA编辑事件的必要条件，不符合该条件的位点将被滤除。

3统计型滤除算法专有子模块

当输入数据同时包括病人的DNA和RNA测序数据时，该子模块生效。

仍以人类最典型的A-I(等效A-G)类型RNA编辑事件为例。经过规则型滤除算法的筛选，保留的RNA编辑事件候选位点在每次测序中只可能有两种结果，即A和G，并且两种结果互斥。由于每个位点在高通量测序中可能属于多个转录本，故在转录阶段发生的RNA编辑事件不一定影响所有的转录本。所以，每个候选编辑位点上同时具有参考位点(A)测序数和变异位点(G)测序数两种数据。本子模块将根据二项分布等原理，通过计算假定该位点确实发生编辑事件(对应的DNA测序没有发生错误)的概率p₁，以及假定该位点事实上没有发生编辑事件(对应DNA测序发生错误)的概率p₂，并设定阈值t，使得p₁/p₂>t，来判断是否在统计学意义上发生了编辑事件

4统计型滤除算法公共子模块

通过之前阶段的滤除算法，已经得到了理论分析的位点结果。这些结果将作为对候选位点结果的“假设”，即认为，已识别出的位点是真实发生了编辑事件的位点。为了作为对比，该子模块将参考人类已知的RNA编辑位点，在统计意义上确定之前结果的“假设”是偶然得到的，还是与真实发生的情况不存在统计意义上的差异。

在本模块中，首先为降低结果的假阳性，提高结果的准确度，将设定拒绝假设的阈值，采用显著性水平分析中的Fisher精确检验方法，计算一个精确的结果，即p值，来作为判断是否偏离原假设的根据。

随后，为了均衡敏感性和特异性，本算法采用了FDR(False Discovery Rate)方法，对已获得的显著性水平进行进一步的控制。该方法同样设置拒绝域，落在拒绝域外的位点将被保留，成为真正识别的RNA编辑位点。

为了实现上述机制，本发明还提出了一种高效率的实现方案。即使用高级编程语言对原始数据进行结构化处理，并将格式化数据导入数据库，使用数据库的查询，更新，插入，联表查询等功能实现上述功能，在涉及统计学计算的模块中，使用自主编写代码和统计学工具相结合的方式完成功能。为简便描述，以人类最典型的A-I(等效A-G)类型RNA编辑事件为例介绍实现原理：

1.信息导入： 

将VCF格式的源数据，通过高级编程语言处理，去除注释信息，将主体测序信息进行格式化处理，导入数据库。

规则型滤除算法公共子模块

2.指定型滤除算法：

导入数据库的信息中，对每一行的位点信息进行筛选，要求参考碱 基列为A碱基，变异碱基列为G(等效I)的位点，同时该位点的参考基因型不为纯合，符合此标准的位点予以保留，进入步骤3，其余位点被滤除。

3.质量控制滤除算法：

首先设定碱基测序质量和测序深度的最低标准值，大于此标准的位点予以保留，进入步骤4，其余位点被滤除。

4.重复区域滤除算法：

从官方网站下载公共参考的重复区域文件，读取记录了全部重复区域的位点坐标区段，载入数据库；

以步骤3的输出作为本阶段算法输入，如果候选编辑位点不在重复区域中，则进入步骤5；

如果候选编辑位点在重复区域中，但属于SINE/Alu功能区，则进入步骤5；

如果候选编辑位点在重复区域中，但不属于SINE/Alu功能区，则予以滤除。

5.测序连接区域滤除算法

从官方网站下载公共参考的基因区域文件，读取记录了基因区域标识的区段，载入数据库；

当RNA编辑候选位点落在编码区边缘，即映射拼接边缘时，被认为测序不可信，予以滤除，在边缘外位点进入步骤6。边缘范围设定临界值n，即编码区临近点以外、以内各n/2个碱基。

6.核苷酸多态性滤除算法

从官方网站下载公共参考的单核苷酸多态性文件，该文件据库包含了最新的人类基因组变异的数据。读取记录发生单核苷酸多态性位点，载入数据库；

当RNA编辑候选位点，其坐标在DNA上对应的点，属于核苷酸多态性，即DNA变异时，该点予以滤除，否则进入步骤7。

规则型滤除算法专有子模块

7.将分析位点同DNA数据进行比对，当RNA编辑事件候选位点对应的DNA位点属于A碱基纯合时，保留该位点进入步骤8，否则滤除。 

统计型滤除算法专有子模块

8.使用高级编程语言编写相关算法

a)假设该位点确实发生了编辑事件，则根据二项分布，该位点被测为A碱基和测为G碱基是互斥离散型随机事件，设ref(A)代表该位点A碱 基被测到次数，ref(G)代表G碱基被测到的次数，则该位点发生编辑事件概率f(T)可由公式2计算得到。

f(A)＝ref(A)/(ref(A)+ref(G))

f(G)＝ref(G)/(ref(A)+ref(G))    (2)

f(T)＝f(A)^ref(A)f(G)^ref(G)

由公式2可知，当f(A)和f(G)相差较多时，f(T)值较小，该位点发生编辑事件概较低，反映到测序信息上，说明该位点编辑水平偏低(f(A)较大)或者偏高(f(G)较大)，无论哪种情况，都较为偏离平均水平，所以f(T)大小变化，符合了这一判断。

b)假设该位点没有发生编辑事件。如果该位点没有发生编辑事件，却被之前的滤除规则识别为RNA编辑事件发生，说明DNA序列测序错误，即DNA不是纯合的AA。由公式3所示，测序发生错误的概率可以通过测序品质因数计算获得：

f(F)＝10^-Q/10    (3)

由公式3可知，测序品质因数(Q)越高，代表该位点越可信，从而DNA的验证越有效。而f(F)同Q值负相关，Q越大，代表没有发生编辑的越低，符合理论上的推断。

c)本发明认为，若假设发生了编辑事件的概率，是假设没有发生编辑事件概率的一万倍时，则认为确实发生了RNA编辑事件。为方便计算，取对数，可得公式4：

$\begin{array}{c}f\left(T\right)/f\left(F\right)>=10000\⇒\mathrm{lgf}\left(T\right)-\mathrm{lgf}\left(F\right)>=4\\ \mathrm{lgf}\left(T\right)+Q/10>=4\end{array}---\left(4\right)$

统计型滤除算法公共子模块

9.Fisher精确检验

本算法的执行，需要一个2*2形式的列表，其基本形式如表1所示：

表1Fisher精确检验的求解列表

	Num(A)	Num(G)	Row Total
				Expected	a	b	a+b
Observed	c	d	c+d
				Column Total	a+c	b+d	a+b+c+d＝n

在表1中，需要提供两组数据以供检验运算，即Expected(期望)组和Observed(观测)组。Expected组参考人类已经发现的RNA编辑位点数据库darned， 并结合之前滤除阶段获得的结果进行计算，其计算规则为：

a)逐个遍历上步骤保留的RNA编辑事件候选位点，设定期望组计数器cal1和cal2，cal1负责累加未变异位点数量，cal2负责累加已变异位点数量。如果该位点为已知编辑位点，则执行步骤b，否则执行步骤c；

b)如果编辑位点已知，分析该位点未变异(A)和已变异(G)的数量，并使用对应的计数器计数；

c)如果变异位点未知，则对于期望组来说，属于完全变异类型，cal1不变化，cal2累加当前位点覆盖深度。

d)遍历全部候选位点，使用公式5计算计数器的平均值，得到期望值a,b(num代表候选编辑位点数量)。

a＝cal1/num

               (5) 

b＝cal2/num

在表1中所示的观测组数据，即为每个点的实际观测结果，其中c值代表碱基A的测序深度，d代表碱基G的测序深度。

获得了Fisher精确检验必须的2*2联表数据后，显著性水平(P值)计算公式如公式6：

计算获得的p值后，进入步骤10。

10.使用FDR控制显著性水平

a)获取Fisher精确检验的p值数组。将所有位点的p值从小到大进行排序，获得p(1),p(2),...,p(m)p(1),p(2),...,p(m)；

b)设定FDR的阈值(q值)。则对经过排序的数列遍历，寻找最大的正整数i，使得p(i)<＝(i*q)/m成立；

c)输出数列p(1),p(2),…,p(i)对应的位点，作为真正具有统计意义精确的RNA编辑位点。

经过上述步骤，获得的位点即为RNA编辑事件机制识别出的位点。

附图说明

图1.RNA编辑事件识别机制流程图

图2.规则型滤除算法公共子模块算法流程图

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本实施例使用的计算机信息工具包括：(1)负责RNA编辑位点识别过程中，测序数据的输入、处理、输出的MySQL数据库管理系统，(2)负责设计数据接口，使得数据具有私密性，SQL语句可复用，具有良好扩展性的Java程序设计语言，(3)具有强大的统计学计算能力的R语言。

本实施例使用的数据同时包括了RNA测序和成本高昂的DNA测序。

本实施例设定关键参数后，使用Java编写文件处理程序，将源数据导入MySQL数据库。使用MySQL数据库的查询，更新，插入，联表查询等功能实现规则型滤除算法模块，使用Java编写统计型滤除算法专有子模块，并使用Java调用R语言实现统计型滤除算法公共子模块。

本实施例为清晰明了，仍处理A-I(G)的编辑类型。具体实施方案如下：

1.信息导入： 

将VCF格式的源数据，通过Java处理，去除注释信息，将主体测序信息进行格式化处理，并仅将Info一栏为”PASS”的位点导入数据库。

规则型滤除算法公共子模块

2.指定型滤除算法：

条件查询参考碱基列为A碱基，变异碱基列为G(等效I)，且参考基因型杂合的位点，查询到的数据将被插入到新表，进入步骤3。

3.质量控制滤除算法：

首先设定测序质量最低标准20(测序错误概率1％)，测序深度的最低为6，条件查询大于此标准的位点并插入新表，进入步骤4。

4.重复区域滤除算法：

从官方网站下载公共参考的重复区域文件，读取记录了全部重复区域的位点坐标区段，载入数据库；

以步骤3创建的新表与本阶段公共参考的重复区域数据表进行联表查询，如果可能变异位点不在重复区域表中，则插入新表，进入步骤5；

如果可能变异位点在重复区域表中，但属于SINE/Alu功能区，则插 入新表，进入步骤5；

如果可能变异位点在重复区域表中，但不属于SINE/Alu功能区，则不再保留该位点信息。

5.测序连接区域滤除算法

从官方网站下载公共参考的基因区域文件，读取记录了基因区域标识的区段，载入数据库；

设定边缘区域的临界值为4(至少覆盖一组密码子的最小偶数)。

以步骤4创建的新表与本阶段公共参考的基因重复区域数据表进行联表查询，如果可能变异位点在编码区边缘2个碱基以上，则插入新表，进入步骤6；

6.核苷酸多态性滤除算法

从官方网站下载公共参考的单核苷酸多态性文件，该文件据库包含了最新的人类基因组变异的数据。读取记录发生单核苷酸多态性位点，载入数据库；

以步骤5创建的新表与本阶段公共参考的核苷酸多态性数据表进行联表查询，如果可能变异位点不在核苷酸多态性表中，则插入新表，进入步骤7；

规则型滤除算法专有子模块

7.以步骤6创建的新表与DNA测序数据表进行联表查询，如果可能变异位点在DNA测序数据表中属于AA纯合型，则插入新表，进入步骤8；

统计型滤除算法专有子模块

8.使用高级编程语言编写相关算法

本步骤使用Java编写相关算法，实现发明内容的对应部分，并将合格的计算结果插入新表，进入步骤9。

统计型滤除算法公共子模块

9.调用R语言的Fisher精确检验模块，计算出每个可能变异位点的p值，记录所有p值作为一个数列，进入步骤10.

10.调用R语言的p.adjust模块，根据p值数列计算FDR值，将p值和FDR值的拒绝域均设定为5％，拒绝域外的位点将插入结果表。

通过实施例检验，本信息处理实现方案可以得到本发明提及机制的结果位点，可以为后续相关研究者的工作做好坚实的铺垫。

Claims (3)

1.一种RNA编辑事件识别机制，其特征在于：要求源数据使用VCF(variantcalling format)格式，识别机制从测序位点的自然属性和统计学计算后的数学属性两方面出发，识别出潜在的RNA编辑位点。

2.如权利要求书1所述的识别机制，识别机制包括规则型和统计型两个滤除算法模块，支持De novo模式(源数据仅包括RNA测序数据)和DNA-RNA模式(源数据包括RNA和DNA测序数据)两种模式的处理，其特征在于：

(1)规则型滤除算法模块考虑测序的实际特点和可能变异位点在基因中的分布，规则对可能变异位点的约束内容包括：可灵活指定检测编辑类型；测序数据质量和测序深度的最低标准要求；是否属于重复区域；是否属于测序连接区域边缘；是否属于单核苷酸多态性；是否在DNA层面上已经发生了变异；

(2)统计型滤除算法模块考虑已获取位点的统计学特点，从最大似然概率，Fisher精确检验和FDR(false discovery rate)检验等方面对位点进行计算，并根据结果判断是否落入拒绝域，确定可能变异位点是否被滤除；

(3)经过上述步骤识别出的位点，将作为本识别机制的结果编辑位点。

3.基于权利要求书1的机制，设计了一种高效率的识别机制实现方案，其特征在于：

(1)使用高级编程语言对源数据进行结构化处理，并将格式化数据导入数据库；

(2)使用数据库的查询，更新，插入，联表查询等功能实现规则型滤除算法模块；

(3)使用高级编程语言和统计学工具相结合的方式编写代码实现概率型滤除算法模块。