转录组讲师带你读文献-m6A和RNA测序结合

我在我在04-转录组笔记推文任务列表（半年期）里面安排了6个经典综述和10篇转录组应用文献给大家，可惜愿意沉下心了认真苦学的并不多。（https://share.mubu.com/doc/14uneHKvPg）

所以安排转录组讲师给大家做一下领读：

下面是转录组讲师的投稿

Date：2020-08-06 Author：jzhang

1 文章信息

• 标题：Epitranscriptomic profiling of N6-methyladenosine-related RNA methylation in rat cerebral cortex following traumatic brain injury
• 标题：外伤性脑损伤后大鼠大脑皮层N6 -甲基腺苷相关RNA甲基化的表观转录组谱分析
• 作者：Yan Zhao（通讯），单位：武汉大学中南医院
• 杂志：Molecular Brain （IF：4.686）， 发表日期：2020 Jan 28

摘要

Background: m6A是真核mRNA最普遍的转录后修饰。据报道，在哺乳动物中枢神经系统中，对感觉体验、学习和损伤有一种刺激依赖m6A调节。大鼠大脑皮层创伤性脑损伤(TBI)的mRNA m6A甲基化模式尚未被研究。

Results: 在本研究中，我们通过甲基化RNA免疫沉淀测序(MeRIP-Seq)对大鼠皮层的mRNA m6A甲基化进行了全基因组图谱分析。TBI后，大鼠大脑皮层METTL14和FTO表达明显下调。利用MeRIP-Seq，我们一共鉴定出2165个显著变化的峰，其中1062个峰显著上调，1103个峰显著下调。这些m6A峰分布在1850个基因中。通过分析m6A峰和mRNA的表达，发现有175个mRNA的甲基化和表达水平在TBI后发生了显著的改变。此外，我们发现功能性FTO对TBI神经损伤的修复是必要的，但使用FTO抑制剂FB23-2对TBI大鼠的空间学习记忆能力没有影响。

Conclusion：本研究探讨了大鼠脑外伤后mRNA m6A甲基化模式，确定FTO可能是TBI表观遗传修饰的干预靶点。

Keywords: m6A methylation, Epigenetic modification, FTO, Traumatic brain injury, Rat cortex

2 数据和方法

MeRIP-Seq和RNA-Seq：18只老鼠，分为6组，每组三个生物学重复，三个样本合并为一个进行测序。6组为3组TBI，3组对照。

质控：Trimmomatic

去核糖体RNA：the SortMeRNA

比对：HISAT2

去重：Picard

m6A-seq data 质量评估：Trumpet R package

peak calling以及差异分析 ：MeTDiff

差异peak注释：ChIPseeker

这个研究并没有明确提供其测序数据或者说表达量矩阵，我看了看文章的附件有没有：

• Table S1. PCR primers used in this study
• Table S2. The detailed information of raw data
• Table S3. Statistical analysis of reads mapped in
• Table S4. The m6A peak number on each gene
• Table S5. The detailed information of significantly changed m6A peaks
• Table S6. The detailed information of significantly changed mRNA
• Table S7. The detailed information of conjoint analysis between mRNA expression and m6A methylation
• Table S8. The raw data of mNSS tests
• Table S9. The raw data of MWM tests

3 Results

3.1 脑外伤后大鼠皮层METTL14和FTO均下调（qRT-PCR结果）

为了探讨探讨TBI后大鼠皮层m6A甲基化状态是否发生改变，使用qRT-PCR 对六个甲基转移酶和去甲基化酶METTL3, METTL14, WTAP, VIRMA, FTO, and ALKBH5进行了评估。与假手术组相比，TBI组甲基化酶METTL14和去甲基化酶表达均显著下调(p < 0.05)。METTL3、WTAP、VIRMA、ALKBH5 4个基因在脑外伤后大鼠皮层的表达无明显变化。METTL14和FTO表达的改变可导致大鼠脑外伤后皮层m6A甲基化的动态变化，如下图：

Fig. 1 Reverse transcription polymerase chain reaction (RT-PCR) of METTL3, METTL14, WTAP, VIRMA, FTO, and ALKBH5

3.2 甲基化RNA免疫沉淀测序结果概述

为了测量sham组和TBI组的m6A甲基化水平，我们使用MeRIP-Seq进行了全转录组m6Aseq分析。我们发现在MeRIP-Seq样本中平均产生11.67 Gb的测序数据，输入样本产生13.46 Gb的测序数据。去除3 '和5 '端上的接头序列、低质量read和低质量碱基后，在MeRIP-Seq样本中cleandata的平均数据量为10.66 Gb，在输入样本中cleandata的平均数据量为12.28 Gb。将每个样本的clean reads与大鼠基因组进行比较，获得参考基因组上的位置信息，平均95.43%的reads比对到了大鼠基因组上，多比对率5.48%，唯一比对率89.95%并且用于后续分析。

3.3 m6A峰的拓扑分布

