這里是佳奧！

ChIP-Seq分析進(jìn)入關(guān)鍵步驟。讓我們開始吧！

##合并后的bam以及去除PCR重復(fù)的rmdup.bam
(chipseq) root 21:02:33 /home/kaoku/chipseq/mouse_project/mergeBam
$ ls -lh
總用量 18G
-rw-r--r-- 1 root root  836M  8月 11 15:57 Control.merge.bam
-rw-r--r-- 1 root root  743M  8月 11 20:57 Control.merge.rmdup.bam
-rw-r--r-- 1 root root  1.3G  8月 11 15:58 H2Aub1.merge.bam
-rw-r--r-- 1 root root  1.1G  8月 11 20:58 H2Aub1.merge.rmdup.bam
-rw-r--r-- 1 root root  1.6G  8月 11 15:59 H3K36me3.merge.bam
-rw-r--r-- 1 root root  1.5G  8月 11 20:58 H3K36me3.merge.rmdup.bam
-rw-r--r-- 1 root root  1.3G  8月 11 16:00 Ring1B.merge.bam
-rw-r--r-- 1 root root 1006M  8月 11 20:58 Ring1B.merge.rmdup.bam
-rw-r--r-- 1 root root  1.1G  8月 11 16:01 RNAPII_8WG16.merge.bam
-rw-r--r-- 1 root root  984M  8月 11 20:58 RNAPII_8WG16.merge.rmdup.bam
-rw-r--r-- 1 root root  1.5G  8月 11 16:02 RNAPII_S2P.merge.bam
-rw-r--r-- 1 root root  1.2G  8月 11 20:58 RNAPII_S2P.merge.rmdup.bam
-rw-r--r-- 1 root root  1.4G  8月 11 16:03 RNAPII_S5P.merge.bam
-rw-r--r-- 1 root root  775M  8月 11 20:58 RNAPII_S5P.merge.rmdup.bam
-rw-r--r-- 1 root root  215M  8月 11 16:04 RNAPII_S5PRepeat.merge.bam
-rw-r--r-- 1 root root  210M  8月 11 20:57 RNAPII_S5PRepeat.merge.rmdup.bam
-rw-r--r-- 1 root root  713M  8月 11 16:04 RNAPII_S7P.merge.bam

macs2包含一系列的子命令，其中最主要的就是callpeak， 官方提供了使用實(shí)例

macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

介紹一下各個(gè)參數(shù)的意義：

• -t: 實(shí)驗(yàn)組的輸出結(jié)果treat
• -c: 對(duì)照組的輸出結(jié)果control
• -f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一個(gè)。如果不提供這項(xiàng)，就是自動(dòng)檢測(cè)選擇。
• -g: 基因組大小， 默認(rèn)提供了hs, mm, ce, dm選項(xiàng)， 不在其中的話，比如說擬南芥，就需要自己提供了。
• -n: 輸出文件的前綴名
• -B: 會(huì)保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
• -q: q值，也就是最小的PDR閾值， 默認(rèn)是0.05。q值是根據(jù)p值利用BH計(jì)算，也就是多重試驗(yàn)矯正后的結(jié)果。
• -p： 這個(gè)是p值，指定p值后MACS2就不會(huì)用q值了。
• -m: 和MFOLD有關(guān)，而MFOLD和MACS預(yù)構(gòu)建模型有關(guān)，默認(rèn)是5：50，MACS會(huì)先尋找100多個(gè)peak區(qū)構(gòu)建模型，一般不用改，因?yàn)槟愫艽蟾怕噬喜粫?huì)懂。

我們實(shí)戰(zhàn)使用的代碼：

cd  ~/mergeBam 

##合并的bam文件運(yùn)行macs2
ls  *merge.bam |cut -d"." -f 1 |while read id;
do 
    if [ ! -s ${id}_summits.bed ];
    then 
echo $id 
macs2 callpeak -c  Control.merge.bam -t $id.merge.bam -f BAM -B -g mm -n $id --outdir ../peaks  2> $id.log &  
    fi 
done  

##把兩邊文件分開避免干擾
mkdir dup
mv *rmdup* dup/
cd dup/

##合并并去除PCR重復(fù)的.rmdup.bam文件運(yùn)行macs2
ls  *.merge.rmdup.bam |cut -d"." -f 1 |while read id;
do 
    if [ ! -s ${id}_rmdup_summits.bed ];
    then 
echo $id 
macs2 callpeak -c  Control.merge.rmdup.bam  -t $id.merge.rmdup.bam  -f BAM -B -g mm -n ${id}_rmdup --outdir ../peaks 2> $id.log &  
    fi 
