使用 ESM 计算蛋白序列的突变影响分数

ESM（Evolutionary Scale Modeling） 是由 Facebook AI Research 开发的蛋白质语言模型。其基于大规模的蛋白质序列数据进行训练，能够捕捉蛋白质序列中的复杂模式和进化信息，从而在多个生物学下游应用中表现出色。

预训练好的 ESM 可用于多个生物学下游应用，本文将就其突变影响分析方面的功能进行介绍。

关于 ESM

ESM 预测突变影响的方式是无监督的，也就是说 ESM 的训练过程中并未使用到与突变影响相关的标签进行有监督学习。其评估每个氨基酸错义突变影响程度的方式是掩蔽该位置上的氨基酸并计算两种氨基酸（参考氨基酸与突变氨基酸）在该位置出现的对数似然比：

ESM 的预训练阶段使用了类似于 BERT 的掩码任务，因此 ESM 可以给出每个位点里各类氨基酸出现的概率，而参考氨基酸（Ref，即未突变时的氨基酸）和突变氨基酸（Alt）的概率差异即可作为一种 effect score。可以简单理解为，突变氨基酸出现的概率比参考氨基酸低越多，说明这类突变越有害。这种无监督的评估方式也被用于许多其他的深度学习工具，例如 EVE 和 GPN 等。

Reference: Brandes et al. 2023, Nat Genet

对于多氨基酸突变，ESM 采取的策略是简单的加性模型，即对每个位置计算出其对数似然比后，进行加和得到最终评分（Meier et al. 2021, bioRxiv）。其他研究者的 benchmark 测试表明，ESM 在诸多数据集上具有出色性能（Tabet et al. 2024, Genome Biol）。

ESM 使用

在使用下述任何一种方法前需提前安装好需要的 module：

1	pip install tqdm numpy pandas biopython torch fair-esm

其他需要的 module 可在后续使用 script 报错时自行通过 pip 或 conda 补充安装。

原始的 ESM 工作流程

使用 git clone 下载其 github repository：

1	git clone https://github.com/facebookresearch/esm.git

如果你的设备无法连接外网，你也可以尝试克隆 gitee 源：

1	git clone https://gitee.com/li42110/esm.git

在 git clone 下载完成后，运行以下命令跳转到为 variant prediction 提供的示例文件夹中：

1	cd esm/examples/variant-prediction

该文件夹下的 predict.py 即 ESM 团队提供的用于计算 effect score 的 script，其中示例命令如下：

python predict.py \
    --model-location esm1v_t33_650M_UR90S_1 esm1v_t33_650M_UR90S_2 esm1v_t33_650M_UR90S_3 esm1v_t33_650M_UR90S_4 esm1v_t33_650M_UR90S_5 \
    --sequence HPETLVKVKDAEDQLGARVGYIELDLNSGKILESFRPEERFPMMSTFKVLLCGAVLSRVDAGQEQLGRRIHYSQNDLVEYSPVTEKHLTDGMTVRELCSAAITMSDNTAANLLLTTIGGPKELTAFLHNMGDHVTRLDRWEPELNEAIPNDERDTTMPAAMATTLRKLLTGELLTLASRQQLIDWMEADKVAGPLLRSALPAGWFIADKSGAGERGSRGIIAALGPDGKPSRIVVIYTTGSQATMDERNRQIAEIGASLIKHW \
    --dms-input ./data/BLAT_ECOLX_Ranganathan2015.csv \
    --mutation-col mutant \
    --dms-output ./data/BLAT_ECOLX_Ranganathan2015_labeled.csv \
    --offset-idx 24 \
    --scoring-strategy wt-marginals

--model-location：选择模型，此处选择的是 5 个 ESM-1v 模型（这些模型使用不同的随机种子生成训练数据集并进行训练），也可选择其他模型例如 ESM-1b (esm1b_t33_650M_UR50S) 或 ESM2 (esm2_t33_650M_UR50D) 等。
--sequence：输入的蛋白序列，注意过长的序列可能会引发内存溢出或者其他意外状况（由于 ESM 的 token 最大数量限制），如果你输入的序列过长，请直接跳转后文 优化的 ESM 工作框架 部分。
--dms-input：输入的表格，其中需包含 mutant 列（--mutation-col 中指定的列名），mutant 格式为 [ref amino][idx][alt amino]，例如示例文件中的 H24C 表示的为第 24 位上的 H 突变为 C。请在该表格的该列中填入所有你希望计算 effect score 的 mutant。
--mutation-col：指定突变在输入表格中所在的列名，此处即 mutant。

红框标注即 mutant 列，后续的列并不是必需的

