CN116312776B

CN116312776B - 一种检测差异化rna编辑位点的方法

Info

Publication number: CN116312776B
Application number: CN202211574812.8A
Authority: CN
Inventors: 张弛; 邱建华; 张坤明; 梁红远; 瞿爱东
Original assignee: SHANGHAI INSTITUTE OF BIOLOGICAL PRODUCTS CO LTD
Current assignee: SHANGHAI INSTITUTE OF BIOLOGICAL PRODUCTS CO LTD
Priority date: 2022-12-08
Filing date: 2022-12-08
Publication date: 2024-01-19
Anticipated expiration: 2042-12-08
Also published as: CN116312776A; WO2024120496A1

Abstract

本发明公开了一种检测差异化RNA编辑位点的方法。所述方法包括以下步骤：对一组包含两种及以上处理方式，或肿瘤与正常组织对照样本集进行高通量测序，分别获得DNA与RNA测序数据；通过对比DNA与RNA测序数据进行DNA/RNA联合突变检测，获得RNA编辑位点；对RNA候选编辑位点在在各样本集中的突变深度进行测量并组成数据集，并使用统计学手段评估每一个候选位点相对不同样本集的差异化编辑程度，从而筛选出具有统计学意义的差异化RNA编辑位点。使用本方法能以高灵敏度检出高置信度的差异化RNA编辑位点。

Description

一种检测差异化RNA编辑位点的方法

技术领域

本发明属于生物技术领域，具体地，本发明涉及一种差异化RNA编辑位点的检测方法。

背景技术

核糖核酸(RNA)编辑(editing)属于RNA的转录后修饰的一类，是指在RNA水平上所发生的一个或多个碱基序列的变化。通常来讲，生物界存在两类RNA编辑，分别是碱基替换编辑以及碱基插入或丢失编辑。这两种编辑方式的任何一种都能产生异于基因模板的转录本，从而能影响RNA序列与结构稳定性、可变剪切、以及编码蛋白质序列，从而对基因表达调控和基因产物的多样化发挥重要作用。目前研究最为深入的单碱基替换编辑主要由核苷酸脱氨基反应介导。在哺乳动物中，最常见的是腺苷脱氨酶ADAR家族催化的腺苷(adenosine)向肌苷(inosine)编辑；以及胞苷脱氨酶APOBEC家族催化的胞苷(cytidine)向尿苷(uridine)编辑。近期，由于高通量测序技术的发展和普及，为研究者在全转录组层面上研究RNA编辑位点提供了可能。

RNA编辑位点检测的一般性方法有两种。一是通过比对RNA高通量测序结果和同物种参考基因组，并通过过滤已知基因组层面单核苷酸多样性位点(SNP)从而获得RNA测序中异于已知序列的突变或编辑位点；其代表方法是SNPiR方法[4]，但由于检出噪声过大现已较少使用。另外一种方法是通过比对RNA高通量测序结果和同样本DNA高通量测序结果，并通过联合突变检测发现仅在RNA层面上发生的异于基因组序列的编辑位点，其代表方法是VaDiR方法(Neums等人)和专利CN105483210。这两种方法均能较为有效的检测RNA编辑位点。但是这些方法存在两个主要的不足之处：首先，这些方法不涉及分析具体实验条件或样本组别对RNA编辑位点的影响，也不能确定RNA编辑位点是否受到这些因素影响，所以检出的编辑位点不是差异化RNA编辑位点；其次，这些方法检测RNA编辑位点的灵敏度有待改进，由于DNA测序仍然存在测序深度不足问题，若仅仅通过对比DNA和RNA测序结果，对异质化程度较高的样本不能完全排除DNA模板水平上的突变导致的假阳性问题。

差异化RNA编辑位点(Differential Variants on RNA，DVR)是指通过对比不同实验条件或样本组别，发现的存在显著编辑程度差异的RNA单核苷酸位点。检测差异化RNA编辑位点能考察不同条件下或样本组别对RNA编辑产生的影响，对深入研究生物学RNA单个核苷酸编辑的功能有重要意义。目前代表性的检测差异化RNA编辑位点是Wang.J等人设计的rMATS-DVR方法，该方法通过对比两个RNA测序样本集，并采用统计学方法(如对数秩检验法等)检测出潜在的差异化RNA编辑位点。然而，该方法存在的主要缺陷是不涉及DNA/RNA联合突变检测，所以未能考虑到较低程度的基因组DNA水平上的突变导致的RNA转录本单核苷酸多样性。因此，在异质化程度较高的样本，如肿瘤细胞样本，使用该方法存在大量的假阳性结果问题。

本发明涉及一种准确的结合DNA与RNA联合突变检测手段(图1)，并通过统计学方法准确检测差异化RNA编辑位点的方法。该方法能有效解决上述现有方法的不足。其主要创新之处在于首次使用差异化RNA编辑位点分析链特异的DNA/RNA联合检测结果，该方法在显著提升差异化RNA编辑位点灵敏度的同时，大幅度提高了检出位点是RNA编辑位点的置信度。对同一数据集使用本方法对比过往方法，能明显观察到上述的优势。

发明内容

本发明目的在于提供一种检测差异化RNA编辑位点的方法。

在本发明的第一方面，提供了一种检测差异化RNA编辑位点的方法，包括步骤：

(a)提供待检测差异化RNA编辑位点的N个样本的独立样本集；其中，N为≥2的正整数；

其中，对于每个样本而言，其独立样本集包括：经联合突变检测的RNA突变位点数据集Z；其中，所述的经联合突变检测的RNA突变位点数据集Z是将(i)单个样本的RNA测序数据与参考基因组进行序列比对后获得的第一RNA序列比对文件数据，与(ii)所述单个样本的DNA测序数据与参考基因组进行序列比对后获得的第二DNA序列比对文件数据，进行联合突变检测处理后所获得的RNA突变位点数据集；

(b)对于待检测差异化RNA编辑位点的两个或多个样本，将相应样本的各自的所述经联合突变检测的RNA突变位点数据集Z中的RNA突变位点进行合并，从而获得合并的候选RNA编辑位点数据集，记为第三数据集；

(c)将所述经联合突变检测的RNA突变位点数据集Z，进行基于RNA突变位点的测序深度的质量控制处理，从而获得符合测序深度质量标准的编辑位点的数据集，记为第四数据集；

(d)将步骤(b)获得的第三数据集中的合并的候选RNA编辑位点数据，与步骤(c)获得的第四数据集中的符合测序深度质量标准的编辑位点的数据进行合并，从而获得含有候选RNA编辑位点信息数据和编辑位点的测序深度数据的第五数据集；和

(e)对所述第五数据集中的所述候选RNA编辑位点信息与编辑位点的测序深度数据进行统计分析，从而确定差异化RNA编辑位点；

其中，步骤(b)和(c)的步骤可互换、先后、或同时进行。

在另一优选例中，所述的方法是非诊断和非治疗性的。

在另一优选例中，在步骤(b)中，对相应样本的各自的经联合突变检测的RNA突变位点数据集Z中的RNA突变位点进行合并后，还进行过滤处理，从而获得经过滤的合并的候选RNA编辑位点数据。

