參考：

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

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

image.png

• 思路：
  將交互信息分bin后，計算bin之間距離與交互頻率關系。
  
  
  
  

實踐：以HiC-Pro 結果文件作為輸入

Step1:輸入文件介紹，進行ice 之前的輸入文件

image.png

image.png

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

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

## 進行循環(huán)判斷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()

• 結果

  
  
  image.png

思考：

• 使用raw_matrix 就行，不需要iced_matrix.

• 不是直接使用bedpe文件，計算距離與概率關系。
  效果如下：

  
  
  image.png
• 需要去掉染色體之間的交互片段
• 繪圖時候，需要橫縱軸取對數(shù)

• 按照09 science 概念，需要計算1e5~1e6 區(qū)間斜率，判斷分球型與平衡型。

  
  
  image.png
• 部分工具也可以畫出類似結果 hicplotdistvscounts
  
  image.png

歡迎評論交流~??