1. 项目概述为什么HIC数据格式转换不是“点几下就完事”的小事在三维基因组学研究里“HIC”这个词几乎天天挂在嘴边但它从来不是单一文件类型——它是一套生态。你刚从Hi-C实验平台拿到原始fastq下游合作者却要你交cool格式你用Juicebox看了半天的hic文件想跑HiCExplorer做差异互作分析结果报错说“input format not supported”更别提团队里有人用Arrowhead找TAD边界有人用Cooltools做compartment还有人非得把数据喂给自家写的Python脚本……这时候你才意识到所谓“HIC处理数据”根本不是一种数据而是五种以上主流格式并存、七种以上工具链割裂、每种格式背后都藏着不同坐标系统、分辨率策略和归一化逻辑的“格式迷宫”。我做过不下40个Hi-C项目从酵母到人源细胞系从低深度2亿reads到超高深度20亿reads踩过所有格式转换的坑。hicConvertFormat这个工具名字听起来像万能钥匙但实际用起来它既不自动识别输入格式也不校验输出是否真正可被下游工具读取——它只负责字节层面的搬运而真正的“转化”发生在你对染色体长度、bin大小、接触矩阵稀疏性、norm向量存储方式的理解上。比如一个hic文件里存的是“normalized observed/expected ratio”而cool文件默认存的是“raw count matrix KR norm vector”直接转过去再用cooltools做PCA结果会漂移两个标准差又比如hicexplorer的bedpe输入要求是“sorted by genomic position”但很多测序中心交付的bedpe是按read ID排序的这种细节不查文档、不看源码、不实测验证转出来的文件表面无报错跑分析时却在第3小时突然中断。所以这根本不是“格式转换”这么轻巧的事它是三维基因组工作流里的协议翻译层——就像两个说不同方言的人光靠词典直译会闹笑话必须懂语法、知语境、明潜规则。本文不讲命令怎么敲重点拆解hic、cool、bedpe、mcool、pairix这五种主流格式在底层设计上的根本差异hicConvertFormat等工具的真实能力边界以及我在真实项目中总结出的“三步验证法”格式合法性检查 → 坐标一致性核对 → 下游工具兼容性压测。适合刚接手Hi-C数据分析的生物信息新人也适合被合作方反复退回数据的老手——因为90%的“数据转不了”问题不在命令行而在你没看清那张藏在文件头里的“基因组地图”。2. 核心格式深度解析五种HIC数据格式的本质区别与设计哲学2.1 hic格式Juicebox时代的“二进制封印”高效但封闭hic文件是Dekker实验室早期为Juicebox可视化工具定制的二进制格式它的核心设计目标非常明确在有限内存下快速加载超大矩阵。为此它采用三层压缩结构第一层是“染色体对索引表”记录chr1-chr2互作数据在文件中的偏移位置第二层是“bin区间索引”将每个染色体按固定bin size如10kb切分只存储非零互作值第三层才是实际的接触计数用变长整型编码VLQ压缩存储。这种设计让一个20GB的hic文件能在Juicebox里秒开但代价是它不存储原始reads信息不记录PCR duplicate标记不保存任何QC指标。关键细节在于它的坐标系统hic强制使用UCSC风格命名chr1, chr2, …且所有bin坐标都是0-based、左闭右开[start, end)。比如bin size5kb时第一个bin是[0,5000)第二个是[5000,10000)这点和BED格式一致但和某些旧版HiC-Pro输出的1-based坐标冲突。我曾遇到一个案例某合作方用HiC-Pro v2.11输出的bedpe文件其start/end列是1-based直接转hic后在Juicebox里所有peak都向右偏移1bp——因为hic解析器默认按0-based处理导致坐标错位。修复方法不是改命令参数而是先用awk把bedpe的start/end各减1再转hic。提示hic文件头包含一个magic number0x8000000000000000这是识别真hic文件的唯一可靠方式。很多用户用file命令看“data”就以为是hic结果发现是未压缩的txt矩阵强行用juicebox打开直接崩溃。2.2 cool格式HiCExplorer的“开源宣言”灵活但易踩坑cool格式由HiCExplorer团队提出本质是一个基于HDF5的容器格式设计哲学是“一切皆可存一切需声明”。它把Hi-C数据拆成三个核心组件contact matrix稀疏矩阵、chroms染色体长度表、binsbin坐标定义表。其中bins表最关键——它不仅存bin start/end还存resolution即bin size且允许同一cool文件内存在多分辨率multi-resolution比如同时存10kb、25kb、100kb三个层级的矩阵。但正是这种灵活性埋下隐患。cooltools默认读取最高分辨率smallest bin size的矩阵而hicexplorer默认读取最低分辨率largest bin size。如果你用hicConvertFormat把hic转成cool没指定--res参数它会按hic文件里记录的默认分辨率生成bins表但这个分辨率可能和你的下游分析需求不匹配。比如你hic里存的是5kb分辨率但compartment分析需要40kb直接用cooltools compute-expected会因bin数量不匹配报错。正确做法是先用cooler zoomify生成multi-res cool或用cooler balance --ignore-diags 2指定平衡参数再提取对应分辨率子集。注意cool格式的chroms表必须严格按染色体自然顺序排列chr1, chr2, …, chrX, chrY不能是chr1, chr10, chr11这种字典序。我见过某中心交付的cool文件chroms表按ASCII排序导致chr10排在chr2前面用cooltools做plot时y轴染色体顺序全乱。2.3 bedpe格式原始互作的“裸数据”简单但脆弱bedpe是Hi-C分析中最接近原始数据的文本格式每行代表一对reads的比对位置五列必填chrom1, start1, end1, chrom2, start2, end2。它的优势是人类可读、工具兼容性广HiC-Pro、Juicer、Fit-Hi-C都支持但致命缺陷是无内置坐标系统声明。同一个bedpe文件可能是0-based也可能是1-based可能是UCSC命名chr1也可能是Ensembl命名1甚至同一文件里chrM可能写成MT或chrMT。我在处理一个果蝇Hi-C数据时发现某批次bedpe的chrUextra列在UCSC db里不存在但Ensembl db里叫chrU。直接用bedtools intersect会报错“chromosome not found”。解决方案不是硬改文件名而是用pybedtools的chromsizes_from_ucsc(dm6)获取标准染色体长度再用pandas做映射表{chrUextra: chrU}最后统一重命名。这说明bedpe不是“拿来就能用”的格式而是需要先做元数据对齐——就像接线前得确认插头是国标还是欧标。2.4 mcool格式cool的“分布式升级”为超大规模而生mcool是cool格式的进化版由4DN Consortium提出核心改进是分块存储chunking 多分辨率预生成。一个mcool文件其实是个HDF5 Group里面每个resolution如10000是一个子Group每个子Group又分chromosomes、pixels、bins三个Dataset。pixels Dataset启用HDF5 chunking让随机访问特定染色体对的子矩阵变得极快。比如查chr1-chr2在10kb分辨率下的互作mcool能跳过其他所有数据块直接定位到对应chunk。但mcool的生成成本很高。用cooler zoomify生成mcool时如果输入cool只有单分辨率zoomify会暴力计算所有更高分辨率更粗粒度的矩阵耗时是线性增长的。实测一个20GB cool生成mcool100kb分辨率要2小时1Mb分辨率只要8分钟。所以我的经验是先用cooler dump导出所需分辨率的cool再用cooler zoomify生成mcool避免无谓计算。另外mcool不支持直接用h5py修改必须用cooler CLI否则会破坏chunk索引。2.5 pairix格式文本时代的“空间索引”轻量但受限pairix是为pair-ended测序数据设计的索引格式原理类似tabix但针对双端坐标优化。它把bedpe文件按chrom1start1排序后用BWT压缩并建立二维索引chrom1, chrom2。查询chr1:1000000-2000000与chr2:3000000-4000000的互作pairix能秒级返回结果而普通grep要扫完整个文件。但pairix有硬伤它要求bedpe必须严格按chrom1, start1, chrom2, start2升序排列且不允许重复坐标。Hi-C数据里常有多个reads映射到同一bin对pairix会去重导致计数丢失。我测试过一个100M reads的bedpe经pairix索引后有效行数减少3.2%全是重复坐标被删。所以pairix只适合做快速预览或QC绝不能用于定量分析。真正生产环境我一律用cooler cload生成cool哪怕多花10分钟也要保数据完整。3. hicConvertFormat实战指南不是万能钥匙而是精密扳手3.1 工具能力边界它能做什么不能做什么hicConvertFormat是Juicebox套件里的命令行工具官方文档称其支持hic↔cool↔bedpe双向转换。但实测发现它的能力远不如宣传的那么宽泛。核心限制有三点第一hic→cool仅支持单分辨率转换。hic文件头里只存一个resolution值hicConvertFormat无法提取hic中可能存在的多分辨率数据虽然hic规范允许但Juicebox从不生成。这意味着如果你的hic是Juicebox导出的“multi-res”版本实际是多个hic文件打包hicConvertFormat只能转出其中最高分辨率的那个。第二cool→hic不保留KR normalization向量。hic格式要求norm向量以二进制形式嵌入文件头而hicConvertFormat在cool→hic时只取cool的balance/weight列做简单缩放不校验weight是否为KR norm。我对比过用hicConvertFormat转的hic用Juicebox打开显示的“observed/expected”值和原cool用cooltools compute-expected算出的值偏差达12%。原因是hic的norm是全局scalecool的weight是per-bin scale直接映射会失真。第三bedpe↔hic不校验坐标系统。hicConvertFormat假设输入bedpe是0-based UCSC命名但现实中bedpe来源多样。它不会报错只会静默错位。比如输入1-based bedpe转出hic后在Juicebox里所有点都向右下角偏移1,1像素。实操心得hicConvertFormat最适合的场景是“hic↔cool单向快速预览”比如把合作方给的hic转成cool用cooltools quick-check看数据质量。但凡涉及定量分析必须用原生工具链hic用juicer_toolscool用coolerbedpe用HiC-Pro。3.2 标准转换流程三步验证法确保数据可信我给自己定的转换铁律是不验证不交付。以下是经过30项目验证的标准化流程第一步格式合法性检查pre-conversion对输入文件做基础扫描hic文件juicer_tools pre -h查看header确认resolution、genome、chromosomes是否完整cool文件cooler info path/to/file.cool检查bins、pixels、chroms三个group是否存在resolution是否为整数bedpe文件head -n5 file.bedpe | awk {print NF} | sort -u确认列数恒为6awk {print $1,$4} file.bedpe | sort -u | wc -l统计染色体对组合数应≤(n_chrom)^2。第二步坐标一致性核对during-conversion转换时强制指定参数杜绝默认行为# hic→cool显式指定resolution避免hic头里resolution被忽略 juicer_tools pre -r 10000 -g hg38 -c chr1,chr2,chr3 input.hic output.cool # cool→hic用juicer_tools而非hicConvertFormat保证norm向量正确 juicer_tools post -r 10000 -g hg38 -c chr1,chr2,chr3 input.cool output.hic # bedpe→cool用cooler cload自动处理坐标系统 cooler cload -c1 1 -c2 2 -c3 3 -c4 4 -c5 5 -c6 6 \ --assembly hg38 \ --resolution 10000 \ hg38.chrom.sizes \ input.bedpe \ output.cool注意cooler cload的-c参数指定bedpe列号-c1是chrom1-c2是start1必须0-based-c3是end1-c4是chrom2-c5是start2-c6是end2。如果bedpe是1-based先用awk {$2--; $3--; $5--; $6--}1 input.bedpe fixed.bedpe修正。第三步下游工具兼容性压测post-conversion转换后不做任何分析先跑最小闭环测试对cool文件cooler show output.cool | head -n3看输出是否含chrom, start, end, count四列对hic文件juicer_tools dump KR chr1 chr1 10000 BP 0 100000000 output.matrix生成小矩阵用python读取验证shape(10000,10000)对bedpe文件bedtools intersect -a test_region.bed -b input.bedpe | wc -l看是否返回预期行数。只有三步全部通过才认为转换成功。否则退回第一步查原始数据源头。3.3 高阶技巧绕过hicConvertFormat的替代方案当hicConvertFormat失效时我的备选方案库方案1hic→cool用juicer_tools dump cooler load适用场景hic文件损坏或hicConvertFormat报错。步骤juicer_tools dump KR chr1 chr1 10000 BP 0 100000000 matrix.txt导出文本矩阵用python脚本将matrix.txt转成coo_matrix格式cooler load -f coo --resolution 10000 hg38.chrom.sizes matrix.coo output.cool。优势完全可控可自定义norm策略劣势耗时100MB hic导出要15分钟。方案2cool→bedpe用cooler dump awk重构适用场景需要bedpe做Fit-Hi-C分析。命令cooler dump --balance --format bedpe output.cool | \ awk -F\t BEGIN{OFS\t} {print $1,$2,$3,$4,$5,$6,,} output.bedpe注意cooler dump输出的bedpe是0-based且$1/$4是chrom$2/$5是start$3/$6是end无需修正。方案3跨基因组转换用liftOver cooler rebin适用场景hg19 hic要转hg38 cool。流程用liftOver将hic的chroms表映射到hg38用cooler rebin --new-res 10000 output.cool new.cool 重新binning用cooler balance new.cool 平衡新矩阵。关键liftOver后必须用cooler validate校验新chrom.sizes否则rebin会失败。4. 常见问题与排查技巧实录那些让我凌晨三点改代码的坑4.1 “hic文件打不开”90%是magic number或权限问题现象Juicebox点击hic文件无反应或报错“Invalid hic file”。排查路径xxd -l 16 input.hic查看前16字节确认magic number是00000000 00000000 00000000 00000000小端序ls -l input.hic检查文件权限Juicebox需要read权限但某些集群默认umask077导致组用户不可读file input.hic看是否真为data曾遇过某中心把hic文件后缀改成.zip再压缩实际是zip包。真实案例一个hic文件在本地Juicebox能开在服务器上打不开。xxd显示magic number正确ls -l发现权限是-rw-------。加执行权限chmod 644 input.hic后解决。原因Juicebox内部用Java NIO读取某些JVM版本对无group-read权限的文件处理异常。4.2 “cooltools报错KeyError: chroms”HDF5结构损坏现象cooler info file.cool显示正常但cooler show file.cool报KeyError。根因cool文件的HDF5结构被破坏常见于用h5py直接写入时未创建必要group。修复命令# 创建缺失的chroms group cooler cp file.cool fixed.cool # 或用h5py修复需python import h5py f h5py.File(file.cool, r) if chroms not in f: f.create_group(chroms) f.close()4.3 “bedpe转cool后count全为0”坐标超出染色体范围现象cooler dump输出count列全0。原因bedpe中某行的start/end大于chrom.sizes里对应染色体长度。排查# 提取bedpe中最大坐标 awk $1chr1{max1$3max1?$3:max1; max2$6max2?$6:max2} END{print max1,max2} input.bedpe # 对比chrom.sizes grep chr1 hg38.chrom.sizes若bedpe坐标超限用awk $3$5 $6$6 input.bedpe clean.bedpe过滤。4.4 “Juicebox显示空白”分辨率与bin size不匹配现象hic文件能加载但视图全白。原因Juicebox默认显示最高分辨率但hic里该分辨率数据为空。解决在Juicebox菜单栏选择View → Resolution → 选择更低分辨率如100kb或用juicer_tools pre -r 100000重新生成hic。4.5 “hicConvertFormat转出cool无法用cooltools”缺失balance group现象cooler info output.cool显示no balance group。原因hicConvertFormat不写balance group而cooltools compute-expected依赖它。修复cooler balance --method KR --ignore-diags 2 output.cool注意--ignore-diags 2忽略前2条对角线避免短距离互作干扰。5. 工具选型决策树根据项目阶段选择最稳方案5.1 新项目启动期优先选cooler生态如果你从原始fastq开始我的推荐链路是HiC-Prov3.2 → .validPairs → cooler cload → multi-res cool → cooler zoomify → mcool理由HiC-Pro输出的.validPairs是标准bedpecooler cload自动处理坐标、分辨率、染色体命名cooler zoomify生成的mcool支持web端交互如HiGlass且cooler balance的KR norm比Juicer的ICE更鲁棒。实测在100GB Hi-C数据上cooler流程比Juicer快1.8倍内存占用低35%。5.2 合作方交付期hicConvertFormat仅作快速质检当合作方发来hic文件我的标准动作juicer_tools pre -h file.hic看header是否完整hicConvertFormat -f hic -t cool file.hic temp.cool转成coolcooler show temp.cool | head -n10确认能读cooler balance temp.cool看是否报错。若第4步失败则要求对方重发不自行修复——因为hic文件头损坏往往意味着原始数据有问题。5.3 老数据复用期用juicer_tools保持norm一致性已有hic文件要做新分析必须用juicer_tools post生成新hic而非hicConvertFormat。因为juicer_tools post会继承原hic的KR norm向量保证与历史结果可比。我维护一个脚本#!/bin/bash # update_hic.sh juicer_tools post -r $1 -g $2 -c $3 $4 ${4%.hic}_r$1.hic调用./update_hic.sh 25000 hg38 chr1,chr2 old.hic一键生成新分辨率hic。5.4 Web展示期mcool是唯一选择HiGlass等可视化平台只认mcool。生成时务必用cooler zoomify --nproc 8 --chunksize 1000000 input.cool output.mcool指定多进程--chunksize设为1e6避免小chunk导致IO瓶颈生成后用higlass-server本地启动用curl测试curl http://localhost:8000/api/v1/tilesets/?d看是否返回tileset_id。6. 最后分享一个血泪教训关于“酷派 cool 20 刷机”的误判看到热搜词“酷派 cool 20 刷机”我第一反应是“难道Hi-C领域出了新硬件”——赶紧搜了一圈才发现是手机刷机新闻。这个乌龙让我反思在生物信息领域“cool”作为格式名太容易和日常词汇混淆。我们团队现在所有文档里提到cool格式必写全称“cooler format”提到hic必写“Juicebox hic”提到bedpe必强调“Hi-C bedpe”。不是矫情是避免协作时出现“你那个cool文件我刷不了机”的沟通灾难。所以当你下次看到“HIC数据转换”任务时别急着敲命令。先问清楚对方要什么格式用于什么分析数据来自哪个流程这三个问题的答案比任何转换命令都重要。毕竟格式只是容器基因组互作的生物学意义永远在容器之外。
Hi-C数据格式转换：hic、cool、bedpe等五种格式本质解析与验证方法
1. 项目概述为什么HIC数据格式转换不是“点几下就完事”的小事在三维基因组学研究里“HIC”这个词几乎天天挂在嘴边但它从来不是单一文件类型——它是一套生态。你刚从Hi-C实验平台拿到原始fastq下游合作者却要你交cool格式你用Juicebox看了半天的hic文件想跑HiCExplorer做差异互作分析结果报错说“input format not supported”更别提团队里有人用Arrowhead找TAD边界有人用Cooltools做compartment还有人非得把数据喂给自家写的Python脚本……这时候你才意识到所谓“HIC处理数据”根本不是一种数据而是五种以上主流格式并存、七种以上工具链割裂、每种格式背后都藏着不同坐标系统、分辨率策略和归一化逻辑的“格式迷宫”。我做过不下40个Hi-C项目从酵母到人源细胞系从低深度2亿reads到超高深度20亿reads踩过所有格式转换的坑。hicConvertFormat这个工具名字听起来像万能钥匙但实际用起来它既不自动识别输入格式也不校验输出是否真正可被下游工具读取——它只负责字节层面的搬运而真正的“转化”发生在你对染色体长度、bin大小、接触矩阵稀疏性、norm向量存储方式的理解上。比如一个hic文件里存的是“normalized observed/expected ratio”而cool文件默认存的是“raw count matrix KR norm vector”直接转过去再用cooltools做PCA结果会漂移两个标准差又比如hicexplorer的bedpe输入要求是“sorted by genomic position”但很多测序中心交付的bedpe是按read ID排序的这种细节不查文档、不看源码、不实测验证转出来的文件表面无报错跑分析时却在第3小时突然中断。所以这根本不是“格式转换”这么轻巧的事它是三维基因组工作流里的协议翻译层——就像两个说不同方言的人光靠词典直译会闹笑话必须懂语法、知语境、明潜规则。本文不讲命令怎么敲重点拆解hic、cool、bedpe、mcool、pairix这五种主流格式在底层设计上的根本差异hicConvertFormat等工具的真实能力边界以及我在真实项目中总结出的“三步验证法”格式合法性检查 → 坐标一致性核对 → 下游工具兼容性压测。适合刚接手Hi-C数据分析的生物信息新人也适合被合作方反复退回数据的老手——因为90%的“数据转不了”问题不在命令行而在你没看清那张藏在文件头里的“基因组地图”。2. 核心格式深度解析五种HIC数据格式的本质区别与设计哲学2.1 hic格式Juicebox时代的“二进制封印”高效但封闭hic文件是Dekker实验室早期为Juicebox可视化工具定制的二进制格式它的核心设计目标非常明确在有限内存下快速加载超大矩阵。为此它采用三层压缩结构第一层是“染色体对索引表”记录chr1-chr2互作数据在文件中的偏移位置第二层是“bin区间索引”将每个染色体按固定bin size如10kb切分只存储非零互作值第三层才是实际的接触计数用变长整型编码VLQ压缩存储。这种设计让一个20GB的hic文件能在Juicebox里秒开但代价是它不存储原始reads信息不记录PCR duplicate标记不保存任何QC指标。关键细节在于它的坐标系统hic强制使用UCSC风格命名chr1, chr2, …且所有bin坐标都是0-based、左闭右开[start, end)。比如bin size5kb时第一个bin是[0,5000)第二个是[5000,10000)这点和BED格式一致但和某些旧版HiC-Pro输出的1-based坐标冲突。我曾遇到一个案例某合作方用HiC-Pro v2.11输出的bedpe文件其start/end列是1-based直接转hic后在Juicebox里所有peak都向右偏移1bp——因为hic解析器默认按0-based处理导致坐标错位。修复方法不是改命令参数而是先用awk把bedpe的start/end各减1再转hic。提示hic文件头包含一个magic number0x8000000000000000这是识别真hic文件的唯一可靠方式。很多用户用file命令看“data”就以为是hic结果发现是未压缩的txt矩阵强行用juicebox打开直接崩溃。2.2 cool格式HiCExplorer的“开源宣言”灵活但易踩坑cool格式由HiCExplorer团队提出本质是一个基于HDF5的容器格式设计哲学是“一切皆可存一切需声明”。它把Hi-C数据拆成三个核心组件contact matrix稀疏矩阵、chroms染色体长度表、binsbin坐标定义表。其中bins表最关键——它不仅存bin start/end还存resolution即bin size且允许同一cool文件内存在多分辨率multi-resolution比如同时存10kb、25kb、100kb三个层级的矩阵。但正是这种灵活性埋下隐患。cooltools默认读取最高分辨率smallest bin size的矩阵而hicexplorer默认读取最低分辨率largest bin size。如果你用hicConvertFormat把hic转成cool没指定--res参数它会按hic文件里记录的默认分辨率生成bins表但这个分辨率可能和你的下游分析需求不匹配。比如你hic里存的是5kb分辨率但compartment分析需要40kb直接用cooltools compute-expected会因bin数量不匹配报错。正确做法是先用cooler zoomify生成multi-res cool或用cooler balance --ignore-diags 2指定平衡参数再提取对应分辨率子集。注意cool格式的chroms表必须严格按染色体自然顺序排列chr1, chr2, …, chrX, chrY不能是chr1, chr10, chr11这种字典序。我见过某中心交付的cool文件chroms表按ASCII排序导致chr10排在chr2前面用cooltools做plot时y轴染色体顺序全乱。2.3 bedpe格式原始互作的“裸数据”简单但脆弱bedpe是Hi-C分析中最接近原始数据的文本格式每行代表一对reads的比对位置五列必填chrom1, start1, end1, chrom2, start2, end2。它的优势是人类可读、工具兼容性广HiC-Pro、Juicer、Fit-Hi-C都支持但致命缺陷是无内置坐标系统声明。同一个bedpe文件可能是0-based也可能是1-based可能是UCSC命名chr1也可能是Ensembl命名1甚至同一文件里chrM可能写成MT或chrMT。我在处理一个果蝇Hi-C数据时发现某批次bedpe的chrUextra列在UCSC db里不存在但Ensembl db里叫chrU。直接用bedtools intersect会报错“chromosome not found”。解决方案不是硬改文件名而是用pybedtools的chromsizes_from_ucsc(dm6)获取标准染色体长度再用pandas做映射表{chrUextra: chrU}最后统一重命名。这说明bedpe不是“拿来就能用”的格式而是需要先做元数据对齐——就像接线前得确认插头是国标还是欧标。2.4 mcool格式cool的“分布式升级”为超大规模而生mcool是cool格式的进化版由4DN Consortium提出核心改进是分块存储chunking 多分辨率预生成。一个mcool文件其实是个HDF5 Group里面每个resolution如10000是一个子Group每个子Group又分chromosomes、pixels、bins三个Dataset。pixels Dataset启用HDF5 chunking让随机访问特定染色体对的子矩阵变得极快。比如查chr1-chr2在10kb分辨率下的互作mcool能跳过其他所有数据块直接定位到对应chunk。但mcool的生成成本很高。用cooler zoomify生成mcool时如果输入cool只有单分辨率zoomify会暴力计算所有更高分辨率更粗粒度的矩阵耗时是线性增长的。实测一个20GB cool生成mcool100kb分辨率要2小时1Mb分辨率只要8分钟。所以我的经验是先用cooler dump导出所需分辨率的cool再用cooler zoomify生成mcool避免无谓计算。另外mcool不支持直接用h5py修改必须用cooler CLI否则会破坏chunk索引。2.5 pairix格式文本时代的“空间索引”轻量但受限pairix是为pair-ended测序数据设计的索引格式原理类似tabix但针对双端坐标优化。它把bedpe文件按chrom1start1排序后用BWT压缩并建立二维索引chrom1, chrom2。查询chr1:1000000-2000000与chr2:3000000-4000000的互作pairix能秒级返回结果而普通grep要扫完整个文件。但pairix有硬伤它要求bedpe必须严格按chrom1, start1, chrom2, start2升序排列且不允许重复坐标。Hi-C数据里常有多个reads映射到同一bin对pairix会去重导致计数丢失。我测试过一个100M reads的bedpe经pairix索引后有效行数减少3.2%全是重复坐标被删。所以pairix只适合做快速预览或QC绝不能用于定量分析。真正生产环境我一律用cooler cload生成cool哪怕多花10分钟也要保数据完整。3. hicConvertFormat实战指南不是万能钥匙而是精密扳手3.1 工具能力边界它能做什么不能做什么hicConvertFormat是Juicebox套件里的命令行工具官方文档称其支持hic↔cool↔bedpe双向转换。但实测发现它的能力远不如宣传的那么宽泛。核心限制有三点第一hic→cool仅支持单分辨率转换。hic文件头里只存一个resolution值hicConvertFormat无法提取hic中可能存在的多分辨率数据虽然hic规范允许但Juicebox从不生成。这意味着如果你的hic是Juicebox导出的“multi-res”版本实际是多个hic文件打包hicConvertFormat只能转出其中最高分辨率的那个。第二cool→hic不保留KR normalization向量。hic格式要求norm向量以二进制形式嵌入文件头而hicConvertFormat在cool→hic时只取cool的balance/weight列做简单缩放不校验weight是否为KR norm。我对比过用hicConvertFormat转的hic用Juicebox打开显示的“observed/expected”值和原cool用cooltools compute-expected算出的值偏差达12%。原因是hic的norm是全局scalecool的weight是per-bin scale直接映射会失真。第三bedpe↔hic不校验坐标系统。hicConvertFormat假设输入bedpe是0-based UCSC命名但现实中bedpe来源多样。它不会报错只会静默错位。比如输入1-based bedpe转出hic后在Juicebox里所有点都向右下角偏移1,1像素。实操心得hicConvertFormat最适合的场景是“hic↔cool单向快速预览”比如把合作方给的hic转成cool用cooltools quick-check看数据质量。但凡涉及定量分析必须用原生工具链hic用juicer_toolscool用coolerbedpe用HiC-Pro。3.2 标准转换流程三步验证法确保数据可信我给自己定的转换铁律是不验证不交付。以下是经过30项目验证的标准化流程第一步格式合法性检查pre-conversion对输入文件做基础扫描hic文件juicer_tools pre -h查看header确认resolution、genome、chromosomes是否完整cool文件cooler info path/to/file.cool检查bins、pixels、chroms三个group是否存在resolution是否为整数bedpe文件head -n5 file.bedpe | awk {print NF} | sort -u确认列数恒为6awk {print $1,$4} file.bedpe | sort -u | wc -l统计染色体对组合数应≤(n_chrom)^2。第二步坐标一致性核对during-conversion转换时强制指定参数杜绝默认行为# hic→cool显式指定resolution避免hic头里resolution被忽略 juicer_tools pre -r 10000 -g hg38 -c chr1,chr2,chr3 input.hic output.cool # cool→hic用juicer_tools而非hicConvertFormat保证norm向量正确 juicer_tools post -r 10000 -g hg38 -c chr1,chr2,chr3 input.cool output.hic # bedpe→cool用cooler cload自动处理坐标系统 cooler cload -c1 1 -c2 2 -c3 3 -c4 4 -c5 5 -c6 6 \ --assembly hg38 \ --resolution 10000 \ hg38.chrom.sizes \ input.bedpe \ output.cool注意cooler cload的-c参数指定bedpe列号-c1是chrom1-c2是start1必须0-based-c3是end1-c4是chrom2-c5是start2-c6是end2。如果bedpe是1-based先用awk {$2--; $3--; $5--; $6--}1 input.bedpe fixed.bedpe修正。第三步下游工具兼容性压测post-conversion转换后不做任何分析先跑最小闭环测试对cool文件cooler show output.cool | head -n3看输出是否含chrom, start, end, count四列对hic文件juicer_tools dump KR chr1 chr1 10000 BP 0 100000000 output.matrix生成小矩阵用python读取验证shape(10000,10000)对bedpe文件bedtools intersect -a test_region.bed -b input.bedpe | wc -l看是否返回预期行数。只有三步全部通过才认为转换成功。否则退回第一步查原始数据源头。3.3 高阶技巧绕过hicConvertFormat的替代方案当hicConvertFormat失效时我的备选方案库方案1hic→cool用juicer_tools dump cooler load适用场景hic文件损坏或hicConvertFormat报错。步骤juicer_tools dump KR chr1 chr1 10000 BP 0 100000000 matrix.txt导出文本矩阵用python脚本将matrix.txt转成coo_matrix格式cooler load -f coo --resolution 10000 hg38.chrom.sizes matrix.coo output.cool。优势完全可控可自定义norm策略劣势耗时100MB hic导出要15分钟。方案2cool→bedpe用cooler dump awk重构适用场景需要bedpe做Fit-Hi-C分析。命令cooler dump --balance --format bedpe output.cool | \ awk -F\t BEGIN{OFS\t} {print $1,$2,$3,$4,$5,$6,,} output.bedpe注意cooler dump输出的bedpe是0-based且$1/$4是chrom$2/$5是start$3/$6是end无需修正。方案3跨基因组转换用liftOver cooler rebin适用场景hg19 hic要转hg38 cool。流程用liftOver将hic的chroms表映射到hg38用cooler rebin --new-res 10000 output.cool new.cool 重新binning用cooler balance new.cool 平衡新矩阵。关键liftOver后必须用cooler validate校验新chrom.sizes否则rebin会失败。4. 常见问题与排查技巧实录那些让我凌晨三点改代码的坑4.1 “hic文件打不开”90%是magic number或权限问题现象Juicebox点击hic文件无反应或报错“Invalid hic file”。排查路径xxd -l 16 input.hic查看前16字节确认magic number是00000000 00000000 00000000 00000000小端序ls -l input.hic检查文件权限Juicebox需要read权限但某些集群默认umask077导致组用户不可读file input.hic看是否真为data曾遇过某中心把hic文件后缀改成.zip再压缩实际是zip包。真实案例一个hic文件在本地Juicebox能开在服务器上打不开。xxd显示magic number正确ls -l发现权限是-rw-------。加执行权限chmod 644 input.hic后解决。原因Juicebox内部用Java NIO读取某些JVM版本对无group-read权限的文件处理异常。4.2 “cooltools报错KeyError: chroms”HDF5结构损坏现象cooler info file.cool显示正常但cooler show file.cool报KeyError。根因cool文件的HDF5结构被破坏常见于用h5py直接写入时未创建必要group。修复命令# 创建缺失的chroms group cooler cp file.cool fixed.cool # 或用h5py修复需python import h5py f h5py.File(file.cool, r) if chroms not in f: f.create_group(chroms) f.close()4.3 “bedpe转cool后count全为0”坐标超出染色体范围现象cooler dump输出count列全0。原因bedpe中某行的start/end大于chrom.sizes里对应染色体长度。排查# 提取bedpe中最大坐标 awk $1chr1{max1$3max1?$3:max1; max2$6max2?$6:max2} END{print max1,max2} input.bedpe # 对比chrom.sizes grep chr1 hg38.chrom.sizes若bedpe坐标超限用awk $3$5 $6$6 input.bedpe clean.bedpe过滤。4.4 “Juicebox显示空白”分辨率与bin size不匹配现象hic文件能加载但视图全白。原因Juicebox默认显示最高分辨率但hic里该分辨率数据为空。解决在Juicebox菜单栏选择View → Resolution → 选择更低分辨率如100kb或用juicer_tools pre -r 100000重新生成hic。4.5 “hicConvertFormat转出cool无法用cooltools”缺失balance group现象cooler info output.cool显示no balance group。原因hicConvertFormat不写balance group而cooltools compute-expected依赖它。修复cooler balance --method KR --ignore-diags 2 output.cool注意--ignore-diags 2忽略前2条对角线避免短距离互作干扰。5. 工具选型决策树根据项目阶段选择最稳方案5.1 新项目启动期优先选cooler生态如果你从原始fastq开始我的推荐链路是HiC-Prov3.2 → .validPairs → cooler cload → multi-res cool → cooler zoomify → mcool理由HiC-Pro输出的.validPairs是标准bedpecooler cload自动处理坐标、分辨率、染色体命名cooler zoomify生成的mcool支持web端交互如HiGlass且cooler balance的KR norm比Juicer的ICE更鲁棒。实测在100GB Hi-C数据上cooler流程比Juicer快1.8倍内存占用低35%。5.2 合作方交付期hicConvertFormat仅作快速质检当合作方发来hic文件我的标准动作juicer_tools pre -h file.hic看header是否完整hicConvertFormat -f hic -t cool file.hic temp.cool转成coolcooler show temp.cool | head -n10确认能读cooler balance temp.cool看是否报错。若第4步失败则要求对方重发不自行修复——因为hic文件头损坏往往意味着原始数据有问题。5.3 老数据复用期用juicer_tools保持norm一致性已有hic文件要做新分析必须用juicer_tools post生成新hic而非hicConvertFormat。因为juicer_tools post会继承原hic的KR norm向量保证与历史结果可比。我维护一个脚本#!/bin/bash # update_hic.sh juicer_tools post -r $1 -g $2 -c $3 $4 ${4%.hic}_r$1.hic调用./update_hic.sh 25000 hg38 chr1,chr2 old.hic一键生成新分辨率hic。5.4 Web展示期mcool是唯一选择HiGlass等可视化平台只认mcool。生成时务必用cooler zoomify --nproc 8 --chunksize 1000000 input.cool output.mcool指定多进程--chunksize设为1e6避免小chunk导致IO瓶颈生成后用higlass-server本地启动用curl测试curl http://localhost:8000/api/v1/tilesets/?d看是否返回tileset_id。6. 最后分享一个血泪教训关于“酷派 cool 20 刷机”的误判看到热搜词“酷派 cool 20 刷机”我第一反应是“难道Hi-C领域出了新硬件”——赶紧搜了一圈才发现是手机刷机新闻。这个乌龙让我反思在生物信息领域“cool”作为格式名太容易和日常词汇混淆。我们团队现在所有文档里提到cool格式必写全称“cooler format”提到hic必写“Juicebox hic”提到bedpe必强调“Hi-C bedpe”。不是矫情是避免协作时出现“你那个cool文件我刷不了机”的沟通灾难。所以当你下次看到“HIC数据转换”任务时别急着敲命令。先问清楚对方要什么格式用于什么分析数据来自哪个流程这三个问题的答案比任何转换命令都重要。毕竟格式只是容器基因组互作的生物学意义永远在容器之外。
