學(xué)習(xí)目標(biāo)

1. 描述MACS2call peak算法的不同組成部分。
2. 描述運(yùn)行MACS2所涉及的參數(shù)。
3. 列出并描述MACS2的輸出文件。

call peak

call peak是我們工作流程中的下一步，它是一種計(jì)算方法，用于識(shí)別在進(jìn)行ChIP測(cè)序?qū)嶒?yàn)后基因組中比對(duì)reads富集的區(qū)域。

ChIP-seq

對(duì)于 ChIP-seq 實(shí)驗(yàn)，我們從比對(duì)文件中觀察到的是正負(fù)鏈上的 reads 密度存在鏈不對(duì)稱(chēng)性，且以結(jié)合位點(diǎn)為中心。所選片段的 5' 端將在正鏈和負(fù)鏈上形成組。然后使用統(tǒng)計(jì)測(cè)量評(píng)估這些組的分布，并與背景（Input或 IgG 樣本）進(jìn)行比較，以確定富集位點(diǎn)是否可能是真正的結(jié)合位點(diǎn)。

圖片來(lái)源：Wilbanks 和 Facciotti，PLoS One 2010

一個(gè)常見(jiàn)的問(wèn)題是我們?nèi)绾卧诒姸嗟腸all peak選項(xiàng)中進(jìn)行選擇，以及如何確定哪種方法最適合我們的數(shù)據(jù)？Wilbanks 和 Facciotti（2010）在 PLoS ONE 上進(jìn)行的一項(xiàng)研究對(duì) 12 種不同的call peak器在 ChIP-Seq 峰值檢測(cè)中的算法性能進(jìn)行了評(píng)估。結(jié)論是，雖然它們之間存在一定的一致性，但每個(gè)call peak軟件識(shí)別出的峰值數(shù)量差異很大。call peak軟件的選擇也取決于你的數(shù)據(jù)和所研究的蛋白質(zhì)。一種常用的識(shí)別結(jié)合位點(diǎn)的工具叫做基于模型的 ChIP-seq 分析（MACS），這也是我們?cè)诒窘坛讨袑⒁褂玫墓ぞ摺?/span>

CUT&RUN

雖然像 MACS2 這樣的標(biāo)準(zhǔn) ChIP-seq 通常用于從 CUT&RUN 數(shù)據(jù)中call peak，但人們擔(dān)心低讀取深度和低背景水平可能會(huì)使標(biāo)準(zhǔn)call peak器容易出現(xiàn)更多的假陽(yáng)性。為了解決這個(gè)問(wèn)題，Henikoff 小組開(kāi)發(fā)了一個(gè)名為 SEACR（用于 CUT&RUN 的稀疏富集分析）的工具，它提供了一種分析策略，該策略使用背景信號(hào)的全局分布來(lái)校準(zhǔn)一個(gè)簡(jiǎn)單的call peak閾值。

SEACR是如何工作的？

1. 首先，數(shù)據(jù)被解析為信號(hào)塊，這些信號(hào)塊表示由跨片段的 read 對(duì)形成的連續(xù)、非零讀取深度的片段。
2. 通過(guò)對(duì)每個(gè)塊中的 read 計(jì)數(shù)求和來(lái)計(jì)算信號(hào)。
3. 憑經(jīng)驗(yàn)確定一個(gè)閾值進(jìn)行過(guò)濾。 繪制目標(biāo)/IgG 中信號(hào)塊的比例（y 軸），用于確定使目標(biāo)與 IgG 塊的百分比最大化的閾值。
4. 通過(guò)過(guò)濾但與 IgG 塊重疊的富集區(qū)域也被移除。

圖片來(lái)源：“用于 CUT&RUN 染色質(zhì)分析的稀疏富集分析call peak”

ATAC-seq

ATAC-seq 的目標(biāo)是識(shí)別可及染色質(zhì)區(qū)域，并通過(guò)代理識(shí)別調(diào)控元件和轉(zhuǎn)錄因子結(jié)合位點(diǎn)。因此，call peak代表識(shí)別基因組中比對(duì) reads 富集的區(qū)域，這與我們對(duì) ChIP-seq 所做的類(lèi)似。目前，MACS2 是 ENCODE ATAC-seq 流程的默認(rèn)call peak器。與 ChIP-seq 分析工作流程相比，有幾個(gè)參數(shù)需要更改，我們將在本課程的末尾詳細(xì)描述它們。

