HiC | contacts vs distance

日常瞎掰

  在HiC的交互矩阵里，随着基因组距离的增加，染色体片段之间的交互频率降低。通过染色体距离与交互频率间的这样关系，我们可以获知全基因组范围内的pattern。并且，已知在细胞周期、细胞类型等不同情况下会呈现出不同的pattern。因此，通过绘制contacts vs distance的曲线图，我们可以观察不同样本间这种pattern是否有变化。

准备数据

  cooltools结合cooler可以用来快速获取绘图的数据。下面使用软件自带的数据进行演示：

import cooler
import cooltools as ct

cool_file = cooltools.download_data("HFF_MicroC", cache=True, data_dir='./data/')
clr = cooler.Cooler('./test.mcool::/resolutions/100000')
cvd = cooltools.expected_cis(clr=clr, nproc=4)
cvd
     region1 region2  dist  n_valid  count.sum  balanced.sum    count.avg  balanced.avg  balanced.avg.smoothed  balanced.avg.smoothed.agg
0       chr2    chr2     0     2309        NaN           NaN          NaN           NaN                    NaN                        NaN
1       chr2    chr2     1     2290        NaN           NaN          NaN           NaN               0.000773                   0.000837
2       chr2    chr2     2     2280  6796628.0    165.955165  2980.977193      0.072787               0.067210                   0.072680
3       chr2    chr2     3     2271  4343042.0    104.663306  1912.391898      0.046087               0.044613                   0.047050
4       chr2    chr2     4     2269  3035631.0     74.149739  1337.871750      0.032679               0.031651                   0.032792
...      ...     ...   ...      ...        ...           ...          ...           ...                    ...                        ...
3250   chr17   chr17   828        3      146.0      0.006447    48.666667      0.002149               0.000121                   0.000050
3251   chr17   chr17   829        2      133.0      0.006013    66.500000      0.003007               0.000121                   0.000050
3252   chr17   chr17   830        1       55.0      0.002013    55.000000      0.002013               0.000121                   0.000049
3253   chr17   chr17   831        0        0.0      0.000000          NaN           NaN               0.000121                   0.000049
3254   chr17   chr17   832        0        0.0      0.000000          NaN           NaN               0.000122                   0.000049

[3255 rows x 10 columns]

   balanced.avg.smoothed 是用于绘图的数据，该数据是区分不同染色体的，而balanced.avg.smoothed.agg不区分染色体。这两列数据都经过smooth处理，可以绘出平滑的曲线。
  值得注意的是，软件推荐的做法是先去除着丝粒区域，因为该区域pattern与基因组其他区域不同。如果采用此方式需要准备染色体长短臂的位置信息数据框。如果物种是人的话，可以直接用代码获取，其他物种就要自己构建了。

import bioframe
import pandas as pd

hg38_chromsizes = bioframe.fetch_chromsizes('hg38')
hg38_cens = bioframe.fetch_centromeres('hg38')
hg38_arms = bioframe.make_chromarms(hg38_chromsizes,  hg38_cens)
hg38_arms = hg38_arms[hg38_arms.chrom.isin(clr.chromnames)].reset_index(drop=True)
hg38_arms
   chrom     start        end     name
0   chr2         0   93139351   chr2_p
1   chr2  93139351  242193529   chr2_q
2  chr17         0   24714921  chr17_p
3  chr17  24714921   83257441  chr17_q

cvd = ct.expected_cis(clr=clr,view_df=hg38_arms,nproc=4)

绘图

  数据已经拿到，下面就可以绘图了。由于数据跨度范围比较大，绘图的时候一般都是分别对X、Y数据做对数处理。

import numpy as np
import matplotlib.pyplot as plt

cvd['genomic_dist'] = cvd['dist'] * 100000
cvd['balanced.avg.smoothed.agg'].loc[cvd['dist'] < 2] = np.nan

plt.loglog(cvd['genomic_dist'], cvd['balanced.avg.smoothed.agg'])
plt.xlabel('genomic distance')
plt.ylabel('contact frequency')
plt.legend()
plt.show()

结果如下：

  注意，绘图之前还是需要将间隔的bin个数转化为实际的染色体线性距离，如这里用的矩阵分辨率为100kb，所以数据中dist列乘以分辨率得到实际具体。同时，还需要过滤一下数据去除低质量的值。

结束语

  由于交互频率与基因组的线性距有这样的pattern，所以绘制contacts vs distance曲线图用的数据来自于染色体内的交互。cooltools在HiC数据可视化方面功能还是挺多的，而且使用起来也挺方便，不管是命令行还是模块导入，值得推荐。

往期回顾

bed基因注释
scanpy踩坑实录
差异基因密度分布
R绘图配色总结
saddleplot | A/B compartments