易基因｜干货：手把手教你做RNA m5C甲基化测序分析（RNA-BS）

大家好，这是专注表观组学十余年，领跑多组学科研服务的易基因。

本期，我们讲讲m5C RNA甲基化重亚硫酸盐测序测序（RNA-BS）实验怎么做，从技术原理、建库测序流程、信息分析流程等方面详细介绍。

一、m5C RNA甲基化测序（RNA-BS）技术原理

m5C 是 RNA 百余种修饰中研究较多的⼀种。m5C 存在于 tRNA 上时，可以对翻译进⾏调节；存在于 rRNA 上时，可以对核糖体的⽣物合成进⾏调控；存在于 mRNA 上时，则可以影响 mRNA 的结构、稳定性及翻译过程。

早在 1970 年，就已经在 mRNA 上发现 m5C 修饰的存在，但是由于技术限制，mRNA 上 m5C 修饰的研究多年来进展缓慢。近⼏年来，多种 m5C 研究⽅法的出现（如 MeRIP-seq、 miCLIP、Aza-IP、RNA-BS-seq 等）使得 RNA m5C 修饰的研究再次进⼊⼈们的视野。

其中，RNA-BS-seq 是⼀种能够从单碱基分辨率⽔平检测 m5C 的强有⼒的技术。该技术利⽤重亚硫酸盐处理 RNA，使RNA 上没有发⽣修饰的 C 被转化为 U，⽽发⽣ m5C 修饰的 C 碱基则保持为 C。经过 PCR，U 转变成 T，这样便将 m5C 与C 区分开来。结合⾼通量测序，就可以从转录组范围检测 m5C 修饰。

RNA甲基化测序原理示意图如下：

该技术的优势包括：（1）检测范围为 mRNA；（2）单碱基分辨率；（3）⾼准确性。

二、m5C RNA甲基化测序（RNA-BS）建库测序

项⽬最终数据的质量受到样品检测、建库、测序等每⼀环节的影响。为从源头保证测序数据准确可靠，易基因承诺在数据的所有实验室⽣产环节严格把关，确保⾼质量数据的产出。RNA甲基化(调取)建库测序的流程图如下：

（一）Total RNA样品检测

对RNA样品的检测主要包括3种方法：

（1）琼脂糖凝胶电泳分析RNA降解程度以及是否有污染，检测具有明显的18S或28S主带，且条带清晰；

（2）Qubit 2.0对RNA浓度进行精确定量，总RNA 检测总量不低于75ug；

（3）Agilent 2100精确检测RNA的完整性，RIN值不低于7.5。

（二）⽂库构建与质检

（1）文库构建：

① mRNA调取：使⽤Oligo（dT）磁珠，将mRNA从Total RNA⾥⾯调取出来。

② 重亚硫酸盐处理：将mRNA进⾏重亚硫酸盐处理，使⽤EZ RNAMethylationTM kit中的Zymo-SpinTM IC Column柱纯化mRNA。

③ 反转录⼀链合成：加⼊ N6 Pri mer，65℃，5min，置于冰上；加⼊First strand buffer、 dNTP、 DTT, RNase Inhibitor,混匀，加⼊SuperscriptⅡ，混匀后进⾏反转录⼀链的合成。

④ 反转录⼆链合成：加⼊Second strand buffer，dNTP，RNase H，DNA Pol Ⅰ,混匀后置于Thermomixer中16℃反应后，磁珠纯化，EB洗脱。

⑤ 洗脱液通过末端修复、末端加A、接头连接后，进⾏PCR扩增反应，扩增产物即为最后的⽂库。

构建原理图如下：

注：

测序接头：包括P5/P7，index和Rd1/Rd2 SP三个部分（如上图所示）。其中P5/P7是PCR扩增引物及flow cell上引物结合的部分，index提供区分不同⽂库信息的Rd1/Rd2，SP即read1/read2 sequence primer，是测序引物结合区域，测序过程理论上由Rd1/Rd2 SP向后开始进⾏。

（2）文库质检：

⽂库构建完成后，先使⽤Qubit2.0进⾏初步定量，稀释⽂库⾄1ng/ul，随后使⽤Agilent 2100对⽂库的insert size进⾏检测，insert size符合预期后，使⽤qPCR⽅法对⽂库的有效浓度进⾏准确定量（⽂库有效浓度> 2nM），以保证⽂库质量。

（三）上机测序

库检合格后，把不同⽂库按照有效浓度及⽬标下机数据量的需求pooling后在HiSeq平台测序，测序策略为PE150。测序的基本原理是边合成边测序（Sequencing by Synthesis）。在测序的flow cell中加⼊四种荧光标记的dNTP、DNA聚合酶以及接头引物进⾏扩增，在每⼀个测序簇延伸互补链时，每加⼊⼀个被荧光标记的dNTP就能释放出相对应的荧光，测序仪通过捕获荧光信号，并通过计算机软件将光信号转化为测序峰，从⽽获得待测⽚段的序列信息。测序过程如下图所示：