具體來(lái)說(shuō)，我們需要考慮的差異包括： 缺乏Input樣本（陰性對(duì)照）、缺乏雙峰 read 分布（即不需要移動(dòng) reads）。

注意：還有其他 ChIP-seq 工具具有適應(yīng) ATAC-seq 數(shù)據(jù)的功能（即 Genrich，或者有專(zhuān)門(mén)為 ATAC-seq 設(shè)計(jì)的調(diào)用器（即 HMMRATAC）

MACS3

MACS 算法捕捉基因組復(fù)雜性的影響，以評(píng)估富集的 ChIP 區(qū)域的顯著性。雖然它是為檢測(cè)轉(zhuǎn)錄因子結(jié)合位點(diǎn)（窄峰）而開(kāi)發(fā)的，但它也適用于較大的區(qū)域（寬峰）。

我們將在本教程中使用 MACS3。call peak的底層算法與原始 MACS 2相同，但它在功能上有一些增強(qiáng)。MACS/MACS2/MACS3 工作流程如下圖所示。在本課程中，我們將更詳細(xì)地描述這些步驟。

去除冗余

MACS 為處理完全相同位置的重復(fù)標(biāo)簽（即具有相同坐標(biāo)和相同鏈的標(biāo)簽）提供了不同的選項(xiàng)。默認(rèn)是在每個(gè)位置保留一個(gè) read。非常常用的“auto”選項(xiàng)告訴 MACS 根據(jù)二項(xiàng)分布計(jì)算完全相同位置的最大標(biāo)簽數(shù)，使用 1e-5 作為 p 值截止。另一個(gè)選擇是設(shè)置“all”選項(xiàng)，它保留每個(gè)標(biāo)簽。如果指定一個(gè)整數(shù)，則在同一位置最多保留該數(shù)量的標(biāo)簽。這種冗余性在 ChIP 和Input樣本中始終一致地應(yīng)用。 我們不需要擔(dān)心這個(gè)選項(xiàng)，因?yàn)槲覀冊(cè)诒葘?duì)后過(guò)濾步驟中已經(jīng)過(guò)濾掉了重復(fù)項(xiàng)。

ChIP-seq 數(shù)據(jù)的雙峰性質(zhì)

下面，我們展示了感興趣的蛋白質(zhì)和從免疫沉淀中獲得的 DNA 片段（綠色）。

由于這些片段通常從 5' 端進(jìn)行測(cè)序，我們獲得的 reads 不會(huì)給出我們上面圖像中所示的片段堆積。相反，我們?cè)诘鞍踪|(zhì)的兩側(cè)（正鏈和負(fù)鏈上）獲得 read 堆積。

在將 reads 比對(duì)到基因組后，真正結(jié)合位點(diǎn)周?chē)?read 密度應(yīng)該顯示出雙峰富集模式（或成對(duì)的峰）。

shift size

MACS 利用這種雙峰模式憑經(jīng)驗(yàn)建模移動(dòng)大小，從而更好地定位精確的結(jié)合位點(diǎn)。

確定shift size

1. MACS 掃描整個(gè)樣本，搜索所有高度顯著的富集區(qū)域。這僅使用 ChIP 樣本完成！ 這些區(qū)域是由 MACS 在基因組上滑動(dòng)一個(gè) 600bp 的窗口來(lái)找到相對(duì)于隨機(jī)標(biāo)簽基因組分布具有超過(guò) 50 倍富集的標(biāo)簽區(qū)域來(lái)確定的。 注意 1：上面描述的窗口大小和富集倍數(shù)值是默認(rèn)值。雖然有參數(shù)可以讓你修改這些值（即 bw 和 mfold），但不建議調(diào)整它們。 注意 2：MACS2 的默認(rèn)富集倍數(shù)大于上面描述的 MACSv1 工作流程中的值
2. MACS 從在步驟 1 中確定的這些高質(zhì)量峰值中隨機(jī)抽取 1000 個(gè)。
3. 對(duì)于這 1000 個(gè)峰值，MACS 將它們的正鏈和負(fù)鏈標(biāo)簽分開(kāi)，并通過(guò)它們中心之間的中點(diǎn)對(duì)齊它們。對(duì)齊中兩個(gè)峰的模式之間的距離定義為“d”，代表估計(jì)的片段長(zhǎng)度
4. MACS 將樣本中的所有 reads 向 3' 端移動(dòng) d/2，移向最可能的蛋白質(zhì)-DNA 相互作用位點(diǎn)

