三维基因组学——如何用Hi-C数据分析拓扑关联结构域

首先什么是三维基因组学？三维基因组学(Three-Dimensional Genomics, 3D Genomics)，是指在考虑基因组序列、基因结构及其调控元件的同时，对基因组序列在细胞核内的三维空间结构，及其对基因转录、复制、修复和调控等生物过程中功能的研究。

随着基因组学和表观基因组学的发展，基因组学研究的维度经历了从一维到三维的转变：

• 一维：基因组序列的测定与组装
• 二维：调控序列/转录因子—基因互作网络、染色质状态、表观修饰
• 三维：基因组的三维空间结构

基因组不是简单的线性序列，是有三维空间结构的，而且这种三维空间结构可以对DNA辅助、基因转录调控、染色质浓缩分离等生物学过程产生重要的影响。

我们知道，真核生物的染色质结构由低级到高级可以分为4种：

• 一级结构：一系列核小体相互连接成的念珠状结构

• 二级结构：由核小体连接起来的纤维状结构经螺旋化形成中空的螺线管

• 三级结构：由螺线管螺旋化形成的筒状结构，称为超螺线管

• 四级结构：超螺线管进一步螺旋折叠形成染色单体

三维基因组的层级结构类似于真核生物的染色质结构，在不同的分辨率下也可以划分为4个层级结构：

• 染色质环（Chromatin loops），分辨率：<1 Kb - few Mb：染色质在空间中形成的环状结构，使相距很远的染色质区域在三维空间中可以聚集在一起。

• 拓扑关联结构域（Topologically associating Domains, TAD），分辨率：40 Kb - 3 Mb：相互作用相对频繁的染色质区域。

• 染色质区室（A/B compartments），分辨率：1 - 100 Mb：A compartments：开放的染色质、表达活跃、基因丰富、较高的GC含量、激活型的组蛋白标记，通常位于细胞核内部；B compartments：关闭的染色质、表达不活跃、基因缺乏、结构紧凑、沉默基因的组蛋白标记，通常位于细胞核外围。

• 染色质疆域（Chromosome Territory, CT），分辨率：~100 Mb：在真核生物中，细胞核内染色质分布并不是随机的，为了跨越较大距离实现相互作用，这些染色质会在三维空间中靠的更近，这就是染色质疆域。

不同层次分辨率下的研究重点不同，比如最小的loop层次对应的是基因级别甚至更小的元件互作；10kb级别一般就可以鉴定TAD之间的互作关系（也是研究比较多的）；更大一些的染色质区室也是研究比较多的内容；在染色质疆域这个层次主要就是研究染色体之间的互作关系了。相应的，分辨率越高，需要的有效互作数据量也越大。

目前三维基因组结构的检测时基于染色体构象捕获技术，也就是3C（chromosome conformation capture）技术，3C是所有染色体构象捕获技术的基础。根据染色质的互作类型，3C技术又可以大致分为5种方法：

• 1 versus 1：3C的由来，经过酶切、DNA片段绑定、反向交联后，通过qPCR确认互作位点。
• 1 versus Many/All：4C技术，3C技术的升级版。反向交联前和3C一样，经过了第二轮消化和绑定，通过特定位点的反向PCR检测特定位点与全基因组潜在位点互作情况。
• Many versus Many：5C技术，基于3C的另一升级版。5C技术的通量增大。
• Many versus All：Capture-C技术，最大不同是在互作片段对的捕捉技术上，它是利用生物素标记反向互补到限制酶酶切位点，从而进行对所有感兴趣的互作位点和全基因组位点之间互作对之间的捕捉。
• All versus All

其中“1”“Many”以及“All”代表的是在一次实验中所涉及的位点，比如说 “1 versus All” 的意义就是该次实验调查的是一个位点和全基因组中所有可能潜在互作位点之间的互作情况。

近年来all versus All也就是检测全局的互作技术发展比较快，有以下一些技术应用：

