MACS2是peak calling最常用的工具。
callpeak用法
這是MACS2的主要功能，因為MACS2的目的就是找peak，其他功能都是可有可無，唯獨callpeak不可取代。最簡單的用法就是
# 常規(guī)的peak calling
macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
# 較寬的peak calling
macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1
我們先來介紹這個案例里的參數(shù)。首先是常規(guī)的peak calling用到的參數(shù)
-
-t/--treatment FIELNAME和-c/--control FILENAME表示處理樣本和對照樣本輸入。其中-t必須，很好理解，沒有處理組你還找啥Peak。 -
-f/--format FORMAT用來聲明輸入的文件格式，目前MACS能夠識別的格式有 "ELAND", "BED", "ELANDMULTI", "ELANDEXPORT", "ELANDMULTIPET" (雙端測序), "SAM", "BAM", "BOWTIE", "BAMPE", "BEDPE". 除"BAMPE", "BEDPE"需要特別聲明外，其他格式都可以用AUTO自動檢測。 -
-g表示實際可比對的基因組大小。比如說人類是2.7e9，也就是2.7G，而實際人類基因組大概是3.2G左右。這是因為有些地方無法拼接，會用N代替，這部分區(qū)域大概是80%左右。擬南芥根據(jù)NCBI顯示是119,667,750，那么實際能比對大概也就是1.0e8. NBT有一篇文章"PeakSeq enables systematic scoring of ChIP-seq experiments relative to controls"的表1就進(jìn)行了統(tǒng)計。we
-
-n/--name表示實驗的名字, 請取一個有意義的名字。 -
-B/--bdg: 以bedGraph格式存放fragment pileup, control lambda, -log10pvalue 和log10qvale.-
NAME_treat_pileup.bdg: 處理后數(shù)據(jù) -
NAME_control_lambda.bdg： 對照的局部lambda值 -
NAME_treat_pvalue.bdg： 泊松檢驗的P值 -
NAME_treat_qvalue.bdg：Benjamini–Hochberg–Yekutieli處理后的Q值
-
-
-q： q值(最小的FDR)的閾值，默認(rèn)0.05。可以根據(jù)結(jié)果進(jìn)行修正。q值是p值經(jīng)Benjamini–Hochberg–Yekutieli修正后的值。
一般常規(guī)是夠用的，但是如果你需要看那些更加寬的peak，可以按照官方的建議使用如下參數(shù)
-
--broad: broad region最大長度是4d。其中d表示MACS的雙峰模型兩個peak的距離。結(jié)果會得到BED12格式文件，存放著附近高度附近的區(qū)域。由于要足夠的寬，所以需要專門的參數(shù)進(jìn)行統(tǒng)計學(xué)過濾。 -
--broad--cutoff: 用于過濾broad得到的peak，默認(rèn)是q值，如果設(shè)置-p就用p值。
上面的基本參數(shù)可以用在最初的分析。根據(jù)基本的分析結(jié)果，可以有選擇地使用下面的參數(shù)符合特定項目的需求。
比較基礎(chǔ)的參數(shù):
-
-s/--tsize: 二代測序讀長，MACS會用前面10個序列進(jìn)行推測。 -
--outdir： 輸出文件夾 -
--verbose 0/1/2/3: 輸出信息的詳細(xì)度。如果是0就表示不想看到屏幕有輸出。 -
-p/--pvalue： 使用P值，而不是q值，也就是說用未多重矯正的p值進(jìn)行篩選。 -
--to-large： 默認(rèn)是把大樣本縮小和小樣本一樣小，設(shè)置該參數(shù)則是把小樣本放大成大樣本一樣大。
和MACS模型構(gòu)建相關(guān)的參數(shù)。正如MACS的名字所示, Model-based Analysis of ChIP-seq, 它需要先建立模型然后分析。那么問題來了，建立什么模型？模型的目的是什么？這里的模型指的是雙峰模型，建立雙模模型的目的是為了更好的將原始reads朝3'偏移，更好的表示蛋白和DNA的作用位置。這里還要多問一句為什么要偏倚。
這就需要從實驗建庫說起。ChIP-seq目標(biāo)是找到蛋白和DNA的作用位置，所以先要讓蛋白和DNA進(jìn)行交聯(lián)，之后用超聲打碎，再用抗體把與蛋白結(jié)合的DNA收集起來測序。在MACS發(fā)表的2008年，那個時候的測序大多都以單端50bp為主，而超聲破碎的片段肯定大于50 bp（這可以通過電泳圖來了解），也就是說最開始的SE50數(shù)據(jù)比對到參考基因組之后，得到的峰圖并沒有真實反映出原來的文庫情況。但由于比對到基因組正負(fù)鏈的概率是相似的，那么就會形成兩個峰（如下圖），為了更好的還原出最來的文庫片段，就先建立了雙峰模型，以兩峰距離d的一半作為偏倚長度。
如果你的數(shù)據(jù)是SE50或者SE100，為了更準(zhǔn)確地找peak，你需要建立雙峰模型，你可能要調(diào)整--bw, --mfold, --fix-bimodal, --shift, --extsize。 但是對于雙端測序而言，它本身測的就是文庫的兩端，因此建立模型沒有必要，偏倚也沒有必要，你 只需要 設(shè)置參數(shù)--nomodel。
-
--bw: 這個參數(shù)僅僅當(dāng)你知道ChIP實驗中超聲打斷后的條帶長度時才可能需要設(shè)置。用來構(gòu)建雙峰模型。 -
--nomodel: 這個參數(shù)說明不需要MACS去構(gòu)建模型，也就是說下面的參數(shù)除了--shift, --extsize外都會被無視。 -
--extsize： MACS使用這個參數(shù)將read以5'-> 3'衍生至等長片段。比如說你知道你的轉(zhuǎn)錄因子的結(jié)合區(qū)域是200bp，那么參數(shù)就是--extsize 200。當(dāng)且僅當(dāng)--nomodel和--fix-bimodal設(shè)置使用。 -
--shift: 這個參數(shù)是絕對的偏移值，會先于--extsize前對read進(jìn)行整體移動。MACS會通過建模的方式自動計算出read需要偏移的距離，除非你對自己的數(shù)據(jù)非常了解，或者前期研究都表明結(jié)合中心在read后面的那個位置上，你才能比較放心的用這個這個參數(shù)了。正數(shù)表示從5'往3'偏移延長到片段中心，如果是負(fù)數(shù)則是3'往5'偏移延長到片段中心。作者給了幾個例子：- 如果是ChIP-seq數(shù)據(jù)，設(shè)置·--shift 0`
- 如果是
DNase-Seq數(shù)據(jù)：read來自于兩個核小體中間，你想把測序read往兩邊延長用來平滑pileup信號，并且希望用來平滑的窗口是200bp,那么使用`--nomodel --shift -100 --extsize 200'. - 如果是
nucleosome-seq數(shù)據(jù)：因為一個核小體大概有147bp DNA纏繞，于是就需要用半個核小體長度進(jìn)行堆積(pipleup)用于小波分析。參數(shù)為--nomodel --shift 37 --extsize 73.
-
-m/--mfold: 構(gòu)建雙峰模型時使用，默認(rèn)是[5,50]，表示選擇那些倍數(shù)變化在5~10之間的富集區(qū)域。如果找不到100個區(qū)域構(gòu)建模型，并且你還設(shè)置了--fix-bimodal時，它就會用--extsize參數(shù)繼續(xù)分析 -
--nolambda： 設(shè)置這個參數(shù)就意味著不用MACS推薦的動態(tài)lambda，而是使用背景l(fā)ambda作為local lambda，也就是不考慮染色質(zhì)結(jié)構(gòu)等造成的局部偏誤。 -
--slocal, --llocal: 這兩個參數(shù)也是MACS用來計算動態(tài)lambda會用到，分別計算1kb內(nèi)lambda(slocal)和10kb的lambda(llocal)，目標(biāo)是處理類似于開放染色質(zhì)區(qū)域的效應(yīng)。注，如果這兩個參數(shù)太小，輸入數(shù)據(jù)中的尖峰(sharp spike)就可能干掉顯著性的peak。
謹(jǐn)慎使用的參數(shù)：
-
--down-sample：如果你的電腦性能比較差，或者樣本特別大，你希望快點看到一個差不多的結(jié)果，可以使用這個參數(shù)。MACS會對數(shù)據(jù)進(jìn)行隨機(jī)抽樣，所以每次的結(jié)果會不太一樣。如果結(jié)果是要發(fā)文章，不要用這個參數(shù)得到的結(jié)果。 -
--keep-dup: 保留重復(fù)。默認(rèn)MACS(auto)會使用二項分布估計每個位置上是否存在重復(fù)（默認(rèn)是1，也就是每個位置上出現(xiàn)一個read的概率最大）。如果你前期已經(jīng)去重，那就使用all省了這一步.
callpeak結(jié)果文件說明
callpeak會得到如下文件：
- NAME_peaks.xls: 以表格形式存放peak信息，雖然后綴是xls，但其實能用文本編輯器打開，和bed格式類似，但是以1為基，而bed文件是以0為基.也就是說xls的坐標(biāo)都要減一才是bed文件的坐標(biāo)
- NAME_peaks.narrowPeak NAME_peaks.broadPeak 類似。后面4列表示為， integer score for display， fold-change，-log10pvalue，-log10qvalue，relative summit position to peak start。內(nèi)容和NAME_peaks.xls基本一致，適合用于導(dǎo)入R進(jìn)行分析。
- NAME_summits.bed：記錄每個peak的peak summits，話句話說就是記錄極值點的位置。MACS建議用該文件尋找結(jié)合位點的motif。能夠直接載入UCSC browser，用其他軟件分析時需要去掉第一行。
- NAME_peaks.gappedPeak: 格式為BED12+3，里面存放broad region和narrow peaks。
- NAME_model.r，能通過
$ Rscript NAME_model.r作圖，得到是基于你提供數(shù)據(jù)的peak模型。 - .bdg文件能夠用UCSC genome browser轉(zhuǎn)換成更小的bigWig文件。
其他有用的子命令
bdgcmp使用*_treat_pileup.bdg和*_control_lambda.bdg計算得分軌(score track)
bdgpeakcall使用 *_treat_pvalue.bdg 或bdgcmp得到的結(jié)果或begGraph文件進(jìn)行peak calling.bdgbroadcall差不多也是這樣子。
bdgdiff能用來分析4個bedgraph文件,得到treatment1 vs control1, treatment2 vs control2, treatment1 vs control2, treament2 vs control1的得分。
filterdup：過濾重復(fù)，結(jié)果是BED文件
predictd：從比對文件中估計文庫大小或d
randsample： 隨機(jī)抽樣
pileup：以給延伸大小去堆積(pileup)比對得到的reads。這一步不會有去重和測序深度標(biāo)準(zhǔn)化，你需要預(yù)先做這些工作。
