随着测序技术的改进，染色质免疫沉淀后的高通量测序（ChIP-Seq）在研究全基因组蛋白质-DNA相互作用方面越来越流行。为了解决缺乏强大的ChIP-Seq分析方法的问题，我们提出了基于模型的ChIP-Seq分析（MACS），用于识别转录因子结合点。MACS抓住了基因组复杂性的影响来评估富集的ChIP区域的意义，MACS通过结合测序标签位置和方向的信息来提高结合位点的空间分辨率。MACS可以很容易地单独用于ChIP-Seq数据，或与对照样本一起使用，并提高了特异性。此外，作为一个通用的峰值调用器，MACS也可以应用于任何 "DNA富集检测"，如果要问的问题只是：我们在哪里可以找到比随机背景更重要的read覆盖。
macs3是macs2的更新版本
什么是Call Peak
Call Peak其实是一种计算方法，用于识别基因组中由测序得到的比对reads富集的区域。通俗的来说Call peak就是把我们有蛋白富集到核酸时的区域给找出来
比对结果的文件中reads在正负链不均匀分布，但在结合位点聚集。正负链5‘末端的reads各形成一组合，通过统计学的方法评估这些组合的分布并和对照组比较，确定这些结合位点是否是显著的。
为了找个蛋白结合核酸的区域，我们需要进行计算，因为其实在通过测序我们得到的只不过是蛋白结合核酸片段的末端，而非蛋白真正的结合区域。所以就需要在得到的reads进行延伸，
image.png
image.png
理想情况下，我们认为测得的reads最终其实会近似地平均分配到正负链上，这样，对于一个结合热点而言，reads在附近正负链上会近似地形成“双峰”，我们偏移后就可以得到真正的结合区域了。
文档
MACS/callpeak.md at master · macs3-project/MACS (github.com)
image.png
去重：首先软件会对数据进行去重处理；
建模：reads在附近正负链上会近似地形成“双峰”，这个双峰之间的距离称之为d，软件需要进行延伸需要知道这个d值。所以MACS选取了1000个表达10~30倍的regions去建模以计算这个d值。值得注意的是，对组蛋白修饰的ChIP-seq数据peak-calling时，“双峰模型”会建立失败，这是因为组蛋白修饰往往并不是孤立存在的，可能很长一段染色质区间都被同一个组蛋白修饰占据，换句话说，组蛋白修饰的peak并不典型。这也就是下面有分开两者模式去计算的原因。
延伸：得到了d之后，向3'端延伸就可以去进行延伸，延伸长度为d/2；
归一化：把两组数据进行归一化，各组数据可比；
获得候选的peak和背景进行比较；
计算确定λ，代入泊松分布，得到p值；
注意你要找的是转录因子还是组蛋白，一个生成的是narrowpeak，一个生成的是broadPeak，实测后者比前者多上起码百分之二十五以上的peak。
本质上，broad ，计算的是如组蛋白修饰在整个基因body区域的分布；narrow peak，计算如转录因子的结合。narrow peak相对于broad 或者分散的marks更易被检测到。也有一些混合的结合图谱，如PolII包括narrow和broad信号。
image.png
Comparative analysis of commonly used peak calling programs for ChIP-Seq analysis - PMC (nih.gov)
image.png
Callpeak
寻找ChIP-seq的转录因子（TF）的命令：
### 生成的文件是narrowPeak
macs3 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
寻找ChIP-seq的组蛋白（Histone ）Mark的命令：
### 生成的文件是broadPeak
macs3 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
对于ATAC-seq双端数据的操作：
macs3 callpeak -f BAMPE -t ATAC.bam -g hs -n test -B -q 0.01
-f是输入文件的格式，可以是BAM、BED，软件可以自动识别，但是需要注意自动识别是无法区分是PE还是SE，所以建议手动指定；
-g是有效基因组大小，软件预设了模式物种的有效因组大小，如果hs人，mm小鼠；
-n是样本名字；
-B是以bedGraph格式存放fragment pileup, control lambda, -log10pvalue 和log10qvale,使用会增长耗时；
-q就是用q值（FDR）进行筛选，输入得到值为阈值；
broad-cutoff用于过滤 broad得到的peak。
本人使用
注意这里是双端的，如果你要单端请用BAM
macs3 callpeak -t me.rmdup.uniq.bam -c Incontrol.uniq.bam -f BAMPE -n test2_me --broad -g hs --broad-cutoff 0.1
参数
- -t/--treatment FILENAME
这是 MACS 的唯一必需参数。该文件可以是任何受支持的格式——请参阅 --format 选项中的详细信息。如果您有多个对齐文件，您可以将它们指定为 -t A B C。MACS 会将所有这些文件汇集在一起。
-
-c/--control
对照、基因组输入或模拟 IP 数据文件。请遵循与 -t/--treatment 相同的方向。 -
-n/--name
实验的名称字符串。 MACS 将使用此字符串 NAME 创建输出文件，如 NAME_peaks.xls、NAME_negative_peaks.xls、NAME_peaks.bed、NAME_summits.bed、NAME_model.r 等。因此，请避免这些文件名与您现有文件之间的任何冲突。 -
--outdir
MACS3 会将所有输出文件保存到此选项的指定文件夹中。如有必要，将创建一个新文件夹。 -
-f/--format 格式
标签文件的格式可以是 ELAND、BED、ELANDMULTI、ELANDEXPORT、SAM、BAM、BOWTIE、BAMPE 或 BEDPE。默认为 AUTO，这将允许 MACS 自动决定格式。当您组合不同格式的文件时，AUTO 也很有用。
请注意，MACS 无法使用 AUTO 检测 BAMPE 或 BEDPE 格式，您必须指定 BAMPE 和 BEDPE 的格式。（这话就是说如果是双端bam，就要用BAMPE，而不是用BAM或者AUTO，否则你看着结果就感觉日了狗）
如今，最常见的格式是 BED 或 BAM（包括 BEDPE 和 BAMPE）。我们的建议是首先将您的数据转换为 BED 或 BAM。
此外，MACS3 可以检测和读取gzipped 文件。比如.bed.gz文件可以直接使用，不用--format BED解压。
以下是推荐格式的详细说明：
BED
BED 格式可在 UCSC 上找到。
BED 格式输入中的基本列是第 1 列染色体名称、第 2 列起始位置、第 3 结束位置和第 6 条链。
请注意，对于 BED 格式，MACS 需要第 6 列链信息（也就是正负链）。请注意，BED 格式的坐标是从零开始的，并且是半开的。在 UCSC 网站上查看更多详细信息。
image.png
BAM/SAM
如果格式是 BAM/SAM，请查看 (http://samtools.sourceforge.net/samtools.shtml) 中的定义。如果 BAM 文件是为双端数据生成的，MACS 将只保留 left mate(5' end) tag。但是，当指定格式 BAMPE 时，MACS 将使用从比对结果中推断出的真实片段进行读取堆积。
BEDPE or BAMPE
当格式指定为 BAMPE 或 BEDPE 时，将触发特殊模式。这样，MACS3 会将 BAM 或 BED 文件作为双端数据处理。 MACS3 不会构建正负链读取的双峰分布来预测片段大小，而是使用读取对的实际插入大小来构建片段堆积。
BAMPE 格式只是包含双端对齐信息的 BAM 格式，例如来自 BWA 或 BOWTIE 的对齐信息。
BEDPE格式是一种简化的、更灵活的BED格式，它只包含定义染色体名称的前三列，以及来自成对端测序的片段的左和右位置。
请注意，这与bedtools使用的格式不一样，bedtools版本的BEDPE实际上不是标准的BED格式。你可以使用MACS3的子命令randsample将包含双端信息的BAM文件转换成BEDPE格式文件
macs3 randsample -i the_BAMPE_file.bam -f BAMPE -p 100 -o the_BEDPE_file.bed
-
-g/--gsize
这是可映射的基因组大小或有效的基因组大小，它被定义为可测序的基因组大小。由于染色体上的重复特征，实际的可映射基因组大小将小于原始大小，约为基因组大小的90%或70%。对于人类基因组，推荐使用默认的hs-2.7e9。下面是所有预编译的有效基因组大小的参数 -
hs: 2.7e9
-
mm: 1.87e9
-
ce: 9e7
-
dm: 1.2e8
用户可能希望使用 k-mer 工具来模拟 Xbps 长读取到目标基因组的映射，并找到理想的有效基因组大小。然而，通常从整个基因组中去掉简单的重复序列和 Ns，就可以得到一个近似的有效基因组大小。数字上的细微差别不会导致峰值调用的巨大差异，因为这个数字用于估计全基因组噪声水平，与 MACS 建模的局部偏差相比，这个噪声水平通常是最不显著的。
- -s/--tsize
序列标记的大小。如果您没有指定它，MACS 将尝试使用输入治疗文件中的前10个序列来确定标记大小。指定它将重写自动确定的标记大小。
-
-q/--qvalue
调用有效区域的 q 值(最小 FDR)截止值。默认值是0.05。对于广泛的标志，您可以尝试0.05作为截止值。Q 值是使用 Benjamini-Hochberg 程序从 p 值计算出来的。 -
-p/--pvalue
如果指定了 p，MACS3将使用 p-value 而不是 q-value。 -
--min-length, --max-gap
可以使用这两个选项通过指定被调用的峰值的最小长度和允许合并两个邻近区域之间的间隔的最大值来微调峰值调用行为。换句话说，一个被称为峰必须长于最小长度，如果两个附近的峰之间的距离小于最大间隙，那么他们将合并为一个。如果没有设置，MACS3将最小长度的 DEFAULT 值设置为预测的片段大小 d，最大间隙的 DEFAULT 值设置为检测到的读取长度。注意，如果设置的最小长度值小于片段大小，则可能对结果没有影响。对于宽峰调用——宽选项集，合并附近较强峰的 DEFAULT 最大间隙与窄峰调用相同，合并附近较弱(宽)峰的 DEFAULT 最大间隙为4倍。您还可以使用带默认设置的—— cutoff-analysis 选项，并检查不同截止值下的列 avelPeak，以确定合理的最小长度值。
-
-- nolambda
打开这个设置后，MACS 将使用背景 lambda 作为局部 lambda。这意味着 MACS 不会考虑峰值候选区域的局部偏差。 -
--slocal, --llocal
这两个参数控制哪两个水平的区域将检查周围的峰值区域，以计算出最大的局部 lambda。默认情况下，MACS 对于小的局部区域(—— slocal)考虑1000bp，对于大的局部区域(—— llocal)考虑10000bps，其捕获像开放染色质结构域这样的长程效应的偏差。您可以根据您的项目调整这些。记住，如果区域设置得太小，输入数据中的一个尖峰可能会杀死一个重要的峰值。 -
--nomodel
开启时，MACS 将绕过构建shifting model -
--shift
请注意，这不是由 --extsize 取代的传统 --shiftsize 选项！您可以在此处设置 bp 的任意偏移。请谨慎设置默认值 (0) 以外的值。设置 --nomodel 时，MACS 将使用此值移动切割端 (5')，然后应用 --extsize 从 5' 到 3' 方向将它们扩展到片段。当该值为负时，两端将向 3'->5' 方向移动，否则向 5'->3' 方向移动。建议将其保留为 ChIP-Seq 数据集的默认值 0，或 -1 * EXTSIZE 的一半以及 --extsize 选项以检测诸如某些 DNAseI-Seq 数据集之类的丰富切割位点。请注意，如果格式为 BAMPE 或 BEDPE 用于双端数据，则不能设置 0 以外的值。默认值为 0。
以下是一些组合 --shift 和 --extsize 的示例：
寻找丰富的切割位点，例如一些 DNAse-Seq 数据集。在这种情况下，测序读数的所有 5' 端都应在两个方向上延伸，以平滑堆积信号。如果想要的平滑窗口为 200bps，则使用
--nomodel --shift -100 --extsize 200。
对于某些核小体序列数据，我们需要使用核小体一半大小来积聚（pile up）核小体的中心以进行小波分析（例如 NPS 算法）。
由于包裹在核小体上的 DNA 约为 147bps，因此可以使用此选项：
--nomodel --shift 37 --extsize 73
-
--keep-dup
它控制 MACS 对完全相同位置的重复标签的行为——相同的协调和相同的链。默认的自动选项使 MACS 使用 1e-5 作为 p 值截止值，根据二项分布计算完全相同位置的最大标签；并且 all 选项保留每个标签。如果给定一个整数，则最多将这个数量的标签保存在同一位置。默认设置是在同一位置保留一个标签。默认值：1 -
--broad
当此标志打开时，MACS 将尝试通过将附近高度富集的区域放入具有松散截止的广泛区域来合成 BED12 中的广泛区域（一种类似基因模型的格式）。广域由通过--broad-cutoff 参数的截止值控制。请注意，合并附近较弱/较宽峰的 max-gap 值是合并附近较强峰的 max-gap 的 4 倍。后者可以通过 --max-gap 选项控制，默认是 PE 数据中的平均片段/插入长度。默认值：False -
--broad-cutoff
broad region的截止值。除非设置了 --broad，否则此选项不可用。如果设置了 -p，这是一个 p 值截止，否则，它是一个 q 值截止。默认值：0.1 -
--scale-to <large|small>
当设置为large时，将较小的数据集线性缩放到与较大数据集相同的深度。默认情况下或设置为small，较大的数据集将向较小的数据集缩放。请注意，放大小数据会导致更多误报。 -
-B/--bdg
如果此参数打开，MACS 将存储fragment pileup，控制 lambda 在 bedGraph 文件中。两个 bedGraph 文件 ，实验组数据将存储在当前目录中名为 NAME_treat_pileup.bdg，NAME_control_lambda.bdg 用于存储来自控制组的local lambda 值。 -
--call-summits
MACS 现在将重新分析 signal profile的形状（p 或 q-score，取决于截止设置），从而对从通用程序call出来的的每个peak内的亚peak进行去卷积。强烈建议检测相邻的绑定事件。使用时，大peak区域的输出亚peak将具有相同的峰边界，以及不同的分数和峰顶位置。 -
--buffer-size
MACS 使用缓冲区大小来逐步增加内部数组大小，以存储每个染色体或重叠群的读取对齐信息。为了增加缓冲区大小，MACS 可以运行得更快，但如果某些染色体/重叠群只有很少的读取，则会浪费更多的内存。在大多数情况下，默认值 100000 可以正常工作。但是，如果您的比对中有大量染色体/重叠群并且每个染色体/重叠群的读数很少，建议指定较小的缓冲区大小以减少内存使用（但读取比对文件需要更长的时间） .读取对齐文件所需的最小内存约为 # of CHROMOSOME * BUFFER_SIZE * 8 Bytes。默认值：100000
输出文件
NAME_peaks.xls
是一个表格文件，其中包含有关被调用峰的信息。您可以在 excel 中打开它并使用 excel 函数进行排序/过滤。信息包括：
- 染色体名称
- peak起始位置
- peak的结束位置
- peak长度
- 峰顶的位置(absolute peak summit position)
- 峰顶上的堆积高度(pileup height at peak summit)
- -log10(pvalue) for the peak summit （例如 pvalue =1e-10，那么这个值应该是 10）
- 该峰的倍数富集度与当地λ的随机泊松分布相对应
- -log10(qvalue) at peak summit
xls文件 中的坐标是基于 1 的，这与 BED 格式不同。当为peak calling是--broad 时，XLS 文件中的pileup、p value、q value和fold change将是整个峰值区域的平均值，因为在broad peak calling mode 下，不会call peak summit。
NAME_peaks.narrowPeak
是 BED6+4 格式文件，其中包含peak summit position , p-value, and q-value。可以将其加载到 UCSC 基因组浏览器。一些特定列的定义是：
- 5列：整数值。根据是否使用 -p (pvalue) 或 -q (qvalue) 作为分数截止值，它被计算为 int(-10*log10pvalue) 或 int(-10*log10qvalue)
注意，当前此值可能超出 UCSC ENCODE narrowPeak
格式中定义的 [0-1000] 范围。您可以使用以下 1行 awk 让这列值中大于1000的全部改为1000（即 p/q-value = 10^-100）：
awk -v OFS="\t" '{$5=$5>1000?1000:$5 } {print}' NAME_peaks.narrowPeak
- 6列：strand - +/- 表示链或方向（如果适用）。如果没有指定方向,为 ”.”。
- 7列：fold-change at peak summit
- 8列: -log10pvalue at peak Summit
- 9列：-log10qvalue at peak summit
- 10列:relative summit position to peak start(峰顶位置与峰顶起点的相对位置)
NAME_summits.bed
BED 格式，其中包含每个峰的峰顶位置。此文件中的第 5 列与narrowPeak文件中的相同。如果您想在binding sites找到 motifs，建议使用此文件。该文件可以直接加载到 UCSC 基因组浏览器。如果您想通过其他工具对其进行分析，请删除起始轨迹线(beginning track line)。
NAME_peaks.broadPeak
BED6+3 格式，类似于narrowPeak 文件，除了缺少注释峰顶的第 10 列。该文件和 gappedPeak 文件仅在启用 --broad 时可用。由于在宽峰调用模式下，不会调call peak summit，因此第 5 列和第 7-9 列的值是峰区所有位置的平均值。如果要修复第 5 列中的值问题，请参阅narrowPeak。
chr1 1871697 1871949 test2_me1_peak_1 13 . 3.60763 5.19273 1.30816
chr1 1885202 1886229 test2_me1_peak_2 13 . 3.78524 5.22165 1.32138
chr1 2274443 2275382 test2_me1_peak_3 18 . 3.91899 6.14386 1.81172
chr1 2292167 2292821 test2_me1_peak_4 12 . 3.77953 5.14002 1.24934
chr1 2299127 2299497 test2_me1_peak_5 11 . 3.7478 5.02819 1.17133
NAME_peaks.gappedPeak
采用 BED12+3 格式，里面包含宽峰和窄峰。
- 第 5 列是如narrowPeak 中在 UCSC 浏览器上显示灰度的分数。
- 第 7 列是该区域第一个narrowPeak的起点
- 第 8 列是终点
- 第 9 列应该是 RGB 颜色，但是，我们在此处保留 0 以使用默认颜色，因此如果需要，请更改它
- 第 10 列表示有多少块，包括广泛区域的起始 1bp 和结束 1bp
- 第 11 列显示每个块的长度，第 12 列显示每个块的开始
- 第 13 位：倍数变化
- 第 14 位：-log10pvalue
- 第 15 位：-log10qvalue。
该文件可以直接加载到 UCSC 基因组浏览器。如果要修复第 5 列中的值问题，请参阅narrowPeak。
NAME_treat_pileup.bdg 和 NAME_control_lambda.bdg
这两个文件是 bedGraph 格式，可以导入 UCSC 基因组浏览器或转换成更小的 bigWig 文件。
-
NAME_treat_pielup.bdg 包含来自 ChIP/处理样本的pileup signals（根据 --scale-to 选项标准化）。
-
NAME_control_lambda.bdg 包含从对照样本中估计的每个基因组位置的局部偏差，或者当对照样本不存在时从处理样本中估计的局部偏差。
-
子命令 bdgcmp 可用于比较这两个文件，并制作一个包含 p 值、q 值、对数似然和对数倍数变化等分数的 bedGraph 文件。
表观组学：call peak概念和MACS2算法原理，最简明清晰的call peak一文通 (qq.com)
macs3-project/MACS: MACS -- Model-based Analysis of ChIP-Seq (github.com)
Chip-seq处理流程 | Public Library of Bioinformatics (plob.org)
网友评论