MMseqs2 搜到一个 100% 一致的蛋白，为什么还要看覆盖度？

假设你从宏基因组里预测出一条 400 个氨基酸的蛋白，想在参考库里找线索。搜索得到两个候选：一个有 90 个氨基酸与查询序列完全相同；另一个的一致率只有 77.6%，却几乎从头对到尾。该先看哪个？如果要把整个蛋白库去冗余，这三条序列又该不该归到一组？

这正好是 MMseqs2 的两种用途。搜索让你拿着一条或一批查询序列，找到参考库里的匹配对象，并检查每一对比对证据。聚类从一批待整理的序列出发，根据设定的相似程度和覆盖范围，把它们归组并选出代表序列。下面用同一组人工序列跑一遍，两种结果就容易读懂了。MMseqs2 官方仓库与本文使用的版本 18-8cc5c

大库里怎样找到值得细比的序列？

设想把查询蛋白和库里每一条蛋白都从头做一遍带缺口的局部比对。即使两条序列几乎毫无关联，也得付出一次比对的计算代价。MMseqs2 常规 CPU 蛋白搜索的做法是先用成本较低的线索缩小候选范围，再把精细比对留给有希望的序列对。MMseqs2 方法论文

第一条线索是 短词，即由若干氨基酸组成的 k-mer。程序利用目标库的索引，查找查询序列中短词及其相似短词出现的位置。允许相似氨基酸替换很重要：哪怕两个蛋白保留了相近的结构，其中一些残基也可能已经改变。

第二条线索是位置是否对得上。假设查询的第 10 位和第 30 位各有一个短词，分别命中同一目标的第 50 位和第 70 位。两次命中的位置偏移都是 40，在比对矩阵里就落在同一条对角线上。相比孤立命中，这更值得沿着该方向做一次无缺口评分。通过预筛的候选再进入允许插入和缺失的局部比对，得到比对区间、得分和 E 值。预筛是启发式的；很远的同源序列可能进不了候选名单。若本来就在找远缘成员，可以提高搜索灵敏度参数 -s，再观察新增候选及其比对证据。官方手册中的预筛与搜索原理

读一张真实输出表：一致率的分母是谁？

为让数字有出处，我们构造了一条 400 个氨基酸的查询序列 Q，以及三个目标。长目标 full_length 由 Q 替换部分残基得到；短目标 domain_only 直接截取 Q 的一段连续序列；第三条是独立生成的随机序列。短目标的名字用来模拟“只匹配一个局部区域”的场景，这些人工序列没有实验功能注释。下面是这次小例子实际使用的搜索命令，输入为蛋白 FASTA。

mmseqs easy-search query.fa reference.fa hits.tsv tmp_search \
  --threads 2 --alignment-mode 3 \
  --format-output "query,target,pident,alnlen,qlen,tlen,qcov,tcov,evalue,bits"

--alignment-mode 3 让程序直接按比对计算一致率。pident 的分母是局部比对的列数，含内部缺口列。 第二行表示这段比对中的 90 个残基全部相同；用 90 除以查询全长 400，得到的查询覆盖度只有 0.225。目标恰好只有 90 个残基，所以同一行的目标覆盖度是 1.000。覆盖度在输出中用 0 到 1 的小数表示。

E 值回答另一个问题：在当前搜索空间的随机背景下，预期会出现多少个得分至少这么高的偶然匹配？它与数据库规模有关。比特得分越高，表示比对得分越强。在本例中，两条命中都有很小的 E 值，但这并没有把 90 个残基的局部证据延伸到其余 310 个位置。官方输出字段与一致率定义

如果目标是找整条蛋白的功能候选，会优先检查长目标的注释来源、结构域排列和关键残基。如果目标是定位一个保守结构域，短命中同样值得保留，还应导出 qstart,qend,tstart,tend 查看具体位置。面对天然蛋白，短命中也可能来自截断的基因预测。仅凭首位命中，无法确定两个蛋白具有相同的底物特异性。

用同一批序列做聚类，为什么短目标会单独成簇？

现在把 Q 与三个目标合在一个 FASTA 里。假设目的变成整理一份全长蛋白代表集，我们希望长度相近、主体序列相似的蛋白进入同一组，同时保留短片段。下面是本例实际运行的建库、聚类和成员表导出命令。

mmseqs createdb proteins.fa proteinsDB
mmseqs cluster proteinsDB clustersDB tmp_members --threads 2 \
  --min-seq-id 0.7 -c 0.8 --cov-mode 0 --alignment-mode 3
mmseqs createtsv proteinsDB proteinsDB clustersDB members.tsv --threads 1

