日常瞎掰
HiC-seq可以用来研究基因组范围内的染色体片段之间的交互情况，这其中会涉及到染色质区室(A/B compartments)的变化情况。HiC的相关概念这里不再多说，不了解的可以查阅之前的帖子[HiC相关概念总结]。今天我们来说说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。哦了，今天就到这里了~~~
往期回顾
双曲线火山图一键拿捏
ChIP-seq数据质控
ChatGPT！见证AI的力量！
ChIPseeker绘图函数借用
R语言书籍免费领
网友评论