【ChIP-seq分析】超级增强子系列4： 用bedtools来进行共识SE的分析

#策略，非真实运行的代码
#1将所有组内样本peaks信息汇集到一个文件
cat *.bed >allHC.bed 
#排序汇集的peaks文件
，去重并第四列peakID
sort -k1,1 -k2,2n allHC.bed | awk '!seen[$0]++' | awk 'BEGIN{OFS="\t"} {print $0, "peak" NR}' > allHCsorted.bed
#分析这个汇集组内所有样本的peaks与不同样本之间的交集情况
bedtools intersect -a allHC.bed -b HC1.bed  HC2.bed HC3.bed  -wa -wb > allHCacrossallsamples.bed
结果举例如下
chr1	100	200	peak1	chr1	90	190	HC1_peak1	HC1 
chr1	100	200	peak1	chr1	110	210	HC2_peak1	HC2
chr1	100	200	peak1	chr1	105	205	HC3_peak1	HC3
chr1	300	400	peak2	chr1	300	400	HC2_peak2	HC2 
chr1	300	400	peak2	chr1	310	410	HC3_peak2	HC3  
chr1	500	600	peak3	chr1	490	590	HC1_peak2	HC1
chr1	900	1000	peak4	chr1	890	990	HC1_peak3	HC1
chr1	900	1000	peak4	chr1	895	995	HC3_peak3	HC3
# 统计每个 allHC peak 出现的次数（即支持样本数）
cut -f4 allHCacrossallsamples.bed | sort | uniq -c >peakoverlapumber.txt
输出结果如下
3 peak1
2 peak2
1 peak3
2 peak4

bedtools merge HC3sample.bed  -i -d 12500 > HC_SEmerged_regions.bed

# 1.排序：批量排序
for f in *.bed; do sort -k1,1 -k2,2n "$f" > "${f%.bed}.sorted.bed"; done
# 2.然后运行bedtools multiinner进行重叠区域分析
bedtools multiinter -header -names HC1 HC2 HC3 -i *_SE.sorted.bed > HCmulti_SE.txt 

awk '$4 > 1' HCmulti_SE.txt | cut -f1-3 | bedtools merge -i -d 12500 > HC_SEmerged_regions.bed