HiC-HiC分析实战：从互作矩阵构建到顺式反式显著性检验

简介这份资源围绕HiC-HiC分析技术展开面向从事三维基因组、染色质构象研究的科研人员与生物信息学学习者帮助其理解染色体间与染色体内相互作用的推断流程。内容涉及HiC技术原理、实验流程、交互矩阵构建与正则化、结构重构及功能注释等关键环节并延伸至TAD结构、基因调控与疾病机制等应用场景。资源包共11个文件以R脚本为主辅以txt说明文档、md说明与conf配置文件压缩包约29KB涵盖HiCUP配置、GOTHiC分析、交互调用与BED12转换等模块便于按流程查阅与复用。目前已有429人学习下载。对于希望快速上手HiC数据分析、理解从原始读段到染色体三维结构重建完整链路的研究者这份资料提供了可参考的脚本框架与操作说明有助于搭建分析流程并排查常见问题。1. HiC-HiC 分析到底在算什么从一次“互作矩阵对不上”的排查说起HiC-HiC 分析以推断染色体间和染色体内的相互作用这句话第一次看容易误解成“把两个 Hi-C 样本放一起比差异”。我最初也这么以为直到某次帮某实验室复现一套染色质构象流程两个生物学重复的互作矩阵在染色体内区域怎么都对不齐才意识到它真正处理的是同一套 Hi-C 数据内部、两两基因组区间之间的互作关系重点是区分顺式intra-chromosomal同一条染色体内部和反式inter-chromosomal不同染色体之间两类接触。它解决的问题很具体当你手上有一份比对后的 Hi-C 接触数据怎么把它整理成可比较、可统计、可推断的互作矩阵并判断哪些染色体间或染色体内的接触显著偏离背景。适合已经跑过比对、拿到.bam或.pairs文件却卡在“矩阵怎么建、背景怎么定、显著性怎么算”这一步的人。下面按我实际落地的顺序拆开讲参数和坑都会给到。2. 从比对结果到互作矩阵HiC-HiC 分析的数据准备与分箱策略2.1 输入文件到底要什么格式为什么.pairs比.bam更省心HiC-HiC 分析的第一步不是写统计模型而是把接触对整理成规整的表格。常见输入有两类比对后的.bam带 Hi-C 特有的配对信息和 4DN 社区推动的.pairs格式。.pairs本质是制表符分隔的文本每行一对接触列依次是 readID、chr1、pos1、chr2、pos2、strand1、strand2。我一般优先用.pairs原因是它去掉了比对细节后续分箱、过滤、算背景都不用再碰比对软件的黑匣子复现性更好。如果你手上只有.bam可以用pairtools转。下面这段是我常用的转换加过滤命令注意--min-mapq和--max-molecule-size这两个参数直接决定后面矩阵的干净程度。# 从 bam 解析出接触对输出 pairs 格式 pairtools parse \ --chroms-path chrom.sizes \ --min-mapq 30 \ --max-molecule-size 1500 \ --walks-policy 5unique \ --output parsed.pairsam \ input.bam # 去掉未比对、多比对、自连等无效接触 pairtools select \ (pair_type UU) or (pair_type UR) or (pair_type RU) \ --output filtered.pairs \ parsed.pairsam # 按基因组坐标排序为分箱做准备 pairtools sort --nproc 8 --output sorted.pairs filtered.pairs逻辑说明parse把每条 read 的配对信息还原成接触对--min-mapq 30过滤低质量比对Hi-C 里低 MAPQ 往往来自重复序列留着会虚高背景--max-molecule-size 1500限制分子跨度避免把远距离随机碰撞当成真实接触select只保留双端唯一比对UU或一端唯一一端多重UR/RU这是多数分析的标准做法。参数不是死的如果你做的是低深度或单细胞 Hi-C--min-mapq可以降到 10但要在报告里写清楚。2.2 分箱大小怎么选5kb、10kb、25kb 的取舍分箱binning是把连续坐标离散成固定窗口窗口内所有接触对归并成一个矩阵元素。HiC-HiC 分析里 bin 大小直接决定矩阵分辨率和稀疏程度。我一般按测序深度和数据用途选bin 大小适用深度适合的分析注意点5 kb 2 亿接触对启动子-增强子级别稀疏需大量平滑10 kb1–2 亿接触对拓扑结构域边界平衡分辨率与覆盖25 kb 1 亿接触对染色体间互作、区室背景稳定推荐入门100 kb任意全基因组区室丢失局部细节选完 bin 后用cooler建矩阵。下面命令把.pairs转成.cool文件--count-as-float在需要做归一化时有用。# 按 25kb 分箱生成 cool 矩阵 cooler cload pairs \ -c1 2 -p1 3 -c2 4 -p2 5 \ chrom.sizes:25000 \ sorted.pairs \ matrix_25kb.cool # 查看矩阵基本信息 cooler info matrix_25kb.cool-c1 -p1 -c2 -p2指定.pairs里染色体和坐标列的位置不同工具输出的列顺序可能不同这一步错了矩阵会全乱。建完先cooler info看 bin 数和非零元素比例如果非零比例低于 5%说明深度不够或 bin 太小需要合并 bin 或补测序。2.3 归一化为什么不能直接用原始计数比较原始接触计数受 GC 含量、酶切位点分布、可比对性影响直接比较两个区域或两个样本会得出错误结论。HiC-HiC 分析里常用 ICE迭代校正或 KR 归一化。cooler balance是 ICE 的常用实现# 对矩阵做 ICE 归一化生成权重 cooler balance --nproc 8 --cis-only matrix_25kb.cool # 导出归一化后的矩阵用于后续统计 cooler dump --balanced -t 8 matrix_25kb.cool balanced_matrix.tsv--cis-only表示只用染色体内接触估计权重再应用到全矩阵。如果你的分析重点是染色体间互作可以不加这个参数但要注意反式接触本身很稀疏权重估计容易不稳。归一化后矩阵元素才是可比较的“相对接触频率”。3. 区分顺式与反式HiC-HiC 分析里最容易被忽略的统计背景3.1 顺式接触的衰减曲线怎么拟合顺式接触随基因组距离增加而衰减这个衰减曲线是判断实验是否正常的核心指标。HiC-HiC 分析里常用cooltools的expected_cis计算每个距离上的平均接触再拟合幂律。下面代码计算并画出一条染色体的衰减曲线import cooltools import cooler import matplotlib.pyplot as plt # 读取归一化后的 cool 文件 clr cooler.Cooler(matrix_25kb.cool) # 计算顺式期望接触按距离分箱 expected cooltools.expected_cis( clr, nproc8, chunksize1000000, ignore_diags2 # 忽略对角线附近 2 个 bin避免自连和近邻噪声 ) # 取 1 号染色体的衰减曲线 chrom chr1 cis_vals expected[expected[chrom1] chrom] plt.loglog(cis_vals[dist], cis_vals[balanced.avg], labelchrom) plt.xlabel(genomic distance (bp)) plt.ylabel(average contact) plt.legend() plt.savefig(cis_decay.png, dpi150)ignore_diags2是关键参数对角线附近包含大量自连和未消化完全片段不忽略会让曲线头部虚高。正常衰减曲线在 log-log 坐标下应接近直线斜率通常在 -1 到 -1.5 之间。如果斜率明显偏平说明背景噪声大或归一化有问题。3.2 反式接触的背景模型为什么不能简单用全基因组平均反式接触不同染色体之间没有距离衰减概念它的背景取决于两条染色体的可及性和测序深度。HiC-HiC 分析里常用两种背景一是全基因组平均反式接触二是按染色体大小和覆盖度校正的期望。我一般用cooltools的expected_trans# 计算反式期望按染色体对分组 expected_trans cooltools.expected_trans( clr, nproc8, chunksize1000000 ) # 查看 chr1 和 chr2 之间的期望接触 pair_vals expected_trans[ (expected_trans[chrom1] chr1) (expected_trans[chrom2] chr2) ] print(pair_vals[[chrom1, chrom2, balanced.avg]].head())得到期望后观测值除以期望就是富集分数observed/expected。富集分数大于 1 表示接触偏好小于 1 表示排斥。这一步是后续推断染色体间相互作用的基础不能跳过。3.3 用观测/期望矩阵找显著互作区域有了期望就可以逐 bin 算富集并做显著性检验。下面代码计算 chr1 和 chr2 之间的观测/期望矩阵并用泊松检验找显著 bin 对import numpy as np from scipy.stats import poisson # 提取 chr1 和 chr2 的观测矩阵 obs clr.matrix(balanceTrue).fetch(chr1, chr2) # 提取对应的期望值按距离或按染色体对 # 这里假设 expected_trans 已给出 chr1-chr2 的平均期望 exp_val pair_vals[balanced.avg].values[0] # 计算富集分数 enrichment obs / exp_val # 泊松检验观测计数是否显著高于期望 pvals poisson.sf(obs - 1, exp_val) # 校正多重检验取 FDR 0.05 from statsmodels.stats.multitest import multipletests flat_p pvals.flatten() flat_p flat_p[~np.isnan(flat_p)] reject, qvals, _, _ multipletests(flat_p, methodfdr_bh) print(显著 bin 对数量:, reject.sum())poisson.sf(obs - 1, exp_val)算的是观测值大于等于当前计数的概率适合计数数据。如果矩阵已经归一化成浮点可以改用负二项或正态近似。FDR 校正不能省否则几万个 bin 对里会有一堆假阳性。4. 避坑与排查HiC-HiC 分析里我踩过的 5 个坑4.1 矩阵对角线异常高顺式衰减曲线头部翘起现象画出的顺式衰减曲线在近距离处明显高于拟合线log-log 图头部翘起。原因未过滤自连片段和未消化完全的邻近片段这些接触集中在对角线附近。解决在expected_cis里增大ignore_diags一般设 2 到 5同时在pairtools parse阶段用--walks-policy过滤掉连续 walk 的 read。我通常先看对角线前 5 个 bin 的计数占比超过 20% 就说明过滤不够。4.2 反式接触富集分数全小于 1找不到任何显著互作现象所有染色体对的观测/期望都低于 1显著性检验没有阳性。原因期望值被高估常见于用全基因组平均反式接触做背景而你的样本里某些染色体覆盖度偏低。解决改用按染色体对分别估计期望或者用cooltools的expected_trans按染色体大小校正。另外检查归一化权重是否包含反式接触如果--cis-only估计的权重直接套到反式矩阵会系统性偏低。4.3 两个样本的矩阵无法比较bin 边界不一致现象样本 A 和样本 B 的矩阵维度不同或者相同坐标的 bin 对不上。原因建矩阵时用了不同的chrom.sizes或不同的 bin 大小甚至染色体命名不一致chr1vs1。解决统一用同一份chrom.sizes建矩阵前先cooler info对比 bin 数和染色体列表。如果命名不同用sed批量替换后再重新建库。这个坑血泪经验一旦矩阵建错后面所有统计都是白做。4.4 归一化后出现负值或 NaN现象cooler dump --balanced导出的矩阵里有负值或 NaN。原因ICE 迭代不收敛通常是因为某些 bin 覆盖度极低权重估计发散。解决在cooler balance里加--min-nnz 10要求每个 bin 至少有 10 个非零接触或者先合并低覆盖 bin。导出后检查 NaN 比例超过 1% 就要回头处理稀疏 bin。4.5 显著性检验假阳性爆炸FDR 校正后仍有上千个现象泊松检验后 FDR 0.05 的 bin 对数量远超预期。原因没有考虑空间自相关相邻 bin 的接触计数高度相关独立检验假设不成立。解决改用块状置换检验block permutation或对矩阵先做平滑再检验。我一般会把矩阵按 3x3 或 5x5 窗口平滑后再算富集牺牲一点分辨率换稳定性。如果必须用逐 bin 检验把 FDR 阈值收紧到 0.01 并报告效应量。5. 进阶技巧用 HiC-HiC 分析结果做染色体间互作推断的验证方法走到这一步你已经有归一化矩阵、期望曲线和显著 bin 对列表。但怎么确认这些推断不是归一化或背景模型的伪影我常用的验证方法是置换检验加独立数据复核。具体做法把接触对的位置在基因组上随机打乱保持每条染色体的接触总数不变重新建矩阵、算期望、跑显著性重复 100 次统计每次显著 bin 对数量。如果真实数据的显著数量落在置换分布 95% 分位数以内说明你的阳性结果可能只是背景波动。下面是一个简化的置换检验框架import numpy as np import cooler def permute_pairs(pairs_file, chrom_sizes, n_perm100): 打乱接触对位置保持染色体接触总数 # 读取原始 pairs pairs np.loadtxt(pairs_file, dtypestr, usecols(1,2,3,4)) chrom1, pos1, chrom2, pos2 pairs.T perm_counts [] for i in range(n_perm): # 在各自染色体上随机重排坐标 pos1_perm np.random.permutation(pos1) pos2_perm np.random.permutation(pos2) # 重新建矩阵并统计显著 bin 对此处省略建矩阵细节 # ... perm_counts.append(sig_count) return np.array(perm_counts) # 真实数据的显著数量 real_sig 1234 perm_counts permute_pairs(sorted.pairs, chrom.sizes) p_value (perm_counts real_sig).mean() print(f置换检验 p 值: {p_value:.4f})这个框架的关键是保持每条染色体的接触总数和染色体间配对总数不变只打乱位置。如果 p 值小于 0.05说明真实互作信号超出随机背景。另一个验证是用独立样本或独立酶切数据重复同一分析看显著 bin 对是否重叠。我一般要求至少两个独立样本中重叠的 bin 对才进入下游注释。最后说一个我自己的习惯每次跑完 HiC-HiC 分析先把顺式衰减曲线、反式富集分布、置换检验 p 值三张图放在一页上过一遍。如果衰减曲线斜率正常、反式富集不是全小于 1、置换检验 p 值显著才继续往下做注释。这个习惯帮我省掉了很多次回头重跑的时间。希望帮到你。本文还有配套的精品资源点击获取