在另一优选例中，步骤(b)所述的过滤的操作选自下组：

(b1)去除均聚核苷酸序列中的突变；

(b2)去除位于基因组重复序列中的突变；和

(b3)上述b1～b2的组合。

在另一优选例中，执行步骤(b)的软件为SNPiR。

在另一优选例中，步骤(c)具体包括：将所述n个数据集Z按正负链拆分，得到n个正链子集和n个负链子集，在两种子集中分别计算所述候选RNA编辑位点的深度，当所述深度大于编辑质量的阈值，则判定所述编辑位点为可靠的编辑位点。

在另一优选例中，步骤(c)、(d)是使用SAMtools软件的view与mpileup指令执行的。

在另一优选例中，在第五数据集中，包括候选RNA编辑位点信息与编辑深度列表。

在另一优选例中，步骤(e)所述的统计分析还包括：对所述差异化RNA编辑位点作最终的选取。

在另一优选例中，步骤(a)具体包括：

(s1)提供来自n个样本的RNA测序结果，其中n为大于等于2的正整数，和所述n个样本的DNA的全基因组测序(WGS)结果；

(s2)将所述来自n个样本的RNA测序结果与参考基因组分别比对，获得n个RNA序列比对文件并校正；

(s3)将所述参考基因组与所述n个样本的全基因组测序(WGS)结果分别比对，获得n个DNA序列比对文件并校正；和

(s4)结合所述n个校正后的RNA序列比对文件和所述n个校正后的DNA序列比对文件，对每个样本分别执行联合突变检测与过滤，获得经联合突变检测的n个位点集，记为n个数据集Z。

在另一优选例中，执行步骤(s1)、(s2)、(s3)的软件包括：star、BWA、SAMtools、picard、或其组合。

在另一优选例中，执行所述的校正的软件包括picard

在另一优选例中，执行步骤(s4)的软件为GATK，具体地为mutect2。

在另一优选例中，步骤(s2)对所述的RNA序列比对文件的校正具体步骤包括：

(s2a)序列重新排列；

(s2b)标记样本库；

(s2c)标记重复序列

(s2d)分割含N片段；

(s2e)标记重新比对区域；

(s2f)局部重新比对；和

(s2g)BQSR校正。

在另一优选例中，步骤(3)对所述的DNA序列比对文件的校正具体步骤包括：

(s3a)序列重新排列；

(s3b)标记样本库；

(s3c)标记重复序列；

(s3d)标记重新比对区域；

(s3e)局部重新比对；和

(s3f)BQSR校正。

在另一优选例中，步骤(s2)与(s3)为各自独立地执行，可互换、先后、或同时进行。

在本发明的第二方面，提供了一种用于检测差异化RNA编辑位点的设备，所述设备包括：

(M1)输入模块，所述输入模块被配置为输入待测样本的经联合突变检测的RNA突变位点数据集Z；

(M2)处理模块，所述评估模块被配置为执行以下处理操作：将所述数据集Z中的RNA突变位点合并，获得第三数据集；将所述数据集Z进行基于RNA突变位点的测序深度的质量控制处理，获得第四数据集；将第三数据集与第四数据集合并，获得第五数据集；和，对第五数据集中的所述候选RNA编辑位点信息与编辑位点的测序深度数据进行统计分析，从而确定差异化RNA编辑位点；和

(M3)输出模块，所述输出模块被配置为输出所述差异化RNA编辑位点的信息。

在另一优选例中，所述处理模块(M2)包括以下子模块：

(M2a)初步合并子模块，所述评估模块被配置为执行以下处理操作：对于待检测差异化RNA编辑位点的两个或多个样本，将相应样本的各自的所述经联合突变检测的RNA突变位点数据集Z中的RNA突变位点进行合并，从而获得合并的候选RNA编辑位点数据集，记为第三数据集；

(M2b)质量控制子模块，所述评估模块被配置为执行以下处理操作：将所述经联合突变检测的RNA突变位点数据集Z，进行基于RNA突变位点的测序深度的质量控制处理，从而获得符合测序深度质量标准的编辑位点的数据集，记为第四数据集；

(M2c)质控后合并子模块，所述评估模块被配置为执行以下处理操作：将模块(M2a)生成的第三数据集中的合并的候选RNA编辑位点数据，与模块(M2b)生成的第四数据集中的符合测序深度质量标准的编辑位点的数据进行合并，从而获得含有候选RNA编辑位点信息数据和编辑位点的测序深度数据的第五数据集；和

(M2d)统计分析与位点确定子模块，所述评估模块被配置为执行以下处理操作：对所述第五数据集中的所述候选RNA编辑位点信息与编辑位点的测序深度数据进行统计分析，从而确定差异化RNA编辑位点。

在另一优选例中，所述的输出模块包括显示器、打印机、平板电脑(PAD)、智能手机。

在另一优选例中，所述的各模块通过有线或无线方式连接。

在另一优选例中，所述系统还包括(M4)控制模块：所述控制模块被配置为控制上述各模块和/或其子模块的运行。

为了本发明的目标，提供的方案可包含但不限于以下步骤：

对不同组别的样本集中的样本分别集进行链特异全转录组测序，获得RNA数据；对上述样本进行全基因组测序，获得DNA数据；将序列比对数据进行必要的预处理后，根据比对所在链拆分RNA测序数据；通过分析RNA数据和DNA数据，并进行联合突变检测，获得候选RNA编辑位点；使用这些位点，结合RNA数据获得这些位点编辑深度的数据；通过使用统计学手段评估这些位点的编辑程度是否受实验条件或组别划分的影响。

·其中，所述全转录组测序为RNA测序，优选的所述测序样本的制备方法为核糖体RNA去除法，优选的测序类型为链特异性测序(strans-specific sequencing)。

·其中，所述样本集的样本数量≥2个，优选的样本数量≥4个.

·其中，所述组别划分为不同实验处理方式或正常组织/肿瘤组织样本。

·其中，统计学手段包含：对原始或标准化处理后的数据使用T值检验、方差分析、基于负二项模型的差异分析以及基于广义线性模型差异分析，其中优选的统计手段为广义线性模型差异分析。

在优选例中，所述方法包含步骤：

(1)序列比对

1.1)获得不同组别样本集中每个样本的DNA和RNA的原始数据，并对其进行初步的质量检测与过滤获得有效数据；

1.2)结合参考转录组和RNA原始数据，对可变剪切位点进行检测或机器学习；

1.3)使用参考基因组对上述有效数据进行比对；

(2)RNA序列数据预处理

对RNA测序数据进行以下操作：依照突变检测参考基因组对测序数据重新排序；添加样本集标记；标记PCR产生重复序列；分割包含可变剪切位点的测序片段(含N片段)；定位出所有需要进行局部序列重比对的目标序列；对插入或缺失(indel)突变进行局部重新序列比对；矫正上述处理测序数据的碱基质量值。

(3)DNA序列数据预处理