縮放文庫(kù)大小

對(duì)于Input和處理樣本之間序列深度不同的實(shí)驗(yàn)，MACS 線性地將總對(duì)照標(biāo)簽計(jì)數(shù)縮放到與總 ChIP 標(biāo)簽計(jì)數(shù)相同。默認(rèn)行為是將較大的樣本向下縮放

有效基因組長(zhǎng)度

為了計(jì)算λbg（下面討論的一個(gè)參數(shù)），MACS 需要有效基因組大小或可映射的基因組大小?？捎成湫耘c基因組中特定位置的 k-mers 的唯一性有關(guān)。低復(fù)雜性和重復(fù)區(qū)域具有低唯一性，這意味著低可映射性。因此，我們需要提供有效基因組長(zhǎng)度來(lái)校正低可映射區(qū)域中真實(shí)信號(hào)的損失。

我如何獲得有效基因組長(zhǎng)度？ MACS 軟件為常用的生物體（人類(lèi)、小鼠、蠕蟲(chóng)和果蠅）預(yù)先計(jì)算了一些值。如果你愿意，你可以根據(jù)你的生物體和構(gòu)建計(jì)算更準(zhǔn)確的值。deepTools 文檔有針對(duì)更新的構(gòu)建的其他預(yù)先計(jì)算的值，并且也有一些關(guān)于如何計(jì)算它的好材料。

峰值檢測(cè)

在 MACS 將每個(gè)標(biāo)簽移動(dòng) d/2 之后，它然后使用 2d 的窗口大小在基因組上滑動(dòng)以找到候選峰值。基因組上的標(biāo)簽分布可以用泊松分布建模。泊松是一個(gè)單參數(shù)模型，其中參數(shù)λ是該窗口中預(yù)期的 reads 數(shù)量。它僅使用Input對(duì)照樣本計(jì)算。

MACS 為每個(gè)候選峰值計(jì)算一個(gè)λlocal。λlocal 參數(shù)是通過(guò)為不同的窗口大小計(jì)算λ值（如下所示）推導(dǎo)出來(lái)的。從這些值中，保留最大值來(lái)表示λlocal。

λlocal = MAX(λ300bp, λ1kb, λ5kb, λ10kb, λbg)。 λbg 表示使用整個(gè)可映射基因組（即最大窗口大?。┕烙?jì)的背景λ。

通過(guò)這種方式，λ捕捉了局部偏差的影響，并且對(duì)于小局部區(qū)域偶爾的低標(biāo)簽計(jì)數(shù)具有魯棒性。這些偏差的可能來(lái)源包括局部染色質(zhì)結(jié)構(gòu)、DNA 擴(kuò)增和測(cè)序偏差以及基因組拷貝數(shù)變異。

接下來(lái)，基于λ計(jì)算泊松分布 p 值。如果 p 值<1e-5，則認(rèn)為一個(gè)區(qū)域具有顯著的標(biāo)簽富集。任何重疊的富集峰值被合并為一個(gè)單一峰值。

錯(cuò)誤發(fā)現(xiàn)率

每個(gè)峰值被視為一個(gè)獨(dú)立的測(cè)試。因此，當(dāng)我們?cè)谝粋€(gè)樣本中遇到數(shù)千個(gè)檢測(cè)到的顯著峰值時(shí)，我們就有一個(gè)多重測(cè)試問(wèn)題。在 MACSv1.4 中，F(xiàn)DR 是通過(guò)交換 ChIP 和對(duì)照樣本來(lái)憑經(jīng)驗(yàn)確定的。然而，在 MACS2/3 中，p 值現(xiàn)在使用 Benjamini-Hochberg 校正進(jìn)行多重比較校正。

