【生物信息笔记】HOMER 找 DNA motif

What is HOMER？

HOMER is a software for motif discovery and ChIP-Seq analysis

HOMER软件是Linux command line based，常用来查找DNA motif ，偶尔以及一些ChIP-seq的分析（如，peak calling）。

• 其他的DNA motif 查找软件如非常有名的在线tool: MEME。
• 其他的peak calling tool：Macs2 （更常用）

---

感兴趣HOMER其它功能可以到它主页去查找，下载与安装的方法也可以在主页里找到。

安装使用如下：

## Download and install homer (Hypergeometric Optimization of Motif EnRichment)
## // http://homer.salk.edu/homer/
## // http://blog.qiubio.com:8080/archives/3024
## pre-install: Ghostscript，seqlogo,blat
cd ~/biosoft
mkdir homer && cd homer
wget http://homer.salk.edu/homer/configureHomer.pl
perl configureHomer.pl -install
perl configureHomer.pl -install hg19
perl configureHomer.pl -install hg38

如果是对MACS找到的peaks记录文件，还需提取对应的列给HOMER作为输入文件：
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' sample_peaks.bed >sample_homer.bed 如果不熟悉 awk就只好手动改。
findMotifsGenome.pl sample_homer.bed hg19 motifDir -len 8,10,12
最后得到的文件夹里面有一个详细的网页版报告，所以很多人都喜欢用这个软件，而且HOMER 这个软件是一个大杂烩，能解决几乎所有的高通量测序数据的分析。

这里记下的只是DNA motif的查找使用方法：

1. Gene/Promoter-based Analysis:
   findMotifs.pl
   performs motif and gene ontology analysis with lists of Gene Identifiers, both promoter and mRNA motifs (See Gene ID Analysis Tutorial)
   .pl 说明是HOMER里的perl的脚本。

   findGO.pl
   performs only gene ontology analysis with lists of Gene Identifiers (Called by findMotifs.pl, See Gene Ontology Analysis)
   这里是个findGO功能，不过我更常用的是enrichR 或者 DAVID。以上两个脚本都是gene ID based的，只需要准备个文本格式的gene list就也可以使用了。

2. Next-Gen Sequencing/Genomic Position Analysis
   findMotifsGenome.pl
   performs motif analysis from genomic positions (See Finding Motifs from Peaks)
   这个是通过基因组里peak的位置来找DNA motif，比较常用，因为根据测序方法不同，有些peak是在 non-coding promoter 或者 intergenic 等地方（也就是不只在coding gene promoter 的peak）。
   example:

$ cd /Users/ye.liu/Desktop/OA_analysis_06/9_patients_downstream_analysis/2.data_cpm2_p7/DNA_motif/Homer/1.complete_enhancer_promoter_sets/data
$ findMotifsGenome.pl 1.tss_gained_DAPs_gene_189.txt.bed  hg38 ./5.differential_output_size_400_1_to_3/ -bg 3.tss_lost_DAPs_gene_608.txt.bed -S 25 -len 8,10,12,13 -size 400
$ findMotifsGenome.pl 3.tss_lost_DAPs_gene_608.txt.bed  hg38 ./6.differential_output_size_400_3_to_1/ -bg $ 1.tss_gained_DAPs_gene_189.txt.bed -S 25 -len 8,10,12,13 -size 400 
$ findMotifsGenome.pl 2.tss_gained_DAPs_noncoding_91.txt.bed  hg38 ./7.differential_output_size_400_2_to_4/ -bg 4.tss_lost_DAPs_noncoding_509.txt.bed -S 25 -len 8,10,12,13 -size 400 
$ findMotifsGenome.pl 4.tss_lost_DAPs_noncoding_509.txt.bed  hg38 ./8.differential_output_size_400_4_to_2/ -bg 2.tss_gained_DAPs_noncoding_91.txt.bed -S 25 -len 8,10,12,13 -size 400  