对DNA测序数据进行以下操作：依照突变检测参考基因组对测序数据重新排序；添加样本集标记；标记PCR产生重复序列；定位出所有需要进行局部序列重比对的目标序列；对插入或缺失(indel)突变进行局部重新序列比对；矫正上述处理测序数据的碱基质量值。

(4)按照参考基因组的正负链拆分RNA序列比对文件

对每个样本的RNA序列比对文件拆分，依据为文件中每条序列所对应的参考基因组的正链或负链。

(5)联合DNA/RNA突变检测产生候选RNA编辑位点

通过联合突变检测方法，检测仅存在于RNA水平，但不存在于DNA水平的突变位点，即RNA候选编辑位点。

(6)对候选RNA编辑位点的进一步过滤

6.1)合并各样本集中的RNA候选编辑位点。

6.2)过滤检出的突变，仅保留单核苷酸突变。

6.3)过滤检出的突变，去除位于均聚核苷酸重复序列(homopolymer runs)中的突变。

6.4)过滤检出的突变，去除位于基因组重复序列中的突变。

(7)检测候选RNA编辑位点编辑深度

使用检测出的RNA突变坐标，对该坐标下与参考基因组一致的序列和突变的序列进行计数。计算样本集中每一个样本该位点突变序列占总体序列的比值。

(8)差异化RNA编辑位点的统计学分析

使用合适的统计学手段对样本集的数据进行两两对比或多组间对比分析，研究每个候选RNA编辑位点在不同样本集中的编辑程度是否存在统计学差异。

(9)差异化RNA编辑位点的选取

依据统计学分析结果过滤在编辑程度上不具备统计学意义的突变位点获得差异化RNA编辑位点。

应理解，在本发明范围内，本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合，从而构成新的或优选的技术方案。限于篇幅，在此不再一一赘述。

附图说明

图1为显示本发明的技术路线概览。

图2为显示本发明所描述的单个样本集的DNA/RNA联合突变检测RNA编辑位点的方法(步骤1-5)。

图3为显示本发明所描述的差异化RNA编辑位点分析方法(步骤6-9)。

图4为显示APOBEC3B在T-47D细胞株中的诱导表达情况。

图5为显示本发明检出RNA差异化编辑位点的变异等位基因变化分数与错误发现概率分布。

图6为显示本发明检出的A>I差异化RNA编辑位点与(下)与已知REDIportal收录的A>I编辑位点(中)、dbSNP数据库中A>G基因组SNV位点(上)的RNA折叠自由能分布对比。

图7为显示本发明检出的A>I差异化编辑位点上下游7bp核苷酸频率，以及和REDIportal收录的A>I编辑位点对比。

图8为显示本发明检出的C>U差异化RNA编辑位点(下)与已知APOBEC3B介导的C>TDNA编辑位点(中)、dbSNP数据库收录的C>T基因组SNV(上)的RNA折叠自由能分布对比。

图9为显示使用VaDiR方法(中)或本发明所描述方法(右)分别对人工加入差异化RNA编辑位点内部参照(左)的样本集进行RNA差异编辑位点检测结果。

图10为显示使用rMATS-DVR(左)或本发明所描述方法(右)分别对诱导表达APOBEC3B野生型(72小时诱导vs非诱导)样本集进行RNA差异编辑位点检测。

具体实施方式

本发明人经过广泛而深入的研究，首次、意外地开发出一种高灵敏、高置信的检测差异化RNA编辑位点的方法。具体地，经过本发明人的大量尝试、调整、和深入研究，获得的实验结果表明，使用所述方法相较过往捕获RNA编辑位点的方法显著提高了灵敏度；相较捕获差异化RNA编辑位点的方法显著提高了置信度。在此基础上，完成了本发明。

术语

如本文所用，术语“编辑位点的测序深度”与“编辑深度”可互换使用。

样本准备、核酸提取、和测序文库构建

在本发明中，实验设计包含有两个以上组别的样本集，组别可以是不同的实验处理方式、正常/肿瘤组织对照、以及时间序列等。

每个样本集需要两个以上重复，重复类型可以是生物学重复，在不具备收集生物学重复的条件下可采用技术重复。在优选例中，重复数量为4个，重复类型为生物学重复。

核酸提取方法选取上，可用常规的基因组DNA或全RNA提取方法，并达到一定的纯度。在优选例中，核酸的提取方法均为Roche MagnaPure 96仪器以及配套试剂；在优选例中，基因组DNA样本的A260/A280值达1.8-2.0之间，RNA样本的RIN值达8以上。

在全基因组测序文库构建方法的选取上，可选用和测序仪器配套的一般全基因组DNA文库构建方法。

在RNA测序文库的方法选取上，可以选取与所选测序仪器配套的一般文库构建方法。

在本发明中，特别适合的RNA测序文库类型是链特异文库类型，采取的方法可以是first strand法或second strand法。在优选例中采取的方法为first strand法。

RNA测序文库构建过程中，可选用核糖体去除法或者polyA富集法降低rRNA在文库中的含量，在本发明中，所采取的优选方法为核糖体去除法。

文库对应的类型可以是单端测序或双端测序，在优选例中所采取的方法为双端测序。

测序

在本发明中，可以采用测序文库对应的测序平台进行测序。优选的方法包括BGISeq、MGISeq、以及Illumina HiSeq等仪器。

所采取的测序深度应该在选取在合理范围。在优选例中，当采用人类细胞样本时，全基因组测序深度为30X以上，RNA测序深度为六千万reads/每样本。

所采取的测序读长应为100bp以上。在优选例中，在采用双端测序的情况下读长为100bp。

数据处理与差异化RNA编辑位点检测方法

如图1所示，发明的优选例中数据的处理与差异化RNA编辑位点的检测通常包含以下步骤：(1)序列比对；(2)RNA数据预处理；(3)DNA数据预处理；(4)按序列比对所在链拆分RNA测序文件；(5)通过联合突变检测候选RNA编辑位点；(6)候选RNA编辑位点的进一步过滤；(7)检测候选RNA编辑位点的编辑深度；(7)差异化RNA编辑的统计分析；(8)差异化RNA编辑位点的最终选取。

如图1所示，上述步骤可以分为以下两个部分：一、DNA/RNA联合突变检测候选RNA编辑位点(包含上述步骤1-5)；二、差异化RNA编辑位点分析(包含上述步骤6-8)。

其中第一部分的具体步骤如图2所示；第二部分的具体步骤如图3所示。

本发明的数据处理与差异化RNA编辑位点的检测方法，所涉及的软件可以为本领域采用的开源、公开、市售、以及操作者自行编写的软件、代码或脚本进行。

具体实施方法如下：

(1)序列比对：如图2所示，主要目的是由原始测序文件出发生成序列比对文件。

(1.1)原始测序数据的获取：所述原始数据包括多于两个样本集中每一个样本的RNA测序数据以及对应的DNA测序数据。在本发明的优选例中，原始数据包括：对照组RNA，对照组DNA，处理组RNA，处理组DNA的高通量测序数据，数据格式为FASTQ。

(1.2)序列过滤：该步骤包括了序列中测序接头以及低质量碱基的去除。在本发明的优选例中，所采用的软件为cutadapt。

