software

Juicer Tools 实跑指南——接触去重、MAPQ、片段过滤与分辨率,哪些步骤真的改变了矩阵

作者:发布于 2026-11-15约 11 分钟阅读
建好接触图,还要核对计数。Juicer Tools 2.20.00 · pre 与 dump 实跑。01,相同位置会继续累加,3 条输入进入同一像素,不等于去重。02,双端质量取较低值筛选,12 条 → MAPQ 30 后 9 条。03,片段信息与过滤分开,显式开关和有效片段编号都要核验。

拿到一个能够打开的 Hi-C 文件,并不意味着比对、去重和质量过滤已经全部完成。矩阵创建程序可以忠实地汇总输入,却无法替分析者保证输入代表独立分子。另一方面,同一份接触记录在不同分辨率或过滤条件下看起来很不一样,也未必意味着出现了新的染色体结构。

本文下载安装官方发布的 Juicer Tools 2.20.00 文件,实际运行 pre 和 dump,用十二条人工接触记录检查矩阵构建与导出。这里测试的是 Juicer 体系中的矩阵工具,没有执行从 FASTQ 开始的完整 Juicer 流程,也没有进行真实读段比对、重复识别、环检测或性能比较。每条输入都由脚本构造,以便逐项核对输出。

建好接触图,还要核对计数。Juicer Tools 2.20.00 · pre 与 dump 实跑。01,相同位置会继续累加,3 条输入进入同一像素,不等于去重。02,双端质量取较低值筛选,12 条 → MAPQ 30 后 9 条。03,片段信息与过滤分开,显式开关和有效片段编号都要核验。
Juicer Tools 2.20.00 · pre 与 dump 实跑

一、先说清楚这次运行到哪一层

官方快速入门展示了使用 pre 生成接触文件、使用 dump 提取矩阵的基本路径。它与处理原始测序读段的整套流程存在不同输入边界。若输入已经是接触记录,就应说明上游如何产生这些记录,以及哪些质量控制已经执行。

本例使用十一字段的中等长度格式,记录名称、两端方向、染色体、位置、片段编号与两端比对质量。两条虚构染色体分别长五万和四万碱基,使用自编的染色体长度文件。这样不依赖真实参考,也不会把示例位置错误对应到某个基因。

十二条记录中包括位置相同的两条、交换两端位置的一条、落在同一分箱内的一条、三条含低质量端的记录,以及染色体内和染色体间接触。我们还单独构造了片段信息对照和带权重的输入,用来检查它们是否改变计数含义。

创建前先按染色体组合组织输入,并保存字段定义。格式相似不等于列顺序相同,尤其是从其他 pairs 文件转换时,不能只删去表头就当作另一种格式。原始位置的坐标约定、染色体命名和两端质量所在列,都应在转换后用少量记录验证。

二、pre 汇总同一像素,不等于识别重复分子

两条人工记录的两端位置均为一千一百和四千一百,第三条交换了这两个位置。以一千碱基分箱后,三条记录都落在同一个无向接触单元,输出计数为三。它们没有因为位置相同或两端交换而自动变成一次接触。

这个结果有两个含义。首先,染色体内接触的两端交换不会创造另一个独立矩阵单元。其次,单元计数增加并不能说明增加了独立分子,因为矩阵构建过程可以把重复输入一起相加。读段层面的重复识别应由能够使用相应原始信息的步骤完成。

java -Xmx2g -jar juicer_tools.2.20.00.jar pre \
  -n -r 1000,5000 contacts.txt all.hic chrom.sizes
java -Xmx2g -jar juicer_tools.2.20.00.jar dump \
  observed NONE all.hic chr1 chr1 BP 1000 chr1_1000.tsv

这里使用不计算归一化的构建选项,并在导出时明确选择未归一化的观察值。这样可以将输出逐项对应到人工输入,不把权重估计与原始汇总同时混进这个实验。对真实数据,也建议先确认原始矩阵的计数关系,再检查后续归一化结果。

文件名包含“去重复”字样并不是证据。一个接触表是否去重,需要查看产生它的流程和实际记录,而不是根据命名推断。若已经丢失读段位置细节或分子标识,事后仅凭一个较大的矩阵像素也无法判断它由真实高频接触还是重复分子贡献。

