MACS2软件进行Peak calling

MACS2可以检出ChIP-Seq，ATAC-seq 以及 MeRIP-seq (RNA甲基化测序） 等富集的序列。

peaks calling 有不同的方法，MACS2是最常用的call peaks工具。 MACS全称Model-based Analysis of ChIP-Seq，最初的设计是用来鉴定转录因子的结合位点，但是它也可以用于其他类型的富集方式测序。

usage: macs2 [-h] [--version]             {callpeak,bdgpeakcall,bdgbroadcall,bdgcmp,bdgopt,cmbreps,bdgdiff,filterdup,predictd,pileup,randsample,refinepeak}

### macs2检出peaks

## 如果有多个实验组和对照组，先分别merge实验组和对照组的bam文件
samtools merge control.merge.bam c1_sorted.bam c2_sorted.bam
samtools merge chip.merge.bam chip1_sorted.bam chip2_sorted.bam 

## shift model for ATAC-seq
nohup macs2 callpeak -c control.merge.bam -t condition.merge.bam -f BAM -B -g hs -n test --nomodel --shift -100 --extsize 200 --outdir ./peaks 2> test.log & 

nohup macs2 callpeak -c control.merge.bam -t condition.merge.bam -f BAM -B -g hs -n test --outdir ./peaks 2> test.log &

# 设置阈值
macs2 callpeak -t condition.merge.bam -c control.merge.bam -f BAM -g hs -n sample5 --outdir ./peaks -B -q 0.01
 
macs2 callpeak -t condition.merge.bam -c control.merge.bam -f BAM --broad -g hs -n sample4  --outdir ./peaks --broad-cutoff 0.1

macs2 callpeak 主要参数介绍：

输入文件参数：

• -t:实验组，IP的数据文件
• c: 对照组
• f：指定输入文件的格式，默认是自动检测输入数据是什么格式，支持bam,sam,bed等
• g:有效基因组大小，由于基因组序列的重复性，基因组实际可以mapping的大小小于原始的基因组。这个参数要根据实际物种计算基因组的有效大小。软件里也给出了几个默认的-g 值：hs -- 2.7e9表示人类的基因组有效大小(UCSC human hg18 assembly).
  • hs: 2.7e9
  • mm: 1.87e9
  • ce: 9e7
  • dm: 1.2e8

输出文件参数：

• --outdir:输出结果的存储路径 -n:输出文件名的前缀
• -B/--bdg:输出bedgraph格式的文件，输出文件以NAME+'_treat_pileup.bdg' for treatment data, NAME+'_control_lambda.bdg' for local lambda values from control显示。bdg文件适合导入UCSC或IGV进行可视化，或者转换为bigwig格式再进行可视化。

peak calling 参数

• -q/--qvalue 和 -p/--pvalue q value默认值是0.05，与pvalue不能同时使用。
• --broad peak有narrow peak和broad peak, 设置时可以call broad peak 的结果文件。
• --broad-cutoff 和pvalue、以及qvalue相似
• --nolambda: 不要考虑在峰值候选区域的局部偏差/λ

Shift 模型参数：

• --nomodel 这个参数和extsize、shift是配套使用的，有这个参数才可以设置extsize和shift。
• --extsize 当设置了nomodel时，MACS会用--extsize这个参数从5'->3'方向扩展reads修复fragments。比如说你的转录因子结合范围200bp，就设置这个参数是200。
• --shift 当设置了--nomodel，MACS用这个参数从5' 端移动剪切，然后用--extsize延伸，如果--shift是负值表示从3'端方向移动。建议ChIP-seq数据集这个值保持默认值为0，对于检测富集剪切位点如DNAsel数据集设置为EXTSIZE的一半。 示例：

1. 想找富集剪切位点，如DNAse-seq，所有5'端的序列reads应该从两个方向延伸，如果想设置移动的窗口是200bp，参数设置如下： --nomodel --shift -100 --extsize 200
2. 对nucleosome-seq数据，用核小体大小的一半进行extsize,所以参数设置如下： --nomodel --shift 37 --extsize 73 --call-summits

MACS利用此参数重新分析信号谱，解析每个peak中包含的subpeak。对相似的结合图谱，推荐使用此参数，当使用此参数时，输出的subpeak会有相同的peak边界，不同的绩点和peak summit poisitions.

参考：

https://pypi.org/project/MACS2/

https://blog.csdn.net/weixin_44748341/article/details/112050912

https://cloud.tencent.com/developer/article/1346057