一文读懂3D基因组之Hi-C及其数据处理
Hi-C的起源与发展
1953年，沃森和克里克发现DNA双螺旋结构，这标志着生物界进入基因组的研究时代。此后不久，胰岛素的氨基酸序列得以破解，DNA测序逐渐成为生物学家的重要目标。1975年Frederick Sanger等人发明了Sanger测序法，这一突破的性的工作为基因组学的发展甚至是整个生物学发展起到了推动作用。而后人类基因组计划的开启，人类基因草图的完成等一系列的重大科研项目的进行推进了生命科学的发展。
单分子实时测序，离子半导体测序以及454测序等二代测序的发展，使得测序通量下降成本降低，促进了基因组学的发展。
上图主要阐述了线性基因组学主要研究工作，随着研究深入以及科技发展，科学家已经不满足于在二维寻找基因的信息，尤其是发现染色质有较多有趣的现象时，二维基因组并不能很好解释染色质的这些现象，因而人们开始寻找空间结构也能影响基因表达的信息。
近年来，随着线性基因组的发展，研究者利用荧光原位杂交、染色质构象捕获（chromosome conformation capture,3C）等技术，更加深刻了解到细胞核中染色质的三维构象以及内容。之后研究者们将3C技术与基因芯片技术融合产生了4C、5C技术。09年Job Dekker等人将3C技术与二代测序技术结合起来发明Hi-C技术研究染色质DNA与DNA的交互，阮等人发明ChIA-PET技术研究蛋白与DNA交互信息。基因组学的研究也从一维转变为三维。
从一维的基因组信息, 到二维的交互网络, 再到三维的染色质结构. (a) 一维基因组和表观遗传学数据, 以及 ChIA-PET 示例数据, 虚线表示两个 DNA 片段之间有连接; (b) RNAPII 参与的染色质交互网络、相应的双对数曲线图(典型的无尺度网络标志), n 表示染色质交互的数目, f 表示有 n 个交互的节点数目. 在放大的图中, 节点的颜色表示节点来自于不同的染色体; (c) 通过 ChIA-PET 数据重建的一个粗略的染色体在细胞核内的分布结构. 不同的颜色表示不同的染色体. 本图来自于 Sandhu 等人。
HiC的重要应用
01
甲醛交联，固定细胞核内染色质构象
02
酶切，并用生物素标记
03
成环
04
打断
05
钓取带有生物素标记的酶切片段
01
数据过滤，筛选（二代测序通用 ^__^ ）
02
比对；比对的方式主要分两种，一种判断每条reads是否含有酶切位点，有则去掉酶切位点之后的数据在进行bowtie2单端比对；另一种采用单端比对的策略，以25bp为起始长度，每次增加5bp直到该reads比对到基因组具有唯一性。
03
寻找酶切片段；比对寻找到reads pairs在基因组物理位置之后，通过插入片段大小的限制搜索reads pairs两端每条read所对应的最近的酶切片段。酶切片段的位置代表了DNA交互产生的大致位置。
04
筛选fragment pairs；
我们只需要Valid Pairs。
05
Binning；将Valid Pairs的交互信息mapping到 基因组的位置，最终转换成为每两个bin的交互强度。
06
交互矩阵标准化；标准化方法主要分为两类，一类是基于矩阵，进行数学上的标准化，例如迭代等，另一类是基于生物学意义（例如mappingability）上的标准化。
分析代码会在博客（www.zilhua.com）以protocol的形式发表
3D基因组之Hi-C数据分析
（1）需要数据：Hi-C交互矩阵：N * N 的矩阵
由第二步得到的Hi-C数据是这样的格式（中间以tab分割）
（2）软件包：tadbit
https://github.com/zilhua/tadbit
官方网站：http://3dgenomes.github.io/TADbit/
（3）计算代码：
#!/user/bin/python
# -*- coding:UTF-8 -*-'''Created on ${date}@author: www.zilhua.com'''from pytadbit import Chromosome# initiate a chromosome object that will store all Hi-C data and analysis#设置环境,初始化my_chrom = Chromosome(name='chr19', centromere_search=True,species='Homo sapiens', assembly='NCBI36')###my_chrom'''Chromosome chr19:0 experiment loaded:0 alignment loaded:species : Homo sapiensassembly version: NCBI36'''#加载数据以BingRen的数据为例# load Hi-C datamy_chrom.add_experiment('k562', cell_type='wild type', exp_type='Hi-C',identifier='k562',project='TADbit tutorial',hic_data="scripts/sample_data/HIC_k562_chr19_chr19_100000_obs.txt",resolution=100000)my_chrom.add_experiment('gm06690',cell_type='cancer', exp_type='Hi-C',identifier='gm06690',project='TADbit tutorial',hic_data="scripts/sample_data/HIC_gm06690_chr19_chr19_100000_obs.txt",resolution=100000)#参数解释：k562 样本名;cell_type：细胞类型#exp_type:实验类型;hic_data：HiC数据;#resolution:分辨率#对两个样本分别可视化my_chrom.experiments["k562"].view()my_chrom.experiments["gm06690"].view()
#如果这两个样本是对照实验，则可以这样exp = my_chrom.experiments["k562"] + my_chrom.experiments["gm06690"]my_chrom.add_experiment(exp)my_chrom.visualize([('k562', 'gm06690'), 'k562+gm06690'])
#Find Topologically Associating Domainsmy_chrom.find_tad('k562', n_cpus=4)my_chrom.find_tad('gm06690', n_cpus=4)#参数解释：函数find_tad第一个参数k652样本名;#n_cpus线程数，个人电脑一般是1--4,服务器看配置#获取计算结果my_chrom.experiments["k562"].tadsmy_chrom.experiments["gm06690"].tads#将数组转换成txt文本my_chrom.experiments["k562"].write_tad_borders()#结果中，第一列就是TAD的编号，#第二列第三列分别是TAD的起始和终止位置#第四列打分#将TAD的边界在热图上显示my_chrom.visualize('k562', paint_tads=True)
#如果需要对比两个样本的TAD情况my_chrom.visualize([('k562', 'gm06690')],paint_tads=True, focus=(490,620), normalized=True)#参数解释：paint_tads是否描述TAD边界;#focus 显示染色体位置;normalized是否标准化
#TAD densitymy_chrom.tad_density_plot('k562')
#保存和重载数据my_chrom.save_chromosome("some_path.tdb", force=True)
from pytadbit import load_chromosomemy_chrom = load_chromosome("some_path.tdb")print my_chrom.experiments