大家好，这是专注表观组学十余年，领跑多组学科研服务的易基因。
本期，我们讲讲甲基化RNA免疫共沉淀(MeRIP-seq/m6A-seq)实验怎么做，从技术原理、建库测序流程、信息分析流程和研究套路等四方面详细介绍。
表观转录组指RNA序列不发生改变的情况下，由RNA上的化学修饰调节基因表达的现象。胞内RNA的修饰超过100种，其中大部分的表观修饰发生在tRNA及其他非编码RNA上[1]，发生在mRNA上的修饰无论是从种类上还是数量上都较少。下图展示了部分发生在真核生物mRNA上的化学修饰[2]：
图1： mRNA甲基化修饰示意图
mRNA上含量最高的修饰是腺苷酸N6位的甲基化修饰(N6-methyladenosine，m6A)，m6A参与mRNA生命周期的各个方面[2]。m6A修饰主要发生在RRACH (R = G or A and H = A,C, or U)这类motif上，这类motif主要在终止密码子和3’UTR处富集，即m6A主要发生在终止密码子和3’UTR附近。m6A修饰由m6A甲基转移酶（writers）催化产生，被m6A去甲基酶（erasers）去除，被m6A结合蛋白（readers）识别并发挥功能[3]。m6A修饰广泛存在于各类生物中，包括病毒、从酵母、植物和人类等高等动物[1,3]。
研究表明，m6A已知的最主要功能是调控mRNA的稳定性：胞质内经过m6A修饰的mRNA可以被YTHDF2识别，使其富集到Processing body（P-body）从而加速mRNA的降解。此外，m6A修饰还可以改变RNA的二级结构、调控microRNA的靶标识别来调控mRNA的稳定性。在核内，m6A修饰可以调控RNA的剪接、出核过程，从而调控基因表达。m6A还可能与DNA甲基化存在互作关系。下图展示了部分目前已知的m6A与生物功能的关系：
图2：m6A与生物学功能
早期的实验方法已经在rRNA和mRNA中鉴定到m6A的存在，但由于实验技术的限制，对m6A的研究一直没有大的进展。近年来，RNA去甲基化酶FTO的报道，使得对m6A修饰的研究再次进入人们的视野。随后，MeRIP-seq/m6A-seq实验方法的建立使得有机会对该修饰的分布和功能进行深入研究[1]。
利用MeRIP-seq/m6A-seq技术鉴定全基因组范围内具有m6A修饰的区域。其原理为通过特异识别m6A修饰的抗体，对细胞内具有m6A修饰的RNA片段进行免疫共沉淀。对沉淀下来的RNA片段进行高通量测序，结合生物信息学分析，即可在全基因组范围内对m6A修饰的状况进行系统研究。MeRIP-seq/m6A- seq是目前研究m6A修饰使用最广泛的技术之一。
从样品处理到最终数据获得中每一个环节都会对数据质量和数量产生影响，而数据质量又会直接影响后续信息分析的结果。为了从源头上保证测序数据的准确性、可靠性，需对样品处理、建库、测序每一个生成步骤都严格把控，从根本上保证了高质量数据的产出。流程图如下：
（一）Total RNA样品检测
对RNA样品的检测主要包括3种方法：
（1）琼脂糖凝胶电泳分析RNA降解程度以及是否有污染，检测具有明显的18S或28S主带，且条带清晰；
（2）Qubit 2.0对RNA浓度进行精确定量，总RNA 检测总量不低于10ug；
（3）Agilent 2100精确检测RNA的完整性，RIN值不低于7.5；
1. RNA片段化：
取10ug总RNA，加入RNA Fragmentation Reagents（Invitrogen）在Thermomixer中70℃反应10min，将RNA打断成片段，片段大小约100nt，用乙醇法将片段化的RNA沉淀；
2. m6A富集：
将含protein A和protein G的磁珠用IP buffer(150 mM NaCl, 10 mM Tris-HCl pH 7.5)洗涤，然后与5 ug m6A抗体(Millipore) 4℃孵育2h；用 IP buffer洗涤两次，再用 IP buffer将磁珠重悬，加入片段化的RNA，4℃下翻转4h；用IP buffer在4℃下洗涤磁珠3次，再用m6A竞争性洗脱液，在4℃下孵育1h。将含有洗脱m6A RNA的上清收集到新的试管中，用苯酚:氯仿:异戊醇(125:24:1)进行纯化；
3. 文库构建：
用SMARTer® Stranded Total RNA-seq Kit v2 - Pico Input Mammalian User Manual根据说明书对IP与Input样品进行反转录及文库构建；利用AMPure XP bead进行片段大小选择，获得最终文库。
构建原理图如下：
图3：富集流程示意图
文库构建完成后，先使用Qubit2.0进行初步定量，稀释文库至1ng/ul，随后使用Agilent 2100对文库的insert size进行检测，insert size符合预期后，使用qPCR方法对文库的有效浓度进行准确定量（文库有效浓度> 2nM），以保证文库质量。
文库检测合格后，把不同文库按照有效浓度及目标下机数据量的需求pooling后在illumina Nova平台测序，测序策略为PE150。
原始下机数据为FASTQ格式，是高通量测序的标准格式。FASTQ文件每四行为一个单位，包含一条测序序列（read）的信息。该单位第一行为read的ID，一般以@符号开头；第二行为测序的序列，也就是reads的序列；第三行一般是一个+号，或者与第一行的信息相同；第四行是碱基质量值，是对第二行序列的碱基的准确性的描述，一个碱基会对应一个碱基质量值，所以这一行和第二行的长度相同。以下为一条read信息的示例：
图4：FASTQ格式示例
原始下机数据包含建库时引进的接头序列以及质量过低的碱基，这些因素会导致后续比对到基因组的reads较少，从而导致得到的信息较少，因此需要进行过滤。
利用Trimmomatic软件对原始数据进行去除接头序列及低质量碱基等质控步骤，所用参数为：ILLUMINACLIP:MeRIP-PE.fa:2:30:10:1:true SLIDINGWINDOW:30:15 AVGQUAL:15 LEADING:15 TRAILING:15 MINLEN:30
过滤后的数据需要比对至参考基因组，基因组上发生m6A修饰的区段会有大量reads比对上去，从而形成“峰”（peak），根据峰的位置即可判断基因组的哪些位置发生了甲基化。
比对采用hisat2[4]，该软件可快速将短序列比对至参考基因组，并且可以考虑并处理拼接位点，尤其适用于RNA数据的比对。所用参数为默认参数。比对完成后，对结果路径进行去除多重比对、去除低质量比对等过滤，从而得到更精确的比对结果。
MeRIP-seq将RNA上发生m6A修饰的区域富集后进行测序，因此发生m6A修饰的区域，IP文库所覆盖的reads数会显著高于input文库，从而形成“峰（peak）”。检测这些峰的位置即可得到RNA上发生m6A修饰的位置。
鉴定得到m6A修饰后，再对m6A修饰进行注释、分布统计、motif鉴定等分析。
差异m6A修饰（即差异peak）的鉴定采用R包exomePeak[5,6]。该R包先将需要比较的样本的peak进行合并，再计算每个样本的合并后peak内累积的reads数目；再对这些reads进行标准化，使两组样本处于可比较的水平；最后检验两组样本在该peak内的reads数目是否有显著差异。所用参数为：WINDOW_WIDTH=200 SLIDING_STEP=30 DIFF_PEAK_ABS_FOLD_CHANGE = 1.5 FOLD_ENRICHMENT=1.5 FRAGMENT_LENGTH=200。
鉴定差异m6A修饰后，再对差异m6A修饰进行注释、分布统计、motif鉴定等分析。
m6A-seq的input文库相当于RNA-seq文库，可以用于分析基因表达量及鉴定差异表达基因。通过定位到基因组区域或者基因外显子区的测序序列（read）的计数来估计基因的表达水平。Read计数除了与基因的真实表达水平成正比外，还与基因的长度和测序深度成正相关。为了使不同基因、不同实验间计算的基因表达水平具有可比性，利用TPM来对read数量进行标准化，该算法考虑了测序深度和基因长度对计数的影响，是目前最为常见的基因表达水平计算方法之一。
图5：TPM计算公式
注：Ri表示该基因的read count数目；Li表示该基因的长度（kbp)，选取最长的转录本作为该基因的长度。表达量计算软件为StringTie[7]，使用默认参数，用TPM展示基因表达量。
通过对转录本进行组装，可以发现注释文件外的新转录本。对这些新转录本进行一系列严格条件筛选，鉴定出非编码RNA。
（1）使用软件StringTie[7]对转录本进行组装，获得每个样品中所有的转录本，将没有在注释文件中出现的转录本标记为新转录本。对新转录本筛选条件为：a)对于每个样品的组装注释文件，过滤掉转录本exon个数 >=2 但FPKM<=0.5或者exon个数>=1但FPKM<=1的转录本[8];b)过滤单个转录本拼接结果中长度低于200bp的转录本；c)将以上单个过滤后的组装转录本用StringTie[7]合并成最终组装转录本；d）通过Cuffcompare[9]软件分析未知转录本与已知转录本的结构关系。
（2）候选lncRNA编码能力预测。lncRNA是一类非编码RNA，通过预测候选lncRNA的编码能力，进一步鉴定出真正的lncRNA。采用多款主流编码能力预测软件对候选lncRNA进行编码能力预测。预测软件包括CNCI[10]和CPC[11]。
（3）针对基因组有注释的lncRNA基因信息，计算known lncRNA的表达情况。
（4）表达量计算。用TPM标准化lncRNA的表达量。