与MeRIPSeq样本之间的测序数据及其相应输入相比较，我们发现了来自所有样本4900个转录基因的11,789个不同的m6A峰。m6A甲基化位点最多的染色体分别为1、3、10号染色体，对应的m6A甲基化位点个数分别为1386、921和905（图2a）。

Fig.2a The distribution pattern of m6A peaks in different chromosomes.具体来说，m6A修饰位点个数在1 ~ 19个之间，64.31%的基因有1个或2个甲基化位点，35.69%的基因有3个或3个以上甲基化位点。例如，位于4号染色体上的Pclo有最多的甲基化位点（19个点)。此外,我们分析了mRNA基因结构中m6A甲基化峰值的分布模式，我们发现峰大多分布在外显子区域，在3端(终止密码子附近)有显著的富集峰和在5 '端有富集(起始密码子附近)。在CDS区域，输入样本的reads高于IP样本的reads(图2b)。

Fig. 2b The distribution patterns of m6A methylation peaks in gene structures of mRNA. The peaks were mostly distributed on the exon region and there was a distinct enrichment peak at the 3′ terminate (near the stop codon) and an enrichment peak at the 5′ terminate (near the start codon)

我们选择了三个基因来显示m6A甲基化模式。Trmt10c 的峰集中与CDS区域，Coro7的峰位于CDS和3’UTR区域，Plk的峰位于3’UTR、CDS和5’UTR区域(图2c)。

Fig. 2c, IGV plot shows directly the peaks in the genes of Trmt10c, Coro6 and Plk. The peak of Trmt10c was located at the CDS region, the peak of Coro7 was located at the CDS and 3′ UTR regions, and the peak of Plk was located at the 3′ UTR, the CDS, and the 5′ UTR regions

3.4 TBI后m6A甲基化变化显著

为了揭示m6A甲基化在TBI中的作用，我们比较了TBI组和假手术组大鼠皮层的甲基化水平。TBI组平均有15,305个峰值，假手术组平均有16,714个峰值(Fig. 3a) 。

Fig. 3a, The peak number of each group. An average of 15,305 peaks in the TBI group and an average of 16,714 peaks in the sham group were identified.

与输入样本相比，TBI组的对数平均富集倍数为4.11，假手术组的对数平均富集倍数为3.72(图3b)。

Fig. 3b The enrichment of each group. The average logarithmic fold-enrichment of TBI group was 4.11, and the average logarithmic fold-enrichment of sham group was 3.72

TBI组平均峰长为3426.14 bp, sham组平均峰长为3831.73 bp(图3c)。

Fig. 3c The peak length of each group. The average peak length of TBI group was 3426.14 bp, and the average peak length of sham group was 3831.73 bp

共鉴定出2165个显著变化峰(p < 0.01, FDR < 0.05)，其中显著上调峰1062个，显著下调峰1103个（Fig. 4a）。前20个差异甲基化mRNA见表1。如上所示，在4900个基因中，有1580个基因的甲基化峰存在差异。各峰的对数平均富集倍数量为2.53(图4b)，各峰的平均长度为4358.7 bp(图4c)。p值在各峰间的分布如图4d所示。

变化显著的峰主要分布在CDS(49.95%)、内含子(22.86%)、3 ' UTR(19.82%)和5 ' UTR(7.34%)，基因间区仅占0.02%(图5a)。多个基因结构中出现了许多具有显著变化的峰;556个峰分布在CDS和3 ' UTR区，532个峰分布在内含子区，495个峰分布在CDS区(图5b)。

Fig. 5 The distribution of significantly changed peaks after TBI. a, Sector graph shows the ratio of peaks in each region. The peaks with significant changes were mainly distributed in the exon (49.95%), CDS (22.86%), 3′ UTR (19.82%), and 5′ UTR (7.34%) and only 0.02% were found in the intergenic region. b, Showing the exact distribution pattern of significantly changed peaks post-TBI. There were 556 peaks distributed in the CDS and 3′ UTR regions, 532 peaks in the CDS and intron regions, and 495 peaks in the CDS region

图5c中显示了三个具有代表性的基因，其峰值发生了显著变化。TBI后CD14、Dll4和Sox7的m6A水平分别显著上调3.63倍、2.10倍和3.14倍。

Fig. 5c, Three representative genes with significantly changed peaks. The m6A level of CD14, Dll4, and Sox7 were significantly up-regulated after TBI by 3.63 fold, 2.10 fold, 3.14 fold, respectively

3.5 GO和KEGG分析揭示了m6A甲基化背后的生物学信息

GO分析显示m6A甲基化与三部分生物信息相关。分子功能:蛋白结合、poly(A) RNA结合、染色质结合;细胞成分:核质、核、细胞质;生物学过程:轴突束化，小脑发育，在子宫胚胎发育中(图6a)。

通过KEGG分析，我们标注了m6A甲基化水平显著改变的基因。这些基因主要在该途径中富集，包括肾细胞癌、癌症中的microRNAs和癌症中的蛋白多糖(图6b)。

