一种鉴定SL序列反式剪切所形成的基因编码框的方法

技术领域

本发明涉及生物技术领域，特别是涉及一种鉴定SL序列反式剪切所形成的基因编码框的方法。

背景技术

基因编码框(Openreading frame，ORF)在基因组中的注释对于后续的生物学研究和应用至关重要。近年来，研究人员开发了各种算法来预测基因组中的ORF，但是这些算法和工具均是以模式物种的研究为出发点所进行的，虽然在大部分物种的研究中都是适用的，但是一些特殊物种的研究需求则没有被考虑。

在自然界中，大部分真核生物基因转录后形成信使RNA(mRNA)，随后通过顺式剪切去掉内含子(Intron)以后形成最终的成熟mRNA，进行翻译，合成蛋白质。然而，有另外一类真核生物，它们的mRNA成熟过程中除了需要经历内含子的顺式剪切，还需要经历前导(spliced leader，SL)序列的反式剪切。具体来说，在这类生物中，基因转录成mRNA以后，在其5’端进行反式剪切，添加上一段特定的SL序列，然后才能形成成熟的mRNA，进行翻译。由于很多物种中SL序列自身带有启始密码子(如“AUG”、“UUG”等)，在现有技术存在SL序列的添加将有可能引入新的翻译起始位点，从而改变原有的基因编码框，形成新的蛋白质序列，使得原始的基因组注释信息失效，从而引发一系列的错误分析和结论的问题。

发明内容

为了克服现有技术的不足，本发明的目的是提供一种鉴定SL序列反式剪切所形成的基因编码框的方法，本发明解决了现有技术存在SL序列的添加将有可能引入新的翻译起始位点，从而改变原有的基因编码框，形成新的蛋白质序列，使得原始的基因组注释信息失效，从而引发一系列的错误分析和结论的问题。

为实现上述目的，本发明提供了如下方案：

一种鉴定SL序列反式剪切所形成的基因编码框的方法，包括：

获得全部的基因组所有转录本和已知编码框序列；

根据所述转录本和已知编码框序列获取所述基因组的比对文件；

基于转录组测序技术，根据所述比对文件筛选局部比对的序列，并提取未成功比对上的序列；

对所述未成功比对上的序列与指定的SL序列进行比较，若所提取的序列为指定SL序列的末端至少8个碱基序列，则确定该序列为SL剪切位点所产生的序列；

根据所述比对文件中所记录的位置信息推算SL序列的剪切位点，并将所述剪切位点转换成转录组坐标；

根据所述转录组坐标和所述转录本得到完整的mRNA序列；

利用核糖体印迹测序数据，根据所述完整的mRNA序列获得完整的mRNA对应的编码框。

优选地，所述位置信息包括：

第3列的比对染色体信息和第4列的比对坐标。

优选地，所述根据所述转录组坐标和所述转录本得到完整的mRNA序列包括：

将SL剪切位点以前的序列删除，并替换成指定的SL序列以得到完整的mRNA序列。

优选地，根据基因组中所有编码框的密码子使用情况，计算每个密码子在整体基因组编码序列中的出现频率，然后计算每个已知编码框中密码子频率的平均值，并将其转化为Z-score。

优选地，根据所述完整的mRNA序列获得完整的mRNA对应的编码框包括：

基于核糖体印迹测序技术，获取核糖体印迹数据序列；

对所述核糖体印迹数据序列筛选，获得满足条件的核糖体印迹数据序列；

根据所述满足条件的核糖体印迹数据序列获得5’端与P-site之间不同距离的出现频率；

根据每一个碱基位于P-site的概率和所述已知编码框中密码子频率的平均值，对编码框进行预测，得到预测编码框；

对所述预测编码框进行筛选，满足第一条件的预测编码框，进行输出，得到编码框。

优选地，所述每一个碱基位于P-site的概率的包括：

