📚 课程目录第4章多序列比对原理

多序列比对原理：MSA的构建方法、进化信息与保守性分析

聊多序列比对（MSA），我脑子里第一个蹦出来的词是“对齐”。说白了，就是把一大家子序列按同源位置排在一起，看哪些位点保守、哪些位点突变。我刚开始做蛋白质结构预测时，总觉得MSA就是个自动步骤，跑个MAFFT就完事了。后来发现，MSA质量直接决定进化信息——而进化信息又是AlphaFold的“燃料”。这里头坑不少，咱们一个一个说透。

MSA的构建方法

构建MSA的方法可不止一种。我习惯把主流方法分为三类：渐进法、迭代法、一致性法。

方法类别	常见工具	核心思路	我遇到过的场景
渐进比对	ClustalW, Clustal Omega	构建导引树，从最相似序列两两比对起，逐步添加	序列数少（<200）时非常快；但一旦序列多，早期错误会累积
迭代比对	MAFFT, MUSCLE	反复优化，修正渐进比对的局部错误	处理中等规模数据集（几百到几千条）时，我常用MAFFT --auto
一致性比对	T-Coffee, ProbCons	综合多个比对结果或结构信息，生成更可靠的比对	对于高变区或者远源同源，一致性方法能减少“错位”——但跑得慢

举个例子，我现在处理一批细菌组蛋白序列，经常先跑MAFFT，再加个--maxiterate 1000。为什么？因为迭代能修正棘轮效应（ratchet effect）带来的积累偏差。我曾在一个项目中只用了渐进比对，结果一个保守基序被严重错位，导致模型预测完全跑偏——后来换成迭代比对，那块儿就对齐了。

核心公式（打分矩阵）
比对的核心是打分。最常用的是BLOSUM62，它的对数比值：

Score(a,b) = log2( P(a,b) / ( P(a)*P(b) ) )
其中P(a,b)是观测到a与b比对的概率
P(a)是氨基酸a的背景频率

在实际MSA构建时，每个位点的gap开销（gap open & extend）需要在实践中调整。我建议先试试默认参数，再用手动检查几个保守位点是否对得齐。

进化信息：从比对中挖出“隐藏信号”

MSA建好了，怎么提取进化信息？AlphaFold的做法很直接——对每一列（位点）算一个进化特征。我把它拆成两步：序列权重与位置特异性得分（PSSM）。

序列权重

为什么要加权？因为数据库里有些序列高度相似（比如同一物种不同菌株），它们会过度代表。如果不砍权重，保守性分析就会偏。常用的方法叫HHblits权重：

for each sequence i:
    effective_seq_weight[i] = 1 / (1 + sum_{j!=i} (seq_identity(i,j) > threshold ? 1 : 0))

我印象很深——有一次没做序列去冗余，结果MSA里80%都是同一个物种的拷贝，导致保守性分析说“异亮氨酸非常保守”，实际上那是采样偏差。后来我用HHblits自带的权重，这个问题就缓解了。

位置特异性得分矩阵（PSSM）

PSSM就是每列一个20维向量，代表该位点倾向于出现哪类氨基酸。AlphaFold的MSA特征里，每个位点的进化信息就是PSSM+profile+位置熵的混合体。构建PSSM的传统方法：

PSSM[i][aa] = log2( f[i][aa] / background[aa] )
其中f[i][aa]是第i列出现氨基酸aa的加权频率
background[aa]是背景频率（如来自SwissProt）

注意，当f[i][aa]很小时（比如0.01），PSSM会变成负分，表示该氨基酸在该位点不太可能出现。这个信息对AlphaFold的Attention机制特别有用——它知道“这里不能放Glu”。

保守性分析：找到关键残基

保守性分析说白了就是看MSA中每一列有多“稳定”。我常用两种方式：序列Logo和位置熵。

序列Logo

Logo把每个位点的氨基酸按保守程度堆叠起来，越高的字母越保守。构建Logo时每个字母的高度由两部分决定：位置特异性贡献与信息量。公式不复杂：