三、m5C RNA甲基化测序（RNA-BS）信息分析流程

四、m5C RNA甲基化测序（RNA-BS）质控分析

（一）测序数据说明

测序⽚段被⾼通量测序仪测得的图像数据经CASAVA碱基识别转化为序列数据（reads），⽂件为fastq 格式，其中主要包含测序⽚段的序列信息以及其对应的测序质量信息。fastq格式⽂件中每个read由四⾏描述信息组成，如下所示：

图：FASTQ格式示例

上述⽂件中第⼀⾏以“@”开头，随后为Illumina测序标识符(Squence Identifiers)和描述⽂字；第⼆⾏是测序⽚段的碱基序列；第三⾏以“+”开头，随后为Illumina测序标识符(也可为空)；第四⾏是测序⽚段每个碱基相对应的测序质量值，该⾏中每个字符对应的ASCII值减去33，即为该碱基的测序质量值。

测序过程本身存在发⽣机器错误的可能性，测序错误率分布检查可以反映测序数据的质量，序列信息中每个碱基的测序质量值保存在FASTQ⽂件中。如果测序错误率⽤e表示，Illumina的碱基质量值⽤Q 表示，则有：Q =-10log10(e)。Illumina Casava 1.8版本碱基识别与Phred分值之间的简明对应关系⻅下表：

（二）测序数据质控

原始下机数据包含建库时引进的接头序列以及质量过低的碱基，这些因素会导致后续⽐对到基因组的reads较少，从⽽导致得到的信息较少，因此需要进⾏过滤。过滤内容如下：

切除接头(adapter)序列；
如果整条 read 的碱基平均质量值低于15则去掉整条 read;
从序列的开头和末尾截掉碱基质量值⼩于3的碱基；
去除clean后短于20bp⻓度的序列。

（三）数据质量评估

经过原始数据过滤、测序错误率检查、GC含量分布检查数据汇总。

(四) ⽐对质量评估

BS-RNA过程包括三个主要步骤：预处理，⽐对和注释。

第⼀步是对参考基因组序列、测序数据和基因注释⽂件进⾏预处理：

（1）将参考基因组序列并⾏转换两次：（A）胞嘧啶被胸腺嘧啶取代，（B）⻦嘌呤被腺嘌呤取代。这种基因组序列转换只需要在第⼀次使⽤参考基因组序列时进⾏，这意味着它可以重新⽤于使⽤相同参考基因组序列的所有后续分析。

（2）T-rich reads 的胞嘧啶被胸腺嘧啶取代，⽽A-rich reads中的⻦嘌呤被腺嘌呤取代。

（3）修改GTF格式的基因注释⽂件以适应转换后的参考基因组序列。每个注释⾏同时转换两次：“C-T”和“G-A”分别附加到基因注释⽂件中的染⾊体标记上。

接下来，BS-RNA调⽤HISAT2程序，根据修改的注释基因⽂件构建替代剪接，并将预处理的读数与转换后的参考基因组序列对⻬。BS-RNA过滤掉mapping到参考基因组序列的两种类型的reads：（1）mapping到多个位置的reads和（2）mapping到错误链的reads（将富含T的readsmapping到将参考序列转换为T或将参考序列转换为T的反向补体，将A丰富的reads ⽐对到将C转换为T或反向补体将G转换为A的参考序列）。mapping步骤完成后，BS-RNA 将提供原始mapping⽂件（SAM 格式）、过滤mapping⽂件（SAM 格式）和mapping报告⽂件。

⽐对率定义为对于单个样品，将read与参考基因组⽐对，⽤能⽐对上的read数量除以该样品总read数量，即为该样品read的⽐对率。⽐对率计算公式如下：

（五）覆盖率评估

（1）甲基化⽔平计算

甲基化⽔平可根据未转化为 T 的 C 与转化为 T 的 C 的 reads 的⽐例计算得到，即：

Beta-value = C-reads / (C-reads + T-reads) * 100%

其中，Beta-value 即为该胞嘧啶的甲基化⽔平，C-reads 为覆盖该位点的⽀持甲基化的reads 数⽬（测得该位点为 C 的reads），T-reads 为覆盖该位点的不⽀持甲基化的 reads 数⽬（测得该位点为 T 的 reads）。计算原理示意图如下：

（2）C位点数统计

将reads⽐对到基因组后，⽐对到不同位点的reads数（测序深度）不同，测序深度过低会导致计算的甲基化率不可信。因此，统计了所有C位点的测序深度。分别统计甲基化数据三种不同类型的C碱基的测序深度。

五、m5C RNA甲基化测序（RNA-BS）差异甲基化位点(DMC)的鉴定及统计

（1）差异甲基化位点（DMC）的鉴定

（2）DMC的注释

