IDR分析

由于IDR分析比较严格,所以要是打算进行IDR分析的话,call peaks时候 设置阈值时候选择p值,且p < 0.01即可
#设置变量
sample_r1=SRR6322534
sample_r2=SRR6322535
sample=SRR6322534_SRR6322535
input=SRR6322538

#设置基因组大小
gsize=199000000

## call peaks for replicte 1
macs2 callpeak -t ./${sample_r1}.ff.bam -c ./$input.ff.bam -f BAM -n ${sample_r1} -g $gsize --keep-dup all -p 0.01

## call peaks for replicte 2
macs2 callpeak -t ./${sample_r2}.ff.bam -c ./$input.ff.bam -f BAM -n ${sample_r2} -g $gsize --keep-dup all -p 0.01

## call peaks for combined dataset,这里之间把两个bams 文件提供即可,不需要合并bams
macs2 callpeak -t ./${sample_r1}.ff.bam ./${sample_r2}.ff.bam -c ./$input.ff.bam -f BAM -n ${sample} -g $gsize --keep-dup all -p 0.01

## run idr
idr --samples ${sample_r1}_peaks.narrowPeak ${sample_r2}_peaks.narrowPeak --peak-list ${sample}_peaks.narrowPeak --input-file-type narrowPeak --output-file ./${sample}.idr --rank p.value --soft-idr-threshold 0.05 --plot

## get idr produced narrowPeak file,把结果转换成narrow peaks
cut -f 1-10 ${sample}.idr | sort -k1,1 -k2,2n -k3,3n >${sample}.idr.narrowPeak

## filter peaks against blacklist
bedtools intersect -a ${sample}.idr.narrowPeak -b ../ref/blacklist.bed -f 0.25 -v >tmp && mv tmp ${sample}.idr.narrowPeak

参考:基因课------表观基因组学之 ChIP-seq 数据分析

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成,浏览时请结合常识与多方信息审慎甄别。
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

相关阅读更多精彩内容

友情链接更多精彩内容