一种m6A测序的peak功能元件注释的方法、系统、设备和介质技术方案

本发明专利技术公开了一种m6A测序的peak功能元件注释的方法、系统、设备和介质，属于m6A测序技术领域。所述方法包括以下步骤：获得目标物种的m6A测序结果，并利用call

【技术实现步骤摘要】
一种m6A测序的peak功能元件注释的方法、系统、设备和介质


[0001]本专利技术属于m6A测序
，具体地，涉及一种m6A测序的peak功能元件注释的方法、系统、设备和介质。

技术介绍

[0002]近年来，RNA的各种修饰及其功能是表观遗传学研究的热点之一。RNA修饰有很多类型，不同类型的RNA修饰含量和功能也存在很大差异。在各种类型的RNA化学修饰中，甲基化是最主要的修饰形式。7
‑
甲基鸟嘌呤(N7
‑
methylguanosine，m7G)，6
‑
甲基腺嘌呤(N6
‑
methyladenosine，m6A)和5
‑
甲基胞嘧啶(5
‑
methylcytosine，m5C)是最具代表性的几种修饰，又以信使RNA(messenger RNA，mRNA)内部修饰丰度最高的m6A研究最为深入。m6A对RNA的加工代谢及多种生物学过程有着重要的调控作用，其研究成为近几年RNA研究领域的重点和热点。
[0003]RNA修饰研究的重点问题之一是通过RNA修饰检测技术，准确揭示RNA修饰位点的分布和动态调控。m6A RNA免疫共沉淀测序(methylated RNA immunoprecipitation sequencing，MeRIP
‑
seq)技术是最早开发和应用最广的m6A测序技术，该方法结合了MeRIP实验可以广泛且特异的获得m6A修饰的优点和高通量测序(High
‑/>throughput sequencing)通量高、准确度高的优点。MeRIP
‑
seq主要包含MeRIP实验和高通量测序两部分。MeRIP实验中，RNA先被片段化成100～200nt(nucleartide)的片段，再分成两份。第一份片段化后的RNA与带有m6A抗体的磁珠共孵育，筛选出带有m6A修饰的RNA片段，经过建库获得immunoprecipitation(IP)文库。第二份片段化后的RNA不进行抗体孵育步骤，按照一般转录组的建库方式直接建立input文库。IP文库和input文库经过高通量测序，获得原始数据，进入数据分析的过程。原始数据经过质控，比对参考基因组之后，获得reads在基因组各个区域的分布情况。数据分析的关键步骤之一是m6A修饰区域的检测，也称为call
‑
peak。call
‑
peak步骤是以input文库为对照，对IP文库的reads分布结果进行矫正，以去除掉一些由于扩增偏好性等可能的原因造成的假阳性的reads富集区域，尽可能准确地获得真实m6A修饰显著富集的区域，而这些显著富集区域，就是我们理解的“peak”(富集峰)。目前已经有多种call
‑
peak软件被广泛应用，例如exomePeak，MACS2，MeTPeak，RADAR等。call
‑
peak软件的输出结果主要包括peak所在的染色体序号，染色体上的起始位点，染色体上的终止位点，相对基因组的正负链信息，富集倍数，富集显著性等结果，部分软件也会给出peak所在的基因信息。
[0004]由于真核生物中普遍存在可变剪切，一个peak所在的基因可能包含多条转录本，他们的位置和功能也可能存在差异，因此判定一个peak可能来源于哪条转录本是非常重要的。其次，一条成熟的RNA有多个功能区域，包括5'非翻译区(5'UTR)，起始密码子(start codon)，编码区域(CDS)，终止密码子(stop codon)，3'非翻译区(3'UTR)，判断一个peak位于某个功能区域对研究其生物学功能非常有价值。仅有call
‑
peak软件提供的原始信息，研究者不足以对测序结果进行深入的研究，因此需要对每个peak进行更详细的功能元件区域
注释，以及对于各功能元件区域的m6A修饰富集程度进行计算和表示，才能更近一步进行相应的m6A修饰所相关的功能研究。
[0005]目前，现有的区域注释软件大多数都是基于DNA测序而设计，DNA具有双链，这类型的注释软件一般根据区域的位置信息(染色体：起始位置
‑
终止位置)来定位到基因组上与该区域有交集的区间，而不考虑正义链和反义链的信息。而针对RNA类型的测序来说，尤其是链特异性建库相关的测序结果，得到正义链和反义链的信息对于基因/转录本的定位就显得尤为重要，因此产生的对下游研究的重要性也就不言而喻。另外，对于测序结果中同一个待注释区间，根据基因组的结构注释文件信息，其可能得到多个功能元件注释结果，在这种情况下注释信息的选择和注释优先级也是一个很重要的问题。同时，目前现有的注释工具或软件，对于功能元件类型的划分也不尽相同，很多是基于DNA相关技术所关注的类型来重点进行注释，不适用于m6A研究所重点关注的类型。m6A测序peak注释后所得到的功能元件占比，在现有软件，例如ChIPSeeker等，其只考虑了每一个注释到的功能元件的绝对占比，而m6A测序中，由于各功能元件实际区间长度的不尽相同，例如，一个转录本，一般CDS区域范围远大于UTR区域，仅以落在相对功能元件区域的peak数量占比来评估在各个功能元件上的富集程度就显得有失偏颇。

