10个图，教你搞定国自然热点分析！专家都拒绝不了的美图！
探索生信之美，解构每段代码背后的故事
大家好呀，我是风间琉璃，前面两期我们分别介绍了如何使用exomePeak和MACSr两种方法对MeRIP-Seq进行peak鉴定以及定量分析。今天我们将介绍大名鼎鼎的生信大咖Y叔在2015年推出的R包——ChIPseeker用于m6A peak进行注释。之后小伙伴使用ChIPseeker发表文章一定要引用相应的参考文献喔：“G Yu, LG Wang, QY He. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 2015, 31(14):2382-2383. doi:[10.1093/bioinformatics/btv145](http://dx.doi.org/10.1093/bioinformatics/btv145)”。大家可以看到因为同样是根据免疫沉淀原理所衍生出来的技术，所以MeRIP-Seq和ChIP-Seq在分析过程有些步骤可以互相借鉴。
这一期我们以一篇2020年发表在《Molecular Brain》杂志上的案例文献“Epitranscriptomic profiling of N6-methyladenosine-related RNA methylation in rat cerebral cortex following traumatic brain injury”进行参考学习，看看这篇文章怎么应用Chipseeker进行注释的吧。
(最新2021年影响因子已在挑圈联靠公号后台更新，
回复期刊名称，得最新影响因子信息）
研究揭示了在发生创伤性恼损伤的大鼠模型中，METTL14和FTO表达水平在明显下降并揭示了FTO在创伤性脑损伤中的潜在调控机制，在接下来的MeRIP-Seq/RNA-Seq分析中，作者筛选明显发生m6A甲基化水平改变以及转录水平改变的基因，并筛选出里面的核心基因。
在文章的方法学部分，作者阐明使用MeTDiff软件进行peak鉴定以及差异分析。对于Peak的注释，作者则使用的是ChIPseeker包。
好啦，介绍完背景和示例文献之后我们就开始本期的分析吧~
首先我们知道ChIPseeker能够对Peak进行注释以及对peak的覆盖区域进行可视化，但是需要我们输入peak鉴定分析后的结果。格式可以是txt、xls也可以是bed格式，我们以上两期使用MACSr分析得到的m6A_peak.xls文件或者exomePeak分析得到的DiffMod.bed为例。
# 安装包 if (!requireNamespace('BiocManager', quietly = TRUE))# install.packages('BiocManager') BiocManager::install('ChIPseeker') 加载包library(ChIPseeker)library(TxDb.Hsapiens.UCSC.hg19.knownGene)txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene# 读取数据，以exomePeak为例ex_peak <- readPeakFile("exomePeak/DiffMod.bed")ex_peak## GRanges object with 17 ranges and 9 metadata columns:## seqnames ranges strand | V4 V5 V6## <Rle> <IRanges> <Rle> | <character> <numeric> <character>## [1] chr7 150742737-150742836 * | peak_1 9.3551760 +## [2] chr7 150743012-150743111 * | peak_2 0.0880848 +## [3] chr7 150743287-150743661 * | peak_3 1.6419373 +## [4] chr7 150743687-150743761 * | peak_4 3.7446046 +## [5] chr7 150743787-150744061 * | peak_5 3.7446046 +## ... ... ... ... . ... ... ...## [13] chr7 150932208-150933602 * | peak_13 12.98192 +## [14] chr7 150934561-150935735 * | peak_14 9.35518 +## [15] chr7 150776938-150777949 * | peak_15 1.90953 -## [16] chr7 151038922-151042472 * | peak_16 2.31436 +## [17] chr7 150778422-150778721 * | peak_17 1.64194 -## V7 V8 V9 V10 V11 V12## <integer> <integer> <integer> <integer> <character> <character>## [1] 150742736 150742836 0 1 100 0## [2] 150743011 150743111 0 1 100 0## [3] 150743286 150743661 0 1 375 0## [4] 150743686 150743761 0 1 75 0## [5] 150743786 150744061 0 1 275 0## ... ... ... ... ... ... ...## [13] 150932207 150933602 0 2 491,109 0,1286## [14] 150934560 150935735 0 1 1175 0## [15] 150776937 150777949 0 2 72,178 0,834## [16] 151038921 151042472 0 2 11,39 0,3512## [17] 150778421 150778721 0 1 300 0## -------## seqinfo: 1 sequence from an unspecified genome; no seqlengths# 读取数据，以MACSr分析数据为例ma_peak <- readPeakFile("MeRIP/m6A_peaks.xls")ma_peak## GRanges object with 37 ranges and 7 metadata columns:## seqnames ranges strand | length abs_summit pileup## <Rle> <IRanges> <Rle> | <integer> <integer> <integer>## [1] chr7 2-150669467 * | 150669467 150668295 0## [2] chr7 150674833-150675096 * | 265 150674957 3## [3] chr7 150675171-150675319 * | 150 150675245 1## [4] chr7 150677251-150682088 * | 4839 150678657 1## [5] chr7 150683090-150699782 * | 16694 150697589 0## ... ... ... ... . ... ... ...## [33] chr7 150958592-150958740 * | 150 150958666 1## [34] chr7 150977763-151009780 * | 32019 151007778 0## [35] chr7 151011277-151033932 * | 22657 151033931 0## [36] chr7 151038886-151039091 * | 207 151038947 18## [37] chr7 151042309-151042466 * | 159 151042419 13## X.log10.pvalue. fold_enrichment X.log10.qvalue. name## <numeric> <numeric> <numeric> <character>## [1] 1.35752 0.957035 0.000000 m6A_peak_1## [2] 4.73180 3.479330 0.360032 m6A_peak_2## [3] 2.57455 1.860770 0.000000 m6A_peak_3## [4] 3.00963 1.914070 0.000000 m6A_peak_4## [5] 1.35752 0.957035 0.000000 m6A_peak_5## ... ... ... ... ...## [33] 2.29090 1.810360 0.00000 m6A_peak_33## [34] 1.53039 0.970940 0.00000 m6A_peak_34## [35] 1.35752 0.957035 0.00000 m6A_peak_35## [36] 9.26928 4.758400 4.73629 m6A_peak_36## [37] 11.08170 6.837510 6.49158 m6A_peak_37## -------## seqinfo: 1 sequence from an unspecified genome; no seqlengths发现都可以顺畅的读取，非常丝滑~需要解释一下，exomePeak得到的是差异分析之后的数据，而MACSr则是进行鉴定peak的数据，两者并不相同。但是ChIPseeker都能够轻松胜任，并进行后续的分析。接下来我们就直接开始可视化了吧。
通过上游对peak进行鉴定后，我们可以通过在基因组层面上对其进行可视化，展示在哪条染色体上的哪个区间出现peak,peak的大小如何。
covplot(ma_peak, weightCol = "pileup")
接下来我们可以展示每个peak在所有基因启动子区域附近的分布情况。
# 获得所有基因TSS上下游3000bp的区域promoter <- getPromoters(TxDb = txdb, upstream = 3000, downstream = 3000)# 获得peak文件在所有基因上TSS上游的分布情况tagMatrix <- getTagMatrix(ma_peak, windows = promoter)# 绘制热图tagHeatmap(tagMatrix, xlim = c(-3000, 3000), color = "red")
这里面每一行代表一个基因，展示peak在所有基因TSS两侧的分布。除此之外，我们还可以用密度图的方式展示peak文件在所有基因启动子附近信号分布的强度图。
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),#TSS上下游区间 xlab="Genomic Region (5'->3')", #X坐标题目 ylab = "Read Count Frequency") #Y轴坐标题目
#通过重抽样的方式还可以添加置信区间plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)## >> Running bootstrapping for tag matrix... 2021-06-22 23:39:21
对于m6A-Seq分析，我们找到相应的peak之后，我们需要做得是明确这些peak属于哪些mRNA，换句话说，是哪些基因编码的mRNA具有m6A甲基化修饰从而影响基因发挥的功能，这是我们研究者最关心的问题，而ChIPseeker能够对我们鉴定出来的peak进行注释分析。那我们看看吧~
peakAnno <- annotatePeak(ma_peak, tssRegion = c(-3000, 3000), TxDb = txdb, annoDb = "org.Hs.eg.db")## >> preparing features information... 2021-06-22 23:39:22 ## >> identifying nearest features... 2021-06-22 23:39:23 ## >> calculating distance from peak to TSS... 2021-06-22 23:39:23 ## >> assigning genomic annotation... 2021-06-22 23:39:23 ## >> adding gene annotation... 2021-06-22 23:39:33 ## >> assigning chromosome lengths 2021-06-22 23:39:33 ## >> done... 2021-06-22 23:39:33# 查看一下peakAnno## Annotated peaks generated by ChIPseeker## 37/37 peaks were annotated## Genomic Annotation Summary:## Feature Frequency## 5 Promoter (<=1kb) 43.243243## 6 Promoter (1-2kb) 16.216216## 7 Promoter (2-3kb) 10.810811## 2 3' UTR 8.108108## 1 1st Intron 5.405405## 4 Other Intron 10.810811## 3 Distal Intergenic 5.405405# 我们可以将peakAnno的csAnno格式改为GRanges格式或者数据框格式G_Anno <- as.GRanges(peakAnno)G_Anno[, 1:6]## GRanges object with 37 ranges and 6 metadata columns:## seqnames ranges strand | length abs_summit pileup## <Rle> <IRanges> <Rle> | <integer> <integer> <integer>## [1] chr7 2-150669467 * | 150669467 150668295 0## [2] chr7 150674833-150675096 * | 265 150674957 3## [3] chr7 150675171-150675319 * | 150 150675245 1## [4] chr7 150677251-150682088 * | 4839 150678657 1## [5] chr7 150683090-150699782 * | 16694 150697589 0## ... ... ... ... . ... ... ...## [33] chr7 150958592-150958740 * | 150 150958666 1## [34] chr7 150977763-151009780 * | 32019 151007778 0## [35] chr7 151011277-151033932 * | 22657 151033931 0## [36] chr7 151038886-151039091 * | 207 151038947 18## [37] chr7 151042309-151042466 * | 159 151042419 13## X.log10.pvalue. fold_enrichment X.log10.qvalue.## <numeric> <numeric> <numeric>## [1] 1.35752 0.957035 0.000000## [2] 4.73180 3.479330 0.360032## [3] 2.57455 1.860770 0.000000## [4] 3.00963 1.914070 0.000000## [5] 1.35752 0.957035 0.000000## ... ... ... ...## [33] 2.29090 1.810360 0.00000## [34] 1.53039 0.970940 0.00000## [35] 1.35752 0.957035 0.00000## [36] 9.26928 4.758400 4.73629## [37] 11.08170 6.837510 6.49158## -------## seqinfo: 1 sequence from hg19 genome# 数据框格式df_Anno <- as.data.frame(peakAnno)df_Anno[, 1:6]## seqnames start end width strand length## 1 chr7 2 150669467 150669466 * 150669467## 2 chr7 150674833 150675096 264 * 265## 3 chr7 150675171 150675319 149 * 150## 4 chr7 150677251 150682088 4838 * 4839## 5 chr7 150683090 150699782 16693 * 16694## 6 chr7 150700897 150718742 17846 * 17847## 7 chr7 150719744 150720587 844 * 845## 8 chr7 150725488 150725701 214 * 215## 9 chr7 150742495 150744028 1534 * 1535## 10 chr7 150744524 150744672 149 * 150## 11 chr7 150746351 150746644 294 * 295## 12 chr7 150747184 150747383 200 * 201## 13 chr7 150778261 150778725 465 * 466## 14 chr7 150801595 150802188 594 * 595## 15 chr7 150805659 150805807 149 * 150## 16 chr7 150819852 150820124 273 * 274## 17 chr7 150820399 150820584 186 * 187## 18 chr7 150828076 150828315 240 * 241## 19 chr7 150840914 150841435 522 * 523## 20 chr7 150847522 150900137 52616 * 52617## 21 chr7 150909215 150909363 149 * 150## 22 chr7 150909583 150909981 399 * 400## 23 chr7 150910276 150910610 335 * 336## 24 chr7 150910827 150911183 357 * 358## 25 chr7 150924110 150924258 149 * 150## 26 chr7 150929542 150929720 179 * 180## 27 chr7 150930865 150931152 288 * 289## 28 chr7 150932097 150932792 696 * 697## 29 chr7 150934441 150935812 1372 * 1373## 30 chr7 150945902 150946271 370 * 371## 31 chr7 150953850 150953998 149 * 150## 32 chr7 150954848 150954996 149 * 150## 33 chr7 150958592 150958740 149 * 150## 34 chr7 150977763 151009780 32018 * 32019## 35 chr7 151011277 151033932 22656 * 22657## 36 chr7 151038886 151039091 206 * 207## 37 chr7 151042309 151042466 158 * 159通过这个我们得到很多有用的信息：1.peak距离最近的基因是哪些；2.与最近基因的TSS距离；3.位于基因组的哪个区域（promoter、5’UTR、3‘UTR等等）。接下来，我们还可以对注释进行可视化。# 绘制饼图plotAnnoPie(peakAnno)
# 绘制条形图plotAnnoBar(peakAnno)
# 由于不同基因组的注释存在重叠，比如外显子和启动子# 多以我们可以用不同层次的饼图进一步可视化vennpie(peakAnno)# 也可以用UpsetRupsetplot(peakAnno)# 还可以将vennpie以及upsetR合并在一起绘制upsetplot(peakAnno, vennpie = TRUE)
另外我们同样可以根据peak到TSS的距离，通过条图的方式进行可视化。
plotDistToTSS(peakAnno, title = "Distribution of peak relative to TSS")
最后我们还可以对peak修饰的mRNA的功能进行富集分析，明确m6A修饰可能会影响的下游通路与细胞功能。
library(ReactomePA)pathway <- enrichPathway(as.data.frame(peakAnno)$geneId)head(pathway, 2)## [1] ID Description GeneRatio BgRatio pvalue p.adjust ## [7] qvalue geneID Count ## <0 行> (或0-长度的row.names)# 开始绘制 dotplot(pathway)
ChIPseeker一个包能够提供这么多可视化的图片。不愧是Y叔的R包。好啦，本次的分享到这里就结束了，MeRIP-Seq和ChIP-Seq在分析流程中有许多相似的步骤。我们上面讲到很多的可视化图片，这在ChIP-Seq分析中非常实用，对于分析m6A-Seq分析的同学来说，最需要注意的是peak注释，明确peak映射到哪个基因，以及相应的功能注释。好啦，我是风间琉璃，咱们下期见~
往期传送门
—
END—
撰文丨风间琉璃
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
关键词
数据
基因
文件
区域
注释
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明：以上内容为用户推荐收藏至CareerEngine平台，其内容（含文字、图片、视频、音频等）及知识版权均属用户或用户转发自的第三方网站，如涉嫌侵权，请通知[email protected]进行信息删除。如需查看信息来源，请点击“查看原文”。如需洽谈其它事宜，请联系[email protected]。
版权声明：以上内容为用户推荐收藏至CareerEngine平台，其内容（含文字、图片、视频、音频等）及知识版权均属用户或用户转发自的第三方网站，如涉嫌侵权，请通知[email protected]进行信息删除。如需查看信息来源，请点击“查看原文”。如需洽谈其它事宜，请联系[email protected]。