接上回,得到 bedgraph 文件了。
重要秘籍: In-depth-NGS-Data-Analysis
6 Peak calling
- M1 SEACR
mkdir -p $projPath/seacr
seacr="/lustre/home/yfxie/.conda/envs/cuttagg/bin/SEACR_1.3.sh"
cat $config1 | while read line
do
arr=($line)
sample=${arr[0]}
## top0.01和top0.05,因为没有做IgG样本
bash $seacr $projPath/alignment/bedgraph/$sample.fragments.normalized.bedgraph \
0.01 non stringent $projPath/seacr/$sample.seacr_top0.01.peaks
bash $seacr $projPath/alignment/bedgraph/$sample.fragments.normalized.bedgraph \
0.05 non stringent $projPath/seacr/$sample.seacr_top0.05.peaks
done
- M2 MACS2
这里我想直接得出差异peaks,合并rep1和rep2,参考:Build Signal Track,2021-05-23 ChIP-seq数据从头到尾比对及分析汇总(个人分析记录贴)
mkdir -p $projPath/macs2
macs2 callpeak -t ${projPath}/alignment/bam/A_rep1.name.rmDup.bam ${projPath}/alignment/bam/A_rep2.name.rmDup.bam \
-c ${projPath}/alignment/bam/W_rep1.name.rmDup.bam ${projPath}/alignment/bam/W_rep2.name.rmDup.bam\
-g mm -f BAMPE -n A_macs2_q0.1.peaks --outdir $projPath/macs2 -q 0.1 --keep-dup all 2>${projPath}/macs2/A_macs2.log