技术实现思路

[0006]为了解决上述技术问题中的至少一个，本专利技术采用以下技术方案：
[0007]本专利技术第一方面提供一种m6A测序的peak功能元件注释方法，包括以下步骤：
[0008]S1，获得目标物种的m6A测序结果，并利用call
‑
peak软件得到peak信息；
[0009]S2，获得目标物种的基因组，并获得蛋白编码类型的转录本注释信息，包括基因信息、各功能元件的区间位置信息和相对于基因组的正负链信息；
[0010]S3，将步骤S1获得的peak信息，与步骤S2获得的转录本注释信息进行比对，获得每一个peak的功能元件注释结果，具体注释方法为，若某个peak与某个功能元件之间满足：
[0011](1)所述peak所在的染色体与所述功能元件所在的染色体相同；
[0012](2)所述peak所对应的基因与所述功能元件所在的基因相同；
[0013](3)所述peak相对于基因组的正负链信息与所述功能元件相对于基因组的正负链信息相同，
[0014](4)所述peak所在的染色体上的起始位置和终止位置与所述功能元件区间有交集，
[0015]则将所述功能元件的类型作为所述peak的注释。
[0016]在本专利技术的一些实施方案中，所述功能元件包括3'
‑
UTR、startcodon、CDS、stopcodon和5'
‑
UTR。
[0017]在本专利技术的一些实施方案中，所述start codon和所述stop codon的区间范围分别为转录起始位点和转录终止位点的上下游100nt区间，即start codon和stop本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种m6A测序的peak功能元件注释方法，其特征在于，包括以下步骤：S1，获得目标物种的m6A测序结果，并利用call
‑
peak软件得到peak信息；S2，获得目标物种的基因组，并获得蛋白编码类型的转录本注释信息，包括基因信息、各功能元件的区间位置信息和相对于基因组的正负链信息；S3，将步骤S1获得的peak信息，与步骤S2获得的转录本注释信息进行比对，获得每一个peak的功能元件注释结果，具体注释方法为，若某个peak与某个功能元件之间满足：(1)所述peak所在的染色体与所述功能元件所在的染色体相同；(2)所述peak所对应的基因与所述功能元件所在的基因相同；(3)所述peak相对于基因组的正负链信息与所述功能元件相对于基因组的正负链信息相同，(4)所述peak所在的染色体上的起始位置和终止位置与所述功能元件区间有交集，则将所述功能元件的类型作为所述peak的注释。2.根据权利要求1所述的m6A测序的peak功能元件注释方法，其特征在于，所述功能元件包括3'
‑
UTR、startcodon、CDS、stopcodon和5'
‑
UTR。3.根据权利要求2所述的m6A测序的peak功能元件注释方法，其特征在于，所述start codon和所述stop codon的区间范围分别为转录起始位点和转录终止位点的上下游100nt区间。4.根据权利要求1所述的m6A测序的peak功能元件注释方法，其特征在于，如果一个peak注释到同一转录本的多个功能元件，那么根据以下顺序设定功能元件的优先级：Stop codon、3'
‑
UTR、Start codon、5'
‑
UTR、CDS。5.根据权利要求1所述的m6A测序的peak功能元件注释方法，其特征在于，进一步包括：S4，分别计算注释到每一种功能元件的peak数量占比P，以及步骤S2获得的所有功能元件的区...

【专利技术属性】
技术研发人员：陈丽，梁思源，韩斐然，
申请(专利权)人：杭州联川生物技术股份有限公司，
类型：发明
国别省市：