前面给大家介绍了☛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
代码语言:javascript复制$ 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
代码语言:javascript复制# 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了
代码语言:javascript复制# 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相关信息可以参考:https://github.com/crazyhottommy/ChIP-seq-analysis/blob/master/part1_peak_calling.md