CD-HIT 这个工具在生物信息学里属于那种表面上看起来平平无奇、但实际上几乎所有组学分析流程都会默认调用一次的隐形基础件。我第一次真正重视它是在处理一套全基因组重测序数据时BLAST 把所有重复区域和旁系同源基因都原封不动地保留下来下游注释直接跑了一周没出来最后才发现冗余序列导致的组合爆炸问题而 CD-HIT 一条命令就解决了 80% 的麻烦。如果你是做基因组注释、宏基因组、16S 扩增子或者蛋白家族分析的这篇文章就是给你准备的。我会从序列去冗余这个核心诉求出发把 CD-HIT 的原理、参数、输出文件、性能调优到实际工作流里的那些坑一次说清楚。1. 序列去重这件事为什么跑不掉 CD-HIT1.1 它到底解决的是哪一类问题CD-HIT 的核心功能说白了就是一句话给定一批序列按照设定的相似度阈值把那些长得差不多的序列合并成一个簇每条簇保留一个代表序列。这个问题听起来很简单但实际场景里到处都是。比如你做宏基因组组装从样本里拼出来几万条重叠群里面大量序列其实来自同一物种的相同基因只是测序误差和组装断点造成了细微差别又比如你用 BLAST 做蛋白功能注释数据库里一堆直系同源和旁系同源序列如果不先去冗余比对时间和结果噪音都会失控再比如构建参考基因集时来自不同样本的同一个基因序列有 SNP 差异但功能完全一致保留全部信息既浪费存储又让下游分析没法收敛。这些场景共同的需求不是单纯删掉完全相同的序列而是按相似度做聚类把相似但不等同的序列归到一组。这正是 CD-HIT 的强项。1.2 冗余不是垃圾但会让结果失真很多刚入门的同学会有个误解认为冗余序列只是多占点硬盘去不去都能接受。这个想法在数据量小的时候可能勉强成立一旦数据量上来后果是连锁反应先看存储和计算资源。一套普通规模的宏基因组样本去冗余前可能几百 MB 的 FASTABLAST 比对时复杂度是 O(N²) 级别的序列量翻倍比对时间翻四倍。你大概率遇到过比对挂了一晚上还没跑完的情况别怀疑冗余序列十有八九是元凶。再看生物学结论。做基因丰度统计时如果同一个基因在样本 A 里被测序组装出了 10 条几乎一样的 contig在样本 B 里只有 2 条你不去冗余就直接比较丰度得出的差异是完全虚假的。反过来做系统发育树构建冗余序列会让树的拓扑结构被高丰度重复序列绑架单系群的判断很容易出错。所以这里的冗余并不是说这些序列是垃圾而是它们在统计和比较层面造成了信息重复计算。CD-HIT 的聚类恰好是在保留代表信息的条件下把重复计算消除掉这就是它存在于所有标准化流程里的根本原因。1.3 CD-HIT 和 BLAST/MMseqs2 这类工具到底有什么差别很多人会问我为什么不用 BLAST 自己写脚本做聚类或者用 MMseqs2 不是更快吗BLAST 本身是全对全比对工具它的设计目标是找到两两之间的最优匹配输出包含大量匹配关系和打分信息。如果拿它做聚类你需要自己解析比对矩阵、设计聚类策略而且全对全比对的资源开销在序列量上万以后就变得不现实了。CD-HIT 的思路完全不同它采用贪心增量聚类依次读取序列和已有的每个簇的代表序列比对如果相似度达到阈值就归入该簇否则新建簇。这样全程只和代表序列比对不需要全对全矩阵时间复杂度和空间占用都被大幅压缩。MMseqs2 是另一个高性能方案它用了 k-mer 索引和预过滤速度确实比 CD-HIT 快但它的定位更偏快速同源搜索参数体系相对复杂新手上手成本高。CD-HIT 在做聚类去冗余这个特定任务上参数简单直接输出格式直观对上游分析来说已经是足够好的选择。实际项目中我更倾向于 CD-HIT 做常规去冗余只有在序列量极大百万级且对时间极度敏感时才会考虑 MMseqs2 的 linclust 模式。2. 安装与跑通其实很简单但有几个容易忽略的坑2.1 conda 安装是首选源码编译留作备选CD-HIT 的安装方式有两种主流路径我个人的建议是能用 conda 就用 conda。conda install -c bioconda cd-hit一条命令搞定依赖关系主要是 GCC 运行时和 zlib全部自动处理版本也保持稳定。以我现在的使用频率conda 装的版本长期固定在 4.8.1跑各类流程都没出过问题。如果你所在的服务器没有外网或者环境隔离要求严格那就走源码编译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之后二进制文件是生成在当前目录下的不会自动安装到系统路径。很多新手在这里直接敲cd-hit发现 command not found然后就懵了。你需要手动把编译出的二进制文件加进 PATHexport PATH$PWD:$PATH想永久生效的话把这一行写进~/.bashrc或者~/.profile。2.2 检查安装是否成功安装完以后别急着跑正式数据先做一个最小验证cd-hit --help如果能看到一长串参数说明说明基础安装没问题。另外验一下并行版本cd-hit -T 4 --help能正常执行就说明多线程支持没问题。有时候 conda 装的版本因为 OpenMP 编译配置缺失-T参数会被静默忽略这种情况在后续处理大数据时很头疼提前验证能省很多事。2.3 最小示例跑通建立信心我习惯用一个很小的测试 FASTA 来验证全流程。假设test.fa里有 5 条序列seq1 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQ seq2 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQ seq3 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVA seq4 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVE seq5 MKTGYIAKQRQISFVKSHFSRQLEERLGLIEVQ运行cd-hit -i test.fa -o test_out -c 0.9 -n 5如果一切正常你会得到test_out代表序列 FASTA和test_out.clstr聚类信息再跑一个cat test_out看输出内容。看到聚合成几条代表序列这个工具你就已经用得起来了。别小看这一步我遇到过不下三次因为环境问题导致程序貌似正常但输出为空的情况提前用 5 条序列验证能快速暴露环境缺陷。3. 核心参数逐个拆解每个参数背后都是算法逻辑3.1 -c 相似度阈值整个工具的灵魂-c指定的是聚类所用的序列相似度阈值默认值是 0.9表示两条序列的一致性identity达到 90% 就归为一类。这个值的设定完全取决于你的下游分析目的做蛋白家族去冗余或去除近同源序列通常用-c 0.9甚至-c 0.95保证簇内序列高度一致做属或科级别的宏基因组分类有时会用-c 0.8甚至更低允许更多序列汇聚到一族研究直系同源基因聚类时我见过有人用到-c 0.7但说实话低于 0.7 以后 CD-HIT 本身的灵敏度会下降因为 k-mer 表的设计就没法覆盖太远的同源关系。不用纠结哪个值是标准答案你需要理解的是阈值越高簇越紧代表序列越能代表簇里所有成员阈值越低聚类越宽松去冗余力度越大但簇内的异质性也跟着上升。3.2 -n 的 k-mer 选择一个经常被忽略但直接决定结果的参数-n指定 k-mer 的长度它和-c是绑定的CD-HIT 官方文档给了一个参考关系相似度阈值蛋白序列推荐 -n核酸序列CD-HIT-EST推荐 -n0.7510 或 110.84 或 59 或 100.93 或 48 或 91.038 或 9这里要解释一下机制CD-HIT 不是直接对每条序列做全比对来确定聚类关系而是先通过短 k-mer 片段快速过滤只有共享足够多 k-mer 的序列才会进入下一步真正的比对。-n设得越大k-mer 越特异但能匹配上的同源序列就越少容易漏-n设得越小k-mer 越普通灵敏度高但计算量变大误报也变多。我自己的经验是蛋白序列用-c 0.9时-n 5是一个安全选择既不会漏掉真实的同源簇计算速度也让人满意。有人会想用-n 2来提高灵敏度结果运行时间翻了好几倍聚类结果却差不多完全不划算。3.3 -g 参数控制最先收录还是最适收录-g默认值是 0含义是谁先出现谁当代表序列。这对于数据的顺序非常敏感如果你把高质量参考序列放在前面它们更可能成为代表序列如果数据顺序是混乱的代表序列的选择就变成随机的。把-g 1打开以后CD-HIT 会优先寻找和簇内其他成员相似度最高的那条作为代表序列代价是计算量会上升不少。我在做宏基因组非冗余基因集时强烈建议开-g 1否则代表序列可能是某条半残的 contig导致后续功能注释出现偏差。但要注意-g 1只能用于蛋白序列模式CD-HIT核酸模式CD-HIT-EST不支持这个参数我第一次在 EST 模式里试着开-g 1工具直接报错退出查了文档才发现这个限制。3.4 -aS 和 -aL长度覆盖率的隐形约束-c只衡量匹配区域的相似度但它不关心匹配区域占了整条序列的多大比例。举个例子两条长度差 5 倍的序列短序列和长序列的一部分达到 95% 相似按-c 0.9的标准它们会聚成一类但长序列有很大一段是短序列没有的把它们强行归到一类在生物学上常常不合理。-aS表示短序列长度必须覆盖到其自身长度的比例-aL表示短序列长度占长序列长度的比例。比如-aS 0.8意思是短序列 80% 以上的长度都得比对得上才能聚类。做全长基因聚类时我会额外开-aL 0.8或-aS 0.8避免把部分缺失序列和全长序列混在一起。3.5 -T 和 -M线程和内存的调度这两个参数在序列量大的时候直接影响你能不能跑完任务。-T指定线程数-M指定最大可用内存以 MB 为单位。cd-hit -i input.fa -o output -c 0.9 -n 5 -T 8 -M 16000上面的命令让 CD-HIT 使用 8 个线程最多 16 GB 内存。这里有个容易被忽略的点CD-HIT 的写入阶段输出排序不是完全多线程的你开再多的线程最后一步可能还是单线程跑所以当你看到 CPU 占用率突然掉下来时不要惊慌它可能正在写文件。4. 蛋白和核酸模式的区别CD-HIT 和 CD-HIT-EST4.1 算法层面的差异决定了它们的适用边界CD-HIT 最早设计是基于氨基酸序列的它的 k-mer 统计和比对策略天然适配蛋白质序列的 20 种字母表。后来才加入了对核酸序列的支持也就是cd-hit-est。两者的核心算法差异在于蛋白质序列的 k-mer 是 3~5 的长度而核酸序列因为只有 4 种碱基同样的 k-mer 长度随机匹配概率高出很多所以核酸需要更长的 k-mer来保证特异性。这就是为什么上表里核酸的-n推荐值明显更大。另外核酸序列还有正负链问题CD-HIT-EST 在比对时会把反向互补链和正链同等对待确保来自基因两侧不同方向的序列也能被正确定义聚类关系。如果用 CD-HIT 直接处理核酸序列不考虑反向互补很多实际同源的序列会被拆成两族。4.2 参数体系完全不同用核酸模式时命令名是cd-hit-est而不是cd-hit参数表也变了-n推荐值在 8~11 之间同时没有-g参数前面提过-aS仍然有效但-aL的行为略有差异。cd-hit-est -i contigs.fa -o contigs_nr -c 0.95 -n 10 -aS 0.9 -T 8 -M 16000这个命令是我处理细菌基因组 contig 去重时经常用的。需要特别说明cd-hit-est支持的阈值-c最低是 0.8低于这个值不是不能用而是聚类的可靠性急剧下降因为核酸序列本身随机相似的噪音会被误判为同源。如果你非得对核酸做 0.7 这种低阈值聚类建议先用cd-hit把序列翻译成蛋白再做别硬用 EST 模式。4.3 具体场景怎么选这里直接给出我自己的选择习惯场景用哪个工具典型参数蛋白序列去除冗余cd-hit-c 0.9 -n 5 -g 1contig 序列去冗余cd-hit-est-c 0.95 -n 10 -aS 0.9非编码 RNA 序列聚类cd-hit-est-c 0.8 -n 11宏基因组非冗余基因集cd-hit翻译后蛋白-c 0.95 -n 5 -g 1 -aS 0.916S/ITS 扩增子 ASV 聚类cd-hit-est-c 0.97 -n 10但更推荐 vsearch 的 UPARSE 流程后面说注意最后一行16S 扩增子的 OTU 聚类虽然技术上 CD-HIT 能做但它的聚类策略和 UCLUST/UNOISE 这类专门算法在生物学定义上有差异现在主流的扩增子流程更推荐 vsearch 或 QIIME2 内置的聚类方式。CD-HIT 在这个领域的角色更多是作为参考去重工具而不是 OTU 定义工具。5. 一个标准的实战流程从原始 FASTA 到非冗余序列集5.1 输入文件预处理这一步省不了CD-HIT 对输入格式要求相对宽松但有几个规则必须遵守第一FASTA 的序列名不能有重复否则程序在输出时可能会产生冲突。尤其当你从多个样本合并 FASTA 时不同样本里同样的序列名会导致聚类结果混乱。一个稳妥的做法是在合并前给序列名加样本前缀sed s/^/SampleA_/ SampleA.fa SampleA_prefix.fa第二序列长度过滤要在 CD-HIT 之前做好。CD-HIT 本身没有过滤短序列的参数垃圾短序列会直接参与聚类白白浪费 k-mer 表空间。用 seqkit 或 awk 预处理一下seqkit seq -m 100 -M 5000 input.fa filtered.fa5.2 核心命令这样写信息量最大我以一个构建宏基因组非冗余基因集的真实场景为例。假设你已经从多个样本的组装结果中预测出了蛋白序列文件是all_proteins.faa有 80 万条序列cd-hit -i all_proteins.faa -o nr_proteins -c 0.9 -n 5 -g 1 -aS 0.9 -T 16 -M 64000拆解一下-i输入文件-o输出文件前缀实际会生成nr_proteins和nr_proteins.clstr-c 0.9相似度 90%-n 5对应 0.9 阈值的推荐 k-mer-g 1选择最优代表序列-aS 0.9短序列 90% 以上长度要覆盖到-T 16 -M 6400016 线程64 GB 内存跑完以后统计一下grep -c ^ nr_proteins这个数字就是去冗余后的非冗余基因数量。通常原始 80 万条可以去掉 30%~50%具体取决于样本来源和阈值。5.3 验证聚类结果的合理性聚类跑完以后还需要做质量验证。我通常看三个指标代表序列数量是否符合预期和样本多样性、数据库规模对照每条序列的聚类簇大小分布是否合理是否存在单序列簇过多的情况随机抽几条簇内序列做多序列比对确认它们确实满足设定的相似度阈值如果单序列簇占了 80% 以上说明阈值设得太严或者数据本来就高度多样如果所有序列都汇成少数几个巨簇说明阈值过松。没有绝对的正确标准但至少要和你的生物学预期对得上。6. 输出文件的读取和进阶玩法6.1 .clstr 文件的信息量拿到nr_proteins.clstr以后你会看到类似这样的内容Cluster 0 0 115aa, seq1... at 90.43% 1 112aa, seq2... at 91.07% Cluster 1 0 120aa, seq3... at 100.00%每行开头的0、1表示该簇内序列编号115aa是序列长度seq1是该条序列的名字at 90.43%表示该序列在比对区和代表序列的一致度。簇编号不代表任何生物学含义完全取决于序列输入顺序和聚类顺序。处理这个文件时最常用的操作是从.clstr里提取所有代表序列的名字或者每个簇的成员列表。用一个小脚本就能搞定awk /^Cluster/ {cluster$2} /\*/ {print cluster, $3} nr_proteins.clstr | head不过更优雅的做法是直接用序列重命名的方式处理把每个簇的编号映射回代表序列。6.2 提取特定簇序列做下游分析有时候你只关心特定簇比如某个感兴趣的蛋白家族。可以先用 grep 找到簇编号再用序列名从 FASTA 中提取。借助 seqkit 和自写一个小 while 循环就能做。实际场景里我建议写一个简单的 Python 脚本把.clstr解析成 TSV 格式一劳永逸import re clstr_file nr_proteins.clstr with open(clstr_file) as f: current_cluster for line in f: if line.startswith(): current_cluster line.strip().replace(Cluster , ) else: seqname re.search(r(\S), line).group(1) identity re.search(rat (\d\.\d)%, line).group(1) print(f{current_cluster}\t{seqname}\t{identity})这个 TSV 文件可以直接导入 R 或 Python 做后续的簇规模统计、物种丰度矩阵构建。6.3 把 CD-HIT 嵌入到 Snakemake/Nextflow 流程里一旦涉及几十个样本的批量处理手敲命令就不现实了。我在 Snakemake 里的 rule 通常是这样写的rule cdhit: input: data/{sample}.faa output: results/{sample}_nr.faa, results/{sample}_nr.clstr threads: 16 params: c0.9, n5, g1, as0.9 shell: cd-hit -i {input} -o results/{sample}_nr -c {params.c} -n {params.n} -g {params.g} -aS {params.as} -T {threads} -M 64000注意输出文件名的前缀参数传递CD-HIT 的-o决定了两个输出文件的命名-o传入results/{sample}_nr输出就是_nr和_nr.clstr这个逻辑在 workflow 里常见但容易写错。7. 我踩过的坑和性能调优的实战经验7.1 大数据量下内存不够CD-HIT 的内存消耗和输入序列的 k-mer 种类数直接相关序列多样性越高内存占用越大。零散的小序列会急剧增加聚类索引的大小有时候 100 万条蛋白序列就能吃掉 30 GB 内存。解决方案有两个方向第一个方向是从源头减少序列数量比如先用 0.95 的高阈值快速去一遍重再用 0.9 做正式的宽松聚类。两轮 CD-HIT 跑出来的结果通常和直接一次 0.9 聚类差异很小但峰值内存能下降 40%。第二个方向是适当调低-M参数但这里有个权衡-M设置过低时 CD-HIT 会频繁读写临时文件时间成本上升非常明显。我的经验是直接把-M设为服务器物理内存的 80%让操作系统去调度效果最好。7.2 输入序列的污染对聚类结果的影响这是个隐藏很深的坑。如果你的 FASTA 里混有大量不完全蛋白序列比如从基因组预测时把假基因也算上了CD-HIT 会忠实地把它们当作独立的序列成员参与聚类。由于假基因序列通常比正常序列短-aS参数设置不合理时它们会躲过覆盖率的约束被错误地并入真实基因的簇。解决办法是在进入 CD-HIT 之前对输入文件做一轮严格的长度和完整性过滤awk /^/ {if (seq) print seq; seq; print; next} {seqseq$0} END {print seq} input.faa | seqkit seq -m 50 --max-len 5000 clean.faa另外帧移错误frameshift导致的蛋白序列截断是宏基因组注释里最常见的污染源这类序列在用 CD-HIT 聚类后很可能形成大量伪簇。如果发现聚类结果里单序列簇异常多强烈建议回顾上游的基因预测参数。7.3 并行度不是越高越好CD-HIT 的多线程加速效果在 4~16 线程之间提升明显但超过 16 线程后加速比会急剧下降。这和它的算法设计有关预过滤阶段可以并行但聚类的核心决策是串行依赖的——后一个簇是否成立依赖前面所有簇的代表序列做比对。我实际测过一个 50 万条序列的数据集8 线程耗时 25 分钟16 线程 18 分钟32 线程反而 20 分钟。所以盲目堆线程不仅浪费资源还可能是负优化。7.4 聚类统计的误用CD-HIT 的簇不是 OTU有些做微生物多样性的同学试图直接用 CD-HIT 的聚类结果当作 OTU 表做多样性分析。这里我必须泼一盆冷水CD-HIT 作为一个通用聚类工具它的聚类边界没有针对 16S 序列进行物种水平的生物学校准直接用默认的 97% 相似度聚出来的 OTU跟 QIIME2 的 VSEARCH 或 UNOISE 聚出来的 OTU 在数目上可能差出 20% 以上。如果你想做扩增子 OTU 聚类还是老老实实走专门的扩增子分析流程。CD-HIT 更合适的定位是高质量参考数据集的构建工具比如构建非冗余参考基因集、清理比对前的同源序列、压缩 blast 数据库等。最后再分享一个我自己的小习惯跑完 CD-HIT 之后我永远会先看一眼代表序列的数量再随机抽 20 条代表序列回比对到原始数据集里做反向验证。这一步成本很低但能发现 90% 的参数设置错误。序列聚类这类任务往往是数据里的单碱基变异、测序错误加上参数边界效应叠加在一起光看命令行没问题是不够的得用真实数据说话。这套预处理 → 参数合理性 → 输出验证的思路放到任何序列聚类工具上都适用。 CD-HIT序列去冗余详解：原理、参数与实战应用
CD-HIT 这个工具在生物信息学里属于那种表面上看起来平平无奇、但实际上几乎所有组学分析流程都会默认调用一次的隐形基础件。我第一次真正重视它是在处理一套全基因组重测序数据时BLAST 把所有重复区域和旁系同源基因都原封不动地保留下来下游注释直接跑了一周没出来最后才发现冗余序列导致的组合爆炸问题而 CD-HIT 一条命令就解决了 80% 的麻烦。如果你是做基因组注释、宏基因组、16S 扩增子或者蛋白家族分析的这篇文章就是给你准备的。我会从序列去冗余这个核心诉求出发把 CD-HIT 的原理、参数、输出文件、性能调优到实际工作流里的那些坑一次说清楚。1. 序列去重这件事为什么跑不掉 CD-HIT1.1 它到底解决的是哪一类问题CD-HIT 的核心功能说白了就是一句话给定一批序列按照设定的相似度阈值把那些长得差不多的序列合并成一个簇每条簇保留一个代表序列。这个问题听起来很简单但实际场景里到处都是。比如你做宏基因组组装从样本里拼出来几万条重叠群里面大量序列其实来自同一物种的相同基因只是测序误差和组装断点造成了细微差别又比如你用 BLAST 做蛋白功能注释数据库里一堆直系同源和旁系同源序列如果不先去冗余比对时间和结果噪音都会失控再比如构建参考基因集时来自不同样本的同一个基因序列有 SNP 差异但功能完全一致保留全部信息既浪费存储又让下游分析没法收敛。这些场景共同的需求不是单纯删掉完全相同的序列而是按相似度做聚类把相似但不等同的序列归到一组。这正是 CD-HIT 的强项。1.2 冗余不是垃圾但会让结果失真很多刚入门的同学会有个误解认为冗余序列只是多占点硬盘去不去都能接受。这个想法在数据量小的时候可能勉强成立一旦数据量上来后果是连锁反应先看存储和计算资源。一套普通规模的宏基因组样本去冗余前可能几百 MB 的 FASTABLAST 比对时复杂度是 O(N²) 级别的序列量翻倍比对时间翻四倍。你大概率遇到过比对挂了一晚上还没跑完的情况别怀疑冗余序列十有八九是元凶。再看生物学结论。做基因丰度统计时如果同一个基因在样本 A 里被测序组装出了 10 条几乎一样的 contig在样本 B 里只有 2 条你不去冗余就直接比较丰度得出的差异是完全虚假的。反过来做系统发育树构建冗余序列会让树的拓扑结构被高丰度重复序列绑架单系群的判断很容易出错。所以这里的冗余并不是说这些序列是垃圾而是它们在统计和比较层面造成了信息重复计算。CD-HIT 的聚类恰好是在保留代表信息的条件下把重复计算消除掉这就是它存在于所有标准化流程里的根本原因。1.3 CD-HIT 和 BLAST/MMseqs2 这类工具到底有什么差别很多人会问我为什么不用 BLAST 自己写脚本做聚类或者用 MMseqs2 不是更快吗BLAST 本身是全对全比对工具它的设计目标是找到两两之间的最优匹配输出包含大量匹配关系和打分信息。如果拿它做聚类你需要自己解析比对矩阵、设计聚类策略而且全对全比对的资源开销在序列量上万以后就变得不现实了。CD-HIT 的思路完全不同它采用贪心增量聚类依次读取序列和已有的每个簇的代表序列比对如果相似度达到阈值就归入该簇否则新建簇。这样全程只和代表序列比对不需要全对全矩阵时间复杂度和空间占用都被大幅压缩。MMseqs2 是另一个高性能方案它用了 k-mer 索引和预过滤速度确实比 CD-HIT 快但它的定位更偏快速同源搜索参数体系相对复杂新手上手成本高。CD-HIT 在做聚类去冗余这个特定任务上参数简单直接输出格式直观对上游分析来说已经是足够好的选择。实际项目中我更倾向于 CD-HIT 做常规去冗余只有在序列量极大百万级且对时间极度敏感时才会考虑 MMseqs2 的 linclust 模式。2. 安装与跑通其实很简单但有几个容易忽略的坑2.1 conda 安装是首选源码编译留作备选CD-HIT 的安装方式有两种主流路径我个人的建议是能用 conda 就用 conda。conda install -c bioconda cd-hit一条命令搞定依赖关系主要是 GCC 运行时和 zlib全部自动处理版本也保持稳定。以我现在的使用频率conda 装的版本长期固定在 4.8.1跑各类流程都没出过问题。如果你所在的服务器没有外网或者环境隔离要求严格那就走源码编译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之后二进制文件是生成在当前目录下的不会自动安装到系统路径。很多新手在这里直接敲cd-hit发现 command not found然后就懵了。你需要手动把编译出的二进制文件加进 PATHexport PATH$PWD:$PATH想永久生效的话把这一行写进~/.bashrc或者~/.profile。2.2 检查安装是否成功安装完以后别急着跑正式数据先做一个最小验证cd-hit --help如果能看到一长串参数说明说明基础安装没问题。另外验一下并行版本cd-hit -T 4 --help能正常执行就说明多线程支持没问题。有时候 conda 装的版本因为 OpenMP 编译配置缺失-T参数会被静默忽略这种情况在后续处理大数据时很头疼提前验证能省很多事。2.3 最小示例跑通建立信心我习惯用一个很小的测试 FASTA 来验证全流程。假设test.fa里有 5 条序列seq1 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQ seq2 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVQ seq3 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVA seq4 MKTAYIAKQRQISFVKSHFSRQLEERLGLIEVE seq5 MKTGYIAKQRQISFVKSHFSRQLEERLGLIEVQ运行cd-hit -i test.fa -o test_out -c 0.9 -n 5如果一切正常你会得到test_out代表序列 FASTA和test_out.clstr聚类信息再跑一个cat test_out看输出内容。看到聚合成几条代表序列这个工具你就已经用得起来了。别小看这一步我遇到过不下三次因为环境问题导致程序貌似正常但输出为空的情况提前用 5 条序列验证能快速暴露环境缺陷。3. 核心参数逐个拆解每个参数背后都是算法逻辑3.1 -c 相似度阈值整个工具的灵魂-c指定的是聚类所用的序列相似度阈值默认值是 0.9表示两条序列的一致性identity达到 90% 就归为一类。这个值的设定完全取决于你的下游分析目的做蛋白家族去冗余或去除近同源序列通常用-c 0.9甚至-c 0.95保证簇内序列高度一致做属或科级别的宏基因组分类有时会用-c 0.8甚至更低允许更多序列汇聚到一族研究直系同源基因聚类时我见过有人用到-c 0.7但说实话低于 0.7 以后 CD-HIT 本身的灵敏度会下降因为 k-mer 表的设计就没法覆盖太远的同源关系。不用纠结哪个值是标准答案你需要理解的是阈值越高簇越紧代表序列越能代表簇里所有成员阈值越低聚类越宽松去冗余力度越大但簇内的异质性也跟着上升。3.2 -n 的 k-mer 选择一个经常被忽略但直接决定结果的参数-n指定 k-mer 的长度它和-c是绑定的CD-HIT 官方文档给了一个参考关系相似度阈值蛋白序列推荐 -n核酸序列CD-HIT-EST推荐 -n0.7510 或 110.84 或 59 或 100.93 或 48 或 91.038 或 9这里要解释一下机制CD-HIT 不是直接对每条序列做全比对来确定聚类关系而是先通过短 k-mer 片段快速过滤只有共享足够多 k-mer 的序列才会进入下一步真正的比对。-n设得越大k-mer 越特异但能匹配上的同源序列就越少容易漏-n设得越小k-mer 越普通灵敏度高但计算量变大误报也变多。我自己的经验是蛋白序列用-c 0.9时-n 5是一个安全选择既不会漏掉真实的同源簇计算速度也让人满意。有人会想用-n 2来提高灵敏度结果运行时间翻了好几倍聚类结果却差不多完全不划算。3.3 -g 参数控制最先收录还是最适收录-g默认值是 0含义是谁先出现谁当代表序列。这对于数据的顺序非常敏感如果你把高质量参考序列放在前面它们更可能成为代表序列如果数据顺序是混乱的代表序列的选择就变成随机的。把-g 1打开以后CD-HIT 会优先寻找和簇内其他成员相似度最高的那条作为代表序列代价是计算量会上升不少。我在做宏基因组非冗余基因集时强烈建议开-g 1否则代表序列可能是某条半残的 contig导致后续功能注释出现偏差。但要注意-g 1只能用于蛋白序列模式CD-HIT核酸模式CD-HIT-EST不支持这个参数我第一次在 EST 模式里试着开-g 1工具直接报错退出查了文档才发现这个限制。3.4 -aS 和 -aL长度覆盖率的隐形约束-c只衡量匹配区域的相似度但它不关心匹配区域占了整条序列的多大比例。举个例子两条长度差 5 倍的序列短序列和长序列的一部分达到 95% 相似按-c 0.9的标准它们会聚成一类但长序列有很大一段是短序列没有的把它们强行归到一类在生物学上常常不合理。-aS表示短序列长度必须覆盖到其自身长度的比例-aL表示短序列长度占长序列长度的比例。比如-aS 0.8意思是短序列 80% 以上的长度都得比对得上才能聚类。做全长基因聚类时我会额外开-aL 0.8或-aS 0.8避免把部分缺失序列和全长序列混在一起。3.5 -T 和 -M线程和内存的调度这两个参数在序列量大的时候直接影响你能不能跑完任务。-T指定线程数-M指定最大可用内存以 MB 为单位。cd-hit -i input.fa -o output -c 0.9 -n 5 -T 8 -M 16000上面的命令让 CD-HIT 使用 8 个线程最多 16 GB 内存。这里有个容易被忽略的点CD-HIT 的写入阶段输出排序不是完全多线程的你开再多的线程最后一步可能还是单线程跑所以当你看到 CPU 占用率突然掉下来时不要惊慌它可能正在写文件。4. 蛋白和核酸模式的区别CD-HIT 和 CD-HIT-EST4.1 算法层面的差异决定了它们的适用边界CD-HIT 最早设计是基于氨基酸序列的它的 k-mer 统计和比对策略天然适配蛋白质序列的 20 种字母表。后来才加入了对核酸序列的支持也就是cd-hit-est。两者的核心算法差异在于蛋白质序列的 k-mer 是 3~5 的长度而核酸序列因为只有 4 种碱基同样的 k-mer 长度随机匹配概率高出很多所以核酸需要更长的 k-mer来保证特异性。这就是为什么上表里核酸的-n推荐值明显更大。另外核酸序列还有正负链问题CD-HIT-EST 在比对时会把反向互补链和正链同等对待确保来自基因两侧不同方向的序列也能被正确定义聚类关系。如果用 CD-HIT 直接处理核酸序列不考虑反向互补很多实际同源的序列会被拆成两族。4.2 参数体系完全不同用核酸模式时命令名是cd-hit-est而不是cd-hit参数表也变了-n推荐值在 8~11 之间同时没有-g参数前面提过-aS仍然有效但-aL的行为略有差异。cd-hit-est -i contigs.fa -o contigs_nr -c 0.95 -n 10 -aS 0.9 -T 8 -M 16000这个命令是我处理细菌基因组 contig 去重时经常用的。需要特别说明cd-hit-est支持的阈值-c最低是 0.8低于这个值不是不能用而是聚类的可靠性急剧下降因为核酸序列本身随机相似的噪音会被误判为同源。如果你非得对核酸做 0.7 这种低阈值聚类建议先用cd-hit把序列翻译成蛋白再做别硬用 EST 模式。4.3 具体场景怎么选这里直接给出我自己的选择习惯场景用哪个工具典型参数蛋白序列去除冗余cd-hit-c 0.9 -n 5 -g 1contig 序列去冗余cd-hit-est-c 0.95 -n 10 -aS 0.9非编码 RNA 序列聚类cd-hit-est-c 0.8 -n 11宏基因组非冗余基因集cd-hit翻译后蛋白-c 0.95 -n 5 -g 1 -aS 0.916S/ITS 扩增子 ASV 聚类cd-hit-est-c 0.97 -n 10但更推荐 vsearch 的 UPARSE 流程后面说注意最后一行16S 扩增子的 OTU 聚类虽然技术上 CD-HIT 能做但它的聚类策略和 UCLUST/UNOISE 这类专门算法在生物学定义上有差异现在主流的扩增子流程更推荐 vsearch 或 QIIME2 内置的聚类方式。CD-HIT 在这个领域的角色更多是作为参考去重工具而不是 OTU 定义工具。5. 一个标准的实战流程从原始 FASTA 到非冗余序列集5.1 输入文件预处理这一步省不了CD-HIT 对输入格式要求相对宽松但有几个规则必须遵守第一FASTA 的序列名不能有重复否则程序在输出时可能会产生冲突。尤其当你从多个样本合并 FASTA 时不同样本里同样的序列名会导致聚类结果混乱。一个稳妥的做法是在合并前给序列名加样本前缀sed s/^/SampleA_/ SampleA.fa SampleA_prefix.fa第二序列长度过滤要在 CD-HIT 之前做好。CD-HIT 本身没有过滤短序列的参数垃圾短序列会直接参与聚类白白浪费 k-mer 表空间。用 seqkit 或 awk 预处理一下seqkit seq -m 100 -M 5000 input.fa filtered.fa5.2 核心命令这样写信息量最大我以一个构建宏基因组非冗余基因集的真实场景为例。假设你已经从多个样本的组装结果中预测出了蛋白序列文件是all_proteins.faa有 80 万条序列cd-hit -i all_proteins.faa -o nr_proteins -c 0.9 -n 5 -g 1 -aS 0.9 -T 16 -M 64000拆解一下-i输入文件-o输出文件前缀实际会生成nr_proteins和nr_proteins.clstr-c 0.9相似度 90%-n 5对应 0.9 阈值的推荐 k-mer-g 1选择最优代表序列-aS 0.9短序列 90% 以上长度要覆盖到-T 16 -M 6400016 线程64 GB 内存跑完以后统计一下grep -c ^ nr_proteins这个数字就是去冗余后的非冗余基因数量。通常原始 80 万条可以去掉 30%~50%具体取决于样本来源和阈值。5.3 验证聚类结果的合理性聚类跑完以后还需要做质量验证。我通常看三个指标代表序列数量是否符合预期和样本多样性、数据库规模对照每条序列的聚类簇大小分布是否合理是否存在单序列簇过多的情况随机抽几条簇内序列做多序列比对确认它们确实满足设定的相似度阈值如果单序列簇占了 80% 以上说明阈值设得太严或者数据本来就高度多样如果所有序列都汇成少数几个巨簇说明阈值过松。没有绝对的正确标准但至少要和你的生物学预期对得上。6. 输出文件的读取和进阶玩法6.1 .clstr 文件的信息量拿到nr_proteins.clstr以后你会看到类似这样的内容Cluster 0 0 115aa, seq1... at 90.43% 1 112aa, seq2... at 91.07% Cluster 1 0 120aa, seq3... at 100.00%每行开头的0、1表示该簇内序列编号115aa是序列长度seq1是该条序列的名字at 90.43%表示该序列在比对区和代表序列的一致度。簇编号不代表任何生物学含义完全取决于序列输入顺序和聚类顺序。处理这个文件时最常用的操作是从.clstr里提取所有代表序列的名字或者每个簇的成员列表。用一个小脚本就能搞定awk /^Cluster/ {cluster$2} /\*/ {print cluster, $3} nr_proteins.clstr | head不过更优雅的做法是直接用序列重命名的方式处理把每个簇的编号映射回代表序列。6.2 提取特定簇序列做下游分析有时候你只关心特定簇比如某个感兴趣的蛋白家族。可以先用 grep 找到簇编号再用序列名从 FASTA 中提取。借助 seqkit 和自写一个小 while 循环就能做。实际场景里我建议写一个简单的 Python 脚本把.clstr解析成 TSV 格式一劳永逸import re clstr_file nr_proteins.clstr with open(clstr_file) as f: current_cluster for line in f: if line.startswith(): current_cluster line.strip().replace(Cluster , ) else: seqname re.search(r(\S), line).group(1) identity re.search(rat (\d\.\d)%, line).group(1) print(f{current_cluster}\t{seqname}\t{identity})这个 TSV 文件可以直接导入 R 或 Python 做后续的簇规模统计、物种丰度矩阵构建。6.3 把 CD-HIT 嵌入到 Snakemake/Nextflow 流程里一旦涉及几十个样本的批量处理手敲命令就不现实了。我在 Snakemake 里的 rule 通常是这样写的rule cdhit: input: data/{sample}.faa output: results/{sample}_nr.faa, results/{sample}_nr.clstr threads: 16 params: c0.9, n5, g1, as0.9 shell: cd-hit -i {input} -o results/{sample}_nr -c {params.c} -n {params.n} -g {params.g} -aS {params.as} -T {threads} -M 64000注意输出文件名的前缀参数传递CD-HIT 的-o决定了两个输出文件的命名-o传入results/{sample}_nr输出就是_nr和_nr.clstr这个逻辑在 workflow 里常见但容易写错。7. 我踩过的坑和性能调优的实战经验7.1 大数据量下内存不够CD-HIT 的内存消耗和输入序列的 k-mer 种类数直接相关序列多样性越高内存占用越大。零散的小序列会急剧增加聚类索引的大小有时候 100 万条蛋白序列就能吃掉 30 GB 内存。解决方案有两个方向第一个方向是从源头减少序列数量比如先用 0.95 的高阈值快速去一遍重再用 0.9 做正式的宽松聚类。两轮 CD-HIT 跑出来的结果通常和直接一次 0.9 聚类差异很小但峰值内存能下降 40%。第二个方向是适当调低-M参数但这里有个权衡-M设置过低时 CD-HIT 会频繁读写临时文件时间成本上升非常明显。我的经验是直接把-M设为服务器物理内存的 80%让操作系统去调度效果最好。7.2 输入序列的污染对聚类结果的影响这是个隐藏很深的坑。如果你的 FASTA 里混有大量不完全蛋白序列比如从基因组预测时把假基因也算上了CD-HIT 会忠实地把它们当作独立的序列成员参与聚类。由于假基因序列通常比正常序列短-aS参数设置不合理时它们会躲过覆盖率的约束被错误地并入真实基因的簇。解决办法是在进入 CD-HIT 之前对输入文件做一轮严格的长度和完整性过滤awk /^/ {if (seq) print seq; seq; print; next} {seqseq$0} END {print seq} input.faa | seqkit seq -m 50 --max-len 5000 clean.faa另外帧移错误frameshift导致的蛋白序列截断是宏基因组注释里最常见的污染源这类序列在用 CD-HIT 聚类后很可能形成大量伪簇。如果发现聚类结果里单序列簇异常多强烈建议回顾上游的基因预测参数。7.3 并行度不是越高越好CD-HIT 的多线程加速效果在 4~16 线程之间提升明显但超过 16 线程后加速比会急剧下降。这和它的算法设计有关预过滤阶段可以并行但聚类的核心决策是串行依赖的——后一个簇是否成立依赖前面所有簇的代表序列做比对。我实际测过一个 50 万条序列的数据集8 线程耗时 25 分钟16 线程 18 分钟32 线程反而 20 分钟。所以盲目堆线程不仅浪费资源还可能是负优化。7.4 聚类统计的误用CD-HIT 的簇不是 OTU有些做微生物多样性的同学试图直接用 CD-HIT 的聚类结果当作 OTU 表做多样性分析。这里我必须泼一盆冷水CD-HIT 作为一个通用聚类工具它的聚类边界没有针对 16S 序列进行物种水平的生物学校准直接用默认的 97% 相似度聚出来的 OTU跟 QIIME2 的 VSEARCH 或 UNOISE 聚出来的 OTU 在数目上可能差出 20% 以上。如果你想做扩增子 OTU 聚类还是老老实实走专门的扩增子分析流程。CD-HIT 更合适的定位是高质量参考数据集的构建工具比如构建非冗余参考基因集、清理比对前的同源序列、压缩 blast 数据库等。最后再分享一个我自己的小习惯跑完 CD-HIT 之后我永远会先看一眼代表序列的数量再随机抽 20 条代表序列回比对到原始数据集里做反向验证。这一步成本很低但能发现 90% 的参数设置错误。序列聚类这类任务往往是数据里的单碱基变异、测序错误加上参数边界效应叠加在一起光看命令行没问题是不够的得用真实数据说话。这套预处理 → 参数合理性 → 输出验证的思路放到任何序列聚类工具上都适用。