【推荐】如何使用MACS获得ChIP-Seq富集的区域(峰区域)

【声明】本文未经许可,禁止转载!
 一 简介 
随着测序技术的不断改进,染色质免疫共沉淀的高通量测序(ChIP-Seq)已广泛应用于研究基因组蛋白和DNA之间的相互作用关系。基于缺乏强大的ChIP-seq分析方法,提出了一种新的算法——基于模型的ChIP-seq算法分析(MACS),来识别转录因子的结合位点。MACS考虑到了基因组复杂性对评估显著富集的ChIP区域的影响,MACS通过结合tag测序的位置和方向信息,提高了结合位点的空间分辨率。MACS适用于单个的ChIP-seq数据,也可以通过设对照样本来提高特异性和精度。
 
二 使用方法及参数 
 
1. 用法:
下载:http://liulab.dfci.harvard.edu/MACS/Download.html
macs14  -t treat.bed  -c control.bed --tsize tag_len --bw bw --gsize gsize --wig --wigextend 300  --space 50 --diag Ture --mfold 10,30 --name outfile_name 2>outfile_name.log
2. 参数
-t         ChIP-seq treatment files. (输入样本文件,如*.bed)
-c Control files. (对照文件,和输入样本文件格式相同)
--name Used to generate output file names.(输出文件名称)
--format Format of tag file, "ELAND" or "BED" or "ELANDMULTI" or
"ELANDMULTIPET". DEFAULT: "BED"(选择输入文件的格式)--gsize Effective genome size, default:2.7e+9(基因组大小,默认为人)
--tsize Tag size. DEFAULT: 25(基因reads的长度)
--bw Band width. This value is used while building the shifting model.
DEFAULT: 300(扫描窗口大小的一半)
--pvalue Pvalue cutoff for peak detection. DEFAULT:1e-5(确定峰的阀值)
--mfold Select the regions with MFOLD high-confidence enrichment
ratio against background to build model. Default:10,30
(建模时选双峰的阀值区间,需要不断调试,若找不到足够多的峰
来建模,可以适当扩大区间)
--wig Whether or not to save shifted raw tag count at every bp into a
wiggle file. WARNING: this process is time/space consuming!!
(选择是否生成*.wig文件)
--wigextend If set as an integer, when MACS saves wiggle files, it will extend
tag from its middle point to a wigextend size fragment. By
default it is modeled d. Use this option if you want to increase
the resolution in wiggle file. It doesn't affect peak calling.
(保存wiggle文件时由中心延伸的距离,能增加wiggle 文件的分
辨率且不影响找峰)
--space The resoluation for saving wiggle files, by default, MACS will save
the raw tag count every 10 bps. Usable only with '--wig' option.
(值越小,wiggle 文件的分辨率越高)
--nolambda If True, MACS will use fixed background lambda as local lambda
for every peak region.(是否设局部的lambda)
--lambdaset Three levels of nearby region in basepairs to calculate dynamic
lambda, DEFAULT: "1000,5000,10000"(算局部lambda时选的
三个区域)
--nomodel Whether or not to build the shifting model. If True, MACS will not
build model. by default it means shifting size = 100, try to set
shiftsize to change it.(是否建模,不建,则默认shifsize = 100)
--shiftsize The arbitrary shift size in bp. When nomodel is true, MACS will
regard this value as 'modeled' d.(承上,若不建模,可以设
Shiftsize为任意值,默认为100)










 三 找峰过程

1. 去重复
由于在DNA扩增或文库制备中会导致背景偏移,所以MACS首先利用二项分模型对tags进行去重复处理。假设tags落在基因组每个位置的概率是相同的(概率为$1/gsize$),所以以某确定个位置为起点的tags个数服从二项分布,即:
$$\sum _{ k=0 }^{ n }{ { C }_{ n }^{ k }{ p }^{ k }{1-p}^{n-k} }>1-P $$ 
$n$代表tags总数,$p=1/gsize$,P为设定的p-value值。上述公式的解$k$为以特定位置为起点的tags数的最大值。

2. 建模

如果使用nomodel参数或是没有找到足够的双峰(>100),MACS将不会建模,而使:
d=shiftsize*2(shiftsize默认为100)
scanwindow=2*d
如果没有nomodel参数计算,MACS将估计d值。

(1) 令: 
peaksize=2*bw。
min_tags=total_tags*lmfold*peaksize/gsize/2
max_tags=total_tags*umfold*peaksize/gsize/2
min_tags、max_tags 是每个peaksize长度范围内tags数的最小值、最大值。lmfold、umfold是-mfold参数值,total_tags是unique tags总数,bw为-bw参数值。

(2) 确定正负链上tags富集区
把tags按照染色体分开并排序,并按照正负链分开。在染色体范围内扫描并找到在peaksize长度内tags数在min_tags和max_tags之间的区域,合并有overlap的区域,计算此区域的最高点(tags深度最大),记录此区域中的tags。我们得到了所有正链所有的峰值pp1、pp2…,负链的所有峰值mp1,mp2…

(3) 确定每一对双峰的中心位置
正链的每一个最高点向左延伸peaksize,向右延伸peaksize要求延伸后的值符合下面的公式:
$pp_i - peaksize < mp_j$,或$pp_i + peaksize > mp_j$
($pp_i$ 为正链的最高点位置,$mp_j$为负链最高点的位置),把符合上面公式的$pp_i$,$mp_j$认为是一对双峰的最高点。从而得到一个双峰的中心$pp_i + mp_j / 2$。