(1.3)RNA序列比对：将RNA序列比对至参考基因组。在在本发明的优选例中，所采用的软件为star，优选方法为两轮法(star 2-pass)。

(1.4)DNA序列比对：将DNA序列比对至参考基因组。在本发明的优选例中，所采用的软件为Burrow Wheeler Aigner(BWA)，优选方法为BWA MEM法。

(2)RNA数据预处理：如图2所示，主要的目的是对序列比对结果进行必要的校正以适应后续的突变检测工作。这部分工作可以使用Picard Tools程序集进行。

(2.1)依照参考基因组对序列进行重排序。

(2.2)添加样本集标记。

(2.3)标记PCR产生重复序列。

(2.4)分割包含可变剪切位点的测序片段(含N片段)。

(2.5)定位出所有需要进行局部序列重比对的目标序列。

(2.6)对插入或缺失(indel)突变进行局部重新序列比对。

(2.7)对上述处理测序数据进行碱基质量评分校正(Base quality scorerecalibration,BQSR)。

(3)DNA数据预处理:如图2所示，主要的目的是对序列比对结果进行必要的校正以适应后续的突变检测工作。这部分工作可以使用Picard Tool进行。

(3.1)依照参考基因组对序列进行重排序。

(3.2)添加样本集标记。

(3.3)标记PCR产生重复序列。

(3.4)定位出所有需要进行局部序列重比对的目标序列。

(3.5)对插入或缺失(indel)突变进行局部重新序列比对。

(3.6)对上述处理测序数据进行BQSR校正。

(4)按链拆分测序比对文件：如图2所示，该过程旨在链特异拆分RNA测序比对文件，将每个样本分为正链(相对于参考基因组)比对文件和负链比对文件。

(5)联合突变检测：如图2所示，主要目的是仅检出存在于RNA序列中，却不存在于DNA序列中的突变。

(5.1)使用GATK的Mutect2工具对每个样本拆分的序列比对文件分辨进行突变检测，检测过程分为两组：分别是以处理组的DNA序列为“正常”参照，对处理组的RNA序列集进行突变检测；以及以对照组的DNA序列为“正常”参照，对对照组的RNA序列集进行突变检测。

(5.2)使用GATK的FilterMutectCalls工具对上述的突变检测结果进行分类。

(5.3)过滤未满足GATK的“PASS”条件的突变。发明的优选例采用的方法为bcftools软件中的view命令。

(5.4)过滤非单核苷酸突变。优选例中采用GATK工具包中的SelectVariant工具实现。

(6)候选RNA编辑位点的进一步过滤：如图3所示，这个过程旨在去除非单核苷酸突变以及容易产生假阳性的突变位点。

(6.1)去除均聚核苷酸重复序列中的突变，优选例中采用SNPiR程序包中的filter_homopolymer_nucleotides.pl实现。

(6.2)去除位于基因组重复序列中的突变。优选例中采用SNPiR程序包中的pblat_candidates_ln.pl实现。

(7)检测编辑深度：如图3所示，该过程旨在链特异地检测各样本集中每个样本在候选RNA编辑位点上的编辑深度。

(7.1)获取候选RNA编辑位点所对应的突变序列与正常序列的计数，并计算编辑深度，发明的优选例采用的方法为samtools软件中的mpileup命令。

(7.2)收集各样本信息并编制候选RNA编辑位点信息与编辑深度列表，该列表应包含至少以下信息：候选RNA编辑位点所在染色体、碱基位置、该位置参考基因组碱基类型(REF base)、该位置被编辑后的碱基类型(ALT base)、RNA编辑位点所在链、以及各样本集中每个样本在该编辑位点的编辑情况。其中后者包括参考碱基深度、碱基编辑深度、以及含碱基编辑序列占该位点所有序列的比例(碱基编辑比例)。

(8)差异化编辑统计分析：如图3所示，该过程旨在使用edgeR等程序包对RNA突变位点的差异化编辑分析以及高置信度位点的筛选。

(8.1)将RNA编辑位点信息与编辑深度列表信息导入edgeR软件包，并且使用对数似然比检验(log likelihood ratio test,LRT)方法比较候选RNA编辑位点的编辑比例在各样本集的对比中是否存在统计学差异，并计算单个位点的错误发现率(false discoveryrate,FDR)。

(9)差异化RNA编辑位点的最终选取：如图3所示，该过程的选取条件一般是测序深度大于10，并且统计结果错误发现率(FDR)小于5％的RNA编辑位点作为最终确认的差异化RNA编辑位点。

计算机指令与参数等的示例

上述处理中，步骤(1.2)使用开源软件cutadapt，下载地址为

https://github.com/marcelm/cutadapt/releases，命令行范式如下：

·cutadapt-g ADAPTER-O 5-e 0-o sample.trimmed.fastq sample.fastq

--minimum-length 35--discard-untrimmed--info-file＝reads.adapter.txt

·cutadapt-q 10-o output.fastq input.fastq

步骤(1.3)使用开源软件star，下载地址为https://github.com/alexdobin/STAR，命令行范式如下：

·STAR--runThreadN 20--runMode genomeGenerate--genomeDir$GENOME_DIR--genomeFastaFiles$REFERENCE.FA--sjdbGTFfile$GENCODE.GTF--sjdbOverhang 99

·STAR--runThreadN 20--genomeDir$GENOME_DIR--outFileNamePrefix$PREFIX--outSJfilterReads Unique–readFilesIn$DIR

·STAR--runThreadN 20--runMode genomeGenerate--genomeDir$GENOME_DIR--genomeFastaFiles$REFERENCE.FA--sjdbGTFfile$GENCODE.GTF--sjdbFileChrStartEnd$PREFIX_SJ.out.tab--sjdbOverhang 99--limitSjdbInsertNsj 1583333

·STAR--runThreadN 5--genomeDir$GENOME_DIR--outFileNamePrefix

$PREFIX_RNA--readFilesIn 1.fq.gz 2.fq.gz--outSJfilterReads Unique--readFilesCommand zcat--outFilterMultimapNmax 1

步骤(1.4)采用开源软件BWA和SAMtools，下载地址为https://github.com/lh3/bwa和www.htslib.org/download，命令行范式如下：

·bwa mem-t 20$REFERENCE.FA 1.fq 2.fq|samtools view-bT$REFERENCE.FA>

GENOMIC_DNA.bam

