整理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
完。