done  

##查看結(jié)果 合并的bam（運(yùn)行macs2的結(jié)果會(huì)跑到peaks目錄下，這里我改名了，為了分開兩邊的結(jié)果）
(chipseq) root 22:05:35 /home/kaoku/chipseq/mouse_project/peaksMerge
$ ls -lh
總用量 17G
-rw-r--r-- 1 root root 1006M  8月 11 21:53 Control_control_lambda.bdg
-rw-r--r-- 1 root root   77K  8月 11 21:48 Control_model.r
-rw-r--r-- 1 root root     0  8月 11 21:53 Control_peaks.narrowPeak
-rw-r--r-- 1 root root  1013  8月 11 21:53 Control_peaks.xls
-rw-r--r-- 1 root root     0  8月 11 21:53 Control_summits.bed
-rw-r--r-- 1 root root  737M  8月 11 21:53 Control_treat_pileup.bdg
-rw-r--r-- 1 root root  926M  8月 11 21:54 H2Aub1_control_lambda.bdg
-rw-r--r-- 1 root root   77K  8月 11 21:49 H2Aub1_model.r
-rw-r--r-- 1 root root   79K  8月 11 21:54 H2Aub1_peaks.narrowPeak
-rw-r--r-- 1 root root   91K  8月 11 21:54 H2Aub1_peaks.xls
-rw-r--r-- 1 root root   53K  8月 11 21:54 H2Aub1_summits.bed
-rw-r--r-- 1 root root  1.1G  8月 11 21:54 H2Aub1_treat_pileup.bdg
略

##查看比較的統(tǒng)計(jì)信息
$ wc -l *bed
       0 Control_summits.bed
    1115 H2Aub1_summits.bed
   40830 H3K36me3_summits.bed
   26053 Ring1B_summits.bed
   41864 RNAPII_8WG16_summits.bed
   20042 RNAPII_S2P_summits.bed
   38663 RNAPII_S5PRepeat_summits.bed
   62765 RNAPII_S5P_summits.bed
   72640 RNAPII_S7P_summits.bed
  303972 總用量

##查看結(jié)果 合并并去除PCR重復(fù)的bam（peaks目錄是我新建的）
(chipseq) root 22:22:28 /home/kaoku/chipseq/mouse_project/peaks
$ ls -lh
總用量 17G
-rw-r--r-- 1 root root 1006M  8月 11 22:20 Control_rmdup_control_lambda.bdg
-rw-r--r-- 1 root root   77K  8月 11 22:15 Control_rmdup_model.r
-rw-r--r-- 1 root root     0  8月 11 22:20 Control_rmdup_peaks.narrowPeak
-rw-r--r-- 1 root root  1.1K  8月 11 22:20 Control_rmdup_peaks.xls
-rw-r--r-- 1 root root     0  8月 11 22:20 Control_rmdup_summits.bed
-rw-r--r-- 1 root root  737M  8月 11 22:20 Control_rmdup_treat_pileup.bdg
-rw-r--r-- 1 root root  926M  8月 11 22:20 H2Aub1_rmdup_control_lambda.bdg
-rw-r--r-- 1 root root   77K  8月 11 22:16 H2Aub1_rmdup_model.r
-rw-r--r-- 1 root root   86K  8月 11 22:20 H2Aub1_rmdup_peaks.narrowPeak
-rw-r--r-- 1 root root   98K  8月 11 22:20 H2Aub1_rmdup_peaks.xls
-rw-r--r-- 1 root root   59K  8月 11 22:20 H2Aub1_rmdup_summits.bed
-rw-r--r-- 1 root root  1.1G  8月 11 22:20 H2Aub1_rmdup_treat_pileup.bdg
略

##查看比較的統(tǒng)計(jì)信息，去除后結(jié)果略有差別
$ wc -l *bed
       0 Control_rmdup_summits.bed
    1115 H2Aub1_rmdup_summits.bed
   40830 H3K36me3_rmdup_summits.bed
   26053 Ring1B_rmdup_summits.bed
   41841 RNAPII_8WG16_rmdup_summits.bed
   20020 RNAPII_S2P_rmdup_summits.bed
   38663 RNAPII_S5PRepeat_rmdup_summits.bed
   62765 RNAPII_S5P_rmdup_summits.bed
   72577 RNAPII_S7P_rmdup_summits.bed
  303864 總用量

下一步便是可視化操作，我們下一篇再見！