其他call peak軟件： 還有許多其他工具能夠處理兩種類(lèi)型的分布（例如窄峰、寬峰）；每個(gè)工具都有特定的子命令和/或模式來(lái)實(shí)現(xiàn)。因此，在選擇要運(yùn)行的call peak器時(shí)，了解你期望的結(jié)合分布類(lèi)型是很有好處的。    

HOMER：一套用于call peak和基序發(fā)現(xiàn)的工具    

SPP：一個(gè) R 包，在 ENCODE 處理流程中實(shí)現(xiàn)。最適合窄峰call peak。 使用滑動(dòng)窗口根據(jù)上下游側(cè)翼窗口中的片段計(jì)數(shù)計(jì)算分?jǐn)?shù)。    

epic2：理想的寬峰調(diào)用（對(duì)一個(gè)較舊的工具 SICER 的重新實(shí)現(xiàn)）    

haystack bio：表觀遺傳變異性和基序分析管道。    

運(yùn)行MACS2

MACS2參數(shù)

MACS2 中有七個(gè)主要功能可用作子命令。在本課程中我們僅涵蓋 callpeak，但如果你感興趣，可以使用 macs3 COMMAND -h 了解更多信息。

callpeak 是 MACS2 中的主要功能，可以通過(guò)輸入 macs2 callpeak 來(lái)調(diào)用。如果你輸入這個(gè)命令而不指定參數(shù)，你將看到命令行選項(xiàng)的完整描述。以下是常用選項(xiàng)的較短列表：

輸入文件選項(xiàng)

• -t：ChIP 數(shù)據(jù)文件（這是 MACS 唯一必需的參數(shù)）
• -c：對(duì)照或模擬數(shù)據(jù)文件
• -f：輸入文件的格式；默認(rèn)是“AUTO”，這將允許 MACS 自動(dòng)決定格式
• -g：可映射的基因組大小，定義為可以測(cè)序的基因組大小（1.0e+9 或 1000000000 都是可接受的格式）。提供了一些預(yù)編譯的值（例如，“hs”表示人類(lèi)（2.7e9），“mm”表示小鼠（1.87e9），“ce”表示秀麗隱桿線蟲(chóng)（9e7），“dm”表示果蠅（1.2e8））。 注意：雖然 MACS 可以在沒(méi)有Input對(duì)照的情況下調(diào)用峰值，但我們不建議這樣做。對(duì)照樣本增加了峰值調(diào)用的特異性，沒(méi)有它，你會(huì)發(fā)現(xiàn)許多假陽(yáng)性峰值被識(shí)別出來(lái)。

輸出參數(shù)

• --outdir：MACS2 將把所有輸出文件保存到這個(gè)選項(xiàng)指定的文件夾中
• -n：輸出文件的前綴字符串
• -B/--bdg：將片段堆積、對(duì)照 lambda、-log10pvalue 和 -log10qvalue 分?jǐn)?shù)存儲(chǔ)在 bedGraph 文件中

移位模型參數(shù)

• --nomodel：是否構(gòu)建移位模型。對(duì)于 ATAC-seq 峰值調(diào)用設(shè)置為 True
• --bw：帶寬，僅用于模型構(gòu)建時(shí)掃描基因組。不建議調(diào)整這個(gè)參數(shù)
• --mfold：模型構(gòu)建的上限和下限（默認(rèn)為 5 和 50）。不建議調(diào)整這個(gè)參數(shù)

峰值調(diào)用參數(shù)

• -q：峰值檢測(cè)的 q 值（最小錯(cuò)誤發(fā)現(xiàn)率）截止值。默認(rèn)設(shè)置為 0.05，不考慮 p 值
• -p：p 值截止值。這與-q 互斥。在需要更寬松的閾值時(shí)使用（見(jiàn)下面的注釋?zhuān)?。如果設(shè)置了 p 值截止值，q 值將不會(huì)被計(jì)算，并在最終的.xls 文件中報(bào)告為-1
• --nolambda：在峰值候選區(qū)域不考慮局部偏差/lambda
• --broad：寬峰調(diào)用

注意：在這種情況下，放寬 q 值并不像預(yù)期的那樣起作用，因?yàn)樗糠峙c峰值寬度相關(guān)。理想情況下，如果你放寬閾值，你只會(huì)得到更多的峰值。但是對(duì)于 MACS3，放寬閾值也會(huì)導(dǎo)致更寬的峰值。

