【ChIP-seq 实战】六、使用macs2找peaks

作者: 佳奥 | 来源:发表于2022-08-12 10:20 被阅读0次

【ChIP-seq 实战】六、使用macs2找peaks
ChIP-seq数据分析实战训练（三）
第3篇：用MACS2软件call peaks
Peak annotation
ATAC-seq(4) -- 使用macs3找peaks
MACS2 进行peak calling
ChIP-seq（三）MACS2 peak calling
m6A上游 | 5 MACS2 peak calling
ATAC-seq 分析（上）
2022-09-13 calling peaks问题

这里是佳奥！

ChIP-Seq分析进入关键步骤。让我们开始吧！

##合并后的bam以及去除PCR重复的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，官方提供了使用实例

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

介绍一下各个参数的意义：

-t: 实验组的输出结果treat
-c: 对照组的输出结果control
-f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一个。如果不提供这项，就是自动检测选择。
-g: 基因组大小，默认提供了hs, mm, ce, dm选项，不在其中的话，比如说拟南芥，就需要自己提供了。
-n: 输出文件的前缀名
-B: 会保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
-q: q值，也就是最小的PDR阈值，默认是0.05。q值是根据p值利用BH计算，也就是多重试验矫正后的结果。
-p：这个是p值，指定p值后MACS2就不会用q值了。
-m: 和MFOLD有关，而MFOLD和MACS预构建模型有关，默认是5：50，MACS会先寻找100多个peak区构建模型，一般不用改，因为你很大概率上不会懂。

我们实战使用的代码：

cd  ~/mergeBam 

##合并的bam文件运行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重复的.rmdup.bam文件运行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

##查看结果 合并的bam（运行macs2的结果会跑到peaks目录下，这里我改名了，为了分开两边的结果）
(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
略

##查看比较的统计信息
$ 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 总用量

##查看结果 合并并去除PCR重复的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
略

##查看比较的统计信息，去除后结果略有差别
$ 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 总用量

下一步便是可视化操作，我们下一篇再见！