CD-HIT 序列聚类详解三步装好并跑通去冗余全流程【免费下载链接】cdhitAutomatically exported from code.google.com/p/cdhit项目地址: https://gitcode.com/gh_mirrors/cd/cdhit数据库从一万条序列涨到一千万条后跑一遍 BLAST 要多久答案很大程度取决于库里还剩多少冗余——真实生物数据里高度相似的序列占比远超直觉冗余直接拖慢后续所有比对分析。CD-HIT 干的就是这件事按序列相似度把冗余序列聚类分组每个簇只留一条代表序列把数据库规模压下去。它是一组 C 写的序列去冗余与聚类程序编译后就是cd-hit蛋白、cd-hit-est核酸等几个可执行文件命令行几行即可上手。本文从安装、首次运行、原理讲起再带你完成两个真实任务——蛋白库去冗余和 16S rRNA OTU 聚类最后给出多核与多机器的扩容策略。三分钟装好并跑通编译只依赖 g默认启用 OpenMP 多线程和 zlib用于读.gz输入没有第三方库要装git clone https://gitcode.com/gh_mirrors/cd/cdhit cd cdhit make # 生成 cd-hit、cd-hit-est 等 6 个程序 cd cd-hit-auxtools make # 生成 cd-hit-dup 等辅助工具环境要求如下绝大多数 Linux 服务器开箱即用项目最低推荐操作系统Linux / macOS任意 Linux内存4 GB16 GB 以上CPU双核8 核以上首次运行用一段手工构造的蛋白序列验证。seq_a与seq_b只差 1 个残基seq_c与它们完全不同cat demo_protein.fa EOF seq_a MSQVKESQDILQELQELSAELHESQELQELQELQELQELQEL seq_b MSQVKDSQDILQELQELSAELHESQELQELQELQELQELQEL seq_c MFGLLPVTGQLLQVLLALFIVGILFV EOF ./cd-hit -i demo_protein.fa -o demo90 \ -c 0.90 # 相似度阈值 90% -n 5 # 词长 k5适用于 0.7~1.0 -d 0 # 结果中只保留 fasta 头第一列 ID -T 2 # 用 2 个线程 -M 800 # 内存上限 800 MB运行结束会打印Total CPU time并生成两个文件demo90代表序列 fasta和demo90.clstr簇清单Cluster 0 0 45aa, seq_a... * 1 45aa, seq_b... at 98% Cluster 1 0 17aa, seq_c... **标记簇的代表序列at 98%是该序列与代表序列的一致率。含义很直观seq_b被归入seq_a的簇输入 3 条序列输出库只剩 2 条。用./cd-hit -h可以随时查看全部参数。它是怎么聚类的——机制讲人话CD-HIT 的速度不是靠更精确的比对而是靠少做比对。三个关键设计长序列优先的贪心策略。输入先按长度从长到短排序第一条序列直接成为代表序列之后每条序列只与已有的代表序列比较相似就并入、不相似就自己开簇。这样比对对象从全部 N 条降到代表序列那一小撮量级差别是平方关系。默认快速模式把序列归入第一个达标的簇加-g 1则逐一比较、归入最相似的簇慢但准。索引表 短词过滤。CD-HIT 给每个 k-mer 建一张索引表蛋白 k 取 2~5核酸 k 取 8~12全部 k-mer 组合恰好能装进内存。两条序列若相似度达到阈值它们共有的相同 k-mer 数必有下限于是先数公共词数量不够直接放弃根本不启动比对数量够时相同词的位置还能标出比对大概落在哪一段只需沿一条窄带做带状动态规划而不是整条序列全局比对。覆盖度约束。相似不等于能归簇两条序列对齐后还要看覆盖情况短序列对齐部分至少占代表序列的-s比例对齐至少覆盖长、短序列的-aL、-aS比例。这些参数就是下面这张图里各段的含义。CD-HIT 参数速查必需参数参数作用建议取值-i输入文件fasta 或.gz4.8.1必填-o输出前缀生成-o和-o.clstr两个文件必填常用参数参数作用建议取值-c相似度阈值蛋白 0.916S rRNA 0.97-n词长 k蛋白 0.7~1.0 用 50.6~0.7 用 40.5~0.6 用 30.4~0.5 用 2核酸 0.95~1.0 用 100.90~0.95 用 8~90.88~0.90 用 7-T线程数0 自动用满全部核-M内存上限MB0 自动调优默认 800-d.clstr中描述信息长度0只留 ID偶尔用的参数参数作用建议取值-s短序列长度至少占代表序列的比例0.9-l短于该长度的序列直接丢弃10默认大数据可用 20-g1 归入最相似簇精确慢模式0-sc/-sf按簇大小排序输出0-t冗余容忍度允许少量错归换速度2默认同一家族的子命令一个一行cd-hit-est核酸聚类词长默认 10支持-P 1 -i R1 -j R2处理双端数据cd-hit-2d/cd-hit-est-2d两个库之间交叉比对-i db1 -i2 db2找出 db2 中相似于 db1 的序列cd-hit-454454 平台读长去重复cd-hit-div把超大库切分成子库服务分层聚类流程psi-cd-hit蛋白相似度低于 40% 时的聚类需另装 BLAST 作为后端两个完整任务演练任务一给蛋白数据库做去冗余目标把一份蛋白库按 90% 相似聚类成非冗余库供后续功能注释或 BLAST 使用。输入protein.fasta标准 fasta 或.gz。命令./cd-hit -i protein.fasta -o nr90 \ -c 0.90 -n 5 \ -d 0 -T 8 -M 16000 \ -s 0.9 # 防止过短序列混入长序列的簇结果解读nr90是可直接替换原库使用的代表序列文件。簇规模用仓库自带脚本统计perl clstr_size_stat.pl nr90.clstr # 各大小簇的数量分布 perl clstr_rep.pl nr90.clstr rep.tsv # 簇编号、代表 ID、簇大小clstr_rep.pl的输出是制表符分隔的三列文本方便用awk或表格软件继续处理。若发现大量单条序列组成的簇先检查-c是否设得过高再看-s是否过严。任务二16S rRNA 测序数据聚 OTU目标把 MiSeq 双端 16S 测序数据聚成 OTU并与全长参考库联合聚类一步完成分群加注释。输入每个样本一个目录内含R1.fq、R2.fq一份样本清单SAMPLE_file每行样本名 R1.fq R2.fq一份 16S 全长参考库。前置安装 Trimmomatic质控用并把usecases/Miseq-16S/NG-Omics-Miseq-16S.pl顶部的CD_HIT_dir和 Trimmomatic jar 路径改成你自己的。命令先把参考库拼成与样本同形的 PE 格式-p 150 -q 100分别指定 R1、R2 参与聚类的高质量区长度perl usecases/Miseq-16S/16S-ref-db-PE-splice.pl \ -i s1/R1.fq -j s1/R2.fq \ -d Greengene-13-5-99.fasta \ -o gg_PE99 -p 150 -q 100 -c 0.99再生成所有样本的质控与 OTU 脚本并执行0.97是 OTU 聚类阈值0.0001是丰度过滤阈值perl usecases/Miseq-16S/NG-Omics-WF.pl \ -i usecases/Miseq-16S/NG-Omics-Miseq-16S.pl \ -s SAMPLE_file -j otu \ -T otu:150:100:0.97:0.0001:/完整路径/gg_PE99-R1:/完整路径/gg_PE99-R2:75 \ -J write-sh生成的 shell 脚本分两步跑先全部qc.*.sh再全部otu.*.sh。每个样本的otu/OTU.clstr列出每个 OTU 的成员读段chimeric-small-clusters-list.txt记录被剔除的嵌合序列和低丰度簇。多样本项目最后合并perl usecases/Miseq-16S/pool_samples.pl -s SAMPLE_file -o pooledpooled/OTU.txt是每个样本中各 OTU 的计数表加注释可直接作为下游差异分析的输入。再往上走——性能与规模策略先把单机资源用满-T 0让程序自动占满全部核心配合-M 0让内存分配自动调优这两个组合是大数据集的第一选择./cd-hit -i big_db.fasta -o out -c 0.9 -n 5 -T 0 -M 0预过滤短序列-l 20直接丢弃短于 20 aa/nt 的序列它们多数没有生物学意义过滤后内存占用和运行时间都会明显下降。.gz输入4.8.1 支持则省去解压步骤和磁盘占用。超大库分而治之单库大到内存吃紧时用cd-hit-div把大库切成若干子库分别聚类再用cd-hit-2d做子库间交叉比对、合并结果即先粗后细的分层流程跨机器并行多台机器时改用perl cd-hit-para.pl -i huge.fa -o out -c 0.9 -n 5脚本自动把输入切块、分布到各节点聚类后再归并适合亿级读长场景。踩坑速查现象低相似度下聚类结果异常很多该聚的没聚上原因-n超出了-c的适用范围比如-c 0.5仍用-n 5短词过滤的统计前提失效。 解法对照速查表换词长0.5~0.6 用-n 30.4~0.5 用-n 2低于 40% 的蛋白相似聚类超出本算法边界改用psi-cd-hit。现象进程中途被系统杀掉输出不完整原因-M默认只有 800 MB大数据集装不下触发系统 OOM。 解法按物理内存调大-M或设 0 自动调优-T 0摊到更多核心仍不够就-l提高丢弃阈值或按上一节的分层流程切库分步跑。现象代表序列看着不像簇内最像的那条原因默认是快速模式序列归入第一个达标的簇-t 2还允许少量错归换速度。 解法加-g 1切换精确模式归入最相似簇速度慢对长度一致性有要求时再补-s 0.9。CD-HIT 的适用边界是明确的蛋白相似度阈值 40%~100% 区间内、序列能装进内存的场景都能跑单机聚类数亿条蛋白序列在一天内完成是其常规能力低于 40% 相似或基因组尺度的超长序列聚类请改用psi-cd-hit。它的定位始终是把冗余换成速度——先聚类再做一切下游分析。引用说明使用 CD-HIT 发表成果时请引用Li W, Godzik A. CD-HIT: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658-1659.【免费下载链接】cdhitAutomatically exported from code.google.com/p/cdhit项目地址: https://gitcode.com/gh_mirrors/cd/cdhit创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
