关于查找motif，目前有很多种软件可以进行预测。我所在的实验室通常使用FIMO（MEME套件里的一个），但是有很多文献里也提到了HOMER这个软件，并且不乏一些影响因子很高的paper，所以感觉有必要学习一下这个软件。这个官网的信息很多，写的也很详细，但是感觉思维比较混乱，所以这篇学习笔记是我自己整理的官网信息，顺序和官网不一样，有兴趣的同学可以直接去官网学习：homer官方网站：http://homer.ucsd.edu/homer/motif/
HOMER有很多功能，这里我只集中学习用HOMER查找motif。

（一）HOMER是什么？

HOMER(Motif富集的超几何优化)是一套用于Motif发现和ChIP-Seq分析的工具。它是一个命令行程序的集合，用于unix操作系统，主要用perl和c++编写。HOMER最初是作为一个从头（de novo）motif发现算法编写的，非常适合在大规模基因组数据中发现8-12 bp的motif。

硬件要求（括号里是推荐的要求）： 2+ Gb 内存 (4-8+ Gb), 10+ Gb 硬盘空间 (50+ Gb)
软件要求：Unix compatible OS (or cygwin), perl, gcc, make, wget, ghostscript, weblogo, blat

HOMER包含了一个新颖的motif发现算法，设计用于在基因组学应用里的调控元件分析(这里的元件指DNA，而非蛋白质)。这是一种不同的motif发现算法，这意味着它采取两套序列，并试图识别在其中一套序列上相对于另一套的特定富集的调控元件（比如说我们的peaks序列相对于整个基因组序列）。它使用ZOOPS评分(0或1发生在每条序列上)与超几何富集计算(或二项式)来确定motif富集。HOMER也尽力考虑数据集里的排序偏差。它的设计用于ChIP-Seq和启动子分析，但可以应用于几乎任何核酸序列的motif发现。

（二）HOMER的安装

安装homer（下面的安装方法来自：生信技能树视频chipseq 17 使用homer找motif）因为我使用的是服务器上已经安装好的软件，所以这一步需要大家自己摸索一下~:

#install homer
conda install -c bioconda homer

下载configureHomoer.pl:here，并把这个文件放在一个单独的文件夹里，比如我用的就是名为homer的文件夹。

（三）HOMER基因组准备

HOMER使用的基因组和我们平时用的基因组不一样，要预处理一下的，可以直接下载你需要的信息，比如：

#这里我需要hg38的基因组
#同样下载到homer文件夹里
$ perl ./homer/configureHomer.pl -install hg38

HOMER支持的基因组有：
Human (hg18, hg19, hg38), Mouse (mm8, mm9, mm10), Rat (rn4, rn5, rn6), Frog (xenTro2, xenTro3), Zebrafish (danRer7), Drosophila (dm3), C elegans (ce6, ce10), S. cerevisiae (sacCer2, sacCer3), pombe (ASM294v1), Arabidopsis (tair10), Rice (msu6), 你也可以使用自己定义的基因组fasta格式和GTF格式的文件。

下载基因组后，homer文件夹里会多出来一个子文件夹，叫“data”，里面有：

$ ll
total 24
drwx------ 2 fangy04 fangy04 8192 Dec 10 14:08 GO
drwx------ 2 fangy04 fangy04 8192 Dec 10 14:08 accession
drwx------ 3 fangy04 fangy04 8192 Dec 10 14:09 genomes

其中，genomes文件夹里会有一个hg38的子文件夹，是我刚才下载的基因组文件，在hg38文件夹里有这些文件：

