一、不求甚解系列

軟件下載及配置

1. conda安裝： conda install -c bioconda homer

2. 使用configureHomer.pl完成HOMER軟件的配置

# 下載配置文件
wget http://homer.salk.edu/homer/configureHomer.pl 
?
# 使用配置文件進(jìn)行軟件配置
perl configureHomer.pl -install
?
# 下載 hg19 人的參考基因組
# 將hg19替換為mm10可下載mm10 小鼠的參考基因組
# 因?yàn)檎n題組人和小鼠的研究居多，這里只以這兩者舉例
perl configureHomer.pl -install hg19 

Motif 分析及peak注釋

此處以 MACS2 軟件call peak得到的bed文件舉例：

HOMER軟件接受兩種文件格式，這里不展開介紹文件的具體格式，細(xì)節(jié)可看使用詳解系列

findMotifsGenome.pl H3K4Me3.bed hg19 H3K4Me3_motif -len 8,10,12
# peak文件：ChIP-Seq_H3K4Me3_1_homer.bed
# 基因組版本：hg19
# 輸出路徑：ChIP-Seq_H3K4Me3_1_motifDir
# motif長度：-len 8,10,12 
# motif的軟件默認(rèn)長度為8，10，12
?
annotatePeaks.pl H3K4Me3.bed hg19 1> ChIP-Seq_H3K4Me3_1_peakAnn.xls 2>H3K4Me3_annLog.txt </pre>

結(jié)果展示

軟件運(yùn)行后的結(jié)果全部放在我們制定的輸出文件路徑下，將該目錄下所有文件下載到本地，打開homerResults.html 文件，即可得到如下結(jié)果展示：

隨便找的一個(gè)例子，不一定是H3K4me3的motif

二、使用詳解

HOMER 一個(gè)常用的motif分析軟件。它通過比較兩個(gè)序列集，并使用ZOOPS scoring和超幾何分布（或者負(fù)二項(xiàng)分布）進(jìn)行motif的富集分析。它主要用于ChIP-seq和promoter分析，但也可以用于核酸序列的motif分析問題。

HOMER軟件可以進(jìn)行多種類型的motif分析，如 promoter motif analysis ，基因組位置motif分析(ChIP-seq分析中的motif分析)，利用自定義的fasta文件進(jìn)行motif分析，RNA序列的motif分析（分析CLIP-seq數(shù)據(jù)中的RNA binding elements）

HOMER進(jìn)行motif分析時(shí)，需要兩個(gè)數(shù)據(jù)集：

1. 感興趣的目標(biāo)序列，如ChIP-seq分析中的peak文件

2. 背景序列，如ChIP-seq分析中的物種全基因組序列

HOMER包的介紹與安裝

HOMER包分為四種：

1. SOFTWARE： 軟件包，里面包括所有的代碼以及一些常見的數(shù)據(jù)文件，如motif矩陣

2. ORGANISMS：物種特異性數(shù)據(jù)包，包括一些指定物種的accession conversion數(shù)據(jù)，gene description數(shù)據(jù)，GO 分析數(shù)據(jù)。大部分?jǐn)?shù)據(jù)都是以NCBI的gene數(shù)據(jù)庫為基礎(chǔ)的

3. PROMOTERS：啟動(dòng)子數(shù)據(jù)包，對promoter進(jìn)行motif分析時(shí)所需要的相關(guān)數(shù)據(jù)和promoter序列信息。大部分?jǐn)?shù)據(jù)都是基于RefSeq 轉(zhuǎn)錄組。

4. Genomes：基因組數(shù)據(jù)包，包括一些基因組序列和注釋信息

HOMER軟件的安裝可看不求甚解部分。HOMER軟件剛安裝之后，是不包含任何序列信息的。我們可以使用 configureHOMER.pl 腳本下載所需數(shù)據(jù)。 我們可以使用 perl configureHOMER.pl -list 命令得到已有的HOMER包。

有些軟件包的名稱后面會(huì)跟有 -o, -p, -g等參數(shù)，這是HOMER為了區(qū)分同一物種不同類型數(shù)據(jù)所添加的。如 human-p。其中 -o 表示organism，-p 表示promoter，-g 表示genome。 這里給出幾個(gè)軟件包安裝與刪除的例子：

