企业官网建设流程全解析

做生信的人，基本都有被CD-HIT救过命的经历。我最早接触它，是很多年前要处理一批从测序组装结果里预测出来的蛋白序列，差不多三万多条，里面有大量近乎一样的重复。当时想直接两两比对去冗余，结果发现这是O(n²)级别的计算，根本跑不动，后来换成CD-HIT，几分钟就搞定，代表序列一下压缩到几千条。

CD-HIT，全称是Cluster Database at High Identity with Tolerance，简单说就是一个序列聚类和去冗余工具。它的核心能力只有一个：把彼此相似度达到一定阈值的序列聚成一簇，每簇保留一条代表序列，最终输出一套“瘦身”之后还保留主要信息的数据集。这个能力听着简单，但几乎贯穿所有组学分析流程——从构建非冗余蛋白库、宏基因组基因集去冗余，到转录组组装后压缩序列，再到进化分析前去除冗余，几乎都绕不开它。

这篇内容适合谁？如果你刚刚接触生信，手上有一批fasta序列要去冗余，或者想搞明白-c、-n这些参数怎么选，可以照着操作；如果你已经在用CD-HIT，但对-n、-aS这些参数背后的逻辑有点含糊，这篇文章也能帮你把参数和原理打通。我会把原理、安装、命令、参数、案例和排错一次性讲透，全程用实际经验说话。

1. 先搞清楚CD-HIT在替你省哪一步：核心原理与应用场景

1.1 为什么要靠聚类来去冗余

先想想如果没有CD-HIT，你会怎么去冗余。最直觉的做法就是两两比对：每一条序列和另外所有序列做一次全局比对，相似度超过阈值就算冗余。但序列数量一旦上万，两两比对的次数就是n(n-1)/2，三万条序列就是将近5亿次比对。这个计算量用单机跑，基本是灾难，而且大多数比对结果都是“不相似”，浪费得离谱。

另一个麻烦是数据库膨胀。比如从宏基因组里预测基因，一个样本就可能出来几十万条基因序列，人群队列加起来上亿条。如果不去冗余，下游做功能注释、丰度统计、差异分析时，计算量会成倍增长，还会因为近缘序列过多导致统计偏差——某个基因家族的感觉丰度被高估，但实际只是同一个基因的不同拷贝。

CD-HIT解决的正是这个问题。它把“两两比较所有序列”变成“每条序列只和已有的簇代表比较”，并且用短词预过滤把绝大多数不相关序列在比对前就排除掉，因此速度能比暴力比对快几个数量级。典型应用场景包括：

• 构建非冗余蛋白数据库，比如把UniProt里几千万条序列压缩成某个阈值下的代表集
• 宏基因组/宏转录组基因集去冗余，生成非冗余基因目录
• 转录组组装后的转录本压缩，去除等位基因和组装碎片带来的冗余
• 蛋白家族分析前去除冗余序列，避免某些高丰度家族在统计里被过度放大
• 病毒、细菌基因组近缘序列清理，去掉同一株系反复测序造成的重复

可以说，CD-HIT是那种“单看功能很简单，但几乎所有流程里都能用上”的工具。

1.2 贪心增量聚类和短词过滤：CD-HIT快在哪儿

CD-HIT之所以快，靠的不是更聪明的比对算法，而是两个策略：贪心增量聚类和短词过滤。

先说说贪心增量聚类。它的流程是这样的：

1. 把输入的所有序列按长度从长到短排序
2. 取当前最长的序列作为第一个簇的代表序列
3. 读入下一条序列，只和已有的簇代表序列比较
4. 如果这条序列和某个代表的相似度达到阈值，就归入该簇，不再比较其他代表
5. 如果所有代表都不满足阈值，就把它当作新的簇代表
6. 重复直到所有序列处理完

整个过程的关键在于“每条新序列只和代表序列比，不和簇里所有成员比”。这意味着比较次数接近序列数乘以簇数，而不是序列数两两相乘。对高度冗余的数据集，簇数远小于序列数，节省的计算量非常可观。

短词过滤则是另一道加速关卡。CD-HIT在精确比对之前，会先把序列切成固定长度的短词（word），比如长度为5的氨基酸片段。如果两条序列之间连足够多的相同短词都没有，程序直接判定它们不可能达到相似度阈值，连比对都不用做。只有通过了短词过滤，才会进入条带比对，计算精确的相似度。