三、MAPQ 门槛作用于两端中的较低值

我们分别放入第一端质量为五、第二端质量为五,以及染色体间接触中一端质量为五的记录,另一端质量均为六十。另有一条记录的两端质量恰好都是三十,用来检查门槛边界。

未设置质量门槛时,全部十二条记录进入本例矩阵。将门槛设为三十后,三条含低质量端的记录被排除,留下九条。两端都恰好为三十的记录保留,符合本次版本的大于等于门槛行为。不能只看一端质量很高,就认为整对记录能够通过。

java -Xmx2g -jar juicer_tools.2.20.00.jar pre \
  -n -r 1000,5000 -q 30 contacts.txt q30.hic chrom.sizes

这项检查适合迁移到格式转换中。如果误把第二端质量写到了其他列,过滤数量可能仍然看似合理,却已经不再按预期筛选。设计一条只有第一端低质量、一条只有第二端低质量的对照,比只准备两端都很高或都很低的记录更容易发现这种错误。

质量阈值还依赖上游比对器如何定义和输出比对质量。本文的质量数字是人工写入,不能据此认定真实数据中三十就是普遍最优门槛。正式流程应结合读段类型、比对策略和研究目标确定标准,再记录被排除数量及其在样本间的差异。

四、只保留染色体内接触,会改变比例的分母

在质量门槛三十基础上,再启用仅保留染色体内接触的选项,计数从九变成七,因为剩下两条染色体间接触被排除。这里选项中的“对角”指染色体组合的对角块,并不意味着只保留分箱矩阵的主对角线。

本次构建设置 染色体内计数 染色体间计数 合计
无额外质量门槛 9 3 12
两端质量至少 30 7 2 9
质量至少 30,且仅染色体内 7 0 7

按这些输出分别计算,染色体内接触占比从百分之七十五变为约百分之七十七点七八,再变成百分之百。最后的百分之百完全来自主动排除染色体间记录,不能将它解释为文库质量提升。比较不同样本的这类指标时,必须使用相同过滤阶段和分母。

类似问题也会出现在短距离与长距离接触占比中。先删除某类记录,再拿剩余数据计算该类之外的比例,会机械地改变数字。质控指标应尽量对应清楚的输入阶段,保留筛选前后的计数,使读者能够区分数据本身的性质与分析者施加的条件。

五、片段文件与片段过滤不是同一个开关

pre 的官方说明同时介绍了限制性片段字段和片段图选项,但其中关于同片段记录的提示不能代替对安装版本的核验。我们最初按“同片段默认被排除”设计预期结果,实际输出却保留了那条记录,因此进一步增加了对照。

首先,把全部人工片段编号都改为零,但不启用片段过滤,输出仍然是十二次接触。其次,构造每五千碱基一个切点的虚拟片段表,并按位置重新计算两端片段编号。提供这个片段文件后,输出也仍然保留十二次接触。说明在本次二点二零版本和这些命令下,提供片段信息本身没有完成同片段排除。

随后显式启用同片段排除选项。在按五千碱基切点计算编号的输入中,五条染色体内同片段记录被移除,剩下七条。染色体间两端即使片段编号相同,也不会因这个条件被当作同一片段。这些结果均由实际导出矩阵逐项核对。

java -Xmx2g -jar juicer_tools.2.20.00.jar pre \
  -n -r 1000,5000 \
  -f restriction.txt --skip-intra-frag \
  fragment_contacts.txt fragment_skip.hic chrom.sizes

我们还对最初那份人工编号运行同一排除选项,只去掉一个编号相同的染色体内记录,留下十一条。它与按切点重新计算编号后的七条不同,说明过滤依据依赖输入片段编号。真实项目若片段编号不可信,即使开关设置正确,也无法得到可信的片段过滤结果。

官方参数解析实现包含这一独立开关。本文记录的是固定发布文件的实跑行为,不将它无条件推广到所有历史版本。迁移旧流程时,应把软件文件校验和与小规模对照一起保留,避免只沿用一段可能已经过时的说明。

六、分辨率变粗,像素变化不等于接触变化