# install 表示安裝
perl configureHomer.pl -install mouse # 安裝小鼠的promoter數(shù)據(jù)集
perl configureHomer.pl -install hg19 # 安裝人的hg19版本的基因組數(shù)據(jù)集
?
# remove 表示刪除
perl configureHomer.pl -remove mm10 # 刪除小鼠mm10版本的基因組數(shù)據(jù)

此外，需要注意的是，當(dāng)我們下載基因組數(shù)據(jù)時(shí)，同時(shí)會(huì)自動(dòng)下載物種包的數(shù)據(jù)（比如，下載hg19數(shù)據(jù)時(shí)，會(huì)自動(dòng)下載human數(shù)據(jù)）。每次我們下載某個(gè)物種的基因組數(shù)據(jù)或者promoter數(shù)據(jù)時(shí)，它都會(huì)自動(dòng)替我們檢查是否下載了對應(yīng)的物種數(shù)據(jù)。

motif分析

正如前面所說，HOMER可以進(jìn)行多種類型的motfi分析。這里我們以最常見的ChIP-seq流程中的peaks motfi分析舉例。該分析需要使用 findMotifsGeome.pl 腳本。

該腳本的輸入文件支持兩種類型的文件，一種是常見的bed文件格式(MACS2軟件call peak得到的bed文件即可直接使用)，另一種是HOMER軟件指定使用的peak文件格式。 下面詳細(xì)講一下這兩種文件格式:

1. bed文件格式：findMotifsGeome.pl 腳本用到的信息為bed文件的前六列，分別是 chr, start, end, peak ID, score(HOMER用不到該信息，但是必須有), strand 。

2. peak文件格式：該文件格式需要五列信息（使用Tab分隔），分別是 peak ID , chr , start , end ,strand 。

當(dāng)我們準(zhǔn)備分析一個(gè)peak/bed文件的motif時(shí)，需要運(yùn)行以下命令：

findMotifsGenome.pl <peak/BED file> <genome> <output directory> -size # [options]
# 比如： findMotifsGenome.pl ERpeaks.txt hg18 ER_MotifOutput -size 200 -mask

• -mask : 該參數(shù)告訴motif分析程序，在得到一個(gè)可能的motif之后，在后續(xù)的motif分析中是否排除該motif的影響。有點(diǎn)類似于抽樣調(diào)查中的無放回抽樣。

• -size: 指定用于motif分析的片段長度，默認(rèn)為200； -size given 設(shè)置片段大小為目標(biāo)序列長度。越大，motif分析所需要的計(jì)算資源越多。

下面介紹一些常用的其他參數(shù)：

• -p : 設(shè)置線程數(shù)-S : 所尋找的motif數(shù)目，默認(rèn)為25。25已經(jīng)不算少了，如果自定義數(shù)目，推薦設(shè)置更少而不是更多。

• -mis : 所允許的錯(cuò)配數(shù)，默認(rèn)為2

• -cpg : 在對背景序列和目標(biāo)序列進(jìn)行normalization時(shí)，使用CpG%矯正而非GC%。

• norevopp : 只在peak所在的鏈搜索motif，不進(jìn)行反義鏈搜索。HOMER默認(rèn)會(huì)在兩條鏈上都進(jìn)行搜索。

• -rna : 搜索RNA上的motif（如使用CLIP-seq數(shù)據(jù)進(jìn)行分析的時(shí)候）

• -bg : 指定自定義的背景序列。HOMER默認(rèn)隨機(jī)選取基因組序列作為背景。

• -h: findMotifsGenome.pl腳本默認(rèn)使用二項(xiàng)分布進(jìn)行motif富集分析，這在背景序列多于目標(biāo)序列時(shí)是十分有用的。但是有時(shí)我們使用 -bg 參數(shù)自定義背景序列時(shí)，其數(shù)目可能會(huì)小于目標(biāo)序列，此時(shí)使用 -h 參數(shù)選擇超幾何分布會(huì)更加適合。

• -len : motif的長度設(shè)置，默認(rèn)8，10，12，越大越消耗計(jì)算資源。