5’端与P-site之间不同距离的出现频率和所述满足条件的核糖体印迹数据序列比对位置。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供了一种鉴定SL序列反式剪切所形成的基因编码框的方法，本发明通过确定SL剪切点所产生的序列，通过确定SL剪切点所产生的序列获得完整的mRNA序列，从而确定所对应的基因编码框，本发明提升了SL序列添加的准确性，降低了因SL序列添加使得得到编码框不准的可能性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的鉴定SL序列反式剪切所形成的基因编码框的方法流程图；

图2为本发明实施例提供的鉴定SL序列反式剪切所形成的基因编码框的方法原理图；

图3为本发明实施例提供的秀丽隐杆线虫slORF的示意图；

图4为本发明实施例提供的布氏锥虫slORF的示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本文中提及“实施例”意味着，结合实施例描述的特定特征、结构或特性可以包含在本申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是指相同的实施例，也不是与其它实施例互斥的独立的或备选的实施例。本领域技术人员显式地和隐式地理解的是，本文所描述的实施例可以与其它实施例相结合。

本申请的说明书和权利要求书及所述附图中的术语“第一”、“第二”、“第三”和“第四”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤、过程、方法等没有限定于已列出的步骤，而是可选地还包括没有列出的步骤，或可选地还包括对于这些过程、方法、产品或设备固有的其它步骤元。

本发明的目的是提供一种一种鉴定SL序列反式剪切所形成的基因编码框的方法，本发明解决了现有技术存在SL序列的添加将有可能引入新的翻译起始位点，从而改变原有的基因编码框，形成新的蛋白质序列，使得原始的基因组注释信息失效，从而引发一系列的错误分析和结论的问题。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1所示，本发明提供了一种鉴定SL序列反式剪切所形成的基因编码框的方法，包括：

步骤100：获得全部的基因组所有转录本和已知编码框序列；

步骤200：根据所述转录本和已知编码框序列获取所述基因组的比对文件；

步骤300：基于转录组测序技术，根据所述比对文件筛选局部比对的序列，并提取未成功比对上的序列；

步骤400：对所述未成功比对上的序列与指定的SL序列进行比较，若所提取的序列为指定SL序列的末端至少8个碱基序列，则确定该序列为SL剪切位点所产生的序列；

步骤500：根据所述比对文件中所记录的位置信息推算SL序列的剪切位点，并将所述剪切位点转换成转录组坐标；

步骤600：根据所述转录组坐标和所述转录本得到完整的mRNA序列；

步骤700：利用核糖体印迹测序数据，根据所述完整的mRNA序列获得完整的mRNA对应的编码框。

进一步的，所述位置信息包括：

第3列的比对染色体信息和第4列的比对坐标。

进一步的，所述根据所述转录组坐标和所述转录本得到完整的mRNA序列包括：

将SL剪切位点以前的序列删除，并替换成指定的SL序列以得到完整的mRNA序列。

进一步的，根据基因组中所有编码框的密码子使用情况，计算每个密码子在整体基因组编码序列中的出现频率，然后计算每个已知编码框中密码子频率的平均值，并将其转化为Z-score。

进一步的，根据所述完整的mRNA序列获得完整的mRNA对应的编码框包括：

基于核糖体印迹测序技术，获取核糖体印迹数据序列；

对所述核糖体印迹数据序列筛选，获得满足条件的核糖体印迹数据序列；

根据所述满足条件的核糖体印迹数据序列获得5’端与P-site之间不同距离的出现频率；

根据每一个碱基位于P-site的概率和所述已知编码框中密码子频率的平均值，对编码框进行预测，得到预测编码框；

对所述预测编码框进行筛选，满足第一条件的预测编码框，进行输出，得到编码框。

进一步的，所述每一个碱基位于P-site的概率的包括：

5’端与P-site之间不同距离的出现频率和所述满足条件的核糖体印迹数据序列比对位置。

本实施例还公开了鉴定SL序列反式剪切所形成的基因编码框的方法具体步骤：

(1)获取基因组中mRNA和已知编码框的序列

通过读取基因组注释文件信息，提取基因转录本位置信息，并获得其序列信息。并根据注释文件中基因编码区的坐标信息，提取已知编码框的序列和位置信息，进而获得全基因组所有转录本和已知编码框序列。