可以这么理解：短词过滤像是先看两个人有没有共同好友，没共同好友就直接跳过；有共同好友，才坐下来仔细对一遍。这样能过滤掉大量明显不相关的序列。这条原理也解释了为什么-n参数（词长）会影响聚类灵敏度——词长设得越大，预过滤条件越严，很多低相似度但真实同源的序列可能连比对机会都拿不到。

2. 装好CD-HIT，先把最常用的几个命令跑通

2.1 安装方式对比：conda最省心

CD-HIT的安装没什么坑，三种方式都试过，最推荐conda。

用conda装最省心，一条命令搞定：

conda install -c bioconda cd-hit

如果你有apt，也能直接装：

sudo apt install cd-hit

不过apt源里的版本可能比较旧，不建议在生产环境用。想自己编译源码也不难，去GitHub上拉最新的release，解压后直接make：

wget https://github.com/weizhongli/cdhit/releases/download/V4.8.1/cd-hit-v4.8.1-2019-0228.tar.gz tar -xzf cd-hit-v4.8.1-2019-0228.tar.gz cd cd-hit-v4.8.1-2019-0228 make

编译完会生成cd-hit、cd-hit-est、cd-hit-2d、cd-hit-est-2d等可执行文件，建议把目录加进PATH，或者把可执行文件软链到/usr/local/bin。我个人用下来，源码编译最稳，因为不依赖conda环境，也不会有动态库缺这缺那的问题。装完之后跑一下cd-hit -h验证版本和参数列表是否正常。

2.2 五个高频命令：从蛋白到核酸一次说清

CD-HIT不是单个命令，而是一族工具。日常最常用的有五个：

蛋白序列去冗余，用cd-hit：

cd-hit -i proteins.fa -o proteins_nr90.fa -c 0.9 -n 5

核酸序列去冗余，用cd-hit-est：

cd-hit-est -i contigs.fa -o contigs_nr95.fa -c 0.95 -n 10

注意，核酸序列一定不能拿cd-hit去跑，蛋白序列也不要用cd-hit-est。虽然命令长得很像，但内部处理的word长度、替换矩阵和比对逻辑都不一样，混用会让结果完全不可信。

如果你有两个库，想去掉第一个库里和第二个库相似的序列，用cd-hit-2d：

cd-hit-2d -i lib1.fa -i2 lib2.fa -o lib1_nr.fa -c 0.9 -n 5

这个命令在污染序列清理、去背景序列时很好用。处理逻辑是：保留第一个库中与第二个库相似度低于阈值的序列，相当于把lib1里被lib2“覆盖”的冗余序列剔除。核酸版本对应的是cd-hit-est-2d。

还有一个容易被忽略的命令是cd-hit-div，它可以按不同阈值分层聚类，适合构建多级别非冗余数据库，但日常用到的少，这里不展开。

2.3 第一次运行要看懂哪些日志

第一次跑CD-HIT，别急着看结果文件，先看屏幕输出。以蛋白去冗余为例，运行时会打印类似这样的日志：

Program: CD-HIT, V4.8.1 Command: cd-hit -i proteins.fa -o proteins_nr90.fa -c 0.9 -n 5 -T 8 -M 8000 total seq: 30000 longest and shortest sequences: 856 aa, 32 aa clustering at 90% identity, word length = 5 ... 17125 seqs. have been clustered. output written to proteins_nr90.fa

这里最值得看的是clustering at 90% identity, word length = 5这一行，它会明确告诉你当前用的阈值和词长。如果显示的词长和你预期的不一致，说明命令写法有误或者程序自动调整了参数，要回头检查。

日志末尾的xxx seqs. have been clustered表示有多少条序列被归并到了已有簇里。用总序列数减去这个数，就是最终的代表序列数。比如3万条序列，17125条被聚类，剩下12875条就是代表序列，也就是非冗余库的大小。

3. 关键参数逐一拆解：不同阈值到底怎么选

3.1 相似度阈值-c和词长-n是绑定关系

-c是最核心的参数，表示序列相似度阈值，取值范围0到1。它在蛋白去冗余里最常见的取值是0.9，对应UniProt90那种标准；宏基因组基因集去冗余常用0.95或者0.9；如果只是粗略清理，可以放到0.7到0.8。

但-c不能单独拍脑袋定，它和词长-n必须配套。-n作用于短词过滤，决定了预过滤的灵敏度。如果词长和阈值不匹配，会出现两种情况：词长太大，低相似度的同源序列被提前过滤掉，聚类结果偏少；词长太小，预过滤形同虚设，计算量暴涨。

