AiCE 是一款基于蛋白质逆折叠（inverse folding）模型的高适配度突变（high-fitness mutation）预测工具。它把 ProteinMPNN 等 AI 逆折叠模型输出的结构兼容序列，结合结构约束与进化约束，筛选出能优化蛋白质功能的候选突变。本文带你只读scripts/目录下的脚本，从零读懂这套完整预测流水线的实现逻辑，帮你快速掌握单点突变提名、LD 矩阵、SCA 矩阵与组合突变评分四大环节是如何串起来的。

📁 scripts 目录全景：先看目录怎么组织

scripts/是整个 AiCE 的执行中枢，采用"编号主脚本 + Python 原子功能 + 第三方子模块"的三层结构。理解这个组织方式，是读懂流水线的前提：

• 编号主脚本：01.single_mut_prediction.sh→04.com_mut_prediction.sh，对应流水线的四个大步骤，负责"编排调度"；
• 原子 Python 功能：predict_dssp.py、count_residue_freq.py、predicted_single_HF_mutations.py、02.caculated_ld.py、com_mut_prediction.py等，每个只干一件事；
• 第三方子模块：scripts/ProteinMPNN/ 与 scripts/pySCA/，分别承担逆折叠建模和进化耦合分析。

💡 关键设计：每个主脚本都用**文件前缀（prefix）**作为"数据总线"。比如结构文件SpCas9_af3.pdb，会依次生成SpCas9_af3.fa、.ss、.freq、.ld、.sca_matrix.tsv、.comb等同前缀文件，后续步骤靠前缀自动"认领"输入。

输入 SpCas9_af3.pdb │ 01 单点突变提名 ├─→ .fa .ss .freq .re.freq .comb .mut │ ├─ 02 LD矩阵 ──→ .ld .vcf │ ├─ 03 SCA矩阵 ─→ .sca_matrix.tsv .db │ └─ 04 组合突变 ─→ .sca.result .ld.result

🧬 第一步：单点突变提名（01.single_mut_prediction.sh）

scripts/01.single_mut_prediction.sh 是流水线的"入口"，它把逆折叠、二级结构、频率统计、突变预测四件事打包成一条命令。调用链是：

01.single_mut_prediction.sh └─ inverse_MPNN.sh （跑 ProteinMPNN 逆折叠，产出 .fa 序列） └─ predict_dssp.py （DSSP 二级结构 → .ss） └─ count_residue_freq.py （残基频率统计 → .freq / .re.freq） └─ predicted_single_HF_mutations.py （筛选高适配度单点突变 → .comb / .mut）

具体分工如下：

• 逆折叠：scripts/inverse_MPNN.sh 先解析 PDB（parse_multiple_chains.py），再调用 scripts/ProteinMPNN/protein_mpnn_run.py，用--num_seq_per_target 1000、--sampling_temp 0.5采样出上千条结构兼容序列，写成 FASTA（.fa）。
• 二级结构：scripts/predict_dssp.py 调用mkdssp对每个残基标注二级结构，输出.ss。
• 频率统计：scripts/count_residue_freq.py 统计每个位点上各氨基酸的出现频率，找出"高频氨基酸"，输出.freq与重排后的.re.freq。
• 突变提名：scripts/predicted_single_HF_mutations.py 把"频率"和"二级结构"两路信息合并，用beta（全局发生阈值，推荐 0.8）和gamma（柔性区阈值，推荐 0.5）两个参数做"AiCE 过滤"，筛出高适配度（HF）单点突变，写出.comb与.mut。

一句话理解这一步：用 AI 生成"结构允许的所有序列"，再用统计+结构规则挑出最可能提升适应度的单点突变。

📊 第二步：构建 LD 关联矩阵（02.caculated_ld.py）

scripts/02.caculated_ld.py 负责计算连锁不平衡（Linkage Disequilibrium, LD）矩阵，用来刻画"不同位点变异之间在序列空间里的共现依赖"。它有一个巧妙的桥接思路——把蛋白质序列"伪反转译"成 DNA，再借用群体遗传学工具 plink 计算：

