Hi-C 数据分析流程¶
一句话概述¶
Hi-C 是一种捕获染色质三维空间相互作用的测序技术，能揭示基因组的三维折叠结构，包括 TAD（拓扑关联域）和染色质环。
核心知识点表格¶
| 知识点 | 说明 |
|---|---|
| Hi-C | 基于邻近连接的染色质构象捕获技术，全基因组尺度 |
| TAD（拓扑关联域） | 染色质自我互作频率高的区域，通常 200kb-2Mb，是基因调控的基本单元 |
| A/B Compartment | A区室=活跃染色质（开放），B区室=沉默染色质（紧缩） |
| 染色质环（Loop） | 两个基因组位点的直接物理接触，通常由 CTCF/Cohesin 介导 |
| 接触矩阵 | Hi-C 数据的核心表示，对称方阵记录任意两个基因组区域的接触频率 |
| 分辨率 | bin 的大小，常用 1kb、5kb、10kb、25kb、50kb、100kb |
| ICE 标准化 | 迭代校正法，消除技术偏差（GC含量、可mappability、片段长度等） |
| HiC-Pro | 经典 Hi-C 数据处理流程，从 FASTQ 到接触矩阵 |
| Juicer | Aiden Lab 开发的一站式 Hi-C 分析流程 |
| cooltools | Open2C 生态的 Python 分析工具，2024年发表于 PLOS Comp Bio |
白话解释原理¶
想象一下： DNA 在细胞核里不是一条直线，而是像毛线团一样折叠缠绕。有些看起来离得很远的基因，在3D空间中其实靠得很近。
Hi-C 做的事情就像是： 1. 先用"胶水"（甲醛）把空间中靠近的 DNA 片段粘在一起 2. 用"剪刀"（限制性内切酶）把 DNA 剪碎 3. 把粘在一起的断头重新"焊接"（连接） 4. 测序，看看哪些原本不相邻的片段被连在了一起 5. 连接频率越高 → 在3D空间中越近
最终得到一个"接触热图"，像一张温度地图，颜色越深表示两个区域在空间中接触越频繁。
各步骤详解¶
第一步：Hi-C 数据预处理（HiC-Pro）¶
白话解释： HiC-Pro 是最经典的 Hi-C 数据处理工具，从原始 FASTQ 文件处理到接触矩阵。
# 1. 安装 HiC-Pro（推荐 conda）
conda create -n hicpro python=3.9 # 创建环境
conda activate hicpro # 激活环境
conda install -c bioconda hicpro # 安装HiC-Pro
# 2. 准备配置文件 config-hicpro.txt
cat > config-hicpro.txt << 'EOF'
# 基因组参考
BOWTIE2_IDX_PATH = /data/genome/bowtie2_index # Bowtie2索引路径
REFERENCE_GENOME = hg38 # 参考基因组名
GENOME_SIZE = /data/genome/chrom_hg38.sizes # 染色体大小文件
# 限制性内切酶
LIGATION_SITE = AAGCTAGCTT # 连接位点序列（MboI）
GENOME_FRAGMENT = /data/genome/MboI_resfrag_hg38.bed # 酶切位点文件
# 比对参数
BOWTIE2_GLOBAL_OPTIONS = --very-sensitive -L 30 # 全局比对参数
BOWTIE2_LOCAL_OPTIONS = --very-sensitive -L 20 # 局部比对参数
# 分辨率
BIN_SIZE = 5000 10000 40000 100000 500000 1000000 # 生成多个分辨率的矩阵
# 标准化
MAX_ITER = 100 # ICE迭代次数
FILTER_LOW_COUNT_PERC = 0.02 # 过滤低计数bin百分比
FILTER_HIGH_COUNT_PERC = 0 # 过滤高计数bin百分比
EPS = 0.1 # 收敛阈值
EOF
# 3. 准备输入数据目录结构
mkdir -p rawdata/sample1 # 创建数据目录
# 将 FASTQ 文件放入: rawdata/sample1/sample1_R1.fastq.gz, sample1_R2.fastq.gz
# 4. 运行 HiC-Pro
HiC-Pro \
-i rawdata/ \ # 输入目录
-o results/ \ # 输出目录
-c config-hicpro.txt \ # 配置文件
-p # 生成运行脚本
# 5. 提交任务（集群环境）
cd results/
bash HiCPro_step1.sh # 第一步：比对和配对
bash HiCPro_step2.sh # 第二步：构建矩阵和标准化
第二步：使用 Juicer 处理（一站式方案）¶
白话解释： Juicer 是 Aiden Lab 的工具，操作更简单，自动完成比对、过滤、矩阵构建和特征识别。
# 1. 安装 Juicer
git clone https://github.com/aidenlab/juicer.git # 克隆Juicer代码
cd juicer/ # 进入目录
# 2. 准备参考文件
# 生成限制性内切酶位点文件
python misc/generate_site_positions.py \
MboI \ # 使用MboI酶
hg38 \ # 基因组版本
/data/genome/hg38.fa # 基因组FASTA文件
# 3. 准备目录结构（Juicer要求固定结构）
mkdir -p juicer_project/fastq # 创建fastq目录
# 将FASTQ文件放入 fastq/ 目录
# 4. 运行 Juicer
bash scripts/juicer.sh \
-d juicer_project/ \ # 项目目录
-g hg38 \ # 基因组
-s MboI \ # 限制性内切酶
-z /data/genome/hg38.fa \ # 基因组FASTA
-p /data/genome/hg38.chrom.sizes \ # 染色体大小
-y restriction_sites/hg38_MboI.txt # 酶切位点文件
# 5. 输出的 .hic 文件包含:
# - 多分辨率接触矩阵（自动生成 1kb 到 1Mb）
# - KR 标准化矩阵
# - 自动调用的 TAD（Arrowhead算法）和 Loop（HiCCUPS算法）
第三步：使用 Open2C 生态分析（Python，推荐新手）¶
白话解释： cooltools 是 Python 生态的 Hi-C 分析工具包，模块化设计，灵活且代码清晰。
# 安装 Open2C 工具链
# pip install cooler cooltools pairtools
import cooler # Hi-C 数据格式处理
import cooltools # Hi-C 分析工具
import numpy as np # 数值计算
import matplotlib.pyplot as plt # 绑图
# 1. 加载 .cool 文件（Hi-C 接触矩阵）
clr = cooler.Cooler("sample.mcool::resolutions/10000") # 加载10kb分辨率
print(f"分辨率: {clr.binsize}") # 打印分辨率
print(f"染色体: {clr.chromnames}") # 打印染色体列表
# 2. 提取特定区域的接触矩阵
matrix = clr.matrix(balance=True).fetch("chr1:50000000-55000000") # 提取chr1一段区域
# balance=True 使用 ICE 标准化后的数据
# 3. 可视化接触热图
plt.figure(figsize=(8, 8)) # 设置图片大小
plt.imshow( # 绘制热图
np.log1p(matrix), # log转换使信号更均匀
cmap="YlOrRd", # 黄-橙-红配色
origin="upper" # 原点在左上角
)
plt.colorbar(label="log(contact)") # 添加颜色条
plt.title("Hi-C Contact Map (chr1)") # 标题
plt.savefig("contact_map.pdf") # 保存图片
# 4. 计算 A/B Compartment（区室分析）
# 使用特征向量法
eigs = cooltools.eigs_cis( # 计算顺式特征向量
clr, # cooler对象
n_eigs=3, # 前3个特征向量
phasing_track=gene_density # 用基因密度定向
)
# eigs[0] > 0: A区室（活跃）
# eigs[0] < 0: B区室（沉默）
# 5. 计算 Insulation Score（TAD边界检测）
insulation = cooltools.insulation( # 计算隔离分数
clr, # cooler对象
window_bp=[100000, 200000, 500000], # 窗口大小列表
verbose=True # 显示进度
)
# insulation score 的局部最小值就是 TAD 边界
# 6. 寻找 TAD 边界
boundaries = insulation[insulation["is_boundary_100000"]].copy() # 提取100kb窗口的边界
print(f"检测到 {len(boundaries)} 个 TAD 边界") # 打印边界数量
第四步：Loop 检测¶
# 使用 Juicer 的 HiCCUPS 检测染色质环
java -jar juicer_tools.jar hiccups \
--cpu \ # 使用CPU（无GPU时）
-r 5000,10000,25000 \ # 检测分辨率
sample.hic \ # 输入 .hic 文件
loops_output/ # 输出目录
# 输出文件: merged_loops.bedpe
# 包含: 环的两个锚点坐标、接触频率、显著性等信息
# 也可以使用 cooltools 检测 Loop
# python 代码:
# dots = cooltools.dots(clr, expected, max_loci_separation=10000000)
第五步：可视化¶
# 使用 HiCExplorer 可视化
# pip install hicexplorer
# 命令行绘制三角热图
# hicPlotMatrix \
# --matrix sample_10kb.cool \ # 输入矩阵
# --region chr1:50000000-55000000 \ # 区域
# --log1p \ # log转换
# --colorMap YlOrRd \ # 颜色方案
# --outFileName contact_triangle.pdf # 输出文件
# 使用 pyGenomeTracks 绘制多轨道图
# 可以同时展示: Hi-C 热图 + TAD + Loop + ChIP-seq 信号 + 基因注释
常见报错与解决¶
| 报错 | 原因 | 解决方案 |
|---|---|---|
Low mapping rate (<30%) | 基因组版本不对或数据质量差 | 确认参考基因组版本，检查 FastQC 质量 |
Low valid pairs ratio | 限制性内切酶位点文件错误 | 确认酶切位点与实验方案一致 |
Memory error | 高分辨率矩阵太大 | 降低分辨率，或使用 cooler 的 out-of-core 模式 |
No TADs detected | 分辨率太低或参数不合适 | 使用 25-50kb 分辨率，调整 insulation window |
Too many/few loops | HiCCUPS 参数需调整 | 调整 FDR 阈值和局部背景窗口大小 |
ICE normalization failed | bin 过滤过于严格 | 放宽低计数过滤阈值 |
速查表¶
┌──────────────────────────────────────────────────────────┐
│ Hi-C 数据分析速查 │
├──────────────────────────────────────────────────────────┤
│ 三大处理流程: │
│ HiC-Pro → allValidPairs → 稀疏矩阵 │
│ Juicer → merged_nodups → .hic 文件 │
│ distiller → .pairs → .cool/.mcool 文件 │
│ │
│ 常用分辨率: │
│ A/B Compartment: 100kb-1Mb │
│ TAD: 10kb-50kb │
│ Loop: 5kb-10kb │
│ │
│ 标准化方法: │
│ ICE（迭代校正）— 最常用 │
│ KR（Knight-Ruiz）— Juicer 默认 │
│ VC（Vanilla Coverage）— 简单列总和归一化 │
│ │
│ 特征检测: │
│ TAD: hicFindTADs / Insulation Score / Arrowhead │
│ Loop: HiCCUPS / Fit-Hi-C / cooltools dots │
│ A/B Compartment: 特征向量分解（PCA第一主成分） │
│ │
│ 可视化: │
│ Juicebox — .hic 文件交互式查看 │
│ HiCExplorer — 命令行批量绘图 │
│ pyGenomeTracks — 多轨道整合图 │
│ HiGlass — 在线交互式浏览 │
└──────────────────────────────────────────────────────────┘
面试高频问题¶
什么是 TAD？它在基因调控中有什么作用？ 答：TAD（拓扑关联域）是基因组中自我互作频率高于与外部互作的结构域，大小通常为 200kb-2Mb。TAD 把基因和其调控元件（如增强子）限制在同一个"邻居"里，使得增强子只能调控同一 TAD 内的基因。TAD 边界通常富集 CTCF 和 Cohesin，边界的破坏可能导致增强子"越界"激活其他基因（如在某些肿瘤中看到的）。
A/B Compartment 和 TAD 的区别？ 答：A/B Compartment 是更大尺度的染色质区室划分（Mb级），A 区室对应活跃转录的开放染色质，B 区室对应沉默的紧缩染色质。TAD 是比 Compartment 更小的结构单元（百kb级），嵌套在 Compartment 内部。检测方法也不同：Compartment 用 PCA/特征向量法，TAD 用 insulation score 或方向性指数。
Hi-C 数据为什么需要标准化？常用什么方法？ 答：Hi-C 原始接触频率受多种技术偏差影响——GC含量、可比对性（mappability）、限制性内切酶位点密度等。ICE（Iterative Correction and Eigenvector decomposition）是最常用的标准化方法，通过迭代使每行/列的总和相等。KR（Knight-Ruiz）是另一种矩阵平衡方法。标准化后才能准确比较不同区域的接触频率。