【Hi-C】交互頻率與距離衰減曲線繪制

參考：

發(fā)現(xiàn)AB區(qū)室 09 science Comprehensive mapping of long range interactions reveals folding principles of the human genome

我們經(jīng)常在HiC文章中可以看到類似衰減曲線(交互頻率與距離變化關(guān)系),如下圖：
于是嘗試?yán)L制這個曲線(膽小的我走了些彎路)

image.png

思路：
將交互信息分bin后季俩，計算bin之間距離與交互頻率關(guān)系犯眠。

實(shí)踐：以HiC-Pro 結(jié)果文件作為輸入

Step1:輸入文件介紹腊尚，進(jìn)行ice 之前的輸入文件

image.png

Step2: 計算分bin后蟀淮，每條染色體bin number 區(qū)間掰盘，為了后面去掉染色體之間的交互片段做準(zhǔn)備

$ 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中牙勘，染色體之間的交互片段

## 進(jìn)行循環(huán)判斷chr1..chrY內(nèi)部交互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：計算距離與交互頻率的關(guān)系

## 合并所有染色體的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

結(jié)果文件

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()

結(jié)果

image.png

思考：

使用raw_matrix 就行散劫，不需要iced_matrix.
不是直接使用bedpe文件，計算距離與概率關(guān)系牍戚。
效果如下：

image.png
需要去掉染色體之間的交互片段
繪圖時候侮繁，需要橫縱軸取對數(shù)
按照09 science 概念，需要計算1e5~1e6 區(qū)間斜率翘魄，判斷分球型與平衡型鼎天。

image.png
部分工具也可以畫出類似結(jié)果 hicplotdistvscounts

image.png

歡迎評論交流~??

最后編輯于：2020.10.07 00:41:37

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者