用AKtoolbox做协同进化分析：从多序列比对到显著位点对挖掘

简介：AK工具箱是一款面向蛋白质多序列比对协同进化分析的 Matlab 开源软件，遵循简化 BSD 许可证分发，适合生物信息学研究者、结构生物学家及相关专业学生使用。它不依赖 Matlab 自带的生物信息工具箱，集中实现了统计耦合分析、直接耦合分析、互信息等七种主流协同进化算法，并提供多种生物信息学矩阵存档，方便用户对比不同方法的分析结果。压缩包内共包含 116 个文件，其中 103 个为 Matlab 源码文件，覆盖序列读取、二级结构处理和算法主程序等核心模块；同时附带 Windows、Linux、macOS 多平台预编译的动态库，以及蛋白质结构文件、二级结构文件、序列文件和说明文档，解压后即可在常用平台直接运行，省去手动编译环境。整个压缩包大小仅为 264KB，非常轻量。目前已有 13897 人学习。借助完整源码、真实示例数据和跨平台编译文件，读者能快速理解协同进化分析的计算流程，并将其灵活嵌入到自己的蛋白质序列研究项目中，有效减少从零搭建工具链的时间成本。 我最早接触到AKtoolbox，是在分析一批细菌效应蛋白序列的时候。当时已经有多序列比对结果，但光看保守位点远远不够——我更想知道，哪些位置之间存在“你变我也变”的协同信号，也就是协同进化分析。这种信号在蛋白质功能位点挖掘、耐药突变分析和结构预测里特别有价值。AKtoolbox就是一个基于Matlab的开源工具箱，专门把多序列比对文件变成可解释的共变矩阵和显著性热图。它最吸引我的地方是：不需要额外的Python环境，也不用折腾R包依赖，只要熟悉Matlab基本操作，就能跑完从序列读取、矩阵计算到置换检验的完整流程。这篇文章我把算法原理、上手指南和踩过的坑一起写出来，给想用它的朋友一个参考。

1. 这个工具箱到底解决什么问题

1.1 从多序列比对到协同进化分数

很多做序列分析的朋友一开始容易混淆“保守性分析”和“协同进化分析”。保守性分析看的是单个位点在进化中是否保持稳定，比如某一位点基本不变，说明它可能是功能关键位点；而协同进化分析看的是两个位点之间是否出现联合变异——一个位点突变了，另一个位点也随之突变。这件事在蛋白质研究中很重要，因为如果两个残基在空间上靠近，它们在没有直接接触的情况下往往通过互相补偿来维持结构稳定性。比如一个位点体积变大，隔壁位点体积变小，两者一起突变才不会破坏整体折叠。

AKtoolbox做的事情，本质上就是把“两两位点之间是否存在非独立变异”这个问题变成一个可计算、可统计检验的数值问题。它的输入是已经比对好的序列，输出是位点对之间的协同进化得分矩阵和对应的P值。拿到这个结果之后，你可以筛选显著位点对，再映射到三维结构上做进一步分析。

1.2 为什么选择Matlab而不是R或Python

在生物信息领域，R和Python确实是主流，但Matlab有自己的独特优势。矩阵运算是Matlab的看家本领，而协同进化分析的核心恰恰是构建二维矩阵、计算相关性、做置换检验，这些操作在Matlab里几乎不需要写复杂循环，几条矩阵命令就能完成。

另外，Matlab内置的统计工具箱提供了大量现成的概率分布函数和假设检验函数，比如卡方分布、置换检验相关的随机抽样函数，省去了很多底层实现。可视化也是我比较看重的，热图、网络图、散点图在Matlab里有成熟的原生接口，配合AKtoolbox自带绘图函数，出一张论文可用级别的图非常快。

我理解有人会质疑：Python也有Biopython和scikit-learn啊。这话没错，但“够用”和“顺手”是两回事。如果你本身就是Matlab用户，并且只做常规规模的序列数据集分析，开着一个Matlab环境把流程跑通，比再搭一套Python环境高效得多。AKtoolbox开源免费，有需要还能直接改内部算法，这种灵活性也是我选择它的原因之一。

