整理ChIP-seq / CUT & Tag 分析時(shí)用到的工具矾柜。本文只對(duì)使用的工具用法進(jìn)行簡(jiǎn)單介紹荣瑟。
MACS2 是一款常用的peak calling軟件，可以鑒定ChIP-seq/Cut&Tag數(shù)據(jù)的peaks零远，本文簡(jiǎn)單介紹MACS2 的安裝及peak calling的用法叙凡。
安裝
MACS2 可以通過(guò)conda、pip的方法進(jìn)行安裝咒精，也可以下載其源文件進(jìn)行安裝镶柱。以下展示conda的安裝方法
# --prefix 指定安裝到`ChIPseq`環(huán)境中
$ conda install -c bioconda macs2 --prefix=ChIPseq
用法
MACS2 有多個(gè)子命令，這里我們只介紹用于peak calling的callpeak
macs2 callpeak的示例用法如下
macs2 callpeak -t treatment.bam \
-c control.bam \
-f BAM \
-g hs \
-n <output_prefix>
--outdir <out_dir> \
--min-length 100 \
--nomodel --extsize 200 2>macs.log
-t：這是MACS唯一需要的參數(shù)模叙，指明處理組的文件歇拆，如果有多個(gè)可以用空格隔開(kāi)輸入-t A B C.如果沒(méi)有對(duì)照組，也可以單獨(dú)使用-t進(jìn)行peak calling向楼，但假陽(yáng)性率會(huì)有所提升查吊。
-c：對(duì)照組的文件，通常是input或者mock IP
-n：輸出結(jié)果的前綴
-f：輸入文件的類型湖蜕，可以是：ELAND, BED, ELANDMULTI, ELANDEXPORT, SAM, BAM, BOWTIE, BAMPE, or BEDPE格式逻卖。默認(rèn)自動(dòng)檢測(cè)，但如果是BAMPE或BEDPE格式需要手動(dòng)聲明昭抒。
-g：有效基因組大小评也，即目前可被測(cè)到的基因組大小炼杖。
MACS給出了幾個(gè)內(nèi)置的基因組大小，如果callpeak的物種不在此列需要給出相應(yīng)物種的有效基因組大小盗迟。
- hs: 2.7e9
- mm: 1.87e9
- ce: 9e7
- dm: 1.2e8
--outdir：輸出的文件夾
--min-length：peak最短的長(zhǎng)度坤邪，只有大于該長(zhǎng)度的peak才被認(rèn)為是一個(gè)peak。
--max-gap：定義了兩個(gè)相鄰peaks最大的距離罚缕，如果兩個(gè)peaks之間的距離小于該值則被merge為一個(gè)peak艇纺。
--nomodel：輸入這個(gè)參數(shù)后，MACS不會(huì)使用shifting model邮弹。
--extsize：使用--nomodel參數(shù)后黔衡，MACS會(huì)使用該參數(shù)對(duì)reads從5'->3' 延伸到指定的長(zhǎng)度。例如腌乡，TF的結(jié)合區(qū)域是200bp左右時(shí)盟劫，可以通過(guò)設(shè)置--nomodel --extsize 200來(lái)指定reads延伸的長(zhǎng)度。
-B/--bdg：設(shè)置該參數(shù)后与纽，MACS會(huì)用bedGraph存儲(chǔ)片段累積值（NAME_treat_pileup.bdg）和local lambda（NAME_control_lambda.bdg）侣签。
關(guān)于shifting model的一個(gè)解釋可以參考https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html
當(dāng)需要進(jìn)行寬峰的peak calling時(shí)，可以使用以下參數(shù)
--broad：使用該參數(shù)后急迂，MACS進(jìn)行寬峰的callpeak影所。
--broad-cutoff：寬峰的統(tǒng)計(jì)值cutoff。默認(rèn)使用q值過(guò)濾不顯著的峰袋毙，如果加了-p就用p值過(guò)濾型檀。
示例用法
# 常規(guī)轉(zhuǎn)錄因子ChIP-seq的peak calling:
$macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
# 組蛋白修飾類寬峰ChIP-seq的peak calling:
$macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
# 雙端測(cè)序的ATAC-seq的peak calling:
$macs2 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01
輸出結(jié)果
非寬峰
以非寬峰callpeak結(jié)果為例，通常MACS2 會(huì)輸出以下文件
1. NAME_peaks.xls
該文件可以通過(guò)excel打開(kāi)听盖，包括以下信息：
- 染色體名稱
- peak起始的染色體坐標(biāo)
- peak結(jié)束的染色體坐標(biāo)
- peak的長(zhǎng)度
- 峰頂?shù)慕^對(duì)坐標(biāo)
- 峰頂處的reads累積數(shù)
- -log10(pvalue) for the peak summit (e.g. pvalue =1e-10, then this value should be 10)
- 峰頂處reads相對(duì)于local background的富集倍數(shù) (fold enrichment for this peak summit against random Poisson distribution with local lambda)
- -log10(qvalue) at peak summit
2. NAME_peaks.narrowPeak
NAME_peaks.narrowPeak是一個(gè)BED6+4格式的文件胀溺，其中每列的含義如下：
以一個(gè)test.narrowPeak 為例
$head test.narrowPeak -n 2
GL000008.2 3142 3399 test_peak_1 50 . 3.81679 5.1 0.995095 23
GL000194.1 1917 2116 test_peak_2 56 . 4.41176 5.65309 1.289 106
3. NAME_summits.bed
該文件為BED5格式，包含了每個(gè)peaks峰頂?shù)奈恢眯畔⒔钥矗捎糜趍otif和binding sites的預(yù)測(cè)仓坞。
$head test_summits.bed -n 2
GL000008.2 3165 3166 test_peak_1 5.1
GL000194.1 2023 2024 test_peak_2 5.65309
注意的是"NAME_summits.bed:的第五列score與"NAME_peaks.narrowPeak"的第8列（峰頂處的-log10pvalue）一致。
寬峰
如果是寬峰的callpeak腰吟，則產(chǎn)生以下的peak calling結(jié)果：
NAME_peaks.broadPeak
該文件為BED6+3格式无埃，沒(méi)有了 .narrowPeak文件的第10列，即峰頂?shù)奈恢妹汀Ｓ捎赽road peak并不存在峰頂?shù)暮x嫉称，所以這個(gè)文件的第5，7-9列都是整個(gè)peak的平均值（而narrowPeak是峰頂處的值）灵疮。
$head test.broadPeak
chr2L 5129 9425 test_peak_1 67 . 2.56825 6.77836 5.42387
chr2L 365789 388188 test_peak_2 115 . 3.23359 11.5768 10.0304
Blacklist regions
需要注意的是由于NGS測(cè)序讀長(zhǎng)的限制织阅，目前的參考基因組有一些測(cè)不準(zhǔn)的區(qū)域，即黑名單區(qū)域（Blacklist regions）震捣。這些區(qū)域的有時(shí)候會(huì)具有高信號(hào)的富集荔棉，影響我們peak calling的結(jié)果闹炉。
為了提高peak calling的質(zhì)量，我們可以在peak calling結(jié)束后润樱，手動(dòng)去除這些區(qū)域渣触。通過(guò)bedtools這一工具就可以實(shí)現(xiàn)對(duì)blaklist region的去除。
$bedtools intersect -v -a NAME_peaks.narrowPeak -b BLACKLIST.BED > FILTERED.narrowPeak
Blacklist region的bed 文件可以在這里下載：https://github.com/Boyle-Lab/Blacklist/tree/master/lists
目前v2版好像只提供了人壹若、小鼠嗅钻、果蠅和線蟲(chóng)的。
對(duì)MACS2 Peak calling的簡(jiǎn)單介紹到此結(jié)束了舌稀。MACS3已經(jīng)更新了啊犬，大體參數(shù)和使用方法與MACS2區(qū)別不大，但具體的區(qū)別在哪還沒(méi)有深究壁查，以后有機(jī)會(huì)再作更新。
ref
MACS3 callpeak doc: https://github.com/macs3-project/MACS/blob/master/docs/callpeak.md
https://pypi.org/project/MACS2/#description
MACS2 peak calling tutorial: https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html
Blacklist regions: https://www.biostars.org/p/184537/
ENCODE paper about blacklist regions: https://www.nature.com/articles/s41598-019-45839-z
完剔应。
