简介AKtoolbox是一款面向生物信息学研究者的开源Matlab工具箱基于简化的BSD许可证分发用于蛋白质多序列比对MSA的协同进化分析且不依赖Matlab生物信息工具箱。适合从事蛋白质结构预测、功能位点识别和分子进化研究的科研人员使用尤其适合需要在自有数据上快速验证协同进化信号的中高级用户。包内已实现统计耦合分析SCA、互信息MI、直接耦合分析DCA、ELSC、OMES、MCBASC和logR等多种主流协同进化计算方法。整个资源共116个文件以103个m格式函数源码为主体另含面向Windows、Linux与Mac平台的预编译mex可执行文件、fasta/dssp/pdb示例数据、txt说明文档及生物信息学矩阵存档压缩包仅264KB轻量易部署。目前已有13958人学习下载可直接调用上述算法或在其基础上扩展是开展MSA协同进化研究的实用工具箱。1. AKtoolbox 是什么用多序列比对做协同进化分析的开源 Matlab 工具箱AKtoolbox 是一套基于 Matlab 的开源工具箱目标很直接给你一条多序列比对multiple sequence alignmentMSA帮你找出哪些位点在进化上是“绑在一起”的。这类协同变化往往对应功能偶联、结构接触甚至变构调控的信号也是预测蛋白质功能位点的重要线索。相比大家熟悉的单列保守性打分协同进化分析看的是“列与列之间的关系”所以它对输入比对的质量、序列数量、清洗方式都更敏感。如果你手头已经有同源序列的比对文件想做位点共变分析但又不想从零写统计模型AKtoolbox 这类工具箱能把互信息、SCA、背景校正、可视化串成一条完整流程。适合的人群很明确做蛋白质功能位点预测的、做酶工程突变设计的、以及想在实验室已有比对数据上快速筛选候选位点对的研究者。本文会从指标选型讲到输入清洗再到最小跑通和排错按一条能复现的路径走完。2. 协同进化分析的三种核心指标MI、SCA、DCA 与工具箱选型2.1 互信息、SCA、DCA 的差异根据序列量选择方法先理清一个容易混淆的点协同进化不是一个算法而是一类问题的统称。AKtoolbox 这类工具箱内部通常实现了多种统计指标最常见的三种是互信息Mutual Information, MI、统计耦合分析Statistical Coupling Analysis, SCA和直接耦合分析Direct Coupling Analysis, DCA。三者算的都是“位点对分数”但出发点和适用条件差别很大。互信息是最直观的一种。它把每个比对位点看作一个离散变量统计两个位点上氨基酸残基的联合分布与独立分布的偏差。如果两个位点各自独立联合概率应该约等于两个边缘概率的乘积如果存在共变联合频率会系统性偏离这个乘积MI 就是这种偏离程度的度量。MI 的优点是不假设线性关系、实现简单、计算快缺点也突出——样本量偏置明显序列太少时分数系统性虚高而且它分不清“直接耦合”和“通过中间位点传播的间接关联”。SCA 则不是算两两位点对的简单相关而是把每个位点的保守性熵加权后构建一个协方差矩阵再做特征分解从中提取若干个“模式”每个模式对应一组功能上耦合的残基网络。这个方法特别适合找功能相关残基簇比单看位点对更有整体视角。代价是它对输入比对的物种采样和清洗质量要求较高另一条典型工序是它输出的是一个子空间新手刚接触时理解成本略高。DCA 是当前做接触预测最主流的一类方法。它用全局概率模型Potts 模型或马尔可夫随机场拟合整个 MSA 的残基分布模型参数里直接包含了位点之间的“直接”耦合强度从而能把间接传播的信号剥离掉。DCA 在小样本下比 MI 稳在序列数充足时做接触预测的精度明显高于 MI 和 SCA代价是计算量大、参数多需要调正则化系数和迭代次数。指标核心思想序列量要求主要输出典型场景MI两列联合分布 vs 独立分布较低但偏差大位点对分数矩阵快速筛选候选SCA熵加权协方差 特征分解中等需物种采样均匀功能模式/残基网络功能位点网络分析DCA全局 Potts 模型拟合较高建议 ≥500 条直接耦合分数结构接触预测我一般会按数据量决定入口序列少于 50 条时MI 只能当粗筛必须配合置换检验50 到 200 条之间优先用 SCA 或者归一化互信息到了几百上千条且目标是接触预测直接上 DCA。AKtoolbox 里如果同时提供了多种方法起始阶段我建议先用 MI 跑通全流程确认数据质量再切到更重的模型——因为流程上的坑主要在前处理方法本身反而好换。2.2 为什么用开源 Matlab 工具箱落地可读源码、交互调试与矩阵运算先解决一个常被问的问题同样做协同进化为什么不用 PythonPython 生态确实丰富比如常见的生物信息库和 DCA 实现都有现成包但 Matlab 在某些场景下依然更顺手尤其是面向老实验室的既有环境。第一是矩阵操作的表达效率。协同进化分析的本质是大型两两矩阵运算从频率表、协方差矩阵到特征分解Matlab 里都是原生操作写起来接近数学公式本身。热图可视化也是直接imagesc一行出图交互式命令行让“算一步看一步”的调试体验很好。第二是开源的透明度。AKtoolbox 是开源工具箱意味着你可以直接打开每个函数的源码搞清楚某个分数到底是怎么算出来的也可以按自己的数据情况修改背景校正逻辑。相比之下商业封闭工具是一个黑匣子结果异常时只能发邮件问技术支持这在方法探索阶段是很被动的。开源项目的常见做法是先跑通自带示例再单步调试进核心函数把输出和文献中的定义逐一对应。第三是环境门槛。很多做分子进化和结构生物的实验室本来就有 Matlab 授权没有 Python 环境或者不想维护虚拟环境。AKtoolbox 只需要把工具箱目录加入 path 即可运行依赖面很小适合拿到数据马上出结果的工作流。边界也要说清楚AKtoolbox 开源不等于 Matlab 免费你仍然需要一个可用的 Matlab 环境。另外如果你的 MSA 有几万条序列、位点数上千全矩阵 MI 的双重循环在 Matlab 里会明显变慢这时候要考虑分块计算或换用向量化写法这点后面在代码里会专门处理。3. 准备 AKtoolbox 的输入多序列比对的三步清洗与数值化3.1 清洗规则与经验阈值去冗余、gap 列过滤与序列数下限很多人拿到一个 MSA 就直接跑协同进化出来的热图看着很漂亮结果经不起推敲。原因很简单协同进化是对“统计偏差”的度量输入数据里任何系统偏差都会被放大成假信号。所以在进入 AKtoolbox 之前有三步清洗是必须做的。第一步是序列去冗余。如果比对里包含大量近缘序列比如同一个物种亚种或同一个家族内部的近亲这些序列会携带共同的系统发育信号导致某些位点对“看起来共变”实际上只是共享了祖先状态。常见的做法是先计算序列两两之间的不一致比例p-distance再按阈值剔除与代表序列过于相似的序列。经验阈值我一般设在 0.15也就是两条序列的残基差异小于 15% 时只保留一条代表序列如果数据量本身少可以放宽到 0.1保留更多多态性。第二步是 gap 列过滤。多序列比对里不可避免会有插入缺失形成大量 gap。gap 不是氨基酸残基如果直接进入频率统计会把“缺失”误当成一种状态参与共变计算产生虚假关联。常见规则是删除 gap 比例超过 30% 的列即maxGapFrac 0.3。如果比对质量差这个值可以调严到 0.2但如果序列本身长度差异大太严会把有效信息也删掉。第三步是确认序列数的下限。协同进化依赖频率估计样本量不够时任何统计指标都不稳。我的经验是最少 50 条高质量同源序列少于 30 条时结果基本只能当参考不能用来下功能结论。序列数不够时优先做的是回数据库扩充同源序列而不是硬跑分析。提示去冗余和序列数下限是矛盾的。去冗余阈值越大保留序列越少序列数不足时反而要降低去冗余阈值来保留更多样本。先看总量再决定冗余阈值。3.2 用 Matlab 把 fasta 比对转成数值矩阵AKtoolbox 的输入通常是 fasta 格式的比对文件但内部计算需要的是数值矩阵每一行是一条序列每一列是一个比对位点20 种标准氨基酸编码为 1 到 20gap 和其他不确定残基统一编码为 0。下面这个函数读取 fasta 并返回字符矩阵% readFastaAlign.m —— 读取比对后的 fasta返回字符矩阵和序列名 function [align, names] readFastaAlign(filename) fid fopen(filename, r); if fid 0 error(打不开文件%s, filename); end lines {}; while true tline fgetl(fid); if ~ischar(tline), break; end lines{end1, 1} tline; %#okAGROW end fclose(fid); names {}; seqs {}; idx 0; for i 1:numel(lines) line strtrim(lines{i}); if isempty(line), continue; end if line(1) idx idx 1; names{idx, 1} strtok(line(2:end)); %#okAGROW 只取注释行第一个词作为序列名 seqs{idx, 1} ; else seqs{idx, 1} [seqs{idx, 1}, upper(line)]; %#okAGROW 处理多行序列 end end lens cellfun(numel, seqs); if numel(unique(lens)) 1 error(不是合法的比对文件各序列长度不一致); end align vertcat(seqs{:}); end逻辑说明函数逐行读取 fasta开头的行为序列名连续的非空行拼成一条序列。比对文件允许序列人为分段换行所以这里必须做拼接而不是按行直接进矩阵。最后检查所有序列长度是否一致这是 MSA 合法性的基本判断。参数说明filename是比对文件路径输出align是字符矩阵行对应序列列对应比对位点names是对应的序列名元胞数组。如果读入的文件不是比对而是未比对序列长度检查这一步会直接报错这是好事避免后期在错误数据上白跑。接下来把字符矩阵转换成数值矩阵并同时做 gap 列过滤% 把字符比对矩阵编码为整数矩阵 % 20 种标准氨基酸映射到 1-20gap 及其他字符编码为 0 function X encodeAminoAcids(align) [nSeq, nPos] size(align); X zeros(nSeq, nPos); aminoAcids ACDEFGHIKLMNPQRSTVWY; for p 1:nPos for s 1:nSeq thisChar align(s, p); idx strfind(aminoAcids, thisChar); if ~isempty(idx) X(s, p) idx; else X(s, p) 0; end end end end % 删除 gap 占比超过 maxGapFrac 的列 % 返回过滤后的矩阵和列映射表 colMap function [Xf, colMap] filterSites(X, maxGapFrac) nSeq size(X, 1); gapFrac mean(X 0, 1); keep gapFrac maxGapFrac; Xf X(:, keep); colMap find(keep); end逻辑说明encodeAminoAcids用 20 个标准字母的查表方式做映射遇到 B、Z、X、*、- 这类字符一律编码为 0。filterSites按列计算 gap 比例保留不超过阈值的列。colMap是过滤后的列号到原始比对列号的映射后面做结构对照时靠它还原位点编号。参数说明maxGapFrac我默认给 0.3。双层循环在序列数上万时会慢实际使用可以改成先double(align)再做ismember映射效果相同但快一个数量级。注意这里只过滤了列还没做序列去冗余去冗余放在编码之后因为数值矩阵可以直接用pdist算距离% 贪心去冗余按顺序取代表序列把与它距离小于 distThreshold 的序列剔除 % distThreshold 是 p-distance即不一致残基比例阈值 function [Xred, keepIdx] dedupByDistance(X, distThreshold) nSeq size(X, 1); d squareform(pdist(X, hamming)); keep true(nSeq, 1); for i 1:nSeq if ~keep(i), continue; end same d(i, :) distThreshold; same(1:i) false; % 不重复判断自己和更早的序列 keep(same) false; end Xred X(keep, :); keepIdx find(keep); end逻辑说明贪心策略顺序扫描选中第一条未被剔除的序列作为代表把所有与它距离小于阈值的后续序列全部标记为剔除。这样保证任意两条被保留的序列之间的距离都大于等于阈值。pdist(X, hamming)计算的是不一致字符比例编码为 0 的 gap 会被当成一种值参与比较所以这一步必须放在filterSites之后否则 gap 会虚增序列间距离。参数说明distThreshold 0.15意味着两条序列一致度低于 85% 才同时保留想要更严格的去冗余阈值调到 0.1。去冗余会直接减少参与统计的序列数序列总量本来就不够 50 时建议跳过这一步或把阈值放宽到 0.2。4. 跑通 AKtoolbox 最小示例从互信息分数到协同热图4.1 最小可运行脚本与参数设置安装工具箱的常规做法是把整个目录加入 Matlab path然后在脚本里调用。下面是一个完整的最小示例数据准备沿用上一章的清洗流程核心分析先用本地自实现的互信息函数跑通装好 AKtoolbox 后把对应行替换成工具箱自带分析函数即可。% AKtoolbox 最小示例数据清洗 - 互信息 - 热图 addpath(genpath(AKtoolbox)); % 将工具箱目录加入搜索路径 alignFile demo_alignment.fasta; maxGapFrac 0.3; distThreshold 0.15; nPerm 0; % 置换次数先设 0 快速出原始 MI % 第 1 步读取与编码 [align, names] readFastaAlign(alignFile); X encodeAminoAcids(align); % 第 2 步列清洗与序列去冗余 [Xf, colMap] filterSites(X, maxGapFrac); [Xred, idxKeep] dedupByDistance(Xf, distThreshold); % 第 3 步互信息分析 % 使用 AKtoolbox 时把下一行替换为工具箱自带的分析模块 % 参数名以它源码注释为准本示例用本地实现保证可复现 [MI, Z] computeMIwithZscore(Xred, nPerm); % 第 4 步可视化 figure; imagesc(MI); colormap(parula); colorbar; xlabel(过滤后的位点编号); ylabel(过滤后的位点编号); title(Pairwise Mutual Information); % 保存结果colMap/idxKeep 用于还原原始位点编号 save(coevolution_result.mat, MI, Z, colMap, idxKeep, names);逻辑说明脚本按“清洗 → 分析 → 可视化”三阶段组织。第 1、2 步就是上一章的预处理第 3 步计算互信息矩阵第 4 步用imagesc画热图。colMap和idxKeep在保存结果时必须带上否则后续做结构映射时无法知道第 10 列对应原始比对的哪一列。参数说明nPerm 0表示不做置换检验只输出原始 MI适合先确认流程能跑通数据量小的时候这一步几秒就能出图。正式分析时把nPerm设为几百甚至上千但全矩阵置换很慢我会先筛高分层再对候选位点对做置换具体在第 6 章展开。互信息函数的核心实现如下它同时支持置换检验function [MI, Z] computeMIwithZscore(X, nPerm) [nSeq, nPos] size(X); alph 20; % 标准氨基酸字母表大小 MI zeros(nPos, nPos); Z zeros(nPos, nPos); for i 1:nPos-1 col_i X(:, i); valid_i col_i ~ 0; for j i1:nPos valid valid_i (X(:, j) ~ 0); xi col_i(valid); xj X(valid, j); obs mutualInfoPair(xi, xj, alph); MI(i, j) obs; MI(j, i) obs; if nPerm 0 nulls zeros(nPerm, 1); for k 1:nPerm rp randperm(numel(xj)); nulls(k) mutualInfoPair(xi, xj(rp), alph); end mu mean(nulls); sd std(nulls); if sd 0 Z(i, j) (obs - mu) / sd; Z(j, i) Z(i, j); end end end end end function mi mutualInfoPair(x, y, alph) N numel(x); if N 2, mi 0; return; end pxy accumarray([x, y], 1, [alph alph]) / N; px sum(pxy, 2); py sum(pxy, 1); pq px * py; % 独立假设下的期望联合分布 mask pxy 0; mi sum(pxy(mask) .* log(pxy(mask) ./ pq(mask))); end逻辑说明mutualInfoPair用频率表估计联合分布和边缘分布然后计算 KL 散度形式的互信息。accumarray把每一对残基出现的次数累加成二维计数矩阵除以总序列数得到联合频率。px * py是独立假设下的期望联合分布当实际联合分布偏离期望时MI 为正。参数说明alph 20对应 20 种标准氨基酸gap 已经在进入函数前通过valid过滤掉。如果想让 MI 输出单位是比特而不是 nats把log换成log2分数会整体乘以约 1.44 倍。这里实现的是经典频数版 MIAKtoolbox 自带版本通常还会加入小样本校正和背景模型效果更好但接口本质上返回的是一样的得分矩阵。4.2 输出解读MI 值、归一化与热图模式跑通只是第一步更要紧的是会读输出。MI 矩阵是个对称方阵第 i 行第 j 列表示位点 i 和位点 j 的协同进化强度。注意 MI 的绝对数值受序列数和位点变异程度影响直接用原始值跨数据比较没有意义要比较至少先归一化。常见的归一化做法是计算归一化互信息NMI公式为NMI(i,j) MI(i,j) / sqrt(H_i * H_j)其中H_i是第 i 列的香农熵。这样做可以把变异程度高的位点天然的高 MI 拉平突出真正超越单点变异度的共变信号。在自习阶段我会在脚本里加一行NMI MI ./ sqrt(entropyMat * entropyMat)先看 NMI 再看 MI。热图怎么读也有门道。健康的协同进化热图不会全盘高亮而是呈现稀疏的亮点或块状结构。对角线附近的块状高亮往往对应同一结构域内的残基接触远离对角线的点则可能对应长程相互作用或变构通路。如果你看到热图上大片连续的红色区域第一反应不该是“发现重要功能网络”而应该怀疑输入数据有问题——这在第 5 章详细展开。建议保存结果时附带一张参数清单方便追溯输出变量含义用途MI原始互信息矩阵排序、筛选候选位点对ZZ-score 矩阵置换检验版判断显著性Z3 可作候选colMap过滤后列到原始比对列的映射结果映射到序列/结构idxKeep去冗余后保留的序列索引回溯样本构成names保留序列名核对物种采样是否均匀5. AKtoolbox 实战排查5 个高频坑的定位与修正5.1 序列数太少得分矩阵全是噪声现象不到 30 条序列的比对跑出 MI 热图结果整张图都偏高随机挑两个位点分数都很大重复抽几次序列算出的结果完全不稳定。原因互信息是频率估计联合概率表有 20×20400 个格子序列太少时每个格子里的样本少得可怜小样本偏差让 MI 系统性虚高。热图越“红”越说明统计不够而不是信号强。解决先把序列数补到 100 条以上再分析。如果没法扩充就把去冗余阈值从 0.15 放宽到 0.2保留更多近缘序列来撑样本量同时用置换检验给每个位点对算 Z-score而不是直接信原始 MI。请记住一个经验判断序列量小时任何协同进化分数都带有“美丽噪声”结构验证比分数本身更有说服力。5.2 冗余序列制造虚假协同簇现象结果里出现一个非常漂亮的高亮“模块”这些位点在功能注释上完全说不上关联但分数极高。仔细看输入序列发现大部分序列来自少数几个近缘分支采样严重不均。原因这是系统发育信号在作怪。近缘序列共享同一个祖先状态它们会同时“继承”一组相同的残基组合统计上看起来像协同进化实际只是谱系关联。MI 和 SCA 都无法自动排除这种假象只有靠输入清洗解决。解决严格做序列去冗余保证保留序列之间的差异度在阈值以上更关键的是检查物种背景尽量让采样覆盖不同分支而不是集中在一两个近缘类群里。我用dedupByDistance时会把distThreshold从 0.15 收紧到 0.1并输出保留序列的名单人工看一眼。5.3 gap 列混进统计个别位点虚高现象某几个位点的 MI 异常高拉开原始比对一看这一列大部分是 gap只有少量残基。原因gap 编码为 0 后如果把 0 当成普通状态参与统计缺失就会被当成一种“共同特征”。两条序列在同一个位置同时缺失会贡献一次虚假共变gap 比例越高这个效应越强。解决在列过滤阶段把 gap 比例超标的列删掉maxGapFrac设在 0.3进入mutualInfoPair前再用valid过滤掉当前位点对上有 gap 的序列确保统计完全在真实残基上进行。这两个动作缺一不可列过滤管高频 gap 列行过滤管残留的序列缺失。5.4 清洗后位点编号对不上结构坐标现象热图上排名靠前的位点对拿到结构上验证发现位置完全对不上反复检查后确认是编号错位。原因filterSites删除了若干列矩阵列号已经是重新编号的没有映射回原始比对位置。如果直接用过滤后的列号去结构文件里找残基第 10 列可能实际上是原始比对第 15 位。解决从一开始就保留colMap触碰结果矩阵的任何位点对都先用origPos colMap(i)还原。保存结果时把colMap和idxKeep一起落盘不要只存 MI 矩阵。读结构文件做验证时用原始编号去匹配这一步能节省大量返工时间。5.5 忽略背景校正高相关位点全是假阳性现象所有位点对的 MI 都显著特异性很差和随机挑选的位点对没有区分度。原因原始 MI 受每个位点自身的变异程度影响一个位点本身熵很高它和任何位点算出来的 MI 都会偏高。这不是真实的协同进化而是单点变异的背景叠加。解决至少要做归一化折。先用 NMI 或 Z-score 替换原始 MI 做排序再看 AKtoolbox 是否提供背景模型选项比如把序列随机打乱后计算背景分布然后做差。DCA 类方法内置了背景校正如果数据量允许这类方法在特异性上会明显优于裸 MI。这个坑最隐蔽因为它不报错只是结果“看起来合理但经不起验证”。6. 验证结果置换检验与结构接触对照6.1 置换检验对高分层位点对做显著性评估全矩阵的置换检验计算量太大我一般只对排名靠前的位点对做。具体做法把其中一个位点的残基顺序随机打乱重算 MI重复多次得到零分布再比较观测 MI 在零分布中的位置。% 对高分层位点对做置换检验 load coevolution_result.mat; % 读入前面保存的 MI、Xred 等 topIdx find(triu(MI, 1) prctile(MI(triu(MI, 1)), 95)); nPerm 1000; pvals zeros(numel(topIdx), 1); for k 1:numel(topIdx) [i, j] ind2sub(size(MI), topIdx(k)); xi Xred(:, i); xj Xred(:, j); valid xi ~ 0 xj ~ 0; xi xi(valid); xj xj(valid); obs mutualInfoPair(xi, xj, 20); nulls zeros(nPerm, 1); for t 1:nPerm rp randperm(numel(xj)); nulls(t) mutualInfoPair(xi, xj(rp), 20); end pvals(k) (sum(nulls obs) 1) / (nPerm 1); % 加 1 避免 p0 end逻辑说明triu(MI, 1)取严格上三角避免重复和自相关。候选对按 MI 前 5% 选取。置换检验的原假设是两个位点独立打乱其中一列的残基就破坏了共变关系得到的是独立假设下的 MI 分布。pvals(k)是经验 p 值分母用nPerm1是标准做法防止出现无法解释的 0。参数说明nPerm 1000是精度和耗时的折中正式发表建议至少 5000 次。置换检验对每对候选位点都要循环若干次候选数超过 500 时会比较慢可以先按前 1% 筛。这一步结束后把 p 值小于 0.01 的位点对列为最终候选。6.2 用三维结构接触对做交叉验证高分层算出来之后最硬核的验证是看它是否对应真实的空间接触。如果你的项目里有同源蛋白的三维结构文件提取残基间距离小于 8Å 的位点对作为阳性集合然后检查你的高分层富集了多少真实接触。% 结构接触对照用距离 8Å 的残基对验证候选 load coevolution_result.mat; contacts load(structure_contacts.mat); % 每行两个原始位点编号 for k 1:numel(topIdx) [i, j] ind2sub(size(MI), topIdx(k)); origPos sort([colMap(i), colMap(j)]); isContact(k) ismember(origPos, contacts.pairs, rows); end hitRate mean(isContact); % 随机抽取同样数量的位点对作为对照 randIdx randi([1, numel(MI)], numel(topIdx), 1); % 对比 hitRate 与随机抽取对照的命中率应显著更高逻辑说明colMap在这里体现出价值——把过滤后的位点对还原到结构编号。如果高分层里真实接触的比例明显高于随机抽取对照说明协同进化信号和空间邻近一致结果可信度高。参数说明8Å 是常见的接触定义阈值α 碳之间的距离更保守可以用 6Å。同一个数据集建议两个阈值都试一遍如果结论一致稳定性更好。结构文件里的残基编号体系与比对编号未必相同通常需要先做序列号对应这也是最容易出错的一步建议单独写一个映射表核对。我自己最早跑协同进化的时候最大的翻车就是把公共序列库拉回来的数据直接塞进分析热图上那些漂亮的高亮块后来发现全是谱系效应。从那以后我养成了一个雷打不动的习惯先清理输入、再置换检验、最后拿结构说话三条线缺一不可。希望帮到你。本文还有配套的精品资源点击获取