根据官方文档和实际使用经验，我常用的搭配是下面这个表：

注意这只是经验区间。如果序列普遍很短（比如长度只有词长的两三倍），可以适当调小-n；如果序列很长且要求极高相似度，再考虑增大-n。默认情况下cd-hit的-n是5，cd-hit-est的-n大约是10，适合大多数场景。

还有一个容易踩的坑：有些版本里-n设置过大时，程序会直接报错说word size太长。这是因为有些短序列根本切不出足够多的短词，程序没法处理。遇到这种情况，先过滤掉过短序列，或者降低-n。

3.2 覆盖度限制-aS与-aL：全局相似还是局部相似

很多人只盯着-c，忽略了覆盖度参数-aS和-aL，结果聚类结果和自己预期差很远。这两个参数控制的不是相似度，而是比对的覆盖范围。

-aS是短序列中参与比对的区域占比下限，-aL是长序列中参与比对的区域占比下限。默认都是0，表示不限制。

举个例子。序列A长100个氨基酸，序列B长80个氨基酸，两条序列比对后一致区域长度是60个氨基酸。此时A的覆盖度是60%，B的覆盖度是75%。如果只设-c 0.9，只要这60个一致区域里的相似度够高，两条序列就可能被聚到一起，哪怕B的另外25个氨基酸完全对不上。但如果你希望每条参与聚类的序列都要有足够的全长覆盖，就要加-aS：

cd-hit -i proteins.fa -o proteins_nr90.fa -c 0.9 -n 5 -aS 0.8

这样设定之后，短序列B必须有至少80%的序列参与比对，也就是至少64个氨基酸，否则即使相似度超过0.9也不会被聚到一起。

在宏基因组基因去冗余时，我强烈建议加上-aS 0.8或-aS 0.9。因为宏基因组里很容易出现两条基因只在某个功能域上高度相似，但整体序列差异很大。如果不加覆盖度限制，它们会被错误聚类成同一条基因，导致后续功能注释和丰度统计失真。

3.3 模式选择-g、线程-T、内存-M以及其他收尾参数

-g参数控制一致性计算模式。默认-g 1使用全局序列一致性，速度更快，适合去冗余场景。如果设成-g 0，则使用局部一致性，会更容易把共享某段高相似区域的序列聚到一起，但速度明显变慢，而且结果里可能会混入只共享部分序列的成员。除非你刻意想找“共享结构域”级别的相似，否则我建议保持默认。

-M参数控制内存上限，单位是MB。默认值在某些版本里是800MB，对大数据集完全不够。我用3万条蛋白序列时一般给到-M 8000（8GB）。如果不想让内存成为瓶颈，可以直接设-M 0，让程序不限制内存。但要注意，这可能会把机器内存吃满，特别是输入序列很多的时候，建议提前用free -h看一眼可用内存。

-T参数控制线程数，例如-T 8表示用8个线程。CD-HIT的多线程加速主要集中在条带比对阶段，不是全程都能线性加速，所以线程数不是越高越好。我一般设成物理核心数的一半到全部，比如8核机器用-T 8，16核用-T 8到-T 12。

还有几个收尾参数比较实用。-d控制输出FASTA文件中序列标识的显示长度，默认是20个字符，如果想保留完整序列ID，设-d 0。-sc可以在.clstr聚类文件里输出具体的相似度百分比，方便你事后判断聚类质量。-sf设为1时，输出文件会包含每个簇的所有序列，而不只是代表序列，这在需要拿到完整簇成员时非常有用。

4. 实操案例：把三万条蛋白序列按90%一致性去冗余

4.1 输入数据准备与预期结果

先说一个我实际跑过的场景。有一个蛋白序列集合，大约三万条，来源是某个宏基因组样本的基因预测结果，里面有很多近乎一样的片段，我需要按90%的一致性去冗余。

第一步是检查数据格式。CD-HIT对输入FASTA的要求比较基础但严格：序列ID建议不要有空格，序列里不要有非法字符，短序列最好提前过滤掉。如果ID带空格，后续解析.clstr文件时会非常痛苦。我习惯用seqkit做一步清洗：

seqkit seq -m 30 proteins.fa > proteins_filtered.fa

这段命令会丢弃短于30个氨基酸的序列，同时把多行序列按规范格式整理。为什么要过滤短序列？因为过短的序列在聚类时意义不大，还可能干扰word length判断，尤其是当-n设置较大时，短序列根本切不出足够多的短词。

