MeRIP-seq 数据分析之 homer 富集 motif

点击上方关注“老俊俊”

1引言

在文章里我们经常会见到这样的图：

m6A 的序列特征一般是 RRACH(R = A or G, H = A, C, or U) ，GGACU/A 出现频率较多。callpeak 后，需要对得到的 peak 进行 motif 富集分析，验证是否含有经典的 m6A motif 特征。

同样我们可以使用 homer 进行分析，类似的还有 MEME。

2寻找 motif

准备输入数据

有个不错的教程推荐：Homer 软件的介绍-最全面而详细的找 motif 教程

我们既可以 提供 bed 文件 也可以 从 bed 文件提取对应的 fasta 序列 寻找 motif，这里以 bed 文件为例。

需要把 peak 文件准备成以下格式：

IVT_peak_1      chr1    4828519 4828619 .
IVT_peak_2      chr1    4832257 4832358 .
IVT_peak_3      chr1    4839303 4839415 .
IVT_peak_4      chr1    4867403 4867503 .

$ mkdir 9.motif-data && cd 9.motif-data/
$ awk '{print $4"\t"$1"\t"$2"\t"$3"\t"$6}' ../6.macs-data/M3_KO_ > IVT.macs2.bed
$ awk '{print $4"\t"$1"\t"$2"\t"$3"\t"$6}' ../6.macs-data/WT_macs2_intersect.bed > WT.macs2.bed
$ awk '{print $4"\t"$1"\t"$2"\t"$3"\t"$6}' ../6.macs-data/M3_KO_macs2_intersect.bed > M3_KO.macs2.bed

富集 motif

帮助文档：

(base) junjun@JunJunLab:/mnt/f/MeRIP-seq/9.motif-data$ findMotifs.pl

        Program will find de novo and known motifs in a gene list

                Usage:  findMotifs.pl <input list> <promoter set> <output directory> [additoinal options]

                example: findMotifs.pl genelist.txt mouse motifResults/ -len 10

                FASTA example: findMotifs.pl targets.fa fasta motifResults/ -fasta background.fa

        Available Promoter Sets: Add custom promoters sets with loadPromoters.pl
                mouse   mouse   /home/junjun/miniconda3/share/homer/.//data/promoters/  -2000   2000    refseq

                Try typing "perl /home/junjun/miniconda3/share/homer/.//configureHomer.pl -list" to see available promoter sets
                Typing "perl /home/junjun/miniconda3/share/homer/.//configureHomer.pl -install NNN" to install promoter set NNN
...

findMotifsGenome.pl

findMotifsGenome.pl 命令在基因组区域中寻找富集 Motifs，根据 bed 文件位置信息提取序列再富集：

# 批量操作
$ for i in IVT WT M3_KO;do findMotifsGenome.pl ${i}.macs2.bed mm10 ${i}/ -len 7 -S 15 -rna -p 10;done

参数说明：

• mm10： 物种基因组。
• -len： 指定 motif 长度。
• -S： 输出 motif 数量。
• -rna： 输出 RNA motif。
• -p： 多线程。

输出了 3 个文件夹：

生成一个汇总的网页文件，homerResults 文件夹则是具体每个 motif 的信息及 svg 图片：

汇总结果，看看 WT 的：

可以看到第二个 motif 就是 GGACU，第一个和第二个其实是 反向互补 的，还是很显著的。

M3KO 的：

IVT 的：

findMotifs.pl

findMotifs.pl 命令基因 fasta 序列寻找富集的 motif。

使用 exomepeak2 的 bed 文件来提取 fasta 序列，下载小鼠基因组文件：

$ cd ../index/
$ wget https://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/mm10.fa.gz
$ gunzip mm10.fa.gz

批量提取序列：

$ cd ../9.motif-data/
$ conda activate m6aseq
$ for i in IVT WT_reps_intersect M3KO_reps_intersect;do bedtools getfasta -fi ../index/mm10.fa -bed ../7.exomepeak2-data/${i}.bed -fo ./${i}.fa -s;done

motif 富集：

$ conda deactivate m6aseq
$ for i in IVT WT_reps_intersect M3KO_reps_intersect;do findMotifs.pl ${i}.fa fasta ${i}_exomepeak/ -len 7 -S 15 -rna -p 10;done

3提取 motif 绘制 logo

在 homerResults 含有每个 motif 的 motif 矩阵，svg 图片可以在 AI 里面修图，也可以提取 motif 矩阵重新绘制 Logo：

这里使用 ggseqlogo 来绘制：

# 安装
install.packages("ggseqlogo")

# 加载R包
require(ggplot2)
require(ggseqlogo)
library(tidyverse)

# 提取pvalue
row1 <- read.table('../9.motif-data/WT/homerResults/motif2.motif',nrows = 1,sep = '_')
pval <- strsplit(row1$V1,split = '\\,')[[1]][4]

# 读取motif矩阵
mo <- read.table('../9.motif-data/WT/homerResults/motif2.motif',skip = 1) %>% t(.)
rownames(mo) <- c('A','C','G','U')

# 绘图
ggseqlogo(mo,method = 'bits') +
  theme_bw(base_size = 20) +
  ggtitle(pval)

不错，效果还是可以的，和文章里差不多了。

如果有多个 motif 矩阵需要重新绘制，那就批量绘图，这里直接用 WT 样本的 12 个 motif 矩阵来绘图：

批量绘图：