这里是用的Differential ATAC-Peak （DAP）进行的motif查询，两组测序样品比较以后会得到gained DAPs和lost DAPs（样品组／对照组）。在DAP annotation的时候会有peak在coding／noncoding gene promoter （TSS）附近（上下1kb以内）就称它是gene associated with DAP=DAG，我用的是FANTOM CAT data set （2017 Nature） 进行的annotation，因为里面不但覆盖了coding gene 信息还同时有 noncoding gene 的信息。Intergenic 的DAP在这里我没有使用。所以我有四个bed file分别是:

	gained	lost
coding	file1_189	file3_608
non-coding	file2_91	file4_509

然后分别查找只在 gained DAG 里的 de novo DNA motif 和只在 lost DAG 里的 de novo DNA motif。关于background，我分别用对应的bed file来做背景peaks。
所以，
file1 比 file3 得到了 file5: DNA motif 只在 gained coding DAP而不在 lost coding DAP里。（反之得到 file6）
file2 比 file4 得到了 file7: DNA motif 只在 gained non-coding DAP而不在 lost non-coding DAP里。 （反之得到 file8）
file1-4 是指的bed file 5-8是HOMER的output。

---

接下来想要想要比较的只有DAP gain 与 DAP lost，不包括coding 和 noncoding。
所以需要做的事情是把file1 与 file 2结合起来变成 DAP gain
file3 与 file 4 结合起来就是 DAP lost。
之前会用比较笨的方法，bed file的 .bed改名成 .txt，打开复制粘贴到excel然后合并，保存称为.txt （用mac的要保存为windows的txt格式），再改名.bed,还会用到命令 changeNewLine.pl不然是个假的bed文件。
后来知道还有其他方法，linux command line：

$ cat 1.tss_gained_DAPs_gene_189.txt.bed  2.tss_gained_DAPs_noncoding_91.txt.bed > gained_DAP.bed

这么快吗？ rbind了？
检查一下，看看file1 和file2 分别有多少行（row）

$ cat 1.tss_gained_DAPs_gene_189.txt.bed |wc -l
 188
$ cat 2.tss_gained_DAPs_noncoding_91.txt.bed |wc -l
  90

那么合并后的文件应该就是188+90，这么多行了

$ cat gained_DAP.bed |wc -l
 278
#另一种方法
$ wc -l < gained_DAP.bed
 278

再不放心就检查一下，在terminal里查看下bed file。
方法1: cat file 全部输出
方法2: head -n 5 file or tail -n 6 file局部输出

$ head -n 10 gained_DAP.bed 
chr10   110460031   110460730   ENSG00000273143.1       RP11-525A16.4
chr20   58622490    58623170    ENSG00000268941.1       MGC4294
chr5    174750778   174752030   ENSG00000266890.1       MIR4634
chr17   18985476    18985916    ENSG00000263045.1       RP11-28B23.1
chr16   1163540 1164037 ENSG00000259910.1       RP11-616M22.2
chr12   46524079    46525144    ENSG00000257496.1       RP11-474P2.4
chr9    129740242   129741064   ENSG00000255824.1       AL590369.1
chr11   132874516   132875011   ENSG00000255371.1       OPCML-IT2
chr8    27901080    27902479    ENSG00000253615.1       RP11-597M17.2
chr8    66176914    66178013    ENSG00000253138.1       LINC00967

接下来同样办法得到 lost_DAP.bed

$ cat 3.tss_lost_DAPs_gene_608.txt.bed 4.tss_lost_DAPs_noncoding_509.txt.bed > lost_DAP.bed
$ wc -l < lost_DAP.bed                       
    1017

准备好了bed file后，开始进行motif查找，

$ pwd
/Users/ye.liu/Desktop/OA_analysis_06/9_patients_downstream_analysis/2.data_cpm2_p7/DNA_motif/Homer/1.complete_enhancer_promoter_sets/data/test
$ findMotifsGenome.pl gained_DAP.bed hg38 ./Gained_DAP_specific_motif_size_400/ -bg lost_DAP.bed -S 25 -len 8,10,12,13 -size 400
$ findMotifsGenome.pl lost_DAP.bed hg38 ./Lost_DAP_specific_motif_size_400/ -bg gained_DAP.bed -S 25 -len 8,10,12,13 -size 400

每一个会用掉30-40min这样。