步骤(2.1)采用开源软件picard tools，下载地址为(https://broadinstitute.github.io/picard),命令行范式如下：

·java-Xmx4g-jar picard.jar ReorderSam INPUT＝$RNA_INPUT.BAM OUTPUT＝$RNA_REORDERED.BAM S＝true R＝$REFERENCE.FA

步骤(2.2)采用开源软件picard tools，命令行范式如下：

·java-Xmx4g-jar picard.jar AddOrReplaceReadGroups

INPUT＝$RNA_REORDERED.BAM OUTPUT＝$RNA_ADDRG.BAM RGID＝$LABEL_RNARGLB＝$LABEL_RNA RGPL＝ILLUMINA RGPU＝lane1 RGSM＝$LABEL_RNA

步骤(2.3)采用开源软件picard tools，命令行范式如下：

·java-Xmx4g-jar picard.jar MarkDuplicates INPUT＝$RNA_ADDRG.BAMOUTPUT＝$RNA_DEDUP.BAM CREATE_INDEX＝true

VALIDATION_STRINGENCY＝SILENT READ_NAME_REGEX＝null

METRICS_FILE＝$LABEL_RNA_METRICS.TXT

步骤(2.4)采用开源软件GATK,下载地址为(https://gatk.broadinstitute.org)，命令行范式如下：

·java-Xmx4g-jar GenomeAnalysisTK.jar-T SplitNCigarReads-R$REFERENCE.FA-I$RNA_DEDUP.BAM-o$RNA_SPLIT.BAM-rf ReassignOneMappingQuality-RMQF 255

-RMQT 60-U ALLOW_N_CIGAR_READS'

步骤(2.5)采用开源软件GATK,结合已知Indel数据库(下载地址为https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle)信息,命令行范式如下：

·java-Xmx4g-jar GenomeAnalysisTK.jar-T RealignerTargetCreator-R

$REFERENCE.FA-I$RNA_SPLIT.BAM-o$REALIGNER_INTERVALS-known

$KNOW_INDELS

步骤(2.6)采用开源软件GATK，命令行范式如下：

·java-Xmx4g-jar GenomeAnalysisTK.jar-T IndelRealigner-R$REFERENCE.FA-I

$RNA_SPLIT.BAM-o$RNA_INDEL.BAM-targetIntervals$REALIGNER_INTERVALS-known$KNOW_INDELS

步骤(2.7)采用开源软件GATK，命令行范式如下：

·java-Xmx4g-jar GenomeAnalysisTK.jar-T BaseRecalibrator-I$RNA_INDEL.BAM-R$REFERENCE.FA-o$RECALIBRATION_REPORT.GRP-knownSites$KNOW_INDELS

·java-Xmx4g-jar GenomeAnalysisTK.jar-T PrintReads-R$REFERENCE.FA-I

$RNA_INDEL.BAM-BQSR$RECALIBRATION_REPORT.GRP-o

$RNA_RECALIBRATED.BAM-U ALLOW_N_CIGAR_READS

步骤(3.1)采用开源软件picard tools，下载地址为命令行范式如下：

·java-Xmx4g-jar picard.jar ReorderSam INPUT＝$DNA_INPUT.BAM OUTPUT＝$DNA_REORDERED.BAM S＝true R＝$REFERENCE.FA

步骤(3.2)采用开源软件picard tools，命令行范式如下：

·java-Xmx4g-jar picard.jar AddOrReplaceReadGroups

INPUT＝$DNA_REORDERED.BAM OUTPUT＝$DNA_ADDRG.BAM RGID＝$LABEL_DNARGLB＝$LABEL_DNA RGPL＝ILLUMINA RGPU＝lane1 RGSM＝$LABEL_DNA

步骤(3.3)采用开源软件picard tools，命令行范式如下：

·java-Xmx4g-jar picard.jar MarkDuplicates INPUT＝$DNA_ADDRG.BAMOUTPUT＝$DNA_DEDUP.BAM CREATE_INDEX＝true

VALIDATION_STRINGENCY＝SILENT READ_NAME_REGEX＝null

METRICS_FILE＝$LABEL_DNA_METRICS.TXT

步骤(3.4)采用开源软件GATK,结合已知Indel数据库(下载地址为https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle)信息,命令行范式如下：

·java-Xmx4g-jar GenomeAnalysisTK.jar-T RealignerTargetCreator-R

$REFERENCE.FA-I$DNA_SPLIT.BAM-o$REALIGNER_INTERVALS-known

$KNOW_INDELS

步骤(3.5)采用开源软件GATK，命令行范式如下：

java-Xmx4g-jar GenomeAnalysisTK.jar-T IndelRealigner-R$REFERENCE.FA-I

$DNA_SPLIT.BAM-o$DNA_INDEL.BAM-targetIntervals$REALIGNER_INTERVALS-known$KNOW_INDELS

步骤(3.6)采用开源软件GATK，命令行范式如下：

·java-Xmx4g-jar GenomeAnalysisTK.jar-T BaseRecalibrator-I$DNA_INDEL.BAM-R$REFERENCE.FA-o$RECALIBRATION_REPORT.GRP-knownSites$KNOW_INDELS

·java-Xmx4g-jar GenomeAnalysisTK.jar-T PrintReads-R$REFERENCE.FA-I

$DNA_INDEL.BAM-BQSR$RECALIBRATION_REPORT.GRP-o

$DNA_RECALIBRATED.BAM-U ALLOW_N_CIGAR_READS

步骤(4.1)采用开源软件samtools，下载地址为www.htslib.org/download，命令行范式如下：

·samtools view-f 64-b$RNA_RECALIBRATED.BAM>

$RNA_RECALIBRATED_FIRST_END.BAM-@20

·samtools view-f 128-b$RNA_RECALIBRATED.BAM>

$RNA_RECALIBRATED_SECOND_END.BAM-@20

步骤(5.1)采用开源软件GATK，结合gromAD人群生殖系资源(下载地址为https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle)，命令行范式如下：

·gatk Mutect2-R$REFERENCE.FA-I$LIST_OF_RNA_RECALIBRATED.BAM-I

$DNA_RECALIBRATED.BAM-normal$LABEL_DNA--germline-resource$GNOMAD.VCF-O$OUTPUT.VCF--native-pair-hmm-threads 20

步骤(5.2)采用开源软件GATK，命令行范式如下：

·gatk FilterMutectCalls-V$OUTPUT.VCF-R$REFERENCE.FA-O

$OUTPUT_WITH_LABEL.VCF

步骤(5.3)采用开源软件bcftools,下载地址为(https://samtools.github.io/bcftools/bcftools.html)，命令行范式如下：

·bcftools view-f PASS$OUTPUT_WITH_LABEL.VCF>$OUTPUT_PASS.VCF

步骤(5.4)采用开源软件GATK，命令行范式如下：

gatk SelectVariants-R$REFERENCE.FA-V$OUTPUT_PASS.VCF--select-type-to-include SNP-O$OUTPUT_SNV.VCF

步骤(6.1)采用开源软件SNPiR，下载地址为(https://github.com/rpiskol/SNPiR)，命令行范式如下：

·revised_convertVCF.sh$OUTPUT_SNV.VCF$OUTPUT_SNPIR.VCF

·perl filter_homopolymer_nucleotides.pl-infile$OUTPUT_SNPIR.VCF-outfile

$OUTPUT_HOMO.VCF-refgenome$REFERENCE.FA

步骤(6.2)采用开源软件SNPiR、pBLAT(下载地址为https://github.com/icebert/pblat)和bedTools(下载地址为https://github.com/arq5x/bedtools2)，命令行范式如下：

·perl pblat_candidates_ln.pl-infile$OUTPUT_HOMO.VCF-outfile

$OUTPUT_PBLAT.VCF-bamfile$RNA_MERGED.BAM-refgenome$REFERENCE.FA-pblatpath$PATH_TO_PBLAT-threads 20>pblat_error.txt 2>&1'

·bedtools intersect-a$OUTPUT_SNV.VCF-b$OUTPUT_PBLAT.VCF-wa-header>

$OUTPUT_FINAL.VCF

步骤(7.1)采用开源软件samtools，命令行范式如下：

·samtools mpileup-B-d 100000-f$REFERENCE.FA-l$OUTPUT_FINAL.VCF-q 30-Q 17-a-o$RNA_RECALIBRATED.MPILEUP$LIST_OF_RNA_RECALIBRATED.BAM

其中，输出文件$RNA_RECALIBRATED.MPILEUP包含了样本集中每个候选RNA编辑位点深度信息。通过使用常见的数据处理分析软件，如R，并结合候选RNA编辑位点

$OUTPUT_FINAL.VCF等信息，可编制候选RNA编辑位点信息与编辑深度列表。

步骤(8.1)采用开源软件edgeR，下载地址为(http://bioconductor.org)。使用R语言脚本读取候选RNA编辑位点信息与编辑深度列表，并使用edgeR默认的均一化和负二项分布模型，对样本集之间进行两两对比，或者在此基础上利用广义线性模型对多个样本集进行对比。命令行范式参照edgeR用户手册，调集指令如下：

·library(edgeR)；edgeRUsersGuide()

步骤(9.1)中，研究者通过具体需求选取合适的FDR值以及平均突变深度变化，可对候选差异化RNA编辑位点进行筛选。一般选取FDR≤0.05作为评判差异化RNA编辑位点的先决条件。

上述步骤的典型系统运行环境和软件版本为：Linux CentOS 7；python 2.7.11；perl；gatk3.6；picard tools；samtools；bcftools；bedtools；cutadapt；star；bwa；SNPiR；R-4.2.1；edgeR-3.38.4。

与现有的仅通过DNA/RNA联合突变检测RNA编辑位点的手段(如VaDiR方法)相比，本发明结合近期出现的高通量测序新技术和生物信息学分析新方法,引入了差异化编辑分析手段，提高了RNA编辑位点的检测灵敏度。在实施例中(实施例7、图9)，通过在标准数据中掺入内部参照后分析可发现，本发明对高突变深度的内参位点检出度较现有方法有所提升的同时，对低突变深度的内参位点也保持了良好的灵敏度，因此在RNA编辑位点的检出上，性能较前人方法有显著提升。

与现有的仅通过差异化统计学分析的RNA编辑位点检测手段(如rMATS-DVR方法)相比，本发明引入了通过对比转录组和基因组联合突变检测手段，有效的过滤了基因组突变带来的差异化RNA突变，因此排除了前人方法中主要的假阳性来源，大大的提高了检出结果的置信度，并且进一步过滤了其他常见的差异化RNA编辑位点来源，比如单核苷酸重复序列和简单重复序列，进一步提高了检出结果的准确性。在实施例中(实施例8、图10)，检测所得RNA差异化编辑位点较rMATS-DVR方法避免了因为DNA层面上的突变导致的假阳性。

本发明的主要优点

(1)首次结合了联合突变检测与差异化分析方法检测RNA差异化编辑位点。该手段较仅通过差异化分析的检测手段相比，过滤了DNA突变带来的假阳性，大大提高了检出的置信度。

(2)较过往的通过联合突变检测RNA编辑位点方法相比，具有高灵敏度特征。

(3)提供了一种在全转录组层面上研究实验处理条件、时间序列以及肿瘤形成对RNA特异位点编辑影响的手段。

下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明具体条件的实验方法，通常按照常规条件，例如Sambrook等人，分子克隆：实验室手册(New York:Cold Spring HarborLaboratory Press,1989)中所述的条件，或按照制造厂商所建议的条件。除非另外说明，否则百分比和份数是重量百分比和重量份数。

实施例1：慢病毒介导诱导表达胞苷脱氨酶APOBEC3B细胞株的构建

在本实施例中，构建慢病毒介导诱导表达APOBEC3B以及无活性突变对照APOBEC3B(E68Q/E225Q)细胞株。该酶是已知的RNA胞苷(C)至尿苷(U)编辑酶，能作为考察RNA编辑水平的工具。具体操作方法如下：

使用GeneArt服务(Thermo Fisher)全基因合成APOBEC3B片段(NM_004900.5)，在合成前对编码序列进行人类细胞密码子优化。所以用QuickChange试剂盒(Agilent)对APOBEC3B的E68、E225分别进行突变，得到APOBEC3B(E68Q/E225Q)编码DNA片段。将这两个DNA片段使用AgeI和ClaI(BspDI)分别克隆至pTRIPZ(GE Healthcare)载体中，得到慢病毒诱导表达质粒。

使用辅助质粒pMD2.G和psPAX2与上述表达质粒共转染293TN(Clonetech)细胞，24小时后更换一次培养基。并于转染后48和72小时收集两次培养基上清。经过0.44μm微孔过滤后，使用Peg-IT试剂盒(System Biosciences)的标准流程浓缩慢病毒。

使用浓缩过的慢病毒转导T-47D细胞，于转导后48小时更换培养基，并于24小时后于培养基中加入1μg/ml嘌呤霉素，并保持该选择压力对细胞持续进行培养。

实施例2：胞苷脱氨酶APOBEC3B的诱导表达与核酸样本制备

在本实施例中，介绍诱导表达野生型与突变体APOEBC3B之前与之后的细胞样本的制备方法，以及DNA与RNA样本的制备方法。具体操作如下：

取上述能诱导表达APOBEC3B野生型以及突变体的T-47D细胞，并培养至对数生长期。期间取4个RNA生物学重复样本，1个DNA样本。

使用1μg/ml多西环素诱导APOBEC3B野生型以及突变体在T-47D细胞中表达72小时。期间取4个RNA生物学重复样本，1个DNA样本。诱导蛋白表达效果如图4所示

取样操作如下：使用PBS清洗细胞两次，并使用刮铲将细胞刮下。-80°C保存。

核酸样本制备方法为使用Roche MagNA Pure 96仪器(Roche)。DNA样本使用DNAand Viral NA Small Volume Kit试剂盒制备；RNA样本采用Cellular RNA Large VolumeKit试剂盒制备。制备后的DNA样本OD260/280＝1.85左右，RNA样本RIN值(使用GEBioAnalyzer2100测定)为10左右。

实施例3：测序文库的制备与测序

取上述实施例中的核酸样本，送至华大基因(深圳)公司进行测序文库的制备。文库均使用华大基因公司基于DNA纳米球测序技术的文库构建方法构建。其中RNA文库制备过程中，采用核糖体RNA去除法富集非核糖体RNA，并采用链特异的文库构建手段(first-strand)。

文库测序仪器选用BGIseq-500型测序仪。全基因组测序采用双端测序方法，读长100bp，RNA测序采用双端测序方法，读长100bp。

测序深度为全基因组测序每样本30X以上(人类基因组)，RNA测序为每样本六千万reads。

实施例4：测序数据的预处理与RNA突变位点检测

使用“具体实施方式”中所述的序列比对方法，用cutadapt对DNA与RNA测序数据进行去接头处理后，使用bwa软件对全基因组测序所得序列使用人类GRCh38参考基因组进行比对；并使用star软件中的2-pass方法对RNA测序所的序列使用人类GRCh38参考基因组和Gencode39参考转录组进行比对。其中RNA测序中样本成功比对的比率范围为94-98％。

使用“具体实施方式”中所述的RNA数据预处理方法，结合GRCh38参考基因组，以及1000G_phase1.indels.b37.vcf和

Mills_and_1000G_gold_standard.indels.b37.sites.vcf两个公开文件组成的indel参考数据库对RNA样本的序列进行BQSR校正。使用“具体实施方式”中所述的DNA数据预处理方法，对DNA样本，结合GRCh38参考基因组，以及1000G_phase1.indels.b37.vcf和

Mills_and_1000G_gold_standard.indels.b37.sites.vcf两个公开文件组成的indel参考数据库进行BQSR校正。

使用“具体实施方式”中所述的突变检测与初步方法，结合经BQSR校正后的DNA与RNA测序数据，进行RNA突变位点检测，其中使用GRCh38作为参考基因组，gromAD人群生殖系资源使用af-only-gnomad.hg38.vcf.gz公开文件。

本章节所涉及公开文件与数据资源均从gatk网站获得，地址为(https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle).

根据上述的操作流程，用Mutect2软件对比APOBEC3B野生型和突变体在非诱导条件下和诱导条件下的样本集，得到候选突变集。

对这些候选突变集进行过滤筛选后，结果汇总如下表(表1)：

表1对候选突变集过滤筛选后的结果汇总

实施例5：APOBEC3B介导的差异化RNA编辑位点的筛选

使用“具体实施方式”中所述的样本中RNA突变位点的序列计数方法，构建包含所有样本集的RNA编辑位点信息的矩阵文件，并使用edgeR程序包中的广义线性模型，对每个突变位点所对应所有样本的突变深度进行对数秩检验。

选取经对数秩检验后FDR≤0.05的位点，并且要求RNA编辑位点所在链与RefGene标注基因所在链一致的情况下，为A>G(对应RNA层面上A>I编辑)以及C>T(对应RNA层面上C>U编辑)，结果如图5所示：

对比结果可知，对比失活突变体APOBEC3B，当在T-47D细胞中诱导表达野生型APOBEC3B时，能诱导大量的C>U差异化RNA编辑位点的产生，而这些差异化位点的编辑程度变化几乎全部为显著上升。因此，证明这些差异化的C>U的产生时APOBEC3B的RNA胞苷脱氨酶活性介导的。

综上所述，本发明能有效检出某种实验条件或介入手段导致的差异化RNA编辑位点。

实施例6：差异化RNA编辑位点的验证

本实施例采用评估RNA折叠自由能的方法对所获差异化RNA编辑位点进行评估验证。一般地，处于RNA序列中的单核苷酸位点有较低的折叠自由能，而处于基因组非编码序列中的单核苷酸位点由于并非转录为RNA序列，有较高的折叠自由能。本实施例中RNA折叠自由能由RNAfold软件以及R语言软件包“LncFinder”进行，这两款软件的下载地址为

https://github.com/choener/RNAfold和https://www.bioconductor.org/

通过R语言软件包"BSgenome.Hsapiens.UCSC.hg38"获取差异化RNA编辑位点以及其他位点上下游100bp的序列，其中R语言命令行为：

RNA$V6<-getSeq(Hsapiens,RNA$V1,start＝RNA$V2-100,width＝201,as.character＝T,strand＝RNA$V5)

其中，输入文件为tab分隔的文本文档，数据为bed格式，其中第一列为编辑位点所在染色体，第二列为编辑位点的坐标，第五列为编辑位点所在链。

将上述所得文件转化为FASTA格式后，使用下述R语言命令范式使用”LncFinder”软件包进行RNAfold分析：

SS.seq<-run_RNAfold(RNA_FASTA,RNAfold.path＝RNAfold.path,parallel.cores＝20)

通过比对本发明所获A>I差异化编辑位点，对比RARAR数据库经实验验证的A>I差异化位点，以及dbSNP数据库保存的SNV位点。结果显示，本发明检出的A>I差异化编辑位点与实验验证的A>I位点在RNA折叠自由能上没有显著差别，但是与SNV位点有明显差别(图6)。

其中本实施例采用的dbSNP数据库版本为150，下载地址为：

https://www.ncbi.nlm.nih.gov/projects/SNP/snp_summary.cgi

使用R语言软件包"BSgenome.Hsapiens.UCSC.hg38“获得A>I差异化RNA编辑位点上下游各7bp的序列信息，并使用WebLogo3对sequence logo进行绘制。结果如图7所示。

在发明所获A>I差异化编辑位点的序列特异性上，发现与实验所得的A>I差异化编辑位点基本一致(图7)。

通过比对本发明所获APOBEC3B依赖的C>U差异化编辑位点，对比公共数据库GEO中保存的测定APOBEC3B依赖的C>T SNV突变(GSE193225实验,原始数据下载地址

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi？acc＝GSE193225)，以及dbsnp数据库保存的SNV位点。结果显示，本发明检出的C>U差异化编辑位点与实验验证的APOBEC3B在基因组DNA层面编辑的C>T SNV突变位点存在显著差别，并且与dbsnp数据库中保存的C>T SNV位点存在显著差别图8。

这些结果显示捕获的差异化RNA位点与DNA层面的SNV位点有明显的差别，表明了所捕获的位点具备已知的RNA特征以及A>I序列特异性。

实施例7：与既有RNA编辑位点检测手段的对比

本实施例旨在对比本发明与前人开发的不基于差异化编辑分析的RNA编辑位点检出手段(PMID:29267927)。

使用公共数据库GEO中保存的T-47D细胞的RNA测序原始数据，选取GSE193234实验的GSM5777068-GSM5777072号样本(4个生物学重复)作为一个样本集，以及GSM5777076-GSM5777080号样本(4个生物学重复)作为另一样本集。所有样本集原始数据的下载地址为

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi？acc＝GSE193234。

将上述RNA测序样本使用本发明所述的序列比对方法，比对至GRCh38参考基因组，比对时使用Gencode39转录组注释进行辅助比对。

在包含GSM5777076-GSM5777080号样本集中加入突变内参集。为了作对比，突变内参集分为两组，即高突变深度组和低突变深度组。突变内参加入采用bamsurgeon开源软件，软件的下载地址为

https://github.com/adamewing/bamsurgeon。所采用命令行范式为：pythonaddsnv.py-v$VARIANT.TXT-f$INPUT_BAM.BAM-r$REFERECE.FA-o

$OUTPUT_BAM.BAM--tmpdir$TEMP_DIR--ignorepileup--force

其中突变内参集中记录的位点满足在RNA测序样本中的平均深度大于10个reads的条件。

使用公共数据库GEO中保存的T-47D细胞的全基因组测序数据(样本GSM5777027，下载地址为

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi？acc＝GSM5777027)作为参考基因组，对添加了内参的GSM5777076-GSM5777080号样本分别使用VaDiR方法进行RNA突变检测。其中VaDiR方法于Neums et al.,2018详细介绍，该方法仅对比RNA测序结果和DNA测序结果，不涉及差异化编辑分析。VaDiR软件下载地址为https://scicrunch.org/resolver/RRID:SCR_015797

使用本发明所述的方法，结合T-47D细胞的全基因组测序数据(GSM5777027),对GSM5777068-GSM5777072号样本和GSM5777076-GSM5777080样本进行差异化RNA编辑位点分析。对比VaDiR方法，结果如图9所示。

结果显示，本发明所采用的方法无论是在高突变深度内参组别，还是低突变深度内参组别，均相对VaDiR方法有较低的漏检率，因此相对该仅对比基因组和转录组测序结果的方法，本发明所述的引入差异化分析方法对检出RNA编辑位点有更高的灵敏度。

实施例8：与既有差异化RNA编辑位点检测手段的对比

本实施例旨在对比本发明与前人开发的仅基于RNA测序数据的差异化RNA编辑位点分析手段(PMID:28334241)。

使用rMATS-DVR软件，对实施例1-3所获的，除了DNA测序以外的其他数据进行分析。其中，采用了与标准rMATS-DVR分析手段一致，使用了gatk中的UnifiedGenotyper手段，并且仅仅使用RNA测序数据进行突变检测。本实施例仅对比诱导表达APOBEC3B野生型与非诱导样本之间

检测结果如图10所示。

从结果可知，本发明通过在差异化RNA编辑位点分析中引入基因组DNA作为参照，能排除基因组突变带来的干扰，而且能过滤位于均聚核苷酸序列以及简单重复序列中的假阳性结果。

在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims

1.一种检测差异化RNA编辑位点的方法，其特征在于，包括步骤：

其中，步骤(b)和(c)的步骤可互换或同时进行。

2.如权利要求1所述方法，其特征在于，所述的方法是非诊断和非治疗性的。

3.如权利要求1所述方法，其特征在于，在步骤(b)中，对相应样本的各自的经联合突变检测的RNA突变位点数据集Z中的RNA突变位点进行合并后，还进行过滤处理，从而获得经过滤的合并的候选RNA编辑位点数据。

4.如权利要求3所述方法，其特征在于，步骤(b)所述的过滤的操作选自下组：

(b1)去除均聚核苷酸序列中的突变；

(b2)去除位于基因组重复序列中的突变；和

(b3)上述b1和b2的组合。

5.如权利要求1所述方法，其特征在于，执行步骤(b)的软件为SNPiR。

6.如权利要求1所述方法，其特征在于，步骤(a)具体包括：

7.如权利要求6所述方法，其特征在于，步骤(c)具体包括：将所述N个数据集Z按正负链拆分，得到N个正链子集和N个负链子集，在两种子集中分别计算所述候选RNA编辑位点的深度，当所述深度大于编辑质量的阈值，则判定所述编辑位点为可靠的编辑位点。

8.如权利要求1所述方法，其特征在于，步骤(c)、(d)是使用SAMtools软件的view与mpileup指令执行的。

9.如权利要求1所述方法，其特征在于，在第五数据集中，包括候选RNA编辑位点信息与编辑深度列表。

10.如权利要求1所述方法，其特征在于，步骤(e)所述的统计分析还包括：对所述差异化RNA编辑位点作最终的选取。

11.如权利要求6所述方法，其特征在于，执行步骤(s1)、(s2)、(s3)的软件包括：star、BWA、SAMtools、picard、或其组合。

12.如权利要求11所述方法，其特征在于，执行所述的校正的软件包括picard。

13.如权利要求6所述方法，其特征在于，执行步骤(s4)的软件为GATK，具体地为mutect2。

14.如权利要求6所述方法，其特征在于，步骤(s2)对所述的RNA序列比对文件的校正具体步骤包括：

(s2a)序列重新排列；

(s2b)标记样本库；

(s2c)标记重复序列

(s2d)分割含N片段；

(s2e)标记重新比对区域；

(s2f)局部重新比对；和

(s2g)BQSR校正。

15.如权利要求6所述方法，其特征在于，步骤(s3)对所述的DNA序列比对文件的校正具体步骤包括：

(s3a)序列重新排列；

(s3b)标记样本库；

(s3c)标记重复序列；

(s3d)标记重新比对区域；

(s3e)局部重新比对；和

(s3f)BQSR校正。

16.如权利要求6所述方法，其特征在于，步骤(s2)与(s3)为各自独立地执行，可互换或同时进行。

17.一种用于检测差异化RNA编辑位点的设备，其特征在于，所述设备包括：

(M2)处理模块，所述处理模块被配置为执行以下处理操作：将所述数据集Z中的RNA突变位点合并，获得第三数据集；将所述数据集Z进行基于RNA突变位点的测序深度的质量控制处理，获得第四数据集；将第三数据集与第四数据集合并，获得第五数据集；和，对第五数据集中的候选RNA编辑位点信息与编辑位点的测序深度数据进行统计分析，从而确定差异化RNA编辑位点；和

18.如权利要求17所述设备，其特征在于，所述处理模块(M2)包括以下子模块：

(M2a)初步合并子模块，所述初步合并子模块被配置为执行以下处理操作：对于待检测差异化RNA编辑位点的两个或多个样本，将相应样本的各自的所述经联合突变检测的RNA突变位点数据集Z中的RNA突变位点进行合并，从而获得合并的候选RNA编辑位点数据集，记为第三数据集；

(M2b)质量控制子模块，所述质量控制子模块被配置为执行以下处理操作：将所述经联合突变检测的RNA突变位点数据集Z，进行基于RNA突变位点的测序深度的质量控制处理，从而获得符合测序深度质量标准的编辑位点的数据集，记为第四数据集；

(M2c)质控后合并子模块，所述质控后合并子模块被配置为执行以下处理操作：将模块(M2a)生成的第三数据集中的合并的候选RNA编辑位点数据，与模块(M2b)生成的第四数据集中的符合测序深度质量标准的编辑位点的数据进行合并，从而获得含有候选RNA编辑位点信息数据和编辑位点的测序深度数据的第五数据集；和

(M2d)统计分析与位点确定子模块，所述统计分析与位点确定子模块被配置为执行以下处理操作：对所述第五数据集中的所述候选RNA编辑位点信息与编辑位点的测序深度数据进行统计分析，从而确定差异化RNA编辑位点。

19.如权利要求17所述设备，其特征在于，所述的输出模块包括显示器、打印机、平板电脑(PAD)、智能手机。

20.如权利要求17所述设备，其特征在于，所述的各模块通过有线或无线方式连接。

21.如权利要求17所述设备，其特征在于，所述设备还包括(M4)控制模块：所述控制模块被配置为控制上述各模块和/或其子模块的运行。