数据准备好了，不要急着跑。可以先想清楚预期：三万条来自宏基因组的蛋白序列，在90%阈值下，通常能压缩掉20%到40%。如果冗余度特别高，甚至能压缩掉一半以上。如果压缩率过低，就要怀疑数据是不是本身多样性就高，或者参数设置有问题。

4.2 执行命令与运行日志解读

按前面分析的参数逻辑，90%阈值对应蛋白序列词长5。我用的命令是：

cd-hit -i proteins_filtered.fa -o proteins_nr90.fa -c 0.9 -n 5 -T 8 -M 8000 -d 0 -sc 1

运行时间取决于机器。我当时用的是一台8核服务器，跑了大概不到十分钟就结束。运行日志会像这样：

total seq: 30000 longest and shortest sequences: 856 aa, 32 aa clustering at 90% identity, word length = 5 ... 17125 seqs. have been clustered.

这里longest and shortest sequences可以用来快速判断输入数据有没有异常，比如最短序列只有几个氨基酸，那大概率是数据预处理没做好。17125 seqs. have been clustered的意思是：有17125条序列被归并到了现有簇中，没有成为新的代表序列。也就是说，最终代表序列数量大约是30000减去17125，也就是12875条。

这里还有个细节值得注意：虽然我设了-c 0.9，但实际聚类时不同簇里的序列和代表序列之间的最小相似度不完全等于90%。因为CD-HIT做的是贪心聚类，某条序列可能匹配到第一个达到阈值的代表就停止比较，而不会在所有代表里找最优者。这是所有贪心聚类的共性，也是它快的代价。如果项目对聚类质量要求极其苛刻，比如要保证每个簇内所有序列两两相似度都超过阈值，CD-HIT并不是最合适的选择，你可能要去看更精细的聚类算法。

4.3 解读.clstr文件并统计聚类情况

命令跑完后会生成两个文件：proteins_nr90.fa是代表序列FASTA文件，proteins_nr90.fa.clstr是聚类结果文件。

.clstr文件的结构长这样：

>Cluster 0 0 452aa, >seq_001... * 1 320aa, >seq_002... at 91% 2 180aa, >seq_003... at 95% >Cluster 1 0 378aa, >seq_004... * 1 201aa, >seq_005... at 90%

每个>Cluster开头的是一个簇，后面每一行代表簇里的一个成员。这一行的格式从左到右依次是：簇内成员序号、序列长度、序列标识、相似度信息。带*的是代表序列，也就是每条簇里最长的那条。at 91%表示这条序列和代表序列的一致性达到91%。

如果你想看聚类结果的整体情况，用grep数一下簇的数量：

grep -c "^>Cluster" proteins_nr90.fa.clstr

这个数字应该和代表序列FASTA文件里的序列数一致。如果不一致，说明输出文件有问题，可能和中途中断有关。

我还常用一条awk统计每个簇的序列数分布：

awk '/^>Cluster/{if(c>0) print c; c=0; next} {c++} END{print c}' proteins_nr90.fa.clstr | sort -n | uniq -c

这样可以快速看出来绝大多数簇是只有一条序列，还是有大量几十条序列的大簇。如果出现某个簇包含上千条序列，要留意这个簇是不是过于膨胀，可能是某个重复序列家族或者隐藏的污染。

4.4 核酸序列和宏基因组的场景延伸

蛋白去冗余做完，同样的逻辑可以平移到核酸序列上。唯一要记住的是换命令：用cd-hit-est。

核酸序列按95%阈值去冗余的命令：

cd-hit-est -i genes.fa -o genes_nr95.fa -c 0.95 -n 10 -T 16 -M 0 -d 0 -aS 0.8

这条命令在宏基因组基因集构建里非常常见。为什么阈值选0.95而不是0.9？因为核酸序列的冗余更多来自近缘菌株的等位基因或者同一基因的微小变异，95%能在保留功能多样性的同时有效压缩数据量。加上-aS 0.8，是避免两条只共享一个外显子或功能域的基因被误聚到一起。

宏基因组里还有个常见操作是构建非冗余基因目录，这时候往往是先把所有样本的基因预测结果合并，然后用cd-hit-est按95%聚类。数据集如果特别大，比如超过几百万条序列，CD-HIT可能会跑得比较吃力。这种情况下我一般先用MMseqs2快速粗聚一轮，再用CD-HIT对结果做精细去冗余，两者配合在速度和质量上都能兼顾。

