宏基因组功能注释实战：CAZyme与VFDB数据库搭建及DIAMOND比对全流程

宏基因组测序的价格这几年一路走低一个样本几百块钱就能拿到几十个G的原始数据但真正让很多人卡住的从来不是测序本身而是拿到clean data之后的功能注释环节。我见过太多人把数据跑完组装、基因预测之后就停在那里手里攥着一堆ORF序列却不知道下一步该往哪个数据库上比对或者好不容易找齐了数据库跑了两天两夜发现结果文件是空的。CAZyme和VFDB这两个库一个管碳水化合物活性酶一个管毒力因子是环境微生物组和病原微生物组研究里绕不开的两座大山。这篇内容就是把我自己从零搭建这套流程的完整过程拆开来讲包括数据库怎么下、脚本怎么写、参数怎么调、报错怎么排适合刚接触宏基因组功能注释的研究生也适合想把这套流程标准化、跑批量的生信工程师。1. 功能注释的整体设计与数据库选型逻辑1.1 为什么是CAZyme和VFDB这两个库宏基因组功能注释的数据库少说也有几十个KEGG、COG、GO、Pfam、TIGRFAM、CAZyme、VFDB、CARD、MGE等等每个库都有自己的侧重。你不可能把所有库都跑一遍时间和存储都扛不住。选库的核心逻辑是看你的科学问题。CAZyme全称Carbohydrate-Active enZYmes Database专门收录碳水化合物活性酶包括糖苷水解酶GH、糖基转移酶GT、多糖裂解酶PL、碳水化合物酯酶CE、辅助活性酶AA以及碳水化合物结合模块CBM。如果你做的是肠道微生物组、土壤微生物组、堆肥、厌氧消化、木质纤维素降解这类跟多糖代谢强相关的研究CAZyme是必跑的。它能告诉你这个群落里谁在降解纤维素、谁在合成胞外多糖、谁在修剪糖链。VFDB全称Virulence Factor Database收录的是细菌毒力因子。做病原微生物、临床微生物组、环境耐药与毒力传播、食品安全微生物这类方向VFDB是标配。它把毒力因子分成若干大类比如黏附、侵袭、毒素、铁摄取、分泌系统等注释结果能直接支撑“这个环境里存在哪些潜在致病菌、它们携带什么毒力武器”的推论。这两个库一个偏代谢功能一个偏致病潜力组合起来刚好覆盖了环境微生物组研究里最常见的两类问题。而且它们的注释方式都是基于序列比对流程高度相似可以共用一套比对框架工程上很划算。1.2 比对策略的选择BLAST还是DIAMOND早期大家做功能注释基本都用BLASTP准确是准确但速度实在感人。一个几十万条ORF的宏基因组用BLASTP对CAZyme库跑一遍单机跑几天是常态。后来DIAMOND出现速度提升了几百到上千倍灵敏度在默认参数下跟BLASTP相当成了宏基因组功能注释的事实标准。我自己的选择是DIAMOND为主BLASTP只在个别需要极高灵敏度的验证场景下用。DIAMOND的blastp模式支持--sensitive、--more-sensitive、--very-sensitive几档灵敏度宏基因组注释一般用--sensitive就够了追求召回率可以上--more-sensitive但时间会明显增加。这里有个很多人忽略的点DIAMOND建库和比对是两个阶段。建库makedb只需要做一次之后所有样本的比对都复用这个库。CAZyme库建完大概几百MB到1GB出头VFDB更小建库时间都在几分钟级别。建库时加--in指定fasta加--db指定输出前缀生成的.dmnd文件就是后续比对要用的库。1.3 注释结果的判定标准比对完拿到的是比对结果不是注释结果。从比对到注释中间隔着一层判定逻辑这层逻辑直接决定了你最终结果的可靠性。最常用的判定标准是e-value阈值加identity阈值加覆盖度阈值三件套。e-value一般卡1e-5identity卡30%或40%覆盖度卡70%或80%。但具体卡多少要看你的研究目的。如果你做的是新物种新酶的挖掘阈值可以放宽identity卡30%甚至25%覆盖度卡50%宁可要假阳性也不要漏掉潜在的新功能。如果你做的是群落功能的定量比较阈值要收紧identity卡40%以上覆盖度卡80%保证注释到的都是可靠的同源。还有一个关键参数是比对长度。DIAMOND默认输出的是局部比对一条ORF可能只比对上目标蛋白的一小段。如果你不加长度过滤会出现一条100个氨基酸的ORF只比对上10个氨基酸但identity 100%的情况这种注释基本没有生物学意义。我的做法是要求比对长度至少50个氨基酸或者至少达到query长度的70%。注意CAZyme和VFDB的判定阈值不要用同一套。CAZyme库里的蛋白结构域比较保守identity可以适当放宽VFDB里的毒力因子很多是水平转移来的序列差异大identity卡太严会漏掉很多真实的毒力基因。我一般CAZyme用identity 30%、覆盖度70%VFDB用identity 35%、覆盖度70%。2. 数据库下载与本地化部署的实操细节2.1 CAZyme数据库的获取与格式整理CAZyme数据库的官方发布在dbCAN这个平台上dbCAN把CAZyme的注释做了很好的整合提供了现成的HMM库和DIAMOND库。但如果你想要最原始的CAZyme序列还是得从CAZy官网拿。CAZy官网提供的是按家族分类的序列文件每个家族一个fasta比如GH1.fasta、GT2.fasta这样。你需要把这些分散的文件合并成一个总库。合并的时候有个细节每个家族的序列ID可能重复合并前要加上家族前缀否则后面注释结果里分不清是哪条序列比上的。我常用的合并脚本逻辑是这样的遍历所有家族的fasta文件对每条序列的header加上家族名作为前缀然后追加到一个总文件里。用awk或者python都行python更可控。合并完之后用DIAMOND建库建库命令里加--threads指定线程数一般给到CPU核心数的80%就行留一点给系统。dbCAN平台还提供了一个预构建的DIAMOND库叫dbCAN-HMMdb或者dbCAN-DIAMOND直接下载就能用省去自己合并建库的麻烦。但预构建库的版本更新可能滞后如果你需要最新版的CAZyme注释还是自己从CAZy官网拉序列自己建库更靠谱。2.2 VFDB数据库的下载与去冗余处理VFDB的官方下载在微生物所的网站上提供的是FASTA格式的蛋白序列文件名叫VFDB_setA_nt.fas或者VFDB_setA_pro.fassetA是经过实验验证的毒力因子setB是预测的。做注释一般用setA就够了setB假阳性偏高。VFDB的原始序列有个问题同一个毒力因子在不同菌株里可能有多个拷贝序列高度相似但不完全一样。直接建库比对一条ORF可能比对上十几条高度相似的参考序列结果冗余严重。我的做法是先做一步去冗余用CD-HIT把identity 95%以上的序列聚成一类每类取一条代表序列。这样库的大小能压缩30%到50%比对速度提升明显注释结果也更干净。去冗余的命令大概是cd-hit -i VFDB_setA_pro.fas -o VFDB_setA_nr.fas -c 0.95 -n 5 -M 16000 -T 8。-c 0.95是identity阈值-n 5是词长-M是内存限制-T是线程数。去冗余完之后再建DIAMOND库。提示VFDB的序列header里包含了毒力因子的分类信息比如VFG编号、菌属、毒力因子名称。建库前最好把这些信息整理成一个独立的注释表后面注释结果出来可以直接关联省得再去解析header。2.3 数据库版本管理与路径规划数据库这东西版本一多就容易乱。我见过有人跑完分析半年后想复现结果发现数据库更新了同样的脚本跑出来的结果对不上。所以从第一天起就要做好版本管理。我的做法是给每个数据库建一个带日期的目录比如/db/CAZyme_20240115/、/db/VFDB_20240115/目录里放原始序列、去冗余序列、DIAMOND库文件、版本说明文件。版本说明文件里记录下载日期、来源URL、序列条数、去冗余参数。跑分析的时候在脚本里硬编码数据库路径不要用环境变量或者软链接避免以后路径变了找不到。存储方面CAZyme原始库大概几百MB建完DIAMOND库1GB左右VFDB原始库几十MB去冗余后更小DIAMOND库也就几十MB。两个库加起来不到2GB对存储压力不大。但如果你还要跑KEGG、CARD这些大库就得提前规划好磁盘空间。2.4 依赖工具的安装与版本确认这套流程依赖的核心工具就三个DIAMOND、CD-HIT、以及一个能跑python脚本的环境。DIAMOND用conda装最省事conda install -c bioconda diamond装完用diamond --version确认版本。CD-HIT同样conda install -c bioconda cd-hit。这里有个坑DIAMOND的版本更新比较频繁不同版本之间makedb生成的库文件格式可能不兼容。如果你用A版本建的库用B版本去比对可能报“database version mismatch”之类的错。所以建库和比对一定要用同一个DIAMOND版本。我的习惯是在项目目录里记录DIAMOND版本号比如diamond v2.1.9下次跑之前先确认版本一致。python环境建议用conda单独建一个装好pandas、biopython这些常用库。脚本里读写fasta用biopython的SeqIO最稳解析比对结果用pandas最方便。如果你不想装biopython用awk也能处理但代码可读性差很多。3. 从ORF到注释结果的完整实操流程3.1 输入数据的准备与质量检查这套流程的输入是基因预测之后的蛋白序列通常是Prodigal或者MetaGeneMark跑出来的.faa文件。在开始注释之前有几项检查必须做。第一确认序列条数和总长度。一个典型的宏基因组样本Prodigal预测出来的ORF大概在几十万到几百万条之间总氨基酸数在几千万到几亿。如果序列条数明显偏少可能是基因预测参数不对或者输入contig质量太差。第二检查序列ID是否唯一。Prodigal默认输出的ID是递增编号一般不会重复但如果你把多个样本的faa合并了就可能出现ID冲突。合并前给每个样本的ID加上样本前缀比如sample1_1、sample1_2这样。第三去除过短的序列。长度小于30个氨基酸的ORF比对上的概率很低留着只会增加计算量。用seqkit或者自己写脚本过滤一下seqkit seq -m 30 input.faa filtered.faa简单直接。第四确认没有非蛋白字符。有时候基因预测会输出带*或者X的序列这些字符会干扰DIAMOND比对。用seqkit grep -v -s -p X或者类似命令过滤掉。3.2 DIAMOND比对的核心参数与运行比对是整个流程里最耗时的环节参数设置直接决定运行时间和结果质量。我以CAZyme库为例完整的比对命令是这样的diamond blastp \ --query filtered.faa \ --db /db/CAZyme_20240115/CAZyme.dmnd \ --out CAZyme_blastp.tsv \ --outfmt 6 \ --evalue 1e-5 \ --max-target-seqs 1 \ --sensitive \ --threads 32 \ --block-size 4.0 \ --index-chunks 1逐项解释一下。--outfmt 6是标准tabular输出包含qseqid、sseqid、pident、length、mismatch、gapopen、qstart、qend、sstart、send、evalue、bitscore这12列。--max-target-seqs 1表示每条query只保留最优的一条比对结果避免输出文件过大。如果你需要看多条比对结果做后续分析可以设成5或者10。--sensitive是灵敏度档位比默认的--fast慢一些但召回率更高。--block-size 4.0是内存块大小单位是GB根据你的机器内存调整给到内存的60%到70%比较稳妥。--index-chunks 1在内存充足时能加快速度内存紧张时可以设成2或4。运行时间参考一个50万条ORF的样本32线程对CAZyme库跑--sensitive大概需要30到60分钟。VFDB库小很多同样条件10到20分钟就能跑完。如果你有几十个样本建议写个循环脚本批量跑或者用GNU parallel并行。3.3 比对结果的过滤与注释判定比对完拿到的tsv文件不能直接用需要经过过滤和判定才能变成注释结果。我一般用python脚本处理核心逻辑分三步。第一步读入比对结果按qseqid分组每组取bitscore最高的那条。虽然--max-target-seqs 1已经做了初步筛选但保险起见还是再确认一下。第二步应用阈值过滤。对CAZyme我要求pident 30且length 50且evalue 1e-5。对VFDB要求pident 35且length 50且evalue 1e-5。不满足的比对直接丢弃对应的ORF标记为未注释。第三步关联注释信息。CAZyme的sseqid里包含了家族信息比如GH1|ABC123用split取第一部分就是家族。VFDB的sseqid需要关联之前整理的注释表拿到毒力因子分类和名称。处理完的注释结果输出成两张表一张是ORF级别的详细注释包含ORF ID、比对上的参考序列、identity、覆盖度、e-value、注释分类另一张是汇总表统计每个样本里各CAZyme家族或各VFDB分类的丰度。3.4 丰度计算与结果可视化准备功能注释的最终目的通常是做丰度比较所以注释完之后要算丰度。丰度的计算方式有两种一种是基于ORF的计数就是某个功能注释到多少条ORF另一种是基于ORF的丰度加权需要先把reads回贴到ORF上拿到每个ORF的丰度再按功能汇总。第一种方式简单但粗糙适合快速看趋势。第二种方式准确但多一步回贴适合正式发表的数据。我一般两种都算先看计数结果有没有明显异常再用加权结果做最终分析。回贴用bowtie2或者salmon都行。bowtie2需要先建ORF的索引然后把每个样本的clean reads贴回去用--no-unal去掉未比对的reads最后用samtools或者featureCounts统计每个ORF的reads数。salmon更省事直接对ORF序列建索引然后quant输出TPM和counts。拿到ORF丰度之后按注释表把ORF映射到功能分类汇总得到每个样本的功能丰度谱。这个丰度谱就可以拿去做PCoA、LEfSe、随机森林这些下游分析了。4. 脚本优化与批量处理的工程化实践4.1 把零散命令封装成可复用脚本手动敲命令跑一两个样本还行样本一多就容易出错。我的做法是把整个流程拆成几个独立的脚本每个脚本负责一个环节用主控脚本串起来。第一个脚本负责数据库准备包括下载、去冗余、建库输出是.dmnd文件和注释表。第二个脚本负责比对输入是faa文件和数据库路径输出是tsv。第三个脚本负责过滤和注释判定输入是tsv输出是注释表和汇总表。第四个脚本负责丰度计算输入是注释表和reads输出是丰度谱。每个脚本都接受命令行参数比如-i输入、-o输出、-d数据库路径、-t线程数。这样换一个样本只需要改参数不用改代码。脚本头部加上set -e任何一步出错就停止避免错误累积。4.2 并行化与资源调度宏基因组样本多的时候串行跑太慢。有两种并行策略样本间并行和样本内并行。样本间并行就是同时跑多个样本每个样本分配一定的线程。比如你有32核可以同时跑4个样本每个样本8线程。用GNU parallel实现起来很简单ls *.faa | parallel -j 4 diamond blastp --query {} --db CAZyme.dmnd --out {.}.tsv --threads 8 --sensitive样本内并行就是DIAMOND自身的多线程--threads给到32单个样本跑满所有核心。两种策略选哪种取决于你的样本数量和机器配置。样本少就样本内并行样本多就样本间并行。如果是在集群上跑建议用SLURM或者PBS提交任务每个样本一个job资源申请写清楚。DIAMOND对内存的需求跟库大小和block-size有关CAZyme库跑--sensitive大概需要8到16GB内存申请的时候留够余量。4.3 日志记录与中间文件管理跑批量分析最怕的就是跑到一半出错不知道错在哪也不知道哪些样本跑完了。所以日志和中间文件管理必须做好。我的做法是每个样本一个目录目录里放原始faa、比对tsv、过滤后的注释表、日志文件。日志文件记录开始时间、结束时间、运行的命令、消耗的资源、是否有报错。主控脚本每跑完一个样本就往总日志里追加一行记录样本名和状态。中间文件要不要保留看情况。比对tsv文件比较大一个样本可能几百MB到几GB如果磁盘紧张注释表生成之后可以删掉tsv。但建议至少保留一份压缩的tsv万一后面要重新过滤或者换阈值不用重新比对。4.4 结果校验与异常样本识别批量跑完之后不能直接拿结果去分析要先做一轮校验。校验的内容包括每个样本的ORF条数是否在合理范围、注释率是否正常、注释到的功能分类分布是否合理。注释率是个关键指标。CAZyme的注释率一般在5%到15%之间VFDB的注释率一般在1%到5%之间。如果某个样本的注释率明显偏离这个范围比如CAZyme注释率只有1%或者高达30%就要查原因。注释率过低可能是ORF质量差或者数据库不匹配注释率过高可能是污染或者阈值太松。功能分类分布也要看。CAZyme里GH和GT通常是丰度最高的两类如果某个样本里AA或者CBM异常高可能是样本特性也可能是注释错误。VFDB里如果某个毒力因子分类异常富集要结合样本背景判断是不是真实的生物学信号。注意异常样本不一定要剔除但一定要标记出来在后续分析里单独说明。我遇到过注释率偏低的样本查下来是测序深度不够导致ORF预测不完整这种样本在丰度比较时权重应该降低。5. 常见报错与排查技巧实录5.1 DIAMOND建库与比对阶段的典型报错报错一Error: Database file is corrupt or incomplete这个报错通常是建库过程中断了或者磁盘空间不足导致.dmnd文件没写完整。解决办法是删掉重建建库前确认磁盘剩余空间至少是库大小的两倍。建库时加--verbose能看到进度方便判断是不是卡住了。报错二Error: Sequence type mismatchDIAMOND建库时如果用的是蛋白序列比对时query也必须是蛋白序列。如果你拿核酸序列去比对蛋白库就会报这个错。检查一下faa文件是不是真的蛋白序列有时候基因预测输出的是fna核酸需要先翻译。报错三Error: Out of memoryDIAMOND比对时内存不够通常是--block-size设太大了。CAZyme库跑--sensitiveblock-size设4.0GB实际内存占用可能到10GB以上。如果机器内存只有16GB把block-size降到2.0或者1.0或者用--index-chunks 2减少内存峰值。报错四比对结果为空tsv文件一行都没有说明没有任何ORF比对上。可能的原因数据库路径写错了、query文件是空的、e-value阈值太严、序列类型不匹配。排查顺序是先确认数据库文件存在且大小正常再确认query文件非空然后放宽e-value到1e-3试试最后检查序列类型。5.2 注释判定阶段的逻辑陷阱陷阱一一条ORF比对上多个家族CAZyme里有些蛋白是多功能酶同时属于两个家族。如果你的判定逻辑是取bitscore最高的那条可能会漏掉另一个家族。处理方式是允许一条ORF对应多个注释在汇总时分别计数。但要注意去重同一条ORF对同一个家族的多条参考序列只算一次。陷阱二覆盖度计算方式不统一覆盖度有两种算法一种是比对长度除以query长度一种是比对长度除以subject长度。宏基因组注释一般用前者因为query是我们的ORF我们关心的是这条ORF有多少比例比上了。但有些流程用后者导致结果不可比。写脚本的时候一定要在注释里写清楚用的是哪种。陷阱三e-value的科学计数法解析DIAMOND输出的e-value是科学计数法比如1.2e-45。用python的float()能直接解析但用awk的时候要注意awk对科学计数法的支持因版本而异。稳妥的做法是用python处理或者用printf转换。5.3 批量运行时的资源冲突与解决批量跑几十个样本最容易出的问题是资源冲突。多个DIAMOND进程同时读写同一个数据库文件虽然DIAMOND支持并发读但IO压力会很大。解决办法是把数据库文件放在SSD上或者用--index-chunks控制内存映射方式。另一个问题是临时文件冲突。DIAMOND默认会在当前目录生成临时文件如果多个进程在同一个目录跑临时文件可能互相覆盖。解决办法是每个样本在独立目录里跑或者用--tmpdir指定不同的临时目录。还有磁盘写满的问题。比对tsv文件累积起来很占空间几十个样本可能上百GB。建议跑完一个样本就压缩一个或者定期清理中间文件。用gzip压缩tsv压缩比大概5到10倍能省不少空间。5.4 常见问题速查表问题现象可能原因排查方法解决方案建库报错corrupt磁盘满或中断检查磁盘空间和文件大小删掉重建确保空间充足比对结果为空路径错/类型不匹配/阈值严逐项检查路径、序列类型、e-value修正路径确认蛋白序列放宽阈值内存溢出block-size过大监控内存占用降低block-size或增加index-chunks注释率异常低ORF质量差或库不匹配检查ORF长度分布和数据库版本过滤短ORF更新数据库注释率异常高污染或阈值过松检查样本来源和阈值设置收紧阈值排查污染批量运行卡住IO冲突或资源竞争查看进程状态和IO等待分离目录限制并发数结果无法复现数据库版本变了对比数据库版本记录固定数据库版本记录版本号6. 流程扩展与下游分析衔接6.1 从注释结果到群落功能比较拿到CAZyme和VFDB的丰度谱之后下一步通常是做群落间的功能比较。CAZyme的丰度谱可以按家族汇总也可以按EC号汇总看不同处理组之间哪些酶家族显著差异。VFDB的丰度谱可以按毒力因子分类汇总看不同环境里毒力因子的分布差异。比较的方法跟物种组成分析类似PCoA看整体差异PERMANOVA检验组间差异是否显著LEfSe找标志性功能随机森林找关键预测因子。不同的是功能丰度谱通常更稀疏很多功能在部分样本里为0做统计检验时要考虑零膨胀的问题。我一般先用vegan包的adonis2做PERMANOVA再用microbiomeMarker或者LinDA做差异分析。CAZyme的差异结果可以映射到CAZy官网的家族页面上看这些差异家族对应的底物和反应类型帮助解释生物学意义。6.2 与其它功能数据库的联合分析CAZyme和VFDB只是功能注释的一部分实际研究里往往还要联合KEGG、CARD、MGE等库一起分析。联合分析的关键是统一ORF ID和注释格式让不同库的结果能对应到同一条ORF上。我的做法是建一张总表每行是一条ORF列包括ORF ID、CAZyme注释、VFDB注释、KEGG注释、CARD注释等。这样一条ORF如果同时注释到CAZyme和VFDB就能直接看出来方便做共现分析。比如某些CAZyme基因和毒力因子基因在同一个ORF上可能暗示水平转移事件。联合分析还能做功能模块的共丰度网络。把CAZyme家族和VFDB分类作为节点样本间的丰度相关性作为边构建共丰度网络找功能模块。这种分析能揭示功能之间的协同关系比如某些多糖降解酶和特定毒力因子是否倾向于同时出现。6.3 流程的容器化与可移植性如果你想让这套流程在别的机器上也能跑或者想分享给合作者容器化是最省事的方案。用Docker或者Singularity把DIAMOND、CD-HIT、python环境和脚本打包成一个镜像别人拉下来就能跑不用折腾依赖。Dockerfile的写法大概是基础镜像用continuumio/miniconda3然后conda install装DIAMOND和CD-HITpip install装pandas和biopython把脚本复制进去设置好入口点。构建完的镜像大概1到2GBpush到镜像仓库别人docker pull就能用。如果是在集群上跑Singularity更合适因为集群通常不允许Docker的root权限。Singularity可以直接把Docker镜像转成sif文件用singularity exec运行。数据库文件不建议打进镜像太大了用--bind挂载外部目录就行。6.4 性能调优的几个实战技巧最后分享几个我踩过坑之后总结的性能调优技巧。第一DIAMOND建库时加--no-unlink避免建库过程中反复删除临时文件能快不少。这个参数在磁盘IO慢的机器上效果明显。第二比对时用--compress 1压缩输出tsv文件能小一半以上后续读入也快。但注意压缩后的文件不能直接用文本工具看要用diamond view解压或者用python的gzip模块读。第三如果样本量特别大可以考虑先把所有样本的ORF合并去冗余用非冗余ORF集去比对然后把比对结果映射回各样本。这样比对次数从N次降到1次速度提升巨大。去冗余用CD-HITidentity 100%或者99%保留映射关系就行。第四数据库可以按需拆分。CAZyme库如果只关心GH和GT可以把其他家族去掉库小一半比对快一倍。VFDB如果只关心毒素类也可以只保留相关序列。这种定制化库在特定研究里很实用。这套流程我从最早的BLASTP版本一路迭代到现在的DIAMOND加脚本化版本中间踩过的坑基本都写在这里了。数据库下载和建库是一次性的工作脚本写好之后跑批量就是改改参数的事。真正花时间的是理解每个参数背后的逻辑以及根据你的数据特点调整阈值。我个人的体会是功能注释没有一套放之四海皆准的参数先跑一个小样本试参数看注释率和功能分布是否合理再放大到全部样本比一上来就跑全量要稳妥得多。