1. Call peaks 并去除blacklists
sample=SRR6322534
## callpeak
input=SRR6322538
gsize=199000000 # 190M
macs2 callpeak -t ./$sample.ff.bam -c ./$input.ff.bam -f BAM -n $sample -g $gsize --keep-dup all
# filter peaks against blacklist
bedtools intersect -a ${sample}_peaks.narrowPeak -b ../ref/blacklist.bed -f 0.25 -v >${sample}_peaks.f.narrowPeak
2. 计算文库复杂度 #注意这里采用的是未过滤的bam files。
image.png
thread=2 #设置线程数目，要设置，后面用到负责会报错。
#打印表头
echo -e "#TotalReads\tDistinctReads\tOneRead\tTwoReads\tNRF=Distinct/Total\tPBC1=OneRead/Distinct\tPBC2=OneRead/TwoReads" >$sample.pbc.metrics.txt
#计算one reads two reads 的数目 #主要考虑的是坐标起始是否一致，如果一致则是重复
bedtools bamtobed -i $sample.qf.bam | awk 'BEGIN{OFS="\t"}{if($6=="+"){print $1,$2,$6} else{print $1,$3,$6}}' | sort --parallel=$thread | uniq -c | awk 'BEGIN{mt=0;m0=0;m1=0;m2=0} ($1==1){m1=m1+1} ($1==2){m2=m2+1} {m0=m0+1} {mt=mt+$1} END{printf "%d\t%d\t%d\t%d\t%f\t%f\t%f\n",mt,m0,m1,m2,m0/mt,m1/m0,m1/m2}' >> $sample.pbc.metrics.txt
cat $sample.pbc.metrics.txt
#TotalReads DistinctReads OneRead TwoReads NRF=Distinct/Total PBC1=OneRead/Distinct PBC2=OneRead/TwoReads
1521604 1492475 1464416 27418 0.980856 0.981200 53.410752
image.png
3. 计算FRIPs
image.png
## FRiP
all_c=$(sed -n 's/^# .* in treatment: \([0-9]\)/\1/p' ${sample}_peaks.xls) #从peaks.xls文件中计算总的reads counts
peak_c=$(samtools view -L ${sample}_peaks.narrowPeak $sample.ff.bam | wc -l) #计算bam 文件中落在peaks区域的reads counts
echo -e "#Total reads\tPeaks reads\tFRiP" >$sample.FRiP.metrics.txt #保存表头
echo -e "${all_c}\t${peak_c}" | awk '{print $1"\t"$2"\t"$2/$1}' >>$sample.FRiP.metrics.txt #计算FRIPs值 并保村
cat $sample.FRiP.metrics.txt
#Total reads Peaks reads FRiP
1469342 15114 0.0102862
4. 计算cross co-rrelation
image.png
## cross correlation
# this need to install phantompeakqualtools, see https://github.com/kundajelab/phantompeakqualtools
sample=SRR6322534
gsize=199000000
run_spp.R -rf -c=${sample}.ff.bam -p=4 -s=-500:5:500 -savp=${sample}.cc.pdf -out=${sample}.cc.metrics.txt
#展示结果
cat SRR6322534.cc.metrics.txt
SRR6322534.ff.bam 1469342 115,130,170 0.197323712213789,0.197220698531932,0.197185630044492 80 0.1973325 500 0.1941522 1.016335 0.9972433 0
image.png
一把来说length of fragment 选择第一个即可
注释
image.png
参考：基因课------表观基因组学之 ChIP-seq 数据分析
