做致病菌基因组分析的朋友应该都有过这种体验手上拿到一株菌的组装结果老板第一句话就问“它有哪些毒力因子有没有可能致病”这时候你急需一套快速、可复现、结果还能写进文章的注释流程。毒力因子注释最常用的办法就是把预测出来的蛋白序列跟专门的毒力因子数据库做比对。而在比对工具里DIAMOND凭借“比BLASTP快几百倍、灵敏度还不差”的优势成了我这两年跑批量样本的首选。这篇文章我就以2024年6月整理的这套流程为底子完整讲讲怎么用DIAMOND把致病菌蛋白序列比对到VFDB数据库上完成毒力因子注释包括工具选型、数据库下载、建库比对、结果过滤和常见坑内容比较多但每一步都能直接照着操作。这套流程适合谁用主要面向做细菌基因组、宏基因组、比较基因组或者致病机制研究的同学。你要是只会点Linux基础命令、能连服务器跑跑conda那这个流程你完全吃得下来。哪怕你之前只用过网页版BLAST只要愿意跟着敲命令也能在半天内把整套注释流程跑通。读完你会得到一个纯文本格式的注释结果表每一行是一个基因的比对信息哪个基因命中了哪个毒力因子、相似度多少、覆盖度多少、E值多少清清楚楚后续要画图、要做统计、要写方法部分都有据可依。1. 为什么是DIAMOND加VFDB选型思路与核心优势1.1 VFDB数据库到底存的是什么VFDB全称Virulence Factors of Pathogenic Bacteria是由中国军事医学科学院研发维护的致病菌毒力因子数据库在病原微生物研究领域用得非常多。它把已经验证过的毒力因子按功能分成几大类包括黏附、侵袭、毒素分泌系统、铁摄取系统、荚膜合成、脂多糖修饰、蛋白酶等。每个条目都带着基因名、功能描述、分类信息和对应的蛋白序列直接下载FASTA格式就能本地使用。VFDB官网提供了两套蛋白序列文件一个叫VFDB_setA_pro.fas是全量数据集另一个叫VFDB_setB_pro.fas是核心数据集。这两套的区别很关键setB是经过人工筛选、在公开数据库里有明确实验证据支持的核心毒力因子序列数量更少、去冗余做得更干净setA则把一些推测性、完整基因组预测出来的候选因子也收了进来序列数量明显更多。日常做注释我一般建议默认用setB尤其是写文章的时候核心集的结果更经得起审稿人追问。如果你是想做探索性筛选想尽量多找潜在因子再换成setA跑一版做对比。1.2 为什么不用BLASTP而选DIAMOND经典做法是用BLASTP把预测蛋白跟VFDB蛋白库做比对BLASTP本身结果很可靠但慢也是真慢。一次细菌全基因组预测大概产出4000到6000条蛋白序列跟VFDB里几千条序列做两两比对单线程跑BLASTP常常要几十分钟甚至更久。如果是宏基因组拼接出来的数据蛋白序列动辄几十万条那BLASTP基本要跑到天荒地老。DIAMOND解决的就是这个性能问题。它的核心思路是把查询序列和数据库序列都转成词袋式的种子索引seed index再用双索引double indexing策略快速定位可能比对的区域最后在候选区域上做扩展和动态规划比对。说白了它用“先粗筛、再精对齐”的策略大幅减少了真正需要做动态规划比对的序列对数加上底层利用了SIMD指令集做并行加速所以在保证接近BLASTP灵敏度的前提下比BLASTP快两三个数量级。DIAMOND默认的比对模式其实偏向fast模式追求极致的速度但灵敏度比BLASTP略低。如果你希望结果更接近BLASTP推荐把比对模式调成--sensitive或者--very-sensitive。实际测下来我拿一株肠炎沙门氏菌的蛋白序列跑VFDB注释BLASTP大概需要40分钟左右DIAMOND的sensitive模式只要30秒不到结果几乎能对上这个效率差距在做批量样本时就是天壤之别。速度对比可以看下面这个表数字来自我本地服务器上的实测不同机器会略有差异但趋势基本一致。工具与参数单样本耗时约灵敏度相对BLASTP适用场景BLASTP默认参数30-50分钟基准单样本、追求完全一致DIAMOND默认fast5-10秒略低超大规模初筛DIAMOND --sensitive20-40秒很接近常规注释首选DIAMOND --very-sensitive1-3分钟几乎一致样本量不大但要求高1.3 现成工具和自建流程怎么选现在也有不少集成好的工具能做毒力因子注释比如ABRicate、staramr、VFDB的在线比对服务。ABRicate这类工具确实方便一条命令就能把多个数据库跑完但问题在于它把比对参数、数据库处理方式都封装在黑盒里了你想调整一个阈值、换一个数据库版本就得去翻它的源码。在线服务则受限于上传文件大小和数据安全你也不太可能在网页上批量提交几百个基因组做注释。自建DIAMOND加VFDB流程表面上看要自己写的命令多一些但数据库本地化之后可以反复使用参数也完全可控产出结果还能直接放进自己的流程管道里是性价比最高的方案。2. 实操前准备安装DIAMOND并下载VFDB数据2.1 DIAMOND的安装方式DIAMOND的安装非常简单二进制文件直接下载解压就能跑官方GitHub的release页面提供了编译好的Linux、macOS和Windows版本。如果你服务器上有conda我更推荐用conda装方便管理和切换版本命令是这样的conda create -n diamond python3.9 -y conda activate diamond conda install -c bioconda diamond -y diamond --version不想用conda的话官方二进制下载方式也不难先到GitHub release页面找到diamond-linux64.tar.gz这个文件然后执行wget https://github.com/bbuchfink/diamond/releases/download/v2.1.9/diamond-linux64.tar.gz tar xzf diamond-linux64.tar.gz sudo mv diamond /usr/local/bin/ diamond --version这里有个小提醒DIAMOND版本之间数据库格式不完全兼容。你如果以前建过旧版本的DIAMOND数据库升级软件之后最好把数据库重新建一遍否则会报数据库格式过旧之类的错误。我自己就踩过这个坑升级之后忘了重建库跑一个批量任务到一半全部中断才发现是版本不匹配。2.2 下载VFDB数据库文件VFDB数据库可以通过FTP或HTTP直接下载。官网地址是www.mgc.ac.cn/VFs/下载区有明确链接。因为数据库文件名偶尔会有变化我建议你进官网的Download页面找到“VFDB_setB_pro.fas.gz”和“VFDB_setA_pro.fas.gz”这两个文件的直链。命令行下载可以这样写wget http://www.mgc.ac.cn/VFs/Down/VFDB_setB_pro.fas.gz wget http://www.mgc.ac.cn/VFs/Down/VFDB_setA_pro.fas.gz gunzip *.gz下载完记得检查一下文件完整性看大小是否跟官网标注一致。解压之后把两个文件重命名成好记的名字比如vfdb_core.fa和vfdb_full.fa后面用起来方便。再说一下文件内容VFDB的FASTA头信息里包含了完整的毒力因子描述比如“VF0123 (aerobactin biosynthesis protein IucA) [Escherichia coli]”后面的括号里是基因名和功能描述方括号里是物种来源。这些信息在后面做结果解读的时候直接能用不需要额外查表。2.3 准备待注释的蛋白序列VFDB只有蛋白序列数据库所以DIAMOND比对用的是blastp流程严格意义上DIAMOND对应blastp的算法叫DIAMOND blastp。你拿到的输入序列必须是蛋白序列FASTA不是基因组DNA序列。蛋白序列从哪里来如果是你自己组装的细菌基因组通常是从注释结果里提取用Prokka或Bakta跑完基因组注释后会生成一个后缀为.faa的文件那个就是蛋白序列文件。如果你手里只有基因组FASTA想快速拿到蛋白序列可以装个Prokka一条命令搞定conda install -c conda-forge -c bioconda prokka -y prokka input_genome.fasta --outdir prokka_out --prefix sample1 --cpus 8跑完之后prokka_out/sample1.faa就是拿来注释毒力因子的输入文件。这个步骤的作用是把基因组上的编码基因预测出来并翻译成蛋白没有这一步后面DIAMOND就没法比对。宏基因组样本也一样你先用Prokka里的prodigal或者单独跑Prodigal把拼接contig上的基因预测出来翻译成蛋白序列再做比对。拿到.faa文件后最好先用seqkit统计一下序列条数确认文件不是空的、序列没有明显截断conda install -c bioconda seqkit -y seqkit stats sample1.faa正常细菌基因组注释得到的蛋白序列数量在3000到6000条之间太少或者太多都值得怀疑一下。一次我帮同事排查跑出来的蛋白序列只有几十条后来发现他给的“基因组”其实只是一个质粒序列这种低级错误在数据预处理阶段就能拦住。3. 从建库到比对核心流程一步步来3.1 用VFDB构建DIAMOND格式数据库DIAMOND比对前要把FASTA格式的数据库转成它的二进制索引格式这一步用makedb子命令完成。diamond makedb --in vfdb_core.fa --db vfdb_core--in指定输入的FASTA文件--db指定输出数据库的前缀名。运行结束之后会生成一个后缀为.dmnd的文件比如vfdb_core.dmnd这个就是DIAMOND专用的数据库文件。建库过程很快VFDB核心集几千条序列秒级到十几秒就完成了。这里需要特别注意makedb只接受蛋白序列FASTA。如果你不小心把核酸序列丢进去DIAMOND会报错或者建出一个无意义的库。判断方法也简单看FASTA里的序列是ACGT组成的还是氨基酸字母组成的一眼就能看出来。另外建库前建议把FASTA里的注释头信息处理干净不要有空格以外的特殊符号个别数据库文件头信息里如果带了过多空格和特殊符号可能会影响后面输出结果的解析。VFDB官方文件头信息一般没问题但如果你以后自己构建其他数据库比如从NCBI下载的nr序列建库最好先用seqkit或awk把序列描述统一简化一下免得结果表里sseqid那一列长得无法直视。3.2 DIAMOND blastp比对命令与参数详解数据库建好之后进入核心步骤用DIAMOND把查询蛋白序列比对到VFDB数据库上。diamond blastp \ --db vfdb_core \ --query sample1.faa \ --out sample1_vfdb.tsv \ --outfmt 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore \ --evalue 1e-5 \ --id 80 \ --query-cover 80 \ --subject-cover 50 \ --max-target-seqs 15 \ --sensitive \ --threads 16 \ --block-size 8.0 \ --tmpdir /tmp/这几个参数我逐个解释一下。--outfmt 6是DIAMOND的tabular输出格式每一行是一个比对结果默认格式包含12列。手动指定qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore这12个字段是为了确保输出顺序稳定不至于因为版本升级导致列顺序变化。这个如果你写自动化脚本处理结果建议每次都手动写全字段。--evalue设的是E值阈值1e-5是序列比对里常用的默认值含义是期望值小于10的负5次方代表这个匹配在随机情况下出现的概率极低。对毒力因子注释来说E值阈值卡1e-5基本够了不用太严后面还有identity和coverage来把关。--id和--query-cover分别是序列一致性和查询覆盖率阈值。这里我习惯直接设到80意思是至少80%的序列一致、命中区域覆盖查询序列80%以上才算有效比对。这在细菌毒力因子注释里是比较常见且稳妥的阈值。你如果想保守一点可以提到90。但这里有个讲法不同毒力因子基因存在同源基因同一家族成员之间序列相似度就可能刚好在70-85%之间阈值卡太死会把一些真实但变异的毒力因子漏掉阈值太松又容易混进一堆同源蛋白造成假阳性。最佳策略是先宽松过滤比如identity设60拿到全部候选再根据后续分析需求决定要不要收紧。--subject-cover 50是比对区域在数据库序列上的覆盖率这个参数很多教程不提但VFDB数据库里有些毒力因子只有部分结构域被注释查询序列的比对区域在数据库序列上覆盖很低这类结果其实不太可靠所以加一个50%的subject覆盖度过滤比较合理。--max-target-seqs控制每条查询序列最多输出多少个比对结果。默认是25做注释的话设15够用了因为一条蛋白序列通常只对应一个毒力因子输出太多同源匹配只会让结果表变得冗长。不过要注意这个参数是“最多保留多少个最优的比对”跟后面说的“一条序列被多个数据库条目命中”是两个维度别混淆。--sensitive是运行模式前面也提到这个模式能显著提高灵敏度代价是速度变慢但即便是sensitive模式依然比BLASTP快上百倍。--threads设CPU线程数建议根据服务器核心数来定能用多少给多少。--block-size是内存管理参数单位是GB意思是每轮比对最多用多少内存来装载数据库块设太大容易内存溢出设太小会导致数据库被反复读取、速度下降一般设4到8比较合适默认是2如果你的服务器内存有64GB以上可以调到8甚至更高。--tmpdir指定临时目录默认是/tmp注意这个目录空间要够大比对过程中会产生临时文件如果/tmp满了任务会异常退出服务器上跑的话建议指定到一个空间充足的大盘路径。跑完之后结果文件sample1_vfdb.tsv就生成了后面所有分析都是基于这个文件。这个过程大概消耗多少时间一株沙门氏菌的蛋白序列16线程跑sensitive模式整个过程不到1分钟就能跑完非常爽快。3.3 结果文件长什么样用head命令看一眼结果文件head -5 sample1_vfdb.tsv你会看到类似这样的内容prokka_00001 VF0123 98.754 320 4 0 1 320 45 364 1.2e-180 963.4 prokka_00002 VF0456 76.512 410 88 4 10 419 2 411 3.4e-120 655.2 prokka_00003 VF0789 88.235 85 10 0 5 89 330 414 5.6e-45 268.1这12列的含义我做了一张速查表建议你收藏一下列名含义qseqid查询序列ID一般是你的蛋白序列编号sseqid数据库序列IDVFDB里的毒力因子编号pident序列一致性百分比length比对区域的长度mismatch错配数gapopen缺口开位数qstart查询序列比对起始位置qend查询序列比对结束位置sstart数据库序列比对起始位置send数据库序列比对结束位置evalue期望值bitscore比对得分我想特别提醒的是sstart/send和qstart/qend这四列它们决定了比对区域的位置和方向。如果数据库序列是反向互补链上的基因sstart会大于send这本身不是错误但在计算覆盖率时注意用绝对值。做后续分析的时候很多人会忘记检查这个方向问题导致“基因组区间坐标”算错最后画图时对不上。4. 结果过滤、后处理与毒力因子注释解读4.1 如何科学设定过滤阈值DIAMOND输出的原始结果表虽然已经按你设置过的参数过滤过一轮了但实际分析中建议再做一轮严格的后过滤。原因很简单比对工具的过滤参数是全局通用的但VFDB注释结果的可靠性还受数据库冗余度、输入序列质量的影响只靠比对参数一关可能挡不住一些边缘情况。我自己常用的后置过滤器用awk就能实现筛选条件可以结合pident和query-cover比如这样awk -F \t $3 90 $4 100 你的结果文件.tsv假设$3是pident第三列、$4是length第四列$4在这里代表比对长度如果要按覆盖度过滤就得额外计算覆盖率或者用DIAMOND的outfmt扩展参数直接输出qcovhsp和scovhsp。DIAMOND在outfmt里支持这些扩展字段简化后的命令如下diamond blastp \ --db vfdb_core \ --query sample1.faa \ --out sample1_vfdb_qcov.tsv \ --outfmt 6 qseqid sseqid pident length evalue bitscore qcovhsp scovhsp \ --evalue 1e-5 \ --sensitive \ --threads 16这里qcovhsp是查询序列上最高得分比对区域的覆盖率scovhsp是数据库序列上最高得分比对区域的覆盖率。拿到这两个字段之后过滤命令就变得很直接awk -F \t $3 90 $7 90 $8 50 sample1_vfdb_qcov.tsv sample1_vfdb_filtered.tsv这里$3是pident$7是qcovhsp$8是scovhsp。我给的经验阈值是pident 90、qcovhsp 90、scovhsp 50。为什么scovhsp要求低一点因为VFDB里很多毒力因子序列并不是全长的数据库里的条目可能只覆盖基因的一部分强行要求数据库序列全覆盖会丢掉有效结果。筛选阈值这里多说一句如果分析的是高度保守的毒力因子标志基因比如肠杆菌科的侵袭蛋白AInvA阈值设到90以上完全没有问题。但如果你关注的是那些高频重组的表面蛋白比如某些黏附素、菌毛亚基序列差异本身就比较大这种时候把pident卡到90会漏掉很多东西。所以阈值设置一定要结合你研究的物种和关注的功能来定不要拿一套参数走天下。4.2 一因多果与最优匹配提取DIAMOND输出的结果里一条查询蛋白序列可能会匹配到VFDB里的多条毒力因子。原因有两类一是同一毒力因子在不同菌株里以相近的序列被收录数据库里存在自然冗余二是查询蛋白本身是一个多结构域蛋白不同结构域分别匹配到了不同毒力因子的条目。这时候就需要做“去冗余保留最优匹配”的处理。最简单的做法是用sort加awk按qseqid去重保留每条查询序列的第一次出现记录前提是结果已经按比对得分或者E值排序过。DIAMOND默认输出是按查询序列和得分排序的如果你不放心可以先按bitscore降序排好再取第一行sort -t $\t -k1,1 -k12,12nr sample1_vfdb_filtered.tsv \ | awk -F \t !seen[$1] sample1_vfdb_besthit.tsv这几行的逻辑是先按第一列查询序列ID排序再按第12列bitscore数值降序排序然后用awk做去重同一个查询序列ID第一次出现时打印之后出现则跳过得到的文件里每个查询序列只有一条最优比对记录。还有一种更精细的做法是用python脚本处理尤其是你后面要做大量样本批量统计时python处理起来更灵活from collections import defaultdict best_hits defaultdict(list) with open(sample1_vfdb_filtered.tsv) as fh: for line in fh: cols line.strip().split(\t) qid cols[0] bitscore float(cols[11]) best_hits[qid].append((bitscore, line.strip())) with open(sample1_vfdb_besthit.tsv, w) as out: for qid, hits in best_hits.items(): best sorted(hits, keylambda x: x[0], reverseTrue)[0] out.write(best[1] \n)这段代码会读入过滤后的结果文件按bitscore排序后为每个查询序列输出一个最佳匹配。之后你再把这个文件跟VFDB的注释信息表做关联就能知道每个基因对应什么毒力因子了。4.3 将VFDB编号映射为功能注释信息得到最佳匹配列表之后当前文件里只有VFDB的编号比如VF0123。要做注释解读还得把编号翻译成基因名、功能描述和毒力因子分类。VFDB官网上有下载“VFs.xls”或者“VirulenceFactors.xlsx”这类注释表字段包括VFID、GeneName、Function、Category等。你也可以从VFDB数据库的FASTA头信息里解析比如用grep提取基因名和描述。这里我提供一个简单实用的步骤第一步下载VFDB的注释表第二步用python字典把VFDB编号映射到功能描述第三步输出一份容易读的注释结果。VFDB官网的Download页面上有一个“VFs.xls.gz”或者类似的注释文件下载解压后是Excel或制表符分隔的文本。如果是Excel建议用python的pandas读取然后导出成TSVimport pandas as pd df pd.read_excel(VFs.xls, header1) df.to_csv(vfdb_annot.tsv, sep\t, indexFalse)然后写一个小脚本把比对结果跟注释表合并annot {} with open(vfdb_annot.tsv) as fh: header fh.readline().strip().split(\t) for line in fh: cols line.strip().split(\t) annot[cols[0]] dict(zip(header, cols)) with open(sample1_vfdb_besthit.tsv) as fh, open(sample1_vfdb_annotated.tsv, w) as out: for line in fh: cols line.strip().split(\t) vfid cols[1] if vfid in annot: info annot[vfid] out.write(\t.join([cols[0], vfid, info.get(GeneName, ), info.get(Function, ), info.get(Category, )], ) \n)跑完之后你就得到一份能直接看懂的注释表哪条蛋白是毒力因子、属于哪个功能类别、相似程度多高一清二楚。文章里要写“我们采用DIAMOND将蛋白序列比对到VFDB数据库以identity≥90%和coverage≥90%作为筛选阈值鉴定出XX个毒力因子相关基因”这段话所需的所有信息就都齐了。4.4 多因素统计与可视化一个小案例注释结果拿到手最常见的下游分析是统计不同样本中毒力因子的种类和数量画一个热图或者柱状图。比如我处理过一个项目比较了12株大肠杆菌分离株的毒力基因携带情况。流程就是每个样本单独跑一遍DIAMOND注释生成一个“样本×毒力因子”的0/1矩阵1表示该样本携带这个毒力因子再用R做聚类热图。矩阵生成这一步用python很容易samples [sample1, sample2, sample3] matrix {} for s in samples: vfs set() with open(f{s}_vfdb_besthit.tsv) as fh: for line in fh: vfs.add(line.strip().split(\t)[1]) matrix[s] vfs all_vfs sorted(set().union(*matrix.values())) with open(vf_presence_matrix.tsv, w) as out: out.write(Sample\t \t.join(all_vfs) \n) for s in samples: row [s] [1 if vf in matrix[s] else 0 for vf in all_vfs] out.write(\t.join(row) \n)然后在R里读入矩阵画热图代码非常简单library(pheatmap) mat - read.delim(vf_presence_matrix.tsv, row.names 1, check.names FALSE) pheatmap(as.matrix(mat), color c(white, steelblue), legend TRUE, main Virulence factor presence)类似的分析思路也可以扩展到毒力因子功能分类的统计比如统计每个样本里Ⅲ型分泌系统、铁载体合成、荚膜合成等各类别因子的数量画堆叠柱状图。这个做出来作为文章主要图表非常有说服力审稿人也会觉得你的毒力注释做得比较系统。5. 常见问题与排查技巧实录5.1 DIAMOND命令报错的几个高频原因先列一下我这两年帮人排查时碰到最多的几个报错。“Error: Too many sequences. Maximum is 10,000,000”这种情况通常是输入文件里有异常空行或者文件格式不对导致DIAMOND把每一行都当作一条序列来读检查一下FASTA文件是否规范首行是“”序列行没有包含空格或非法字符。“Error: Failed to open database file”一般就是数据库路径写错了检查--db参数是否指向了之前建库生成的.dmnd文件注意不要写成FASTA文件路径。“Error: Database version mismatch”意思是DIAMOND软件版本和数据库构建时的版本不匹配。这个问题最容易在上游工具升级后出现解决办法就是重新用当前版本的DIAMOND执行一次makedb。内存相关报错比如“Memory allocation failed”或者“Killed”一般是--block-size设得太大超过了服务器可用内存调小到4或者2就行。另一个容易被忽视的坑是--tmpdir所在分区空间不足。比对时DIAMOND会写很多临时文件如果/tmp空间只有几个GB而你的查询序列又很大任务会莫名中断建议单独设一个空间充足的目录比如/work/tmp。5.2 比对结果为空或者命中率偏低出现“结果文件是空的”或者“几百条序列只有一两条命中”这种情况先不要怀疑数据库按顺序排查第一输入序列是不是蛋白序列。很多人会把核苷酸序列直接丢进来DIAMOND blastp对核酸序列比对也能跑出一堆类似于随机匹配的结果但E值和identity都会很糟糕命中数会少得离谱。确认方法就是看FASTA里的序列字母如果只有A、C、G、T那一定是核酸需要先做基因预测翻译成蛋白。第二阈值是不是卡太严。如果你的pident设成99、qcov设成95然后比对的是跟VFDB数据库物种较远的环境菌株很可能什么都筛不出来。建议先用宽松参数跑一遍比如--id 50 --query-cover 50 --very-sensitive看看原始命中情况再逐步收紧这样能定位是阈值问题还是数据问题。第三数据库选型是不是合适。VFDB核心集setB只有部分代表性序列你研究的物种如果在数据库里没有收录或者收录很少命中数自然偏低。这种情况建议换全量setA试试或者结合CARD、ResFinder等其他功能注释数据库一起看。比如某些细菌的毒力因子研究文献本来就少VFDB收录有限这种情况下多库联合注释比单靠VFDB靠谱得多。第四看你的序列是不是前瞻性地污染了。宏基因组样本里如果宿主序列占比高预测出的蛋白里有大量宿主蛋白这些蛋白自然不可能比对到细菌毒力因子上。这里有条件的话可以先对contig做分类学注释或者blast到NT库把非目标物种的contig过滤掉再跑基因预测和毒力注释。5.3 毒力因子注释结果的边界与数据库更新这里说一个审稿时经常被问到的问题DIAMOND比对结果能不能直接写成“该菌株携带某毒力因子”我的建议是DIAMOND比对结果本身是“序列相似性证据”不等于“功能验证证据”。在文章里写的时候最好表述为“该菌株基因组中存在编码某毒力因子同源蛋白的基因序列”或者“预测含有与某毒力因子高度相似的蛋白编码序列”。这样做既准确又稳妥。如果后续需要更强的证据可以再看比对区域是否覆盖完整阅读框、基因是否位于完整的基因簇内、有没有启动子区域突变等。数据库更新的问题也要提一下。VFDB数据库基本每年都会更新不同版本的序列数量、描述信息可能有差异。做研究的时候建议固定使用一个版本并且在文章的“数据可用性”部分写明数据库下载日期和版本比如“VFDB数据库下载于2024年6月版本为VFDB 2024”。这既是科研规范也方便别人复现你的结果。我自己的习惯是每次下载数据库后把下载日期和文件MD5值记录在一个README里放在流程所在目录长期项目里这个习惯能帮你省很多追溯的麻烦。6. 小结一下我的实操心得最后聊一点个人感受。DIAMONDVFDB这套流程我大概是从2022年开始在项目里大规模使用的到现在跑过的样本加起来有上千个。最大的体会是注释本身不难难点全在“过程和结果能不能经得起追溯和复现”。我见过很多同事用在线工具点几下拿到一份结果就开始写文章最后审稿人要求提供比对参数和数据库版本当场拿不出来只能重跑一遍既浪费时间还影响信心。所以我特别建议你把每一步命令、每个数据库的下载日期和版本、每次阈值调整的原因都记录清楚这不只是给自己留条后路更是科研工作者的基本素养。再分享一个小技巧如果你的样本特别多比如一次跑一两百个基因组建议写一个循环脚本把每株菌的蛋白序列文件按统一规则命名然后逐个调用DIAMOND比对生成的结果文件统一放到results目录下。这样后面统计矩阵、画热图、做功能分类分析都会轻松很多。这算是我踩过不少坑之后总结出的最实用经验了。