生信地基系列--MACS3 Call peaks

2022-11-06  本文已影响0人  可能性之兽

随着测序技术的改进，染色质免疫沉淀后的高通量测序（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

参数

请注意，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

用户可能希望使用 k-mer 工具来模拟 Xbps 长读取到目标基因组的映射，并找到理想的有效基因组大小。然而，通常从整个基因组中去掉简单的重复序列和 Ns，就可以得到一个近似的有效基因组大小。数字上的细微差别不会导致峰值调用的巨大差异，因为这个数字用于估计全基因组噪声水平，与 MACS 建模的局部偏差相比，这个噪声水平通常是最不显著的。

序列标记的大小。如果您没有指定它，MACS 将尝试使用输入治疗文件中的前10个序列来确定标记大小。指定它将重写自动确定的标记大小。

可以使用这两个选项通过指定被调用的峰值的最小长度和允许合并两个邻近区域之间的间隔的最大值来微调峰值调用行为。换句话说，一个被称为峰必须长于最小长度，如果两个附近的峰之间的距离小于最大间隙，那么他们将合并为一个。如果没有设置，MACS3将最小长度的 DEFAULT 值设置为预测的片段大小 d，最大间隙的 DEFAULT 值设置为检测到的读取长度。注意，如果设置的最小长度值小于片段大小，则可能对结果没有影响。对于宽峰调用——宽选项集，合并附近较强峰的 DEFAULT 最大间隙与窄峰调用相同，合并附近较弱(宽)峰的 DEFAULT 最大间隙为4倍。您还可以使用带默认设置的—— cutoff-analysis 选项，并检查不同截止值下的列 avelPeak，以确定合理的最小长度值。

以下是一些组合 --shift 和 --extsize 的示例：
寻找丰富的切割位点，例如一些 DNAse-Seq 数据集。在这种情况下，测序读数的所有 5' 端都应在两个方向上延伸，以平滑堆积信号。如果想要的平滑窗口为 200bps，则使用

--nomodel --shift -100 --extsize 200。

对于某些核小体序列数据，我们需要使用核小体一半大小来积聚（pile up）核小体的中心以进行小波分析（例如 NPS 算法）。
由于包裹在核小体上的 DNA 约为 147bps，因此可以使用此选项：

--nomodel --shift 37 --extsize 73

输出文件

NAME_peaks.xls

是一个表格文件，其中包含有关被调用峰的信息。您可以在 excel 中打开它并使用 excel 函数进行排序/过滤。信息包括：

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 基因组浏览器。一些特定列的定义是：

注意，当前此值可能超出 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
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 格式，里面包含宽峰和窄峰。

该文件可以直接加载到 UCSC 基因组浏览器。如果要修复第 5 列中的值问题，请参阅narrowPeak。

NAME_treat_pileup.bdg 和 NAME_control_lambda.bdg

这两个文件是 bedGraph 格式，可以导入 UCSC 基因组浏览器或转换成更小的 bigWig 文件。

表观组学：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)

上一篇 下一篇

猜你喜欢

热点阅读