

实战（二）：利用本地Alphafold3预测蛋白质和金属离子相互作用

Dr.X的基因空间

2024-12-02

导读：利用Alphafold3在本地完成蛋白质和金属离子的相互作用预测

实战（二）：利用本地Alphafold3预测蛋白质和金属离子相互作用

写在前面的
上一期推送初步介绍了如何利用本地搭建的Alphafold3完成蛋白质结构的从头预测。由于Alphafold3的特色并不局限于蛋白质结构的预测，它还有其独特的相互作用预测能力。因为，在DeepMind发布的Alphafold3 github网页中并未给出如何预测蛋白质和其他分子的相互作用的教材。所以，本期的推送将不依赖其官网教程，而是详细给出如何自定义选择蛋白和其他生物分子，并利用Alphafold3在本地完成蛋白-生物分子相互作用的分析。

利用Alphafold3预测钙离子结合蛋白与钙离子结合的相互作用

本次推送我将利用实验室常用的一个具有钙离子结合活性的蛋白和钙离子为示例，介绍如何使用alphafold3预测蛋白质与钙离子的相互作用。这个蛋白具有钙离子结合结构域（Calcium-binding domain），该结构域通常具有特定的三维结构，能够与钙离子形成稳定的配位，从而调节蛋白质的功能。以前分析这个蛋白的相关特性时，是通过序列比对的方式先在PDB数据库中确定与该蛋白具有相似序列的且有已知的和钙离子相互作用的三维结构的其他蛋白质，再通过同源建模的方式得到自己的蛋白和钙离子如何相互作用的信息。这次，在不依赖同源建模的情况下，利用如下文所述的方法即可用Alphafold3一步到位直接从头完成蛋白结构+金属离子相互作用结构预测。由于是自定义蛋白序列，下文的方法可以迁移拓展到其他各种蛋白的预测上。

数据准备

1.蛋白质序列准备

按照如下格式准备好待研究蛋白的fasta格式文件

>C1
MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH

将上面的蛋白质序列文件保存命名为proteins.faa。

2.编写生成用于Alphafold3预测蛋白和分子相互作用的json格式的批量生成程序

#!/usr/bin/perl
use strict;
use warnings;
use JSON;
use Bio::SeqIO;

my $fasta = Bio::SeqIO -> new(-file => $ARGV[0], -format => 'fasta');
while(my $seq = $fasta -> next_seq){
    my ($seq_id, $sequence) = ($seq -> id, $seq -> seq);
    # 创建数据结构
    my $data = [
        {
            "name" => $seq_id,
            "modelSeeds" => [42],
            "sequences" => [
                {
                    "proteinChain" => {
                        "sequence" => $sequence,
                        "count" => 1
                    }
                },
                {
                    "ion" => {
                        "ion" => $ARGV[1],
                        "count" => 1
                    }
                }
            ]
        }
    ];
    # 创建 JSON 对象
    my $json = JSON->new->utf8->pretty->canonical(1); # 打开 canonical 模式
    # 将数据结构编码为 JSON 字符串
    my $json_text = $json->encode($data);
    open(OUT, ">", $seq_id."_AF3_input.json");
    print OUT $json_text;
    close OUT;
}

将上面的代码保存为make-alphafold3-json.pl程序文件。该程序可以批量处理多个蛋白质序列。

3.制作蛋白质-金属离子相互作用输入文件

make-alphafold3-json.pl程序有两个输入参数，第一个参数输入蛋白质的fasta格式文件，第二个参数输入金属离子名称，如钙离子输入为CA。程序具体用法如下：

perl make-alphafold3-json.pl proteins.faa CA

完成运行后会在工作目录下产生一个后缀名为_AF3_input.json的文件，该文件可以用于本地Alphafold3进行蛋白-金属离子相互作用预测。本推送中生成的输入数据如下：

[
   {
      "modelSeeds" : [
         42
      ],
      "name" : "C1",
      "sequences" : [
         {
            "proteinChain" : {
               "count" : 1,
               "sequence" : "MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH"
            }
         },
         {
            "ion" : {
               "count" : 1,
               "ion" : "CA"
            }
         }
      ]
   }
]

使用本地Alphafold3完成蛋白-金属离子相互作用预测

1.运行代码

使用nohup命令通过docker在后台运行alphafold3

nohup \
docker run \
#挂载计算机中输入数据路径为docker空间中的虚拟路径/root/af_input
--volume --volume ~/home/software-test/AF3-test/af_input:/root/af_input \
#挂载计算机中~/home/software-test/AF3-test/af_output路径到docker虚拟空间，作为结果输出路径
--volume ~/home/software-test/AF3-test/af_output:/root/af_output \
#挂载af3.bin模型参数所在路径~/home/software/alphafold3/models为docker空间模型参数路径/root/models
--volume ~/home/software/alphafold3/models:/root/models \
#挂载AF3数据库路径~/home/software/alphafold3/database为docker空间数据路径/root/public_databases
--volume ~/home/software/alphafold3/database:/root/public_databases \
#使用2个GPU完成预测结构
--gpus 2 \
#调用docker空间中的alphafold3
alphafold3 \
python \
/app/alphafold/run_alphafold.py \
#指定输入蛋白质序列数据C1_AF3_input.json
--json_path=/root/af_input/C1_AF3_input.json \
--model_dir=/root/models \
--output_dir=/root/af_output &