在一千碱基分辨率下,那三条重复位置记录位于一千与四千对应的非对角单元,另外一条两端都在一千到两千范围内的接触位于主对角线。改成五千碱基分辨率后,四条记录都进入零到五千的同一个对角单元,计数为四。

原来分散的信号变成一个更强的对角像素,但接触记录没有增加。对于整条虚构第一染色体,不加质量门槛时,一千与五千分辨率的计数总和都为八。第二条染色体与染色体间部分的相应总和也分别保持为一和三。

这类变化说明,不同分辨率下的像素数、非零比例和对角线计数不能直接逐项比较。一个细分辨率下位于非对角线的接触,在粗分辨率下可能进入对角线。若某个后续步骤会排除对角像素,改变分辨率还会进一步改变参与分析的数据集合。

矩阵能生成很细的分箱,并不意味着数据支持同样精细的结构判断。分箱越细,同样数量的接触分散到更多单元,每个单元可用证据可能更少。应结合覆盖、重复间一致性和待回答的问题选择分辨率,而不是只比较图像放大后看起来是否更细。

七、NONE 不保证数值就是原始读段数

我们另外构造三行带分数的短格式输入,同一接触位置分别赋值二点五和一点二五,另一位置赋值四。实际构建再以未归一化方式导出后,第一个单元为三点七五,所有单元总和为七点七五。

这说明导出时不应用归一化,只表示没有在读取时执行所选的矩阵归一化,并不能保证输入最初就是每条记录权重为一的原始接触。若输入已经包含某种分数、缩放值或预处理后的数值,导出仍会保留这些来源特征。

带权输入 汇总到的单元 实际导出值
同一位置的 2.5 与 1.25 一千与四千的分箱对 3.75
另一位置的 4 八千与一万一千的分箱对 4.00
三行合计 两个单元 7.75

因此,交接接触文件时应说明数值列的来源,而不是只写某个归一化名称。特别是在多工具格式转换中,一张已经变换过的矩阵如果被再次当作原始计数归一化,可能改变后续解释。文件容器负责保存数值,数值代表什么仍然需要分析记录说明。

八、这次没有运行的分析,不能从矩阵直接推断

本文所有主要构建都显式关闭归一化,导出使用未归一化观察值,没有计算或验证 KR 等权重,也没有进行显著接触、染色体环或结构域检测。因此,示例中出现一个高计数像素,只能说明人工输入在该单元聚集,不能将它命名为真实染色体环。

Juicer Tools 的优势是提供成熟的接触文件构建与提取接口,便于将数据交给后续浏览和分析。它的使用代价,是必须明确输入格式、过滤阶段、分辨率和数值语义。相同文件扩展名不足以保证两份矩阵可以直接比较,尤其当它们的上游去重或质量筛选不同。

对于实际项目,可以先核对各染色体组合的计数总和,再抽取少量已知接触检查分箱位置,最后进入归一化和结构分析。每增加一个处理阶段,都说明它改变了哪些数据或数值,能让最终发现有清楚的证据链,而不是只依赖一张视觉上鲜明的热图。

九、复现实跑与这次排查的经验

下载本组完整脚本与实跑记录,在本篇目录执行下面命令即可生成接触表与矩阵。脚本使用官方固定发布地址,并校验 JAR 文件摘要。如果运行环境无法连接下载站点,可先从官方地址获取同名文件放到输出目录,再运行脚本完成校验和分析。

python3 run.py output

本次遇到的下载不可达问题通过在可连接的环境下载后传入解决,没有更换软件来源。脚本另外固定了 Java 运行容器,限制容器两核、四 GiB 内存,并保留各条命令与日志。这些限制只是教学执行环境,没有进行大规模耗时比较,也没有据此推断完整 Hi-C 流程的资源需求。

最值得保留的经验是用明确的正反例核对过滤行为。重复位置、单端低质量、门槛边界和同片段记录分别检查,可以定位究竟是哪条规则造成数量变化。只看最后的文件大小、成功提示或接触总数,无法区分这些原因。

当需要批量构建和审阅多样本、多分辨率接触矩阵时,可通过云生信生物分析服务器服务统一保存处理环境与分析记录。

把知识用于当前任务

下一步可以直接开始分析

按研究目标查找在线工具、课程、计算资源或完整分析服务。

按目标找方案仍有疑问,联系客服 ↗