
如何用CD-HIT给百万级序列快速去冗余从安装、聚类到避坑的完整指南【免费下载链接】cdhitAutomatically exported from code.google.com/p/cdhit项目地址: https://gitcode.com/gh_mirrors/cd/cdhit手上几百万条蛋白质序列想直接做功能注释BLAST全量跑下去服务器先扛不住——这正是CD-HIT序列聚类要解决的问题按你指定的相似度阈值先把序列库压成一份非冗余代表序列库再拿小库去跑下游分析。官方实测32核机器聚合一亿条蛋白质序列不到一天。如果你是做宏基因组、转录组或微生物组分析的新手这篇带你一次装好、跑通、调对参数。 先认清它是谁序列查重器不是序列比对器CD-HIT更像给序列库做查重的工具就像相册里有100张几乎一样的照片它挑一张留作代表并告诉你剩下99张各自跟哪张重、重多少——它不评价照片好不好也不做进化分析。适合干的事蛋白质序列聚类去冗余cd-hit核酸/EST/读段聚类cd-hit-est两个库互相比对找出B库里相对A库的新序列cd-hit-2d测序读段的重复去除cd-hit-dup在 cd-hit-auxtools/ 目录编译明确不会干的事能力边界别踩蛋白质相似度低于40%的聚类——底层短词过滤器在这个区间失效官方替代方案是 psi-cd-hit/基因组级别的超长序列——cd-hit-est处理不了这类数据请换工具系统发育/同源关系判定——它只回答谁和谁像树分析需要别的软件 首次运行3条命令装好跑通最小聚类最短安装路径Linux已有ggit clone https://gitcode.com/gh_mirrors/cd/cdhit # 获取源码 cd cdhit make # 编译主程序默认多线程版 cd cd-hit-auxtools make cd .. # 编译辅助工具两个编译注意点4.8.1起支持直接读.gz压缩输入需要zlib编译报zlib错误就先装开发包Ubuntusudo apt install zlib1g-dev极老系统不支持OpenMP时改用make openmpno最小可跑通命令3条测试蛋白、90%相似度聚类cd-hit -i test.fasta -o test_out -c 0.9 -n 5跑完你会得到两个文件这就是全部输出形态test_out代表序列fasta每个簇留的那一张照片test_out.clstr文本簇表Cluster 0开始一个新簇行尾带*的是该簇代表后面的百分比是它与代表的相似度️ 参数不用背按输入输出/准确性/性能三个抽屉归类面对一长串参数不慌关键是知道什么情况下才需要动它。输入输出抽屉——换数据时动。-i/-o是仅有的两个必需参数输入输出文件各一个。生产任务记得加-d 0默认.clstr里只保留20个字符的描述-d 0改成保留header到第一个空格为止的完整标识否则下游数据库回查时对不上号。准确性抽屉——觉得聚类不干净时动。-c相似度阈值蛋白质常用0.9核酸常用0.95~0.97-nk-mer词长相当于先按多长的单词开头做粗筛的筛子规格必须与-c匹配蛋白质-n 5对应阈值0.7~1.0-n 2才能下探到0.4~0.5核酸-n 10,11对应0.95~1.0-n 8,9对应0.90~0.95-aL/-aS比对覆盖长/短序列的比例要求防短序列被长序列一口吞-g 1精确模式序列归入最相似的簇默认是归入第一个达标的簇快但不一定最优上图就是-aL/-aS的几何含义R是代表序列S是待比对的短序列参数管的是比对段必须覆盖各自多少比例。性能抽屉——跑得慢或内存爆时动。-T线程数0表示用满所有CPU-M内存上限MB默认800大数据库建议给到16000以上0为不限制小心OOM-l丢弃短于该长度的序列默认10用来提前砍掉短碎屑。 实战场景一百万级蛋白质库去冗余场景拿到SPROT全量蛋白要做90%非冗余库再跑功能注释。输入输出输入一个fasta或.gz输出代表序列fasta 簇表。cd-hit -i uniprot_sprot.fasta -o uniprot_nr90 -c 0.9 -n 5 -d 0 -T 8 -M 16000结果怎么读uniprot_nr90的行数就是代表序列数——500万压到200万说明冗余率60%打开uniprot_nr90.clstr抽查每段Cluster x是一个簇带*的行是代表。想快速看簇大小分布跑一句perl clstr_size_stat.pl uniprot_nr90.clstr即可。 实战场景二核酸序列聚类cd-hit-est场景EST或转录本数量大要按95%相似度去冗余。输入输出同上两个文件。核酸字母表只有A/C/G/T四种短词随机撞上的概率远高于蛋白质所以-n默认就是10而不是5别照抄蛋白命令。cd-hit-est -i transcripts.fasta -o transcripts_nr95 -c 0.95 -n 10 -d 0 -T 8 -M 16000结果怎么读与场景一相同的两个输出文件。含内含子的转录本会有长gaps稀释局部一致性这类数据建议追加-aL/-aS覆盖约束把只像一小段的配对挡在簇外。 实战场景三16S rRNA的OTU聚类微生物组场景MiSeq双端16S读段要按97%相似度聚出OTU并顺便完成注释。输入输出输入PE read文件输出OTU簇表、注释表和运行日志。仓库自带完整工作流在 usecases/Miseq-16S/ 目录核心脚本是cd-hit-otu-miseq-PE.pl。它的特色是不要求先把PE两端拼接成contig而是直接取R1/R2的5高质量区例如R1前200bp、R2前150bp参与聚类还能把全长16S参考库的目标区如V3-V4剪接出来与样品一起聚聚类和注释一次做完。perl usecases/Miseq-16S/cd-hit-otu-miseq-PE.pl -i sample_PE.fq -o otu_result -c 0.97结果怎么读输出目录下的OTU.log记录全流程-c 0.97是OTU默认阈值对应16S菌种水平划分的行业惯例。完整依赖需要Trimmomatic质控和路径配置方法见usecases/Miseq-16S/README。 避坑指南5个常见现象的原因与对策内存爆了怎么办现象程序中断报memory allocation failed。原因序列默认存内存-M上限不够或数据真的大。对策调低-M比如8000让它走磁盘溢出用-l过滤短序列减小体量仍不行就分阶段聚类。代表序列没少下来怎么办现象跑完发现簇数只比输入少一点点。原因-c定太高或-n与-c不匹配阈值低于0.7却用-n 5会让大量相似序列漏过粗筛。对策对照前面准确性抽屉的词长对照表检查-n把-c降到0.85~0.9或加-s长度差约束。不相关的短序列被并进大簇怎么办现象clstr里出现明显不该在一起的小片段。原因没开覆盖约束短序列在长序列上局部命中就被判冗余。对策加-aL 0.9 -aS 0.9要求比对覆盖90%以上。超大库太慢怎么办现象几千万条序列单机要跑很久。原因单进程瓶颈。对策先-T 0用满所有核真正的大库用下图的分层策略——cd-hit-div先按多样性把大库切成子库各子库内部用cd-hit聚类代表序列之间再用cd-hit-2d跨库比较最后合并。分布式集群可配合cd-hit-para.pl。高频坑速查表现象最可能原因第一优先对策memory allocation failed数据量超过-M调低-M触发磁盘溢出或-l预过滤代表序列数降不下来-c过高、-n与-c不匹配核对词长对照表降低-c短序列乱入大簇未开覆盖约束加-aL 0.9 -aS 0.9编译失败提示zlib相关缺zlib开发包安装zlib开发包后重新make蛋白相似度想低于40%算法边界换psi-cd-hit 收尾速查一页带走常用命令配套工具一句话仓库根目录有一批Perl脚本处理.clstr结果——clstr_rep.pl提取代表序列、clstr_size_stat.pl统计簇大小、clstr2tree.pl生成树文件、clstr_quality_eval.pl评估聚类质量用法都是perl 脚本名 xxx.clstr完整清单在 官方用户手册。一页速查表任务命令骨架安装编译git clone https://gitcode.com/gh_mirrors/cd/cdhit cd cdhit make蛋白90%去冗余cd-hit -i in.fa -o out -c 0.9 -n 5 -d 0 -T 8 -M 16000核酸95%去冗余cd-hit-est -i in.fa -o out -c 0.95 -n 10 -d 0 -T 8 -M 16000两库比较找新序列cd-hit-2d -i db1 -i2 db2 -o out -c 0.9 -n 516S OTU聚类perl usecases/Miseq-16S/cd-hit-otu-miseq-PE.pl -i PE.fq -o otu -c 0.97参数起步推荐值-c蛋白0.9 / 核酸0.95~0.97-n蛋白5对应阈值0.7~1.0/ 核酸10对应0.95~1.0-T0用满CPU-M16000起步按机器内存上调0不限制-d生产任务固定用0引用信息Li W, Godzik A. CD-HIT: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics, 2006.【免费下载链接】cdhitAutomatically exported from code.google.com/p/cdhit项目地址: https://gitcode.com/gh_mirrors/cd/cdhit创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考