2.运行记录

由于这次输入的蛋白质序列较短，总共267个氨基酸。整个预测过程不到40分钟就完成了。下面是预测过程中输出的记录。

I1201 06:24:56.178634 140130585034752 folding_input.py:1027] Detected /root/af_input/alphafold_input.json is an AlphaFold Server JSON since the top-level is a list.
I1201 06:24:56.178923 140130585034752 folding_input.py:1033] Loading 1 fold jobs from /root/af_input/alphafold_input.json
I1201 06:24:56.812045 140130585034752 xla_bridge.py:895] Unable to initialize backend 'rocm': module 'jaxlib.xla_extension' has no attribute 'GpuAllocatorConfig'
I1201 06:24:56.814730 140130585034752 xla_bridge.py:895] Unable to initialize backend 'tpu': INTERNAL: Failed to open libtpu.so: libtpu.so: cannot open shared object file: No such file or directory
I1201 06:25:09.705902 140130585034752 pipeline.py:40] Getting protein MSAs for sequence MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:25:09.707513 139988016559680 jackhmmer.py:78] Query sequence: MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:25:09.708212 139987647460928 jackhmmer.py:78] Query sequence: MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:25:09.708981 139988016559680 subprocess_utils.py:68] Launching subprocess "/hmmer/bin/jackhmmer -o /dev/null -A /tmp/tmp1hukqmu0/output.sto --noali --F1 0.0005 --F2 5e-05 --F3 5e-07 --cpu 8 -N 1 -E 0.0001 --incE 0.0001 /tmp/tmp1hukqmu0/query.fasta /root/public_databases/uniref90_2022_05.fa"
I1201 06:25:09.709141 139987144144448 jackhmmer.py:78] Query sequence: MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:25:09.709293 139987647460928 subprocess_utils.py:68] Launching subprocess "/hmmer/bin/jackhmmer -o /dev/null -A /tmp/tmpew28852j/output.sto --noali --F1 0.0005 --F2 5e-05 --F3 5e-07 --cpu 8 -N 1 -E 0.0001 --incE 0.0001 /tmp/tmpew28852j/query.fasta /root/public_databases/mgy_clusters_2022_05.fa"
I1201 06:25:09.709605 139986942817856 jackhmmer.py:78] Query sequence: MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:25:09.901079 139987144144448 subprocess_utils.py:68] Launching subprocess "/hmmer/bin/jackhmmer -o /dev/null -A /tmp/tmp89sqme3k/output.sto --noali --F1 0.0005 --F2 5e-05 --F3 5e-07 --cpu 8 -N 1 -E 0.0001 --incE 0.0001 /tmp/tmp89sqme3k/query.fasta /root/public_databases/bfd-first_non_consensus_sequences.fasta"
I1201 06:25:09.901253 139986942817856 subprocess_utils.py:68] Launching subprocess "/hmmer/bin/jackhmmer -o /dev/null -A /tmp/tmphf0fmn30/output.sto --noali --F1 0.0005 --F2 5e-05 --F3 5e-07 --cpu 8 -N 1 -E 0.0001 --incE 0.0001 /tmp/tmphf0fmn30/query.fasta /root/public_databases/uniprot_all_2021_04.fa"
I1201 06:29:34.884322 139987144144448 subprocess_utils.py:97] Finished Jackhmmer in 264.983 seconds
I1201 06:38:24.718533 139988016559680 subprocess_utils.py:97] Finished Jackhmmer in 795.009 seconds
I1201 06:44:49.236162 139986942817856 subprocess_utils.py:97] Finished Jackhmmer in 1179.335 seconds
I1201 06:46:08.973908 139987647460928 subprocess_utils.py:97] Finished Jackhmmer in 1259.265 seconds
I1201 06:46:08.992091 140130585034752 pipeline.py:73] Getting protein MSAs took 1259.29 seconds for sequence MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:46:08.992226 140130585034752 pipeline.py:79] Deduplicating MSAs and getting protein templates for sequence MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:46:09.010417 139988016559680 subprocess_utils.py:68] Launching subprocess "/hmmer/bin/hmmbuild --informat stockholm --hand --amino /tmp/tmpzjr7_zgz/output.hmm /tmp/tmpzjr7_zgz/query.msa"
I1201 06:46:09.166867 139988016559680 subprocess_utils.py:97] Finished Hmmbuild in 0.156 seconds
I1201 06:46:09.167667 139988016559680 subprocess_utils.py:68] Launching subprocess "/hmmer/bin/hmmsearch --noali --cpu 8 --F1 0.1 --F2 0.1 --F3 0.1 -E 100 --incE 100 --domE 100 --incdomE 100 -A /tmp/tmpotzhrp5a/output.sto /tmp/tmpotzhrp5a/query.hmm /root/public_databases/pdb_seqres_2022_09_28.fasta"
I1201 06:46:15.128245 139988016559680 subprocess_utils.py:97] Finished Hmmsearch in 5.960 seconds
I1201 06:46:40.539313 140130585034752 pipeline.py:108] Deduplicating MSAs and getting protein templates took 31.55 seconds for sequence MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:46:40.539500 140130585034752 pipeline.py:115] Filtering protein templates for sequence MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:46:40.540956 140130585034752 pipeline.py:124] Filtering protein templates took 0.00 seconds for sequence MGMALPKTGKPTAKQVVDWAINLIGSGVDVDGYYGRQCWDLPNYIFNRYWNFKTPGNARDMAWYRYPEGFKVFRNTSDFVPKPGDIAVWTGGNYNWNTWGHTGIVVGPSTKSYFYSVDQNWNNSNSYVGSPAAKIKHSYFGVTHFVRPAYKAEPKPTPPADKITWNWKGVFYPNPEKAIRVRKTAGLTGTVVEEDSWLYTKDDWVKFDQVIKKDGYWWIRFKYQREGSSTNNFYCAVCRITDKEQKIKNEKYWGTIEWALEHHHHHH
I1201 06:46:45.414394 140130585034752 pipeline.py:165] processing C1, random_seed=42
I1201 06:46:45.440038 140130585034752 pipeline.py:258] Calculating bucket size for input with 268 tokens.
I1201 06:46:45.440350 140130585034752 pipeline.py:264] Got bucket size 512 for input with 268 tokens, resulting in 244 padded tokens.
[06:46:47] UFFTYPER: Unrecognized atom type: Ca+2 (0)
[06:46:48] UFFTYPER: Unrecognized atom type: Ca+2 (0)
/alphafold3_venv/lib/python3.11/site-packages/alphafold3/model/confidences.py:337: RuntimeWarning: invalid value encountered in divide
  return np.nanmean(value * mask_with_nan, axis=axis) / np.nanmean(
Running AlphaFold 3. Please note that standard AlphaFold 3 model parameters are
only available under terms of use provided at
https://github.com/google-deepmind/alphafold3/blob/main/WEIGHTS_TERMS_OF_USE.md.
If you do not agree to these terms and are using AlphaFold 3 derived model
parameters, cancel execution of AlphaFold 3 inference with CTRL-C, and do not
use the model parameters.
Found local devices: [CudaDevice(id=0), CudaDevice(id=1)]
Building model from scratch...
Processing 1 fold inputs.
Processing fold input C1
Checking we can load the model parameters...
Running data pipeline...
Processing chain A
Processing chain A took 1290.93 seconds
Processing chain B
Processing chain B took 0.00 seconds
Output directory: /root/af_output/c1
Writing model input JSON to /root/af_output/c1
Predicting 3D structure for C1 for seed(s) (42,)...
Featurising data for seeds (42,)...
Featurising C1 with rng_seed 42.
Featurising C1 with rng_seed 42 took 3.59 seconds.
Featurising data for seeds (42,) took  8.34 seconds.
Running model inference for seed 42...
Running model inference for seed 42 took  62.31 seconds.
Extracting output structures (one per sample) for seed 42...
Extracting output structures (one per sample) for seed 42 took  0.32 seconds.
Running model inference and extracting output structures for seed 42 took  62.63 seconds.
Running model inference and extracting output structures for seeds (42,) took  62.63 seconds.
Writing outputs for C1 for seed(s) (42,)...
Done processing fold input C1.
Done processing 1 fold inputs.

使用pymol查看蛋白-金属离子相互作用

Alphafold3输出了5个预测结果，根据ranking score评估，第四个模型是最优结果。

seed,sample,ranking_score
42,0,0.9080788902069712
42,1,0.907647363594914
42,2,0.909683814307021
42,3,0.9097932768353593
42,4,0.9110906368967179

因此下载第四个模型对应的cif文件在本地使用pymol可视化，如下图所示，红色部分标记的氨基酸则是预测出来与钙离子（图中绿色球体）具有相互作用的氨基酸。

对比在线Alphafold3网页预测结果

打开Alphafold3在线官网，将相同的序列输入至官网预测，并下载预测结果。使用pymol比较本地和在线的差异，下图中蓝色的结构是在线预测的结构，绿色的结构是本地预测结果，红色标记是在线和本地预测出来的钙离子。可以看出整体空间结构仍然具有一些差异，总体RMSD在1.9左右，但是无论是在线还是本地均将钙离子预测在了正确的结合位点。我检查了一下在线和本地预测结构有差异的原因，这是因为1.在线预测和本地预测的随机数种子不同，在线预测的随机数种子是228491813，而本地预测的随机数种子是42，我将本地预测的随机数种子调整为和在线预测的一致后，这种差异会显著缩小，2.蛋白质中存在一些无规则卷曲，这些无规则卷曲本身预测的鲁棒性比阿尔法螺旋和贝塔折叠等低，因此尽管使用相同随机数种子，也有可能在无规则卷曲处存在细微差异。建议一次预测中使用多个不同的随机数种子评估一下预测结果的可靠性。