3.6 TBI后基因表达谱

转录组的标准分析，比较容易复现，基本上看我六年前的表达芯片的公共数据库挖掘系列推文即可；

• 解读GEO数据存放规律及下载，一文就够
• 解读SRA数据库规律一文就够
• 从GEO数据库下载得到表达矩阵 一文就够
• GSEA分析一文就够（单机版+R语言版）
• 根据分组信息做差异分析- 这个一文不够的
• 差异分析得到的结果注释一文就够

结果发现，TBI后大鼠大脑皮层中428个mrna表达增加，280个mrna表达减少(p < 0.05, log2FC > 1)(图7b)。

Fig 7b, The expression level of mRNA post-TBI. A total of 428 mRNA expression increased and 280 mRNA expression decreased (p < 0.05, log2FC > 1) in rat cerebral cortex after TBI.

我们在火山图中标记出上调最明显的5个基因(Il1r2、Fosl1、Ptx3、Hsbp8和Msn)和下调最明显的5个基因(Hdac9、Tril、Fam107a、Pcdh20和Prom1)。

3.7 TBI后m6A甲基化与基因表达的联合分析

其中m6A峰值的mRNA有175个，其水平均发生显著变化，其中40个mRNA水平均上调，41个mRNA水平均下调。其中47个基因mRNA表达上调，m6A峰下调，47个基因mRNA表达下调，m6A峰上调。如图8。

3.8 FTO抑制剂FB23-2可以提高TBI大鼠mNSS评分

为了探讨m6A去甲基化酶FTO对TBI术后大鼠神经功能的控制作用，我们在TBI + FB23-2、TBI +DMSO、Sham+ FB23-2、Sham+DMSO四组进行mNSS和MWM实验。sham+ FB23-2组与sham+DMSO组连续4天mNSS评分无差异，说明FB23-2对正常大鼠神经功能无影响。所有TBI组（即TBI + FB23-2、TBI + DMSO组）mNSS得分显著高于假手术组（即Sham+ FB23-2, Sham+DMSO)，连续4天，提示CCI对神经功能造成损害。

Fig. 9 Effects of FTO inhibitior FB23–2 on neurological functions and cognitive and learning abilities in TBI.> a, the score of four groups in mNSS: TBI + FB23–2 (n = 7), TBI + DMSO (n = 8), Sham+FB23–2 (n = 5), and Sham+DMSO (n = 5). > > b and c, the escape latency time and target quadrant time of four groups in MWM, respectively: TBI + FB23–2 (n = 7), TBI + DMSO (n = 8), Sham+FB23–2 (n = 5), and Sham+DMSO (n = 5). “*”, “**”, and “***” indicate p < 0.05, p < 0.01, and p < 0.001 in TBI vs Sham, respectively; “#”, “##”, and “###” indicate p < 0.05, p < 0.01, and p < 0.001 in FB23–2 vs DMSO. Error bars represent SD. NS indicates p > 0.05

4 Conclusion

我们的研究表明，脑外伤后大鼠皮层的METTL14和FTO表达下调。通过MeRIPSeq，我们发现1580 mRNA中m6A甲基化水平发生了显著变化。将m6A峰与mRNA表达联合分析，发现有175个mRNA在TBI后发生显著变化。这些基因可能是干预脑外伤表观遗传调控的关键基因。此外，我们发现FTO在TBI大鼠的神经功能维护中发挥重要作用，可能是减轻TBI损伤的干预靶点。

5 总结

• 从结果中可以看出来。第一个结果查看了6个甲基化转移酶和去甲基化酶，**提出了两个基因METTL14和FTO的表达变化。**而后面的结果2，3，4，5，6，7却没有对这两个基因做任何的描述，都是对表达水平的变化和peak变化做了一个总的描述，并且这些变化中最显著的top20中也没有这两个基因。此外，挑选用来展示差异peak的三个基因是CD14，Dl14和Sox7。
• 后面将peaks和表达联合起来做分析，看peak上调表达下调和peak下调表达上调的基因中，只提到了各自有47个，却没有深入挖掘这47+47中的调控模式。
• 最后一个结果3.8显得非常生硬，没有前面铺叙后面沉淀。

6 扩展部分

6.1 TBI

traumatic brain injury，创伤性脑损伤。美国每年有530多万人因创伤性脑损伤(TBI)而致残，170万人遭受脑损伤，造成600多亿经济损失。创伤性脑损伤的原因包括爆炸伤、运动损伤和颅骨穿透。TBI的病理生理机制有直接机械力引起的原发性损伤和神经元的继发性损伤两方面。撞击产生的机械力直接损害轴突和神经元的一部分。继发性损伤的发生机制较为复杂。

6.2.甲基化机制：

• “writers” are methyltransferases：METTL3，METTL14 ；
• “readers” are m6A binding proteins involved in the translation process;
• “erasers” are demethylases：ALKBH5 and FTO。