2. 协同进化分析的常用算法与Matlab实现思路

2.1 互信息与卡方检验

协同进化分析最朴素也最常用的算法是互信息（Mutual Information, MI）。它的核心思想来自信息论：如果两个位点的变异相互独立，那么它们联合分布应该是边际分布的乘积；如果两者存在关联，联合分布会明显偏离独立假设。

给定两个位点i和j，互信息定义为：

MI(i,j) = Σ P(xi,xj) log( P(xi,xj) / (P(xi)P(xj)) )

这里P(xi)和P(xj)分别是位点i和位点j的氨基酸分布概率，P(xi,xj)是两位点氨基酸组合的联合概率。如果两个位点完全没有关联，MI值近似为0；关联越强，MI值越大。

在Matlab里做这件事有个很方便的函数——histcounts2，它可以直接统计二维分布的频数。我实现的MI计算核心大概长这样：

function mi = calc_mi(col_i, col_j, alpha) % col_i, col_j 是两列氨基酸索引，alpha是氨基酸类别数 N = numel(col_i); % 计算两位点的联合频数 joint = accumarray([col_i, col_j], 1, [alpha alpha]); joint = joint / N; % 计算边际频数 pi = sum(joint, 2); pj = sum(joint, 1); % 避免log(0) joint = joint + eps; pi = pi + eps; pj = pj + eps; % 展开成向量并计算互信息 joint_vec = joint(:); prod_vec = (pi * pj); prod_vec = prod_vec(:); mi = sum(joint_vec .* log(joint_vec ./ prod_vec)); end

我用accumarray替代双重循环，在序列数几千、位点长度几百的情况下，运算效率依然可以接受。单纯MI值只能给变异关联的强弱排序，不能直接判定显著性，所以通常还需要配合卡方检验或者置换检验。

2.2 背景信号去除与置换检验

做协同进化分析有一个不能忽略的干扰因素：系统发育背景。如果一组序列拥有共同的进化历史，那么不同位点会沿着同样的进化树分化，导致位点之间出现大量假阳性的相关性。也就是说，你看到的“协同进化”，可能只是因为所有位点都跟着同一棵树的拓扑结构在变，而不是位点之间真的存在功能耦合。

处理这种背景信号常用的做法是“残差化”：先计算两个位点之间的进化距离或相似度，拟合一个背景模型，再把残差当作真正有意义的协同进化信号。AKtoolbox中相关模块也遵循这个思路——先计算原始MI矩阵，再用回归方法去拟合背景距离，最后输出残差化的得分。

置换检验是用来判断显著性最直观的手段。它的做法是把一个位点的氨基酸顺序打乱，破坏原有的组合关系，然后重新计算互信息。重复很多次之后，你会得到一个“随机情况下互信息可能达到的水平”分布。把真实观察到的MI值和这个分布比较，就能算出P值——如果真实值落在分布的尾端，说明不太可能随机产生。

置换检验的代码模式我写了很多次，核心结构是：

function pval = permutation_pval(obs_mi, col_i, col_j, alpha, nperm) null_mi = zeros(nperm, 1); N = numel(col_i); for r = 1:nperm col_j_perm = col_j(randperm(N)); null_mi(r) = calc_mi(col_i, col_j_perm, alpha); end pval = (sum(null_mi >= obs_mi) + 1) / (nperm + 1); end

注意最后加1做平滑处理，避免出现P=0这种不合理的极端值。置换次数越多，P值越精确，但计算成本也越高，实际使用需要在两者之间平衡。

2.3 Matlab向量化计算的细节

我刚开始写的时候没太注意向量化，用双层循环遍历所有位点对，结果序列长度500的MSA就要算12万个位点对，跑一次要一个通宵。后来优化了两次，速度提升非常明显。

