简介这是一款开源的Matlab工具箱面向蛋白质多序列比对的协同进化分析采用简化的BSD许可证发布并且不依赖Matlab自带的生物信息工具箱。工具箱集成了统计耦合分析、直接耦合分析、互信息、显式似然等多种主流协同进化计算方法每种方法均以独立函数形式封装用户可按需调用便于对不同算法的计算结果进行横向比较同时工具内还附带多种生物信息学矩阵档案为后续的深度分析提供基础。整个压缩包共包含116个文件其中核心为103个Matlab脚本另有支持Windows、macOS与Linux平台的动态链接库、可执行文件以及蛋白质结构、二级结构、序列等示例数据文件并包含C源码以便有需要的用户进行底层修改。压缩包体积仅为264KB轻量小巧下载和部署都非常方便。示例数据中提供了蛋白质序列与二级结构文件能够帮助使用者快速验证整个分析流程。目前已有13958人学习或下载特别适合从事蛋白质结构预测、进化保守性研究以及相关算法开发的科研人员和Matlab开发者。1. 从多序列比对里找“一起变”的残基AKtoolbox 是什么、为什么值得拆我刚接触多序列比对的时候总觉得分析到保守位点就完事了。后来在某个模拟项目X里做功能残基筛选被一个问题卡住两个位点单独看都不够保守但它们的变异模式总是成对出现这种信号在传统的序列Logo图里根本看不见。AKtoolbox就是为解决这个问题而生的开源软件一个纯MATLAB工具箱专门做多序列比对的协同进化分析。它把FASTA读取、位点过滤、互信息/统计耦合计算、置换检验、结果导出串成一条完整链路尤其适合手里没有Python生信环境、却已经在用MATLAB做数据分析的工程师和研究生。如果你手头已经有一份比对好的FASTA从敲第一行命令到拿到残基配对得分矩阵通常用不了一个上午。2. 协同进化分析为什么选 MATLAB先想清楚提取的是哪一种“联动信号”2.1 协同进化信号来自哪里MSA 给我们的原始输入是什么多序列比对本质上是一张字符表每一行是一条序列每一列是比对上的一个位点字符可能是氨基酸残基也可能是-这样的空位。我们看一个位点是否保守通常会统计这个位点上出现频率最高的残基占多大比例但这是单点统计它丢掉了一个关键信息——位点之间是否存在“组合偏好”。举个最简单的例子。假设一个假想MSA里有100条序列位点12和位点80单独看都各有约一半序列出现两类残基单点保守性得分都不高。可当你把两个位点放在一起看时会发现规律非常强只要位点12是K位点80几乎一定是D只要位点12是E位点80几乎一定是R。这说明两个位点根本不是独立变异的它们很可能在三维结构上靠在一起或者共同受到同一选择压力一个变了另一个必须跟着补偿。协同进化分析想提取的就是这种成对非随机联动。它把每个位点当作一个离散变量把每一对位点看作一个二维列联表问题判断观测到的联合分布是否显著偏离两个边缘分布的乘积。理解这一点之后后面所有预处理参数都有了方向我们希望在MSA里尽可能去掉噪声留下真正能反映“一起变”的信号。2.2 工具箱里的三类核心算法MI、SCA 与扰动统计AKtoolbox把主流的协同进化算法封装成了几个统一入口你不需要把每个算法从头推导一遍但必须清楚它们各自的输出含义否则很容易拿错结果去解释生物学问题。第一类是互信息也就是 MI。它比较的是每个位点对联合分布p(x,y)与两个位点边缘分布乘积p(x)p(y)的差异差异越大互信息值越高。MI实现最简单、计算最快适合做初步筛选。但它的缺点也很明显天然偏向保守位点。因为保守位点本身只有少数几个字符类别概率估计相对稳定一旦出现少量共变残差就会被放大而高度可变的位点类别太多概率被摊薄即使存在真实共变也容易被淹没。第二类是统计耦合分析也就是 SCA。它不会直接给出二维的“位点A到位点B”的分数而是先对每个位点做熵加权频率归一化然后对加权频率矩阵做特征分解再通过谱重构得到位点之间的耦合关系。SCA对保守位点有一定惩罚所以更容易找出结构域内部的功能模块适合分析较大范围的残基网络。但它比较“吃数据”要求MSA序列数和序列多样性都足够如果只拿20条序列去跑SCA出来的特征向量基本是噪声。第三类可以叫扰动统计法风格上和ELSC、McBASC这一类方法接近。它的思路是把每条序列当作一次自然扰动用回归或协方差模型分离出其他位点对当前位点变异的贡献然后做显著性检验。这类方法比MI更“条件化”能消除一部分背景相关但计算量最大而且对MSA里的异常序列特别敏感如果里面混进几条和外类群序列差异极大的序列它们可能主导整个协方差估计。下面这张表可以帮助选型算法类型输出形式适合场景主要风险互信息 MI位点对之间的非负分数快速初筛找候选配对偏向保守位点gap处理不当容易假高统计耦合 SCA位点谱权重与重构的耦合矩阵找结构域协同模块、功能界面相关残基需要足够多样本小样本下不稳定扰动统计残差或检验统计量中长序列、需要定量显著性时对异常序列敏感计算耗时长2.3 为什么在 MATLAB 里做这件事以及 AKtoolbox 的定位很多人会先问现在Python生态里做生物信息学已经很成熟为什么还要用MATLAB我在实际使用中的体会是很多实验室和工程团队的老课题代码都是MATLAB写的交互界面上也习惯用MATLAB做热力图和散点图如果为了一个协变分析再搭一套Python环境先要解决BioPython版本、numpy版本这些依赖冲突可能比分析本身更耗时。AKtoolbox把从FASTA读取到统计检验的完整流程都收进了一组前缀为ak_的函数。它不提供多序列比对功能输入必须是已经比对好的MSA这意味着它定位的是一个中段分析组件前面接任何比对软件后面接MATLAB的可视化和统计分析。作为开源软件另一个好处是当你对某个结果存疑时可以直接打开ak_covar这类函数的源码确认里面到底有没有做gap masking、有没有去冗余而不是把算法当黑匣子用。3. 环境准备与数据格式把 FASTA 变成 MATLAB 喜欢的形状3.1 拉源码并登记路径我手头这版 AKtoolbox 是纯命令行函数式架构没有 App Designer 界面目录下主要是ak_前缀的脚本函数、一个测试目录和几个示例数据文件。没有界面反而是好事命令行接口稳定写批处理也方便。拿到源码后第一步不是马上跑demo而是先把工具箱路径加到 MATLAB 搜索路径下并确认核心函数能被识别。% 假设源码解压后放在 ~/src/AKtoolbox addpath(genpath(~/src/AKtoolbox)); % 检查核心函数是否可见 which ak_readfasta这段代码里genpath的作用是递归添加 AKtoolbox 下所有子目录避免漏掉嵌套的工具函数目录。执行which ak_readfasta后如果返回的是完整路径说明路径配置正常如果返回“未找到”先检查是否用了addpath(genpath(...))而不是只addpath外层目录这是最常见的漏配原因。3.2 从 FASTA 到字符矩阵AKtoolbox 的输入函数把 FASTA 文件读进来后会保留每条序列的头部信息和序列内容再用另一个函数把变长序列转成统一的字符矩阵。我不建议自己用fastaread之后直接堆字符串因为序列长度不一致时很容易在拼接处产生隐性错误。% 读入已经比对好的多序列比对文件 fastaData ak_readfasta(A3L1_aln.fasta); % 转换成字符矩阵行序列列MSA位点 seqMat ak_alignment2matrix(fastaData); % 看一下规模序列数和位点数 [ns, np] size(seqMat); fprintf(序列数%d位点数%d\n, ns, np);这里的fastaData在函数内部被处理成一个结构数组每个元素包含头部和序列体。ak_alignment2matrix做的是字面量转换把每一条序列展开成一行字符不足的长度用-补齐如果输入的多序列比对里原本就有gap这些gap会原样保留。拿到seqMat之后后续所有函数都按这个矩阵接收数据这是整个分析链路的统一数据格式。3.3 关键预处理过滤高gap列并对近缘序列去冗余多序列比对不是拿来就能直接算协同进化的。很多MSA里存在这样的列只有少数序列有残基大部分序列在这个位置是-这种列带来的“变异信号”根本不是进化约束而是比对软件对插缺事件的占位。另外如果数据集里包含大量同源度极高的近缘序列比如同一菌株重复测序得到的序列它们会把谱系偏好放大让少数远缘序列的变异信息被淹没。常见的做法是先做一次位点过滤和去冗余。AKtoolbox里对应函数是ak_filter_alignment% 过滤高gap列并对相似度高于90%的序列去冗余 [alnClean, info] ak_filter_alignment(seqMat, ... max_gap_frac, 0.3, ... min_identity, 0.9); % info里保存了保留下来的是原始MSA哪些列 origPos info.positions; fprintf(过滤后保留位点数%d\n, size(alnClean, 2));max_gap_frac表示一列中gap占比超过30%就删除min_identity表示序列相似度高于90%时只保留其中一条代表性序列。我自己的习惯是如果MSA本来序列数就不多min_identity可以放到0.95甚至1.0优先保证样本量如果序列冗余严重则改成0.8。这里还有一个必须养成的习惯info.positions一定要保存下来。一旦过滤函数执行完列号就变成了“过滤后新矩阵的列号”而原始MSA的列号存在info.positions里后面做结构映射时全靠它。3.4 数据规模不同时的参数决策预处理没有一组参数能适配所有数据按照我跑过的几个MSA案例粗估可以按下面的思路选数据情况建议参数理由序列数很多超过200条min_identity设为0.8~0.9去掉重复代表避免近缘序列主导序列数偏少只有30~50条min_identity设为0.95以上或不去冗余尽量保留样本量否则统计不稳定比对覆盖很差gap比例高max_gap_frac设为0.3先砍掉大量无信息的占位列比对覆盖很好gap很少max_gap_frac设为0.5只删除那些几乎全是gap的区域这里不用把这些值当成铁律它们解决的是同一个问题让矩阵里的每一个统计量都尽量来自真实的残基状态而不是来自比对噪声。4. 跑通一次默认流程拿到互信息得分矩阵和显著性估计4.1 最小运行命令数据预处理完之后就可以直接调用协同进化分析入口。AKtoolbox 里最常见的入口函数是ak_covar你不需要分别调用十几个子函数只要把字符矩阵放进去并指定算法即可。res ak_covar(alnClean, ... method, mi, ... normalize, zscore, ... nperm, 1000);method参数指定算法mi表示互信息也可以换成sca或perturb。normalize指定标准化方式zscore会基于置换检验的随机分布把原始分数转成z分数。nperm是置换次数这里设置1000次意味着对每个位点对都做1000次随机打乱的统计量估计用于后面计算经验p值。这里的逻辑是先计算观测数据的原始互信息矩阵然后把每一列内的残基字符随机打乱破坏位点之间的真实联动再重新计算一次随机情况下的互信息。重复nperm次之后就能得到随机背景的均值和标准差再用观测值减去随机均值、除以标准差得到标准化后的得分。这个得分比原始互信息更有可比性因为它已经减掉了背景偏差。4.2 返回结构里有什么ak_covar返回的res是一个 MATLAB struct里面主要字段如下字段含义res.score标准化后的位点对得分矩阵对称矩阵行列数等于输入位点数res.zscore对应的z分数矩阵用于评估偏离随机背景的程度res.pvalue置换检验得到的经验p值矩阵res.positions这些位点对应原始MSA中的列号等价于之前info.positionsres.method本次运行使用的算法名第一次跑完建议先看size(res.score)确认矩阵维度等于预处理后的列数而不是原始MSA列数。如果维度对不上多半前面某一步没有把矩阵传对。4.3 快速看最强的位点对得分矩阵是对称矩阵第i行第j列和j行i列是同一个配对因此只看上三角就够。下面这段代码可以打印出得分最高的前10个位点对score res.score; % 只取上三角排除对角线 [I, J, V] find(triu(score, 1)); % 按得分降序排列 [Vs, idx] sort(V, descend); % 打印前10个候选 for k 1:min(10, numel(idx)) rowIdx I(idx(k)); colIdx J(idx(k)); fprintf(候选配对 bit %d - %d : score %.3f z %.2f p %.3f\n, ... rowIdx, colIdx, Vs(k), ... res.zscore(rowIdx, colIdx), ... res.pvalue(rowIdx, colIdx)); end这段代码里的triu(score, 1)会保留矩阵对角线以上元素第二个参数1表示从对角线右侧开始取直接去掉对角线。find返回的三个输出分别是行列下标和非零值。排序后打印时行的索引和列的索引要一起从I(idx(k))和J(idx(k))取避免排序后错位。看到的结果里如果某个配对的p值特别低但z分数却不高要留意是否来自低频残基组合后面进阶部分会提到怎么处理。4.4 参数组合怎么调实际分析很难一次到位的下面是我通常使用的参数组合参数推荐值说明method先mi后scaMI做快速初筛SCA做模块验证normalizezscore减掉随机背景比原始分数可解释nperm1000追求统计稳定性时使用nperm100调参阶段先用快速模式有一点要特别注意切换method后不同算法的得分量纲完全不一样SCA的权重甚至会出现负值不能拿MI的阈值去直接套SCA结果。调参阶段我会先把nperm降到100甚至0这里nperm0表示只算原始分数不做置换让流程快速跑通确认矩阵尺寸和输出字段没问题后再把置换次数提上去。5. 避坑与常见问题让我返工的五个具体失误5.1 没有去冗余结果看起来“太完美”现象跑完IM后排名靠前的位点对几乎都落在一个非常保守的结构域上并且这些位点本身单独看也很保守和前几个位点的互信息高得像假数据一样。原因MSA里如果包含大量相似度超过95%的近似序列它们在谱系上是重复抽样会让少数几个常见残基组合反复出现互信息被这些重复样本刷高最终主导全局排名。协同进化分析真正需要的是多样化的独立序列而不是同一谱系的一堆副本。解决回到ak_filter_alignment把min_identity从0.9收紧到0.8必要时按构建好的进化谱系进行加权抽样。我一般会先打印去冗余前的序列数size(seqMat,1)和去冗余后size(alnClean,1)如果发现只去掉了不到10%的序列说明数据本身多样性不足即使算出来结果也最好不要直接下强结论。5.2 把 gap 当成一个普通字符现象得分矩阵里有一整块连续列在结果中排名特别高比如位点120到150。单独看它们的原始互信息发现这些列在两两之间几乎都是最高分但和功能没有任何关系。原因这些列在MSA里是一大片gap区许多序列在这个区域整体缺失gap字符因此成对共现。在互信息计算眼里这等价于“有大量序列在同一个位置缺失”于是统计量表现出强关联但这根本不是进化意义上的协同变异。解决用max_gap_frac先把高gap列过滤掉。如果某个区域没有完全过滤干净可以额外保留一个gap掩码矩阵在调用ak_covar时把掩码传入至少不要让-参与最后的字符频次统计。我常用的做法是先把max_gap_frac定在0.3再人工看info.positions是否还有连续成片被保留的gap区域。5.3 过滤后的位点编号和结构编号对不上现象结果显示第42位和第88位是显著共进化残基对但把这两个编号直接对应到三维结构时发现结构里该区域根本不存在这两个残基或者对应到了完全不同位置的残基上。原因MSA本身不是从结构编号开始的。原始比对里可能删掉了一些空位而ak_filter_alignment又基于旧的MSA矩阵删除了高gap列所有后续得到的编号都是“当前矩阵的第几列”而不是原始序列上的残基编号更不是PDB文件里的编号。解决从一开始就保存info.positions每个分析结果里的位点编号都要通过它做一次映射把“当前矩阵列号”转成“原始MSA列号”再结合你读入结构文件时建立的序列到结构编号映射做第二次转换。这个映射我每次都会写进脚本第一段而不是等结果出来后再手工拼。5.4 置换检验时间爆炸现象当位点数到200以上、nperm设为1000时作业跑了几个小时还没结束。中间杀掉进程后检查发现每个位点对都要重新扫一遍矩阵计算量直接膨胀。原因nperm是对每一位点对做置换重算位点数从100增加到200时两两配对数增加到4倍如果再叠加1000次置换总体耗时就是非线性上升。AKtoolbox默认在置换时会重建整棵统计量矩阵保证p值的一致性但代价就是慢。解决调参阶段先用nperm100甚至nperm0。如果最终报告里只关心排名前几十的候选位点对可以只对候选对单独做置换检验这样计算量减少一个量级。另一个做法是使用块置换而不是逐字符置换保持序列内部相关性随机性更贴近真实进化过程。5.5 小样本硬套 SCA还拿谱去下结论现象只有20条序列跑SCA后得到的第一特征向量看起来有明确分层第一特征值解释了很大方差但换用MI方法后这些“强信号位点”大多对不上。原因SCA在构造协方差矩阵时需要足够的独立样本去估计残基之间的协变模式20条序列对L×L维矩阵来说自由度严重不足。少数几条每序列自带的特殊突变组合会被当成稳定特征放大看起来像模块化信号实际上是抽样噪声。解决少于50条序列的MSA不要做SCA先用MI做初筛序列数足够时再跑SCA并且对比两种方法得到的强候选列表是否一致。如果两个方法都指向同一组位点对结论才相对可信。小样本下看到SCA谱再有“规律”也不要直接写进报告除非用置换检验验证过。6. 进阶从得分矩阵到残基配对网络以及一个收尾习惯6.1 给得分矩阵做保守性屏蔽协同进化分数高的位点对不一定都有生物学意义尤其是互信息方法它会把保守位点之间的微弱共变放得很大。一个实用的技巧是先计算每个位点的保守程度把过于保守的位点排除掉只保留有一定变异度、又表现出联动信号的位点。% 计算每个位点的香农熵-1表示完全保守0表示完全均匀 H ak_entropy(alnClean); conserv 1 - H; % 设置一个保守性阈值太保守的位点不进入最终网络 keep conserv 0.2; scoreClean res.score .* (keep keep);这里的ak_entropy返回每个MSA位点的归一化香农熵值conserv越接近1代表位点越保守。keep是一个逻辑向量keep keep会把它扩展成一个逻辑矩阵只有两个位点都足够可变时对应位置上的得分才保留下来其他位置全部置零。这个操作能明显减少保守位点主导的问题让网络更聚焦在有真实进化约束的残基上。6.2 从矩阵到边列表拿到清理后的得分矩阵下一步是把矩阵转成网络分析工具能读的边列表。我一般输出一个三列TSV文件分别是位点A、位点B、得分。% 取上三角转成网络边表 thr prctile(scoreClean(scoreClean eps), 99); [I, J, W] find(scoreClean thr); % 写出边表三列分别表示起始位点、目标位点、权重 edges [I, J, W]; writematrix(edges, covar_edges.tsv, Delimiter, tab);阈值我通常设为得分矩阵非零部分的高分位点第99百分位数优先保留最强信号。scoreClean thr会生成一个逻辑稀疏矩阵find提取出的三列就是边表。你可以在网络可视化软件里用这个边表画残基相互作用网络也可以把它和三维结构中的残基距离做对比看那些显著共进的残基对是否在空间上确实接近。最后说一个让我改变工作习惯的教训以前我拿到新MSA就直接跑协同进化结果总要先浪费一两轮在检查数据质量和编号映射上。从那以后我每次跑AKtoolbox之前都强制自己走一遍“序列数、冗余度、gap比例、编号映射”四项检查确认无误再分析。这个习惯帮我省掉了很多返工时间希望也能帮到你。本文还有配套的精品资源点击获取 AKtoolbox：MATLAB中的多序列比对协同进化分析实践
简介这是一款开源的Matlab工具箱面向蛋白质多序列比对的协同进化分析采用简化的BSD许可证发布并且不依赖Matlab自带的生物信息工具箱。工具箱集成了统计耦合分析、直接耦合分析、互信息、显式似然等多种主流协同进化计算方法每种方法均以独立函数形式封装用户可按需调用便于对不同算法的计算结果进行横向比较同时工具内还附带多种生物信息学矩阵档案为后续的深度分析提供基础。整个压缩包共包含116个文件其中核心为103个Matlab脚本另有支持Windows、macOS与Linux平台的动态链接库、可执行文件以及蛋白质结构、二级结构、序列等示例数据文件并包含C源码以便有需要的用户进行底层修改。压缩包体积仅为264KB轻量小巧下载和部署都非常方便。示例数据中提供了蛋白质序列与二级结构文件能够帮助使用者快速验证整个分析流程。目前已有13958人学习或下载特别适合从事蛋白质结构预测、进化保守性研究以及相关算法开发的科研人员和Matlab开发者。1. 从多序列比对里找“一起变”的残基AKtoolbox 是什么、为什么值得拆我刚接触多序列比对的时候总觉得分析到保守位点就完事了。后来在某个模拟项目X里做功能残基筛选被一个问题卡住两个位点单独看都不够保守但它们的变异模式总是成对出现这种信号在传统的序列Logo图里根本看不见。AKtoolbox就是为解决这个问题而生的开源软件一个纯MATLAB工具箱专门做多序列比对的协同进化分析。它把FASTA读取、位点过滤、互信息/统计耦合计算、置换检验、结果导出串成一条完整链路尤其适合手里没有Python生信环境、却已经在用MATLAB做数据分析的工程师和研究生。如果你手头已经有一份比对好的FASTA从敲第一行命令到拿到残基配对得分矩阵通常用不了一个上午。2. 协同进化分析为什么选 MATLAB先想清楚提取的是哪一种“联动信号”2.1 协同进化信号来自哪里MSA 给我们的原始输入是什么多序列比对本质上是一张字符表每一行是一条序列每一列是比对上的一个位点字符可能是氨基酸残基也可能是-这样的空位。我们看一个位点是否保守通常会统计这个位点上出现频率最高的残基占多大比例但这是单点统计它丢掉了一个关键信息——位点之间是否存在“组合偏好”。举个最简单的例子。假设一个假想MSA里有100条序列位点12和位点80单独看都各有约一半序列出现两类残基单点保守性得分都不高。可当你把两个位点放在一起看时会发现规律非常强只要位点12是K位点80几乎一定是D只要位点12是E位点80几乎一定是R。这说明两个位点根本不是独立变异的它们很可能在三维结构上靠在一起或者共同受到同一选择压力一个变了另一个必须跟着补偿。协同进化分析想提取的就是这种成对非随机联动。它把每个位点当作一个离散变量把每一对位点看作一个二维列联表问题判断观测到的联合分布是否显著偏离两个边缘分布的乘积。理解这一点之后后面所有预处理参数都有了方向我们希望在MSA里尽可能去掉噪声留下真正能反映“一起变”的信号。2.2 工具箱里的三类核心算法MI、SCA 与扰动统计AKtoolbox把主流的协同进化算法封装成了几个统一入口你不需要把每个算法从头推导一遍但必须清楚它们各自的输出含义否则很容易拿错结果去解释生物学问题。第一类是互信息也就是 MI。它比较的是每个位点对联合分布p(x,y)与两个位点边缘分布乘积p(x)p(y)的差异差异越大互信息值越高。MI实现最简单、计算最快适合做初步筛选。但它的缺点也很明显天然偏向保守位点。因为保守位点本身只有少数几个字符类别概率估计相对稳定一旦出现少量共变残差就会被放大而高度可变的位点类别太多概率被摊薄即使存在真实共变也容易被淹没。第二类是统计耦合分析也就是 SCA。它不会直接给出二维的“位点A到位点B”的分数而是先对每个位点做熵加权频率归一化然后对加权频率矩阵做特征分解再通过谱重构得到位点之间的耦合关系。SCA对保守位点有一定惩罚所以更容易找出结构域内部的功能模块适合分析较大范围的残基网络。但它比较“吃数据”要求MSA序列数和序列多样性都足够如果只拿20条序列去跑SCA出来的特征向量基本是噪声。第三类可以叫扰动统计法风格上和ELSC、McBASC这一类方法接近。它的思路是把每条序列当作一次自然扰动用回归或协方差模型分离出其他位点对当前位点变异的贡献然后做显著性检验。这类方法比MI更“条件化”能消除一部分背景相关但计算量最大而且对MSA里的异常序列特别敏感如果里面混进几条和外类群序列差异极大的序列它们可能主导整个协方差估计。下面这张表可以帮助选型算法类型输出形式适合场景主要风险互信息 MI位点对之间的非负分数快速初筛找候选配对偏向保守位点gap处理不当容易假高统计耦合 SCA位点谱权重与重构的耦合矩阵找结构域协同模块、功能界面相关残基需要足够多样本小样本下不稳定扰动统计残差或检验统计量中长序列、需要定量显著性时对异常序列敏感计算耗时长2.3 为什么在 MATLAB 里做这件事以及 AKtoolbox 的定位很多人会先问现在Python生态里做生物信息学已经很成熟为什么还要用MATLAB我在实际使用中的体会是很多实验室和工程团队的老课题代码都是MATLAB写的交互界面上也习惯用MATLAB做热力图和散点图如果为了一个协变分析再搭一套Python环境先要解决BioPython版本、numpy版本这些依赖冲突可能比分析本身更耗时。AKtoolbox把从FASTA读取到统计检验的完整流程都收进了一组前缀为ak_的函数。它不提供多序列比对功能输入必须是已经比对好的MSA这意味着它定位的是一个中段分析组件前面接任何比对软件后面接MATLAB的可视化和统计分析。作为开源软件另一个好处是当你对某个结果存疑时可以直接打开ak_covar这类函数的源码确认里面到底有没有做gap masking、有没有去冗余而不是把算法当黑匣子用。3. 环境准备与数据格式把 FASTA 变成 MATLAB 喜欢的形状3.1 拉源码并登记路径我手头这版 AKtoolbox 是纯命令行函数式架构没有 App Designer 界面目录下主要是ak_前缀的脚本函数、一个测试目录和几个示例数据文件。没有界面反而是好事命令行接口稳定写批处理也方便。拿到源码后第一步不是马上跑demo而是先把工具箱路径加到 MATLAB 搜索路径下并确认核心函数能被识别。% 假设源码解压后放在 ~/src/AKtoolbox addpath(genpath(~/src/AKtoolbox)); % 检查核心函数是否可见 which ak_readfasta这段代码里genpath的作用是递归添加 AKtoolbox 下所有子目录避免漏掉嵌套的工具函数目录。执行which ak_readfasta后如果返回的是完整路径说明路径配置正常如果返回“未找到”先检查是否用了addpath(genpath(...))而不是只addpath外层目录这是最常见的漏配原因。3.2 从 FASTA 到字符矩阵AKtoolbox 的输入函数把 FASTA 文件读进来后会保留每条序列的头部信息和序列内容再用另一个函数把变长序列转成统一的字符矩阵。我不建议自己用fastaread之后直接堆字符串因为序列长度不一致时很容易在拼接处产生隐性错误。% 读入已经比对好的多序列比对文件 fastaData ak_readfasta(A3L1_aln.fasta); % 转换成字符矩阵行序列列MSA位点 seqMat ak_alignment2matrix(fastaData); % 看一下规模序列数和位点数 [ns, np] size(seqMat); fprintf(序列数%d位点数%d\n, ns, np);这里的fastaData在函数内部被处理成一个结构数组每个元素包含头部和序列体。ak_alignment2matrix做的是字面量转换把每一条序列展开成一行字符不足的长度用-补齐如果输入的多序列比对里原本就有gap这些gap会原样保留。拿到seqMat之后后续所有函数都按这个矩阵接收数据这是整个分析链路的统一数据格式。3.3 关键预处理过滤高gap列并对近缘序列去冗余多序列比对不是拿来就能直接算协同进化的。很多MSA里存在这样的列只有少数序列有残基大部分序列在这个位置是-这种列带来的“变异信号”根本不是进化约束而是比对软件对插缺事件的占位。另外如果数据集里包含大量同源度极高的近缘序列比如同一菌株重复测序得到的序列它们会把谱系偏好放大让少数远缘序列的变异信息被淹没。常见的做法是先做一次位点过滤和去冗余。AKtoolbox里对应函数是ak_filter_alignment% 过滤高gap列并对相似度高于90%的序列去冗余 [alnClean, info] ak_filter_alignment(seqMat, ... max_gap_frac, 0.3, ... min_identity, 0.9); % info里保存了保留下来的是原始MSA哪些列 origPos info.positions; fprintf(过滤后保留位点数%d\n, size(alnClean, 2));max_gap_frac表示一列中gap占比超过30%就删除min_identity表示序列相似度高于90%时只保留其中一条代表性序列。我自己的习惯是如果MSA本来序列数就不多min_identity可以放到0.95甚至1.0优先保证样本量如果序列冗余严重则改成0.8。这里还有一个必须养成的习惯info.positions一定要保存下来。一旦过滤函数执行完列号就变成了“过滤后新矩阵的列号”而原始MSA的列号存在info.positions里后面做结构映射时全靠它。3.4 数据规模不同时的参数决策预处理没有一组参数能适配所有数据按照我跑过的几个MSA案例粗估可以按下面的思路选数据情况建议参数理由序列数很多超过200条min_identity设为0.8~0.9去掉重复代表避免近缘序列主导序列数偏少只有30~50条min_identity设为0.95以上或不去冗余尽量保留样本量否则统计不稳定比对覆盖很差gap比例高max_gap_frac设为0.3先砍掉大量无信息的占位列比对覆盖很好gap很少max_gap_frac设为0.5只删除那些几乎全是gap的区域这里不用把这些值当成铁律它们解决的是同一个问题让矩阵里的每一个统计量都尽量来自真实的残基状态而不是来自比对噪声。4. 跑通一次默认流程拿到互信息得分矩阵和显著性估计4.1 最小运行命令数据预处理完之后就可以直接调用协同进化分析入口。AKtoolbox 里最常见的入口函数是ak_covar你不需要分别调用十几个子函数只要把字符矩阵放进去并指定算法即可。res ak_covar(alnClean, ... method, mi, ... normalize, zscore, ... nperm, 1000);method参数指定算法mi表示互信息也可以换成sca或perturb。normalize指定标准化方式zscore会基于置换检验的随机分布把原始分数转成z分数。nperm是置换次数这里设置1000次意味着对每个位点对都做1000次随机打乱的统计量估计用于后面计算经验p值。这里的逻辑是先计算观测数据的原始互信息矩阵然后把每一列内的残基字符随机打乱破坏位点之间的真实联动再重新计算一次随机情况下的互信息。重复nperm次之后就能得到随机背景的均值和标准差再用观测值减去随机均值、除以标准差得到标准化后的得分。这个得分比原始互信息更有可比性因为它已经减掉了背景偏差。4.2 返回结构里有什么ak_covar返回的res是一个 MATLAB struct里面主要字段如下字段含义res.score标准化后的位点对得分矩阵对称矩阵行列数等于输入位点数res.zscore对应的z分数矩阵用于评估偏离随机背景的程度res.pvalue置换检验得到的经验p值矩阵res.positions这些位点对应原始MSA中的列号等价于之前info.positionsres.method本次运行使用的算法名第一次跑完建议先看size(res.score)确认矩阵维度等于预处理后的列数而不是原始MSA列数。如果维度对不上多半前面某一步没有把矩阵传对。4.3 快速看最强的位点对得分矩阵是对称矩阵第i行第j列和j行i列是同一个配对因此只看上三角就够。下面这段代码可以打印出得分最高的前10个位点对score res.score; % 只取上三角排除对角线 [I, J, V] find(triu(score, 1)); % 按得分降序排列 [Vs, idx] sort(V, descend); % 打印前10个候选 for k 1:min(10, numel(idx)) rowIdx I(idx(k)); colIdx J(idx(k)); fprintf(候选配对 bit %d - %d : score %.3f z %.2f p %.3f\n, ... rowIdx, colIdx, Vs(k), ... res.zscore(rowIdx, colIdx), ... res.pvalue(rowIdx, colIdx)); end这段代码里的triu(score, 1)会保留矩阵对角线以上元素第二个参数1表示从对角线右侧开始取直接去掉对角线。find返回的三个输出分别是行列下标和非零值。排序后打印时行的索引和列的索引要一起从I(idx(k))和J(idx(k))取避免排序后错位。看到的结果里如果某个配对的p值特别低但z分数却不高要留意是否来自低频残基组合后面进阶部分会提到怎么处理。4.4 参数组合怎么调实际分析很难一次到位的下面是我通常使用的参数组合参数推荐值说明method先mi后scaMI做快速初筛SCA做模块验证normalizezscore减掉随机背景比原始分数可解释nperm1000追求统计稳定性时使用nperm100调参阶段先用快速模式有一点要特别注意切换method后不同算法的得分量纲完全不一样SCA的权重甚至会出现负值不能拿MI的阈值去直接套SCA结果。调参阶段我会先把nperm降到100甚至0这里nperm0表示只算原始分数不做置换让流程快速跑通确认矩阵尺寸和输出字段没问题后再把置换次数提上去。5. 避坑与常见问题让我返工的五个具体失误5.1 没有去冗余结果看起来“太完美”现象跑完IM后排名靠前的位点对几乎都落在一个非常保守的结构域上并且这些位点本身单独看也很保守和前几个位点的互信息高得像假数据一样。原因MSA里如果包含大量相似度超过95%的近似序列它们在谱系上是重复抽样会让少数几个常见残基组合反复出现互信息被这些重复样本刷高最终主导全局排名。协同进化分析真正需要的是多样化的独立序列而不是同一谱系的一堆副本。解决回到ak_filter_alignment把min_identity从0.9收紧到0.8必要时按构建好的进化谱系进行加权抽样。我一般会先打印去冗余前的序列数size(seqMat,1)和去冗余后size(alnClean,1)如果发现只去掉了不到10%的序列说明数据本身多样性不足即使算出来结果也最好不要直接下强结论。5.2 把 gap 当成一个普通字符现象得分矩阵里有一整块连续列在结果中排名特别高比如位点120到150。单独看它们的原始互信息发现这些列在两两之间几乎都是最高分但和功能没有任何关系。原因这些列在MSA里是一大片gap区许多序列在这个区域整体缺失gap字符因此成对共现。在互信息计算眼里这等价于“有大量序列在同一个位置缺失”于是统计量表现出强关联但这根本不是进化意义上的协同变异。解决用max_gap_frac先把高gap列过滤掉。如果某个区域没有完全过滤干净可以额外保留一个gap掩码矩阵在调用ak_covar时把掩码传入至少不要让-参与最后的字符频次统计。我常用的做法是先把max_gap_frac定在0.3再人工看info.positions是否还有连续成片被保留的gap区域。5.3 过滤后的位点编号和结构编号对不上现象结果显示第42位和第88位是显著共进化残基对但把这两个编号直接对应到三维结构时发现结构里该区域根本不存在这两个残基或者对应到了完全不同位置的残基上。原因MSA本身不是从结构编号开始的。原始比对里可能删掉了一些空位而ak_filter_alignment又基于旧的MSA矩阵删除了高gap列所有后续得到的编号都是“当前矩阵的第几列”而不是原始序列上的残基编号更不是PDB文件里的编号。解决从一开始就保存info.positions每个分析结果里的位点编号都要通过它做一次映射把“当前矩阵列号”转成“原始MSA列号”再结合你读入结构文件时建立的序列到结构编号映射做第二次转换。这个映射我每次都会写进脚本第一段而不是等结果出来后再手工拼。5.4 置换检验时间爆炸现象当位点数到200以上、nperm设为1000时作业跑了几个小时还没结束。中间杀掉进程后检查发现每个位点对都要重新扫一遍矩阵计算量直接膨胀。原因nperm是对每一位点对做置换重算位点数从100增加到200时两两配对数增加到4倍如果再叠加1000次置换总体耗时就是非线性上升。AKtoolbox默认在置换时会重建整棵统计量矩阵保证p值的一致性但代价就是慢。解决调参阶段先用nperm100甚至nperm0。如果最终报告里只关心排名前几十的候选位点对可以只对候选对单独做置换检验这样计算量减少一个量级。另一个做法是使用块置换而不是逐字符置换保持序列内部相关性随机性更贴近真实进化过程。5.5 小样本硬套 SCA还拿谱去下结论现象只有20条序列跑SCA后得到的第一特征向量看起来有明确分层第一特征值解释了很大方差但换用MI方法后这些“强信号位点”大多对不上。原因SCA在构造协方差矩阵时需要足够的独立样本去估计残基之间的协变模式20条序列对L×L维矩阵来说自由度严重不足。少数几条每序列自带的特殊突变组合会被当成稳定特征放大看起来像模块化信号实际上是抽样噪声。解决少于50条序列的MSA不要做SCA先用MI做初筛序列数足够时再跑SCA并且对比两种方法得到的强候选列表是否一致。如果两个方法都指向同一组位点对结论才相对可信。小样本下看到SCA谱再有“规律”也不要直接写进报告除非用置换检验验证过。6. 进阶从得分矩阵到残基配对网络以及一个收尾习惯6.1 给得分矩阵做保守性屏蔽协同进化分数高的位点对不一定都有生物学意义尤其是互信息方法它会把保守位点之间的微弱共变放得很大。一个实用的技巧是先计算每个位点的保守程度把过于保守的位点排除掉只保留有一定变异度、又表现出联动信号的位点。% 计算每个位点的香农熵-1表示完全保守0表示完全均匀 H ak_entropy(alnClean); conserv 1 - H; % 设置一个保守性阈值太保守的位点不进入最终网络 keep conserv 0.2; scoreClean res.score .* (keep keep);这里的ak_entropy返回每个MSA位点的归一化香农熵值conserv越接近1代表位点越保守。keep是一个逻辑向量keep keep会把它扩展成一个逻辑矩阵只有两个位点都足够可变时对应位置上的得分才保留下来其他位置全部置零。这个操作能明显减少保守位点主导的问题让网络更聚焦在有真实进化约束的残基上。6.2 从矩阵到边列表拿到清理后的得分矩阵下一步是把矩阵转成网络分析工具能读的边列表。我一般输出一个三列TSV文件分别是位点A、位点B、得分。% 取上三角转成网络边表 thr prctile(scoreClean(scoreClean eps), 99); [I, J, W] find(scoreClean thr); % 写出边表三列分别表示起始位点、目标位点、权重 edges [I, J, W]; writematrix(edges, covar_edges.tsv, Delimiter, tab);阈值我通常设为得分矩阵非零部分的高分位点第99百分位数优先保留最强信号。scoreClean thr会生成一个逻辑稀疏矩阵find提取出的三列就是边表。你可以在网络可视化软件里用这个边表画残基相互作用网络也可以把它和三维结构中的残基距离做对比看那些显著共进的残基对是否在空间上确实接近。最后说一个让我改变工作习惯的教训以前我拿到新MSA就直接跑协同进化结果总要先浪费一两轮在检查数据质量和编号映射上。从那以后我每次跑AKtoolbox之前都强制自己走一遍“序列数、冗余度、gap比例、编号映射”四项检查确认无误再分析。这个习惯帮我省掉了很多返工时间希望也能帮到你。本文还有配套的精品资源点击获取