• -N: 用于motif分析時(shí)所需的序列數(shù)目，通常當(dāng)我們設(shè)置 -len過大，內(nèi)存不夠時(shí)，會(huì)選擇減小 -N參數(shù)或者 -size 參數(shù)。

三、原理

前面我們說到，HOMER軟件可以進(jìn)行多種類型的motif分析，如 promoter motif analysis ，基因組位置motif分析(ChIP-seq分析中的motif分析)，利用自定義的fasta文件進(jìn)行motif分析，RNA序列的motif分析（分析CLIP-seq數(shù)據(jù)中的RNA binding elements）。 但是，不管我們?nèi)绾握{(diào)用HOMER，對于發(fā)現(xiàn)motif的基本步驟都是一致的，下面主要介紹一下各步驟的原理：

原理理解即可，HOMER已經(jīng)將各個(gè)步驟封裝起來，不用一步一步分別進(jìn)行。

預(yù)處理

1. 序列提取

• 若給出的基因組位置信息，則提取出來的是對應(yīng)的基因組序列；

• 若給出的是基因accession number，給需要選擇適當(dāng)?shù)膒romoter區(qū)域

2. 背景選擇

如果沒有自定義背景序列信息（即： -bg 參數(shù)），HOMER將自動(dòng)為用戶選取。如果使用的是基因組位置，將從整個(gè)基因組序列抓取GC含量一致的序列作為背景序列。若進(jìn)行的是promoter分析，則將所有的promoter作為背景。

3. GC含量矯正

將目標(biāo)序列和背景序列對GC含量進(jìn)行分bin（5%區(qū)間），背景序列通過調(diào)節(jié)權(quán)重得到和目標(biāo)序列同樣的GC含量分布。這可以使得HOMER在分析來自CpG島時(shí)的序列時(shí)，不會(huì)簡單的找到那些GC富集的motif。下圖是一個(gè)ChIP-seq分析中GC含量分bin的例子

4. 自動(dòng)標(biāo)準(zhǔn)化

除了GC含量之外，目標(biāo)序列還會(huì)存在其他的情況影響分析結(jié)果。如外顯子上的密碼子偏好等生物學(xué)現(xiàn)象，由A富集序列優(yōu)先排列引起的實(shí)驗(yàn)偏差。當(dāng)這些偏差顯著時(shí)，將會(huì)影響HOMER的motif分析。HOMER通過給背景序列分配適當(dāng)?shù)臋?quán)重，減少寡核苷酸序列頻率（如AA）在背景與目標(biāo)的差異。

重頭預(yù)測motif

默認(rèn)情況下HOMER調(diào)用homer2版本，若想使用舊版本，在命令行中添加 -homer1 參數(shù)即可

5. 解析輸入序列

根據(jù)設(shè)置的motif長度，將輸入序列解析為寡核苷酸，并生成一個(gè)寡核苷酸頻度表。該表只記錄出現(xiàn)的寡核苷酸和其出現(xiàn)的次數(shù)，可以增加motif搜索時(shí)的效率，但是也會(huì)破壞寡核苷酸與其原始序列的關(guān)系。

6. 寡核苷酸自動(dòng)均一化

HOMER除了在第四步中對全局序列進(jìn)行了自動(dòng)均一化矯正，還在寡核苷酸頻度表中進(jìn)行了矯正。該方法的主要思想是將那些較短的寡核苷酸與較長的寡核苷酸相平衡，但是該方法由于存在許多與motif同長的寡核苷酸，需要分配數(shù)量眾多的權(quán)重。對于那些極端的序列偏好，該方法不會(huì)移除它們。

7. 全局搜索

在構(gòu)建好寡核苷酸頻度表后，HOMER進(jìn)行全局motif搜索。其基本原理就是若某個(gè)motif富集，則其存在的寡核苷酸也同樣富集。為了增加靈敏度，HOMER允許搜索motif時(shí)的錯(cuò)配，同時(shí)跳過多個(gè)錯(cuò)配的寡核苷酸節(jié)省計(jì)算資源。

8. 矩陣優(yōu)化

9. Mask and Repeat

