CD-HIT 这名字在生信圈里基本是绕不开的。只要做过基因组注释、宏基因组分析、蛋白家族筛选或者准备非冗余数据库大概率都会跟它打交道。它的核心作用就一句话对核酸或蛋白序列进行聚类去冗余用一条代表序列代替一堆相似序列从而极大缩小数据集规模。简单说它是个减法工具专治序列数据膨胀。这工具我用了好几年从早期的 4.x 版本用到现在的 4.8.1期间踩过不少坑也总结了一些自己的使用习惯。这篇东西不是单纯把官方文档翻译一遍而是想把我在实际项目中遇到的细节、参数调优的思路、以及那种“明明跑完了但结果不对”的排查过程写出来希望对刚开始接触 CD-HIT 的朋友有帮助。1. 为什么需要 CD-HIT序列去冗余到底在解决什么问题先聊一个最基础但很多人没认真想过的问题序列冗余是怎么来的以及它为什么会成为分析流程里的一个障碍。举个例子你做了一个宏基因组样本的组装得到几百万条 contig然后预测出几十万个基因。但这几十万个基因里有很多是高度相似的拷贝可能来自同一个物种的不同菌株或者同一个基因的不同片段。如果直接拿这几十万个基因去做后续的物种注释、功能注释或者丰度分析计算量会非常可观而且结果里会充斥着重复信息干扰你对样本真实多样性的判断。再比如你在做蛋白质家族分析从多个物种里收集了上万条同源蛋白序列。这些序列里有大量非常相近的成员如果直接拿去建 HMM 模型或者做多序列比对不仅慢而且近缘序列的偏向性会严重影响模型质量。这时候把序列在 90% 或 95% 相似度水平上去冗余保留下来的代表序列就能更均匀地反映整个家族的多样性。1.1 冗余数据的来源与影响冗余序列的来源多种多样不仅仅是生物学上的真实重复。我经常遇到的几种情况测序深度高导致的覆盖冗余。高深度测序下同一个基因组区域会被多条 reads 覆盖组装出来的 contig 之间天然存在大量重叠。同一物种不同菌株间的近缘序列。比如同一个物种的两株菌它们的 16S rRNA 基因序列可能 99% 相似这时候保留一条就够了。不同批次测序数据的合并。把多个样本的组装结果合并在一起跨样本的序列冗余非常惊人有时候去冗余后只剩三分之一。冗余数据带来的影响也相当实际。最直观的是计算资源浪费下游分析变慢磁盘占用变大。更隐蔽的是统计偏差问题比如某个基因家族在样本里有 100 个拷贝但其中 95 个几乎一模一样如果你不去冗余后面做丰度分析时这个家族的信号就会被严重放大误导生物学结论。1.2 去冗余的常见策略与 CD-HIT 的定位去冗余的思路其实不复杂本质上就是“序列两两比较相似度超过阈值就聚类”。但这里有个关键的工程问题如果直接用经典的 Needleman-Wunsch 或 Smith-Waterman 算法做全对全比对复杂度是 O(n²)几万条序列还能勉强跑几十万条序列就彻底歇菜了。所以去冗余工具的核心竞争力在于如何在尽量不损失准确性的前提下把比对速度提上去。目前的方案主要分几类基于 k-mer 索引和贪心聚类的工具代表就是 CD-HIT。基于图聚类和马尔可夫聚类的工具代表是 MMseqs2 的 linclust 模式还有 usearch、vsearch 这些。一些基于长读长组装算法思路的去冗余工具比如 purgedups主要用于基因组组装层级。CD-HIT 的优势在于它足够轻量、使用简单、内存占用可控而且几十万条序列这个规模下表现很稳定。虽然速度上可能不如 MMseqs2 那么激进但它的易用性和结果的可解释性让我一直把它作为首选工具。2. CD-HIT 核心机制拆解贪心聚类和 word 过滤是怎么回事CD-HIT 能快速处理大规模序列靠的不是什么黑魔法而是两个工程上的关键设计短 word 过滤策略和贪心聚类算法。搞清楚这两个机制你才能真正理解它的参数为什么这么设计以及什么场景下该调什么参数。2.1 贪心聚类的基本流程贪心算法的核心思想很简单先对序列按长度排序从最长的序列开始把它作为代表序列然后扫描其他序列看是否满足相似度阈值满足条件的就归入这个类不满足的就作为新的代表序列继续下一轮扫描。排序按长度降序是有讲究的。一条序列被聚类后跟代表序列比对时某些区域可能达不到阈值但更长的那条作为代表通常能代表更多的保守区域后续的类成员覆盖度也更容易达标。我在实际使用中体会是这个“贪心”策略虽然不是全局最优但胜在效率高结果的偏向性在大多数场景下是可以接受的。需要注意的一点是贪心聚类的结果跟序列输入顺序是有关的。这也是为什么 CD-HIT 会先按长度排序再聚类而不是按输入顺序直接处理。你如果做多轮聚类第一轮结果的代表序列进入第二轮最终的代表序列集合可能会因为第一轮的聚类阈值不同而有差异。2.2 word 过滤与“短单词”匹配的加速原理CD-HIT 最核心的加速手段是 word 过滤也叫 k-mer 过滤。它的原理是如果两条序列要满足 90% 的相似度那么它们必然共享很多长度为 w 的完全匹配的短片段word。于是 CD-HIT 会预先为每条序列建立 word 索引聚类时先快速扫描索引只对共享 word 数量达到一定阈值的序列对进行真正的比对。这里的 w 是参数-n它跟相似度阈值-c密切相关。因为相似度阈值越高序列间需要共享的短片段就越长、越严格word 长度可以设置得更大过滤效率也更高阈值越低序列间允许的差异越多word 必须更短否则容易漏掉真正相似的序列对。CD-HIT 官方文档里给出了一个经验对照表比如-c 0.9时建议-n 5-c 0.8时建议-n 4-c 0.7时建议-n 3。这个对照表不是随便写的它是基于大量测试得到的经验值需要记住的是-n太大容易漏掉相似序列-n太小会降低过滤效率、增加计算时间。我见过有人直接不设-n让程序自己判断这种情况 CD-HIT 会按版本内置的默认规则来定多数时候问题不大但如果你追求性能最优化建议手动指定。2.3 为什么说 CD-HIT 的结果是“足够好”而不是“精确”这里想提醒一点CD-HIT 的聚类结果并不保证是全局最优解。贪心策略 word 过滤决定了它是在“速度”和“准确度”之间做了权衡。在大多数应用场景里这种取舍是完全值得的。比如你做蛋白家族去冗余CD-HIT 聚出来的类可能不是最紧凑的但代表序列的生物学代表性没问题后续分析也不会被干扰。但如果你的研究对聚类边界非常敏感比如要严格定义单拷贝直系同源基因家族那 CD-HIT 可能不是最优选择OrthoFinder 这类基于系统发育信号的工具更合适。工具选型没有绝对的对错关键看你的下游分析对精度的要求有多高。3. 安装与基础用法从零开始跑通第一个聚类任务CD-HIT 的安装非常友好不管你是用 conda 的日常用户还是习惯源码编译的老派玩家都能在几分钟内搞定。这里我推荐优先走 conda 路线省心省力。3.1 Conda 安装与可用版本确认用 conda 安装只需要一行命令conda install -c bioconda cd-hit安装完成后可以用cd-hit --version确认版本。目前 bioconda 默认源里最新的稳定版本基本在 4.8.1这个版本很稳定也是我一直在用的版本。有几个需要注意的细节如果之前装过旧版 4.6.x建议升级到 4.8.x因为新版修复了一些内存管理的 bug对大文件处理更稳。CD-HIT 的二进制文件是区分线程版本的cd-hit本身是单线程版本cd-hit-mp才是并行版本且需要单独安装。在 conda 环境里cd-hit和cd-hit-mp可能是分开的包如果你需要并行建议显式安装cdhit和cdhit-mp两个包部分渠道名称不同需要以你所在环境为准。3.2 源码编译方式与依赖说明有些服务器环境没有外网或者集群的安全策略不允许用 conda这时候源码编译就是唯一选择了。CD-HIT 的源码编译非常轻量依赖只有 GNU C 编译器没有其他乱七八糟的库下载源码后进入目录直接执行wget https://github.com/weizhongli/cdhit/archive/refs/tags/v4.8.1.tar.gz tar -zxvf v4.8.1.tar.gz cd cdhit-4.8.1 make make install编译完成后可执行文件在cd-hit-est核酸序列版本、cd-hit蛋白序列版本和cd-hit-mp并行版本等名称下按需调用。如果你的系统里同时装了 GNU 和 Intel 编译器建议用 GNU 编译兼容性更好。3.3 蛋白序列聚类的最小完整示例安装完成后我们先跑一个最小示例确认工具能正常工作。假设你有一个蛋白序列文件proteins.fasta要按 90% 相似度聚类去冗余cd-hit -i proteins.fasta -o proteins_c90 -c 0.9 -n 5 -T 4 -M 2000解释一下每个参数-i指定输入文件必须是 FASTA 格式。-o指定输出文件前缀程序会生成proteins_c90代表序列和proteins_c90.clstr聚类结果文件。-c 0.9表示序列相似度阈值0.9 即 90%。-n 5表示 word size对于 90% 阈值CD-HIT 建议用 5。-T 4表示使用 4 个线程这个参数只在使用并行版本时生效单线程版本会忽略它。-M 2000表示内存上限设置为 2000 MB。运行日志里会显示类似这样的信息Total CPU time 23.45s The core of CD-HIT is done看到The core of CD-HIT is done就说明聚类完成了。这个示例流程如果你能顺利跑通那 CD-HIT 的基本使用就没有障碍了接下来可以深入理解参数。4. 参数选型的实操细节阈值、word size、比对模式怎么搭配CD-HIT 的参数不算多但每一个都值得认真对待。我见过太多人直接拿着默认参数跑完就收工结果出的结果跟预期差别很大然后又回来怀疑工具出了问题。实际上参数选型是 CD-HIT 使用中最需要考虑清楚的部分。4.1 相似度阈值-c的选择策略-c是 CD-HIT 最核心的参数它决定了两条序列被聚到同一类的严格程度。不同应用场景推荐的阈值差别很大应用场景推荐阈值说明蛋白序列去冗余同一蛋白家族0.9 - 0.95高阈值保留更多序列差异核酸序列去冗余菌株水平0.95 - 0.99株水平聚类需要非常严格宏基因组基因目录构建0.9 - 0.95平衡计算量和生物学意义构建非冗余参考数据库0.9蛋白/ 0.99核酸参考库通常希望尽量精简16S/ITS 等标记基因聚类0.97 - 0.99通常对应 OTU 或 ASV 聚类水平在这个基础上我还要分享一个经验阈值到底选多少最终应该由你的下游分析对“错误聚类”的容忍度决定。如果你做的是多序列比对和系统发育树构建建议阈值稍微高一点比如 0.95避免把旁系同源基因混在一起如果你只是想粗筛一下再去做精细聚类0.8 甚至 0.7 都可以。4.2 Word size-n的选择不要盲目跟风-n的选择逻辑我在前面原理部分已经提到了它跟-c强相关。CD-HIT 官方提供了一个建议表我可以直接抄来并用实际经验验证过-c 0.9对应-n 5-c 0.8对应-n 4-c 0.7对应-n 3-c 0.6对应-n 2-c 0.5对应-n 2如果你拿不准最简单的办法就是直接不设-n让 CD-HIT 根据-c自动分配。但如果你追求极致的速度可以尝试比建议值大 1 的-n比如-c 0.9时用-n 6看看会不会漏掉同源序列。我做过测试在 90% 阈值下-n 6偶发会漏掉一些边界情况的序列所以不建议为了提速牺牲召回率。还有一个容易忽略的点-n的范围在 CD-HIT 里是有限制的蛋白序列 1-5核酸序列 1-10超过上限程序会直接报错并提示你检查参数。4.3 比对覆盖度相关的参数-aL 和 -aS-c控制的是相似度但相似度有两种计算维度全局相似度和局部覆盖度。CD-HIT 里有两个参数来分别控制这两种维度-aL控制比对区域占短序列长度的最小比例取值范围 0-1默认 0.0。如果设成 0.8意味着短序列必须有 80% 以上的区域被比对上否则不聚类。-aS控制比对区域占短序列长度的最小比例默认 0.0。如果设成 0.8意味着短序列必须有 80% 以上的长度参与到比对中。这两个参数有什么区别举个例子一条 1000 bp 的序列和一条 300 bp 的序列如果相似度超过 95%但 300 bp 的短序列只有 250 bp 能比对到长序列上那全局相似度是 (250 × 2) / (1000 300) ≈ 0.38达不到 0.9 阈值程序会认为它们不聚类。但如果你设了-aL 0并且只要求短序列的覆盖度-aS就可能把这种部分覆盖的情况聚类到一起。实际操作里如果你做的是全长序列的去冗余-aL和-aS可以保持默认但如果你做的是部分序列、片段化组装结果或宏基因组分箱后的序列建议设置-aS 0.8甚至更高避免把只有部分同源的序列硬凑在一起。4.4 其他重要参数-g、-T、-M、-d这几个参数虽然不像-c和-n那么核心但也非常影响运行效果和资源占用的冲突点。-g这个参数控制是否进行全局最优比对。默认值是 0表示按“先到先得”的方式聚类即只要满足阈值就把序列归入当前类设置为 1 时CD-HIT 会尝试寻找该类最优的代表序列。代价是速度变慢、内存占用增加。我的建议是如果你使用-c 0.9以上做严格聚类可以设置-g 1结果更干净如果是粗筛保持默认即可。-T线程数。CD-HIT 的并行版本支持多线程-T 0表示使用所有可用核心-T 4表示用 4 核。需要注意单线程版本忽略这个参数所以如果你发现设了-T没有效果先检查一下自己用的是不是cd-hit而不是cd-hit-mp。-M内存上限单位是 MB。CD-HIT 对内存的控制比较保守如果你不设这一项它默认会用满所有可用内存。在共享服务器上跑任务建议设置一个合理的内存上限比如-M 4000避免被管理员找上门。-d聚类结果文件.clstr里每条序列名称显示的长度默认 20。这个参数纯粹影响展示效果试过几次后我觉得改成 100 比较舒服否则序列名太长被截断查起来很麻烦。示例cd-hit -i proteins.fasta -o proteins_c90 -c 0.9 -n 5 -g 1 -T 4 -M 4000 -d 1004.5 核酸序列聚类与 cd-hit-est蛋白序列用cd-hit核酸序列则要用cd-hit-est。两者的参数逻辑基本一致但默认 word size 的合理范围不同。核酸序列的-n可以到 10但一般建议在 8-10 之间。cd-hit-est -i contigs.fasta -o contigs_c95 -c 0.95 -n 10 -T 4 -M 4000这个命令用于对基因组 contig 做 95% 相似度的去冗余。这里有一个经验如果输入是基因组 contig 而不是完整的基因序列建议先做一轮更严格的过滤比如先按 99% 去冗余再做一次 95% 的聚类效果会更好。因为 contig 通常比较长一步聚类到 95% 可能会把不同物种的相近片段混在一起。5. 输出文件解读与后续衔接很多人跑完 CD-HIT 之后就只拿了那个代表序列文件完全忽略了.clstr文件这其实很可惜。.clstr文件的信息量远比你想象的大它不仅仅是辅助文件更是你判断聚类质量的依据。5.1.clstr文件的格式详解.clstr文件是纯文本格式用 tab 键分割字段。每一条序列会输出以下信息序列 ID每条序列在输入文件中的标识符默认显示前 20 个字符由-d控制。聚类分组编号从 0 开始递增的整数。序列所属类每个类包含的序列数。长度序列长度bp 或 aa。链方向对于核酸序列可能有正向或反向互补正两条链的相对关系。实际打开一个.clstr文件它看起来是这样的Cluster 0 0 200aa, seq1... 1 198aa, seq2...其中Cluster 0表示第 0 个类每一行前面的数字是该序列在类中的编号后面是长度和序列 ID。如果一个类的代表序列后面带着*表示它是这个类的代表序列。这个文件的用途是你可以根据.clstr文件检查聚类质量比如看一个类里序列长度差异是否过大、代表序列是否合理。如果发现某个类里长度差异非常大说明这个类可能聚类过头了需要调整-aS或-c。5.2 提取代表序列和类成员列表的常用操作基于.clstr文件你可以做很多后处理。这里分享几个我经常用的命令提取每个类的序列数统计grep ^ proteins_c90.clstr | awk -F {print $2} | sort | uniq -c | sort -k2 -n提取某个特定类的所有序列 IDawk /^/{cluster$0; next} {print cluster, $0} proteins_c90.clstr | grep Cluster 5如果你想把聚类结果转成“代表序列到所有成员的映射表”可以用 Python 快速处理这也是我在做下游分析时最常用的方式之一。5.3 与下游分析工具InterProScan、eggNOG-mapper 等的衔接去冗余之后代表序列通常要进入下游功能注释流程。比如我会将proteins_c90作为输入喂给 InterProScan 或 eggNOG-mapper 做功能注释这样能省下不少计算时间同时避免近缘序列重复命中导致的冗余注释。这里有个实践经验对序列做功能注释时建议保留一份“代表序列到原始序列 ID 的映射表”这样当你拿到代表序列的功能注释结果后可以回溯到所有原始序列。否则下游分析做完你会发现根本不知道哪条代表序列对应哪些原始序列那种感觉非常痛苦。建映射表的小技巧是在聚类前就自己做一个 ID 对应关系。比如给每条输入序列改一个简洁的唯一 ID聚类后再从.clstr文件重建完整的 ID 映射。这个方法虽然有点笨但非常可靠。6. 常见问题与排查技巧实录这部分是我想重点写的。CD-HIT 用起来虽然简单但实际跑数据时各种问题层出不穷。我把自己遇到过的高频问题列出来希望能帮你省点时间。6.1 内存不足与段错误这是 CD-HIT 最经典的问题之一尤其在处理几十万条序列的大文件时内存消耗会突然暴涨。表现为程序运行到一半直接退出终端报Segmentation fault或者Killed。排查思路按顺序来先确认输入文件有没有问题。如果 FASTA 文件格式不规范比如序列行有非法字符、注释行没有以开头可能导致解析异常进而内存异常。检查-M参数是否合理。如果你把-M设得太低比如 1000CD-HIT 在构建索引时可能会被迫频繁读写临时文件速度变慢极端情况下也可能报错。-M 0表示不限制内存但这样会占满所有可用内存不建议在共享服务器上这样做。如果文件合法、参数也合理但内存还是暴涨可能是 word index 构建阶段出了问题。可以尝试降低-n比如从 5 降到 4减少索引量。最后如果你用的是平行版本cd-hit-mp注意它有个已知问题线程数设得过高会导致内存占用非线性增长。-T 8和-T 4的内存差别不是 2 倍可能是 3 倍以上。提示如果你用cd-hit-mp遇到内存问题时优先降低线程数而不是内存上限。默认情况下-T 0会用满全部核心这在大文件场景下其实隐患很大我建议显式指定-T 4或-T 8再配合-M限制总内存。6.2 “Out of memory”但服务器内存还很充足这个坑很隐蔽。有时候你会发现服务器明明还有 256 GB 内存但 CD-HIT 报错说内存不够。原因通常是你的编译版本是 32 位的无法寻址超过 4 GB 的内存空间。解决办法很直接重新编译成 64 位版本。如果你用的 conda 安装的二进制一般就是 64 位的这个问题主要出现在源码编译时。编译前检查一下系统架构确保是在 64 位环境下编译的。另外还有一种情况你设置了-M但设置的单位搞错了。-M的单位是 MB不是 GB。把-M 2000当成 2000 GB自然会报内存不足。这个错误很傻瓜但确实碰见过同事掉进去。6.3 聚类结果里出现大量单序列类代表序列数几乎没减少这是很多新手会懵的问题跑了半天去冗余结果输出文件大小几乎没变每个类里就一条序列。这种情况的原因通常有几种阈值设得太高。90% 的阈值本身就会保留很多单序列类尤其是数据本身多样性很高时。序列长度差太大。如果输入里有大量长度差异很大的序列短序列可能达不到覆盖度要求导致不聚类。这时候可以适当调整-aS。输入序列的 ID 里有特殊字符。有些序列 ID 里有空格或|导致解析异常间接影响聚类结果。建议输入前统一把 ID 里的空格替换成下划线。遇到这种情况先不要急着调阈值先检查数据本身。我试过一次把序列 ID 里的空格替换掉之后聚类效果立刻有了明显改善。6.4 并行版本的线程设置无效如果你用的是cd-hit-mp但发现运行时 CPU 占用率始终只有 100%可能是你用了cd-hit的单线程版本。cd-hit和cd-hit-mp是两个独立的可执行文件前者直接忽略-T参数。还有一个容易忽略的点某些 Linux 发行版里conda 会同时安装两个版本但默认的cd-hit软链接指向了单线程版。验证方法很直观which cd-hit which cd-hit-mp如果确认用的是cd-hit-mp但线程还是起不来检查一下你的 shell 环境变量偶尔会遇到旧版本路径遮蔽新版本的情况。6.5 代表序列会比原始序列变短或变少这个“问题”其实是 CD-HIT 的正常行为不算 bug。因为 CD-HIT 聚类时代表序列是从类里选择的一条序列而不是多条序列的保守 consensus。CD-HIT 不会生成共识序列它只是挑了一条最长的或其他规则下选出的序列作为代表。所以代表序列的长度完全取决于被选中那条序列的长度可能比类里某些成员短但不可能凭空变成共识序列。如果你需要每个类生成一条 consensus 序列后续可以用consensus工具对每个类的成员重新比对生成但这已经超出 CD-HIT 的职责范围了它是聚类工具不是多序列比对工具。6.6 输入文件不是标准 FASTA 格式导致的奇怪报错CD-HIT 对 FASTA 格式要求比较严格但报错信息往往不直观。常见问题包括序列行里含有数字或空格比如序列中间有-或.这些字符会导致 CD-HIT 在解析阶段行为异常。注释行开头之后的序列行里包含非标准氨基酸字符蛋白序列里出现了 U 或 *CD-HIT 会保留这些字符但可能会影响后续比对。文件结尾没有换行符有时候会导致最后一条序列解析不完整。这类问题的排查思路也比较简单先用seqkit stats查看输入文件的基本统计信息再用seqkit grep筛掉含异常字符的序列。我自己有个习惯所有输入 CD-HIT 的 FASTA 文件都会先用seqkit seq -g -w 0做一次清洗去掉 gaps 和非法字符。7. 加速与资源优化大规模序列的实战调优CD-HIT 虽然已经很快但当你处理真正的大规模数据时还能通过一些技巧再挤压出不少性能空间。这里分享几个我在实际项目中验证有效的优化手段。7.1 利用序列预过滤减少输入规模在做 CD-HIT 之前先做一轮轻量级的预过滤可以显著减少 CD-HIT 的计算量。预过滤的思路是用seqkit按长度过滤去掉过短或过长的序列具体阈值根据你的研究目标而定。用seqkit rmdup或dedupe工具先去重完全一样的序列减少 CD-HIT 处理冗余信息。比如你有一个 500 万条序列的文件但其中有很大一部分是完全重复的。先去重后可能只剩 300 万条这节省下来的时间非常可观。7.2 多轮聚类策略由粗到细的分步聚类我发现一个比较好用的套路先做一个低阈值的粗聚类比如 0.8然后在每一类内部再做一个高阈值的精细聚类比如 0.95。这样分步聚类的好处是每一步的工作量被分解得更小且每一轮的输入文件都更集中整体速度反而比一步到位更快。比如你的任务是按 95% 相似度对一堆蛋白序列聚类。可以这样操作# 第一轮粗聚类到 80% cd-hit -i all_proteins.fasta -o proteins_c80 -c 0.8 -n 4 # 提取代表序列 # 第二轮对代表序列细聚类到 95% cd-hit -i proteins_c80 -o proteins_c95 -c 0.95 -n 5这样两轮跑下来最终的代表序列集合可能跟直接一步-c 0.95得到的结果很接近但整体运行时间往往更低。7.3 用环境变量提高文件系统 IO 性能CD-HIT 在运行过程中会创建大量临时文件。如果你的/tmp目录是普通的机械硬盘而你的项目目录在 SSD 上可以考虑把临时目录指向更快的位置export TMPDIR/path/to/ssd/tmp cd-hit -i input.fasta -o output -c 0.9虽然这个操作看起来微不足道但在处理百万条序列时IO 性能的差距会被放大得非常明显。另外输出文件最好写到项目本地避免写到网络文件系统上。8. 实战示例从原始 FASTQ 到非冗余蛋白序列的完整流水线理论说了那么多不如直接看一个完整的实战案例。这里我用自己做过的一个宏基因组项目为例展示从原始测序 reads 到非冗余蛋白序列目录的完整流水线。8.1 案例背景与分析设计这个项目的目标是构建一个人类肠道宏基因组基因目录。原始数据来自 50 个样本的宏基因组测序每个样本约 10 Gb。分析流程是质控 宿主序列过滤。样本单独组装。对所有样本的 contig 做基因预测。将所有预测得到的蛋白序列合并成一个文件。用 CD-HIT 去冗余得到非冗余基因目录。整个流程最耗时的步骤就是第 5 步。合并后的蛋白序列文件大概有 600 万条序列如果不做策略优化直接丢给 CD-HIT 很可能会跑半天甚至直接内存爆掉。8.2 CD-HIT 命令与参数选择我先做了一步预处理# 将所有样本的蛋白序列合并 cat sample*/prodigal/*.faa all_orfs.faa # 清除非法字符和过短序列 seqkit seq -g -w 0 all_orfs.faa all_orfs_clean.faa seqkit seq -m 30 -M 5000 all_orfs_clean.faa all_orfs_filtered.faa然后用 CD-HIT 做两轮聚类# 第一轮粗聚类90% 相似度 cd-hit -i all_orfs_filtered.faa -o gene_catalog_c90 -c 0.9 -n 5 -T 8 -M 8000 -d 100 # 提取代表序列 # 第二轮细聚类95% 相似度 cd-hit -i gene_catalog_c90 -o gene_catalog_c95 -c 0.95 -n 5 -T 8 -M 8000 -d 100两轮跑完后最终得到约 180 万条代表蛋白序列。这个规模跟类似研究的基因目录大小一致说明聚类效果是合理的。8.3 分析结果的质控与可视化拿到结果后我通常会做以下几步质控统计聚类前后序列数量计算去冗余率。检查.clstr文件里每个类的序列数分布确认没有异常巨大的类。随机抽样一些类人工检查代表序列是否合理。用 eggNOG-mapper 对代表序列做功能注释评估功能覆盖度是否足够。其中第 2 点很关键。如果出现一个类包含了几万条序列的情况说明这个类可能是由重复序列或低复杂度区域导致的需要特别关注。9. 最后聊几句实操心得CD-HIT 这个工具最让我满意的一点是它的稳定性和可预测性。不像一些新潮的工具版本更新频繁算法动不动就大变样CD-HIT 的架构非常稳定很多年前写的流程现在依然能流畅运行这种可维护性在生信工具里相当难得。我手头就有几条两年前跑的流程现在原封不动拿出来重跑结果还是一模一样。踩过几次坑之后我的体会是CD-HIT 使用最大的门槛不是工具本身而是你得清楚自己的数据适合什么阈值、什么参数。多试几次不同的-c和-n组合观察聚类结果的变化这是最快的学习路径。工具本身不会告诉你阈值设多少合适只有对数据的理解才是最终的判断依据。如果你刚开始接触 CD-HIT建议先从一个小规模的测试集入手比如选 1 万条序列分别跑-c 0.8、-c 0.9、-c 0.95三个阈值对比聚类结果的差异感受一下参数对结果的影响。这个过程花不了多少时间但对后面处理大规模数据非常有帮助。最后再分享一个小技巧CD-HIT 的.clstr文件一定要保留好不要因为占空间就随手删掉。很多时候你后面需要回溯某个代表序列到底包含了哪些原始序列这时候.clstr就是唯一线索。我自己会在项目目录里专门建一个clustering/子目录把输入、输出、参数和日志文件都整理好这样即使隔了很久再回头看项目也能快速搞清楚当时的分析逻辑。 
长河君 · 2026/10/4 12:44:52 · 阅读 1.2k CD-HIT序列去冗余实用指南：参数调优与实战排错
CD-HIT 这名字在生信圈里基本是绕不开的。只要做过基因组注释、宏基因组分析、蛋白家族筛选或者准备非冗余数据库大概率都会跟它打交道。它的核心作用就一句话对核酸或蛋白序列进行聚类去冗余用一条代表序列代替一堆相似序列从而极大缩小数据集规模。简单说它是个减法工具专治序列数据膨胀。这工具我用了好几年从早期的 4.x 版本用到现在的 4.8.1期间踩过不少坑也总结了一些自己的使用习惯。这篇东西不是单纯把官方文档翻译一遍而是想把我在实际项目中遇到的细节、参数调优的思路、以及那种“明明跑完了但结果不对”的排查过程写出来希望对刚开始接触 CD-HIT 的朋友有帮助。1. 为什么需要 CD-HIT序列去冗余到底在解决什么问题先聊一个最基础但很多人没认真想过的问题序列冗余是怎么来的以及它为什么会成为分析流程里的一个障碍。举个例子你做了一个宏基因组样本的组装得到几百万条 contig然后预测出几十万个基因。但这几十万个基因里有很多是高度相似的拷贝可能来自同一个物种的不同菌株或者同一个基因的不同片段。如果直接拿这几十万个基因去做后续的物种注释、功能注释或者丰度分析计算量会非常可观而且结果里会充斥着重复信息干扰你对样本真实多样性的判断。再比如你在做蛋白质家族分析从多个物种里收集了上万条同源蛋白序列。这些序列里有大量非常相近的成员如果直接拿去建 HMM 模型或者做多序列比对不仅慢而且近缘序列的偏向性会严重影响模型质量。这时候把序列在 90% 或 95% 相似度水平上去冗余保留下来的代表序列就能更均匀地反映整个家族的多样性。1.1 冗余数据的来源与影响冗余序列的来源多种多样不仅仅是生物学上的真实重复。我经常遇到的几种情况测序深度高导致的覆盖冗余。高深度测序下同一个基因组区域会被多条 reads 覆盖组装出来的 contig 之间天然存在大量重叠。同一物种不同菌株间的近缘序列。比如同一个物种的两株菌它们的 16S rRNA 基因序列可能 99% 相似这时候保留一条就够了。不同批次测序数据的合并。把多个样本的组装结果合并在一起跨样本的序列冗余非常惊人有时候去冗余后只剩三分之一。冗余数据带来的影响也相当实际。最直观的是计算资源浪费下游分析变慢磁盘占用变大。更隐蔽的是统计偏差问题比如某个基因家族在样本里有 100 个拷贝但其中 95 个几乎一模一样如果你不去冗余后面做丰度分析时这个家族的信号就会被严重放大误导生物学结论。1.2 去冗余的常见策略与 CD-HIT 的定位去冗余的思路其实不复杂本质上就是“序列两两比较相似度超过阈值就聚类”。但这里有个关键的工程问题如果直接用经典的 Needleman-Wunsch 或 Smith-Waterman 算法做全对全比对复杂度是 O(n²)几万条序列还能勉强跑几十万条序列就彻底歇菜了。所以去冗余工具的核心竞争力在于如何在尽量不损失准确性的前提下把比对速度提上去。目前的方案主要分几类基于 k-mer 索引和贪心聚类的工具代表就是 CD-HIT。基于图聚类和马尔可夫聚类的工具代表是 MMseqs2 的 linclust 模式还有 usearch、vsearch 这些。一些基于长读长组装算法思路的去冗余工具比如 purgedups主要用于基因组组装层级。CD-HIT 的优势在于它足够轻量、使用简单、内存占用可控而且几十万条序列这个规模下表现很稳定。虽然速度上可能不如 MMseqs2 那么激进但它的易用性和结果的可解释性让我一直把它作为首选工具。2. CD-HIT 核心机制拆解贪心聚类和 word 过滤是怎么回事CD-HIT 能快速处理大规模序列靠的不是什么黑魔法而是两个工程上的关键设计短 word 过滤策略和贪心聚类算法。搞清楚这两个机制你才能真正理解它的参数为什么这么设计以及什么场景下该调什么参数。2.1 贪心聚类的基本流程贪心算法的核心思想很简单先对序列按长度排序从最长的序列开始把它作为代表序列然后扫描其他序列看是否满足相似度阈值满足条件的就归入这个类不满足的就作为新的代表序列继续下一轮扫描。排序按长度降序是有讲究的。一条序列被聚类后跟代表序列比对时某些区域可能达不到阈值但更长的那条作为代表通常能代表更多的保守区域后续的类成员覆盖度也更容易达标。我在实际使用中体会是这个“贪心”策略虽然不是全局最优但胜在效率高结果的偏向性在大多数场景下是可以接受的。需要注意的一点是贪心聚类的结果跟序列输入顺序是有关的。这也是为什么 CD-HIT 会先按长度排序再聚类而不是按输入顺序直接处理。你如果做多轮聚类第一轮结果的代表序列进入第二轮最终的代表序列集合可能会因为第一轮的聚类阈值不同而有差异。2.2 word 过滤与“短单词”匹配的加速原理CD-HIT 最核心的加速手段是 word 过滤也叫 k-mer 过滤。它的原理是如果两条序列要满足 90% 的相似度那么它们必然共享很多长度为 w 的完全匹配的短片段word。于是 CD-HIT 会预先为每条序列建立 word 索引聚类时先快速扫描索引只对共享 word 数量达到一定阈值的序列对进行真正的比对。这里的 w 是参数-n它跟相似度阈值-c密切相关。因为相似度阈值越高序列间需要共享的短片段就越长、越严格word 长度可以设置得更大过滤效率也更高阈值越低序列间允许的差异越多word 必须更短否则容易漏掉真正相似的序列对。CD-HIT 官方文档里给出了一个经验对照表比如-c 0.9时建议-n 5-c 0.8时建议-n 4-c 0.7时建议-n 3。这个对照表不是随便写的它是基于大量测试得到的经验值需要记住的是-n太大容易漏掉相似序列-n太小会降低过滤效率、增加计算时间。我见过有人直接不设-n让程序自己判断这种情况 CD-HIT 会按版本内置的默认规则来定多数时候问题不大但如果你追求性能最优化建议手动指定。2.3 为什么说 CD-HIT 的结果是“足够好”而不是“精确”这里想提醒一点CD-HIT 的聚类结果并不保证是全局最优解。贪心策略 word 过滤决定了它是在“速度”和“准确度”之间做了权衡。在大多数应用场景里这种取舍是完全值得的。比如你做蛋白家族去冗余CD-HIT 聚出来的类可能不是最紧凑的但代表序列的生物学代表性没问题后续分析也不会被干扰。但如果你的研究对聚类边界非常敏感比如要严格定义单拷贝直系同源基因家族那 CD-HIT 可能不是最优选择OrthoFinder 这类基于系统发育信号的工具更合适。工具选型没有绝对的对错关键看你的下游分析对精度的要求有多高。3. 安装与基础用法从零开始跑通第一个聚类任务CD-HIT 的安装非常友好不管你是用 conda 的日常用户还是习惯源码编译的老派玩家都能在几分钟内搞定。这里我推荐优先走 conda 路线省心省力。3.1 Conda 安装与可用版本确认用 conda 安装只需要一行命令conda install -c bioconda cd-hit安装完成后可以用cd-hit --version确认版本。目前 bioconda 默认源里最新的稳定版本基本在 4.8.1这个版本很稳定也是我一直在用的版本。有几个需要注意的细节如果之前装过旧版 4.6.x建议升级到 4.8.x因为新版修复了一些内存管理的 bug对大文件处理更稳。CD-HIT 的二进制文件是区分线程版本的cd-hit本身是单线程版本cd-hit-mp才是并行版本且需要单独安装。在 conda 环境里cd-hit和cd-hit-mp可能是分开的包如果你需要并行建议显式安装cdhit和cdhit-mp两个包部分渠道名称不同需要以你所在环境为准。3.2 源码编译方式与依赖说明有些服务器环境没有外网或者集群的安全策略不允许用 conda这时候源码编译就是唯一选择了。CD-HIT 的源码编译非常轻量依赖只有 GNU C 编译器没有其他乱七八糟的库下载源码后进入目录直接执行wget https://github.com/weizhongli/cdhit/archive/refs/tags/v4.8.1.tar.gz tar -zxvf v4.8.1.tar.gz cd cdhit-4.8.1 make make install编译完成后可执行文件在cd-hit-est核酸序列版本、cd-hit蛋白序列版本和cd-hit-mp并行版本等名称下按需调用。如果你的系统里同时装了 GNU 和 Intel 编译器建议用 GNU 编译兼容性更好。3.3 蛋白序列聚类的最小完整示例安装完成后我们先跑一个最小示例确认工具能正常工作。假设你有一个蛋白序列文件proteins.fasta要按 90% 相似度聚类去冗余cd-hit -i proteins.fasta -o proteins_c90 -c 0.9 -n 5 -T 4 -M 2000解释一下每个参数-i指定输入文件必须是 FASTA 格式。-o指定输出文件前缀程序会生成proteins_c90代表序列和proteins_c90.clstr聚类结果文件。-c 0.9表示序列相似度阈值0.9 即 90%。-n 5表示 word size对于 90% 阈值CD-HIT 建议用 5。-T 4表示使用 4 个线程这个参数只在使用并行版本时生效单线程版本会忽略它。-M 2000表示内存上限设置为 2000 MB。运行日志里会显示类似这样的信息Total CPU time 23.45s The core of CD-HIT is done看到The core of CD-HIT is done就说明聚类完成了。这个示例流程如果你能顺利跑通那 CD-HIT 的基本使用就没有障碍了接下来可以深入理解参数。4. 参数选型的实操细节阈值、word size、比对模式怎么搭配CD-HIT 的参数不算多但每一个都值得认真对待。我见过太多人直接拿着默认参数跑完就收工结果出的结果跟预期差别很大然后又回来怀疑工具出了问题。实际上参数选型是 CD-HIT 使用中最需要考虑清楚的部分。4.1 相似度阈值-c的选择策略-c是 CD-HIT 最核心的参数它决定了两条序列被聚到同一类的严格程度。不同应用场景推荐的阈值差别很大应用场景推荐阈值说明蛋白序列去冗余同一蛋白家族0.9 - 0.95高阈值保留更多序列差异核酸序列去冗余菌株水平0.95 - 0.99株水平聚类需要非常严格宏基因组基因目录构建0.9 - 0.95平衡计算量和生物学意义构建非冗余参考数据库0.9蛋白/ 0.99核酸参考库通常希望尽量精简16S/ITS 等标记基因聚类0.97 - 0.99通常对应 OTU 或 ASV 聚类水平在这个基础上我还要分享一个经验阈值到底选多少最终应该由你的下游分析对“错误聚类”的容忍度决定。如果你做的是多序列比对和系统发育树构建建议阈值稍微高一点比如 0.95避免把旁系同源基因混在一起如果你只是想粗筛一下再去做精细聚类0.8 甚至 0.7 都可以。4.2 Word size-n的选择不要盲目跟风-n的选择逻辑我在前面原理部分已经提到了它跟-c强相关。CD-HIT 官方提供了一个建议表我可以直接抄来并用实际经验验证过-c 0.9对应-n 5-c 0.8对应-n 4-c 0.7对应-n 3-c 0.6对应-n 2-c 0.5对应-n 2如果你拿不准最简单的办法就是直接不设-n让 CD-HIT 根据-c自动分配。但如果你追求极致的速度可以尝试比建议值大 1 的-n比如-c 0.9时用-n 6看看会不会漏掉同源序列。我做过测试在 90% 阈值下-n 6偶发会漏掉一些边界情况的序列所以不建议为了提速牺牲召回率。还有一个容易忽略的点-n的范围在 CD-HIT 里是有限制的蛋白序列 1-5核酸序列 1-10超过上限程序会直接报错并提示你检查参数。4.3 比对覆盖度相关的参数-aL 和 -aS-c控制的是相似度但相似度有两种计算维度全局相似度和局部覆盖度。CD-HIT 里有两个参数来分别控制这两种维度-aL控制比对区域占短序列长度的最小比例取值范围 0-1默认 0.0。如果设成 0.8意味着短序列必须有 80% 以上的区域被比对上否则不聚类。-aS控制比对区域占短序列长度的最小比例默认 0.0。如果设成 0.8意味着短序列必须有 80% 以上的长度参与到比对中。这两个参数有什么区别举个例子一条 1000 bp 的序列和一条 300 bp 的序列如果相似度超过 95%但 300 bp 的短序列只有 250 bp 能比对到长序列上那全局相似度是 (250 × 2) / (1000 300) ≈ 0.38达不到 0.9 阈值程序会认为它们不聚类。但如果你设了-aL 0并且只要求短序列的覆盖度-aS就可能把这种部分覆盖的情况聚类到一起。实际操作里如果你做的是全长序列的去冗余-aL和-aS可以保持默认但如果你做的是部分序列、片段化组装结果或宏基因组分箱后的序列建议设置-aS 0.8甚至更高避免把只有部分同源的序列硬凑在一起。4.4 其他重要参数-g、-T、-M、-d这几个参数虽然不像-c和-n那么核心但也非常影响运行效果和资源占用的冲突点。-g这个参数控制是否进行全局最优比对。默认值是 0表示按“先到先得”的方式聚类即只要满足阈值就把序列归入当前类设置为 1 时CD-HIT 会尝试寻找该类最优的代表序列。代价是速度变慢、内存占用增加。我的建议是如果你使用-c 0.9以上做严格聚类可以设置-g 1结果更干净如果是粗筛保持默认即可。-T线程数。CD-HIT 的并行版本支持多线程-T 0表示使用所有可用核心-T 4表示用 4 核。需要注意单线程版本忽略这个参数所以如果你发现设了-T没有效果先检查一下自己用的是不是cd-hit而不是cd-hit-mp。-M内存上限单位是 MB。CD-HIT 对内存的控制比较保守如果你不设这一项它默认会用满所有可用内存。在共享服务器上跑任务建议设置一个合理的内存上限比如-M 4000避免被管理员找上门。-d聚类结果文件.clstr里每条序列名称显示的长度默认 20。这个参数纯粹影响展示效果试过几次后我觉得改成 100 比较舒服否则序列名太长被截断查起来很麻烦。示例cd-hit -i proteins.fasta -o proteins_c90 -c 0.9 -n 5 -g 1 -T 4 -M 4000 -d 1004.5 核酸序列聚类与 cd-hit-est蛋白序列用cd-hit核酸序列则要用cd-hit-est。两者的参数逻辑基本一致但默认 word size 的合理范围不同。核酸序列的-n可以到 10但一般建议在 8-10 之间。cd-hit-est -i contigs.fasta -o contigs_c95 -c 0.95 -n 10 -T 4 -M 4000这个命令用于对基因组 contig 做 95% 相似度的去冗余。这里有一个经验如果输入是基因组 contig 而不是完整的基因序列建议先做一轮更严格的过滤比如先按 99% 去冗余再做一次 95% 的聚类效果会更好。因为 contig 通常比较长一步聚类到 95% 可能会把不同物种的相近片段混在一起。5. 输出文件解读与后续衔接很多人跑完 CD-HIT 之后就只拿了那个代表序列文件完全忽略了.clstr文件这其实很可惜。.clstr文件的信息量远比你想象的大它不仅仅是辅助文件更是你判断聚类质量的依据。5.1.clstr文件的格式详解.clstr文件是纯文本格式用 tab 键分割字段。每一条序列会输出以下信息序列 ID每条序列在输入文件中的标识符默认显示前 20 个字符由-d控制。聚类分组编号从 0 开始递增的整数。序列所属类每个类包含的序列数。长度序列长度bp 或 aa。链方向对于核酸序列可能有正向或反向互补正两条链的相对关系。实际打开一个.clstr文件它看起来是这样的Cluster 0 0 200aa, seq1... 1 198aa, seq2...其中Cluster 0表示第 0 个类每一行前面的数字是该序列在类中的编号后面是长度和序列 ID。如果一个类的代表序列后面带着*表示它是这个类的代表序列。这个文件的用途是你可以根据.clstr文件检查聚类质量比如看一个类里序列长度差异是否过大、代表序列是否合理。如果发现某个类里长度差异非常大说明这个类可能聚类过头了需要调整-aS或-c。5.2 提取代表序列和类成员列表的常用操作基于.clstr文件你可以做很多后处理。这里分享几个我经常用的命令提取每个类的序列数统计grep ^ proteins_c90.clstr | awk -F {print $2} | sort | uniq -c | sort -k2 -n提取某个特定类的所有序列 IDawk /^/{cluster$0; next} {print cluster, $0} proteins_c90.clstr | grep Cluster 5如果你想把聚类结果转成“代表序列到所有成员的映射表”可以用 Python 快速处理这也是我在做下游分析时最常用的方式之一。5.3 与下游分析工具InterProScan、eggNOG-mapper 等的衔接去冗余之后代表序列通常要进入下游功能注释流程。比如我会将proteins_c90作为输入喂给 InterProScan 或 eggNOG-mapper 做功能注释这样能省下不少计算时间同时避免近缘序列重复命中导致的冗余注释。这里有个实践经验对序列做功能注释时建议保留一份“代表序列到原始序列 ID 的映射表”这样当你拿到代表序列的功能注释结果后可以回溯到所有原始序列。否则下游分析做完你会发现根本不知道哪条代表序列对应哪些原始序列那种感觉非常痛苦。建映射表的小技巧是在聚类前就自己做一个 ID 对应关系。比如给每条输入序列改一个简洁的唯一 ID聚类后再从.clstr文件重建完整的 ID 映射。这个方法虽然有点笨但非常可靠。6. 常见问题与排查技巧实录这部分是我想重点写的。CD-HIT 用起来虽然简单但实际跑数据时各种问题层出不穷。我把自己遇到过的高频问题列出来希望能帮你省点时间。6.1 内存不足与段错误这是 CD-HIT 最经典的问题之一尤其在处理几十万条序列的大文件时内存消耗会突然暴涨。表现为程序运行到一半直接退出终端报Segmentation fault或者Killed。排查思路按顺序来先确认输入文件有没有问题。如果 FASTA 文件格式不规范比如序列行有非法字符、注释行没有以开头可能导致解析异常进而内存异常。检查-M参数是否合理。如果你把-M设得太低比如 1000CD-HIT 在构建索引时可能会被迫频繁读写临时文件速度变慢极端情况下也可能报错。-M 0表示不限制内存但这样会占满所有可用内存不建议在共享服务器上这样做。如果文件合法、参数也合理但内存还是暴涨可能是 word index 构建阶段出了问题。可以尝试降低-n比如从 5 降到 4减少索引量。最后如果你用的是平行版本cd-hit-mp注意它有个已知问题线程数设得过高会导致内存占用非线性增长。-T 8和-T 4的内存差别不是 2 倍可能是 3 倍以上。提示如果你用cd-hit-mp遇到内存问题时优先降低线程数而不是内存上限。默认情况下-T 0会用满全部核心这在大文件场景下其实隐患很大我建议显式指定-T 4或-T 8再配合-M限制总内存。6.2 “Out of memory”但服务器内存还很充足这个坑很隐蔽。有时候你会发现服务器明明还有 256 GB 内存但 CD-HIT 报错说内存不够。原因通常是你的编译版本是 32 位的无法寻址超过 4 GB 的内存空间。解决办法很直接重新编译成 64 位版本。如果你用的 conda 安装的二进制一般就是 64 位的这个问题主要出现在源码编译时。编译前检查一下系统架构确保是在 64 位环境下编译的。另外还有一种情况你设置了-M但设置的单位搞错了。-M的单位是 MB不是 GB。把-M 2000当成 2000 GB自然会报内存不足。这个错误很傻瓜但确实碰见过同事掉进去。6.3 聚类结果里出现大量单序列类代表序列数几乎没减少这是很多新手会懵的问题跑了半天去冗余结果输出文件大小几乎没变每个类里就一条序列。这种情况的原因通常有几种阈值设得太高。90% 的阈值本身就会保留很多单序列类尤其是数据本身多样性很高时。序列长度差太大。如果输入里有大量长度差异很大的序列短序列可能达不到覆盖度要求导致不聚类。这时候可以适当调整-aS。输入序列的 ID 里有特殊字符。有些序列 ID 里有空格或|导致解析异常间接影响聚类结果。建议输入前统一把 ID 里的空格替换成下划线。遇到这种情况先不要急着调阈值先检查数据本身。我试过一次把序列 ID 里的空格替换掉之后聚类效果立刻有了明显改善。6.4 并行版本的线程设置无效如果你用的是cd-hit-mp但发现运行时 CPU 占用率始终只有 100%可能是你用了cd-hit的单线程版本。cd-hit和cd-hit-mp是两个独立的可执行文件前者直接忽略-T参数。还有一个容易忽略的点某些 Linux 发行版里conda 会同时安装两个版本但默认的cd-hit软链接指向了单线程版。验证方法很直观which cd-hit which cd-hit-mp如果确认用的是cd-hit-mp但线程还是起不来检查一下你的 shell 环境变量偶尔会遇到旧版本路径遮蔽新版本的情况。6.5 代表序列会比原始序列变短或变少这个“问题”其实是 CD-HIT 的正常行为不算 bug。因为 CD-HIT 聚类时代表序列是从类里选择的一条序列而不是多条序列的保守 consensus。CD-HIT 不会生成共识序列它只是挑了一条最长的或其他规则下选出的序列作为代表。所以代表序列的长度完全取决于被选中那条序列的长度可能比类里某些成员短但不可能凭空变成共识序列。如果你需要每个类生成一条 consensus 序列后续可以用consensus工具对每个类的成员重新比对生成但这已经超出 CD-HIT 的职责范围了它是聚类工具不是多序列比对工具。6.6 输入文件不是标准 FASTA 格式导致的奇怪报错CD-HIT 对 FASTA 格式要求比较严格但报错信息往往不直观。常见问题包括序列行里含有数字或空格比如序列中间有-或.这些字符会导致 CD-HIT 在解析阶段行为异常。注释行开头之后的序列行里包含非标准氨基酸字符蛋白序列里出现了 U 或 *CD-HIT 会保留这些字符但可能会影响后续比对。文件结尾没有换行符有时候会导致最后一条序列解析不完整。这类问题的排查思路也比较简单先用seqkit stats查看输入文件的基本统计信息再用seqkit grep筛掉含异常字符的序列。我自己有个习惯所有输入 CD-HIT 的 FASTA 文件都会先用seqkit seq -g -w 0做一次清洗去掉 gaps 和非法字符。7. 加速与资源优化大规模序列的实战调优CD-HIT 虽然已经很快但当你处理真正的大规模数据时还能通过一些技巧再挤压出不少性能空间。这里分享几个我在实际项目中验证有效的优化手段。7.1 利用序列预过滤减少输入规模在做 CD-HIT 之前先做一轮轻量级的预过滤可以显著减少 CD-HIT 的计算量。预过滤的思路是用seqkit按长度过滤去掉过短或过长的序列具体阈值根据你的研究目标而定。用seqkit rmdup或dedupe工具先去重完全一样的序列减少 CD-HIT 处理冗余信息。比如你有一个 500 万条序列的文件但其中有很大一部分是完全重复的。先去重后可能只剩 300 万条这节省下来的时间非常可观。7.2 多轮聚类策略由粗到细的分步聚类我发现一个比较好用的套路先做一个低阈值的粗聚类比如 0.8然后在每一类内部再做一个高阈值的精细聚类比如 0.95。这样分步聚类的好处是每一步的工作量被分解得更小且每一轮的输入文件都更集中整体速度反而比一步到位更快。比如你的任务是按 95% 相似度对一堆蛋白序列聚类。可以这样操作# 第一轮粗聚类到 80% cd-hit -i all_proteins.fasta -o proteins_c80 -c 0.8 -n 4 # 提取代表序列 # 第二轮对代表序列细聚类到 95% cd-hit -i proteins_c80 -o proteins_c95 -c 0.95 -n 5这样两轮跑下来最终的代表序列集合可能跟直接一步-c 0.95得到的结果很接近但整体运行时间往往更低。7.3 用环境变量提高文件系统 IO 性能CD-HIT 在运行过程中会创建大量临时文件。如果你的/tmp目录是普通的机械硬盘而你的项目目录在 SSD 上可以考虑把临时目录指向更快的位置export TMPDIR/path/to/ssd/tmp cd-hit -i input.fasta -o output -c 0.9虽然这个操作看起来微不足道但在处理百万条序列时IO 性能的差距会被放大得非常明显。另外输出文件最好写到项目本地避免写到网络文件系统上。8. 实战示例从原始 FASTQ 到非冗余蛋白序列的完整流水线理论说了那么多不如直接看一个完整的实战案例。这里我用自己做过的一个宏基因组项目为例展示从原始测序 reads 到非冗余蛋白序列目录的完整流水线。8.1 案例背景与分析设计这个项目的目标是构建一个人类肠道宏基因组基因目录。原始数据来自 50 个样本的宏基因组测序每个样本约 10 Gb。分析流程是质控 宿主序列过滤。样本单独组装。对所有样本的 contig 做基因预测。将所有预测得到的蛋白序列合并成一个文件。用 CD-HIT 去冗余得到非冗余基因目录。整个流程最耗时的步骤就是第 5 步。合并后的蛋白序列文件大概有 600 万条序列如果不做策略优化直接丢给 CD-HIT 很可能会跑半天甚至直接内存爆掉。8.2 CD-HIT 命令与参数选择我先做了一步预处理# 将所有样本的蛋白序列合并 cat sample*/prodigal/*.faa all_orfs.faa # 清除非法字符和过短序列 seqkit seq -g -w 0 all_orfs.faa all_orfs_clean.faa seqkit seq -m 30 -M 5000 all_orfs_clean.faa all_orfs_filtered.faa然后用 CD-HIT 做两轮聚类# 第一轮粗聚类90% 相似度 cd-hit -i all_orfs_filtered.faa -o gene_catalog_c90 -c 0.9 -n 5 -T 8 -M 8000 -d 100 # 提取代表序列 # 第二轮细聚类95% 相似度 cd-hit -i gene_catalog_c90 -o gene_catalog_c95 -c 0.95 -n 5 -T 8 -M 8000 -d 100两轮跑完后最终得到约 180 万条代表蛋白序列。这个规模跟类似研究的基因目录大小一致说明聚类效果是合理的。8.3 分析结果的质控与可视化拿到结果后我通常会做以下几步质控统计聚类前后序列数量计算去冗余率。检查.clstr文件里每个类的序列数分布确认没有异常巨大的类。随机抽样一些类人工检查代表序列是否合理。用 eggNOG-mapper 对代表序列做功能注释评估功能覆盖度是否足够。其中第 2 点很关键。如果出现一个类包含了几万条序列的情况说明这个类可能是由重复序列或低复杂度区域导致的需要特别关注。9. 最后聊几句实操心得CD-HIT 这个工具最让我满意的一点是它的稳定性和可预测性。不像一些新潮的工具版本更新频繁算法动不动就大变样CD-HIT 的架构非常稳定很多年前写的流程现在依然能流畅运行这种可维护性在生信工具里相当难得。我手头就有几条两年前跑的流程现在原封不动拿出来重跑结果还是一模一样。踩过几次坑之后我的体会是CD-HIT 使用最大的门槛不是工具本身而是你得清楚自己的数据适合什么阈值、什么参数。多试几次不同的-c和-n组合观察聚类结果的变化这是最快的学习路径。工具本身不会告诉你阈值设多少合适只有对数据的理解才是最终的判断依据。如果你刚开始接触 CD-HIT建议先从一个小规模的测试集入手比如选 1 万条序列分别跑-c 0.8、-c 0.9、-c 0.95三个阈值对比聚类结果的差异感受一下参数对结果的影响。这个过程花不了多少时间但对后面处理大规模数据非常有帮助。最后再分享一个小技巧CD-HIT 的.clstr文件一定要保留好不要因为占空间就随手删掉。很多时候你后面需要回溯某个代表序列到底包含了哪些原始序列这时候.clstr就是唯一线索。我自己会在项目目录里专门建一个clustering/子目录把输入、输出、参数和日志文件都整理好这样即使隔了很久再回头看项目也能快速搞清楚当时的分析逻辑。