流程跑了那么多，不看一点软件原理吗（以MACS为例） / 四六文摘

有很多粉丝在婉拒我邀请他写笔记的时候最常用的一个理由是，生信技能树已经有了那么多好的笔记，他自己实在是想不出可以写啥。比如NGS流程系列：

但实际上，这仅仅是借口，我把自己七八年的流程经验分享给大家，就必然会漏掉很多细节，这都是大家可以发挥的地方啊！比如我们的转录组讲师在学习我的《ChIP-seq数据分析》流程就针对其中的MACS软件进行细致的讲解：

下面是转录组讲师的投稿

学习目的

描述MACS2 peak calling 算法的不同组成部分
描述运行MACS2所涉及的参数
列出并描述MACS2的输出文件

Peak Calling流程

Peak calling，用于识别用ChIP测序实验产生的数据比对到基因组中的reads富集的区域。

对于ChIP-seq实验，我们可以从比对文件中观察到，以结合位点为中心，read密度在+/-链上的分布不对称。所选片段的5 '端将在正链和负链上形成group。然后用统计方法评估这些group的分布，并与背景(输入或模拟IP样本)进行比较，以确定富集位点是否可能是一个真正的结合位点。

有各种工具可用于peak calling，而MACS2是最常用的程序之一。我们将在本节中演示它。注意，在这个例子中术语“tag”和序列“read”是可以互换使用的。

注意：

我们使用的数据的目的是研究两个转录因子，因此我们的重点是鉴定作为点状结合位点（ punctate binding sites）的简并短序列（short degenerate sequences）。ChIP-seq分析算法专门用于识别两种富集类型中的一种(或对每种有特定的方法):

宽峰或宽域（broad peaks）：即覆盖整个基因体的组蛋白修饰
窄峰（narrow peaks）：即转录因子结合。窄峰更容易检测，因为我们寻找的是具有更高振幅的区域，与宽或分散的标记（marks）相比，更容易从背景中区分出来。
“混合”（mixed）：还有一些“混合”（mixed）结合谱很难被算法识别。这方面的一个例子是PolII的结合特性，它在启动子结合，跨越了基因的长度，导致混合信号(narrow and broad)。

MACS2

一种常用的识别转录因子结合位点的工具叫做ChIP-seq模型分析(Model-based Analysis of ChIP-seq，MACS)。mac算法捕捉基因组复杂性的影响，以评估腹肌的ChIP区域重要性。虽然它是为检测转录因子结合位点而开发的，但它也适用于更大的区域。

MACS结合了测序tags的位置和方向信息，提高了结合位点的空间分辨率。mac可以很容易地单独用于ChIP样本，或者与增加peak calls特异性的对照样本一起使用。mac的工作流程如下所示。在这一课中，我们将更详细地描述这些步骤：

Removing redundancy

MACS提供了不同的选项来处理完全相同位置的重复tags，即具有相同协调和相同sttand的tags。默认情况下，在每个位置保留一个read。auto选项，这是非常常用的，告诉mac基于二项分布计算在完全相同的位置的最大tags数，使用1e-5作为pvalue阈值。另一种方法是设置all选项，它保留每个tags。如果指定了integer，那么许多tags将保存在相同的位置。这种冗余性适用于ChIP和对照样本。

为什么需要考虑重复片段？

具有相同起始位点的reads被认为是重复片段。这些重复片段可以来自实验，但是能产生ChIP信号。

不良类型的重复：

如果测序核酸分子起始量低，这可能导致在测序前reads过度扩增。PCR中的任何偏差都将加剧这一问题，并可能导致人为富集区域。此外，被列入黑名单(重复)的超高信号区域，也能产生高重复。在分析之前屏蔽这些区域可以帮助解决这个问题。

良好的重复：

ChIP-seq也会存在一定的生物重复reads，因为你只是测了基因组的一小部分。如果你的覆盖深度过深或者你的蛋白质只结合到少数位点，这个数字会增加。如果存在一定的生物学重复reads，去除这些重复可能导致低估ChIP信号。

The take-home