當(dāng)最優(yōu)oligo被優(yōu)化成motif后，motif 對應(yīng)的序列從要分析的數(shù)據(jù)中移除，接下來再分析最優(yōu)的…..直到設(shè)定個(gè)數(shù)的motif被發(fā)現(xiàn)。

搜索已知的motif（homer2版本）

10. 加載moitf數(shù)據(jù)庫

HOMER從以前的數(shù)據(jù)中加載一個(gè)已知的motif列表從而搜索已知的motif還可以通過在命令行中指定的motif( -mknown)或通過編輯homer包中的文件( data/knownTFs/known.motifs )來添加motifs。

11. 掃描所有的motif

HOMER通過掃描所有的序列，確定每個(gè)motif的富集度，并計(jì)算其顯著性。

motif輸出的結(jié)果分析

12. Motif文件

HOMER的輸出有許多后綴名為 motif 的文件，其內(nèi)容如下：

每個(gè)motif信息是一塊，均以>開頭。>所在行的信息以tab分隔。 motif首行信息解釋：

• > + 一致性序列：如圖上的>ASTTCCTCTT

• Motif名稱：如圖上的1-ASTTCCTCTT

• 比值檢測概率的log值：8.059752

• P-value得log值：-23791.535714

• 占位符：如上圖得0，不具有任何信息

• 逗號分隔得富集信息，如：T:17311.0(44.36%),B:2181.5(5.80%),P:1e-10317

  • T表示帶有該motif的目標(biāo)序列在總的目標(biāo)序列（target）中得百分比

  • B表示帶有該motif的背景序列在總的背景序列（background）中得百分比

  • P表示最終富集的p-value

• 逗號分隔的motif統(tǒng)計(jì)信息，如：Tpos:100.7,Tstd:32.6,Bpos:100.1,Bstd:64.6,StrandBias:0.0,Multiplicity:1.13

  1. Tpos：motif在目標(biāo)序列中的平均位置

  2. Tstd：motif在目標(biāo)序列中位置的標(biāo)準(zhǔn)差

  3. Bpos：motif在背景序列中的平均位置

  4. Bstd：motif在背景序列中位置的標(biāo)準(zhǔn)差

  5. StrandBias: 鏈偏好性，在正義鏈上的motif數(shù)與反義鏈motif數(shù)的比值的log

  6. Multiplicity：具有一個(gè)或多個(gè)結(jié)合位點(diǎn)的序列中每個(gè)序列的平均出現(xiàn)次數(shù)

13 denovo預(yù)測的motif

首先對產(chǎn)生的motif進(jìn)行去冗余，然后將motif轉(zhuǎn)換為向量值從而計(jì)算皮爾遜相關(guān)系數(shù)。HOMER產(chǎn)生的motif結(jié)果使用網(wǎng)頁展示，該文件名稱為 homerResults.html ，而目錄homerResults/ 存放著所有該網(wǎng)頁依賴的文件和圖片。生成的motif首先去冗余。以下是一個(gè)簡單的輸出示例：

shi

輸出結(jié)果按照p-value排序，最后一列是一個(gè)鏈接到motif文件的超鏈接，可以從這個(gè)文件中找到包含此motif的其他序列。在Best Match/Details 列中，HOMER將會(huì)展示與denovo motif最匹配的已知motif（不可全信，最匹配不一定是最優(yōu)）。該列還包含有一個(gè)名為 More information 的超鏈接，點(diǎn)擊后會(huì)出現(xiàn)以下頁面：

該頁面包含motif的一些基本信息，如鏈接到motfi文件的超鏈接，且可以生成motif logo的pdf格式。

接著是與已知motif的match。這部分主要看看denovo motif和已知的motif的相似性，需要檢查一下是否合理，如下：

如果點(diǎn)擊 Simiar motifs 超鏈接，將會(huì)展示在denovo motif搜尋過程中與該motif相似的其他motifs（這些motif的enrichment value相對較低）。通常我們會(huì)點(diǎn)擊查看一下是否有一些motifs被不正確的分類，比如幾個(gè)motif具有幾個(gè)相同的堿基。具體界面如下：

14. 已知motif的輸出

該輸出也是以網(wǎng)頁形式展現(xiàn)，頁面結(jié)果根據(jù)enrichement進(jìn)行排序，具體如下：