一、 基本介绍
HOMER（Hypergeometric Optimization of Motif EnRichment），用于Motif的发现与查找。homer软件找motif整合了两个方法，包括依赖于数据库的查询和de novo的推断。HOMER通过比较目标和背景序列，使用ZOOPS scoring (zero or one occurrence per sequence)和超几何检验进行富集分析。HOMER也可以对peaks进行注释。
二、 使用方法
(1) 下载数据
刚安装的Homer实际没有包含参考序列或注释数据 ，但是可以使用 configureHomer.pl下载数据。homer包分为4种：
- SOFTWARE：homer工具包，包含一些常用数据。
- ORGANISMS：物种特异性的数据，包含Gene accession、gene description、GO analysis信息。大多数数据来自于NCBI Gene database。下载promoter 或 genome 数据时，会自动下载对应Organism 包。
- PROMOTERS：Promoter 序列信息和Promoter 富集分析的文件。大多数数据来自RefSeq transcript。
- GENOMES：基因组序列及其注释信息。
#使用homer附带的配置脚本来下载数据库。
perl /home/analysis/miniconda3/envs/bioinf3/share/homer/configureHomer.pl -list #查看Homer数据包
ORGANISMS：fly-o v6.3 Drosophila melanogaster (fly) accession and ontology information
PROMOTERS：fly-p v5.5 fly promoters (fly)
GENOMES：dm6 v6.4 fly genome and annotation for UCSC dm6
GENOMES：dm3 v6.0 fly genome and annotation for UCSC dm3
#安装homer数据包
perl /home/analysis/miniconda3/envs/bioinf3/share/homer/configureHomer.pl -install dm6
perl /home/analysis/miniconda3/envs/bioinf3/share/homer/configureHomer.pl -install fly-o
perl /home/analysis/miniconda3/envs/bioinf3/share/homer/configureHomer.pl -install fly-p
#安装位置
ls /home/analysis/miniconda3/envs/bioinf3/share/homer/data/
configureHomer.pl的选项：
-list：列出可用的包。
-install：安装homer需要用到的数据包，安装在homer文件目录的/data/genomes/里。
-version：安装homer包时，可以指定包版本。
-remove：移除包。
-update：更新所有包到最新版本。
-reinstall：强制重新安装所有已经安装过的包。
-all：安装所有包。
-url：安装时，使用的资源地址，默认http://homer.ucsd.edu/homer/。
(2) 寻找motifs
findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size # [options]
findMotifsGenome.pl TF_homer.bed dm6 homer/TF_motif -p 4 -len 8,10,12
findMotifsGenome.pl TF_homer.bed dm6 homer/m6A_motif -p 4 -len 6,8,10 -rna
peak/BED file：bed格式输入文件。HOMER peak文件应该至少有5列，用制表符分隔，其他列将被忽略。
Column1: Unique Peak ID
Column2: chromosome
Column3: starting position
Column4: ending position
Column5: Strand (+/- or 0/1, where 0="+", 1="-")
genome：参考基因组。
output directory：结果文件的路径。
-p：线程数。
-len ("-len <#>" or "-len <#>,<#>,...", default 8,10,12)：motif长度设置，默认8,10,12。越大需要的计算资源越多。
-S：结果输出多少motifs，默认25。
-size ("-size <#>", "-size <#>,<#>", "-size given", default: 200)：用来进行motif发现的片段长度，默认为200。例如：-size -300,100表示peak上游100bp，下游300bp区域。如果想在提供的peak中寻找motifs，使用参数-size given。然而，对于转录因子peaks，大多数motifs 被发现位于peak 中心 +/- 50-75 bp的范围内，所以最好根据peak 的大小将寻找motif的区域设为固定值。如果从一个转录因子中分析ChIP-seq峰，推荐用50bp建立一个给定转录因子结合的主motif；200bp用于寻找一个转录因子的主motif和共同富集motif；对于组蛋白标记区域，500-1000bp更合适。
-rna: 输出RNA motif，使用RNA motif数据库。
-norevopp：只在peak所在的链搜索motif，不进行反义链搜索。HOMER默认会在两条链上都进行搜索。
-nomotif：don't search for de novo motif enrichment.
-mset <vertebrates|insects|worms|plants|yeast|all>：check against motif collects, default: auto.
-basic：just visualize de novo motifs, don't check similarity with known motifs.
-noknown：don't search for known motif enrichment, default: -known.
-mknown <motif file>：known motifs to check for enrichment.
-bg：自定义背景序列。如果将某种实验的peak与另一种实验peak相比较，可以再创建一个peak/BED文件（参数："-bg <peak/BED file>"），将会对背景进行移除GC-bias操作和自动标准化。例如，当尝试寻找特定于某一细胞类型特定峰的motif时，指定两种细胞类型的共同峰做背景。
-N：用于motif寻找的背景序列数目，HOMER 使用的motif 寻找算法需要使用背景序列区域作为对照，可能选择50000 或 peaks总数两倍的随机背景序列。default=max(50k, 2x input)；耗内存参数。HOMER 会选择和目标数据一致GC 含量分布的序列作为背景序列。例如，目标序列是GC高含量的，那么背景序列也会如此。
-cpg：Normalize CpG% content instead of GC% content. 考虑到HOMER 可能卡在CGCGCGCG这样的motifs。
-gc：使运算按总GC含量而不是CpG岛含量进行标准化，是默认的。
-mis：motif错配碱基数，默认2bp。如果寻找12-15 bp Motif，可以设置3-4bp的错配。
-h：使用超几何检验代替二项式分布。Motif 富集分析使用超几何分布和二项式分布。一般情况下，序列较多或者背景序列远远多于目标序列，二项式分布计算比较快，因此findMotifsGenome.pl默认使用二项式分布。当自定义背景序列时，这时序列较少，使用超几何检验比较好("-h")。findMotifs.pl用于启动子分析，并且默认使用超几何检验。
#预处理，将MACS2得到的peaks.narrowPeak或summits.bed文件转换成homer需要的bed文件。
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' cbx7_peaks.narrowPeak > sample_homer.bed
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' TF_summits.bed > TF.bed
#由于dm6和ensembl的基因组格式不太相同，所以需要进行转换。
cat TF.bed | cut -f 2 | sort | uniq -c
cat TF.bed | awk '{if($2=="2L"||$2=="2R"||$2=="3L"||$2=="3R"||$2=="4"||$2=="X"||$2=="Y") {print $1"\tchr"$2"\t"$3"\t"$4"\t"$5}}' > TF_homer.bed
cat TF_homer.bed | cut -f 2 | sort | uniq -c
#寻找富集的motifs
findMotifsGenome.pl TF_homer.bed dm6 homer/TF_motif -p 4 -len 8,10,12 2>>log/homer_findMotifsGenome_log.txt &
findMotifsGenome.pl TF_homer.bed dm6 homer/m6A_motif -p 4 -len 5 -rna 2>>log/homer_findMotifsGenome_log.txt &
(3) motifs结果
① homerMotifs文件
分为motifs8、motifs10、motifs12和all.motifs。每个motif信息以 > 起始，首行是motif各种统计信息，其他行对应motif各个位点A/C/G/T的频率。
>TCTCTCTC 1-TCTCTCTC 7.868595 -30.087445 0 T:94.0(13.09%),B:2725.2(5.68%),P:1e-13 Tpos:96.5,Tstd:56.2,Bpos:97.7,Bstd:70.6,StrandBias:-0.2,Multiplicity:1.53
0.001 0.001 0.001 0.997
0.012 0.986 0.001 0.001
0.001 0.001 0.001 0.997
0.032 0.888 0.061 0.019
0.044 0.087 0.188 0.681
0.016 0.561 0.130 0.293
0.001 0.001 0.031 0.967
0.014 0.971 0.014 0.001
首行内容：
序列 —— TCTCTCTC
motif名称 —— 1-TCTCTCTC
log odds检测阈值 —— 7.868595：用于确定结合的vs未结合位点。
富集P-value对数值 —— -30.087445
用于老版本格式的占位符 —— 0
T:94.0(13.09%),B:2725.2(5.68%),P:1e-13 ——
- T:#(%) - 包含motif的目标序列数 / 目标序列总数
- B:#(%) - 包含motif的背景序列数 / 背景序列总数
- P:# - 富集显著性 p-value
用逗号分隔Motif统计量 ——
Tpos:96.5,Tstd:56.2,Bpos:97.7,Bstd:70.6,StrandBias:-0.2,Multiplicity:1.53：
- Tpos: motif在目标序列中的平均位置 (0 = start of sequences)
- Tstd: motif在目标序列中位置的标准偏差
- Bpos: motif在背景序列中的平均位置 (0 = start of sequences)
- Bstd: motif在背景序列中位置的标准偏差
- StrandBias: log (+ strand occurrences / - strand occurrences)
- Multiplicity: 在一个peak中motif出现的平均次数
② html文件
低复杂度的motif序列的核苷酸倾向于都是同一种核苷酸，从而导致GC含量异常。当目标序列和背景库中序列之间存在系统性偏差时会导致这样的结果。通它们的GC含量非常高。在这种情况下，可以在motif分析命令中添加参数“-gc”，从而使运算按总GC含量而不是CpG岛含量进行标准化。HOMER非常敏感，所以如果序列的组成有偏差，HOMER 很可能会发现。新版本中的Autonormalization可以尽量减小这个问题的发生。
motifs 有时候会出现一些序列模式的重复。这种motifs一般会有数个差不多序列的motifs。除非有充分的理由相信这些可能是真实的，否则背景序列可能有问题。如果你的目标序列在外显子和其他类型的序列上高度富集，就会出现这种情况；并且如果"-gc"参数也不能改善结果，用户就需要考虑自己正在分析序列的类型以及怎么去匹配他们。
低质量和低重复的Motifs。这种发生在motif看起来很靠谱，但是在序列中出现的百分比缺失很低的。例如，寡核苷酸和重复序列在用户序列中出现从而导致极高的显著性。统计上使显著的，但是事实上却并不是。一些调节基因的启动子序列会发生这样的事情。原则上，motif 在不到5%的靶序列中存在的话，这个motif不太可信。
一些高质量的motifs 可能会出现在结果的后面。如果一个motif在序列中高度富集， HOMER 会发现他们，然后继续寻找新的motif。后续的motifs 可能会掩盖先前找的motif。常常发生于ChIP-Seq数据中，免疫沉淀的蛋白高表达以及与大量的结合位点紧密结合。这些motifs 可能结合第一位的motif，但是亲和性不高。处理这种情况的方法是重复motif 分析的过程，但是丢掉 the top motif（排名靠前的），添加参数-mask <motif file>就可以在motif分析过程中忽略这些top motif 。
(4) peaks注释
annotatePeaks.pl <peak file | tss> <genome version> [additional options...]
annotatePeaks.pl TF_homer.bed dm6 1>homer/TF_peakAnn.xls 2>>log/homer_annotatePeaks_log.txt &
工作原理：注释peaks/regions的过程分为两个主要部分。第一种是确定与最近的TSS的距离，并将峰分配给该基因。第二种是确定峰/区域中心所占据区域的基因组注释。
The process of annotating peaks/regions is divided into two primary parts. The first determines the distance to the nearest TSS and assigns the peak to that gene. The second determines the genomic annotation of the region occupied by the center of the peak/region.
※ 到最近的 TSS 的距离：
annotatePeaks.pl 在"/homer/data/genomes/<genome>/<genome>.tss"中加载一个文件，该文件包含 RefSeq 转录起始位点的位置。使用这些位置来确定最近的TSS，报告距离Distance to TSS（负值表示TSS的上游，正值表示下游），以及与位点相关的各种注释信息。
※ 基因组注释：
为了根据重要的基因组特征注释给定峰的位置，annotatePeaks.pl 调用一个单独的程序（assignGenomeAnnotation）来有效地将峰分配给基因组范围内数百万个可能的注释之一。提供两种类型的输出。第一个是Basic Annotation，包括峰是否在TSS (transcription start site), TTS (transcription termination site), Exon (Coding), 5' UTR Exon, 3' UTR Exon, Intronic, or Intergenic，这些都是许多研究人员感兴趣的常见注释。第二个是Detailed Annotation还包括更详细的注释，也考虑了重复元件和CpG岛。由于某些注释重叠，因此优先级基于以下内容进行分配（在并列的情况下，它是随机的[即如果有两个重叠的重复元素注释]）：
1. TSS (by default defined from -1kb to +100bp)
2. TTS (by default defined from -100 bp to +1kb)
3. CDS Exons
4. 5' UTR Exons
5. 3' UTR Exons
6. **CpG Islands
7. **Repeats
8. Introns
9. Intergenic
** Only applicable for the "Detailed Annotation".
尽管 HOMER 不允许显式更改 TSS 区域的定义（-1kb 到 +100bp），但可以通过在 EXCEL 中按"到最近的 TSS 的距离"列对注释输出进行排序，然后选择感兴趣的区域内的注释输出。
三、 bedtools
(1) bedtools介绍
bedtools可以处理基因组广泛使用的BAM、BED、GTF、VCF等格式的文件，进行取交集、并集、补集、计数、格式转变等操作。bedtools总共有二三十个工具/命令来处理基因组数据。比较典型而且常用的功能举例如下：
- 格式转换：bam转bed（bamToBed），bed转其他格式（bedToBam，bedToIgv）
- 对基因组坐标的逻辑运算：交集（intersectBed，windowBed），邻集（closestBed），补集（complementBed），并集（mergeBed），差集（subtractBed）
- 计算覆盖度（coverage）：（coverageBed，genomeCoverageBed）
注意事项：
- bedtools默认输入文件的分隔符为TAB，除了bam格式的文件。
- 如果未使用-sorted参数，则bedtools默认不支持大于512M的染色体。
- -sorted参数和-g参数必须存在一个。
- 当进行多个文件比较时，染色体的命名方式必须统一，chrX和X不可以同时存在。
- bed是0-based，gff、vcf是1-based，bedtools自身可以识别gff、vcf文件的正确位置，所以不需要再对gff文件初始位置减去1转变成bed格式。
- 当比较两个文件时，fileB一般会被加载到内存中，fileA会被遍历各行与内存中的fileB进行比较。所以一般将比较小的文件作为fileB。
(2) bedtools window
#窗口，例如寻找lncRNA上下游10kb以内的蛋白编码基因
bedtools window -a lnc.bed -b pr.bed -l 10000 -r 10000 > lnc_pr_window.txt
#把lnc.bed比对到pr.bed寻找overlap，寻找范围，上游10000bp，下游10000bp
(3) bedtools closest
#最近，例如寻找距离lncRNA最近的蛋白编码基因
sort -k1,1 -k2,2n lnc.bed > lnc.sorted.bed
sort -k1,1 -k2,2n pr.bed > pr.sorted.bed
bedtools closest -a lnc.sorted.bed -b pr.sorted.bed > lnc_pr_closest.txt
(4) bedtools intersect
#取交集
bedtools intersect -a cpg.bed -b exons.bed > result.bed
#从cpg.bed中取出与exons.bed不重叠的区域
bedtools intersect -a cpg.bed -b exons.bed -v > result.bed
#取至少有50%区域重叠的区域
bedtools intersect -a WT_rep1_peaks.narrowPeak -b WT_rep2_peaks.narrowPeak -f 0.5 -F 0.5 -e > WT_macs2_intersect.bed
-a：查询文件，可以是BAM/BED/GFF/VCF格式，A中的每一行都会与B进行比较以找出重叠的部分。
-b：参考文件，依照-b输入的文件在A里面找重叠的部分，可以是BAM/BED/GFF/VCF格式。-b可以输入多个文件，用空格分开即可。现在可以用通配符（*）识别多个输入文件。
-filenames：标明重叠区域是源于和哪个B文件发生重叠。
-f：指定两个区域最小重叠多少bp才算与A发生重叠，默认1bp，也可以是一个比例。
-F：指定两个区域最小重叠多少bp才算与B发生重叠，默认1bp，也可以是一个比例。
-e：要求A或B满足最小分数。换句话说，如果-e与-f 0.90和-F 0.10一起使用，则要求覆盖90%的A或10%的B。如果没有-e，两个分数都必须满足。
-wa：在交集中输出A的完整区间，而不仅是重叠的部分。
-wb：在交集中输出B的完整区间，而不仅是重叠的部分。
-v：输出A unique的区域。
-wo：输出A、B重叠区域重叠的bp。
-wao：输出A、B重叠区域重叠的bp，假如B中不存在A的这个区域则输出0。
-s：在找交集的过程中需要考虑链特异性。
-S：强制不考虑正负链，而只考虑基因组区域是否重叠。
-split：Treat "split" BAM or BED12 entries as distinct BED intervals.
网友评论