在諸如：ChIP-seq幽纷、ATAC-seq、MeRIP-seq(m6A)博敬、Cut&Tag中我們常常聽到一個(gè)詞Call Peak霹崎。Call Peak究竟是什么東西，得到的結(jié)果又有什么生物學(xué)意義呢冶忱？
概念
Call Peak其實(shí)是一種計(jì)算方法尾菇，用于識(shí)別基因組中由測序得到的比對reads富集的區(qū)域。
其實(shí)通俗的來說Call peak就是把我們有蛋白富集到核酸時(shí)的區(qū)域給找出來囚枪。
- 對于ChIPseq/Cut&Tag來說我們這個(gè)區(qū)域就是我們的轉(zhuǎn)錄因子所可以和DNA互作的區(qū)域派诬；
- 對于ATAC-seq來說這個(gè)就是特定條件下的開放染色質(zhì)的區(qū)域；
- 對于m6A測序來說就是發(fā)生m6A甲基化修飾的RNA區(qū)域链沼。
方法
為了找個(gè)蛋白結(jié)合核酸的區(qū)域默赂，我們需要進(jìn)行計(jì)算，因?yàn)槠鋵?shí)在通過測序我們得到的只不過是蛋白結(jié)合核酸片段的末端括勺，而非蛋白真正的結(jié)合區(qū)域缆八。所以就需要在得到的reads進(jìn)行延伸，如下圖：
理想情況下疾捍，我們認(rèn)為測得的reads最終其實(shí)會(huì)近似地平均分配到正負(fù)鏈上奈辰，這樣，對于一個(gè)結(jié)合熱點(diǎn)而言乱豆，reads在附近正負(fù)鏈上會(huì)近似地形成“雙峰”奖恰，我們偏移后就可以得到真正的結(jié)合區(qū)域了。
常用軟件
常用的Call Peak軟件很多宛裕，但總體上根據(jù)它們的推測peak的方法總體上可以分為3類：
1.基于馬爾科夫模型：HMMRATAC瑟啃；
2.基于count：MACS2、HOMER揩尸、F-seq等蛹屿；
3.基于峰形狀：PICS、PolyaPeak岩榆、CLC等错负。
不同的peak caller軟件有各自適用的范圍，但是理論上這些軟件都可以去做蛋白-核酸的call peak朗恳，對于各個(gè)軟件最適合的使用場景如下圖：
MACS2原理
MACS2是最常用的call peaks工具湿颅，是ENCODE的ChIP-seq和ATAC-seq流程采用的工具载绿。 MACS全稱Model-based Analysis of ChIP-Seq粥诫，該工具是大名鼎鼎的Liu tao大神開發(fā)。最初的設(shè)計(jì)是用來鑒定轉(zhuǎn)錄因子的結(jié)合位點(diǎn)崭庸，但是理論上它也可以用于其他的DNA-蛋白結(jié)合類型的富集方式測序分析怀浆。
軟件的計(jì)算流程如下：
- 去重：首先軟件會(huì)對數(shù)據(jù)進(jìn)行去重處理谊囚；
- 建模：前面在圖1中，我們談到了reads在附近正負(fù)鏈上會(huì)近似地形成“雙峰”执赡，這個(gè)雙峰之間的距離稱之為d镰踏，軟件需要進(jìn)行延伸需要知道這個(gè)d值。所以MACS2選取了1000個(gè)表達(dá)10~30倍的regions去建模以計(jì)算這個(gè)d值沙合。
- 延伸：得到了d之后奠伪，向3'端延伸就可以去進(jìn)行延伸，延伸長度為d/2首懈；
- 歸一化：把兩組數(shù)據(jù)進(jìn)行歸一化绊率，各組數(shù)據(jù)可比；
- 獲得候選的peak和背景進(jìn)行比較究履；
- 計(jì)算確定λ滤否，代入泊松分布，得到p值最仑；
補(bǔ)充知識(shí)：
（1）假設(shè)TF在基因組上的分布是沒有任何規(guī)律藐俺，測序得到的read在基因組上的分布也必然是隨機(jī)的，某個(gè)堿基上覆蓋的read的數(shù)目應(yīng)該服從二項(xiàng)分布泥彤。當(dāng)n很大欲芹，p很小時(shí)，二項(xiàng)分布可以近似用泊松分布替代吟吝。
（2）在MACS2中耀石，它并不是使用一個(gè)固定的λ值，而是采用了一個(gè)動(dòng)態(tài)規(guī)劃的思路爸黄，服從λlocal = max(λBG, λ1k, λ5k, λ10k)滞伟，泊松分布的公式和模型如下圖：
image.png
image.png
λ是滑動(dòng)窗口的期望值， n是測序得到的read總數(shù)目炕贵， l是單個(gè)read的長度梆奈，s是基因組的大小
從上述的公式中我們可以看到，n称开、l亩钟、s都是確定的，所以只需要確定了λ我們就可以得到P值鳖轰。如果只是隨機(jī)分布的片段他們會(huì)落在上述密度分布圖封頂?shù)膬啥饲逅郑瑑蛇叺母怕适菢O小的，MACS2使用了單邊檢驗(yàn)蕴侣，當(dāng)落在右邊的P達(dá)到了閾值時(shí)候焰轻，我們就認(rèn)為這個(gè)是我們要找的peak了。
- 計(jì)算FDR：當(dāng)有對照組存在時(shí)候昆雀，就可以得到兩個(gè)P辱志，F(xiàn)DR=P(對照組)/P(實(shí)驗(yàn)組)
實(shí)操
其實(shí)MACS2已經(jīng)更新到了MACS3了蝠筑，簡單介紹下使用方法：
尋找ChIP-seq的轉(zhuǎn)錄因子（TF）的命令：
macs3 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
尋找ChIP-seq的組蛋白（Histone ）Mark的命令：
macs3 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
對于ATAC-seq雙端數(shù)據(jù)的操作：
macs3 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01
-f是輸入文件的格式，可以是BAM揩懒、BED什乙，軟件可以自動(dòng)識(shí)別，但是需要注意自動(dòng)識(shí)別是無法區(qū)分是PE還是SE已球，所以建議手動(dòng)指定臣镣；
-g是有效基因組大小，軟件預(yù)設(shè)了模式物種的有效因組大小智亮，如果hs人退疫，mm小鼠；
-n是樣本名字鸽素；
-B是以bedGraph格式存放fragment pileup, control lambda, -log10pvalue 和log10qvale,使用會(huì)增長耗時(shí)褒繁；
-q就是用q值（FDR）進(jìn)行篩選，輸入得到值為閾值馍忽；
broad-cutoff用于過濾 broad得到的peak棒坏。
最后：如果想了解更多和生信或者精品咖啡有關(guān)的內(nèi)容歡迎關(guān)注我的微信公眾號(hào)：生信咖啡，更多精彩等你發(fā)現(xiàn)遭笋！
