刘小泽写于2020.6.12
来自：https://hbctraining.github.io/Intro-to-ChIPseq/lessons/05_peak_calling_macs.html
看到好资料顺手整理一下
-
MACS全称是：Model-based Analysis of ChIP-seq (MACS)：https://github.com/taoliu/MACS
-
它开发的初衷是检测转录因子结合位点（transcription factor binding sites，TFBS），后来也出了broad peaks和mix模式（比如PolII的结合存在于启动子和跨越基因长度，broad和narrow并存）
-
通过整合测序reads位置和来源信息，增强了结合位点的空间分辨率
-
可以针对某一个样本，也可以结合一个control对照（可以增加鉴定peaks的特异性）
-
整个流程是：
去重复
-
重复的定义是：同一链的同一坐标，多个reads
-
提供了多个参数：默认（auto）是一个位置只保留一条read。另外还有：
- all：保留所有重复reads
- integer：指定一个位点最多包含几个重复reads
-
这个操作对ChIP和control样本同时处理
-
重复也分好坏
-
不容忍的重复：起始底物浓度低，容易导致过度扩增（overamplification），产生的是PCR重复；另外blacklisted (repeat) regions with ultra high signal也具有很高的重复，因此在分析前去掉这部分区域也是有帮助的
-
可以考虑的重复：如果测序深度过高或者蛋白只是结合很少几个位点，那么这其实是真的生物学重复，去掉会丢失一些信号
-
如何处理？
- Consider your enrichment efficiency and sequencing depth. 用IGV查看：真正的峰会有多个有偏移的重叠reads，而只有PCR重复的样品在没有偏移的情况下也能完美地叠加
- 如果要进行differential binding analysis，可以保留
- 如果知道是在基因组重复区域结合的，就需要使用双端测序，并要保留
- 除此以外，最好在call peaks之前去掉重复
-
Modeling the shift size
一个真实的结合位点处的reads分布应该是：双峰（bimodal or paired peaks）。MACS就会利用这个双峰进行建模，估计shift size，来定位预测的结合位点
如何建模？
- 基于ChIP样本，在全基因组范围扫描数据，寻找置信度高的富集区域
- 需要指定参数bandwidth（超声破碎的大小）、mfold（阈值high-confidence fold-enrichment）
- 然后沿着基因组滑动窗口
shift size：
MACS会随机抽样1000个高质量的peaks，将正负链的reads区分开来，并按照reads的中心点进行映射，得到两种peaks，峰越高说明这里的reads越多。最后，两种peaks的中心点之间的距离定义为d（表示estimated fragment length）。之后MACS将所有reads向3‘端移动d/2的距离，得到最可能的蛋白-DNA结合区域
Scaling libraries
如果input和treat样本之间的测序深度不同，MACS会对文库进行矫正，让input向treat靠拢。结果就是：total control read count to be the same as the total ChIP read count
中间计算过程还需要考虑
- Effective genome length
-
Peak detection：
The read distribution along the genome can be modeled by a Poisson distribution 泊松分布
The Poisson is a one parameter model, where the parameter λ is the expected number of reads in that window.
MACS uses a dynamic parameter, λlocal, defined for each candidate peak.
Estimation of false discovery rate
A region is considered to have a significant tag enrichment if the p-value < 10e-5
MACS 1.4版本设定FDR是根据经验
MACS2：多重比较 Benjamini-Hochberg correction
MACS2的功能
这里重点关注核心功能——callpeak
输入参数
- -t：IP样本（表示treat）
- -c：control或mock样本
- -f：数据类型，默认是AUTO，MACS自动判断数据类型
可以指定：ELAND, BED, ELANDMULTI, ELANDEXPORT, SAM, BAM, BOWTIE, BAMPE, or BEDPE【Eland：是过时了的带比对信息的序列文件，目前基本已经被淘汰】 - -g：有效基因组大小。由于染色体存在重复区域，实际比对过程中使用的基因组，比实际的基因组范围要小（可能也就是真实基因组的90%或70%）
默认是hs -- 2.7e9，另外还有：mm: 1.87e9、ce: 9e7、dm: 1.2e8
输出参数：
- --outdir：输出文件夹
- -n：输出文件的前缀
- -B/--bdg：把fragment pileup, control lambda, -log10pvalue and -log10qvalue这些值存到bedGraph文件，命名为NAME_treat_pileup.bdg或NAME_control_lambda.bdg
Shifting model参数：
- -s/--tsize：reads的长度。如果不设置，MACS会使用前10条序列自己判断
- --bw：建模过程中使用的bandwidth值，即破碎后DNA fragment片段大小
- --mfold：建模过程使用的阈值
Peak calling参数：
- -q：minimum FDR cutoff
- -p：p-value cutoff（替代q-value）
- --broad：调用broad peak calling
举个例子：
macs2 callpeak -t bowtie2/H1hesc_Nanog_Rep1_aln.bam \
-c bowtie2/H1hesc_Input_Rep1_aln.bam \
-f BAM -g 1.3e+8 \
-n Nanog-rep1 \
--outdir macs2
整个运行过程会打印很多日志，所以可以使用重定向将日志保存
输出的文件
一般一个样本会得到6个文件：
-
_peaks.narrowPeak: BED6+4 format file which contains the peak locations together with peak summit, pvalue and qvalue -
_peaks.xls: a tabular file which contains information about called peaks. Additional information includes pileup and fold enrichment -
_summits.bed: peak summits locations for every peak. To find the motifs at the binding sites, this file is recommended -
_model.R: an R script which you can use to produce a PDF image about the model based on your data and cross-correlation plot -
_control_lambda.bdg: bedGraph format for input sample -
_treat_pileup.bdg: bedGraph format for treatment sample
了解下文件格式：
- narrowPeak：ENCODE计划使用的BED 6+4格式，前6列是标准格式，后4列是附加
-
WIG：全称是Wiggle format (WIG)，表示基因组上一个区域的信号。详情在http://genome.ucsc.edu/goldenPath/help/wiggle.html
它的二进制版本是：bigwig，可以用UCSC的工具wigToBigWig 转换一下例如：
# 来自：http://www.chenlianfu.com/?p=2366 track type=wiggle_0 name="sampleA1" description="RNA-Seq read counts of species A" variableStep chrom=chr01 span=10 10001 13 10011 15 10021 12 fixedStep chrom=chr01 start=100031 step=10 span=10 17 15 20 1. 第一行必须如理示例中格式。只有name和description这两个参数的值可以随意填写。 2. 有两种方法进行数据描述。分别是variableStep和fixedStep。前者数据内容用2列表示，后者数据部分仅用1列表示。 3. 这两种方法的几个参素意义为： chrom 设置序列名 start fixStep中Locus的起始位置 step fixStep中Locus的步进 span 一个数据对应碱基数目 -
BedGraph：用窗口的方式代替原始的每个碱基的测序深度，文件分为track line、data line，详情在：http://genome.ucsc.edu/goldenPath/help/bedgraph.html
bedgraph在原始depth的基础上合并了相同测序深度的连续碱基
统计一下每个样本中有多少peaks：
wc -l *.narrowPeak
如果要看shift size的图：
同时还会生成一个cross-correlation plot（Pearson相关性图）
Rscript Nanog-rep1_model.r
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论