第一个优化是减少重复计算。对于MI这类对称统计量，矩阵是对称的，只需要计算上三角部分，直接省掉一半计算量。第二个优化是把字符比较全部改成整数索引比较。Matlab处理char类型虽然方便，但在大规模矩阵计算上，数值类型的速度明显更好。我的习惯是把氨基酸字母映射成1到20的整数索引，gap处理为0，后面所有计算都用整数矩阵。

第三个优化是使用parfor并行。Matlab的并行计算工具箱在独立位点对的计算上几乎可以线性加速，因为每个位点对的计算互不依赖。我一般会在机器CPU比较充足的时候开8个worker，把置换检验的循环丢进去跑。如果你的机器内存不大，记得把大矩阵用single类型存储，精度下降有限，但内存占用几乎减半。

3. AKtoolbox上手实操

3.1 获取代码与运行环境准备

AKtoolbox既然是开源项目，最直接的获取方式就是从代码托管平台克隆仓库。下载解压之后，建议把整个工具箱目录添加到Matlab路径中，这样在任何目录都能直接调用函数。命令行操作如下：

addpath(genpath('D:/AKtoolbox')); savepath;

我建议用genpath一次性递归添加所有子目录，避免漏掉某些辅助函数。savepath保存路径设置，这样下次启动Matlab就不用重新添加了。

版本方面，我使用Matlab 2020b到2023b都测试过，核心函数没有遇到兼容性问题。如果你的Matlab版本比较老，主要留意histcounts2这个函数是否可用——它是2015b之后引入的，太老的版本需要换成hist2或者手动分箱。

3.2 输入文件的准备与检查

AKtoolbox的输入是多序列比对文件，而不是原始未比对序列。很多第一次用的人在这步踩坑：拿一堆未比对的序列直接丢进去，结果程序报错或者结果完全不可信。

如果你想分析一批同源序列，需要先用MAFFT、Muscle或Clustal Omega做多序列比对，导出为FASTA格式。比对好的FASTA长这样：

>seq1 MSTNPKPQRKTKTV >seq2 MSTNPKPQRKTKSV >seq3 MSTNPKPQRKTKTV

所有序列的长度必须一致——因为比对之后每一列代表同一个进化位置。导入之后建议先检查一遍数据质量。我的判断标准有三个：一看序列中重复序列是不是过多，二看gap占比，三看有没有非标准氨基酸字符。gap占比太高的列对MI计算干扰很大，比如某一位点一半序列都是gap，算出来的联合分布会很稀疏，容易产生虚假的高MI值。遇到这种情况，我一般会直接过滤掉gap比例超过20%的列，再做分析。

3.3 核心调用流程和参数选择

AKtoolbox的接口逻辑比较清晰，核心步骤大致是：读入序列、转成数值矩阵、计算协同进化得分矩阵、做置换检验、画图。下面这段代码是我根据实际仓库的demo改编的典型调用过程：

% 读取比对好的FASTA文件 seqs = fastaread('my_alignment.fasta'); % 把结构体数组转成char矩阵，这一步要求所有序列等长 aln = char(seqs.Sequence); % 将字母映射为数值索引，自定义函数做清洗和过滤 alnNum = AK_prepare_alignment(aln, 'remove_gap_cols', true); % 计算协同进化得分矩阵 miMat = AK_mutual_information(alnNum); % 置换检验，500次置换，得到P值矩阵 pMat = AK_permutation_test(alnNum, 'nperm', 500); % 画出热图 AK_plot_coev(miMat, pMat);

关于置换次数，我给个实用建议：初次探索数据集的时候，500次就够用了，几分钟能跑完；如果你打算用结果支撑实验验证或者投稿，建议至少跑1000次。置换次数太少，P值分辨率太低，很多临界显著的位点对会被判定为不显著。

关于参数选择，还有一个很容易被忽略的问题：序列数量。协同进化分析本质上需要足够的样本量来估计联合分布。如果你的序列数量少于50条，MI估计会非常不稳定，此时即使跑出来结果，也建议只当线索，不要当结论。

3.4 结果解读与可视化

AKtoolbox输出的核心结果是一个位点对协同进化得分矩阵，行列都是位点编号，数值越高表示协同进化信号越强。P值矩阵则用来筛掉不可信的得分。

