当purge_dups“用力过猛”用Hi-C数据手动校正基因组重复片段的实战策略最近在整理几个高杂合度物种的组装项目时我反复遇到一个令人头疼的问题使用purge_dups进行重复序列过滤后BUSCO评估里的缺失率Missing会突然飙升有时甚至能翻两到三倍。看着那些辛苦组装出来的contig被成批标记为“冗余”而删除而后续的基因注释又提示一些关键的单拷贝同源基因不见了这种滋味确实不好受。如果你也正为此困扰觉得自动化的去冗余工具像一把不够精细的剃刀那么这篇文章或许能给你提供一个不同的思路——我们不必完全依赖算法“黑箱”的判决可以借助Hi-C数据提供的三维空间信息亲自上手像一位基因组“编辑”一样对重复片段进行可视化的审视和手动校正。这种方法的核心在于将Hi-C互作信号作为一种独立的、在物理空间上真实的约束条件来验证和修正基于序列相似度和覆盖度所做的判断。当两个contig在序列上高度相似被purge_dups判定为重复但它们的Hi-C互作模式却清晰地指向基因组上不同的、独立的位置时我们就有了充分的理由去质疑自动过滤的结果并将这些contig“抢救”回来。整个过程融合了Juicebox的可视化交互、Hi-C数据挂载Scaffolding的逻辑以及一些基于经验的判断技巧虽然比一键运行脚本更耗时但对于追求高质量、尤其是关注特定基因区域完整性的参考基因组构建来说往往是值得的。1. 理解问题根源为什么purge_dups会“误伤”在直接动手操作之前我们有必要先拆解一下purge_dups的工作原理以及它可能在哪些场景下“用力过猛”。这能帮助我们在后续手动校正时更有针对性地去检查哪些区域。purge_dups的判断逻辑主要基于两个维度测序深度Coverage和序列相似性Similarity。它的理想假设是在一个二倍体基因组中纯合区域的测序深度大约是杂合区域的两倍而由于单倍型分离或嵌合组装产生的冗余contig其深度会接近或低于杂合区域并且会与其他contig存在高度的序列相似性。然而现实中的基因组复杂性常常打破这个理想模型高度重复区域与近期串联重复一些基因家族如抗病基因R基因、某些转录因子家族存在大量的串联复制这些区域序列相似性极高但它们在基因组上是真实存在的不同拷贝承担着功能分化。purge_dups很容易将这些区域整体误判为单倍型冗余而删除。杂合度分布极度不均有些物种的基因组存在大片段的结构变异或杂合度“岛屿”导致局部区域的覆盖度分布偏离全局模型。purge_dups基于全局覆盖度分布设定的阈值low cutoff, mid cutoff, high cutoff在这些区域可能不再适用。组装错误导致的假相似性偶尔组装算法可能会在复杂区域产生一些局部的、错误的重复序列模式这些“假重复”也会被purge_dups捕获并过滤。注意不要完全抛弃purge_dups的结果。它的输出文件如dups.bed和中间文件覆盖度分布图PB.cov.png仍然是极有价值的诊断工具。我们的策略是将其作为一个“高灵敏度”的预警系统然后用Hi-C数据这个“高特异性”的验证工具去进行二次审查。为了更直观地理解哪些contig容易被误判我们可以看下面这个简化的对比表格Contig 特征purge_dups 可能标记Hi-C 互作模式提供的证据手动校正倾向序列相似度高覆盖度处于“杂合峰”HAPLOTIG(疑似单倍型冗余)互作信号集中于基因组一个特定区域且与“主”contig互作模式不同保留可能是独立基因拷贝覆盖度异常高 high cutoffREPEAT(高重复)互作信号分散与多个染色体区域有弱连接可能为真实串联重复或重复序列需结合注释判断覆盖度异常低 low cutoffJUNK(垃圾序列)几乎没有清晰的Hi-C互作信号很可能为污染或组装错误考虑删除与另一contig末端高度相似OVLP(末端重叠)两个contig的Hi-C互作图谱能无缝拼接可能是未正确连接的同一 scaffold尝试合并这个表格为我们后续在Juicebox中的观察提供了一个初步的“侦查指南”。2. 准备工作获取Hi-C互作矩阵与构建可视化基础手动校正的前提是获得一个可靠的、基于Hi-C数据的基因组挂载结果并将其加载到交互式可视化工具中。这里我们以常用的Juicer流程和Juicebox工具为例。2.1 生成Hi-C接触矩阵.hic文件首先你需要将Hi-C双端测序数据比对到你的contig集合上即purge_dups处理前的primary contigs并生成Juicebox可读的.hic文件。我通常使用Juicer的流程它封装了从比对到矩阵生成的全过程。# 假设你的基因组fasta文件为 hifiasm.primary.fa Hi-C数据为 hic_R1.fastq.gz 和 hic_R2.fastq.gz # 1. 创建参考基因组索引 bwa index hifiasm.primary.fa # 2. 运行Juicer脚本需要提前配置好Juicer的目录和依赖 # 这里展示一个简化版的命令核心实际需根据Juicer的README配置运行脚本 ./juicer.sh -g 基因组ID -d 工作目录 -p 染色体尺寸文件 \ -s 限制性内切酶 -y 限制性位点文件 \ -z hifiasm.primary.fa \ -q 队列名 -l 资源需求 \ -D juicer主目录 \ hic_R1.fastq.gz hic_R2.fastq.gz运行成功后在输出目录的aligned文件夹里你会找到最终生成的merged_nodups.txt文件。接下来使用Pre脚本juicer_tools.jar来生成.hic文件。# 3. 生成 .hic 文件 java -jar juicer_tools.jar pre -n \ aligned/merged_nodups.txt \ output.hic \ (awk {print $1\t$2} hifiasm.primary.fa.fai) # 从.fai文件生成染色体尺寸文件得到的output.hic文件就是包含了所有contig间互作强度的矩阵是后续可视化的基础。2.2 导入Juicebox并初步观察下载并打开Juicebox现在已集成到Juicebox Assembly Tools (JBAT)中。将.hic文件和你的基因组fasta文件导入。加载 .hic 文件: 在Juicebox中通过File - Open...选择你的.hic文件。加载 Assembly (fasta): 接着通过File - Load Assembly...加载你的hifiasm.primary.fa文件。这时你会看到基于Hi-C互作热图主对角线为contig内部互作颜色越深互作越强和contig列表。初始视图通常是全基因组级别的互作信号看起来可能像一团“毛球”。我们的第一步是寻找明显的染色体尺度chromosome-scale的互作块blocks。这些块表现为热图上远离主对角线的、明亮的正方形或矩形区域意味着属于不同contig的片段在三维空间上紧密相邻很可能位于同一条染色体上。3. 手动校正实战从可疑contig到决策现在进入核心环节。我们将结合purge_dups的输出在Juicebox中逐一审查被标记的contig。3.1 定位与审查被标记的contig首先从purge_dups输出的dups.bed或purge_dups.log中提取出所有被标记为HAPLOTIG或REPEAT的contig ID列表。在Juicebox中使用搜索框Search box直接输入contig ID快速定位到该contig在热图和列表中的位置。观察这个contig的内部互作模式主对角线上的信号是否清晰、连续这能反映其自身组装质量。最关键的一步观察这个contig与基因组其他所有contig的互作模式。你需要拖动视图仔细查看以这个contig所在行和列为起点的整个行和列。判断逻辑如下情况A该contig与另一个或少数几个contig有极强的、特异的互作信号。这强烈暗示它们位于同一条染色体的相邻位置。如果这两个contig恰好序列相似度很高被purge_dups标记为重复那么它们很可能是等位区域allelic regions即来自同源染色体的相同区段。由于Hi-C数据来自细胞群体等位区域在三维空间上也可能被拉近从而产生互作信号。这时你需要判断保留哪一个比较两者的长度和连续性N50。检查BUSCO基因是否只存在于其中一个上。通常保留更长、更完整、包含更多单拷贝BUSCO基因的那个。情况B该contig与基因组上多个不相邻的位点都有较弱的互作。这更符合重复序列元件如转座子的特征它们分散在基因组各处。如果该contig本身很短且被purge_dups标记为REPEAT那么它很可能就是真实的重复序列可以考虑过滤。情况C该contig几乎没有清晰的互作信号像一个“孤岛”。如果它同时被purge_dups标记为JUNK且覆盖度极低这很可能是外源污染或严重的组装错误。倾向于删除。3.2 利用“比对视图”进行精细比对Juicebox的一个强大功能是能够加载序列比对信息。你可以将purge_dups运行中的关键中间文件——基因组自比对结果*.self.paf——进行转换并加载进来。# 将 minimap2 自比对的 PAF 文件转换为 Juicebox 可读的 .pairs 或 .syn 格式 # 这里需要借助一些脚本例如 paftools.js (来自minimap2) 或自定义脚本进行格式转换 # 假设你已转换得到一个文件 duplicates.syn在Juicebox中通过File - Load Annotation...加载这个比对文件。它会以弧线arcs的形式显示在热图上方直观地展示出哪些contig之间存在大范围的序列相似性。当你在热图上看到一个可疑的互作信号时抬头看一眼是否有比对弧线连接着这两个contig这能立刻帮你确认这种空间邻近性是否源于序列相似性从而区分是真实的染色质交互还是重复序列造成的假象。3.3 做出编辑决策并导出新基因组在Juicebox中你可以直接对assembly进行编辑删除Contig在contig列表右键选择“Remove”。合并Contig如果你发现两个被标记为重复的contig实际上是同一个scaffold被错误打断的两部分它们的Hi-C互作图谱能完美衔接你可以尝试将它们合并Merge。但这需要谨慎最好有末端序列比对的支持。重排或调整方向对于复杂的嵌合体你可能还需要进行手动切割和重排。所有的编辑操作都会记录在Juicebox的“编辑历史”中。完成一轮审查和编辑后务必导出新的基因组fasta文件File - Export Assembly...。同时强烈建议将当前的“视图状态”包括你添加的注释、高亮等保存为.jbx文件方便日后回溯检查。4. 效果验证与迭代优化导出了手动校正后的基因组工作只完成了一半。必须进行系统的质量评估并与自动化过滤的结果进行对比。基础组装指标使用QUAST或BUSCO快速查看contig数量、N50、BUSCO完整性C%和重复率D%的变化。我们的目标是在BUSCO完整度C%基本不下降的前提下显著降低重复率D%同时严格控制缺失率M%的上升。busco -i manually_curated_genome.fa -m genome -l eukaryota_odb10 -o busco_manual -c 20Hi-C挂载质量评估将校正后的基因组重新进行Hi-C挂载可以使用3D-DNA或ALLHiC并计算挂载到染色体水平的比例、scaffold N50以及使用HiCExplorer的hicPlotMatrix观察校正后的热图是否更干净、染色体内部的交互矩阵矩阵块对角线是否更清晰。一个成功的校正应该能提升挂载的连续性和准确性。基因空间完整性检查如果你有该物种的RNA-seq数据或近缘物种的蛋白质序列运行一个基因注释流程如Braker3比较校正前后预测到的基因数量、尤其是那些单拷贝直系同源基因的完整性。这是检验是否误删了功能区域的最直接证据。这个过程往往是迭代式的。你可能需要根据验证结果再次打开Juicebox调整之前的某些决策。例如如果发现某个重要基因家族成员大量丢失就回去重新检查那些被删除的、包含该基因的重复contig看看它们是否真的有独特的Hi-C互作模式而被误杀了。手动校正基因组重复片段本质上是在计算自动化的效率与生物学判断的精度之间寻找平衡点。它要求研究者不仅熟悉工具流程更要理解数据背后的生物学意义。虽然耗时但对于构建一个真正可靠、能为后续功能研究提供坚实基础的参考基因组而言这种投入是至关重要的。每一次在Juicebox中放大热图仔细分辨那些模糊的互作信号时你都是在直接与基因组的三维结构对话这种体验本身也是生物信息学分析中一种独特的乐趣。