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格式。默认自动检测,但如果是BAMPEBEDPE格式需要手动声明。
-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

完。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,222评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,455评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,720评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,568评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,696评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,879评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,028评论 3 409
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,773评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,220评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,550评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,697评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,360评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,002评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,782评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,010评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,433评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,587评论 2 350

推荐阅读更多精彩内容