4. 多序列比对与同源搜索实战

在蛋白结构预测流程中，同源搜索是决定后续建模质量的关键前置步骤。无论是为 AlphaFold2 准备 MSA 输入，还是用模板法寻找结构模板，本质上都依赖同一件事：从海量序列库中快速、准确地找出与查询序列同源的候选。本节围绕三个主力工具——BLAST、MMseqs2、HMMER——讲清它们的定位差异、安装方式、关键参数含义，以及如何组织一次可复现的大规模搜索任务。

三个工具的定位与选择

三者并非互相替代，而是覆盖不同的搜索场景。BLAST 是最通用的启发式比对工具，blastp 适合单条或少量查询对中等规模库的搜索，结果可解释性强；MMseqs2 以极高的速度见长，通过预过滤和 k-mer 匹配大幅压缩搜索空间，适合数万到数百万条序列的批量搜索，是构建大规模 MSA 的常用选择；HMMER 则基于隐马尔可夫模型，用 hmmsearch 或 jackhmmer 做谱搜索，对远缘同源更敏感，适合寻找低相似度但结构相关的家族成员。

一个实用的判断原则是：需要速度与规模，选 MMseqs2；需要单条查询的详细比对与 E 值解释，选 BLAST；需要检测远缘同源或做迭代搜索，选 HMMER。实际工程中常组合使用，例如先用 MMseqs2 快速筛出候选，再用 HMMER 对候选做精细验证。

安装与环境准备

三个工具都可通过 conda 安装，建议为每个项目建立独立环境，避免依赖冲突。以下命令为示例，实际版本号以当前发布为准：

# 创建独立环境
conda create -n homolog_search -c bioconda -c conda-forge \
  blast mmseqs2 hmmer -y
conda activate homolog_search

# 验证安装
blastp -version
mmseqs version
hmmsearch -h | head -n 3

数据库是搜索的另一半。BLAST 需要 makeblastdb 建库，MMseqs2 用 mmseqs createdb，HMMER 则直接读取 FASTA 或 HMM 文件。以 UniRef50 为例，下载后建库：

# BLAST 建库
makeblastdb -in uniref50.fasta -dbtype prot -out uniref50_db

# MMseqs2 建库
mmseqs createdb uniref50.fasta uniref50_db

# HMMER 可直接用 FASTA，或先建 HMM
hmmbuild family.hmm family_alignment.sto

关键参数调优

参数决定搜索的灵敏度与速度平衡。BLAST 中，-evalue 控制显著性阈值，默认 10 对蛋白搜索偏宽松，构建 MSA 时常收紧到 1e-3 或 1e-5；-word_size 默认 3，调小到 2 可提高灵敏度但变慢；-max_target_seqs 限制返回条数，注意它并非“最优 N 条”，而是“最先满足阈值的 N 条”，大规模搜索时建议设大一些再后处理。

MMseqs2 的核心参数是 -s（灵敏度，1–7.5，越高越敏感越慢）和 --max-seqs。构建 MSA 常用 -s 5.7 或 -s 7.5，配合 --max-seqs 1000。其 easy-search 封装了建库、搜索、格式转换全流程，适合快速上手：

mmseqs easy-search query.fasta uniref50_db result.m8 tmp \
  -s 5.7 --max-seqs 1000 --format-output "query,target,evalue,pident,alnlen"

HMMER 的 hmmsearch 中，-E 是序列层面的 E 值阈值，--domE 是结构域层面阈值，--incE 和 --incdomE 控制哪些命中被纳入结果。做远缘搜索时，jackhmmer 的迭代次数 -N 很关键，通常 3–5 轮即可，过多会引入噪声。

大规模搜索任务的组织

当查询序列达到数千条时，逐条调用 BLAST 会非常低效。推荐用 MMseqs2 的批量模式，或把查询拆分成多个 FASTA 分片并行提交。以下是一个可复现的批量搜索流程示例：

# 1. 建库（若尚未建）
mmseqs createdb uniref50.fasta uniref50_db

# 2. 批量搜索，输出为 m8 格式
mmseqs search queries_db uniref50_db result_db tmp \
  -s 5.7 --max-seqs 1000 --threads 16

# 3. 转换为可读结果
mmseqs convertalis queries_db uniref50_db result_db result.m8 \
  --format-output "query,target,evalue,pident,alnlen,qcov"

结果判读时，重点关注 evalue、pident（一致性百分比）、alnlen（比对长度）和 qcov（查询覆盖度）。构建 MSA 时，通常保留 E 值小于 1e-3、覆盖度大于 50% 的命中，再送入 MAFFT 或 Clustal Omega 做多序列比对。若命中数过少，可放宽 E 值或提高 MMseqs2 灵敏度；若命中数过多且噪声大，则收紧阈值并检查是否有重复序列。

搜索策略与结果流向

下图概括了从查询序列到 MSA 输入的完整数据流，以及三个工具在其中的分工与关键参数位置。图中标注了各环节的输入输出格式和判断条件，便于对照实际操作。

整个流程的核心判断点在于过滤环节：E 值和覆盖度阈值直接决定 MSA 的深度与质量。若后续结构预测效果不佳，可回到搜索环节调整灵敏度或更换数据库，而不是盲目增加迭代次数。把搜索参数、数据库版本和过滤条件记录在项目配置中，是保证结果可复现的基本要求。