5. 常见问题与排错实录

5.1 运行报错类

用CD-HIT这么多年，最常遇到的问题基本集中在内存、序列格式和词长这几个点上。我从实际踩坑经验里整理了下面这个表格。

我自己印象最深的一次内存问题，是在处理一条包含上百万条蛋白序列的蛋白库时，忘了设置-M，结果程序跑到一半被系统直接杀掉，前面的计算全部白费。从那以后我养成的习惯是：跑大规模数据前先跑一条小命令，比如samtools view -H或者seqkit stats，确认机器可用内存和输入文件大小，再决定-M给多少。

5.2 结果不对类

比报错更麻烦的是结果错而不自知。下面这几个情况是我帮别人排查时见过最多的问题。

第一条是核酸序列用错程序。cd-hit和cd-hit-est长得太像，很容易搞混。如果你把DNA序列用cd-hit跑，程序不一定报错，但聚类结果会有问题，因为word length和比对模型都不匹配。这个坑尤其隐蔽，因为我见过有人跑完还拿着结果做了下游分析，直到怎么都解释不通才发现。

第二条是覆盖度参数没设置导致误聚类。前面提过，不加-aS和-aL时，CD-HIT允许局部相似度够高就把整条序列聚到一个簇。这在某些场景下是bug，但在另一些场景下是特性。关键是你要明确自己到底想按“全长相似”还是“局部相似”来聚类。想做严格去冗余就加-aS 0.8以上，别偷懒。

第三条是-n设置错误导致聚类结果异常偏少或偏多。有个朋友做蛋白序列，-c设了0.7，但-n没改仍是5，结果聚类结果明显偏少。后来把-n改成4，结果合理很多。这说明-c和-n的绑定关系不是纸面上的理论，而是会影响实际结果的操作细节。

第四条是贪心聚类导致的“次优”结果。CD-HIT的贪心策略决定了，如果两条高度相似的序列在比较时，先碰到了另一条也满足阈值的序列，它们会被分到不同的簇。这在很多下游分析中无所谓，但如果你要做精细的直系同源基因分析，最好在CD-HIT聚类之后再结合系统发育做二次确认，不要完全依赖聚类结果。

5.3 超大数据量的处理思路与替代工具

CD-HIT适合的规模是几十万条序列这个量级。到了百万甚至千万条序列，即使能跑完，时间和内存都会非常吃紧。这时候不要硬撑，可以考虑下面这些思路。

思路一是用更快的工具做粗聚类。MMseqs2是目前最成熟的大规模序列聚类工具之一，它的linclust模式专门为超大数据库设计，速度能比CD-HIT快一个甚至两个数量级。我的工作流是：先用MMseqs2按95%粗聚类，把序列压到一个合理规模，再用CD-HIT按需要的阈值做精细去冗余。这样既拿到了CD-HIT的精确控制，又避免了长时间等待。

思路二是切成小批次分别运行，最后合并簇。这个操作要小心，因为跨批次的冗余关系在合并时可能丢掉。我一般是先按样本或者染色体区域把序列分开，每个分区单独跑，最后再用cd-hit-2d做一次跨分区去冗余，把不同批次之间的冗余也消掉。

思路三是考虑VSEARCH或者UCLUST这类工具。它们在OTU聚类、16S分析里用得很广，尤其是核酸序列上的表现值得关注。不过如果目标是蛋白序列的高精度去冗余，我个人还是首选CD-HIT，毕竟它在蛋白领域的参数调教和结果可解释性上更成熟。

回到我自己的经验，这些年用下来，最顺手的CD-HIT工作流其实是这样的：拿到输入数据后，先用seqkit stats看一下序列长度分布，过滤掉过短序列；然后抽一小批序列，比如1万条，在不同参数组合下各跑一遍，看看代表序列数量的变化趋势；确定好阈值和词长之后，再整库跑。整个过程其实没多少技术难度，真正值钱的判断力都在参数选择的理解上。比如什么时候该加-aS 0.9，什么时候可以放心用-c 0.7，这些判断来自你对数据本身的认识，而不是单纯背参数。最后一个我常提醒自己的细节：无论多信任CD-HIT，跑完后都去看一眼.clstr文件里的簇大小分布，异常的大簇往往意味着数据里有污染物或者重复区域，这时候先别急着往下游走，回头检查数据本身才是最省时间的做法。