1. Hi-C技术。近年最火的全基因组 3D 基因组测序技术，不仅可以用于检测全基因组的三维基因组结构和染色质互作，同时还可以用于辅助基因组组装（同样用的很多，现在组参考基因组都要测HiC，以组装到染色体水平）等。
2. ChIA-PET技术（ chromatin interaction analysis paired-end tag sequencing）：与3C基础上发展的其他技术不同，区别在于第一步对于互作位点的 DNA 破碎不是通过限制酶进行消化，而是利用超声波击碎。然后应用抗体对特定蛋白参与的互作区段进行富集，并对这些互作区段进行消化连接，提取含有接头的双端序列（ paired end tag， PET）进行互作检测。
3. DLO Hi-C技术（ digestion-ligation-only Hi-C）：该技术相对于传统的全基因组染色体构象捕获技术 Hi-C 而言更加高效简单，仅需要两轮的消化连接过程，无须生物素标记，未连接的 DNA 也可以被有效地去除，极大地提高了染色体构想捕获效率。 DLO Hi-C 技术更像是 Hi-C 技术的一个升级优化。
4. 原位Hi-C (in situ Hi-C)：使用完整的细胞核进行后续反应，减少了Hi-C中随机连接造成的背景噪音。使用四碱基酶Mbol进行酶切，提高了分辨率，实验时间由Hi-C的7天缩短至3天。
5. HiChIP (in situ Hi-C followed by chromatin immunoprecipitation)：该技术是一种利用原位Hi-C原理和转座酶介导构建文库来解析染色质构象的方法。

现在Hi-C测序是应用最广泛的三维基因组学技术，对Hi-C数据的处理，构建Hi-C互作图谱和鉴定TAD结构域是最关键的问题，也发展了一大批专门的生物信息学算法软件，我这里用经典的软件HiC-Pro、HiCPlotter和R包HiTC跑一遍Hi-C数据分析的流程。

为了快速得到结果，这里选择了Gossypium hirsutum四倍体陆地棉的两条染色体参考序列，和一部分Hi-C测序结果（50万条reads的双端测序结果）跑下流程。后续制作Hi-C互作矩阵，和分析Compartment和TAD的.bed后缀文件和.matrix后缀文件来自尤师姐（前面的这点数据做不出互作图，太少了，跑完所有数据又很慢……）。

两个软件安装过程就不说了，可以参考github官方，如果HiC-Pro很难安装依赖的话可以用官方提供的singularity镜像：

nservant/HiC-Pro: HiC-Pro: An optimized and flexible pipeline for Hi-C data processing (github.com)

当前路径文件结构如下：

3.1 获得染色质互作矩阵（HiC-Pro）

1
2
3
4
5
6
7
8

# 1.准备酶切片段
## dpnii是选择的限制性核酸内切酶，可以用别的，查看digest_genome.py源码
python digest_genome.py Gh_genome.fa –r dpnii –o Gh_dpnii.bed
# 2.统计酶切片段大小
awk ’{print $3-$2;}’ Gh_dpnii.bed > Gh_dpnii_size.txt
# 3.构建参考基因组索引，生成基因组大小的文件
samtools faidx Gh_genome.fa
awk ‘{print $1”\t”$2;}’ Gh_genome.fa.fai > Gh_size.txt

这里酶切参考基因组生成的bed文件如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

# 4.修改config文件
## 主要修改基因组索引文件路径、基因组大小文件路径、酶切片段文件路径和分辨率
## BOWTIE2_IDX_PATH、GENOME_SIZE、GENOME_FRAGMENT、BIN_SIZE

vim config.txt


#######################################################################
## Alignment options
#######################################################################

FORMAT = phred33
MIN_MAPQ = 0

BOWTIE2_IDX_PATH = /home/Bioinfor/Gh_index
BOWTIE2_GLOBAL_OPTIONS = --very-sensitive -L 30 --score-min L,-0.6,-0.2 --end-to-end --reorder
BOWTIE2_LOCAL_OPTIONS =  --very-sensitive -L 20 --score-min L,-0.6,-0.2 --end-to-end --reorder

#######################################################################
## Annotation files
#######################################################################

REFERENCE_GENOME = Gh
GENOME_SIZE = /home/Bioinfor/Gh_size.txt

#######################################################################
## Digestion Hi-C
#######################################################################

GENOME_FRAGMENT = /home/Bioinfor/Gh_dpnii.bed
LIGATION_SITE = AAGCTAGCTT
MIN_FRAG_SIZE = 100
MAX_FRAG_SIZE = 160000
MIN_INSERT_SIZE = 200
MAX_INSERT_SIZE = 600

#######################################################################
## Contact Maps
#######################################################################

BIN_SIZE = 5000 20000 100000
MATRIX_FORMAT = upper

主要是修改以上内容，但也要注意下Hi-C双端测序文件的后缀，要让软件能匹配上。

1
2

# 5.运行HiC-Pro，获得染色质互作矩阵
HiC-Pro -c config.txt -i Gh/ -o HiC_Pro_out

HiC-Pro的主要结果放置在目录HiC_Pro_out/ hic_results/matrix/Gh/下，为5000、20000、100000 分辨率下的_abs.bed 以及_iced.matrix后缀文件。其他结果文件和分析图片可以在hic_results文件夹里查看，这里不展示了。