現(xiàn)在我們對(duì)可以修改命令的不同方式有了一定的了解，讓我們?yōu)槊總€(gè)野生型重復(fù)樣本設(shè)置命令：

macs2 callpeak \
-t WT_1_sorted_final_filter.bam \
-c WT_1_Input_sorted_final_filter.bam \
-f BAM -g mm \
-n WT_1 \
--outdir macs2 2> WT_1.log

作為一個(gè)通用的峰值調(diào)用器，如果要問(wèn)的問(wèn)題僅僅是：“在哪里我們可以找到比隨機(jī)背景更顯著的 reads 覆蓋度？”那么 MACS2 可以應(yīng)用于任何 DNA 富集分析。下面，我們對(duì) CUT&RUN 和 ATAC-seq 數(shù)據(jù)的峰值調(diào)用所需的變化進(jìn)行說(shuō)明。

CUT&RUN 的參數(shù)如何變化？

對(duì)于 CUT&RUN-seq 數(shù)據(jù)的峰值調(diào)用，幾乎不需要進(jìn)行太多改變。唯一值得注意的區(qū)別是 CUT&RUN 測(cè)序數(shù)據(jù)通常是雙端的。為了考慮到這一點(diǎn)，你可以添加格式參數(shù)。    

• -f BAMPE

MACS2 中的雙端分析模式。在這種模式下，MACS2 正確地解釋測(cè)序 DNA 片段的完整范圍，并丟棄未正確配對(duì)的比對(duì)。 當(dāng)以單端模式分析雙端數(shù)據(jù)集時(shí)，MACS2 會(huì)去除每對(duì)中的第二個(gè) read（“R2” read），然后將剩余的“R1” reads 視為單端 reads。它從“單端”R1 reads 中建模片段長(zhǎng)度，然后將 read 長(zhǎng)度擴(kuò)展到模式的平均值。使用這種模式處理雙端數(shù)據(jù)可以使用實(shí)際的片段長(zhǎng)度，以獲得更準(zhǔn)確的最終結(jié)果。

ATAC-seq 的參數(shù)如何變化？

為了識(shí)別基因組中的可及區(qū)域，我們需要在過(guò)濾后獲得的無(wú)核小體 BAM 文件上調(diào)用峰值。目前，MACS2 是 ENCODE ATAC-seq 流程的默認(rèn)峰值調(diào)用器，因此下面我們提供如果使用 ATAC-seq 數(shù)據(jù)作為輸入時(shí)推薦的參數(shù)變化。

• -f BAMPE：MACS2 中的雙端分析模式  
• --nomodel：繞過(guò)構(gòu)建移位模型。read 堆積不代表雙峰模式，因?yàn)槲覀儧](méi)有在測(cè)定特定的蛋白質(zhì)-DNA 相互作用。開(kāi)放區(qū)域本質(zhì)上是單峰的，不需要對(duì) reads 進(jìn)行任何移位。
• --keep-dup all：保留所有 reads，因?yàn)槲覀円呀?jīng)從 BAM 文件中過(guò)濾了重復(fù)項(xiàng)。
• --nolambda：MACS2 將使用背景 lambda 作為局部 lambda（因?yàn)槲覀儧](méi)有用于 ATAC-seq 的input對(duì)照樣本）。

MACS輸出文件

對(duì)于每個(gè)樣本（2 個(gè)重復(fù)），應(yīng)該有 4 個(gè)文件輸出到結(jié)果目錄中，所以總共有 8 個(gè)文件：

• _peaks.narrowPeak：BED6+4 格式文件，其中包含峰值位置以及峰值頂點(diǎn)、p 值和 q 值
• _peaks.xls：一個(gè)表格文件，其中包含有關(guān)調(diào)用的峰值的信息。其他信息包括片段堆積和倍數(shù)富集（ChIP-seq 標(biāo)簽計(jì)數(shù)與λlocal 的比率）
• _summits.bed：峰值中片段堆積最高的位置。這些是預(yù)測(cè)的精確結(jié)合位置，建議用于基序發(fā)現(xiàn)
• _model.R：一個(gè) R 腳本，你可以使用它根據(jù)你的數(shù)據(jù)生成關(guān)于模型的 PDF 圖像和互相關(guān)圖。