$ ll
total 4539344
drwxr-xr-x 5 fangy04 fangy04       8192 Oct 19  2019 annotations
-rw-r--r-- 1 fangy04 fangy04      11672 Oct 19  2019 chrom.sizes
-rw-r--r-- 1 fangy04 fangy04 3273481150 Jan 16  2014 genome.fa 
-rw-r--r-- 1 fangy04 fangy04    3270764 Oct 19  2019 hg38.aug
-rw-r--r-- 1 fangy04 fangy04   43703292 Oct 19  2019 hg38.basic.annotation
-rw-r--r-- 1 fangy04 fangy04  705525298 Oct 19  2019 hg38.full.annotation
-rw-r--r-- 1 fangy04 fangy04     164707 Oct 19  2019 hg38.miRNA
-rw-r--r-- 1 fangy04 fangy04  528510489 Oct 19  2019 hg38.repeats
-rw-r--r-- 1 fangy04 fangy04   24471159 Oct 19  2019 hg38.rna
-rw-r--r-- 1 fangy04 fangy04   29505106 Oct 19  2019 hg38.splice3p
-rw-r--r-- 1 fangy04 fangy04   29473833 Oct 19  2019 hg38.splice5p
-rw-r--r-- 1 fangy04 fangy04    3270627 Oct 19  2019 hg38.stop
-rw-r--r-- 1 fangy04 fangy04    3270740 Oct 19  2019 hg38.tss
-rw-r--r-- 1 fangy04 fangy04    3270352 Oct 19  2019 hg38.tts
drwxrwxr-x 2 fangy04 fangy04       8192 Dec 10 14:53 preparsed

