三维基因组分析TAD¶
一句话概述¶
基于Hi-C数据识别拓扑关联域（TAD）边界（Arrowhead/insulation score方法）和染色质环（Loop，HiCCUPS方法），解析三维基因组结构对基因表达调控的影响。
核心知识点总览¶
| 知识点 | 关键内容 | 重要程度 |
|---|---|---|
| Hi-C数据处理 | 原始数据→接触矩阵→归一化 | ⭐⭐⭐⭐⭐ |
| TAD定义与生物学 | 自相互作用增强的基因组域 | ⭐⭐⭐⭐⭐ |
| TAD边界识别 | Insulation score/Arrowhead/TopDom | ⭐⭐⭐⭐⭐ |
| 染色质环(Loop) | HiCCUPS/HICCUPS检测焦点互作 | ⭐⭐⭐⭐ |
| A/B compartment | PCA识别活性/抑制性染色质区室 | ⭐⭐⭐⭐ |
| 归一化方法 | ICE/KR/VC矩阵归一化 | ⭐⭐⭐ |
| 差异TAD分析 | 条件间TAD边界/强度变化 | ⭐⭐⭐ |
| 与基因表达整合 | TAD边界与基因调控的关联 | ⭐⭐⭐ |
各步骤详解¶
第一步：Hi-C数据与三维基因组基础¶
白话解释： Hi-C是测量基因组三维结构的技术。它通过交联→酶切→连接→测序，捕获在三维空间中互相靠近的基因组区域。结果是一张"接触矩阵"——矩阵中每个格子代表两个基因组区域互相接触的频率。TAD（拓扑关联域）是接触矩阵上的"方块"——域内的互作远多于域间，就像城市中不同社区的人更多在社区内部走动。
技术细节： Hi-C数据的层级结构： - Compartments (A/B)：百万碱基(Mb)尺度，对应活性/抑制性染色质 - TADs：数百千碱基(100kb-1Mb)尺度的自关联域 - Loops：特定位点间的点对点接触（如启动子-增强子环） - Sub-TADs/Insulating neighborhoods：TAD内更精细的结构
# Hi-C接触矩阵示意
#
# pos1 pos2 pos3 pos4 pos5 pos6 pos7 pos8
# pos1 [■■■■] . . . . . . .
# pos2 [■■■■] . . . . . . .
# pos3 [■■■■] . . . . . . .
# pos4 [■■■■] . . . . . . .
# pos5 . . [■■■■] . . . . .
# pos6 . . [■■■■] . . . . .
# pos7 . . . . [■■■■] . . .
# pos8 . . . . [■■■■] . . .
#
# [■■■■] = TAD (高频互作域)
# TAD边界 = 域之间互作急剧下降的位置
第二步：Hi-C数据处理流程¶
白话解释： Hi-C原始数据是paired-end reads，每条read对代表一次"接触"。处理流程包括：比对到参考基因组→过滤无效reads→构建接触矩阵→归一化校正系统偏差。
技术细节：
# === Hi-C数据处理（HiC-Pro）===
# HiC-Pro是最常用的Hi-C处理流水线
# 1. 配置文件
cat > config.txt << 'EOF'
BOWTIE2_IDX_PATH = /path/to/bowtie2_index
REFERENCE_GENOME = hg38
GENOME_SIZE = chrom_hg38.sizes
GENOME_FRAGMENT = HindIII_resfrag_hg38.bed
LIGATION_SITE = AAGCTAGCTT
MIN_FRAG_SIZE = 100
MAX_FRAG_SIZE = 100000
MIN_INSERT_SIZE = 100
MAX_INSERT_SIZE = 600
BIN_SIZE = 5000 10000 25000 40000 100000 500000 1000000
MATRIX_FORMAT = upper
N_CPU = 16
EOF
# 2. 运行HiC-Pro
HiC-Pro -i rawdata/ -o hicpro_output/ -c config.txt
# === 或使用Juicer ===
# Juicer是4DN/Aiden Lab开发的流水线
juicer.sh \
-g hg38 \
-s HindIII \
-d /path/to/fastq/ \
-D /path/to/juicer/ \
-t 16
# 输出：inter_30.hic文件（.hic格式，用于Juicebox可视化）
# === 或使用cooler/hicstuff（Python生态）===
# hicstuff处理
hicstuff pipeline \
-g genome.fa \
-e DpnII \
--threads 16 \
-o hicstuff_output \
R1.fq.gz R2.fq.gz
# 输出：.cool格式（用于cooler/cooltools）
# === 接触矩阵归一化 ===
# 使用cooler进行ICE归一化
cooler balance sample.cool
# 或KR归一化
cooler balance --cis-only --max-iters 200 sample.cool
# 使用Juicer Tools进行KR归一化
java -jar juicer_tools.jar dump observed KR \
sample.hic chr1 chr1 BP 10000 chr1_10kb_KR.txt
# HiC-Pro的ICE归一化
python iced/scripts/ice --filter_low_counts_perc 0.02 \
--filter_high_counts_perc 0 \
--max_iter 100 \
--output_bias bias.txt \
raw_matrix.txt
第三步：TAD边界识别¶
白话解释： TAD边界就是两个相邻TAD之间的"墙"。在这里，互作频率急剧下降。识别TAD边界的方法有多种，最常用的是Insulation Score（绝缘分数）和Arrowhead算法。
技术细节：
# === Insulation Score方法（cooltools）===
import cooler
import cooltools
import numpy as np
import matplotlib.pyplot as plt
# 加载cool文件
clr = cooler.Cooler("sample.mcool::/resolutions/25000")
# 25kb分辨率适合TAD分析
# 计算insulation score
insulation_table = cooltools.insulation(clr, window_bp=[100000, 200000, 500000])
# window_bp: 绝缘窗口大小，通常用多个值
# 识别TAD边界（insulation score局部极小值）
# cooltools自动在insulation结果中标注边界
# 边界 = insulation score曲线的局部极小值点
# boundary_strength = 该极小值的拓扑突出度(topographic prominence)
# 查看insulation score
chr1_ins = insulation_table[insulation_table['chrom'] == 'chr1']
plt.figure(figsize=(15, 4))
plt.plot(chr1_ins['start'] / 1e6, chr1_ins['log2_insulation_score_200000'],
linewidth=0.5)
plt.xlabel("Position (Mb)")
plt.ylabel("Insulation Score")
plt.title("chr1 Insulation Score (200kb window)")
# 标记边界
boundaries = chr1_ins[chr1_ins['is_boundary_200000'] == True]
for _, b in boundaries.iterrows():
plt.axvline(x=b['start']/1e6, color='red', alpha=0.3, linewidth=0.5)
plt.savefig("insulation_score_chr1.pdf")
# === Arrowhead算法（Juicer Tools）===
# Arrowhead基于角标变换检测TAD
java -jar juicer_tools.jar arrowhead \
-k KR \
-r 10000 \
sample.hic \
arrowhead_output
# 输出：arrowhead_output/10000_blocks.bedpe
# 格式：chr1 start1 end1 chr2 start2 end2 score ...
# === TopDom（R包）===
# TopDom也是常用的TAD识别方法
# TopDom TAD识别
library(TopDom)
# 准备输入矩阵（n×n归一化接触矩阵）
# 格式：3列bed + n列矩阵值
# TopDom需要的输入格式为domain文件
# 运行TopDom
result <- TopDom(matrix.file = "chr1_40kb_matrix.txt",
window.size = 5) # 窗口大小
# 结果包含：domain（TAD域）和boundary（边界）
tads <- result$domain
boundaries <- result$domain[result$domain$tag == "boundary", ]
cat("Number of TADs:", sum(result$domain$tag == "domain"), "\n")
cat("Number of boundaries:", nrow(boundaries), "\n")
cat("Mean TAD size:", mean(tads$size[tads$tag == "domain"]) / 1000, "kb\n")
第四步：染色质环（Loop）检测——HiCCUPS¶
白话解释： 染色质环是接触矩阵中的"亮点"——两个特定位置之间有异常高频率的互作。最典型的环是CTCF介导的环和启动子-增强子环。HiCCUPS是检测这些环的标准方法。
技术细节： HiCCUPS对接触矩阵中每个像素，比较其与周围背景（donut/bottom-left/horizontal/vertical四种背景滤波器）的信号强度，用泊松检验识别显著富集的焦点互作。四种背景分别排除不同类型的伪信号：donut排除中心像素本身的贡献，bottom-left排除TAD角落效应，horizontal/vertical排除条纹状信号。
# === HiCCUPS Loop Detection ===
# 使用Juicer Tools的HiCCUPS
java -jar juicer_tools.jar hiccups \
-k KR \
-r 5000,10000,25000 \
--threads 16 \
sample.hic \
hiccups_output
# 输出：merged_loops.bedpe
# 格式：chr1 x1 x2 chr2 y1 y2 color observed_bl expected_bl ...
# 查看loop数量
wc -l hiccups_output/merged_loops.bedpe
# === 使用cooltools检测loops ===
# cooltools dots (loops) detection
import cooltools
# 计算expected（用于背景估计）
expected = cooltools.expected_cis(clr, nproc=8)
# 检测loops（dots）
# cooltools.dots返回的DataFrame中，显著性通过各kernel的q值列判断
# 列名格式为 la_exp.<kernel>.qval（如 la_exp.donut.qval）
dots_df = cooltools.dots(
clr,
expected,
max_loci_separation=10_000_000, # 最大搜索距离10Mb
nproc=8
)
# 过滤显著loops（基于donut kernel的q值）
# 注意：返回的已经是通过统计阈值筛选的候选点
print(f"Candidate dots: {len(dots_df)}")
第五步：A/B Compartment分析¶
白话解释： 在更大的尺度上（兆碱基级别），基因组分为A区室（活跃转录、开放染色质）和B区室（基因沉默、紧密染色质）。通过对归一化接触矩阵做PCA，第一主成分（PC1）的正负值就对应A/B区室。
技术细节：
# === A/B Compartment分析（cooltools）===
import cooltools
# 计算compartment（需要基因密度/GC含量作为参考信号来确定A/B方向）
# 方法：对O/E矩阵做PCA
# 注意：cooltools Python API中没有cooltools.genome.fetch_chromsizes_and_gc函数
# 需要使用bioframe库获取GC含量作为phasing track
import bioframe
# 获取染色体大小和bins
bins = clr.bins()[:]
chromsizes = bioframe.fetch_chromsizes('hg38')
centromeres = bioframe.fetch_centromeres('hg38')
arms = bioframe.make_chromarms(chromsizes, centromeres)
# 计算GC含量（需要下载参考基因组fasta）
# genome = bioframe.load_fasta('/path/to/hg38.fa')
# gc_cov = bioframe.frac_gc(bins[['chrom', 'start', 'end']], genome)
# 这里假设已有gc_cov DataFrame，包含chrom, start, end, GC列
gc_cov = bins[['chrom', 'start', 'end']].copy()
# gc_cov['GC'] = bioframe.frac_gc(gc_cov, genome) # 实际使用时取消注释
# 计算eigenvectors
# eigs_cis的phasing_track参数接受一个DataFrame（含chrom, start, end及数值列）
# 它用该数值列与eigenvector做相关来确定A/B方向
eig_values, eig_vectors = cooltools.eigs_cis(
clr,
phasing_track=gc_cov, # 参数名是phasing_track，不是位置参数+phasing_track_col
n_eigs=3
)
# PC1 > 0: A compartment (active)
# PC1 < 0: B compartment (repressed)
# 可视化
import matplotlib.pyplot as plt
chr1_eig = eig_vectors[eig_vectors['chrom'] == 'chr1']
plt.figure(figsize=(15, 3))
colors = ['red' if x > 0 else 'blue' for x in chr1_eig['E1']]
plt.bar(chr1_eig['start']/1e6, chr1_eig['E1'], color=colors, width=0.1)
plt.xlabel("Position (Mb)")
plt.ylabel("PC1 (Eigenvector)")
plt.title("A/B Compartments - chr1")
plt.savefig("compartments_chr1.pdf")
第六步：差异TAD分析与条件比较¶
白话解释： 比较不同条件（如正常vs肿瘤、分化前后）的TAD结构变化：哪些TAD边界消失了（边界弱化）？哪些新出现了？TAD结构的改变可能导致增强子-启动子互作异常，从而改变基因表达。
技术细节：
# === 差异TAD分析 ===
# 1. 比较两个条件的insulation score
import pandas as pd
import numpy as np
ins_ctrl = pd.read_csv("control_insulation.tsv", sep="\t")
ins_tumor = pd.read_csv("tumor_insulation.tsv", sep="\t")
# 合并
merged = ins_ctrl.merge(ins_tumor, on=['chrom', 'start', 'end'],
suffixes=('_ctrl', '_tumor'))
# 计算insulation score差异
merged['delta_insulation'] = (merged['log2_insulation_score_200000_tumor'] -
merged['log2_insulation_score_200000_ctrl'])
# 边界变化
# 边界消失（ctrl有边界，tumor无）
lost_boundaries = merged[(merged['is_boundary_200000_ctrl'] == True) &
(merged['is_boundary_200000_tumor'] == False)]
# 新出现边界
gained_boundaries = merged[(merged['is_boundary_200000_ctrl'] == False) &
(merged['is_boundary_200000_tumor'] == True)]
print(f"Lost boundaries: {len(lost_boundaries)}")
print(f"Gained boundaries: {len(gained_boundaries)}")
# 2. TAD边界处的基因表达变化
# 边界消失可能导致原本分隔的增强子-基因互作增强
# === TAD边界与CTCF/cohesin的关联 ===
library(GenomicRanges)
# 加载TAD边界
boundaries <- import("tad_boundaries.bed")
# 加载CTCF ChIP-seq peaks
ctcf_peaks <- import("CTCF_peaks.narrowPeak")
# 计算TAD边界与CTCF的重叠
overlap <- findOverlaps(boundaries, ctcf_peaks)
pct_with_ctcf <- length(unique(queryHits(overlap))) / length(boundaries) * 100
cat(sprintf("%.1f%% TAD boundaries overlap with CTCF peaks\n", pct_with_ctcf))
# 通常>70%的TAD边界有CTCF结合
第七步：Hi-C数据与基因调控整合¶
白话解释： 三维基因组结构最终服务于基因表达调控。将TAD/Loop信息与RNA-seq、ChIP-seq数据整合，可以理解：哪些增强子通过Loop与哪些基因互作？TAD边界如何限制增强子的活动范围？
技术细节：
# === 整合分析 ===
# 1. Loop锚点（anchor）处的功能注释
# 检查loop两端是否分别对应增强子和启动子
loops <- read.table("merged_loops.bedpe", header = TRUE)
promoters <- import("promoters_2kb.bed")
enhancers <- import("H3K27ac_enhancers.bed")
# 找增强子-启动子Loop
loop_anchor1 <- GRanges(loops$chr1, IRanges(loops$x1, loops$x2))
loop_anchor2 <- GRanges(loops$chr2, IRanges(loops$y1, loops$y2))
ep_loops <- loops[
(countOverlaps(loop_anchor1, enhancers) > 0 &
countOverlaps(loop_anchor2, promoters) > 0) |
(countOverlaps(loop_anchor1, promoters) > 0 &
countOverlaps(loop_anchor2, enhancers) > 0),
]
cat(sprintf("Enhancer-Promoter loops: %d / %d total\n", nrow(ep_loops), nrow(loops)))
# 2. TAD内基因共调控
# 同一TAD内的基因表达相关性应高于TAD间
# 计算intra-TAD vs inter-TAD基因对的表达相关性
实战命令速查¶
# Hi-C处理快速流程
# HiC-Pro
HiC-Pro -i raw/ -o output/ -c config.txt
# Juicer
juicer.sh -g hg38 -s HindIII -d fastq/ -t 16
# TAD检测
java -jar juicer_tools.jar arrowhead -k KR sample.hic output/
# Loop检测
java -jar juicer_tools.jar hiccups -k KR sample.hic output/
# Compartment
cooltools eigs-cis sample.cool --output compartments.tsv
面试常问点¶
Q1: TAD的生物学意义是什么？¶
A: TAD是基因表达调控的功能单元——TAD内的增强子主要调控TAD内的基因，TAD边界限制了调控元件的作用范围。CTCF和cohesin是TAD边界的主要结构蛋白。TAD边界在细胞类型间高度保守（~80%共享），但TAD内部的环路互作是细胞类型特异的。TAD边界破坏可导致增强子错误激活邻近基因（enhancer hijacking）。
Q2: Insulation score和Arrowhead方法的原理区别？¶
A: Insulation score计算每个位点周围方形窗口内的平均互作强度——TAD内部的insulation高，边界处急剧下降形成"山谷"。Arrowhead基于角标变换（corner score）——在接触矩阵中检测"箭头"形的互作模式，直接输出TAD域而非边界点。前者是连续评分可调阈值，后者是离散的域检测。
Q3: Hi-C数据的分辨率取决于什么？¶
A: 分辨率取决于测序深度（有效reads数）。Rao et al. (2014)的定义：分辨率是使≥80%的bin至少有1000条reads的最小bin大小。粗略估算：所需reads数 ≈ (基因组大小/bin大小) × 1000。人类基因组（~3Gb）：1kb分辨率需要~30亿reads，10kb需要~3亿，100kb需要~3000万。TAD分析通常需要25-100kb分辨率（~3000万-1亿reads），Loop检测需要5-10kb分辨率（>3亿reads）。注意这是理想均匀分布的估计，实际需要更多reads。
Q4: Hi-C数据归一化为什么重要？哪些偏差需要校正？¶
A: Hi-C原始信号受多种系统偏差影响：(1) 酶切位点分布不均匀；(2) GC含量影响PCR效率；(3) 可比对性（mappability）差异；(4) 片段长度偏差。ICE（Iterative Correction and Eigenvector decomposition）和KR（Knight-Ruiz）归一化通过矩阵平衡消除这些偏差。不归一化的矩阵会产生大量虚假的高/低互作区域。
Q5: TAD边界在疾病中如何被破坏？¶
A: 常见机制：(1) CTCF结合位点的点突变/甲基化改变导致边界蛋白无法结合；(2) 结构变异（缺失/倒位/易位）破坏TAD边界物理结构；(3) cohesin突变影响环路挤出机制。后果是"enhancer hijacking"——原本被边界隔离的增强子越界激活邻近TAD的原癌基因。典型例子：T-ALL中的TAL1 enhancer hijacking。
易错点¶
1. 分辨率选择不当¶
用100kb分辨率检测Loop（需要5-10kb）或用5kb分辨率分析Compartment（需要100kb-1Mb），都是不恰当的。不同层级的结构需要匹配的分辨率。
2. 忽略Hi-C数据质量指标¶
有效配对比例（valid pairs ratio）过低(<20%)说明建库质量差。长程互作比例（>20kb）过低说明连接效率差。这些质量问题会严重影响下游分析。
3. 直接比较不同深度样本的接触频率¶
不同测序深度的样本即使归一化后，接触频率的绝对值也不可直接比较。条件比较应使用标准化方法（如计算观测/期望比O/E）或专用差异分析工具。
4. 混淆TAD和compartment¶
TAD（~100kb-1Mb）和compartment（~1-10Mb）是不同尺度的结构。Compartment switching（A→B）不等于TAD边界变化。两者的识别方法和分辨率需求完全不同。
5. 过度依赖单一TAD caller¶
不同TAD识别方法对参数和分辨率敏感，结果可能不完全一致。建议用2-3种方法（如insulation score + Arrowhead + TopDom）取共识边界。
补充知识¶
Hi-C技术变体¶
- Micro-C：使用MNase代替限制酶，核小体分辨率
- HiChIP/PLAC-seq：ChIP + Hi-C，特异性检测特定蛋白介导的互作
- Capture Hi-C：靶向捕获特定区域的互作
- scHi-C：单细胞Hi-C，揭示细胞间异质性
分析工具生态¶
- Juicer/Juicebox：Aiden Lab全套工具（处理+分析+可视化）
- cooler/cooltools：Python生态，灵活高效
- HiC-Pro：处理流水线
- HiGlass：交互式Hi-C可视化
- FAN-C：全面的Hi-C分析框架
引用推荐¶
- TAD discovery: Dixon et al., Nature, 2012
- Arrowhead: Rao et al., Cell, 2014
- Insulation score: Crane et al., Nature, 2015
- cooltools: Abdennur et al., PLoS Computational Biology, 2024