附图说明" src="/d/file/p/2020/11-23/ad46c4f7ec3386b613049dc53debed3f.gif" />
RNA 5mC重亚硫酸盐测序的文库构建方法及其应用
技术领域t
本发明属于基因组测序技术领域，具体的，涉及一种通过利用只含ACT三碱基的随机六聚引物对重亚硫酸盐处理的RNA片段进行逆转录后的文库构建方法及其在5mC高通量测序中的应用。t
背景技术t
RNA同时具有调控和信息分子的双重功能，在众多细胞机制中发挥着核心作用。RNA转录后的修饰为RNA功能的多样化奠定了化学基础。自然界中的RNA修饰广泛存在于A、U、C、G四类核苷酸上。截止目前，RNA修饰数据库RNAMDB共收录了109种RNA的修饰形式，其中甲基化修饰约占RNA修饰总量的80％，这类甲基化修饰主要发生在碱基基团的氮原子N上，以及嘌呤和嘧啶的C原子上，或2’-OH氧原子等特殊位置上。RNA核苷上的甲基化反应主要依赖于生物体内的甲基转移酶及甲基基团供体完成。6-甲基腺嘌呤(N6-methyladenosine,m6A)，是发生在碱基A第六位N原子上的甲基化修饰，作为真核生物中最常见的一种RNA转录后修饰，因其丰富的含量和高度的保守型，近年来得到了广泛的关注和研究。t
除m6A外，5-甲基胞嘧啶(5-methylcytosine，5mC)是另一种在RNA中广泛存在的甲基化修饰形式。虽然5mC修饰在DNA研究中已被广泛报道，而在RNA中的研究尚处于起步阶段。早在20世纪70年代5mC就已被发现存在于仓鼠细胞的mRNA中，其修饰位点的分布特征以及生物学功能等仍不清楚，直到近期对于RNA5mC的基因组水平上的单位点检测及生物学功能的研究才有了初步探索。t
和m6A一样，RNA5mC修饰也应该是动态可逆的，甲基转移酶以S-腺苷甲硫氨酸(SAM)作为甲基供体，将甲基转移到胞嘧啶C形成5-甲基胞嘧啶(5mC)。RsmB是第一个被发现的RNA5mC甲基转移酶，主要催化细菌rRNA上的甲基化形成。随后30多种RNA上的5mC甲基转移酶陆续被发现，这些甲基转移酶主要可分为NOP2/NOL1,YebU/Trm4,RsmB/Yn1022c和PH1991/NSUN四类，而且这些酶在真核生物中有很高的保守性，近年来NSUN蛋白家族被广泛深入研究。人的NSUN蛋白家族共有9个蛋白，该家族的多个成员都具有潜在的5mC甲基转移酶功能结构域，其中NSUN2的催化活性已被证实。除NSUN2外，DNMT2也被认为可能是哺乳动物的RNA5mC甲基转移酶，它和NSUN2的催化位点有一定的交集。其它的酶如NSUN1和NSUN3-7被预测可以催化一些保守的甲基化位点的甲基化过程。NSUN1、NSUN2、NSUN5已被证实能结合RNA，但这些酶的特异性结合底物并不是很清楚。t
利用5mC抗体免疫沉淀结合亚硫酸氢盐处理测序的结果发现了古生菌mRNA中的多个5mC修饰，其保守序列为AU(m5C)GANGU，和古生菌的rRNA上的保守序列一致，说明mRNA和rRNA上的5mC修饰可能通过同一种甲基转移酶催化形成。t
DNA中5mC的去甲基化反应由TET家族蛋白催化介导，而近期一项研究表明TET蛋白也能介导RNA5mC去甲基化形成5hmC，而且相较于TET1对DNA的高度活性，TET3对RNA的选择性似乎更强。针对RNA5mC去甲基化酶的研究还有待更深入的探索。t
相比m6A，5mC在mRNA中的功能研究还并不深入。基于现有研究推测，5mC的修饰可能会调控mRNA的选择性剪接，5mC修饰水平可能会影响外显子的保留水平以及转录本的组装形式。除此之外，5mC还可能与蛋白质翻译以及mRNA的稳定性有较为密切的关系。总之，5mC修饰在RNA中广泛存在而且其功能可能涉及到包括细胞内信号转导，组织发育分化和癌症等许多方面，而对mRNA中5mC位点的检测和分布规律的探索，5mC修饰酶及结合蛋白的发现将有助于阐述5mC修饰对mRNA加工代谢的调控机制。t
RNA整体水平上的5mC修饰可用质谱技术检测，而对于5mC的单位点鉴定迄今为止已有四种方法被报道，分别为5mC重亚硫酸盐测序(Bisulfitesequencing,BS-seq)、5mC-RIP、Aza-IP、miCLIP。其中RNA重亚硫酸盐测序是目前为止最为理想的方法，重亚硫酸盐处理能够将核苷酸序列中未发生甲基化的胞嘧啶C(占绝大部分)脱氨基转换成尿嘧啶U，而甲基化的胞嘧啶则保持不变，进行PCR扩增后尿嘧啶U全部转换成胸腺嘧啶T，因此可与原本具有甲基化修饰的C碱基区分开来，该方法可以检测到RNA经重亚硫酸盐处理后转换为尿嘧啶的非甲基化胞嘧啶，最后，对PCR产物进行测序，并且与未经处理的序列进行比较，可判断胞嘧啶位点是否发生了甲基化。已有研究证实了tRNA和rRNA中均存在5mC修饰，早期研究由于受到检测技术的限制，在mRNA中仅发现少数的5mC修饰，而近期通过重亚硫酸盐测序方法发现在HeLa细胞的mRNA和非编码RNA中5mC均存在着广泛的修饰，并且指出mRNA中的5mC修饰主要富集在非编码区(UTR)。t
发明内容t
本发明人发现已报道的RNA5mC重亚硫酸盐测序在文库构建的逆转录阶段采用含有ATCG四个碱基的随机六聚引物，其存在一定的局限性。RNA中甲基化C在C碱基中所占比例很低，利用超高效液相色谱串联质谱检测到人细胞mRNA中的5mC水平在所有C碱基中的占比不足千分之一(如图1所示，在HeLa细胞中mRNA的5mC的水平是0.03918％，293T细胞中为0.072401％)。因此，RNA经重亚硫酸盐处理后，绝大部分(例如99.7％以上，甚至99.9％以上)的非甲基化C都将转换成为U，导致RNA序列中A、U、G的富集以及C的显著减少，按照四个碱基ATCG的随机排列，有(4^6)4096种组合，其中包含G的组合有(4^6-3^6)3367种，即超过80％的序列不能被ATCG四碱基随机六聚引物有效识别，造成了扩增过程中包含G的ATCG四碱基随机六聚引物无法有效配对原来的非甲基化C区域，导致随机六聚引物逆转录效率降低和扩增的偏好性。其次，现有的RNA5mC重亚硫酸盐测序采用了基于“双碱基编码原理”的SOLID测序平台，与Illumina平台直观的碱基序列不同，SOLID测序将reads利用颜色空间进行编码，将每一个碱基与它邻近的碱基用一种颜色表示,但在荧光解码阶段，鉴于其是双碱基确定一个荧光信号，因而一旦发生错误就容易产生连锁的解码错误，并不适用于5mC单碱基转化水平的重亚硫酸盐测序。因此，现有的测序结果无法反映出RNA5mC修饰的真实的分布规律。t
本发明人发现，利用只含有ACT的三碱基随机六聚引物对经重亚硫酸盐处理的RNA进行逆转录PCR，三碱基随机六聚引物中的A能够高效匹配原来的非甲基化C区域(重亚硫酸盐处理之后为U)，克服了原来含有G的四个碱基ATCG的随机六聚引物无法有效配对由非甲基化C转变而来的U区域的缺陷，可以提高经重亚硫酸盐转化的RNA与逆转录引物的结合效率，其逆转录效率显著高于传统的ACTG四碱基随机六聚引物，从而有利于5mC位点的检测以及对5mC的分布规律和功能机制的进一步分析。其原理如图2所示，从总RNA中分离纯化mRNA或者前体mRNA(pre-mRNA)，经片段化处理及重亚硫酸盐转化后，绝大部分C转变为U，ACTG四碱基随机六聚引物由于其中存在G，仅能与少量重亚硫酸盐转化后的片段匹配，因此逆转录效率较低，仅能逆转录扩增部分片段，无法全面覆盖RNA(或相应的cDNA)的全部序列，而利用ACT三碱基随机六聚引物对已转化的C(即非甲基化的C)的匹配效率要远远高于普通的ACTG四碱基随机六聚引物，经II链合成和PCR扩增后，能有效扩增绝大多数含非甲基化C的RNA片段，由此构建的cDNA文库基本能完全覆盖测序样品的所有序列，能得到更优的更真实的数据。t
在此基础上，本发明涉及5mC甲基化RNA的重亚硫酸盐测序文库构建、高通量测序及甲基化检测方法，包括以下步骤：t
(1)将RNA样品片段化；t
通过物理或化学方法将RNA样品打断成为适合于测序的长度的片段。优选将RNA样品打断为约100nt大小的片段。可以使用商购的试剂盒，例如AmbionRNA片段化试剂盒实现RNA样品的片段化。t
“RNA样品”是分离纯化的待测RNA，可以来自于人体、动物体、植物体或其器官、组织或细胞。所述RNA可以是mRNA、tRNA、rRNA、总RNA或其它包含或可能包含5mC甲基化位点的RNA。例如可以使用Trizol法从组织或细胞中提取总RNA。还可以使用商购的试剂盒提取(分离纯化)RNA，例如使用Ambion的mRNA提取试剂盒提取mRNA。t
为了鉴定重亚硫酸盐处理RNA方法的可靠性，RNA样品中可以包含一定比例的非甲基化RNA标准品作为鉴定重亚硫酸盐转化效率的标准品，例如RNA样品可以是非甲基化RNA标准品与待测RNA以1:100(重量比)的比例混合的混合物。非甲基化RNA标准品序列不含5mC甲基化位点，因此在用重亚硫酸盐处理时全部胞嘧啶理论上都将转换为尿嘧啶，在重亚硫酸盐处理后对RNA样品中非甲基化RNA标准品序列中胞嘧啶转化为尿嘧啶的情况(例如百分比)进行分析，可以反映重亚硫酸盐转化效率。本发明中可以使用小鼠DHFR基因的RNA序列作为非甲基化RNA标准品。DHFR的RNA序列可以用T7高效RNA合成试剂盒(NEB、E20405)从携带小鼠DHFR基因全长的质粒pcDNA3-HA-mDHFR上转录出。t
(2)对片段化的RNA样品进行重亚硫酸盐处理；t
将片段化的RNA样品用重亚硫酸盐例如重亚硫酸钠进行处理，以使非甲基化胞嘧啶转换为尿嘧啶。用重亚硫酸盐处理RNA的方法是本领域公知的。重亚硫酸盐处理可以包括重亚硫酸盐温育、脱盐和脱磺酸化步骤。重亚硫酸盐处理的示例性方法包括用新鲜配制的重亚硫酸盐溶液(40％重亚硫酸钠，600μM氢醌,pH5.1)溶解mRNA沉淀，置于75℃PCR仪上处理4小时；然后对重亚硫酸盐处理后的样品进行脱盐处理，例如使用MicroBio-spin6chromatographycolumns(Bio-Rad)；随后加入1MTris(pH9.0)于75℃PCR仪上1小时进行脱磺酸处理。最后加入糖原和3M醋酸钠(pH5.2)及预冷的纯乙醇，放置-80℃过夜沉淀mRNA。t
(3)利用ACT三碱基随机六聚引物对重亚硫酸盐处理后的RNA样品进行逆转录，合成cDNA；t
为了提高经重亚硫酸盐处理的RNA的逆转录效率，以重亚硫酸盐处理后的RNA样品为模板，使用ACT三碱基随机六聚引物和反转录酶合成cDNA。例如使用SuperscriptIIReverseTranscriptaseKit(Invitrogen)进行逆转录合成cDNA，并用自行合成的ACT三碱基随机六聚引物代替试剂盒中的AGCT四碱基随机六聚引物。t
“ACT三碱基随机六聚引物”是指由A、C、T三种碱基随机合成的六聚引物的混合物，该混合物中任一引物均由且仅由A、C、T组成。t
(4)用逆转录合成的cDNA构建测序文库。t
用逆转录合成的cDNA进行II链合成，纯化产物后，在cDNA片段的末端加上“A”碱基(Atailing)，再连接接头(Adapter)并纯化回收，以纯化回收的DNA片段为模板进行PCR扩增，并纯化回收PCR产物，获得RNA5mC重亚硫酸盐测序文库。t
可以由逆转录合成的cDNA开始，使用KAPAStrandedmRNA-SeqKit(Cat.No.KK8420)构建测序文库，具体地，可以按照试剂盒的操作说明进行II链合成及之后的建库(即构建测序文库)步骤。t
(4)进行高通量测序。t
利用上述构建的测序文库进行高通量测序。例如在二代测序平台例如Illumina的Hiseq测序平台上采用边合成边测序的方法进行序列测定。所述Illumina的Hiseq测序平台包括例如Hiseq2000/2500/3000/4000测序平台。t
(5)对测序结果进行数据分析，获得RNA的5mC甲基化信息。t
重亚硫酸盐处理会将非甲基化的胞嘧啶C转变为U，而甲基化的胞嘧啶保持不变，然后PCR扩增使得U变成T。因此，在将测序得到的序列与参考基因组进行比对时，比对结果为C－C(参考基因组上在某个位置上是C，测得的reads在该位置上也是C)的是甲基化的胞嘧啶，比对结果为C－T的是非甲基化的胞嘧啶。具体的分析过程包括下面三个步骤：t
质量控制：对测序得到的序列(reads)进行碱基质量控制、修剪引物等处理。例如可以通过FASTX-Toolkit，Cutadapt，Trimmomatic等软件[1-3]进行测序数据的质量控制。t
序列比对：将质量控制后的测序序列与人参考基因组上的原始序列进行比对。在一些情况下，考虑到转录组测序数据会包含一些跨内含子的序列，因此，对于那些在与人参考基因组原始序列的比对中未比对上的测序序列，再与人参考基因组的转录组数据以及外显子-外显子连接构成的连接序列进行比对。可以使用hg19版本的人参考基因组进行比对。比对可以通过Bismark软件[4]进行。t
甲基化水平分析：对于人参考基因组原始序列中的每个胞嘧啶位置，确定测序序列中在该位置表现为甲基化的序列数(即比对结果为C－C，以下称甲基化序列数)和表现为非甲基化的序列数(即比对结果为C－T，以下称非甲基化序列数)。对于与人参考基因组的转录组数据以及外显子-外显子连接序列比对的结果，与上述方法相同确定每个胞嘧啶位置的甲基化序列数和非甲基化序列数，并将其对应的位置信息转换为人参考基因组原始序列位置信息，并综合起来计算甲基化水平。使用公式M/(U+M)计算每个胞嘧啶位置的甲基化水平，其中U和M分别是在选定胞嘧啶位置上的非甲基化序列数和甲基化序列数。在一些情况下，对于比对上的测序序列，如果在一个测序片段(即测序文库构建过程中打断形成的片段)中所包含的胞嘧啶位置中有超过30％的位置测序结果为C，鉴于RNA甲基化的低水平，如此高比例的甲基化胞嘧啶位点很可能为处理不充分导致，因此在计算甲基化水平之前先将这些测序片段过滤掉。t
本发明中使用只含有ACT三碱基的随机六聚引物对重亚硫酸盐处理的RNA片段进行逆转录，能够高效匹配5mC修饰区域，可以大幅提高逆转录和PCR扩增的效率，更有利于逆转录时与随机引物配对区域潜在的5mC位点的检测。t
进一步，本发明中使用不含有5mC甲基化位点的DHFR基因的RNA作为标准品，可以鉴定重亚硫酸盐转化效率。t
进一步，本发明的重亚硫酸盐处理方法，可以达到99.7％的转化效率。t
进一步，本发明中使用Hiseq2000测序平台进行高通量测序，能够高效、准确地用于5mC重亚硫酸盐测序，规避了SOLID测序平台错误率高的弱势。t
利用本发明的方法检测RNA的5mC甲基化位点，可以提高逆转录合成效率，检测更加全面和完整。利用本发明方法构建的cDNA文库经由Hiseq2000测序平台所得的测序数据中，检测到了14691个5mC修饰位点，且大多倾向于分布在mRNA的编码区，尤其是翻译起始位点附近，进一步研究结果表明5mC还分布在剪切位点附近的外显子中，且5mC分布存在一定的序列偏好性，相对于CpG和CHG，5mC甲基化位点主要分布在CHH区域。t
附图说明t
图1显示了利用超高效液相色谱(UHPLC)检测的两种人细胞系HeLa和293T细胞的mRNA中5mC占所有胞嘧啶总数的比例。t
图2显示了本发明RNA5mC重亚硫酸盐测序的原理示意图。t
图3显示了实施例1中对于DHFRRNA标准品的逆转录效率的鉴定结果。A图显示经重亚硫酸盐处理前后的mRNA样品利用不同的逆转录引物进行逆转录合成，并利用DHFR特异性引物进行PCR扩增的产物的琼脂糖凝胶电泳检测结果。B图显示对琼脂糖凝胶电泳检测结果的定量分析。t
图4显示了5mC修饰的类型分布。A图是显示所有含有5mC修饰基因的基因类型的饼图，包括编码蛋白的mRNA和非编码RNA。B图是显示含有5mC修饰的各种非编码基因的个数的柱形图。t
图5显示了Hela细胞系中5mC在mRNA上各区间的分布特征。A图显示了Hela细胞中5mC在不同区间的比例分布。B图显示了实际和预期的甲基化位点在各个区间的分布比例。t
图6是5mC在mRNA5’UTR、CDS以及3’UTR区间的分布曲线。t
图7显示了5mC修饰附近的序列特征。A图显示了5mC位点及下游碱基CG,CHG,CHH的比例，其中H代表A、C、U。B图显示了5mC附近5个碱基的序列组成。t
图8显示了利用重亚硫酸盐处理结合Sanger测序方法对实施例1中利用ACT随机六聚引物测序所得数据中的3个基因PLOD3、COL4A5和FAM129B的5mC位点的验证。左图为重亚硫酸盐未处理和处理的Sanger测序峰图和序列比对图；右图为高通量测序所得对应基因的甲基化水平，其中，空心圆圈表示非甲基化即转化了的胞嘧啶，实心圆圈表示甲基化即未转化的胞嘧啶。t
具体实施方式t
可以理解的是，在此描述的特定实施方式通过举例的方式来表示，其并不作为对本发明的限制。在不偏离于本发明范围的情况下，本发明的主要特征可以用于各种实施方式。本领域的技术人员将会意识到或能够确认，仅仅使用常规实验，许多等同物都能应用于本文所描述的特定步骤中。这些等同物被认为处在本发明的范围之内，并且被权利要求所覆盖。t
本发明一方面提供了一种提高重亚硫酸盐处理后的RNA样品的逆转录效率的方法，其特征在于：t
利用ACT三碱基随机六聚引物对重亚硫酸盐处理后的RNA样品进行逆转录，合成cDNA。t
本发明另一方面提供了RNA5mC重亚硫酸盐测序文库构建方法，其特征在于包括以下步骤：t
(1)将RNA样品片段化；t
(2)对片段化的RNA样品进行重亚硫酸盐处理；t
(3)利用ACT三碱基随机六聚引物对重亚硫酸盐处理后的RNA样品进行逆转录，合成cDNA；t
(4)用逆转录合成的cDNA构建测序文库。t
本发明另一方面提供了5mC甲基化RNA的测序方法，其特征在于包括以下步骤：t
(1)将RNA样品片段化；t
(2)对片段化的RNA样品进行重亚硫酸盐处理；t
(3)利用ACT三碱基随机六聚引物对重亚硫酸盐处理后的RNA样品进行逆转录，合成cDNA；t
(4)用逆转录获得的cDNA构建测序文库；t
(5)进行高通量测序。t
本发明又一方面提供了RNA5mC甲基化检测方法，其特征在于包括以下步骤：t
(1)将RNA样品片段化；t
(2)对片段化的RNA样品进行重亚硫酸盐处理；t
(3)利用ACT三碱基随机六聚引物对重亚硫酸盐处理后的RNA样品进行逆转录，合成cDNA；t
(4)用逆转录合成的cDNA构建测序文库；t
(5)进行高通量测序；t
(6)对测序结果进行数据分析，获得RNA的5mC甲基化信息。t
本文中使用的术语“5mC”指5-甲基胞嘧啶。t
本文中的术语“重亚硫酸盐处理”是指使用重亚硫酸盐处理RNA，使RNA中所包含的未甲基化胞嘧啶转变为尿嘧啶。t
在一个具体实施方案中，上述RNA样品是待测RNA或包含待测RNA。待测RNA可以包含5mC。在优选的实施方案中，所述待测RNA是分离纯化的。在一个具体实施方案中，RNA样品中包含的待测RNA可以来自于人体、动物体、植物体或其器官、组织或细胞。在一个具体实施方案中，所述RNA可以是mRNA、tRNA、rRNA、总RNA或其它包含或可能包含5mC甲基化位点的RNA。在一个具体实施方案中，RNA样品是分离纯化的待测RNA。在一个具体实施方案中，RNA样品是非甲基化RNA标准品与分离纯化的待测RNA以1:100(重量比)的比例混合的混合物。在优选的实施方案中，所述非甲基化RNA标准品是小鼠DHFR基因的RNA序列。t
在一个具体实施方案中，所述“片段化”是将RNA样品打断，优选地打断为约100nt大小的RNA片段。t
在一个具体实施方案中，所述重亚硫酸盐处理包括重亚硫酸盐温育、脱盐和脱磺酸化步骤。t
在优选的实施方案中，所述重亚硫酸盐处理基本上完全将RNA样品中的未甲基化的胞嘧啶转化为尿嘧啶。其中“基本上完全”是指经重亚硫酸盐处理后，RNA样品中未甲基化胞嘧啶转化为尿嘧啶的转化效率大于等于99％，更优选大于等于99.6％，更优选大于等于99.7％，最优选大于等于99.9％。其中“转化效率”是指RNA样品中未甲基化胞嘧啶经重亚硫酸盐处理转化为尿嘧啶的比例。t
在一个具体实施方案中，所述ACT三碱基随机六聚引物是由A、C、T三种碱基随机合成的六聚引物的混合物。t
在一个具体实施方案中，用逆转录合成的cDNA构建测序文库包括II链合成、产物纯化、末端加“A”碱基(Atailing)、连接接头(Adapter)、纯化回收、PCR扩增和纯化回收PCR产物的步骤。t
在一个具体实施方案中，所述高通量测序是利用构建的测序文库在测序平台，优选第二代测序平台，更优选Hiseq2000测序平台上进行测序。t
在一个具体实施方案中，所述数据分析是对测序数据进行质量控制之后将其与人的全基因组序列(hg19)及转录组和外显子-外显子连接序列进行比对，并基于比对结果计算每个测到的胞嘧啶位点的甲基化水平。在一个具体实施方案中，使用公式M/(U+M)计算某个胞嘧啶的甲基化水平，U和M分别是测序数据中这个胞嘧啶位点上的非甲基化序列数和甲基化序列数。在一个具体的实施方案中，对于比对上的测序片段，如果该测序片段所包含的胞嘧啶位置超过30％的位置测序结果为C，在计算甲基化水平之前先将这些序列过滤掉。t
本发明又一方面提供了根据本发明的RNA5mC重亚硫酸盐测序文库构建方法构建的测序文库。t
本发明的又一方面提供了根据本发明的RNA5mC重亚硫酸盐测序文库构建方法构建的测序文库在RNA5mC甲基化检测中的应用。t
本发明的又一方面提供了ACT三碱基的随机六聚引物在RNA5mC甲基化检测中的应用，所述RNA5mC甲基化检测包括将RNA片段化、重亚硫酸盐处理、逆转录合成cDNA、高通量测序和数据分析的步骤，其特征在于：在对重亚硫酸盐处理后的RNA进行逆转录合成cDNA时，使用ACT三碱基的随机六聚体作为引物。t
下面将结合实施例对本发明的实施方案进行详细描述，但本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应当视为对本发明范围的限制。t
除非特别指明，以下实施例中所用的ACT三碱基随机六聚引物由Invitrogen公司合成。t
除非特别指明，以下实施例中所用的PCR仪为Veriti96孔热循环仪，购自AppliedBiosystems公司。t
除非特别指明，以下实施例中所用的试剂均为分析纯级的试剂，且可从常规渠道商购获得。t
实施例1t
1)RNA提取及片段化处理t
利用T7体外转录试剂盒(NEB、E2040S)从携带小鼠DHFR基因全长的质粒pcDNA3-HA-mDHFR上转录出DHFR的RNA序列，并以此作为鉴定重亚硫酸盐转换效率的标准品。利用Trizol法从人HeLa细胞中提取总RNA。根据mRNAPurificationKit(Ambion,61006)说明，从总RNA中分离纯化mRNA。将2μgmRNA与20ng体外转录得到的DHFR-RNA标准品混匀成mRNA样品，根据AmbionRNA片段化试剂盒的方法说明，将该mRNA样品片段化成100nt大小的片段。片段化条件为每个反应2μgmRNA/9μl，加入1μl片段化试剂，90℃反应1分钟后再加入1μl终止缓冲液。向反应后的mRNA产物中加入糖原和3M醋酸钠(pH5.2)及预冷的纯乙醇，放置-80℃过夜沉淀mRNA。t
2)重亚硫酸盐处理、脱盐及脱磺酸化t
离心回收mRNA沉淀，用100μl新鲜配制的重亚硫酸盐溶液(40％重亚硫酸钠，600μM氢醌,pH5.1)溶解mRNA沉淀,分成2管，每管50μl，置于75℃PCR仪上处理4小时。然后用MicroBio-spin6chromatographycolumns(Bio-Rad)对重亚硫酸盐处理后的样品进行脱盐处理。随后加入1MTris(pH9.0)于75℃PCR仪上1小时进行脱磺酸处理。最后加入糖原和3M醋酸钠(pH5.2)及预冷的纯乙醇，放置-80℃过夜沉淀mRNA。t
3)cDNA合成及文库构建t
离心回收mRNA沉淀，溶解后分出300ng重亚硫酸盐处理后的样品，利用SuperscriptIIReverseTranscriptasekit(Invitrogen)进行逆转录cDNA合成，其中使用自行合成的ACT三碱基随机六聚引物代替试剂盒中原有的随机六聚引物。然后根据KAPA公司的KAPAStrandedmRNA-SeqKit(platform)(Cat.No.KK8420)链特异性mRNA建库试剂盒说明进行后续的II链合成、产物纯化、A-tailing、加接头、纯化、PCR扩增、PCR产物纯化、文库检测过程-。t
为了鉴定ACT三碱基随机六聚引物的逆转录效率，将经重亚硫酸盐处理的和未经重亚硫酸盐处理的上述mRNA，分别使用不同的随机六聚引物进行逆转录，并将逆转录合成的cDNA分别使用DHFR特异性引物进行PCR扩增，对PCR扩增产物进行琼脂糖凝胶电泳分析(图3A)，并对电泳条带进行定量统计(图3B)，检测不同引物对DHFRRNA标准品的逆转录效率，所使用的引物分别是ACT三碱基随机六聚引物和传统的AGCT四碱基随机六聚引物。结果显示，对于未经重亚硫酸盐处理的mRNA样品，利用两种逆转录引物逆转录合成并扩增后DNA的量相当，说明其逆转录扩增效率相当；而对于重亚硫酸盐处理过的mRNA样品，利用ACT三碱基随机六聚引物进行逆转录合成并进行PCR后，DNA的量比ACTG四碱基随机六聚引物高3倍以上，说明其逆转录扩增效率较普通ACTG四碱基随机六聚引物高3倍以上。t
(4)高通量测序t
使用所构建的文库在Hiseq2000测序平台上进行测序，获得测序数据。t
(5)甲基化位点鉴定及甲基化水平计算t
质量控制：用FASTX-Toolkit软件对测序数据进行碱基质量控制、修剪引物等处理。t
序列比对：对于进行质量控制之后获得的序列，通过Bismark软件与人的参考基因组(hg19版本)上的原始序列进行比对。对于通过hg19原始序列比对未比对上的数据，再通过Bismark软件与hg19的转录组数据以及外显子-外显子连接构成的连接序列进行比对。t
甲基化水平分析：对于人参考基因组原始序列中的每个胞嘧啶位置，确定测序序列中在该位置表现为甲基化的序列数(即比对结果为C－C，以下称甲基化序列数)和表现为非甲基化的序列数(即比对结果为C－T，以下称非甲基化序列数)。对于与人参考基因组的转录组数据以及外显子-外显子连接序列比对的结果，与上述方法相同确定每个胞嘧啶位置的甲基化序列数和非甲基化序列数，并将其对应的位置信息转换为人参考基因组原始序列位置信息，并综合起来计算甲基化水平。使用公式M/(U+M)计算每个胞嘧啶位置的甲基化水平，其中U和M分别是在选定胞嘧啶位置上的非甲基化序列数和甲基化序列数。对于比对上的测序序列，如果在一个测序片段(即测序文库构建过程中打断形成的片段)中所包含的胞嘧啶位置中有超过30％的位置测序结果为C，鉴于RNA甲基化的低水平，如此高比例的甲基化胞嘧啶位点很可能为处理不充分导致，因此在计算甲基化水平之前先将这些测序片段过滤掉。t
为分析重亚硫酸盐处理的转化效率，计算DHFR标准品的重亚硫酸盐转化效率，将测序数据进行质量控制之后与DHFR标准品的序列进行比对，并基于比对结果计算标准品的转化效率，使用公式U'/(U'+M')计算标准品的转化效率，U'和M'分别是测序数据中DHFR的所有胞嘧啶位点上的非甲基化序列数之和和甲基化序列数之和。同样，对于比对上的DHFR的测序片段，如果该测序片段所包含的胞嘧啶位置超过30％的位置测序结果为C，在计算甲基化水平之前先将这些序列过滤掉。根据测序结果，U'＝8337，M'＝3097496，DHFR标准品的重亚硫酸盐转化效率为0.997(即99.7％)，即DHFRRNA标准品的转化效率达99.7％，代表了所测mRNA样品的重亚硫酸盐转化效率。t
对获得的人Hela细胞mRNA5mC甲基化位点的分布规律进行进一步分析，如以下实施例2-5所示。t
实施例2t
从Ensembl数据库下载人的基因组注释文件(版本号72)。对测序分析得到的甲基化位点,用Bedtools软件基于下载的注释文件进行注释。将注释结果中含有5mC修饰的基因基于ensembl中对基因的归类进行分类。结果表明含有5mC甲基化修饰的基因中87％属于蛋白编码基因(图4A)，其余的修饰基因为非编码基因，包含pseudogene，lincRNA,antisense,processedtranscript等。其中，5mC修饰的pseudogene有202个，lincRNA62个，antisense58个，processedtranscript39个(图4B)。t
实施例3t
为了检查5mC修饰在转录本上的各个区域的分布情况，首先将转录本按照Ensembl数据库的注释信息分为四个区域，5’非编码区(5’UTR)，蛋白编码区域(CDS)，内含子(Intron)和3’非编码区(3’UTR)。将实施例1中mRNA上的5mC修饰位点位置信息提取出来，通过Bedtools软件与Ensembl数据库下载的注释信息进行分析，检测这些甲基化位点在mRNA的四个区域(5’UTR、CDS、Intron以及3’UTR)的数量分布，结果表明这些甲基化位点在CDS中的分布为35％，内含子中为28％，3’非编码区为24％，5’非编码区为13％(图5A)。统计测序结果中测到的mRNA上C位点(覆盖次数大于等于1)在四个区域(5’UTR、CDS、Intron以及3’UTR)的分布比例，将其作为期望的甲基化C位点的分布，用期望的甲基化位点的分布和实际甲基化位点的分布来检测5mC位点在四个区域分布的偏好性，结果表明，和期望值相比，5mC修饰位点显著富集在蛋白编码区域和5’非编码区以及3’非编码区(图5B)。t
实施例4t
为了检查5mC修饰在转录本上的分布特征，将mRNA的5’非编码区、蛋白编码区域以及3’非编码区三个区间分别按照总长度标准化成总成为100的长度，并通过Bedtools软件与其进行位置关联，统计在1-100各个区间上的甲基化位点个数，并将每个小区间的甲基化位点的个数除以所有甲基化位点的个数之和，得到甲基化位点在mRNA整体水平上的百分比分布曲线(图6)。结果表明，在CDS起始位点附近，5mC位点分布相对于其他区域有着显著的富集。t
实施例5t
5mC甲基化位点可以根据其附近的碱基信息分为CpG，CHG及CHH，其中H为A、C、或U。对筛选得到的甲基化位点，提取其附近的2-3nt的序列，分析其中CpG，CHG及CHH的分布数量及比例。结果表明，相对于CpG和CHG，5mC甲基化位点主要分布在CHH区域(图7A)。同时，我们以甲基化位点为中心，将其上下游各延长2nt得到以甲基化位点为中心的5nt序列,通过fastaFromBed程序对序列进行提取，并利用Weblogo展示甲基化位点附近5nt的序列(图7B)。t
实施例6t
为了进一步证明测序结果的可靠性，从实施例1的测序数据中选取了含有5mC修饰的3个PLOD3、COL4A5和FAM129B，利用重亚硫酸盐处理并结合Sanger测序方法对其中的5mC甲基化位点进行验证，即将上述基因的对应甲基化位点两侧设计PCR引物，以经重亚硫酸盐未处理和处理的RNA为模板进行逆转录PCR，通过琼脂糖凝胶电泳检测PCR扩增产物，并将对应位置的条带切下回收DNA片段，进行Sanger测序。图8显示了验证结果。图8的左图为重亚硫酸盐未处理和处理的Sanger测序峰图和序列比对图。图8的右图示意性地以高通量测序结果中10条序列(reads)的测序结果为代表显示了测得的甲基化胞嘧啶的次数，代表了该位点的甲基化水平；其中用圆圈表示不同序列中甲基化位点左右两侧胞嘧啶的甲基化状态，空心圆圈表示非甲基化即重亚硫酸盐转化后测序结果为T的胞嘧啶，实心圆圈表示甲基化即未转化的胞嘧啶。结果表明三个基因中分别存在一个经重亚硫酸盐处理后未发生转化的甲基化胞嘧啶，与高通量测序结果吻合，进一步证实了本发明方法的可靠性和准确性。t
参考文献t
[1]GordonA.FASTX-Toolkit.http://hannonlab.cshl.edu/fastx_toolkit/.t
[2]MartinM.Cutadaptremovesadaptersequencesfromhigh-throughputsequencingreads.EMBnetJ.2011；17(1):10-12.t
[3]BolgerAM,LohseM,UsadelB.Trimmomatic:aflexibletrimmerforIlluminasequencedata.Bioinformatics.2014；30(15):2114-20.t
[4]Krueger,F.Andrews,S.R.Bismark:aflexiblealignerandmethylationcallerforBisulfite-Seqapplications.Bioinformatics.2011；27(11):1571-2.t
