m6A-seq 数据处理及图表复现交流群

前面我们分享了 跟着Nature Medicine学MeDIP-seq数据分析，数据和代码都是公开，这个2G的压缩包文件，足以学习3个月，写60篇教程。同时也分享了 全套MeRIP-seq文章图表复现代码，其实MeRIP-seq其实就是RNA水平的又叫做m6a测序。

但很多粉丝留言表示这些英文教程看不懂，数据也很分散，没有中文解说实在是很难跟下来，希望我们出一个手把手系列教程。

这个全套 MeRIP-seq 图表复现代码在GitHub：https://github.com/al-mcintyre/merip_reanalysis_scripts 这个也是接近2G的压缩包！

其实很早以前我就在《生信技能树》发布过教程：新的ngs流程该如何学习（以CUT&Tag 数据处理为例子），提到了我自己是不太可能去把所有的ngs流程全部录制视频的，只能说是更好的传达学习方法给到大家。其实如果你看过我表观组学，比如《ChIP-seq数据分析》 和 《ATAC-seq数据分析》 就会发现其实这个m6A数据处理大同小异的，当然了，肯定是会有一些细微差异是需要注意的。

虽然我没有时间，但是我们的两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，还是有幸招募到了愿意花时间给大家做整理的小伙伴。所以有了这个MeRIP-seq 图表复现交流群。

01-文献解读|在人的四种不同组织中m6A甲基化的遗传驱动因素

1.文献信息

标题：Genetic drivers of m6A methylation in human brain, lung, heart and muscle 发表时间：2021年7月1日 杂志：Nature Genetics（2020IF=27.605） 作者：麻省理工学院计算机与人工智能实验室Manolis Kellis课题组

数据和方法

主要数据

image-20210711235626527

可以看到，样本量还是蛮可观的！

m6A meRIP–seq across human tissues

文章主要使用的链特异性建库和双端45bp的测序策略，第一次遇到测这么短的啊：

We used the SMARTer Stranded RNA-Seq Kit from Clontech/Takada, which is optimized to work with 100pg of starting material. The libraries were sent for 2×45-base-pair paired-end sequencing.

甲基化数据提到的数据库为eGTEx，可以参考：https://cloud.tencent.com/developer/news/397979 进行此数据库的了解。

分析流程：

• 1.去除tRNA和rNRA：使用bowtie2将read比对到tRNA和rRNA，保留 the unmapped reads
  • tRNA：downloaded from the University of California, Santa Cruz Table Browser
  • rRNA：downloaded from the National Center for Biotechnology Information Nucleotide database
• 2.比对到参考基因组：使用hisat2 比对到the hg38 human genome：GENCODE (v26; downloaded from https://gtexportal.org/home/datasets
• Peak calling：MACS2

我们后续就按照这个流程来进行图表复现！

主要结果

这是一篇研究quantitative trait loci（QTL）与m6A修饰关联的文章，我们这次主要关注m6A的地方，即结果1。

结果1：m6A variation across tissues and individuals

107个病人共176个样本，在QC之后，剩余91个病人129个样本做分析：53 brain, 12 heart, 32 muscle

and 32 lung samples，见下图A。

Peak calling共得到>278,000个peak位点（每个位点至少在两个病人中存在），平均每个样本约20,000个位点。

这些Peaks中，与以前的结果相比，有很多都是以前未检测到的peak。比较的对象是来自数据库RMBase v2.0中发表的peaks。

RMBase v2.0: deciphering the map of RNA modifcations from epitranscriptome sequencing data

并且，检测到的Peaks位点的主要集中在终止密码子附近。

作者还将以前没有检测到的m6A位点序列特征分布也绘制了一个图，如下，这个图5’UTR位置的信息要比终止密码子的信号值高：

image-20210712011454191

保守序列特征：GGACH

image-20210712011135022

利用m6A谱对样本进行相似性分析，发现组织类型是主要的差异来源，与基于RNA表达谱的样本相似性分析结果类似。

基于m6A谱的样本相似性聚类：

基于RNA表达谱的样本相似性聚类：

具有组织特异m6A修饰位点的基因的功能富集结果显示，组织功能相关功能富集：

此外，作者还对具有组织特异性m6A修饰位点的基因进行了展示：

这些具有m6A修饰位点信息的基因可以在数据库进行查询展示：https://www.gtexportal.org/home/browseEqtls?location=chr1%3A750000-850,000

• 大脑特异m6A修饰的POU3F2基因：
• 肺特异的m6A修饰EGFR 基因

最后，不同组织总都有m6A修饰的基因在不同组织中的表达不差异，具有m6A修饰组织特异的表达 也特异。这表明m6A可以导致其他广泛表达的转录本的组织特异性功能。

其他内容结果如下，详细版本可以前往看文献：

结果2：m6A genetic driver discovery and validation

结果3：Tissue specificity of m6A QTLs.

结果4：m6A QTLs and eQTLs sometimes overlap but are mostly independent

结果5：m6A QTLs help interpret GWAS loci

结果6：Tissue-specific m6A-QTL enrichments of GWAS variants

结果7：Novel m6A regulator prediction

非常好的资源：

此篇文章贡献了作者全部的分析代码，如下：

下载链接

http://compbio.mit.edu/m6AQTLs/

可惜的是下载不到数据啊，即使是m6A的bed文件也没有公开，只能在GTEx中进行相关位点信息查询。

如有万能的网友能搞到数据，我就可以给你复现整个文章了。

不过，我还有plan B：准备找一套其他的数据，利用此篇文献的所有代码进行数据分析。后续即将更新~