MACS2软件call m6A甲基化peaks

前面给大家介绍了☛m6A甲基化数据分析流程，提到过两个peak calling的软件

Peak calling

• 用R包exomePeak call peak
• 用MACS2 call peak

其实目前可用的peak calling工具很多，如下表所示

更多信息也可以参考

http://wodaklab.org/nextgen/data/peakfinders.html

MACS2这个软件是目前用的比较多，也是大家比较认可的一款软件。MACS2是一个用python2.7写的工具，所以当你同时在使用python3.6和python2.7时，使用前请务必激活python2.7（将 python2.7/anaconda2 的安装目录添加到环境变量中），安装方法为：

下载安装anaconda2

$ wget -c -P basic_tool/ https://repo.continuum.io/archive/Anaconda2-5.0.1-Linux-x86_64.sh
$ sh Anaconda2-5.0.1-Linux-x86_64.sh
$ echo 'export PATH=../anaconda2/bin:$PATH' >>~/.bashrc

下载安装MACS2

# 1. 用源码安装
$ wget -c -P biosoft/ https://pypi.python.org/packages/9f/99/a8ac96b357f6b0a6f559fe0f5a81bcae12b98579551620ce07c5183aee2c/MACS2-2.1.1.20160309.tar.gz
$ cd biosoft  && tar zxvf MACS2-2.1.1.20160309.tar.gz
$ cd MACS2-2.1.1.20160309 && python setup.py install
$ echo 'export PATH=../MACS2-2.1.1.20160309/bin:$PATH' >>~/.bashrc

# 2. 用bioconda安装
$ conda install -c bioconda macs2

安装成功后就可以直接使用MACS2进行peak calling了

# MACS首先的工作是要确定一个模型，这个模型最关键的参数就是峰宽d，这个d就是bw(band width)，而它的一半就是shiftsize
$ macs2 callpeak -c controlFile.bam -t treatmentFile.bam -m 10 30 -p pvalue -f BAM -g gsize -B -n filename.preffix --outdir ChIP_seq/CallPeak 2>ChIP_seq/CallPeak/filename.macs2.log

输入文件参数：

• -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显示。

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的一半。

MACS2软件call m6A甲基化