考虑富集效率和测序深度。尝试通过基因组浏览器区分非重复数据。真实的peak会有多个有offsets的重叠reads，而只有PCR重复的样本会在没有offsets的情况下完美叠加。区分生物重复和PCR伪产物的一个可能的解决方案是在你的实验设置中引入UMIs。
为差异binding分析保留重复。
如果你希望在重复的区域结合，使用双端测序并保留重复。

否则，最佳实践是在peak calling之前删除重复reads

Modeling the shift size

真实结合位点周围的tags密度应呈现双峰富集模式(或成对峰)。MACS利用这种双峰模式对移动的尺寸进行经验模拟，以更好地定位精确的结合位点。

为了找到成对的peaks来建立模型，MACS首先扫描整个数据集，寻找高度显著的富集区域。这只是使用ChIP样本完成的！给定声波大小(带宽bandwidth)和高可信的折叠富集(fold-enrichment，mfold)， mac在基因组中滑动两个带宽窗口，找到tags富集超过mfold的区域，相对于随机的标签基因组分布。

MACS随机抽取1000个这样的高质量peaks，分离它们的正链和负链tags，并通过它们中心之间的中点对它们进行对齐。对准过程中两个peaks的模之间的距离定义为“d”，表示估计的fragment长度。MACS将所有的tags向3 '端移动d/2，以找到最可能的蛋白质-DNA相互作用位点。

Scaling libraries

在input和处理样本之间测序深度不同的实验中，MACS线性缩放控制tags总数，使其ChIP tag总数相同。默认的做法是为了更大的样本被缩小。

Effective genome length

从tags count计算λBG，MAC2需要有效基因组大小或用来比对的基因组的大小。比对性能与在基因组特定位置的k-mers唯一性有关。低复杂性和重复区域具有低唯一性，即低的比对性能。因此，我们需要提供有效的基因组长度，以纠正低比对性能区域的真实信号丢失。

那么，如何获取有效参考基因组长度呢？

MACS2软件有一些预先计算好的常用生物(human, mouse, worm and fly)的有效参考基因组长度。你可以基于你的物种和参考基因组计算一个更精确的值。deepTools工具有计算好的比较新的参考基因组有效长度，也可以直接用来计算。

Peak detection

当mac将每个tag按d/2移动后，它会使用2d大小的窗口在基因组中滑动，以找到候选peaks。tags沿基因组的分布可以用泊松分布来模拟。泊松是一个单参数模型，其中参数λ是该窗口的预期read数。

公式看起来略微有一点点复杂哦！主要是背后的一个简单的统计学概念呢，泊松分布！

MACS没有使用从整个基因组估计的统一的λ，而是使用了一个动态参数，λlocal，即每个候选peak。lambda参数是根据对照样本估计的，并通过在各种窗口大小上取最大值来推导的。

λlocal = max(λBG, λ1k, λ5k, λ10k).

通过这种方式，lambda捕获了局部偏倚的影响，并且在小的局部区域低tag计数是稳健的。这些偏差的可能来源包括局部染色质结构、DNA扩增和测序偏倚，以及基因组拷贝数变异。

**p值：**如果基于λ的泊松分布p值< 10e-5(可以从默认值更改)，则该区域具有显著的tags富集。

fold enrichment：将重叠的富集peaks进行合并，每个tag位置从中心向外扩展'd’个碱基。最高fragment堆积峰的位置，以后称为顶点（summit），被预测为精确的结合位置。ChIP-seq tag数和λlocal之间的比值被报告为fold富集（fold enrichment）。

Estimation of false discovery rate

每个peak都被认为是一个独立的检验，因此，当我们遇到在一个样本中检测到的数千个有效peak时，我们就有了多个检验的问题。在MACSv1.4中，FDR是通过交换ChIP和对照样品进行经验测定的。在MACS2中，现在使用Benjamini-Hochberg校正对p值进行多重比较。

总结

从这篇文章中，我们可以看到：

peak是如何检验出来的
在进行peak calling之前是否去重的问题的讨论
后续软件中需要重点关注的几个参数：

-g（有效参考基因组长度的大小，设置不正确会影响λ值的计算）
p值是什么统计模型计算得到
fold enrichment值的理解

参考资料

https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html

流程跑了那么多，不看一点软件原理吗（以MACS为例）