ChIP-seq peak，我想call就call

目前最常用的鉴定转录因子结合位点的工具是刘小乐老师开发的「基于模型的芯片序列分析（Model-based Analysis of ChIP-seq, MACS）」。

MACS算法通过捕捉基因组复杂度的影响来评估富集的ChIP区域的显著性。MACS通过结合测序标签（ChIP-tag）的「位置」和「方向」信息来提高结合位点的空间分辨率。

MACS可直接用于单ChIP样本分析，当然最好是有对照样本一起使用。它的原理图相信很多人都看过，如下图所示，蓝色椭圆形部分为蛋白-DNA结合为止，由于超声打断DNA的位点还有大小不是固定，因此得到的DNA片段也不是一样的，在测序后得到的序列和基因组相比，是形成峰分布的样子。由于双端测序，则有双峰（即使是单端测序，也会做成双端来处理）：

1. MACS 安装

因为macs2用的python2，会和一些用python3的软件冲突，所以我们得为macs2单独建立一个环境

conda create -n macs2 -c conda-forge -c bioconda macs2 -y
source activate macs2

2. peak calling

2.1 MACS2用法

macs2有很多子命令，其中我们最关心的是callpeak

# 一个例子: for regular peak calling on TF ChIP-seq:
macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

# 本文中的代码，因为文献里指明q value 是0.05，所以这里我们设置为0.05

INDIR=03-read-align
OUTDIR=04-peak-calling

mkdir -p $OUTDIR
# 以ESC-1为例
conda activate macs2
macs2 callpeak -t $INDIR/ESC_TP53_1.bam -c $INDIR/ESC_input_1.bam -f BAM -g hs --outdir $OUTDIR -n ESC_TP53_1 -B -q 0.05

参数说明：

❝

(1) -t : t是treatment组的意思，意为实验组，在本文中即为P53 ChIP组，如果有多个样本，可以写为 -t A B C。

(2) -c ：c是control组的意思，即为实验组，如果有多个样本，同-t一样的用法。

(3) -n ：n是name的首字母，就是给输出文件设置名字。

(4) -f ：是format，即文件格式（ELAND，BED，SAM，BAM）等等。

(5) -g ：gsize，可用图表示的基因组大小，不同种属有不同的设定，一般人类的是 hs: 2.7e9。

(6) -q ：可设置qvalue阈值。

(7) -B ：bdg，以bedGraph格式存放fragment pileup, control lambda, -log10pvalue 和log10qvale。

❞

因为样本比较少，偷懒不写循环了，还是手动处理。

2.2 MACS2 peak calling输出数据

检查输出文件可见6个不同的文件：

其中：

（1）_control_lambda.bdg：一个不怎么重要的文件，之所以大，主要是记录各个区间的lambda值

（2）_model.R：是R代码文件，可用于绘图，得到的是基于所提供的数据的peak模型。

（3）_peaks.narrowPeak：BED6+4格式文件，包含峰的位置以及峰的峰高、pvalue和qvalue，这是另外一种作者常喜欢上传的数据格式；

（4）_peaks.xls ：包含关于被称为峰值的信息的xls文件。和narrowPeak的区别是它的坐标是从1开始的，而narrowPeak是从0开始的。

（5）_summits.bed ：包含每个峰定点的位置，如果要找在binding site的motif，建议使用这个文件。

（6）_treat_pileup.bdg：实验组的bedGraph格式。

2.3 查看输出文件

「xls文件」

拿到这么多文件，我第一个感兴趣的是**.xls「以及」_peaks.narrowPeak「文件，原因是这两个文件很小而且都见过。于是乎我先打开」.xls「文件查看，发现了一个很有意思的地方：有一个人胚胎干细胞组（ESC）样本没有peaks，而其余三个则都有。先是检查代码没有问题，再检查数据质量都是好的。之后对比作者上传的bigWig文件，ESC-2样本确实无峰，是此可以」排除我们数据处理**的问题。

经过思考，个人猜测可能是人胚胎干细胞（ESC）ChIP-seq样本制备问题，ESC样本制备困难是已知的，因此还有相应的优化procedure发表，感兴趣可查：*An Optimized Protocol for ChIP-Seq from Human Embryonic Stem Cell Cultures (STAR Protocols)*。

「peaks.narrowPeak」文件

narrowPeak文件是我们在Chip-seq公共数据中常看到的数据类型。可以检查narrowPeak文件里peak的数目，可见如我们xls文件缩减，ESC-2样本是一个peak都没call到的。

用IGV直接打开后如下图，我们明确知道CDKN1A基因是P53蛋白的经典下游基因，因此我们使用该基因作为参考，可见ESC-2样本仍旧无peak，则从我个人的角度来看，这个样本或是不合格样本，其数据不可用于统计或者后续分析，此时作者设置重复样本的优势就体现出来了。除了发现不合格样本之外，还可以看到P53蛋白在神经前体细胞组（NPC）的结合位点多了一个，这意味着什么，我们还需要进一步思考。

「BED文件」

该文件一样可以使用IGV直接打开，可看到一个非常小大约一个碱基大小的小块，对应着某个碱基。这文件记录的是峰顶的位置，只有一个碱基的大小是可以理解的。具体的图就不放出来了。

至此，我们peak calling基本结束，可以看到从ESC到NPC分化，我们找到600-1000个peak，如果一个基因平均有1-3个peak，那么最后注释出来的基因会有200-400个左右，通过比较P53蛋白在ESC和NPC中binding基因的交集或者子集，或许能找到对NPC分化关键基因，如此集合转录组或者其他基因组学信息，能够进一步缩小下一步研究的范围。