第4章：MSA（多序列比对）特征构建原理

咱们直接进入正题。MSA（多序列比对）在AlphaFold里到底扮演什么角色？我个人的理解是——它就像是给模型看的“进化纪录片”。所有来自不同物种的同源蛋白质序列，被排列在一起，算出一个profile。这玩意儿记录了每个位置上，氨基酸能怎么变、不能怎么变。

说白了，AlphaFold能预测出蛋白的三维结构，很大程度上就是因为读懂了这张进化史。没有MSA，模型就失去了判断哪些位点保守、哪些位点可以自由晃动的依据。

4.1 特征构建的第一站：序列数据库与搜索工具

先说说入口。MSA不可能凭空生成，你得有东西喂进去。早期我刚开始做结构预测时，手头就一个单序列，跑个PSI-BLAST就以为完成了任务。后来跟AlphaFold接触久了，才意识到事情远没那么简单。

AlphaFold官方流程中，输入序列首先会被扔进两个数据库：一个是Uniclust30（现在更多用BFD 或 MGnify），另一个是UniRef90。用到的搜索工具是jackhmmer 和 hhblits。我建议在实际工作中，至少交替使用这两种工具，因为它们各有侧重。

• jackhmmer：迭代搜索，每一步用上一步的profile去扫数据库。适合发现远亲同源序列。
• hhblits：基于HMM（隐马尔可夫模型）的比对方法，速度更快，对远程同源检测更敏感。

为什么要有两个工具？原因很简单——单一搜索容易有偏差。我遇到过这样的情况：只用jackhmmer，结果里99%都是人类的同源基因，家族多样性很低；换成hhblits后，突然多了很多从细菌到植物的序列，MSA的“进化广度”一下子就出来了。

4.2 比对结果：不是所有序列都值得留下

搜索完了，你手里可能有几万条序列。但别急着全部塞进模型。相信我，全部保留只会让矩阵爆炸，而且很多低质量的匹配会引入噪声。我个人习惯做几轮过滤：

1. 去除冗余：序列相似度高于99%的，只留一条。
2. 覆盖度筛选：只保留比对上至少50%长度的序列。
3. E-value阈值：默认我选1e-3，实在序列数少可以放松到1e-1。

嗯，这里要注意：AlphaFold模型本身对MSA的深度是有阈值的。我记得官方论文里提到，他们训练时，MSA深度通常限制在1,024或512条序列。如果你拿到的MSA达到了上万条，可以随机抽一个子集，或者按聚类权重采样。

# 一条典型的AlphaFold MSA构建命令（简化版）
jackhmmer \
  -N 3 \
  --noali \
  --cpu 64 \
  -A query.a3m \
  query.fasta \
  uniref90.fasta

# 再跑一次hhblits，增加序列多样性
hhblits \
  -i query.a3m \
  -d bfd_metaclust.fasta \
  -oa3m query_merged.a3m \
  -n 2 \
  -cpu 64

4.3 从比对到特征：Raw MSA如何被编码

拿到a3m格式的比对后，AlphaFold不是直接把它当作特征塞进网络的。这里有很多处理步骤。我先说大体的流向：

好，接下来说说具体的编码逻辑。原始MSA中的每一列（对应一个残基位置），其实有N条同源序列在该位置的氨基酸种类。AlphaFold会做以下几件事：

4.4 具体特征维度拆解

你想想看，模型拿到这些数据后，首先看到的是每个位置上20种氨基酸的出现频率。如果一个位置在进化中只出现“亮氨酸”和“异亮氨酸”两种，模型就知道这个位置很可能是疏水核心的一部分。如果某个位置出现20种氨基酸都有，那模型直接判定这个位置是loop，今天绑这里，明天绑那里，不固定。

4.5 如何处理超长序列与超大MSA

蛋白序列超过1,500个残基怎么办？或者你搜出来的MSA有10万条序列？这些情况我在实际项目中经常遇到。AlphaFold官方方案里有一个硬限制：

