HiC-seq
可以用来研究基因组范围内的染色体片段之间的交互情况,这其中会涉及到染色质区室(A/B compartments
)的变化情况。HiC的相关概念这里不再多说,不了解的可以查阅之前的帖子[HiC相关概念总结:https://www.jianshu.com/p/5129bddde3ea]。今天我们来说说saddleplot
的事,一个可以展示染色质区室开放及交互情况的热图。先一睹为快,如下图。
saddleplot
主体本身就是一个普通的热图,没有什么特别的地方,只要得到了绘图矩阵用绘图就很简单了。其实,关键就在于绘图数据的获得,想要获得数据得先知道应该用哪些数据以及如何计算。不过,好在这些耗费精力的事情已经有人做了,咱们有时候更多的是站在巨人的肩膀上,使用现成的软件探索数据本身的意义。当然,毕竟现在厉害的人很多,人家解决问题的方式可能是造一个可以解决问题的工具,所以,可以用来绘制saddleplot
软件也有不少。工具不在于多,能解决问题就行,咱们今天就来说说如何用cooltools
作图。这个工具是用python
编写,可以在命令行直接调用,或者在python
里面导入使用。
https://osf.io/3h9js/download
,如果用window浏览器下载会得到名为test.mcool
的文件,而在linux里面用wget
下载最好重新命名一下,否则文件名就叫download
。这个文件格式是mcool
,里面包含1k
、10k
、100k
、1M
四个分辨率的矩阵数据。注意,使用的时候需要指定使用哪个分别率的数据,例如下面使用了10k
分别率的数据。# 下载测试数据
wget --no-check-certificate -O test.mcool https://osf.io/3h9js/download
# 绘图过程
cooltools eigs-cis -o test.gc.eigs data/test.mcool::resolutions/100000
cooltools expected-cis -p 5 -o test.obs_exp.tsv data/test.mcool::resolutions/100000
cooltools saddle --contact-type cis --strength --out-prefix test.saddleplot -n 38 --qrange 0.025 0.975 --fig pdf data/test.mcool::resolutions/100000 test.gc.eigs.cis.vecs.tsv test.obs_exp.tsv
结果如下:
使用命令行绘图还是挺省事的,三行命令即可得到图。不过,不出意外的话意外就会出现了。仔细一看,绘出图中主体热图部分没有问题,可是上面和左边两个直方图好像有问题呀。直方图是体现A/B compartments
各自成分内的交互频率,也就是AA
、BB
交互的频率。可是,明明从热图上看差异还是挺明显的,怎么直方图却没有波澜呢?原因在哪里呢?经过本人的各种检查数据,最后得出一个结论:这该不会是软件自身的BUG
吧!
从最后一个绘图命令可知,绘图需要三个输入文件:test.mcool
、test.gc.eigs.cis.vecs.tsv
、test.obs_exp.tsv
。不过,确切点说,原始输入只有test.mcool
文件,后面两个都是由其计算得到。
重点来了,其实在使用cooltools eigs-cis
命令得到test.gc.eigs.cis.vecs.tsv
时,可以提供一个额外的相位文件通过参数--phasing-track
指定。该文件是用来矫正A/B compartments
信息的,为什么要矫正呢?这是A/B compartments
信息是由交互矩阵经过PCA降维得来,具体哪些是A哪些是B,需要明确一下,一般认为值>0为A成分,值<0为B成分。所以,为了数据的准确性,最好提供一下相位文件,相位信息可以来源于基因组的GC含量,或者真实的ChIP-seq
数据。相位文件格式如下:
chrom start end GC
chr2 0 100000 0.4358666666666667
chr2 100000 200000 0.40953
chr2 200000 300000 0.42189
chr2 300000 400000 0.43187
python
导入cooltools
包方式来绘图。采用这种方式的好处就是可以更好的控制输出,如果需要可以做一些自定义的修改。import pandas as pd
import numpy as np
import cooltools
import cooler
import bioframe
import warnings
from cytoolz import merge
#绘图函数定义,函数代码过多,可以直接到官网拷贝:
#https://cooltools.readthedocs.io/en/latest/notebooks/compartments_and_saddles.html
def saddleplot()
# 读取交互矩阵
clr = cooler.Cooler('test.mcool::resolutions/100000')
# 根据基因组文件获取相位信息
bins = clr.bins()[:]
hg38_genome = bioframe.load_fasta('./hg38.fa')
gc_cov.to_csv('hg38_gc_cov_100kb.tsv',index=False,sep='\t')
gc_cov = bioframe.frac_gc(bins[['chrom', 'start', 'end']], hg38_genome)
# PCA特征向量
view_df = pd.DataFrame({'chrom': clr.chromnames, 'start': 0, 'end': clr.chromsizes.values, 'name': clr.chromnames})
cis_eigs = cooltools.eigs_cis(clr, gc_cov, view_df=view_df, n_eigs=3)
eigenvector_track = cis_eigs[1][['chrom','start','end','E1']]
cvd = cooltools.expected_cis(clr=clr,view_df=view_df)
# 得到绘图的数据
Q_LO = 0.025
Q_HI = 0.975
N_GROUPS = 38
interaction_sum, interaction_count = cooltools.saddle(clr, cvd, eigenvector_track, 'cis', n_bins=N_GROUPS, qrange=(Q_LO,Q_HI), view_df=view_df)
#绘图
saddleplot(eigenvector_track, interaction_sum/interaction_count, N_GROUPS, qrange=(Q_LO,Q_HI), cbar_kws={'label':'average observed/expected contact frequency'})
结果如下:
使用代码绘图也不是很麻烦,而且还方便控制输出结果。不过,最重要的是这种方式不会出现上面的BUG
,图里面的直方图是正确的。
saddleplot
作为展示A/B compartments
的图,可以很轻松的知道AA、BB、AB、BA四种类型的交互情况,即交互数量和交互强度的估计,可谓一眼就能知道染色质的开放情况。如果多个条件的样本放在一起比较,便可以知道不同条件对染色质开放的影响。一般来说,热图里面从左到右为B->A(从上到下为B->A),不同软件可能会有所不同,对于这点需要留意以免搞反了。顺便提一下,saddleplot
还有cis
和trans
的区别,这里展示的是cis
结果。最后,附上一篇讲诉saddleplot
计算原理的帖子,方便想要学习的同学:https://blog.sciencenet.cn/home.php?mod=space&uid=2970729&do=blog&id=1367103。哦了,今天就到这里了~~~
--------- The End ---------
基因组 | 转录组 | 表观组 | 单细胞
寻找灵魂的工具人
长按扫码加关注