安裝默認使用Conda
macs2的基本命令
macs2 callpeak -t treatment.bam -c control.bam -g hs -B -f BAM -n prefix -q 0.00001
macs2 callpeak -t treatment.bam -c control.bam -g hs -B -f BAMPE -n prefix -q 0.00001
參數(shù)解讀：
-g 基因組的選擇
-B 輸出bgd文件，下游bigwig文件生成所需
-f 雙端測序使用BAMPE出爹，單端的話不需要加參數(shù)像棘，默認是auto識別，但是識別不了BAMPE
-q 使用閾值蜓谋，根據(jù)需要選擇
輸出文件包括
prefix_model.r
prefix_peaks.narrowPeak
prefix_peaks.xls
prefix_summits.bed
prefix_treat_pileup.bdg
prefxi_ _treat_pvalue.bdg
prefix_model.r 為模型建立示意圖，使用Rscript prefix_model.r 命令可以在Linux環(huán)境下作圖
prefix_peaks.xls 為包含了peak信息的表格丽已，可以用excel打開，其中包括：
染色體
peak起始位點
peak終止位點
peak區(qū)域的長度
peak定點的絕對位置
堆積信號值及其-log10(pvalue)
fold enrichment及其-log10(pvalue)
#注意：XLS文件里面的起始位點和GFF文件相似從1起始买决，bed文件從0開始
prefix_peaks.narrowPeak 為包含了peak位置及summit位置沛婴，p值吼畏，q值的文件，可以直接上傳到UCSC browser嘁灯，其中特定列的信息為：
5th：整數(shù)信號值
7th：fold-change
8th：-log10 pvalue
9th：-log10 qvalue
10th：summit距離peak起始位點的距離
prefix_summits.bed 包含每個peak峰頂summit的位置泻蚊，等同于從narrowPeak文件里面剝離出來的，方便用來尋找motif的binding丑婿，同樣可以載入UCSC
prefix_peaks.broadPeak &prefix_peaks.gappedPeak 文件都是為選擇 -broad 參數(shù)后得到的文件性雄，本質上和narrowPeak文件一致，不過增加了broadregion羹奉，當然summit的定義也有區(qū)別秒旋，需要自己去指定
prefix_treat_pileup.bdg 文件為輸入文件treat組的bedGraph file， prefxi_ _treat_pvalue.bdg為對照組的bedGraph file
為了方便在IGV上查看ChIP-seq的結果和后期的可視化展示诀拭，所以我們需要把macs2的結果轉化為bw提供給IGV(bdg file to wig file transformation)
一共分為三步
1.首先使用bdgcmp得到FE或者logLR轉化后的文件(Run MACS2 bdgcmp to generate fold-enrichment and logLR track)
macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_FE.bdg -m FE
macs2 bdgcmp -t H3K36me1_EE_rep1_treat_pileup.bdg -c H3K36me1_EE_rep1_control_lambda.bdg -o H3K36me1_EE_rep1_logLR.bdg -m logLR -p 0.00001
# 參數(shù)解釋
-m FE 計算富集倍數(shù)降低噪音
-p 為了避免log的時候input值為0時發(fā)生error迁筛，給予一個很小的值
2. 預處理文件與對應參考基因組染色體長度文件下載
使用conda安裝以下三個處理軟件：bedGraphToBigWig/bedClip/bedtools
下載染色體長度文件：http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/chromInfo.txt.gz 并解壓（針對human，其余物種的皆可以按照類似網(wǎng)址下載）
寫一個小小的sh腳本方便一步轉化name.sh：
#!/bin/bash
# check commands: slopBed, bedGraphToBigWig and bedClip
which bedtools &>/dev/null || { echo "bedtools not found! Download bedTools: <http://code.google.com/p/bedtools/>"; exit 1; }
which bedGraphToBigWig &>/dev/null || { echo "bedGraphToBigWig not found! Download: <http://hgdownload.cse.ucsc.edu/admin/exe/>"; exit 1; }
which bedClip &>/dev/null || { echo "bedClip not found! Download: <http://hgdownload.cse.ucsc.edu/admin/exe/>"; exit 1; }
# end of checking
if [ $# -lt 2 ];then
echo "Need 2 parameters! <bedgraph> <chrom info>"
exit
fi
F=$1
G=$2
bedtools slop -i ${F} -g ${G} -b 0 | bedClip stdin ${G} ${F}.clip
LC_COLLATE=C sort -k1,1 -k2,2n ${F}.clip > ${F}.sort.clip
bedGraphToBigWig ${F}.sort.clip ${G} ${F/bdg/bw}
rm -f ${F}.clip ${F}.sort.clip
chmod +x name.sh
3. 最后就是生成bw文件
./name.sh H3K36me1_EE_rep1_FE.bdg hg19.len
./name.sh H3K36me1_EE_rep1_logLR.bdg hg19.len
最后得到產(chǎn)物耕挨，至于的使用哪一個作為輸入文件大家就根據(jù)需要來吧
H3K36me1_EE_rep1_FE.bw
H3K36me1_EE_rep1_logLR.bw
參考文獻：
