企业官网建设流程全解析

做致病菌基因组分析的朋友应该都有过这种体验：手上拿到一株菌的组装结果，老板第一句话就问“它有哪些毒力因子？有没有可能致病？”这时候你急需一套快速、可复现、结果还能写进文章的注释流程。毒力因子注释最常用的办法，就是把预测出来的蛋白序列跟专门的毒力因子数据库做比对。而在比对工具里，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秒不到，结果几乎能对上，这个效率差距在做批量样本时就是天壤之别。速度对比可以看下面这个表，数字来自我本地服务器上的实测，不同机器会略有差异，但趋势基本一致。

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 python=3.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列的含义，我做了一张速查表，建议你收藏一下：

我想特别提醒的是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里很多毒力因子序列并不是全长的，数据库里的条目可能只覆盖基因的一部分，强行要求数据库序列全覆盖会丢掉有效结果。

筛选阈值这里多说一句，如果分析的是高度保守的毒力因子标志基因，比如肠杆菌科的侵袭蛋白A（InvA），阈值设到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, key=lambda x: x[0], reverse=True)[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读取，然后导出成TSV：

import pandas as pd df = pd.read_excel("VFs.xls", header=1) df.to_csv("vfdb_annot.tsv", sep="\t", index=False)

然后写一个小脚本，把比对结果跟注释表合并：

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. 小结一下我的实操心得

最后聊一点个人感受。DIAMOND+VFDB这套流程，我大概是从2022年开始在项目里大规模使用的，到现在跑过的样本加起来有上千个。最大的体会是，注释本身不难，难点全在“过程和结果能不能经得起追溯和复现”。我见过很多同事用在线工具点几下，拿到一份结果就开始写文章，最后审稿人要求提供比对参数和数据库版本，当场拿不出来，只能重跑一遍，既浪费时间还影响信心。所以我特别建议你把每一步命令、每个数据库的下载日期和版本、每次阈值调整的原因都记录清楚，这不只是给自己留条后路，更是科研工作者的基本素养。

再分享一个小技巧：如果你的样本特别多，比如一次跑一两百个基因组，建议写一个循环脚本，把每株菌的蛋白序列文件按统一规则命名，然后逐个调用DIAMOND比对，生成的结果文件统一放到results目录下。这样后面统计矩阵、画热图、做功能分类分析，都会轻松很多。这算是我踩过不少坑之后总结出的最实用经验了。