日常瞎掰
在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
