【Hi-C】交互频率与距离衰减曲线绘制

参考：

发现AB区室 09 science Comprehensive mapping of long range interactions reveals folding principles of the human genome

我们经常在HiC文章中可以看到类似衰减曲线(交互频率与距离变化关系),如下图：
于是尝试绘制这个曲线(胆小的我走了些弯路)

image.png

思路：
将交互信息分bin后，计算bin之间距离与交互频率关系。

实践：以HiC-Pro 结果文件作为输入

Step1:输入文件介绍，进行ice 之前的输入文件

image.png

Step2: 计算分bin后，每条染色体bin number 区间，为了后面去掉染色体之间的交互片段做准备

$ cat B1_L7_A001.inter.clean.40000.raw_abs.bed | bedtools groupby -i - -g 1 -c 4 -o min >min.txt
$ cat B1_L7_A001.inter.clean.40000.raw_abs.bed | bedtools groupby -i - -g 1 -c 4 -o max >max.txt
$ paste min.txt max.txt | cut -f 1,2,4 >chr_40kb_bin.txt

image.png

Step3: 去掉raw_matrix中，染色体之间的交互片段

## 进行循环判断chr1..chrY内部交互matrix
cat chr_40kb_bin.txt | while read id;do 
echo $id;
arr=( $id )
chr=${arr[0]};
start=${arr[1]};
end=${arr[2]}
echo $end;

( nohup cat B1_L7_A001.inter.clean.40000.raw.matrix | awk -v start=$start -v end=$end  'BEGIN{FS=OFS="\t"}{if($1>start && $1<end && $2>start &&  $2<end){print $0}}' >test_${chr} & )
done

Step4：计算距离与交互频率的关系

## 合并所有染色体的bin -bin-count 文件（上三角）
$ cat test* | awk 'BEGIN{FS=OFS="\t"}($2>$1){print $0}' > test_merge_bin


## 计算距离与交互频率列表
$ cat test_merge_bin | awk 'BEGIN{FS=OFS="\t"}{print ($2-$1)*40000,$3}' > test_merge_bin2

## 相同距离合并
$ cat test_merge_bin2 | awk 'BEGIN{FS=OFS="\t"}{s[$1]+=$2}END{for (i in s){print i,s[i]}}' | LC_ALL=C sort -k 1,1n > test_merge_bin3

## 除以交互片段之和，计算交互频率
$ sum=`cat test_merge_bin3 | awk '{sum+=$2}END{print sum}'`
$ cat test_merge_bin3 | awk -v Sum=$sum 'BEGIN{FS=OFS="\t"}{print $1,$2/Sum}' > test_merge_bin4

结果文件

image.png

Step5: 画图

library(tidyverse)
library(ggplot2)
data=read_delim("test_merge_bin4",delim="\t",col_names=F)

coefficients( lm(log10(X2)~log10(X1),data = data %>% filter(X1>1e5 & X1 <1e6)))[2]

ggplot(data,aes(x=X1,y=X2)) + #+geom_point()+
  geom_vline(xintercept = 1e5,linetype="dashed")+
  geom_vline(xintercept = 1e6,linetype="dashed")+
  scale_x_continuous(trans= 'log10')+
  scale_y_continuous(trans= 'log10')+
  stat_smooth(method = 'loess', span = 0.1, se = FALSE, level = 0.95)+
  annotate("text",x=1e8,y=1e-2,label="slope : -1.02")+
  labs(x="Distance",y="contact probability")+
  theme_classic()