--dms-output：输出的表格，其中最后 n 列（n 的具体数字取决于你输入的模型数量）为各模型的评分。
--offset-idx：突变位置的偏倚，例如在 ESM 示例中第一个 mutant 为 H24C，但 H 其实位于序列的开端，因此将此参数设置为 24 使其从 24 开始位置计数。如果你从 1 开始，则将该项设置为 1 即可。
--scoring-strategy：计算 effect score 的策略，分为三种。wt-marginals 直接将 wt 序列（未突变序列）输入模型得到每个 token 的 embedding 并计算其出现各个氨基酸的概率；masked-marginals 将 wt 序列中每个 token 逐个掩码并计算掩码位置上出现各个氨基酸的概率；pseudo-ppl 将 mutant 序列中每个 token 逐个掩码并计算整条序列的伪困惑度作为其 effect score。

关于 --scoring-strategy 参数，masked-marginals 的计算方式是最符合直觉和逻辑的，但由于 masked-marginals 计算中要对每个 token 都掩码计算一次，因此时间开销很高（特别是在序列很长的情况下）。相比之下，wt-marginals 虽然向模型暴露了原始序列信息，但仅需计算一次便可得到所有位置上各个氨基酸的出现概率。已有研究分析表明 wt-marginals 也具有很好的性能表现（虽然其会高估原始氨基酸的概率），建议此处选择 wt-margins 即可。

运行以上命令后，如果你缺乏这些模型的文件，脚本会自行开始下载并存储到特定位置（例如 Linux 为 home 目录下的 .cache/torch/hub/checkpoints/）。该步骤同样可能遇到网络问题，此时你可以选择自行从其他途径下载到指定路径中。如果要预测其他序列，请修改 --sequence 参数并自行准备 --dms-input 的输入文件。

下载的 ESM model

优化的 ESM 工作框架

使用 ESM github repository 中的 script 虽然可以达成目的，但是它有以下一些缺点：

不适用过长的蛋白序列。
定制化程度较高，迁移到自己的需求上时有一定的难度（如果数据量大），门槛相对更高。

基于此，这里推荐一个由 Brandes 等人开发的流程，该流程经过了以下改进：

只需输入 fasta 文件即可批量计算这些序列中所有氨基酸位点上各个突变的 effect score。
对于过长的序列，该流程使用滑动窗口方法将长序列分解成具有 overlap 的多条序列，并对重叠位置上的 effect score 进行加权平均计算。
对于序列中的插入或缺失同样可以进行计算（衡量标准是整条序列的似然比差异）。

使用 git clone 配置：

1 2	git clone https://github.com/ntranoslab/esm-variants.git cd esm-variants

该文件夹下的 esm_score_missense_mutations.py 使用示例：

1	python esm_score_missense_mutations.py --input-fasta-file [fasta] --output-csv-file [output csv] --model-name esm1b_t33_650M_UR50S

[fasta] 处填入你所需要分析的 fasta 序列文件路径，[output csv] 处填入输出的 csv 表格文件路径（含文件名）。
可通过 --model-name 指定使用的 ESM 模型，默认为 ESM-1b。

以下是一个运行示例，展示如何使用，首先建立一个 test.fa，填入以下内容：

>seq1
FISHWISHFQRCHIPSTHATARECRISP
>seq2
RAGEAGAINSTTHEMACHINE

运行脚本，指定输出文件为 test.csv：

1	python esm_score_missense_mutations.py --input-fasta-file test.fa --output-csv-file test.csv

运行结束后，test.csv 中的第一列为 fasta 文件中对应的序列 id，第二列为该序列所有可能的氨基酸突变，第三列为基于 ESM1b 计算得到的 ESM score：

关于 InDel 的 score 计算此处不进行详细阐述，如感兴趣可参考其 github repository README：

https://github.com/ntranoslab/esm-variants?tab=readme-ov-file#scoring-multi-residue-mutations

参考材料

Meier et al. 2021, bioRxiv: Language models enable zero-shot prediction of the effects of mutations on protein function
Tabet et al. 2024, Genome Biol: Benchmarking computational variant effect predictors by their ability to infer human traits
Brandes et al. 2023, Nat Genet: Genome-wide prediction of disease variant effects with a deep protein language model
EVE (Nature 2021): Disease variant prediction with deep generative models of evolutionary data
GPN (PNAS 2023): DNA language models are powerful predictors of genome-wide variant effects

题外话

esm-variant repository 中有一个 issue 针对其函数命名提出了质疑：

https://github.com/ntranoslab/esm-variants/issues/12

发起该 issue 的人是 ColabFold 的创建者，大致内容是 esm-variant 算出来的 score 不应该称作 pll (pseudo-likelihood)。基于这一点，本文也避开了 pll 这一说法直接使用 score 一词进行说明。各位朋友如果要使用其他关于 score 的称呼请自行斟酌。