拿到结果后，我建议按以下顺序处理：先设置一个P值阈值（比如0.05），筛选显著的位点对；再按得分从高到低排序，优先关注Top 20的位点对；最后把这些显著位点对映射到蛋白结构上，看它们是否在空间上靠近。如果两个位点在一级序列上距离很远，但三维结构上距离很近，这种协同进化信号就值得重点研究——它很可能指向一个功能相关的远程相互作用网络。

可视化这块，AKtoolbox自带的热图功能相当实用。你可以在图上叠加显著性标记，比如用星号标出P值小于0.01的位点对。不过要注意：热图适合展示全局模式，真的要看具体位点对，最好以数值表格为准。我通常会把显著的位点对导出成CSV，方便后续用PyMOL做结构映射。导出命令很简单：

[i, j, pval] = find(pMat < 0.05); T = table(i, j, miMat(sub2ind(size(miMat), i, j)), pval, ... 'VariableNames', {'Pos_i', 'Pos_j', 'MI', 'P_value'}); writetable(T, 'significant_pairs.csv');

4. 实战踩坑记录与排查方法

4.1 序列矩阵化容易出错的几个点

第一个坑是字符编码不一致。多个测序文件拼接后，有的序列用大写字母，有的用小写字母，导致在映射索引的时候匹配不上，结果全是错误码。我的习惯是读入之后立刻统一字符：

aln = upper(aln);

第二个坑是非标准氨基酸字符。B（天冬酰胺或天冬氨酸）、Z（谷氨酰胺或谷氨酸）、J（亮氨酸或异亮氨酸）、X（未知氨基酸）在真实数据里并不少见。如果不处理，它们会被映射成错误索引，干扰后续计算。我一般把这些字符统一视为gap或者直接丢弃该列，具体看它们在整个序列中出现的比例。

第三个坑是gap的编码方式。不同比对工具导出的gap可能用-、.或~表示，读进来之后一定要统一替换成同一个符号再处理。如果你在计算MI时把gap当作一种普通字符参与统计，结果可能完全跑偏，因为gap的含义是“缺失位置”，和真实的氨基酸状态本质不同。

4.2 置换检验太慢怎么办

这是我自己遇到过最头疼的问题。序列数量500，位点长度600，置换500次，相当于要计算500×(600×599/2) ≈ 9000万次位点对MI。如果不做优化，在普通电脑上跑几天都算不完。

我后来总结出三个有效的提速方案。第一，只对上三角做置换检验，因为MI矩阵是对称的；第二，把置换循环改成parfor并行；第三，如果真的算力紧张，不要对全部位点对做检验，先按MI得分从高到低排序，只对得分最高的前5%位点对做置换检验——这样不仅能大幅缩短时间，还能让检验更聚焦在潜在信号上。这个策略在文献里也有应用，逻辑上说得通：置换检验本质是筛显著性，没信号的位置做再多置换是浪费时间。

4.3 不同工具的效果对比

我除了AKtoolbox，也用过R里的Bio3D和Python生态中的EVcouplings。简单做个对比，方便你按实际需要选型：

如果你手头只有几十条同源序列，其实不太建议直接上Potts模型类工具——参数多、数据需求大，容易过拟合。这种情况下，用AKtoolbox这类轻量互信息方法更稳妥。反过来，如果序列量达到几千甚至几万条，Potts模型的精度优势就体现出来了，MI类方法可能会低估复杂耦合关系。

我在实际操作中更愿意把AKtoolbox当作快速筛查环节：先跑一遍互信息矩阵和置换检验，筛出强共变位点对，再根据序列量决定是否需要进一步上更复杂的模型。输入比对清理干净，置换检验次数给足，剩下的交给工具去跑，大部分情况都能得到可靠线索。最后还是要提醒一句：协同进化信号是很好的提示，但它不等同于物理上的直接相互作用，把它当线索而不是结论，后续一定要结合结构信息和实验证据验证。

本文还有配套的精品资源，点击获取