(2)已知编码框中密码子使用频率特征的提取

根据基因组中所有编码框的密码子使用情况，计算每个密码子在整体基因组编码序列中的出现频率，然后计算每个已知编码框中密码子频率的平均值，并将其转化为Z-score。

(3)确定SL剪切位点

读取RNA-Seq转录组测序技术、与基因组的比对文件(sam或bam格式)，根据比对文件第6列的详细比对信息，筛选局部比对的序列(标识符‘S’)，提取未成功比对上的序列，并与指定的SL序列进行比较，如果所提取的序列为指定SL序列的末端8个(或以上)碱基序列，则确定该序列为SL剪切位点所产生的序列。根据比对文件中所记录的位置信息(第3列：比对染色体信息；第4列：比对坐标)，推算SL序列的剪切位点，并将其转换成转录组坐标。

(1)获取成熟mRNA序列信息

根据步骤(3)中所获得的SL剪切位位点转录组坐标信息，结合步骤(1)中所获得的每一个转录本序列信息，将SL剪切位点以前(5’端方向)的序列删除后，替换成指定的SL序列，从而形成完整的、成熟的mRNA序列。

(2)核糖体印迹数据(RPF)质量评估

读取Ribo-Seq数据与基因组序列的比对文件(sam或bam格式)，根据所记录的位置信息(第3列：比对染色体信息；第4列：比对坐标)提取每一条序列在基因组上的比对位置坐标，并将其转换为转录组坐标。使用multitaper算法对不同长度RPF序列的3碱基周期性，对完全不具有周期性的噪音数据进行过滤。各长度的3碱基周期性通过进行评估，频率显示为3.33或0.34Hz，P值≤0.01的长度得以保留，用于后续分析。

(3)核糖体印迹数据特征训练

通过提取比对到已知编码框启始或终止密码子的RPF比对信息，计算每条RPF 5’端与启始密码子(P-site)或终止密码子(A-site)的距离，统计各个长度的RPF的5’端与P-site之间不同距离的出现频率。

(4)核糖体印迹数据与密码子频率在编码框预测中的权重分配

根据各RPF在相位0，1和2位置出现的频率计算其分布集中度。分布集中度由复杂度(熵)进行描述，公式如下：

其中，i表示不同的相位(0，1和2)，P

(5)候选编码框的选取与搜索

依据(4)中所有包含有SL序列的成熟转录本的序列信息，提取所有候选编码框序列，依据标准为，拥有启始密码子(NUG)和终止密码子(UAG，UAA，UGA)并且其长度为3的倍数。优先搜索AUG起始的候选编码框，由长到短，逐一进行计算，AUG起始的候选编码框全部搜索完全且不满足输出条件后，再进行NUG编码框的搜索和计算，搜索结果输出，作为候选编码框用于下一步处理。

(6)编码框的预测

根据Ribo-Seq数据中各RPF的比对位置以及第(6)步中获得的其5’端与P-site之间的距离信息，计算各转录本上每一个碱基正好位于P-site的概率，并将其转化为Z-score。

结合每个碱基位于P-site的概率以及步骤(2)中计算所得的其作为密码子使用的概率，进行四组统计检验，分别是(a)位于相位0上的Z-score值极显著大于(单尾检验)位于相位1上的Z-score；(b)位于相位0上的Z-score值极显著大于(单尾检验)位于相位2上的Z-score；(c)位于相位0上的密码子的使用频率值极显著大于(单尾检验)位于相位1上的密码子频率；(d)位于相位0上的密码子的使用频率值极显著大于(单尾检验)位于相位2上的密码子频率。

以上统计所得的4个P值经加权卡平方算法(Weighted chi-square method)合并成最终P值，计算方法如下，

首先按照步骤(7)中分配的权重，将P值转化为卡平方值，公式如下：

其中M表示合并后的卡方值，i为第i个检验，Pi为第i个检验的P值，wi为第i个P值的权重，由于wi之和须为1，且RPF和密码子使用频率各进行了两次检验，因此，相对应P值的权重为上一步中计算所得RPF/密码频率权重的一半。