这里选择 70% 序列一致率和 80% 覆盖度来演示全长去冗余，两项阈值都应随研究目的调整。--cov-mode 0 要求比对同时充分覆盖两端。对于 400 对 90 的例子，即使短序列整段都对上，最多也只能覆盖长序列的 22.5%，达不到这一门槛。官方覆盖模式与聚类说明

原始 members.tsv 每行两列，第一列为代表 ID，第二列为成员 ID。要统计簇大小，按第一列分组计数即可。代表是聚类算法选出来的序列；如果要用它做功能注释，仍要核对注释质量。尤其在真实大库里，同簇也不能当作“任意两成员都已逐对验证超过阈值”的证明。想比较簇内两条具体蛋白，就把它们单独对齐，查看覆盖范围与差异位置。

>Q
GMNLNLMFDGSAIKNKVWNFVVLSVYFFTGGASFEPCSCYDDLSLRFVFKTICRCKHLDWTMHIGKLRCWVMMCFADREQSKQNPNPEYFDFTYYDISDKNGGLDVLELIPCICNLSPMPHMILKAFNDITLSHAKSLIVIWREPPNQIRNSSPMFCLLILVVVDFVTFSAIPIVFEKGRQDECYFEGKERMHYMRHKWEYGAITPLAMYTTCEHSNHCHPPWNQAHDGSMRVIKANPYHGYAGPYFKSDAGRMVGSPYSPFKCDPYNGRTYWDIFWAFMHVDMVIITWMRNHIFRVSPHYAQEMTVNYLMCDLSSIWCMGEIRDPQLNGIDMHRRPDQMMKPDWERCMPPQHMWKSCHCDVMLQNQTPRFRKWTGHDGCQTRLTIGAVTYNAWQWAYDG
>domain_only
AFNDITLSHAKSLIVIWREPPNQIRNSSPMFCLLILVVVDFVTFSAIPIVFEKGRQDECYFEGKERMHYMRHKWEYGAITPLAMYTTCEH
>full_length
KMNPDLMFDWSQIWTKVWNFVVLSVYEFTGGASFEPCSPFGDLALGFVFKYIHRMKHLDWTDYIGKLRCWVMMCFADDEQYKENPNPELFCVEYYDISDQNISWDVLELTPFICNVNPMPTMMLKNFNDITLSHAPVLIVIWCEPQNHIRPSSPMFMLLPLVVVDFVNFSAYPIVFEKGAQDECYFEKVERMHYMRHKPEYIAITPLGMTFPCEHSNHCHSNVAQAHDGSMRVVKANPYPGYFGPYVKSDASNVNGSPYSPRKCDPYMGRTYWDIFWCFNEVDMVNITWMCNHEFRVSPHYSQEMTVNYLMCDLSSIWCMGEIRDPQLNGEDMWRRPDQMMKPDWGCCMPPQHMQKSCHTDCMDQNQGPRFRKWTGHDKCQNRLTIGAYTYIAWQWKYDG
>unrelated
IPEPANLWPSLSHQDDWKTVNHHAIPWCLRSNKAWCQPRHIARFWECATTMDKDTWIPFTMGKLGCVDGGLHNWEGTDIAKGLREIHSGGYVAKYMIQCHRMDKWMMGTGFIAYGSPTDLWIHQIGRMVFVCVRFTGWNHCSWQGLSVSFPIKNNMECINIQFSAIKEMCPHQAKKEFNPKQMGLPYVMSTWEPWCLTLEPKEVCLTRVWICNCYNESSCDTPPKDCCDMTYWFIWEATESKKIVFCPPMWSWKCNLTNFEWVPASMYFRQHYFSKGKEFRSETRVCKRQSPWNCYLNASSISERVNHLQIEPNYHYNLLEIKWAHHWYPFAMRDCLKHEKNDYRFFYEICITANATWTRQRPETNRTWTGYVRYMNHGVGTSMYSCDDLFIWATGLNLG

下次拿到一个命中，先问自己要保留什么

给一条未知蛋白找线索，保存命中序列、双向覆盖和比对坐标，之后才有依据区分全长候选、共享结构域与片段。给整库去冗余，保存输入库版本、聚类参数和代表—成员映射，否则只留下代表 FASTA，以后很难追溯某条样本序列去了哪里。本例中的两张表之所以给出不同答案，正是因为“能找到局部相似”和“适合用一条全长代表替代”是两个不同的问题。

实际蛋白库若超出本地处理能力，可根据目标库大小、索引占用和试跑时临时目录的峰值，按需查看云生信科研服务器租赁配置。