主要展示后续用的文件，以下为Gh_5000_abs.bed文件：

这个文件将染色体划分为5000 bp的bin，并且在第四列进行编号。

以下为Gh_5000_iced.matrix文件：

这个文件第一列和第二列都是bin编号，第三列为两个bin之间归一化之后的互作强度。

3.2 绘制染色质互作图（HiCPlotter）

这里使用的是尤师姐提供的跑完了A01染色体的Gh_20000_iced.matrix与Gh_20000_abs.bed文件（我的示例Hi-C数据量太少）。

1
2
3
4
5
6
7
8
9
10
11

# HiCPlotter绘制染色体互作图
python HiCPlotter.py -f Gh_20000_iced.matrix -o Ghir_A01 -r 20000 -tri 1 -bed Gh_20000_abs.bed -n Ghir_A01 -chr Ghir_A01
## HiCPlotter绘制染色质互作几个常用参数：
-chr 染色体名称，如果染色体名称不是Chr*，-chr参数需传入对应的值
-o 输出文件名
-tri 默认值为0，1代表着传入的matrix和bed文件是HiC-Pro运行的结果
-r 分辨率
-n 任务名

# ZModem协议传输文件
sz Ghir_A01-Ghir_A01.ofBins\(0-5887\).20K.png

染色体Ghir_A01内部的互作图Ghir_A01-Ghir_A01.ofBins(0-5887).20K.png如下所示。

3.3 鉴定TAD（HiCPlotter）

这里使用的数据与3.2一样。

1
2
3
4
5
6
7
8
9
10

# HiCPlotter鉴定TAD
python HiCPlotter.py -f Gh_20000_iced.matrix –bed Gh_20000_abs.bed –o TAD -r 20000 -n Ghir_A01 -chr Ghir_A01 -tri 1 -fh 0 -s 600 -e 900 -ptd 1 -pi 1
## HiCPlotter绘制TAD几个常用参数：
-ptd 默认0，输入1，调用绘制TAD算法
-fh 输入文件中抬头需要删除的行数，没有header lines就是0
-s 起始位点（第几个bin）
-e 结束位点（第几个bin）

# ZModem协议传输文件
sz TAD-Ghir_A01.ofBins\(600-900\).20K.png

染色体Ghir_A01的12Mb到18Mb之间TAD鉴定结果图TAD-Ghir_A01.ofBins(600-900).20K.png如下：

3.4 鉴定染色质区室和TAD（HiTC包）

这里使用的数据来自尤师姐提供的LG02_50000_abs.bed与LG02_50000_iced.matrix文件，也是HiC-Pro跑出来的结果文件。

因为用到了R包，集群R有点问题，我暂时传到本地用自己电脑跑了下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# 下载和加载HiTC包
> BiocManager::install("HiTC")
> setwd("D:/zhuomian/D5_50Kb/D5_50Kb")
> library(HiTC)

# 导入数据（importC函数读取HiC-pro的结果）
> mydata <- importC("D:/zhuomian/D5_50Kb/D5_50Kb/LG02_50000_iced.matrix", xgi.bed = "D:/zhuomian/D5_50Kb/D5_50Kb/LG02_50000_abs.bed")

# 主成分分析和鉴定Compartment
> pc <- pca.hic(mydata$LG02LG02, npc = 1, asGRangesList = TRUE)
> write.csv(pc, file = "D:/zhuomian/D5_50Kb/D5_50Kb/LG02_50k.csv")
> mydata2 <- read.table("D:/zhuomian/D5_50Kb/D5_50Kb/LG03_50k.csv",header = TRUE,sep = ",")
> barplot(mydata2$score)

# 热图展示
> mapC(mydata$LG02LG02, trim.range = 0.94, col.pos = c("white", "red"))

这里主成分分析得到的LG02_50k.csv，统计第一主成分score做条形图就可以鉴定染色质区室Compartment A/B了，以下链接是提出者的文章。

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2858594/

跑出来的互作热图如下：

也可以选取染色体的一部分（一般是自己感兴趣的部分），做热图：

1
2

> data1 <- extractRegion(mydata$LG02LG02,chr = "LG02", from = 1000, to = 5000000)
> mapC(data1, trim.range = 0.94,col.pos = c("white", "red"))

关于如何用HiTC包分析Hi-C数据，也可以看bioconductor的一篇文章：

Analyzing Hi-C data with the HiTC BioC package (bioconductor.org)