计算自由度(k)

k＝2{E(M)}

其中，

s

其中，w

其中，

0.75ρ

最后可以求解ρ的近似值为-2.167+(10.028-4q

根据计算所得的自由度k和合并后的卡方值，根据卡方分布2χ

(1)编码框输出错误发现率(FDR)控制

输出P值≤0.001的修选编码框并根据Benjamini和Hochberg法控制FDR≤0.0001，满足这一标准的候选编码框进行最后的结果输出。

本实施还公开了关于秀丽隐杆线虫(Caenorhabditis elegans)slORF的鉴定：

从NCBI数据库(https://www.ncbi.nlm.nih.gov/)下载所需的实验数据(编号：PRJNA208993)，该数据由Stadler M等于2013年发表于PLoS Genet，文章名称为“Conservedtranslatome remodeling in nematode species executing ashared developmentaltransition”，共包含48个样本的测序数据，其中包括6个转录组测序数据和6个核糖体印迹测序数据。从以下地址中下载本实施例所选用的秀丽隐杆线虫(C.elegans)参考基因组下载地址为http://ftp.ensemblgenomes.org/pub/metazoa/release-54/fasta/caenorhabditis_elegan s/dna/Caenorhabditis_elegans.WBcel235.dna.toplevel.fa.gz，同时下载基因组注释文件(http://ftp.ensemblgenomes.org/pub/metazoa/release-54/gtf/caenorhabditis_elegan s/Caenorhabditis_elegansWBcel235.54.gtf.gz)，SL序列为GGUUUAAUUACCCAAGUUUGAG。利用本发明方法，首先对秀丽隐杆线虫的转录组测序数据和核糖体印迹测序数据进行预处理，调用STAR软件将其与基因组序列进行比对，比对结果输出为bam格式，随后将转录组数据和核糖体印迹数据比对结果，以及秀丽线虫全基因组序列文件、注释文件和上述SL序列输入到本发明中开发的软件slORFfinder中，进行slORF预测。结果显示，在此数据集中，成功鉴定得到了334个slORF，图2所示为其中一个例子。

本实施例还公开了布氏锥虫(Trypanosomabrucei)slORF的鉴定：

从NCBI数据库(https://www.ncbi.nlm.nih.gov/)下载所需的测序数据(编号：PRJNA246300),该数据由Jensen BC等于2014年发表于BMC Genomics，文章名称为“Extensive stage-regulation of translation revealed by ribosome profilingofTrypanosomabrucei”。该数据集中共包含22个样本数据，我们使用了其中9个转录组测序数据和9个核糖体印迹测序数据。下载本实施例所需的布氏锥虫(T.brucei)参考基因组序列(https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/445/GCF_000002445.2_ASM244v1/GCF_000002445.2_ASM244v1_genomic.fna.gz)和基因组注释文件(https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/445/GCF_000002445.2_ASM244v1/GCF_000002445.2_ASM244v1_genomic.gtf.gz)，该实施例中SL序列为AACUAACGCUAUUAUUGAUACAGUUUCUGUACUAUAUUG。我们首先对布氏锥虫的转录组测序数据和核糖体印迹测序数据进行预处理，调用STAR软件将其与基因组序列进行比对，比对结果输出为bam格式，根据本发明方法的描述，随后将转录组数据和核糖体印迹数据比对结果，以及布氏锥虫全基因组序列文件、注释文件和上述SL序列输入到本发明中开发的软件slORFfinder中，进行锥虫slORF预测。结果显示，在此数据集中，成功鉴定得到了586个slORF，图3所示为其中一个例子。

本发明的有益效果如下：

本发明通过分析这类生物中转录组测序(RNA-Seq)和核糖体印迹测序(Ribo-Seq)数据，搜索SL序列在各基因mRNA上剪切位点，进而推断基因在SL反式剪切以后所形成的成熟mRNA序列，借助核糖体印迹测序数据在成熟mRNA上的分布特征，确定该mRNA的翻译相位，从而预测新的编码框。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。