• MSA的序列数被截断到最大dapth，默认是max_msa_clusters = 256, 加上随机采样达max_extra_msa = 512。
• 如果总序列数超过这个值，会按hhblits或jackhmmer提供的聚类信息，优先选聚类中心；其余序列则有概率被忽略。

# 模拟AlphaFold对超大MSA的截断策略
total_sequences = 100000
max_clusters = 256
max_extra_sequences = 512

# 第一步：选256个聚类中心
cluster_centers = select_cluster_centers(msa, n=max_clusters)
# 第二步：从剩余序列中随机抽512条
extra_sequences = random_sample(remaining_sequences, n=max_extra_sequences)
# 最终输入序列数 = 256 + 512 = 768

4.6 为什么MSA中的“gap”也重要

别忽略那些横杠。MSA里的gap不是空白，它意味着某些同源序列在此处插入了额外的残基（或者目标序列有插入）。模型通过读取gap的位置，能识别哪些区域是“弹性插入区”。

我遇到过这样一个案例：一个多结构域蛋白，在预测两个domain之间的linker时 AlphaFold总是失败。后来检查MSA特征发现，linker区域在其他同源序列中被多个gap替代。原来那些同源序列根本没有这个linker结构。模型学习到“这个区域可以不存在”，于是预测时就倾向于把它当作一个无序结构。嗯，这才对上了。

所以gap信息可以直接影响模型对柔性区域和结构域边界的判断，千万别把它当作无意义特征扔掉。

4.7 MSA特征质量的自检清单

每次运行完MSA构建后，我建议先做个快速自检。下面是我自用的检查表：

1. 有效序列数：过滤后是否大于50条？低于这个数，模型基本等于在猜结构。
2. 序列多样性：序列是否来自不同物种门类？全是一类物种的MSA容易导致过拟合。
3. 覆盖度：目标序列的N端和C端是否都有覆盖？如果两端序列空洞，模型预测两端通常会糊掉。
4. gap比例：整条MSA的gap比例不宜超过30%。gap太多说明同源性很差，可能混入了不完全匹配的序列。

4.8 从raw MSA到one-hot的实操代码

我写一个简短的demo，模拟AlphaFold内部的特征编码逻辑。这当然不是真正的源代码（官方代码要复杂得多），但核心思想就在这了。

import numpy as np

def encode_msa(msa_seqs, max_len=1024, alphabet='ACDEFGHIKLMNPQRSTVWY-X'):
    """
    模拟对MSA进行one-hot编码
    msa_seqs: 列表，每个元素是一条比对后的序列（长度=target_len）
    max_len: 最大序列数，不足padding
    """
    n_seq = min(len(msa_seqs), max_len)
    target_len = len(msa_seqs[0])
    msa_feat = np.zeros((n_seq, target_len, len(alphabet)), dtype=np.int32)

    aa_to_idx = {aa: i for i, aa in enumerate(alphabet)}

    for i in range(n_seq):
        seq = msa_seqs[i][:target_len]
        for j, ch in enumerate(seq):
            if ch in aa_to_idx:
                msa_feat[i, j, aa_to_idx[ch]] = 1
            else:
                # 遇到未知字符，设为X的索引
                msa_feat[i, j, aa_to_idx['X']] = 1

    return msa_feat

# 假设我们有3条比对的序列
raw_msa = [
    'ACDEFG--HIJKLM',
    'A-DEFGSHIJKLM',
    'AC-EFG-HIJKLM',
]
features = encode_msa(raw_msa, max_len=10)
print(features.shape)  # (3, 13, 23)

这段代码其实把MSA里所有信息都压缩进去了：氨基酸类型、gap、以及残基位置。你发现没有，模型其实并不知道每一条序列背后对应哪个物种，也不知道进化树的结构。它看到的，就是一堆纯字符的统计规律。

所以调优MSA特征，说白了就是在调整“模型相信哪些共进化关系”。MSA覆盖率越高、多样性越好、同源性越明确，模型读出来的进化约束就越准，结构预测结果也就越可信。我建议大家在每次训练前，先把这个步骤跑通，再做其他花活。