(4) 求d值
在每一个双峰的中心位置向左、右分别延伸peaksize长度,统计正链tags在这个区间的每一点的深度,然后再统计负链tags在这个区间的每一点的深度,最后把所有的双峰中心附近peaksize区域每一点深度统计结果叠加,就得到了一对双峰。计算此双峰的两个最高值。d值为两个最高点的差值。如图:
chip-seq-peak.png


(5)shiftsize=d/2
  scan_window=2*max(tsize,d)
 
3. 获得富集区域(Call peaks)
 
1) lambda_bg = scan_window*total_tags/gsize
  rate=treat_total_tags/contral_total_tags
  要求rate要在0.5~2之间。


2)利用泊松分布计算出扫描窗口中最小tags数目
$$
\sum _{ k=0 }^{ n }{ \frac { { \lambda  }^{ k } }{ k! } { e }^{ -\lambda  } } \le P\quad 且\quad \sum _{ k=0 }^{ n+1 }{ \frac { { \lambda  }^{ k } }{ k! } { e }^{ -\lambda  } } \ge P
$$
$\lambda = \lambda_{bg},P = P_{value}$,上述公式所得的n值即为扫描窗口中应有的最小tags。

3)平移合并
  把正链的tags向又平移d/2长度,负链的tags向左平移d/2,然后把所有平移后的tags放在一起考虑。
4)以scan_windown长度扫描每条染色体,记录下窗口中tags大于n的窗口,合并有overlap的窗口,合并后的区域为候选的peak并统计出peak长度,peak最高点,tags数等。
5)如果有control样本,则对control样本进行3),4)处理,得到候选的假阳性peak。

4. peak筛选
 
1)有对照情况
a) 以每个候选peak的高点为中心,分别延伸lregion(默认200),sregion(默认1000)长度,如下:
left_lregion = peak_summit-lregion/2
right_lregion = peak_summit+lregion/2
left_sregion = peak_summit-sregion/2
right_sregion = peak_summit+sregion/2






在treat样本中统计出peak长度范围内的tags数num_peak,在control样本中统计出left_lregion ~right_lregion范围内的tags数cnum_l,left_sregion ~right_sregion范围内的tags数cnum_s:
tlambda_peak=num_peak/peak_length* max(peak_length,scan_window)
lambda_tbg=lambda_bg/scan_window*max(peak_length,scan_window)
clambda_lregion=cnum_l/lregion*rate* max(peak_length,scan_window)
clambda_sregion=cnum_s/sregion*rate* max(peak_length,scan_window)
lambda_lcoal=max(lambda_tbg,clambda_lregion,clambda_sregion)
认为此peak范围内的tags总数应该符合$\lambda = \lambda_{lcoal}$的泊松分布,所以:
$$1-\sum _{ k=0 }^{ n }{ \frac { { \lambda_{local}  }^{ k } }{ k! } { e }^{ -\lambda_{local}  } } =\quad P_{tmp}$$
$n =\lambda_{peak}$。如果$P_{tmp}$ 大于 $P_{value}$ 则认为此peak为不可信,过滤掉此peak。$P_{tmp}$ 小于 $P_{value}$ 认为此peak可信。最终得到所有可信的peak。

b) 假阳性计算

把control样本看成treat样本,把treat样本看成control样本并进行a)步处理,那么就可以得到所有假阳性的peak,所以FDR定义为:
FDR=100 * Num_control/Num_teat
Num_teat为在peak中大于此peak的P值的peak的数目
Num_control为在假阳性peak中大于此peak的P值的peak的数目
c) 富集度的计算
fold_enrichment = peak_height  /  local_lambda *  max(peak_length,scan_window) / d
peak_height为peak中最高的深度值。
2) 如果没有control样本
以每个候选peak的高点为中心,延伸tregion(默认5000)长度。
left_tregion = peak_summit-tregion/2
left_tregion = peak_summit-tregion/2
tlambda_tregion=tnum_t/tregion*max(peak_length,scan_window)
tnum_t为treat 样本中lregion范围内tags的总数。
lambda_lcoal=max(lambda_bg,tlambda_tregion),与有control样本的方法相同,计算出P_tmp的值并与P_Value比较判断此peak是否可信。没有control样本则无法计算FDR。

四 输出文件
 
输出文件格式包括 *.xls  *.bed  (峰的具体信息如位置、峰值、pvalue值以及富集度等等) *.log (日志文件,记录MACS找峰的具体过程信息) *.r(MACS自动生成的建模用的R脚本)以及一个保存tag序列的压缩包文件夹。
下面列出一个.xls的peak信息文件来加以说明:
Peak-info.png

前面加“#”号的号为说明文档,下面开始的各列的最上方那排标题代表peak的各信息含义:
chr:             染色体位置
start: peak的起始位点
end: peak的终止位点
length: peak的长度
summit: peak的最高点位置
tags: peak区域的tag数
-10*log10(pvalue):pvalue值开以10为底的对数后乘以负10的值
fold_enrichment: peak富集度。它是tag数与peak区域平均tag水平的比值
FDR(%): control peak 数 /ChIP peak 数[size=13] [/size]
参考链接:

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3120977/


0 个评论

要回复文章请先登录注册