MACS2 进行peak calling
整理ChIP-seq / CUT & Tag 分析时用到的工具。本文只对使用的工具用法进行简单介绍。
MACS2 是一款常用的peak calling软件，可以鉴定ChIP-seq/Cut&Tag数据的peaks，本文简单介绍MACS2 的安装及peak calling的用法。
安装
MACS2 可以通过conda、pip的方法进行安装，也可以下载其源文件进行安装。以下展示conda的安装方法
# --prefix 指定安装到`ChIPseq`环境中
$ conda install -c bioconda macs2 --prefix=ChIPseq
用法
MACS2 有多个子命令，这里我们只介绍用于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唯一需要的参数，指明处理组的文件，如果有多个可以用空格隔开输入-t A B C.如果没有对照组，也可以单独使用-t进行peak calling，但假阳性率会有所提升。
-c：对照组的文件，通常是input或者mock IP
-n：输出结果的前缀
-f：输入文件的类型，可以是：ELAND, BED, ELANDMULTI, ELANDEXPORT, SAM, BAM, BOWTIE, BAMPE, or BEDPE格式。默认自动检测，但如果是BAMPE或BEDPE格式需要手动声明。
-g：有效基因组大小，即目前可被测到的基因组大小。
MACS给出了几个内置的基因组大小，如果callpeak的物种不在此列需要给出相应物种的有效基因组大小。
- hs: 2.7e9
- mm: 1.87e9
- ce: 9e7
- dm: 1.2e8
--outdir：输出的文件夹
--min-length：peak最短的长度，只有大于该长度的peak才被认为是一个peak。
--max-gap：定义了两个相邻peaks最大的距离，如果两个peaks之间的距离小于该值则被merge为一个peak。
--nomodel：输入这个参数后，MACS不会使用shifting model。
--extsize：使用--nomodel参数后，MACS会使用该参数对reads从5'->3' 延伸到指定的长度。例如，TF的结合区域是200bp左右时，可以通过设置--nomodel --extsize 200来指定reads延伸的长度。
-B/--bdg：设置该参数后，MACS会用bedGraph存储片段累积值（NAME_treat_pileup.bdg）和local lambda（NAME_control_lambda.bdg）。
关于shifting model的一个解释可以参考https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html
当需要进行宽峰的peak calling时，可以使用以下参数
--broad：使用该参数后，MACS进行宽峰的callpeak。
--broad-cutoff：宽峰的统计值cutoff。默认使用q值过滤不显著的峰，如果加了-p就用p值过滤。
示例用法
# 常规转录因子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
# 双端测序的ATAC-seq的peak calling:
$macs2 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01
输出结果
非宽峰
以非宽峰callpeak结果为例，通常MACS2 会输出以下文件
1. NAME_peaks.xls
该文件可以通过excel打开，包括以下信息：
- 染色体名称
- peak起始的染色体坐标
- peak结束的染色体坐标
- peak的长度
- 峰顶的绝对坐标
- 峰顶处的reads累积数
- -log10(pvalue) for the peak summit (e.g. pvalue =1e-10, then this value should be 10)
- 峰顶处reads相对于local background的富集倍数 (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是一个BED6+4格式的文件，其中每列的含义如下：
以一个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格式，包含了每个peaks峰顶的位置信息，可用于motif和binding sites的预测。
$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，则产生以下的peak calling结果：
NAME_peaks.broadPeak
该文件为BED6+3格式，没有了 .narrowPeak文件的第10列，即峰顶的位置。由于broad peak并不存在峰顶的含义，所以这个文件的第5，7-9列都是整个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测序读长的限制，目前的参考基因组有一些测不准的区域，即黑名单区域（Blacklist regions）。这些区域的有时候会具有高信号的富集，影响我们peak calling的结果。
为了提高peak calling的质量，我们可以在peak calling结束后，手动去除这些区域。通过bedtools这一工具就可以实现对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版好像只提供了人、小鼠、果蝇和线虫的。
对MACS2 Peak calling的简单介绍到此结束了。MACS3已经更新了，大体参数和使用方法与MACS2区别不大，但具体的区别在哪还没有深究，以后有机会再作更新。
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
完。