*.fa or *.fa.masked是unmasked/masked基因组文件
genome.tss (参考基因转录起始位点的位置)
genome.tts (参考基因转录终止位点的位置)
genome.splice3p (参考基因3'剪切位点的位置)
genome.splice5p (参考基因5'剪切位点的位置)
genome.aug (参考基因翻译起始密码子的位置)
genome.stop (参考基因翻译终止密码子的位置)
genome.rna (refseq RNA definition file)
genome.repeats.rna (repeat RNA definition file)
genome.basic.annotation (exon/intron/TSS/TTS/intergenic区域注释)
genome.full.annotation (CpG岛和重复区域注释)
annotation/ subdirectory (包含GO的注释文件)

（四）准备HOMER输入文件

可以使用HOMER peak files文件或者任何BED文件。
HOMER的peak文件应该是最少5列的（用TAB分隔）：
第一列: Unique Peak ID（peak的ID）
第二列: chromosome（染色体）
第三列: starting position（起始位置）
第四列: ending position（结束位置）
第五列: Strand (+/- or 0/1, where 0="+", 1="-")（链）

那么问题来了：How many sequences can HOMER handle?
从理论上讲，可以很多(即数百万)。该算法被设计用于10k条目标序列和50k背景序列。如果你使用大量的序列来运行findMotifs.pl时，你可能希望使用-b参数，该参数切换到cumulative binomial distribution进行motif评分，这可以使它更快地计算，并在使用大量序列时给出基本相同的结果。在findMotifsGenome.pl中默认使用二项式。

（五）Finding motifs

在某些情况下，用户可能要考虑什么情况下用哪一个工具，所以HOMER提供了以下选项:

启动子：使用findMotifs.pl分析启动子（这个功能不是本文探讨的主要内容，请参考官网：这里），如果你希望使用一组特定的启动子作为背景，请将它们放在一个文本文件中(第一列是ID)，并使用-bg <background IDs file>参数。在目标和背景中发现的基因将被从背景中移除，这样它们就不会相互抵消。

基因组区：比如说enhancers，或其他调控元件。当分析峰/区域，使用findMotifsGenome.pl，你可以指定适当的背景区域的基因组区域，通过将它们放置在它们自己的峰文件中，并使用-bg <background peak file>。例子：当尝试寻找特定于某一细胞类型特定峰的motif时，指定两种细胞类型的共同峰作为背景——这将有助于消除主要motif，并揭示共富集motif。如果峰在外显子附近，指定外显子上的区域作为背景，以消除triplet偏倚。

这里我的BED峰文件是来自ATAC-seq的峰，且不在启动子区域，所以应该使用findMotifsGenome.pl。代码非常简单，就一行，但是建议在运行前先看一下下面的参数介绍部分，因为根据你的实验分析需要，可能要更改参数设置，这里我只用了默认值，只作为练习：

# finding motif
# 基本格式
# findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size # [options]
$ findMotifsGenome.pl /gpfs/home/fangy04/peaks.bed /gpfs/home/fangy04/homer/data/genomes/hg38 /gpfs/home/fangy04/pks_MotifOutput_2/ -size 200

（六）findMotifsGenome.pl的一些重要参数

（1）Masked vs. Unmasked Genome (-mask or hg18 vs. hg18r)
实际上，这通常并不重要。由于HOMER是一种差异motif发现算法，重复序列通常在目标序列和背景序列很常见。然而，一个转录因子与某一类重复序列结合的情况并不少见，这可能导致处理几个大范围的相似序列，导致结果偏倚的原因。通常使用masked会更安全。

（2）Region Size (-size <#>, -size <#>,<#>, -size given, 默认: 200，这个参数是必须输入的)
用于寻找motif的区域的大小是重要的。如果从一个转录因子中分析ChIP-Seq峰，推荐用50 bp建立一个给定转录因子结合的主motif，200 bp用于寻找一个转录因子的主motif和“共富集”motif。对于组蛋白标记区域，500-1000 bp可能更合适(即H3K4me或H3/H4乙酰化区域)。理论上，HOMER可以处理非常大的区域(即10kb)，但是区域越大，序列就越多，执行时间也就越长。这些区域将以峰的中心为基础。如果你更愿意用偏移量，可以指定-size -300,100来搜索大小为400的区域，该区域位于峰中心上游100 bp(如果在假定的TSS区域进行motif查找，这是很有用的)。如果你的区域长度不一定，使用-size given选项HOMER将使用准确的区域作为input。

（3）Motif length (-len <#>或-len <#>,<#>,..., default： 8,10,12)
指定要找到的motif长度。HOMER会分别找到每个motif大小的结果，然后在最后把结果结合起来。随着长度的增加，寻找motif所需的时间也大大增加。一般来说，在尝试更长的长度之前，最好先尝试短的motif(即小于15)。HOMER也可以找到更长的motif，但最好使用较小的数据集来找长motif，否则可能会花费太长时间(或占用太多内存)。减少总资源消耗的另一个技巧是减少背景序列的数量(-N <#>)。

（4）Mismatches allowed in global optimization phase (-mis <#>, default: 2)
HOMER寻找有可能的候选motif时，首先检查原始oligos富集，允许错配。你允许的错配越多，算法就越敏感，特别是对于较长的motif。然而，这也会降低算法的速度。如果搜索长度超过12-15 bp的motif，最好将这个值增加到至少3甚至4。

（5）Number of CPUs to use (-p <#>, default: 1)
HOMER现在是多核兼容的。它不是完全并行的。一般来说，motif的长度越长，加速的效果越好。

（6）Number of motifs to find (-S <#>, default: 25)
指定要查找的每个长度的motifs的数量。25已经足够多了。建议减少这个数量，特别是对于长motif，以减少总执行时间。

（7）Normalize CpG% content instead of GC% content ("-cpg")
如果HOMER在寻找类似“CGCGCGCG”的motif时卡住了，你可以尝试禁用GC/CpG 标准化(-noweight)。

（8）Region level autonormalization (-nlen <#>, default 3, -nlen 0 是禁用)
（9）Motif level autonormalization (-olen <#>, default 0 ,你还可以disabled)
自动标准化试图消除低排序核苷酸的序列偏倚。(1-mers, 2-mers…< # >)。区域水平自动标准化(默认为1/2/3 mers)尝试通过调整背景区域的权重对其进行标准化。如果这不能完成任务(自动标准化不能保证消除所有序列偏差)，你可以尝试motif水平的自动标准化(-olen <#>)。在发现新motif时，在oligo表上执行自动标准化程序。

（10）User defined background regions (-bg <peak file of background regions>)
你可以自己选择背景区域!!这非常有用，因为荷马是一个差异motif发现算法。例如，你可以给HOMER一组peaks(由另一个因子共同结合的峰），并将它们与其他峰进行比较（比较特定一种细胞类型的ChIP-Seq峰与特定于另一种细胞类型的峰）。为了实现这一功能，你要准备第二个peak/BED文件，并使用参数-bg。HOMER会使用mergePeaks自动检查背景峰是否与目标峰重叠，并放弃重叠区域。

（11）Hypergeometric enrichment scoring (-h)
默认情况下，findMotifsGenome.pl使用二项分布来给motifs评分。当背景序列的数量大大超过目标序列时，这种方法很好；但是，如果你使用上面的-bg选项，并且背景序列的数量小于目标序列，那么使用超几何分布(-h)比较好。

（12）Find enrichment of individual oligos (-oligo)
这将在名为oligo.length.txt的输出目录中创建输出文件。

（13）Search for RNA motifs (-rna)
如果是RNA数据(例如Clip-Seq)，该选项只会限制HOMER搜索+链(相对于峰),并将输出RNA motif logo(例如U代替T)。它还将试图把发现的motif与RNA motif数据库比较，但是RNA motif数据库现在只包含microRNA。

（14）Optimize motifs (-opt <motif file>)
与寻找新的motif不同，HOMER将试图优化所提供的motif。当你试图改变一个motif的长度，或找到一个非常长的给定motif时，这个参数就非常有用。例如，如果你指定-opt <file>和-len 50，它将尝试将motif扩展到50bp并对其进行优化。

（七）findMotifsGenome.pl是如何进行分析的？

findMotifsGenome.pl是通过一系列步骤来寻找高质量的motifs的:

1. 验证峰/BED文件
HOMER确保你有有效的峰，并检查确保你有唯一的峰标识符（就是峰的ID，也就是BED文件里的第一列）。如果有重复，它将通知你，并将为峰名添加数字，以确保它们对于下游分析是唯一的。

2. 从基因组中提取与输入文件区域相对应的序列，过滤包含>70%“N”的序列
这一步很容易解释。HOMER也会丢弃以“N”为主的序列。如果你觉得你扔掉了太多的序列，试着在一个unmasked的基因组上运行findMotifsGenome.pl。

3.计算峰序列的GC/CpG含量
CpG岛是哺乳动物基因组中序列含量偏差最大的来源，不幸的是，它在转录起始位点附近，所有的活动都在那里。。。默认情况下，HOMER跟踪GC%(使用-cpg来使用CpG%)。

4. 准备基因组序列作为背景序列
这个步骤只有在你第一次从给定大小的区域(-size <#>)去找motif的时候才会完成。HOMER取基因TSS附近的区域(+/- 50kb)，并将其分割成指定大小的区域。然后计算它们的GC/CpG%，并将其存储起来，以便下次在类似大小区域中搜索motifs时加速执行。

5. 随机选择背景区域进行motif查找
由于HOMER是差异motif发现算法，它必须使用背景序列区域作为对照。默认情况下，HOMER选择了足够多的随机背景区域，这样区域总数为50000，或者是峰总数的2倍，甚至更多的数量(-N <#>)。使用的总序列越多，程序运行的速度就越慢，但你要确保有足够的背景区域来可靠地估计motif频率。HOMER试图选择与输入序列的GC含量分布相匹配的背景区域(以5%的增量)。例如，如果你的输入区域富含GC, HOMER会从基因组中GC丰富的区域中随机选择区域作为对照。如果提供了自定义背景区域(-bg <peak/BED file>)，HOMER会自动确保这些区域不与目标区域重叠(使用mergePeaks)。自定义区域仍将按照GC含量进行标准化。

6. 序列偏差的自动标准化（Autonormalization）
Autonormalization是HOMER提供的一个独特的程序，试图消除由lower-order oligo（这个不会翻译）序列引入的偏差。它是通过假设你的目标区域和背景区域在1-mers, 2-mers, 3-mers中是不平衡的来工作的。Autonormalization的oligo的最大长度由-nlen <#>指定(默认为3，禁用-nlen 0)。例如，目标序列中不应该比背景中有更多的A。在计算出每个oligo的不平衡后，对每个背景序列的权重进行少量的调整，以帮助对不平衡进行标准化。如果目标序列具有丰富的A，那么含有较多A的背景序列会被赋予较高的权重，而含有较少A的背景序列会被赋予较低的权重。权重只增加了少量，并且这个过程在爬山算法优化（hill climbing optimization）中重复了很多次。这一程序有助于消除与某些基因组区域相关的序列偏倚，或可能由偏倚(如偏倚测序)引入的的实验结果偏倚。

7. 检查已知motifs的富集程度
HOMER根据目标和背景序列筛选可靠的motif库进行富集，返回p值小于0.05的富集motifs。先进行已知的motif富集，因为它通常比较快，可以更快地看到你的目标区域富集了什么。已知motif的富集将被输出一个名为“knownResults.html”的文件。

8. de novo寻找新的motif
默认情况下，HOMER会搜索len 8、10和12 bp的motifs(使用-len <#，#，#>，之间没有空格，即-len 6,10,15,20)。从de novo motif查找的输出将显示在“homerResults.html”文件中。

（八）findMotifsGenome.pl的输出文件

上面运行后，大概运行了1个小时，生成了以下文件：

homerMotifs.motifs8/10/12：这些是de novo 查找motif的输出文件，由motif长度分隔，并代表算法的单独运行。
homerMotifs.all.motifs：简单地由所有homerMotifs.motifs组成的连接文件。
motifFindingParameters.txt：用于执行findMotifsGenome.pl的命令，包含你使用的参数
knownResults.txt：包含关于已知motifs丰富的统计信息的文本文件(在EXCEL中打开)。
seq.autonorm.tsv：用于lower-order oligo标准化的autonormalization统计。
homerResults.html：de novo查找motif的格式化输出。

homerResults子文件夹：包含homerResults.html网页文件，包括motif<#>.motif文件，用于寻找每个motif的具体实例。
knownResults.html：查找已知motifs的格式化输出。
knownResults子文件夹：包含knownResults.html网页文件，包括known<#>.motif文件。

（九）理解输出的Motif文件

HOMER的真正输出是“*.motif”文件，其中包含必要的信息，以识别motif的未来实例。它们在findMotifs.pl和findMotifsGenome.pl的输出目录中。一个典型的motif文件看起来是这样的（下面这个是我运行后生成的一个motif文件）:

>NNATGASTCATH   Fra1(bZIP)/BT549-Fra1-ChIP-Seq(GSE46166)/Homer  7.011739    -10.259778  0   T:4.0(57.14%),B:3353.3(3.23%),P:1e-4
0.227   0.228   0.279   0.266
0.285   0.199   0.287   0.229
0.461   0.212   0.256   0.070
0.017   0.005   0.016   0.962
0.015   0.009   0.745   0.231
0.949   0.010   0.001   0.040
0.097   0.363   0.411   0.130
0.028   0.013   0.026   0.933
0.189   0.789   0.009   0.013
0.963   0.012   0.001   0.024
0.074   0.253   0.198   0.475
0.244   0.347   0.172   0.237

第一行以一个“>”开头，后面跟着各种信息，其他行是每个位置的各个核苷酸具体概率(A/C/G/T)。标题行实际上是用制表符分隔的，并包含以下信息:

1. “>”+序列(实际上不用于任何东西，可以是空的)示例:>NNATGASTCATH
2. motif名称(如果几个motif在同一个文件中，应该是唯一的)。例子:Fra1(bZIP)/BT549-Fra1-ChIP-Seq(GSE46166)/Homer
3. log odds检测阈值，用于确定结合的vs未结合位点。示例:7.011739
4. 富集的log-p值
5. 0(表示向后兼容性的占位符，在旧版本中用于描述“gapped” motif，结果证明它并不是很有用)
6. 发生的信息，用逗号分隔，例如:T:4.0(57.14%),B:3353.3(3.23%),P:1e-4
   T:带motif的靶标序列数，占靶序列总数的%
   B:带motif的背景序列数目，占总背景的%
   P:最终富集P值
7. 用逗号分隔的Motif统计信息（这一个信息是与链有关的，因为我的BED文件里链的哪一列没有指定是+/-，我输入的是.，所以这里没有第7个信息），官网的例子：Tpos:100.7,Tstd:32.6,Bpos:100.1,Bstd:64.6,StrandBias:0.0,Multiplicity:1.13
   Tpos: motif在目标序列中的平均位置(0 =序列开始)
   Tstd:目标序列中位置的标准偏差
   Bpos:motif在背景序列中的平均位置(0 =序列开始)
   Bstd:背景序列中位置的标准差
   StrandBias:在+链出现与在-链出现的log比例。
   Multiplicity:在具有一个或多个结合位点的序列中，每个序列平均出现的次数。

（十）如何评价你找到的这些motif质量？

这一部分的图我就用官网上的了，因为比较有代表性，你可以用下面的方法来评价你自己找到的motif的好坏与否，和可信度。

WARNING:因为这部分是最难理解的，所以官网在这里又说一遍。HOMER会在motif结果的旁边打印出最佳的motif猜测，但在你确定这个motif在你的序列里富集之前，强烈建议你看看比对的情况!!下面是一个可能发生的例子:

在这种情况下，HOMER已经确定YY1是“最好的猜测”匹配这个新motif。点击“More Information”看看有什么:

在这个例子中，你可以看到，motif比对在已知的YY1 motif的边缘，而不是YY1 motif 的核心(CAAGATGGC)。这并不意味着YY1 motif在你的数据中没有得到富集，但除非有其他motif结果显示YY1 motif在其他地方得到富集，否则YY1 motif不太可能在你的数据集中得到富集（有点绕，多读几遍就懂了）。

和往常一样，记住HOMER是一个de novo的查找motif的工具!! 即使HOMER会猜出最好的匹配，如果它是一个新的motif，无论如何你都不想相信那个匹配。因此，你可以看到比对的重要性，并得到证据来支持或反对这个匹配。

在很多情况下，HOMER会找到p值很低的motif，但这些motif可能看起来“不太靠谱”。质量差的motif可以大致分为以下几类:

（1）Low Complexity Motifs（低复杂度的motif）:
（在v3版本里比较少见）这些类型的motif倾向于在每个位置显示相同的1、2、3或4个核苷酸的集合。例如:

当目标序列和背景序列之间存在系统性偏差时，这些motif通常会出现。通常它们的GC含量非常高，在这种情况下，你可以尝试在motif finding命令中添加-gc，以总GC含量而不是CpG含量进行标准化。

其他时候，例如，在一些生物体中将启动子序列与随机基因组背景序列进行比较时，会显示出对嘌呤或嘧啶的偏好。HOMER非常敏感，所以如果序列的组成有偏差，HOMER很可能会发现。新版本中的Autonormalization解决了这个问题。

（2）Simple Repeat Motifs（简单重复的motifs）:
(在v3.0+版本里很少出现) 有时motifs还会显示出重复的序列：

通常像这样的motif将伴随着几个其他的motif，看起来高度相似。除非有充分的理由相信这些可能是真实的，否则最好假设背景可能有问题。如果你的目标序列在外显子和其他类型的序列上高度富集，就会出现这种情况，如果-gc不起作用，你可能必须仔细考虑你要分析的序列类型并试图匹配它们。(例如启动子vs.启动子，外显子vs.外显子等)。你也可以尝试使用-olen <#>在寡核苷酸水平上自动标准化序列偏差。

（3）Small Quantity Motifs / Repeats（小数量的motif）:
这个解释起来有点难。这些看起来像真正的motif，但在目标序列里只有一个非常低的百分比。例如，一个oligo或一个repeat，在一些目标序列中出现，看起来是一个重要motif。从统计学上讲，它们是富集的，但可能不是真的。当从一个小的调控基因列表中寻找启动子的motif时，这是最大的问题。原则上，motif在不到5%的靶序列中存在，这可能是个问题。

（4）Leftover Junk（剩下的“垃圾”？我并不想这么翻译。。。）
这些motif出现在你的结果列表比较靠后的位置。如果一个元件在你的序列中高度富集，HOMER会找到它，掩盖它，然后继续寻找motif。在这种情况下，许多HOMER发现的其他motif将被高度富集的motif抵消。例如(另一个pu1例子):

排名很高的motif:

排名靠后的motif：

这并不一定是阴性的结果，但是应该把它们放在特定的环境中。这通常发生在ChIP-Seq数据集里，并与大量的结合位点紧密结合。这些“其他”motif很可能也能结合PU.1，只是可能代表低亲和力结合位点，但不建议在这种情况下关注它们。一个更安全的方法是重复motif寻找步骤，但是用不包含高富集motif的那些序列。或通过添加-mask <motif file>到motif finding 命令再尝试一下。