MACS的原理

泊松分布

泊松分布是统计与概率中重要的离散分布之一，泊松分布表示在一定的时间或空间内出现的事件个数，比如某一服务设施在一定时间内受到的服务请求的次数、DNA序列的变异数、汽车站台的等候人数。根据MACS的论文中所描述的，Chip-Seq实验中全基因组的reads分布恰好是符合泊松分布的。

泊松分布的概率分布为
$P(X = k)=\frac{e^{-\lambda}\lambda^{k}}{k!}$
其中e代表的是自然常数，而是单位时间（或单位面积）内随机事件的平均发生率，比如在一定时间内某一服务设施受到的请求次数是5次。

Poisson.png

另外，泊松分布实际上只有一个参数，即 $\lambda$ ，其方差和期望也是 $\lambda$ 。同时，随着 $\lambda$ 的增加，图像分布会趋于对称。

参考资料

二项分布、泊松分布、正态分布的关系
 泊松分布和指数分布：10分钟教程
 wiki_泊松分布

MACS的算法概览

Adjusting read position based on fragment size distribution

Chip-Seq的主要过程为：交联——超声破碎——特异性识别——测序。所以我们测序得到的片段就是我们转录因子结合位点周围的片段。需要注意的一点是，MACS软件出现的年代是2008年，那时候的测序读长都很短，大约50bp左右，且以单端测序为主，并没有真实反应DNA-蛋白结合片段的长度。所以说，我们如果拿测得的50bp去做reads数目的堆积，势必会与真实的结合位置有一定的偏移。事实上，测序的短reads会在真实的结合位置两侧形成双峰，如下图所A示。这也是MACS双峰模型构建的理论基础。

bimodal.png

值得一提的是，像转录因子一类的蛋白与DNA，其结合位点比较narrow，所以双峰模型的构建是比较合理的。但像图B所示的，一些蛋白与DNA会产生较宽的结合区域（诸如一些组蛋白修饰），这时候双峰就不那么显著了。

更为麻烦的是，有时候会有一些混合的结合位点模式，比如Polll蛋白，其会在启动子区域结合，也会覆盖整个基因区域。

为了衡量真实的测序片段大小，d，MACS会粗略地以2倍的超声破碎片段长度作为window来鉴定初步的富集区域。为了避免重复区域或者PCR导致极端富集区域的影响，MACS会随机挑选1000个区域作为模型peak构建区域。这些区域的reads富集程度是基因组背景的10-30倍。对于每个区域的模型peak，MACS都会分离出对比到正链和负链上的reads，然后分别计算出这些reads的位置。从而分别构建出这个区域内的正负链上的模型peak，正负链上模型peak顶点之间距离就记为d。在d确定之后，所有的reads都会朝着3'的方向横移（shift）d/2的距离，从而更好地模拟出蛋白-DNA结合位点。

在2012年的Identifying ChIP-seq enrichment using MACS这篇文章中，作者也提到对于一些过度破碎或者有着很宽的结合位点情况，可能会造成算出来的d很小。对于这种情况，我们一般建议用一个特定的片段长度，而非是预测出来的d。

注意shift和extend的区别，在2008年原始的MACS文章中，作者用的是shift，而到了12年的文章，作者写错，写成了extend。当然，在MACS2中，这两种情况都存在了。

Calculate peak enrichment using local background normalization

基于先前已经调整位置的reads，MACS会在全基因组范围内以2d长度的window来寻找那些有显著富集的区域。有重叠的window会融合成一个候选区域。因为会有许多因素影响不同范围内的reads富集程度，所以MACS用了动态的参数 $\lambda_{local}$ 来对于reads数目的富集进行泊松分布的建模。即MACS并不会用一个常数 $\lambda$ ，而是用一个会在不同区域有变化的 $\lambda_{local}$ 。动态参数值定义为
$\lambda_{local}=max(\lambda_{BG},[\lambda_{region},\lambda_{1k}],\lambda_{5k},\lambda_{10k})$

$\lambda_{BG}$ 来自于全基因组的计算， $\lambda_{region}$ 则来自在control中的对应区域，剩下的 $\lambda_x$ 则来自control中，以得到的候选区域为中心，1k，5k，10k范围内的区域计算。见下图

lambda.png

如果control不在，则local值只是在Chip的样本中计算，而region和1k值也会被舍弃。同时如果Chip-Seq和control的样本测序深度不同，MACS会默认地把测序深度更深的样本缩放。

关于以及p值这一步的计算可能需要看源代码才可以了解了。

但根据MACS2的wiki来说，似乎p值和的计算都是以单个碱基为单位考虑的。

基于泊松分布的模型，我们就可以以单尾检验，计算出p值了。MACS默认以p=1 x 10-5为阈值。

Estimating the empirical false discovery rate by exchanging ChIP-seq and control samples

这里MACS用的Chip和control的置换，从而检验出FDR值我并没有看懂。不过MACS2用的已经是Benjamini-Hochberg方法了，还是比较好懂的。

参考资料：

Evaluation of Algorithm Performance in ChIP-Seq Peak Detection
Model-based Analysis of ChIP-Seq (MACS)
Identifying ChIP-seq enrichment using MACS
In-depth-NGS-Data-Analysis-Course

MACS2中的一些参数介绍

-f/--format FORMAT

可以接受多种格式参数，默认使用AUTO来检测格式。但并不能检测“BAMPE”或者“BEDPE”格式，即双端测序格式。所以，当你的数据是双端测序数据时，你应该用BAMPE或者BEDPE参数。当你设置成双端参数的时候，MACS2就会跳过建模计算d的那一步，而是直接用片段的insert size来建立堆积。

--extsize

如果使用这个参数，那么MACS就会使用你设置的数值，来把reads从5‘—3’补齐到你指定的数值。这个参数只有当--nomodel参数设置了，或者MACS建模失败，--fix-bimodal开启的时候才可以用。

--shift

shift参数会先于extsize参数执行。如果你设置的数值为正，reads会从5‘—3’偏移，而数值为负，reads会从3'—5‘偏移。当格式为BAMPE或者BEDPE的时候，不能设置参数。

--broad

会放宽cutoff的阈值，然后把临近的区域结合起来，形成较宽的peak区域。与broad-cutoff参数是一起的，broad-cutoff参数默认为q-value的参数，为0.1。

有趣的是，shift后面数值如果为正，则正负链的reads会朝着中心偏移，如果后面数值为负，则正负链的reads会各自远离中心，即正链reads向左，负链reads向右。

给个例子：

chr1    500    550    read1    .    (+)
chr1    700    750    read2    .    (-)

--shift -100
chr1    400    450    read1    .    (+)
chr1    800    850    read2    .    (-)

--extsize 200
chr1    400    600    read1    .    (+)
chr1    650    850    read2    .    (-)