MeRIP-seq 数据分析之 homer 注释 peaks

原创 JunJunLab 老俊俊的生信笔记 2022-08-15

收录于合集

点击上方关注老俊俊

1引言

callpeak 步骤得到的是每个 m6A 位点的 位置信息，包括染色体、起始位置、终止位置、正负链 等。我们具体想知道该位点位于哪个基因上，属于哪一个位置，基因属于哪种类型等信息就需要 对 peaks 进行注释，exomePeak2 的结果文件会有基因名注释，macs2 结果就没有了。

可以做这个事情的软件由余老师的 Chipseeker、ChIPpeakAnno 等包，今天介绍一下怎么使用 homer 进行 peak 注释。

2安装

使用 conda 安装：

$ conda install -c bioconda homer

下载配置文件，把内容复制保存到 configureHomer.pl 里即可，一个 perl 脚本：

地址：
https://github.com/IGBIllinois/HOMER/blob/master/configureHomer.pl

下载参考基因组，手机开热点贼快：

#下载 mm10小鼠参考基因组
$ perl /path-to-homer/configureHomer.pl -install mm10
#下载 hg19人的参考基因组
$ perl /path-to-homer/configureHomer.pl -install hg19

3注释 peaks

$ cd /mnt/f/MeRIP-seq/
$ mkdir 8.annotate-peaks-data

## annotate macs2 data
$ annotatePeaks.pl 6.macs-data/IVT_peaks.narrowPeak mm10 > 8.annotate-peaks-data/IVT_peaks.annotated.tx
$ annotatePeaks.pl 6.macs-data/WT_macs2_intersect.bed mm10 > 8.annotate-peaks-data/WT_macs2_intersect.annotated.txt
$ annotatePeaks.pl 6.macs-data/M3_KO_macs2_intersect.bed mm10 > 8.annotate-peaks-data/M3KO_macs2_intersect.annotated.txt

## annotate exomePeak2 data
$ annotatePeaks.pl 7.exomepeak2-data/IVT.bed mm10 > 8.annotate-peaks-data/IVT.annotated.txt
$ annotatePeaks.pl 7.exomepeak2-data/WT_reps_intersect.bed mm10 > 8.annotate-peaks-data/WT_reps_intersect.annotated.txt
$ annotatePeaks.pl 7.exomepeak2-data/M3KO_reps_intersect.bed mm10 > 8.annotate-peaks-data/M3KO_reps_intersect.annotated.txt

注释也是比较全面的：

主要由 注释位置、Entrez ID、Esembl ID、gene name、基因描述等：

4绘图

我们根据 Annotation 列的分类进行绘图，读取 exomepeak2 的注释结果。

新建 R project，创建 step5_annotate_peaks.R 脚本：

# 设置工作路径
setwd('../8.annotate-peaks-data/')

# 加载R包
library(tidyverse)
library(ggplot2)
library(reshape2)

# 读取注释文件
IVT <- read.delim('IVT.annotated.txt',header = T)
WT <- read.delim('WT_reps_intersect.annotated.txt',header = T)
M3ko <- read.delim('M3KO_reps_intersect.annotated.txt',header = T)

# 保存到list
df <- list(IVT = IVT$Annotation,WT = WT$Annotation,M3ko = M3ko$Annotation)

# 切分保存
lapply(1:3, function(x){
  tmp = sapply(strsplit(df[[x]],split = '\\('),'[',1)
  res = table(tmp)
}) %>% do.call('rbind',.) %>% as.data.frame() -> type_comb

# 添加名称
type_comb$name = c('IVT','WT','M3ko')

# 宽数据转长数据
final <- melt(type_comb)

# 因子化
final$name <- factor(final$name,levels = c('IVT','WT','M3ko'))

# 绘制饼图
ggplot(final,aes(x = '',y = value,fill = variable)) +
  geom_col(position = position_fill()) +
  theme_void() +
  theme(legend.position = 'bottom',
        strip.text.x = element_text(size= 14)) +
  facet_wrap(~name) +
  coord_polar(theta = 'y') +
  scale_fill_brewer(palette = 'Set1',name = 'Region types')

看到注释在 intron 里是最多的，其次是 exon 和 3'UTR。可以看到 METTL3 KO 后，exon 和 3'UTR 的 peak 明显减少了。

我们把输入文件换成 macs2 的，后面代码不变走一遍：

# 读取注释文件 macs2
IVT <- read.delim('IVT_peaks.annotated.txt',header = T)
WT <- read.delim('WT_macs2_intersect.annotated.txt',header = T)
M3ko <- read.delim('M3KO_macs2_intersect.annotated.txt',header = T)

...省略

啥玩意，可以看到 exon 所占比例是最多的，其次是 intron，METTL3 KO 后，3'UTR 的 peak 明显减少了，exon 上感觉差的不大。

这样对比，感觉 exomepeak2 callpeak 的区域在 intron 会多很多，相对于 macs2 的话，具体什么原因和算法及模型有关，还是其它原因？

绘制条形堆积图展示：

# 条形堆积图
# 宽数据转长数据
final_macs2 <- melt(type_comb)
final_exomepeak2 <- melt(type_comb)

# 添加分组
res <- rbind(final_macs2,final_exomepeak2)
res$sample <- c(rep('Macs2',24),rep('ExomePeak2',24))

# 绘图
ggplot(res,aes(x = name,y = value,fill = variable)) +
  geom_col(position = position_fill()) +
  theme_bw(base_size = 16) +
  scale_y_continuous(labels = scales::percent_format()) +
  theme(legend.position = 'right',
        strip.text.x = element_text(size= 14)) +
  facet_wrap(~sample) +
  scale_fill_brewer(palette = 'Set1',name = 'Region types') +
  xlab('') + ylab('Percent')