了解一下call peaks的MACS

刘小泽写于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也具有很高的重复，因此在分析前去掉这部分区域也是有帮助的

  • 可以考虑的重复：如果测序深度过高或者蛋白只是结合很少几个位点，那么这其实是真的生物学重复，去掉会丢失一些信号

    

  • 如何处理？

  1. Consider your enrichment efficiency and sequencing depth. 用IGV查看：真正的峰会有多个有偏移的重叠reads，而只有PCR重复的样品在没有偏移的情况下也能完美地叠加
  2. 如果要进行differential binding analysis，可以保留
  3. 如果知道是在基因组重复区域结合的，就需要使用双端测序，并要保留
  4. 除此以外，最好在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!