CD-HIT快速上手：从一条命令到百万级序列聚类

CD-HIT快速上手从一条命令到百万级序列聚类【免费下载链接】cdhitAutomatically exported from code.google.com/p/cdhit项目地址: https://gitcode.com/gh_mirrors/cd/cdhitCD-HIT是一个蛋白质与核酸序列聚类/比较工具核心任务是对大规模序列库做去冗余给定相似度阈值它把相似序列归成簇每簇留一条代表序列。如果你手上有一批UniRef、转录本或16S读长需要快速算出非冗余版再交给下游比对、注释或网络构建这篇内部向的上手笔记就是给你看的。⚡ 30秒跑通一条命令看到聚类结果先在当前目录造一个3条序列的最小输入跑一遍最简聚类预期得到demo.clstr和demo两个输出文件。# 造一个3条序列的测试库seq001与seq002仅1个残基不同seq003差异大 cat demo.fa EOF seq001 MAAAAKLLCCGDDDEEEFFGGHHHIIJJJ seq002 MAAAAKLLCCGDDDEEEFFGGHHHIIKKK seq003 PPPQQQRRRSSSTTTUUUVVVWWWWWXXX EOF # 90%相似度阈值聚类-n 5为蛋白质默认词长 ./cd-hit -i demo.fa -o demo -c 0.9 -n 5打开demo.clstr你看到的是这样的结构每行以开头的行是一个新簇带*号的成员是该簇代表序列at XX%是成员与代表的相似度。Cluster 0 0 27aa, seq001 1 27aa, seq002 at 96% Cluster 1 0 27aa, seq003 *seq001和seq002归为一簇、长序列当代表seq003单独成簇——这就是CD-HIT的全部行为模型。 原理速览只有做参数决策时需要懂的部分不用看论文记住下面四条参数就不会乱填按长度降序聚类首条即代表。序列先排序最长的先当某簇代表后续序列只要与已有某簇代表达到阈值就并入该簇。所以哪条留下来当代表基本由长度决定想要最短的留下是反直觉的。k-mer预筛 精确比对两段走。先用词-n控制词长筛掉明显不相似的序列对再对剩下的做局部比对。-n不是越大越快词长必须匹配你的阈值区间官方对应关系是 蛋白质-n 5用于 0.7~1.0、-n 4用于 0.6~0.7核酸cd-hit-est默认-n 10。阈值设 0.9 却用-n 3速度和结果都不可靠。相似度的默认口径是全局。-G 1默认下相似度 相同残基数 / 较短序列全长。也就是说短序列挂到长序列代表上时只要它全身都对上就算数长序列的尾巴不管。这正是-aL、-aS这类覆盖度参数存在的原因。图1-aL/-AL/-aS/-AS 对齐覆盖度控制示意图来自用户手册快速模式 vs 精确模式。默认-g 0是碰到第一个达标的簇就入簇快但簇成员可能略杂-g 1会选相似度最高的簇慢但更干净。⚙️ 环境与安装make一遍加一条验证命令依赖只有C编译器和zlib4.8.1起支持.gz输入需要zlib大多数Linux自带。老系统编译不了多线程可退回make openmpno装不上zlib用make zlibno。# 获取源码当前版本4.8.1 git clone https://gitcode.com/gh_mirrors/cd/cdhit cd cdhit # 编译主程序默认带多线程 make # 编译辅助工具cd-hit-dup、cd-hit-lap等16S脚本依赖它 cd cd-hit-auxtools make验证命令不带参数运行./cd-hit会打印全部选项说明能打印出来即编译成功。配置项最低要求推荐配置操作系统Linux/Unix有gUbuntu 20.04内存4GB16GB聚类时序列默认驻留内存CPU双核8核以上-T多线程才划算磁盘2GB源码输出50GB大库输出会再产生一套代表序列 参数速查表必用、常用、进阶三档档位参数作用推荐值范围必用-i/-o输入fasta / 输出前缀必填必用-c序列相似度阈值蛋白0.9~0.95核酸0.95~0.98常用-n词长k-mer蛋白5核酸10随阈值下调常用-T线程数物理核数0表示用满所有CPU常用-M内存上限(MB)用于自动调词表4000~320000为自动常用-l短于此长度直接丢弃蛋白10/核酸50常用-s短序列长度须达代表的比例0.7~0.9默认0不限制进阶-aL/-aS对齐须覆盖长/短序列的比例0.9表示几乎全长对齐进阶-AL/-AS覆盖不足的绝对长度豁免配合-aL用如60进阶-t允许的冗余度容忍少量重复换速度默认2越大越快进阶-g1选最相似簇慢而准对结果敏感时设1进阶-B1序列存磁盘而非内存超大库必开进阶-d.clstr里描述截断长度0完整保留id去冗余分析常用两个高频坑-G 0局部相似度单独使用没有意义官方明确提示要配覆盖度参数才用-n与-c不匹配是结果不对的第一嫌疑。 实战工作流三个真实场景场景一蛋白质库去冗余问题UniProt SP有60万条直接拿去建BLAST库太慢需要90%非冗余版本。# 90%全局相似度聚类-s 0.8限制成员长度差不超过20%-d 0保留完整序列号 ./cd-hit -i sprot.fasta -o sprot90 -c 0.9 -n 5 \ -T 8 -M 16000 -s 0.8 -d 0得到sprot90代表序列fasta和sprot90.clstr簇归属。解读方法grep -c ^Cluster sprot90.clstr就是簇数单个簇想看成员用grep -A20 ^Cluster 0$ sprot90.clstr。后续若按簇做统计clstr_size_stat.pl可直接处理clstr文件需要代表序列回查全库时用clstr_rep.pl。场景二转录本聚类找异构体问题RNA-seq组装出的转录本里同一基因有多个剪接变体要按95%相似度归簇避免把不同异构体误并成一簇。# cd-hit-est是核酸版本-n 10为核酸默认词长 # -aS 0.9要求对齐覆盖短序列90%防止尾巴短一点就并入 ./cd-hit-est -i transcripts.fasta -o tr95 -c 0.95 -n 10 \ -T 4 -M 8000 -aS 0.9 -G 1注意这里-G 1是默认行为相似度按短序列全长算。转录本场景如果你更关心可变区是否不同可以改用-G 0并配套加-aL/-aS约束覆盖度——这是少数值得动-G的场合。场景三16S rRNA的OTU聚类MiSeq双端问题MiSeq下机的是双端原始reads需要先拼接、去重、分阶段聚类到OTU而不是直接拿cd-hit-est硬跑。项目自带一条完整流水线usecases/Miseq-16S/cd-hit-otu-miseq-PE.pl内部就是cd-hit-dup100%去重嵌合体过滤→ 99.25% → 97%三阶段递进聚类再比对参考库。# -i/-i2为双端reads-c 0.97为OTU相似度阈值输出到otu目录 perl usecases/Miseq-16S/cd-hit-otu-miseq-PE.pl \ -i R1.fasta -i2 R2.fasta -o otu -c 0.97前置条件cd-hit-auxtools已 make脚本会检查cd-hit-dup是否存在缺了直接报错退出。产物在输出目录内OTU.log记录各阶段序列数。图2CD-HIT-OTU-MiSeq 16S rRNA双端读数聚类流程来自usecases/Miseq-16S 性能调优先判断瓶颈再动参数内存还是磁盘默认序列驻留内存-B 0。报错或机器开始swap时第一反应是加-B 1把序列落盘而不是升内存。-M的真正作用是给程序预算去自动分配词表大小——把它设为物理内存的70%左右如16G机器给-M 12000比盲目调大更稳。线程何时有效聚类是串行主循环多线程加速主要吃在序列量大十万级以上且CPU核多时。经验值4条线程对4核机器的加速接近线性超过物理核数无收益还占内存。-T 0会用满所有逻辑核跑批时留2核给系统。单库跑不动就拆。官方方案是cd-hit-para.pl把输入按长度降序切成N块第一块用cd-hit自聚其余块逐块用cd-hit-2d打到已聚结果上最后合并。单机多块并行用cd-hit-2d-para.pl集群则传--B hosts主机列表。图3cd-hit-para.pl 分块并行聚类流程来自用户手册两轮聚类换精度。第一轮低阈值如0.85粗聚对代表序列再跑第二轮高阈值0.95能显著减少第二轮数据量适合大类归并类内细分的数据库构建需求。 排错三板斧症状运行中报内存分配失败或被OOM杀掉。原因序列驻留内存词表随输入规模膨胀。 解法按顺序试-M调小预算 →-B 1落盘 → 输入先-l过滤短序列 → 彻底跑不动就按上节拆库。症状簇数量与预期严重不符过多或过少。原因90%的情况是-c与-n不匹配或阈值口径用错全局 vs 局部。 解法回到-n与阈值对应表核对确认-G取值过聚时加-aL 0.9收紧覆盖度欠聚时检查-s是否把本该入簇的短序列挡在门外。症状编译报zlib相关错误。原因4.8.1起读.gz输入依赖zlib头文件。 解法Ubuntusudo apt install zlib1g-devCentOSsudo yum install zlib-devel或干脆make zlibno放弃gz支持。症状用16S流水线时脚本报no .../cd-hit-dup然后退出。原因cd-hit-auxtools没单独make。 解法进cd-hit-auxtools目录执行make后重跑。症状.clstr文件里序列号被截断拿它回查全库对不上。原因-d默认20描述在20字符或第一个空格处截断。 解法重跑时加-d 0保留完整defline到第一个空格。 下一步CD-HIT的用法收敛想清楚全局还是局部口径、阈值配什么词长、内存给多少预算剩下的都是重复劳动。参数全集以不带参数运行./cd-hit打印的说明为准版本演进看 ChangeLog完整手册见 doc/cdhit-user-guide.tex同目录有PDF版16S与miRNA-seq两条完整流水线在usecases/下可直接复用。引用请标注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），仅供参考