1. 蛋白 → DNA：用内置的最优密码子表OPTIMAL_CODONS把每条.fa序列反转译成 DNA FASTA；
2. DNA → VCF：以第一条序列为参考，逐位点生成分型（genotype），写出标准.vcf；
3. plink 算 LD：调用../scripts/plink/plink --r2 square得到 r² 关联矩阵，输出.ld与.vcf；
4. 后处理：用sed把制表符换成逗号，并清理中间文件。

⚠️ 注意：LD 环节默认绑定 plinkv1.9.0-b.7.7，plink v2.0 需自行适配。这一步的产物.ld、.vcf会交给第四步做组合评分。

🧩 第三步：构建 SCA 进化耦合矩阵（03.caculated_sca.sh + pySCA）

SCA（统计耦合分析）刻画的是位点之间的进化耦合，是 AiCE 引入"进化约束"的核心。scripts/03.caculated_sca.sh 是一个薄封装，它遍历每个.fa，依次调用scripts/pySCA/下的两个核心模块：

• scaProcessMSA.py：预处理多序列比对（MSA），生成中间.db数据库文件；
• scaCore.py：执行 SCA 核心计算（SVD 分解、sector 划分、耦合矩阵），输出.sca_matrix.tsv。

pySCA 本身源自 scripts/pySCA/ 的开源实现，AiCE 作者针对本项目做了改造（见 scripts/pySCA/README）。

SCA 的核心思想是：蛋白质家族里协同进化的位点往往在功能上相互依赖。下图展示了 SCA 把 DHFR 蛋白的整体"扇区（sector）"进一步分解到 IC1–IC6 各独立耦合模式的立体视图，红色区域就是耦合强度高的残基簇：

SCA 还会把耦合模式组织成层级化的"扇区（sector）"，用于识别功能模块。下面的层级树说明了 β-内酰胺酶中不同 IC（独立耦合）模式如何归并成 sector 1 / sector 2：

在教程数据里，SCA 还会输出各独立分量（V、U）的散点分布，用来直观检查耦合的稳定性——这也是理解.sca_matrix.tsv内部结构的钥匙：

🎯 第四步：组合突变预测（04.com_mut_prediction.sh）

最后一步 scripts/04.com_mut_prediction.sh 把前序产物"合体"，对组合突变打分。它会为每个.fa前缀同时跑两次 scripts/com_mut_prediction.py：

• 一次用.sca_matrix.tsv（进化耦合视角）；
• 一次用.ld+--vcf（序列共现视角）。

两种矩阵各自给出：平均两两得分、对数得分，以及一个0/1 推荐标志（默认阈值：LD 分 0.5、SCA 分位 0.9），最终写入.sca.result和.ld.result。参数很灵活——

• 传数字（如2）：枚举该阶数的所有突变组合；
• 传列表文件：只对你指定的位点组合打分。

🧰 辅助脚本与工具一览

🔗 串起来：用"前缀数据流"读懂整条流水线

读 AiCE 的scripts/，最省力的方式是记住一条主线——同一个文件前缀在四步里像接力棒一样传递数据：

1. 01把 PDB 变成序列 + 结构 + 频率三件套（.fa/.ss/.freq）并提名单点突变（.comb/.mut）；
2. 02借"蛋白→DNA→VCF→plink"算出序列共现矩阵（.ld）；
3. 03借 pySCA 算出进化耦合矩阵（.sca_matrix.tsv）；
4. 04用两张矩阵给组合突变交叉打分，给出是否推荐的 0/1 标志。

结构约束来自第 1 步的 DSSP + ProteinMPNN，进化约束来自第 2、3 步的 LD 与 SCA，四步叠加正是论文"结构约束 + 进化约束"方法论在代码里的落地。理解了这条前缀数据流，你就已经掌握了整个 AiCE 高适配度突变预测流水线的实现骨架。

📌 官方演示见 AiCE_demo.ipynb，完整用法与参数说明见 README.md。