R(i) = log2(20) - H(i)
其中H(i) = - sum_{aa} p(i,aa) * log2(p(i,aa))
p(i,aa)是位点i出现aa的频率

这个R(i)越高，说明该位点越保守。我记得做丝氨酸蛋白酶催化三联体时，位点His、Asp、Ser三个字母在Logo里几乎占满整列——绝对保守。这给我省了很多突变实验的时间。

位置熵H(i)

位置熵和R(i)是一体两面：熵越低，保守性越高。AlphaFold在输入特征中直接用了H(i)作为权重的补充。我建议在检查MSA质量时，可以输出每个位点的熵：

import numpy as np
# 假设msa_matrix是形状 (N_seq, L) 的整数编码矩阵
def positional_entropy(msa):
    N, L = msa.shape
    entropy = []
    for i in range(L):
        counts = np.bincount(msa[:,i], minlength=21)  # 20氨基酸+gap
        freqs = counts / N
        freqs = freqs[freqs > 0]
        ent = -np.sum(freqs * np.log2(freqs))
        entropy.append(ent)
    return entropy

这个函数可以快速定位哪些区域“乱哄哄”——那些熵高的区域通常对应loop或无序区，不适合做模板。

一个重要警告
不要只跑一个MSA工具就交差。我曾在不同数据集上试过Clustal Omega和MAFFT，结果保守位点有5-10%的差异。建议至少跑2个工具（如MAFFT + HHblits），然后取一致性高的区域做后续分析。
另外，MSA的gap处理很关键。有些研究表明gap插得太少（过对齐）或太多（欠对齐）都会让AlphaFold的置信度下降——我习惯用MAFFT的--ep 0.123参数控制gap惩罚。

图：MSA构建与保守性分析的典型流程。注意每一步都可能引入偏差，需要人工检查。

实践避坑指南

我踩过的坑，列出来给大家参考：

选工具要看序列规模：100条以下用Clustal Omega很快；几千条用MAFFT iterative；上万条必须用HMM-based（如HHblits或HMMER）。我见过有人用ClustalW跑5万条序列，跑了三天还没完——太不值了。
gap惩罚别图省事：默认参数往往偏向于“尽量少插gap”。但对高变区，适当容忍gap能提升比对的生物学意义。我处理RNA结合蛋白时，把gap open从默认的-2.0改成-1.5，结果更合理了。
别迷信一个MSA：AlphaFold论文里说他们用HHblits + jackhmmer构建两个MSA，然后合并。你在自己的任务中也可以这么做——多源比对能补假阴性。我试过HHblits单独跑，有些远源同源会漏；加上jackhmmer之后，多找到几个关键家族。
检查序列同源性阈值：很多工具默认序列相似性阈值是0.62，但这个对全局比对不够严格。我建议在做领域特异性数据集时，手动设--id 0.8去除冗余。

一个小技巧
当你手头有结构数据（比如已知的晶体结构），可以用结构比对辅助MSA。工具如Expresso（T-Coffee的变种）会把结构信息作为约束，效果非常稳。我一次做跨膜蛋白比对时，纯序列比对全乱了，加结构约束后一下就对了。

总结：MSA是AlphaFold的“命根子”

我个人经验，AlphaFold预测的置信度（pLDDT）在保守区域通常很高，在进化信息少的区域就会掉。所以，花时间把MSA建好，比后期调参更值。

嗯，这儿注意一点：不同数据库（UniRef90 vs BFD）得到的MSA覆盖度不一样。我建议先用UniRef90跑一遍，再用BFD或MGnify补充。如果你做的是细菌或病毒蛋白，UniRef可能不够，BFD通过宏基因组数据能捡到很多远源同源。有一次我处理一个未知功能蛋白，UniRef90只比对出30条同源，BFD一下子给到700条，模型直接飞升。

保守性分析出来后，你可以用这些信息指导突变试验分析。比如某个位点保守到全是甘氨酸，那你换成脯氨酸试试——我敢打赌结构肯定会崩。这也是AlphaFold帮助设计突变实验的经典玩法，后面章节会深入讲（但今天先到这儿）。