【HiC】绘制TAD边界强度

背景

经常在HiC文章中看到TAD边界强度变化图，大体上可以反映TAD之间隔绝情况,绝缘系数越低，TAD之间交互越弱。
如下图：

image.png

画图的方法：

• 有点类似ChIP-seq peak profiler图,可以借助deeptools来绘制。
  

  image.png
  image.png

• 同时我们也可以用bedtools进行统计，用R画图。

---

实践：

TAD边界位点

image.png

1.deeptools 绘图

注：deeptools默认binsize为10bp,得到的矩阵会非常很大(数据为G级别)，binsize需要根据自己需要进行修改，一定不要默认参数跑哦。

## Step1: boundary 边界位点
$ cat *boundaries.bed | cut -f 4,5 | grep -v "track" | cut -d "|" -f 3 | tr ":-" "\t" | sort -Vk 1 -k2,2n > sample_boundary.bed

## Step2: insulation score 转换成bigwig格式
$ cat *bedGraph| grep -v "track"  | sort -Vk 1 -k2,2n | grep -v "chrY" | grep -v "NA" > sample_insulationscore.bedgraph

$ bedGraphToBigWig sample_insulationscore.bedgraph hg19.ChrSize.bed sample_insulationscore.bw

## Step3：deeptools画图
computeMatrix reference-point  --referencePoint center --binSize 10000 \
-b 500000 -a 500000 \
-R sample_boundary.bed -S sample_insulationscore.bw \
--skipZeros -o sample_matrix_boundary_500Kb_bin10Kb.gz

plotProfile -m sample_matrix_boundary_500Kb_bin10Kb.gz \
              -out sample_matrix_boundary_500Kb_bin10Kb.pdf \
              --plotTitle "boundary strength"

画图效果如图：

image.png

2.bedtools+R 绘图

注：用bedtools overlap 函数，取bed与bedgraph交集，最终输出结果不是按顺序输出。如下图所示。

image.png

• bedtools

## bedgraph区域保留NA位置
$ cat *bedGraph| grep -v "track"  | sort -Vk 1 -k2,2n | grep -v "chrY" >  bedtools_sample_insulationscore.bedgraph

## bed文件拓展
$ 

$ cat sample_boundary.bed | awk 'BEGIN{FS=OFS="\t"}{$2=$2-480000;$3=$3+480000}{print $0}' > bedtools_sample_boundary.bed 

## 统计matrix
$ cat bedtools_sample_boundary.bed | while read id; do 
echo $id| tr " " "\t" > bedtools_tmp;
intersectBed -a bedtools_tmp -b bedtools_sample_insulationscore.bedgraph -wa -wb|sort -k 6,6n | cut -f 8 | xargs| tr " " "\t" >> bedtools_matrix
done

• R 画图

## R 画图
x <- readline()
setwd(x)

library(ggplot2)
library(patchwork)

boundary_matrix <- read.table("bedtools_matrix2")
dim(boundary_matrix)
vector_test <- apply(boundary_matrix,2 ,function(x){median(x,na.rm = T)})
# plot(1:25,vector_test,type = "l")

dat=as.data.frame(vector_test)
dat$num=c(1:25)
p1<-ggplot(data = dat,aes(x=num,y=vector_test))+
  #geom_line()+
  geom_smooth(se = F)+
  scale_x_continuous(breaks =c(1,12.5,25),
                     labels = c("-500Kb","0","500Kb"))+
  theme_classic()+
  xlab("TAD boundary")+ylab("Mean insulation score")
p1
p2 <- ggplot(data = dat,aes(x=num,y=vector_test))+geom_line()+
  scale_x_continuous(breaks =c(1,12.5,25),
                     labels = c("-500Kb","0","500Kb"))
p1+p2

Left Fig: 显示折线图拟合之后的线图；Right Fig：显示原始折线图

---

思考：

• 画此图需要搞清楚，TAD boundary边界如何取点
• 选取合适的工具及其参数，统计boundary附近的InsulationScore分布
• 可以进一步提高计算InsulationScore的精确度，采用10Kb为binsize进行统计，得到的结果不需要进行拟合直接画图

image.png

欢迎评论留言，交流学习😀~