鉴定出的DMC包含染⾊体、起始位置、终⽌位置等信息。根据DMC的位置信息，结合基因组注释信息中所有基因的位置信息及各个基因元件（5utr, cds, intron, 3utr, ncRNA, tRNA）等位置信息，鉴定DMC与哪些基因的哪些基因元件有重叠，以此来判断DMC修饰哪些基因的哪些基因元件。

（3）DMC修饰基因的统计

根据DMC的注释⽂件，提取出DMC修饰的基因及其的信息，以更加⽅便地查看DMC修饰的基因。

（4）DMC在染⾊体上的分布

根据DMC的位置信息，统计DMC落在哪些染⾊体上，并⽤图形展示，以了解DMC在染⾊体上的分布有偏好性。

（5）DMC在基因元件上的分布

同样地，根据DMC的位置信息，分别统计Hyper DMC及Hypo DMC 落在哪些基因元件上。

（6）DMC修饰基因的功能富集分析

基因本体（ Gene Ontology， GO）是基因功能国际标准分类体系，提供了⼀套动态更新的标准词汇表来描述⽣物体中基因和基因产物的属性，可以挖掘出⼀些⽣物学相关的途径。 GO分为三个Ontology，分别是：分⼦功能（MolecularFunction， MF）、细胞组分（ Cellular Component， CC）和⽣物过程（ Biological Process， BP）。

KEGG（ Kyoto Encyclopedia of Genes and Genomes，京都基因与基因组百科全书）是基因组破译⽅⾯的数据库。在给出染⾊体中⼀套完整基因的情况下，它可以对蛋⽩质交互⽹络在各种细胞活动起的作⽤做出预测。KEGG Pathway显著性富集分析应⽤超⼏何检验，找出与整个基因组背景相⽐，在差异甲基化修饰的基因中显著富集的Pathway。

将鉴定出的DMC所修饰的基因，利⽤GO和KEGG数据库进⾏功能富集分析。

对位于cds区的Hyper-DMC和Hypo-DMC修饰的基因进⾏功能富集分析：

① 对实验组相对于对照组甲基化⽔平升⾼的DMC（All-Hyper DMC）修饰的基因做功能富集分析；

② 对实验组相对于对照组甲基化⽔平降低的DMC（All-Hypo DMC）修饰的基因做功能富集分析；

富集分析采⽤Fisher检验，结合BH校正。富集分析结果包括表格和图⽚两部分，其中，表格为所有富集到的GO/KEGG条⽬，包括显著和不显著的。

易基因RNA m5C甲基化测序(RNA-BS)技术介绍

易基因提供适用于不同科研需求的m5C甲基化测序技术：

常规mRNA m5C甲基化测序（RNA-BS）：
mRNA分离后首先通过亚硫酸盐处理，非甲基化的C转变为U，m5C修饰的碱基保持不变，结合高通量测序，可以对RNA上的每一个C碱基修饰进行定位与定量。
常规mRNA +lncRNAm5C甲基化测序（全转录组RNA-BS）：

易基因科技建立的升级版m5C RNA甲基化测序服务，去除人rRNA后，剩余RNA经重亚硫酸盐处理后，结合高通量NGS策略，可在全转录组范围内单碱基分辨率地检测基因m5C甲基化修饰分布。

技术优势：

高深度：超高深度重亚硫酸盐处理，检测准确性极高；
高通量：结合高通量NGS，全转录组范围内检测；
单碱基：单碱基分辨率，快速检测和分析RNA中的m5C。
高准确：精确的检测mRNA等每一个C碱基的的修饰水平。

研究方向：

与m6A甲基化类似，m5C甲基化已被证明与肿瘤、神经系统紊乱、代谢性疾病、病毒感染以及个体发育等密切相关。
此外，RNA甲基化（m5C）与人类疾病密切相关，其功能涉及调控干细胞应激、细胞毒性应激、mRNA出核和植物细胞发育及基因表达等方面。

以上就是关于甲基化RNAm5C甲基化测序（RNA-BS）实验流程和分析思路的介绍，易基因科技提供全面的RNA甲基化研究整体解决方案，技术详情了解请致电易基因0755-28317900。

参考文献：

[1] Ashburner, M. and C. A. Ball, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium.Nat Genet, 2000, 25 (1): 25-9.

[2] Yang X, et al.5-methylcytosine promotes mRNA export - NSUN2 as the methyltransferase and ALYREF as an m5C reader. Cell Res. 2017 May;27(5):606-625.

[3] Frank Jühling et al. metilene: Fast and sensitive calling of differentially methylated regions from bisulfite sequencing data.Genome Research, 2016, 26: 256-262.

[4] Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic acids research, 2000,28(1): 27-30.

[5] A F L , A L H , C J W B , et al. BS-RNA: An efficient mapping and annotation tool for RNA bisulfite sequencing data[J].Computational Biology & Chemistry, 2016, 65:173.

[6] Bohnsack K E, Höbartner C, Bohnsack M T. Eukaryotic 5-methylcytosine (m5C) RNA methyltransferases: mechanisms,cellular functions, and links to disease[J]. Genes, 2019, 10(2): 102.