使用 conda 安装 Homer
默认安装 homer 最新版
conda install homer -y # -y 表示直接安装，不询问确认与否
查看 HOMER 的数据库文件
which homer# 返回 homer 所在路径，我这里是~/miniconda3/envs/chipseq/bin/homercd打不开
Linux 查找软件安装路径的五种方法
-
find -name homer# 最推荐 -
locate homer# 相当于find -name,速度更快。不搜索具体目录，而是搜索一个数据库（/var/lib/locatedb） -
whereis homer#只能用于程序名的搜索，而且只搜索二进制文件（参数-b）、man说明文件（参数-m）和源代码文件（参数-s） -
which homer#在PATH变量指定的路径中，搜索某个系统命令的位置，并且返回第一个搜索结果 -
type -p homer# 相当于which homer
参考：HOMER安装和使用攻略，如何获取Motif?
我这里**在主目录**使用了 `find -name homer`
结果:
./miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer
./miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/cpp/backup/homer
./miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/cpp/homer
./miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/bin/homer
./miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/bin/homer
./miniconda3/envs/chipseq/bin/homer
./miniconda3/envs/chipseq/share/homer
./miniconda3/envs/chipseq/share/homer/bin/homer
./miniconda3/envs/chipseq/share/homer/cpp/backup/homer
./miniconda3/envs/chipseq/share/homer/cpp/homer
1 cd./miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer
2 perl configureHomer.pl -install hg38/mm10 # 下载参考数据库
3 cd ./data/genomes # 如果下载成功，可以看到有一个 hg38/mm10 的文件夹
4 du -sh hg38/mm10 # 查看文件夹大小，为5.7G/。
zhu:下载很慢
我是使用conda (conda 使用基础命令安装的homer，但是我不知道版本于是 查一下怎么看软件版本
1 conda --version #查看conda的版本
我现在的版本：==》conda 4.9.2
2 conda list homer #查看软件版本
Name Version Build Channel
homer 4.11 https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
HOMER isn't part of bioconda (partially due to its unorthodox code base/organization, although there is a HOMER package that may or may not work), but regardless, bioconda can help install the required packages and other NGS software even if HOMER itself must be installed manually.
homer官网说明
关于homer
- HOMER 基于 C++ 和 Perl 语言
- 用于 motif 查找和 数据分析的工具 需要两个序列作为参数：
- （1）参考序列：hg19、mm10 等基因组序列、promoter 序列、自定义的 FASTA 序列
- （2）所要分析的序列：DNA 或 RNA 序列
HOMER 适用于在大规模数据中寻找 DNA 或 RNA 序列的 motif
关于Motif以及实战
motif是比较有特征的短序列，会多次出现的，一般认为它的生物学意义重大，做完CHIP-seq分析之后，一般都会寻找motif 。查找有两种，一种是de novo的，要求的输入文件的fasta序列，一般是根据peak的区域的坐标提取好序列 。另一种是依赖于数据库的搜寻匹配，会将现有的ChIP-seq数据进行整合，提供更全面，更准确的motif数据库。
两个工具：
（1）Homer
（2）R包bioconductor
(3) MEME 或
参考文献A survey of motif finding Web tools for detecting binding site motifs in ChIP-Seq data
这里主要介绍使用HOMER寻找motif
1 提供包含基因组坐标的文件 比如 peak文件或BED文件
（1）分析peak文件中富集的motif：
findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size
简单例子：
findMotifsGenome.pl ERpeaks.txt hg18 ER_MotifOutput/ -size 200 -mask
-mask 使用repeated-mask序列
-size 设置motif长度
homerMotifs.motifs <#> : these are the output files from the de novo motif finding, separated by motif length, and represent separate runs of the algorithm.
homerMotifs.all.motifs : Simply the concatenated file composed of all the homerMotifs.motifs<#> files.
motifFindingParameters.txt : 记录执行findMotifsGenome.pl的命令
knownResults.txt : 记录motif的统计数据，text file(open in EXCEL).
seq.autonorm.tsv : autonormalization statistics for lower-order oligo normalization.
homerResults.html : de novo motif finding的格式化输出
实践
1 下载数据库，并锁定位置
find -name mm10
./miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/data/genomes/mm10
2 先找到自己处理好的bed文件位置
/home/yangjy/project/epi/align/mergeBam
2.* 发现在之前处理过后有两套_summits.bed文件，决定批量移动一下，只处理一种
使用mv命令：
2.1 查看当前目录
/home/yangjy/project/epi/align/mergeBam/
2.*2 查看要移动到的目录
/home/yangjy/project/epi/peaks
先cd到文件当前所在目录
使用mv *.log_summits.bed /home/yangjy/project/epi/peaks
3.0 所以处理好的bed文件所在目录为：
/home/yangjy/project/epi/peaks
3 来到处理目录下
cd ~/project/epi/motif
3.1 写批量脚本，先echo一下
for id in /home/yangjy/project/epi/peaks/*.bed;
do
echo $id
done
成功的话：出现：
/home/yangjy/project/epi/peaks/Control.log_summits.bed
/home/yangjy/project/epi/peaks/H2Aub1.log_summits.bed
/home/yangjy/project/epi/peaks/H3K36me3.log_summits.bed
/home/yangjy/project/epi/peaks/mm10.tss.bed
/home/yangjy/project/epi/peaks/Ring1B.log_summits.bed
/home/yangjy/project/epi/peaks/RNAPII_8WG16.log_summits.bed
/home/yangjy/project/epi/peaks/RNAPII_S2P.log_summits.bed
/home/yangjy/project/epi/peaks/RNAPII_S5P.log_summits.bed
/home/yangjy/project/epi/peaks/RNAPII_S7P.log_summits.bed
好的，下一步：
加上找文件名的代码
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
于是：==》
在
for id in /home/yangjy/project/epi/peaks/*.bed;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
done
有：
/home/yangjy/project/epi/peaks/Control.log_summits.bed
Control
/home/yangjy/project/epi/peaks/H2Aub1.log_summits.bed
H2Aub1
/home/yangjy/project/epi/peaks/H3K36me3.log_summits.bed
H3K36me3
/home/yangjy/project/epi/peaks/mm10.tss.bed
mm10
/home/yangjy/project/epi/peaks/Ring1B.log_summits.bed
Ring1B
/home/yangjy/project/epi/peaks/RNAPII_8WG16.log_summits.bed
RNAPII_8WG16
/home/yangjy/project/epi/peaks/RNAPII_S2P.log_summits.bed
RNAPII_S2P
/home/yangjy/project/epi/peaks/RNAPII_S5P.log_summits.bed
RNAPII_S5P
/home/yangjy/project/epi/peaks/RNAPII_S7P.log_summits.bed
RNAPII_S7P
继续 对bed文件进行处理
首先我们要取重要的信息，第4,1,2,3列的信息，使用awk,print命令：
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >homer_peaks.tmp
使用less命令打开任意一个文档看一下
less /home/yangjy/project/epi/peaks/H2Aub1.log_summits.bed
看到：
chr1 4492667 4492668 H2Aub1.log_peak_1 4.02077
chr1 4493157 4493158 H2Aub1.log_peak_2 8.94022
chr1 4493467 4493468 H2Aub1.log_peak_3 4.02077
chr1 4496553 4496554 H2Aub1.log_peak_4 2.72721
chr1 4497102 4497103 H2Aub1.log_peak_5 4.02077
chr1 4497409 4497410 H2Aub1.log_peak_6 4.02077
chr1 5916554 5916555 H2Aub1.log_peak_7 5.07861
chr1 5917032 5917033 H2Aub1.log_peak_8 5.96592
chr1 9546061 9546062 H2Aub1.log_peak_9 1.92244
第4,1,2,3列的信息应该分别表示,eg：H2Aub1.log_peak_1 chr1 4492667 4492668
验证一下：
for id in /home/yangjy/project/epi/peaks/*_summits.bed;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >homer_peaks.tmp
done
查看命令：head homer_peaks.tmp
可以看到Homer Peak/Positions file有4列：
第一列: 染色体
第二列: 起始位置
第三列: 终止位置
第四列: 链的方向(+/- or 0/1, where 0="+", 1="-")
RNAPII_S7P.log_peak_1 chr1 4085494 4085495 +
RNAPII_S7P.log_peak_2 chr1 4559682 4559683 +
RNAPII_S7P.log_peak_3 chr1 4559912 4559913 +
RNAPII_S7P.log_peak_4 chr1 4560134 4560135 +
RNAPII_S7P.log_peak_5 chr1 4560750 4560751 +
RNAPII_S7P.log_peak_6 chr1 4769884 4769885 +
RNAPII_S7P.log_peak_7 chr1 4770731 4770732 +
RNAPII_S7P.log_peak_8 chr1 4771139 4771140 +
RNAPII_S7P.log_peak_9 chr1 4771620 4771621 +
RNAPII_S7P.log_peak_10 chr1 4772002 4772003 +
为什么要这样，因为Homer软件的说明书要求这样的格式；才能给做注解；
接下来开始做注释：
为什么给了mm10就可以自动做注释；
下载数据之前先查看数据列表，看Homer已经有哪些数据包：
语法公式：
perl /path-to-homer/configureHomer.pl -list
具体事例：perl /home/yangjy/miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/configureHomer.pl -list
等同于：perl /home/yangjy/miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/.//configureHomer.pl -list
主要是给系统一个configureHomer.pl路径
save to :/home/yangjy/miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer//update.txt
启用安装：
同前面描述，不介意再来一遍：（Homer数据包安装使用-install参数）
perl /path-to-homer/configureHomer.pl -install mouse#下载老鼠启动子数据
perl /path-to-homer/configureHomer.pl -install mm10#下载 mm10小鼠参考基因组
perl /path-to-homer/configureHomer.pl -install hg19#下载 hg19人的参考基因组
实例：perl /home/yangjy/miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/configureHomer.pl -install mm10
很慢，可以挂载后台：
nohupperl /home/yangjy/miniconda3/pkgs/homer-4.11-pl526hc9558a2_3/share/homer/configureHomer.pl -install mm10&
综上：
for id in /home/yangjy/project/epi/peaks/*_summits.bed;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >homer_peaks.tmp
findMotifsGenome.pl homer_peaks.tmp mm10 ${sample}_motifDir -len 8,10,12
annotatePeaks.pl homer_peaks.tmp mm10 1>${sample}.peakAnn.xls 2>${sample}.annLog.txt
done
==》输出文件
其实control的bed为0 没有必要做注释；
下载看看文件如何：pwd---/home/yangjy/project/epi/motif
cd motifDir;会看到以下文件
homerMotifs.all.motifs
homerMotifs.motifs10
homerMotifs.motifs12
homerMotifs.motifs8
homerResults
homerResults.html
knownResults
knownResults.html
knownResults.txt
motifFindingParameters.txt
seq.autonorm.tsv
打开后缀为.html的文件，可以看到
单独分析 寻找富集Motifs
findMotifsGenome.pl命令用于在基因组区域中寻找富集Motifs
命令使用格式为：
findMotifsGenome.pl <Homer Peak/Positions file> <genome> <output directory> -size # [options]
findMotifsGenome.pl homer_peaks.tmp mm10 ${sample}_motifDir -len 8,10,12
annotatePeaks.pl homer_peaks.tmp mm10 1>${sample}.peakAnn.xls 2>${sample}.annLog.txt
#参数解释
- -输入文件：awk处理好的Homer Peak/Positions file
- -参考基因组：这里是mm10
- -输出文件：给一个路径和输出文件的名字
- -len：motif大小设置，默认8,10,12；越大需要的计算资源越多
参考：
- http://homer.ucsd.edu/homer/motif/index.html
- https://mp.weixin.qq.com/s?src=11×tamp=1611902775&ver=2857&signature=tPyh2rni2VaBUUqMK-0wdP0iC8yTxY0zTQgXokd7pARHxx8NrzB1b2jJntEEvL3UzXWoRABKLn2h2gHUU1yQ0WLc8hvYe0Ow196f75LoKvWmf5cpzRKU6urXtC9ztPAf&new=1
补充
通过meme去来找motif，需要bed格式的peaks的坐标来获取fasta序列。
MEME，链接：http://meme-suite.org/
