计算解析并可视化Alphafold3结构预测结果

先放结论

本篇博客所讨论的结构数据处理、以及可视化内容，我按照自己的想法已经实现、整合成了1个工具，并做了一定的开发。

详细使用以及文档可以参考：

https://github.com/MaybeBio/AlphaFold3-SeqVisToolkit

欢迎PR和Issue！

背景

我有了alphafold3预测的结构数据，这个数据大多是PDB格式或者是mmCIF格式，也就是生物大分子晶体结构数据。

现在我们想要解析、挖掘这个结构预测文件中的空间坐标信息，希望能够得到一些有用的信息，

但是，怎么访问、处理、解析以及可视化这个文件？怎么整合一维的序列特征信息，整合到三维的结构可视化中？

二，走向3D：BioPython的PDB模块

Bio.PDB是Biopython中处理生物大分子晶体结构的模块。除了别的类之外，Bio.PDB包含PDBParser类，此类能够产生一个Structure对象，以一种较方便的方式获取文件中的原子数据。只是在处理PDB文件头所包含的信息时，该类有一定的局限性。

这里提一句题外话, 大家感兴趣的，可以用我之前博客中提到的、我开发的工具LibInspector去check一下Bio.PDB这个库。相关的文件我已经测试过，并放在仓库https://github.com/MaybeBio/LibInspector/blob/main/test/Bio_PDB.md中了。

另外注意：PDB 文件格式不再被修改或扩展以支持新内容，PDBx/mmCIF 在 2014 年成为标准的 PDB 存档格式。全球蛋白质数据库（wwPDB）的所有站点都使用大分子晶体学信息文件（mmCIF）数据词典来描述 PDB 条目的信息内容。mmCIF 使用灵活且可扩展的键值对格式来表示大分子结构数据，并且对单个 PDB 条目中可以表示的原子数、残基数或链数没有限制（不拆分条目！）。

所以我接下来的演示还是以mmCIF格式数据为主：

1，坐标文件读取

关于文件的读取访问部分，先上结论：

# 1,读取PDB文件: # 创建一个 PDBParser 对象 from Bio.PDB.PDBParser import PDBParser p = PDBParser(PERMISSIVE=1) # 接着通过 PDBParser 解析PDB文件, 就产生了Structure对象, 实际上就是面向对象编程那一套 structure_id = "1fat" filename = "pdb1fat.ent" s = p.get_structure(structure_id, filename) # 2,读取mmCIF文件: # 与PDB文件的情形类似，先创建一个 MMCIFParser 对象 from Bio.PDB.MMCIFParser import MMCIFParser parser = MMCIFParser() # 然后用这个解析器从mmCIF文件创建一个结构对象： structure = parser.get_structure('1fat', '1fat.cif')

这里我举一个我自己的实际例子：

我手头上有一个CTCF_noDNA.cif的结构预测文件，是取的Alphafold3预测结果5个文件中的model 0，因为官方声明model 0是置信度最高的。

现在我想访问并解析一下这个数据，但是手动解析太麻烦了，所以借助计算机：

from Bio.PDB import MMCIFParser parser = MMCIFParser() ctcf_noDNA = parser.get_structure('CTCF_noDNA', 'CTCF_noDNA.cif')

这里注意一下，初始化一个解析器的示例的时候，MMCIFParser()的参数参考官网：

https://biopython.org/docs/1.75/api/Bio.PDB.MMCIFParser.html

有人初始化喜欢warning都不看，当然默认是False。

如果是读取BinaryCIF文件，同理如下：

还有一些不常见的文件格式读取需求，可以参考：https://biopython.org/docs/latest/Tutorial/chapter_pdb.html

2，文件内容访问（主要是字典数据）

其中比较重要的一点是结构对象的header属性：这是一个将头记录映射到其相应值的Python字典。

有了键值对，我们就可以按照字典处理逻辑来访问这个cif文件中的具体结构数据了，

但是，我们一般操作的时候其实很少直接访问这个cif文件，因为键值对其实可以直接保存到字典数据中，减少直接访问cif文件的开销。

为了尽量少访问mmCIF文件，可以用MMCIF2Dict 类创建一个Python字典来将所有mmCIF文件中各种标签映射到其对应的值上。若有多个值（像 _atom_site.Cartn_y 标签，储存的是所有原子的y坐标值），则这个标签映射到一个值列表。从mmCIF文件创建字典如下：

from Bio.PDB.MMCIF2Dict import MMCIF2Dict ctcf_noDNA_dict = MMCIF2Dict('CTCF_noDNA.cif')

我们可以看到，这个字典结构其实会更加复杂：

访问就很简单，比如说我要获取包含所有原子y坐标的列表：

3，写入PDB/mmCIF文件

有了文件读取，那当然少不了文件输出了。

就好比我们想要在pymol或者是chimerax中导入这些结构数据，在进行数据处理+美化+select之后，我们有将部分数据导出的需要。

这个可以用PDBIO类实现。当然也可很方便地输出一个结构的特定部分。

比如说前面1-小节文件读取中我读进去的1个结构对象ctcf_noDNA

from Bio.PDB import PDBIO # 注意需要导入这个类 io = PDBIO() io.set_structure(ctcf_noDNA) io.save('ctcf_noDNA.pdb')

可以看到，这就是一个简单的pdb文件，我们已经拿到了一个PDB文件了，

其实这里也解决了一个问题，就是mmCIF文件转化到PDB文件，怎么操作，

就可以用这个！

当然，我们主要还是以写入mmCIF文件为主：

MMCIFIO 类可用于将结构写入 mmCIF 文件格式，

from Bio.PDB import MMCIFIO io = MMCIFIO() io.set_structure(ctcf_noDNA) io.save("1.cif")

现在我已经将我的CTCF_noDNA.cif文件格式转换为了1.cif文件：

Select 类可以与前面的 PDBIO 类以类似的方式使用。使用 MMCIF2Dict 读取的 mmCIF 字典也可以写入。

比如说我前面2-小节中构建的一个字典ctcf_noDNA_dict，我同样可以输出为mmcif格式

io = MMCIFIO() io.set_dict(ctcf_noDNA_dict) io.save("2.cif")

当然，上面只是简单的输出全部结构数据，实际操作过程中我们肯定是只想输出部分特殊区域或者是指定范围残基的结构数据，也就是写出结构的一部分。

可以用 Select 类（也在 PDBIO 中）来实现。 Select 有如下四种方法：

在默认情况下，每种方法的返回值都为1（表示model/chain/residue/atom被包含在输出结果中）。通过子类化 Select 和返回值0，我们可以从输出中排除model、chain等。

比如说下面的演示，我将只输出我这个结构数据中的Pro脯氨酸：

注意，这里是继承了基类Select

此处参考官方文档：https://biopython.org/docs/latest/Tutorial/chapter_pdb.html

但是注意官方文档中有些内容没有更新过来，也就是residue.get_name()。

现在更新之后residue已经没有get_name这个方法了，我的方法get_resname如下：

from Bio.PDB import Select class ProSelect(Select): def accept_residue(self, residue): if residue.get_resname()=='PRO': return True else: return False io = PDBIO() io.set_structure(ctcf_noDNA) io.save('pro_only.pdb', ProSelect())

现在我已经拿到了这个结构文件上的所有Pro信息：

官方文档还提供了其他方法：

4，结构表示

一个 Structure 对象的整体布局遵循称为SMCRA（Structure/Model/Chain/Residue/Atom，结构/模型/链/残基/原子）的体系架构：

这是许多结构生物学家/生物信息学家思考结构的方式，并提供了一种简单而有效处理结构的方法。额外的信息基本上是按需添加的。这种数据结构不一定最适合表示结构中的大分子内容，但它绝对必要，用于解释描述结构的数据文件（通常是 PDB 或 MMCIF 文件）中的数据。

• 如果这种层次结构无法表示结构文件的内容，那么可以相当肯定该文件包含错误或至少没有明确描述结构。
• 如果无法生成 SMCRA 数据结构，就有理由怀疑存在问题。
• 解析 PDB 文件可用于检测可能的问题。

Structure 对象的 UML 图，参考官网：

我们可以先忽视disorder部分，整体逻辑链就是SMCRA。

这个图其实可以从class的继承与引用角度来看，用来表示大分子结构的 Structure 类的SMCRA体系的UML图。带方块的实线表示集合，带箭头的实线表示引用，带三角形的实线表示继承，带三角形的虚线表示接口实现。

结构，模型，链，残基都是实体基类的子类。原子类仅仅（部分）实现了实体接口（因为原子类没有子类）。

对于每个实体子类，我们可以用该子类的一个唯一标识符作为键来提取子类（比如，可以用原子名称作为键从残基对象中提取一个原子对象；用链的标识符作为键从域对象中提取链）。

无序原子和残基用DisorderedAtom和DisorderedResidue类来表示，二者都是DisorderedEntityWrapper基类的子类。它们隐藏了无序的复杂性，表现得与原子和残基对象无二。

一般地，一个实体子类（即原子，残基，链，模型）能通过标识符作为键来从父类（分别为残基，链，模型，结构）中提取。

格式类似如下：

child_entity = parent_entity[child_id]

我们还可以从一个父实体对象获得所有子实体的列表。需要注意的是，这个列表以一种特定的方式排列（例如根据在模型对象中链对象的链标识符来排序）。

child_list = parent_entity.get_list()

我们也可以从子类得到父类

parent_entity = child_entity.get_parent()

在SMCRA的所有层次水平，我们还可以提取一个完整id 。

完整id是包含所有从顶层对象（结构）到当前对象的id的一个元组。

其实就是对每一个残基一个可以溯源的完整的身份证，比如说一个残基对象的完整id示例如下：

这个残基id表示该残基不是异质残基（也不是水分子），因为其异质值为空；而序列标识符为10，插入码为"A"。

要得到实体的id，用 get_id 方法即可：

entity.get_id()

可以用 has_id 方法来检查这个实体是否有子类具有给定id：

entity.has_id(entity_id)

实体的长度等于其子类的个数：

nr_children = len(entity)

对于从父实体得到的子实体，可以删除，重命名，添加等等，但这并不包含任何完整性检查（比如，有可能添加两个相同id的残基到同一条链上）。这就真的需要包含完整性检查的装饰类（Decorator）来完成了，但是如果我们想使用原始接口的话可以查看源代码（Entity.py，在官方biopython仓库中可以找到）。

下面我们来正式分析介绍SMCRA体系，

4.1 结构

结构对象是层次中的最高层。其id是用户指定的一个字符串。结构包含一系列子模型（S-M）。大部分晶体结构（但不是全部）含有一个单一模型，但是NMR结构通常由若干模型构成。晶体结构中大部分的无序也能导致多个模型（因为无序一般是按照构象集合的角度来研究的，也就是ensemble，所以一般都是集合形式）。

所以我们一般的代码中访问cif等坐标文件，都是从最高层structure开始的，

我们看类型也知道：

然后看子实体也就是model（因为S下一层就是M），我们想看子实体数目实际上就是看实体长度。

确实是一个，至少对于我们本篇博客分析的对象来说，也就是alphafold3预测的结构结果而言，就是1个结构中只含有1个单体model。

4.2 模型

结构域对象的id是一个整数，源自该模型在所解析文件中的位置（自动从0开始）。晶体结构通常只有一个模型（id为0），而NMR文件通常含有多个模型。然而许多PDB解析器都假定只有一个结构域， Bio.PDB 中的 Structure 类就设计成能轻松处理含有不止一个模型的PDB文件。

所以我们其实可以像从可迭代数据中获取切片数据一样，依据索引index获取一个实体的子实体数据。

比如说，从一个结构对象中获取其第一个模型

first_model = structure[0]

以我前面的结构对象为例，index单一不会报错，但是访问超过数目的下标就会报错。

4.3 链

模型对象存储着子链的列表。那就是同样的访问方式：
链对象的id来自PDB/mmCIF文件中的链标识符，是个单字符（通常是一个字母）。模型中的每个链都具有唯一的id。例如，从一个模型对象中取出标识符为“A”的链对象：

chain_A = model["A"]

对照前面的例子就是：ctcf_noDNA[0]是我们唯一的model，

然后我们可以从中获取各个链，比如说：

我们同样可以先看一下有多少条链：

4.4 残基

链对象储存着残基对象的列表。所以是同样访问逻辑，

比如说A链，727，其实就是CTCF蛋白质的全长了（也就是说这个例子中，A链其实就是蛋白质对象链）

但是要注意，残基存储的数据格式比较特殊，不是简单的list嵌套，其实顶层是嵌套，越往下越是底层实现的存储数据结构就会越复杂。

一个残基id是一个三元组：

采用该方案表示非标准field的原因见https://biopython.org/docs/latest/Tutorial/chapter_pdb.html#sec-hetero-problems

• 序列标识符 （resseq），一个描述该残基在链上的位置的整数（如100）；
• 插入码 （icode），一个字符串，如“A”。插入码有时用来保存某种特定的、想要的残基编号体制。一个Ser 80的插入突变（比如在Thr 80和Asn 81残基间插入）可能具有如下序列标识符和插入码：Thr 80 A, Ser 80 B, Asn 81。这样一来，残基编号体制保持与野生型结构一致。

也就是三元组（hetfield，resseq，icode）。

比如说，上面的葡萄酸残基id，可以是(’H_GLC’, 100, ’A’)；

如果异质标签和插入码为空，那么可以只使用序列标识符：

这其实就是一个正常的残基了

# Full id residue=chain[(' ', 100, ' ')] # Shortcut id residue=chain[100]

异质标签的起因是许许多多的PDB文件使用相同的序列标识符表示一个氨基酸和一个异质残基或一个水分子，如果不使用异质标签的话，这会产生一个很明显的问题。

比如说，举一个例子：
插入码为空的Asn 10具有残基id (’ ’, 10, ’ ’) ；Water 10，残基id (’W’, 10, ’ ’)；一个序列标识符为10的葡萄糖分子（名称为GLC的异质残基），残基id为 (’H_GLC’, 10, ’ ’) 。在这种情况下，三个残基（具有相同插入码，都是空和序列标识符10）可以位于同一条链上，因为它们的残基id是不同的，也就是说如果我们不设置三元组中的第一个项，那么这3种情况将无法区分。

也就是为了将3种情况（10, ’ ’)区分开来，我们才定义了hetfield，表明10号位置上这3种情况分别是：原始的残基、水分子、杂原子。

因为理论上来讲，（10, ’ ’)可以作为一个key，value我们可以填写该残基真实的字符串，也就是说按照正常逻辑来讲，我们可以下面这样：

一个残基对象存储着一个子原子集，它还包含一个表示残基名称的字符串（如 “ASN”）和残基的片段标识符（这对X-PLOR的用户来说很熟悉，但是在SMCRA数据结构的构建中没用到）。 就是此处的片段标识符的延伸。

大多数情况下，hetflag和插入码均为空，如 (’ ’, 10, ’ ’) 。在这些情况下，序列标识符可以用作完整id的快捷方式：

就像前面展示的那样，因为这个就是正常的10号位的残基，也就是正常的序列标识符，所以就按照序列下标取切片一样方便。

# use full id res10 = chain[(" ", 10, " ")] # use shortcut res10 = chain[10]

一个链对象中每个残基对象都应该具有唯一的id。但是对含无序原子的残基，要以一种特殊的方式来处理，详见https://biopython.org/docs/latest/Tutorial/chapter_pdb.html#sec-point-mutations。

一个残基对象还有大量其它方法：

residue.get_resname() # returns the residue name, e.g. "ASN" residue.is_disordered() # returns 1 if the residue has disordered atoms residue.get_segid() # returns the SEGID, e.g. "CHN1" residue.has_id(name) # test if a residue has a certain atom is_aa(residue) # 来检验一个残基对象是否为氨基酸

is_aa是在Polypeptide的方法中，

主要就是第2个参数，用于检查是否是标准氨基酸字母

总体，举例来说：

residue = ctcf_noDNA[0]["A"][727] print(residue.get_resname()) # 残基名称 print(residue.get_id()) # 残基id print(residue.is_disordered()) # 是否无序 print(residue.get_segid()) # SEGID

SEGID暂时不清楚？

for atom in residue: print(atom.get_name(), atom.get_coord()) # 原子名称和坐标

print(residue.get_full_id()) # 完整ID信息 print(residue.get_parent().get_id()) # 链ID print(residue.get_parent().get_parent().get_id()) # 模型ID print(residue.get_parent().get_parent().get_parent().id) # 结构ID

residue.has_id('CA') # 检查是否有CA原子

4.5 原子

原子对象储存着所有与原子有关的数据，它没有子类。原子的id就是它的名称（如，“OG”代表Ser残基的侧链氧原子）。在残基中原子id必需是唯一的。

此外，对于无序原子会产生异常，这部分需要参考官网https://biopython.org/docs/latest/Tutorial/chapter_pdb.html#sec-disordered-atoms。

原子id就是原子名称（如 ’CA’ ）。在实践中，原子名称是从PDB文件中原子名称去除所有空格而创建的。

但是在PDB文件中，空格可以是原子名称的一部分。通常，钙原子称为 ’CA…’ 是为了和Cα原子（叫做 ’.CA.’ ）区分开。在这种情况下，如果去掉空格就会产生问题（如统一个残基中的两个原子都叫做 ’CA’ ），所以保留空格。

在PDB文件中，一个原子名字由4个字符组成，通常头尾皆为空格。为了方便使用，空格通常可以去掉（在PDB文件中氨基酸的Cα原子标记为“.CA.”，点表示空格）。为了生成原子名称（然后是原子id），空格删掉了，除非会在一个残基中造成名字冲突（如两个原子对象有相同的名称和id）。对于后面这种情况，会尝试让原子名称包含空格。这种情况可能会发生在，比如残基包含名称为“.CA.”和“CA…”的原子，尽管这不怎么可能。

所存储的原子数据包括原子名称，原子坐标（如果有的话还包括标准差），B因子（包括各向异性B因子和可能存在的标准差），altloc标识符和完整的、包括空格的原子名称。较少用到的项如原子序号和原子电荷（有时在PDB文件中规定）也就没有存储。

为了处理原子坐标，可以用 ’Atom’ 对象的 transform 方法。用 set_coord 方法可以直接设定原子坐标。

1个Atom对象有如下附加方法：

a.get_name() # atom name (spaces stripped, e.g. "CA") a.get_id() # id (equals atom name) a.get_coord() # atomic coordinates a.get_vector() # atomic coordinates as Vector object a.get_bfactor() # isotropic B factor a.get_occupancy() # occupancy a.get_altloc() # alternative location specifier a.get_sigatm() # standard deviation of atomic parameters a.get_siguij() # standard deviation of anisotropic B factor a.get_anisou() # anisotropic B factor a.get_fullname() # atom name (with spaces, e.g. ".CA.")

此处，我以CTCF蛋白第727号位也就是C端最后一位的Arg残基为例进行展示：

residue = ctcf_noDNA[0]['A'][727] # 获取模型0，链A，残基727 for atom in residue: # 每一个属性都可以通过相应的方法获取 print(f"Atom Name原子名字: {atom.get_name()}") print(f"Atom Id原子编号: {atom.get_id()}") print(f"Atom Coordinates原子坐标: {atom.get_coord()}") print(f"Atom Vector原子向量: {atom.get_vector()}") print(f"Atom B-factor原子B因子: {atom.get_bfactor()}") print(f"Atom Occupancy原子占有率: {atom.get_occupancy()}") print(f"Atom Altloc原子替代位置: {atom.get_altloc()}") print(f"Atom Sigatm原子参数标准偏差: {atom.get_sigatm()}") print(f"Atom Siguij各向异性B因子标准偏差: {atom.get_siguij()}") print(f"Atom Anisou各向异性B因子: {atom.get_anisou()}") print(f"Atom Fullname原子全名: {atom.get_fullname()}")

为了表示原子坐标，使用 siguij、各向异性 B 因子和 sigatm 的 Numpy 数组。

涉及到numpy数组，其实就可以进行一系列的向量运算了，

举个Bio.PDB的 Vector 模块功能的例子：

下面这个是Gly甘氨酸，

假设我们要查找Gly残基的Cβ原子的位置，如果存在的话（因为Gly其实是没有Cβ原子，只有Cα原子）。将Gly残基的N原子沿Cα-C化学键旋转-120度（这个度数的数值，可以回顾有机化学中大概正四面体的空间构想进行想象），能大致将其放在一个真正的Cβ原子的位置上。

怎么做呢？就是下面这样使用 Vector 模块中的rotaxis方法（能用来构造一个绕特定坐标轴的旋转）：

# get atom coordinates as vectors n = residue['N'].get_vector() c = residue['C'].get_vector() ca = residue['CA'].get_vector() # center at origin n = n - ca c = c - ca # find rotation matrix that rotates n # -120 degrees along the ca-c vector rot = rotaxis(-pi * 120.0/180.0, c) # apply rotation to ca-n vector cb_at_origin = n.left_multiply(rot) # put on top of ca atom cb = cb_at_origin+ca

此处我还是以我的例子进行解释，CTCF蛋白的第3号位正好就是Gly，可以拿来演示：

我们先获取每N、C、以及Cα原子的空间三维坐标，就把它当成numpy数组进行操作运算即可。

然后我们再以Cα作为坐标原点，也就是将参考坐标系以Cα为参考，也就是得到其他原子相对于Cα的相对坐标向量。

from Bio.PDB.vectors import rotaxis import math rot = rotaxis(-math.pi * 120.0/180.0, c)

这个rot矩阵，如果要深究的话，其实我们可以探讨线性代数里的坐标系变换，其实就是一个变换，此处不表。

总之我们现在拿到了变换，然后再左乘变换矩阵，对 N 原子的相对向量执行旋转

然后，将虚拟 Cβ 的坐标映射回全局坐标系，

现在我们已经拿到了旋转变换前后的N原子坐标，然后可以在三维坐标系中进行可视化

总体执行下来就是：

# get atom coordinates as vectors 获取原子坐标作为向量 from Bio.PDB.vectors import rotaxis from math import pi n = residue['N'].get_vector() c = residue['C'].get_vector() ca = residue['CA'].get_vector() # center at origin 以ca为参考, 即ca在原点, 获取相对坐标 n = n - ca c = c - ca # find rotation matrix that rotates n # -120 degrees along the ca-c vector 拿到坐标系变换矩阵, 围绕ca-c向量旋转-120度 rot = rotaxis(-pi * 120.0/180.0, c) # apply rotation to ca-n vector 对 N 原子的相对向量执行旋转 cb_at_origin = n.left_multiply(rot) # put on top of ca atom 将虚拟 Cβ 的坐标映射回全局坐标系 cb = cb_at_origin+ca

这个例子展示了在原子数据上能进行一些相当不平凡的向量运算，这些运算会很有用。除了所有常用向量运算（叉积（用 ** ），点积（用 * ），角度， 取范数等）和上述提到的 rotaxis 函数，Vector 模块还有方法能旋转（ rotmat ）或反射（ refmat ）一个向量到另外一个向量上。

4.6 从结构中提取指定的Atom/Residue/Chain/Model

这个其实就是select选择语法，

举些例子如下：

其实就是python数据结构操作的一些简单语法，像切片之类

structure = ctcf_noDNA model = structure[0] chain = model['A'] residue = chain[100] atom = residue['CA']

还可以用一个快捷方式：

atom = structure[0]['A'][100]['CA'] atom

5，无序Disorder

Bio.PDB能够处理无序原子和点突变（比如Gly和Ala残基在相同位置上）。

5.1 一般性方法

无序可以从两个角度来解决：原子和残基的角度。一般来说，我们尝试压缩或封装所有由无序引起的复杂性。如果我们仅仅想遍历所有Cα原子，那么我们不必在意一些具有无序侧链的残基。另一方面，应该考虑在数据结构中完整地表示无序性。因此，无序原子或残基存储在特定的对象中，这些对象表现得就像不存在无序。这可以通过仅表示无序原子或残基的子集来完成。至于挑选哪个子集（例如使用Ser残基的哪两个紊乱OG侧链原子位置），由用户来决定。

5.2 无序原子

无序原子由普通的 Atom 对象表示，但所有表示相同物理原子的 Atom 对象都存储在一个 DisorderedAtom 对象中。在 DisorderedAtom 对象中的每个 Atom 对象都可以使用其 altloc 指定符进行唯一索引。 DisorderedAtom 对象将所有未捕获的方法调用转发到选定的 Atom 对象，默认为表示具有最高占位率的原子。当然，用户可以更改选定的 Atom 对象，利用其 altloc 指定符。通过这种方式，原子无序性得以正确表示，而不会带来太多额外的复杂性。换句话说，如果你对原子无序性不感兴趣，你将不会受到它的困扰。

每个无序原子都有一个特征性的altloc标识符。我们可以设定：一个 DisorderedAtom 对象表现得像与一个指定的altloc标识符相关的 Atom 对象：

5.3 无序残基

• 常见情况：最常见的情形是一个包含一个或多个无序原子的残基。这显然可以通过使用 DisorderedAtom 对象来表示无序原子，并将 DisorderedAtom 对象像普通 Atom 对象一样存储在 Residue 对象中来解决。DisorderedAtom 将表现得和普通原子完全一样（实际上是占有率最高的原子），它将所有未捕获的方法调用转发到它所包含的 Atom 对象（选定的 Atom 对象）之一。
• 点突变：这部分参考https://biopython.org/docs/latest/Tutorial/chapter_pdb.html#point-mutations

6，非标准残基 Hetero residues

7, 浏览Structure对象

解析PDB/mmCIF文件，提取一些Model、Chain、Residue和Atom对象。

此处我还是以CTCF作为例子进行分析，

解析mmCIF文件，提取一些Model、Chain、Residue和Atom对象

from Bio.PDB import MMCIFParser parser = MMCIFParser() ctcf_noDNA = parser.get_structure('CTCF_noDNA', 'CTCF_noDNA.cif') structure = ctcf_noDNA model = structure[0] # 获取第一个模型 chain = model['A'] # 获取链A residue = chain[1] # 获取残基编号1也就是第1个残基 atom = residue['CA'] # 获取CA原子 structure, model, chain, residue, atom

迭代遍历一个结构中的所有原子

如果是看残基：

如果是看原子，其实就是简单的循环语句，

for model in structure: for chain in model: for residue in chain: for atom in residue: print(atom)

有个快捷方式可以遍历一个结构中所有原子：

类似地，遍历一条链中的所有原子，可以这么做：

atoms = chain.get_atoms() for atom in atoms: print(atom)

遍历模型中的所有残基

同样，如果我们想要遍历一遍模型中的残基，

residues = model.get_residues() for residue in residues: print(residue)

我们也可以用 Selection.unfold_entities 函数来获取一个结构的所有残基：

from Bio.PDB import MMCIFParser, Selection res_list = Selection.unfold_entities(structure, 'R') res_list

或者获得链上的所有原子：

atom_list = Selection.unfold_entities(chain, 'A') atom_list

明显的是， A=atom, R=residue, C=chain, M=model, S=structure 。

而且，我们可以用这种标记返回层次中的上层，如从一个 Atoms 列表得到（唯一的） Residue 或 Chain 父类的列表：

从原子倒推残基：

或者从原子倒推链：

更多遍历方法，参考官方的API文档。

从链中提取异质残基（如resseq 10的葡萄糖（GLC）部分）

residue_id = ("H_GLC", 10, " ") residue = chain[residue_id]

打印链中所有异质残基

结合我们前面讲解到的关于残基存储数据结构的三元组，这个问题其实就很好解决了。

就是看第一位的hetfiled是什么内容进行判断，

normal_residue = 0 for residue in chain.get_list(): residue_id = residue.get_id() hetfield = residue_id[0] if hetfield[0]=="H": print(residue_id) elif hetfield[0] == " ": normal_residue += 1 print("Number of normal residues:", normal_residue) print("Ratio of normal residues to total residues:", normal_residue/len(chain.get_list()))

比如说我手头上的这个蛋白质，727个残基都是正常的，

输出一个结构分子中所有B因子大于50的CA原子的坐标

B-factor是一个结构生物学的概念，简单解释：

B 因子反映原子在晶体中的热运动幅度（B 因子越大，原子运动越剧烈，结构稳定性可能越低）。筛选 “B 因子> 50 的 CA 原子”，常用于：

• 识别蛋白质结构中柔性较高的区域（如无规则卷曲）；
• 排查晶体结构解析中可能存在误差的区域（高 B 因子可能伴随坐标精度低）。

某种程度上，每个残基的B-factor数据可以用于IDR建模等。

其实看了前面这么多函数我们都可以发现，每一个对象或者类，其属性的逐层访问，来来去去就只有那么几个函数，

比如说先判断有没有这个属性，has？

如果有这个属性，基本上就是get方法，可以从上往下也可以从下往上倒推。

for model in structure.get_list(): for chain in model.get_list(): for residue in chain.get_list(): if residue.has_id("CA"): ca = residue["CA"] if ca.get_bfactor() > 50.0: print ca.get_coord()

同样的，其实我们可以做一些简单的数据统计，比如说做得更深一点：

ca_number = 0 for model in structure.get_list(): for chain in model.get_list(): for residue in chain.get_list(): if residue.has_id("CA"): ca = residue["CA"] if ca.get_bfactor() > 50.0: ca_number += 1 print(ca.get_coord()) print(f"number of CA with B-factor > 50.0: {ca_number}") print(f"ratio of CA with B-factor > 50.0: {ca_number / (len(Selection.unfold_entities(structure, 'A')))}")

输出的原子坐标我们其实可以在三维空间中查看位置分布，比如说看看是不是在外周还是在某个区域等等。

然后我这里是做了一个统计，就是看看满足这个阈值条件的原子数目和比例有多少？
可以看到，比例才7%左右，还是比较小的。

输出所有含无序原子的残基

for model in structure.get_list(): for chain in model.get_list(): for residue in chain.get_list(): if residue.is_disordered(): resseq = residue.get_id()[1] resname = residue.get_resname() model_id = model.get_id() chain_id = chain.get_id() print(model_id, chain_id, resname, resseq)

我的数据中并没有，当然这个数据中的概念标注，和我们一般分析的IDR区域还是有点区别的，

遍历所有无序原子，并选取所有具有altloc A的原子（如果有的话）

for model in structure.get_list(): ... for chain in model.get_list(): ... for residue in chain.get_list(): ... if residue.is_disordered(): ... for atom in residue.get_list(): ... if atom.is_disordered(): ... if atom.disordered_has_id("A"): ... atom.disordered_select("A")

从 Structure 对象中提取多肽

为了从一个结构中提取多肽，需要用 PolypeptideBuilder 从 Structure 构建一个 Polypeptide 对象的列表，如下所示：

from Bio.PDB.Polypeptide import PPBuilder ppb = PPBuilder() model_nr = 1 polypeptide_list = ppb.build_peptides(structure,model_nr) for polypeptide in polypeptide_list: print(polypeptide)

Polypeptide对象正是Residue对象的一个UserList，总是从单结构域（在此例中为模型1）中创建而来。我们可以用所得 Polypeptide 对象来获取序列作为 Seq 对象，或获得Cα原子的列表。

polypeptide.get_sequence()

前面是通过从structure中构建一个polypeptide，

下面是通用从化学键角度构建多肽对象。

多肽可以通过一个C-N 化学键或一个Cα-Cα化学键距离标准来建立。

from Bio.PDB.Polypeptide import PPBuilder, CaPPBuilder ppb = CaPPBuilder() ?ppb

from Bio.PDB.Polypeptide import PPBuilder, CaPPBuilder # Using C-N ppb=PPBuilder() for pp in ppb.build_peptides(structure): print(pp.get_sequence()) # Using CA-CA ppb=CaPPBuilder() for pp in ppb.build_peptides(structure): print(pp.get_sequence())

需要注意的是，上例中通过 PolypeptideBuilder 只考虑了结构的模型 0。尽管如此，还是可以用 PolypeptideBuilder 从 Model 和 Chain 对象创建 Polypeptide 对象。

获取结构的序列

要做的第一件事就是从结构中提取所有多肽（如上所述）。然后每条多肽的序列就容易从 Polypeptide 对象获得。该序列表示为一个Biopython Seq 对象，它的字母表由 ProteinAlphabet 对象来定义。

按照上面说法，正常处理顺序：从结构中提取所有多肽，从每条多肽polypeptide对象中获取序列数据。

8，分析结构

8.1 度量距离

重载原子的减法运算来返回两个原子之间的距离。

# Get some atoms ca1 = residue1['CA'] ca2 = residue2['CA'] # Simply subtract the atoms to get their distance distance = ca1-ca2

此处我还是以我自己的数据，CTCF蛋白质为例，我想简单看一下N端第1个位置上的残基和C端最后1个位置上的残基之间的距离（在这个结构背景中）。

8.2 度量角度

用原子坐标的向量表示，和 Vector 模块中的 calc_angle 函数可以计算角度。

此处我还是以自己的蛋白为例子，我想随便看第1、566、727这3个位置上的残基构成的角度。

from Bio.PDB.vectors import calc_angle atom1, atom2, atom3 = chain[1]['CA'], chain[566]['CA'], chain[727]['CA'] vector1 = atom1.get_vector() vector2 = atom2.get_vector() vector3 = atom3.get_vector() angle = calc_angle(vector1, vector2, vector3) angle

8.3 度量扭转角

用原子坐标的向量表示，然后用 Vector 模块中的 calc_dihedral 函数可以计算角度。

同样，只不过计算扭转角需要4个点，

from Bio.PDB.vectors import calc_dihedral atom1, atom2, atom3, atom4 = chain[1]['CA'], chain[288]['CA'], chain[566]['CA'], chain[727]['CA'] vector1 = atom1.get_vector() vector2 = atom2.get_vector() vector3 = atom3.get_vector() vector4 = atom4.get_vector() angle = calc_dihedral(vector1, vector2, vector3, vector4) angle

三，多链复合物结构分析：我的一个例子

下面以我手头上的一个多链复合物结构预测文件为例：

2个CTCF蛋白质，22个Zn离子，2条互补的DNA单链——》随便模拟的1个多链复合物结构

我们先导入这个结构

fromBio.PDBimportMMCIFParser,Selection parser=MMCIFParser()CTCF_2_DNA=parser.get_structure('CTCF_2_DNA','/data2/AlphaFold3-SeqVisToolkit/tests/test_data/Multi_protein_complex/fold_ctcf_2_dna_model_1.cif')

同样，只有一个model id，我们就用structure[0]

structure=CTCF_2_DNA model=structure[0]model.child_dict

可以看到，正好26条chain，对应上了（26=2蛋白质+22锌离子+2 DNA单链）

因为每一条chain都是一个实体类，所以我们可以看看具体每一个chain有什么属性，比如说我们以chain A为例，就是一个CTCF单体蛋白；

下面看的是实例属性：

# 遍历 __dict__，打印 属性名: 属性值forattr_name,attr_valueinmodel['A'].__dict__.items():print(f"{attr_name}:{attr_value}")

总体和我们前面介绍的一致，我们主要关心的属性是child_list和child_dict。

当然，光看dict是不够的，因为dict只包含实例属性，不包含方法、类属性或通过slots定义的属性。

想要看更多的类属性、方法，可以参考LibInspector—为小白智能解析、阅读（几乎所有）Python工具库中的第一、二小节，

我们可以先用dir查看一下：

object=model['A']forattrindir(object):ifattr.startswith('_'):continueprint(f"{attr} ——> {getattr(object, attr)}")

我们先看一下这个链的所有残基：

同样的，res.id残基id是三元组的存储数据格式，我们先来检查一下数据

forchaininmodel:forresinchain: print(res.id[0])

我们可以看到输出有Zn锌原子杂原子，

我们回过头来看一下我们的这几条链chain

forchaininmodel: print(f"{chain.id} ——> {len(chain)}")

chain A、B是蛋白质，

C-X是配体Zn原子，Y、Z是两条DNA单链。

我们先看蛋白质结构数据，其中是没有杂原子的：

再看核酸结果，同样是没有杂原子的：

剩余的就是配体，我们可以看到，确实每一个配体原子都是按照杂原子的形式存储的：

forchaininmodel:ifchain.idin['A','B','Y','Z']:continueforresinchain: print(res.id[0])

因为考虑到杂原子在这里指的一般就是配体了，而我们的结构分析结果肯定得保留这一部分配体，

至于水分子，因为我们一般不会考虑水分子在这个蛋白质-DNA互作结构数据中的展示，所以这里我就全过滤掉了。

总之，现在我们的结构数据中的各种chain，

AB是蛋白质，YZ是DNA，其余是Zn2+配体，对于配体我们上面的异原子可以很快检测出来；

那么对于蛋白质和核酸DNA，我们又该如何区分呢？

forchaininmodel:ifchain.idin['A','B']:forresinchain: print(res, res.id)elifchain.idin['Y','Z']:forresinchain: print(res, res.id)

从结果输出中看，我们是能够直接从残基的id层面直接区分的：

蛋白质的话，很简单，可以使用Polypeptide模块：

参考：https://biopython.org/docs/1.75/api/Bio.PDB.Polypeptide.html

要先导入这个多肽类，

from Bio.PDBimportPolypeptideforchaininmodel:ifchain.idin['A','B']:forresinchain: print(res, Polypeptide.is_aa(res,standard=True))elifchain.idin['Y','Z']:forresinchain: print(res, Polypeptide.is_aa(res,standard=True))

效果如下：

当然，我们能够用这个函数确定一个残基是否是氨基酸（以及是否是标准的），但是核酸呢？

我们之所以在这里需要确定这些残基的身份类型，

其中一大原因在于我们需要依据残基的身份类型，来确定相应的操作逻辑，

比如说我要计算残基之间的距离，对氨基酸，我就可以选择α-C，但是核酸呢，一般是选择编号为1的C位，其他类型的残基就更不一样了，所以为了后续计算的需求，我们需要先确定当前残基的身份。

from Bio.PDBimportPolypeptideforresinmodel['A']or model['B']:ifPolypeptide.is_aa(res,standard=True): print(res.get_resname(), res.child_dict)

首先是蛋白质：

我们以氨基酸的CA也就是α-C作为其表征原子（比如说距离计算时，可以用该原子表征该氨基酸）

我们可以先看一下这些表征原子的三维坐标

from Bio.PDBimportPolypeptideforresinmodel['A']or model['B']:ifPolypeptide.is_aa(res,standard=True): print(res.get_resname(), res['CA'], res['CA'].get_coord())

接下来，仿照氨基酸一样原理，我们可以看下核酸残基的每一个原子如何，说不定可以从中找到表征的方法：

forresinmodel['Y']ormodel['Z']:print(res,res.get_atoms())

但是下面这残基：我们根本就没法看出来是什么

也就是我们借助下一层数据（atom类）的思路，如果只是调用get_atoms函数是看不出来什么的

换个正确的遍历方法：

forchain_idin['Y','Z']: chain=model[chain_id]forresinchain: print(res.id, res.get_resname())

我们不调用get_atom方法，我们直接for循环，从原子自身层面上看看有没有什么特征：

forchain_idin['Y','Z']: chain=model[chain_id]forresinchain:foratominres: print(chain_id, res.id, atom.name)

关键的1号位C

如果我们用1号位的C去限制DNA，

forchain_idin['Y','Z']:chain=model[chain_id]forresinchain:if"C1'"inres:print(chain_id,res.id)

我们其实可以发现：正好是将42位长度的残基都整理了出来，正好是DNA双链长度

或者，其实实在没办法我们可以先定义1个核酸规范残基名称的列表或者是字典，然后再用resname判断逻辑；

# 定义常见的核酸残基名称nucleic_acids={'DA','DT','DC','DG','DI','A','C','G','U','I'}res_name=res.resname.strip()# 判断是否可能是核酸is_nucleic=res_nameinnucleic_acids

或者稍微复杂一点，其实我们可以仅仅从残基的resname这个类属性上去解析一切残基分子的身份，只要我们为每一个残基定义的名称列表足够：例如

def obtain_resname(res):ifres.resname[:2]=="CA":resname="CA"elifres.resname[:2]=="FE":resname="FE"elifres.resname[:2]=="CU":resname="CU"else: resname=res.resname.strip()ifresnameinMETAL:return"M"else:returnresname

事实上，很多方法，都是直接从残基名字中进行推断的

# Only count RNA residuesnum_residues=0forresinc:if(res.resname.strip()=='A'or res.resname.strip()=='C'or res.resname.strip()=='G'or res.resname.strip()=='U'): num_residues+=1

总之，除了蛋白质可以用polypeptide函数来判断是否是一个残基之外，其实最广泛最普遍的做法，就是看resname也就是残基名字了，我们大胆点可以像下面这样定义：

protein={'CYS':'C','ASP':'D','SER':'S','GLN':'Q','LYS':'K','ILE':'I','PRO':'P','THR':'T','PHE':'F','ASN':'N','GLY':'G','HIS':'H','LEU':'L','ARG':'R','TRP':'W','ALA':'A','VAL':'V','GLU':'E','TYR':'Y','MET':'M',"UNK":"X"}dna={'DA':'A','DC':'C','DG':'G','DT':'T'}rna={'A':'A','C':'C','G':'G','U':'U'}# 然后我们可以直接抽取第1个残基的name来确定该chain或者是原子的具体类型def get_type(self): first_res=self.child_list[0].resname.strip()# Get the first residue name to see what kind of sequence it isiffirst_resinself.protein:#if the first residue is in the protein dictionary, we have a protein sequencereturn"protein"eliffirst_resinself.dna:#if the first residue is in the dna dictionary, we have a dna sequencereturn"dna"else:return"rna"

我的判断逻辑代码如下：

elif"C1'"inres: rep_atom=res["C1'"]rname=res.resname.strip()ifrnamein{'DA','DC','DG','DT'}: rtype="DNA"elifrnamein{'A','C','G','U'}: rtype="RNA"else: rtype="Nucleic"# Fallback for modified bases

至于其余的，在我们当前这个结构中，其实就是配体了，

然后配体的话，我这里统一取第1个原子

forchaininmodel:ifchain.idin['A','B','Y','Z']:continueforresinchain: print(res, list(res.get_atoms())[0])

四，以AlphaFold的预测结果文件为例

1，Alphafold3输出文件组织

以我手头上的一个预测结果文件为例，

我们下载下来的数据，一般组织如下：

一般是5个model，也就是5次采样扩散结果，我们一般取model 0，置信度最高。

• full_data是我们python数据分析中用于解析的，有各种详细矩阵数据

• model cif文件不用说，用于mmCIF parse解析，就是结构model数据

• summary_confidences 置信数据，主要是一些简短的统计指标

我们数据分析，主要是用cif解析结构+full_data详细注释。

至于request json，就是这个提交job的详细输入设置，具体可以参考我之前的博客使用Json文件批量执行Alphafold3 web server预测；

具体格式如下：

其实就是1个蛋白质序列，1个DNA双链（两条单链），11个Zn离子（ion）

[{"name":"FLZnwitha6cseDNA","modelSeeds":["1960824427"],"sequences":[{"proteinChain":{"sequence":"MEGDAVEAIVEESETFIKGKERKTYQRRREGGQEEDACHLPQNQTDGGEVVQDVNSSVQMVMMEQLDPTLLQMKTEVMEGTVAPEAEAAVDDTQIITLQVVNMEEQPINIGELQLVQVPVPVTVPVATTSVEELQGAYENEVSKEGLAESEPMICHTLPLPEGFQVVKVGANGEVETLEQGELPPQEDPSWQKDPDYQPPAKKTKKTKKSKLRYTEEGKDVDVSVYDFEEEQQEGLLSEVNAEKVVGNMKPPKPTKIKKKGVKKTFQCELCSYTCPRRSNLDRHMKSHTDERPHKCHLCGRAFRTVTLLRNHLNTHTGTRPHKCPDCDMAFVTSGELVRHRRYKHTHEKPFKCSMCDYASVEVSKLKRHIRSHTGERPFQCSLCSYASRDTYKLKRHMRTHSGEKPYECYICHARFTQSGTMKMHILQKHTENVAKFHCPHCDTVIARKSDLGVHLRKQHSYIEQGKKCRYCDAVFHERYALIQHQKSHKNEKRFKCDQCDYACRQERHMIMHKRTHTGEKPYACSHCDKTFRQKQLLDMHFKRYHDPNFVPAAFVCSKCGKTFTRRNTMARHADNCAGPDGVEGENGGETKKSKRGRKRKMRSKKEDSSDSENAEPDLDDNEDEEEPAVEIEPEPEPQPVTPAPPPAKKRRGRPPGRTNQPKQNQPTAIIQVEDQNTGAIENIIVEVKKEPDAEPAEGEEEEAQPAATDAPNGDLTPEMILSMMDR","count":1}},{"dnaSequence":{"sequence":"TTGATTTCTTACCTTTTGGAGCCGCATGATGTCGCTGTCTAC","count":1}},{"ion":{"ion":"ZN","count":11}},{"dnaSequence":{"sequence":"GTAGACAGCGACATCATGCGGCTCCAAAAGGTAAGAAATCAA","count":1}}]}]

2，Full_data的解析

此处我还是以上面例子的full_data为例，进行解析：

importjson struc_file_path="/data2/AlphaFold3-SeqVisToolkit/tests/test_data/CTCF_withDNA/fold_flznwitha6csedna_full_data_0.json"withopen(struc_file_path,"r")asf:struc_matrix=json.load(f)struc_matrix

我们先看一下这里有哪些数据：

atom_chain_ids

1个list，长度为结构中原子数，值为每一个原子对应的链ID，比如说是A chain，或者是B chain等

可以通过集合简单查看其中有多少个chain，

set(struc_matrix['atom_chain_ids']),len(set(struc_matrix['atom_chain_ids']))

比如说我这个例子里一共有14条chain：

符合我们前面的数据（1条蛋白质序列，2条DNA单链序列，11个Zn离子，1+2+11=14），并且我们可以看得出来，DNA是单链单独建模的，离子则是单独一个建模的

另外用字典简单统一一下这个结构中每一条chain有多少个原子：

fromcollectionsimportCounter Counter(struc_matrix['atom_chain_ids'])

最高的那个是蛋白质CTCF，一共有5784个氨基酸原子；然后是两条DNA单链，最后是1个Zn离子单独一条链，单独1个原子

atom_plddts

1个list，长度为结构中原子数目，值为每一个原子具体的pLDDT置信值；

我们可以依据这个指标，将每一个原子划分一下置信类别，然后再上色，就可以绘制出AlphaFold一样的效果

contact_probs

1个list（实际应该视为array数据），其实就是一个矩阵，但是注意，这里是token数xtoken数的矩阵，不是原子数。

importnumpyasnp np.asarray(struc_matrix['contact_probs'])

我们可以看到822，确实不是原子数，这里指的是token数

原子数目的话，我们其实很好理解，就是一个原子一个单位；

那这里的token应该如何理解呢？其实就是一个残基、碱基，或者是一个重原子！

那么这里的822=727（CTCF蛋白有727个残基）+11Zn原子+42x2（DNA双链，共42个碱基对），正好对应上；

所以，这里的token，我们应该理解为残基/生物分子基本单位级别的，一般高于atom原子级别。

那么这里的contact prob，其实就是每一个残基之间距离关系的一个映射，描述的就是残基之间的距离关系。

值都是0-1之间，描述任意两个残基之间的互作概率。

我们可以简单地将其进行可视化

importnumpyasnp np.asarray(struc_matrix['contact_probs'])importmatplotlib.pyplotasplt plt.imshow(np.asarray(struc_matrix['contact_probs']),cmap='hot',interpolation='nearest')plt.colorbar()

可以发现互作概率绝大多数残基之间都是很低的（注意，我这里说的是残基，并不局限于蛋白质序列残基，因为全长822包括727的蛋白质序列，当然我这里可以将727的蛋白质序列单独画1个热图出来），

所以我们很自然地想要看一下AlphaFold对于contact的定义是什么？

结构方面的contact，基本上都是threshold，就是划1个阈值，比如说这里的8A，

当然，常规的想法，肯定是硬阈值，8以内是1，8以外是0？

但是我们看到实际的contact probs数据中其实是有小数的，

所以我们得回到AlphaFold3的原文中进行查看，

参考补充资料：https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-024-07487-w/MediaObjects/41586_2024_7487_MOESM1_ESM.pdf

结合 AlphaFold 3 官方补充资料，对 contact prob（接触概率）的解释其实就很容易理解了：

其实说白了，就是任意两个token之间的距离，这也是一个分布数据，我们不应该视为绝对数值（单一数值当然对于一个硬阈值的逻辑判断就是“是或者不是”），如果是一个分布，那么对于一个分布而言，一个硬阈值也只能划分出来阈值左右两侧数据的分布。

contact prob的准确定义与计算（基于文档4.4节、5.7节）

contact prob（接触概率）是模型预测的“两个token的代表原子距离≤8Å”的概率值，是从distogram（距离分布直方图）中推导的关键权重参数，用于后续全局gPDE评分（公式16）的计算，其核心细节如下：

1. 依赖token代表原子：文档4.4节明确，计算接触概率时，需先确定每个token的“代表原子”（如标准氨基酸用Cα/Cβ、标准核苷酸用C1’、配体用单原子token自身），仅基于这两个代表原子的距离判断是否“接触”。
2. 从distogram提取概率：模型会先输出所有token对(i,j)的distogram——即距离落在不同“区间（bin）”的概率（如3.25-50.75Å分为38个等宽bin，超50.75Å为第39个bin）。contact prob(i,j)是所有“距离≤8Å的bin”的概率总和（文档5.7节公式16注释），例如若“3.25-4Å”概率0.3、“4-6Å”概率0.5、“6-8Å”概率0.14，则contact prob(i,j)=0.3+0.5+0.14=0.94。
3. 核心作用：gPDE评分的权重：如文档5.7节所示，contact prob作为公式16中的权重p_ij，用于凸显“模型确信接触的token对”对全局评分的贡献——接触概率越高（如0.94），该token对的PDE值（预测距离误差）在计算gPDE时权重越大；若概率接近0（如0.02），则该token对基本不影响全局评分。

为什么用8Å作为“接触”的硬性阈值？（基于文档2.5.1节、6.3节）

8Å并非随意设定，而是AlphaFold 3团队基于生物分子互作的物理特性+模型训练目标确定的，核心原因有两点：

1. 匹配生物分子“有效互作”的实际距离范围：
   文档2.5.1节定义“界面（interface）”时，已明确“链间重原子最小距离＜5Å”为真实互作；但在contact prob计算中放宽到8Å，是为了覆盖“弱互作/间接接触”的场景（如配体与蛋白的疏水作用、核酸与蛋白的远程静电作用），避免遗漏潜在的功能相关token对。同时，文档6.3节提到核酸相关LDDT计算用30Å半径、蛋白用15Å，可佐证“不同分子类型的互作距离阈值需适配大小/作用方式”，8Å是兼顾“特异性”与“完整性”的折中值。
2. 模型训练的一致性与鲁棒性：
   若阈值过窄（如5Å），可能因模型对“距离预测的微小偏差”误判为“不接触”，导致权重丢失；若过宽（如10Å），则会引入大量无关token对，稀释有效信息。8Å是团队通过训练验证的“最优阈值”——既能稳定区分“真实接触”与“随机距离”，又能适配不同分子（蛋白、核酸、配体）的接触特性，保证跨数据集（PDB、蒸馏集）训练的一致性。

当然，这些只是我推断的，毕竟原始文献附录中对于距离数据的使用比较多，8A并不是一个唯一的数据标准，但就像超参数调参一样，这些数据很难一句话解释为什么这么设置。

对照我的疑问，为什么输出是小数（如0.94）而非0/1？（基于文档3.7节、4.3节）

这是由AlphaFold 3的概率预测机制和模型不确定性决定的，完全符合文档中的模型设计逻辑：

1. distogram本质是概率分布，而非确定距离：
   模型并非直接预测“token对(i,j)的距离是7.2Å”，而是输出“距离落在各个bin的概率”（如文档4.4节所述，distogram head的输出是binned概率而非确定值）。contact prob是这些bin的概率总和，自然是0-1之间的小数——例如0.94代表“该token对距离≤8Å的总概率为94%”，而非“100%接触”。
2. 体现模型的不确定性：
   小数结果是模型对“接触与否”的置信度体现。例如0.94意味着“模型高度确信接触，但仍有6%的概率距离＞8Å”（可能源于MSA信息不足、模板缺失等）；若为0.5，则代表模型对“是否接触”判断模糊。这种概率化输出比“0/1硬判断”更能反映真实生物结构的不确定性（如部分动态互作的构象波动），也更利于后续gPDE评分的精准加权。

总而言之，其实我们拿到mmCIF数据，真正进行CIF parse解析之后，对于1个model，我们可以获取任意一对token之间的预测距离；当我们seed设置的越多，获取的数据越多，其实类似于构建了一个分布数据之后，我们也同样可以得到一个token（i，j）之间距离分布的数据，我们再用硬性阈值，就会得到一个是否8A以内contact的概率值。当然，这只是一个类比说法，实际计算contact probs还是参考原文。

另外一个自然而然的问题：分布概率与确定距离？

一个很自然的问题，就是一个seed，采样5次，获得5个model，然后一个model的结果中就有contact probs结果，又有mmcif结构（这个结构可以得到一个确切的距离），那么我应该如何理解contact probs是该一次seed一个model里的距离分布，和该model得到的一个确切距离这两者之间的矛盾关系呢？

结合AlphaFold 3官方补充资料（尤其是3.7节扩散模块、4.4节距离分布图预测、5.7节模型选择），这两者并非“矛盾”，而是模型从“概率预测”到“确定结构”的完整输出链路，核心是“概率分布指导采样，采样结果对应单一实例”，具体逻辑如下：

本质：两者属于模型输出的“不同维度”，无矛盾

首先需明确两者的核心定位——它们是模型针对“token对距离”的两种互补输出，服务于不同目标，不存在逻辑冲突：

逻辑上是从“概率分布”到“确定结构”

同一个seed下的1个model（对应1次采样）中，这两个输出是“先后生成、前者指导后者”的关系，具体流程完全匹配文档3.7节扩散模块和4.4节距离分布预测的设计：

1. 第一步：模型先预测distogram（接触概率的来源）
   在模型主干网络（Pairformer Stack）处理后，Distogram Head会输出所有token对的distogram——即“距离落在不同区间（bin）的概率”（如3.25-50.75Å分为38个bin，超50.75Å为第39个bin，文档4.4节）。
   contact probs就是从这个distogram中提取“距离≤8Å的所有bin的概率和”（文档5.7节），本质是模型对“该token对距离范围”的概率判断，反映“不确定性”。
2. 第二步：基于概率分布，通过扩散模型采样得到确定结构
   模型的Diffusion Module（扩散模块）会以“主干网络输出的单表示/对表示”“distogram隐含的距离概率”为条件，通过200步去噪过程（文档3.7.1节），从纯噪声中采样生成唯一的原子坐标集合{x_pred_l}（文档18算法SampleDiffusion）。
   这个坐标集合就是mmCIF结构的来源——将其按token代表原子（如Cα、C1’）计算距离，得到的就是“确切距离”，本质是从概率分布中采样出的“一个具体实例”，反映“确定性结果”。

对照我的疑问，为什么同一model中会同时存在“概率”和“确定距离”？

这是AlphaFold 3“概率建模+生成式预测”框架的必然结果，核心目的是同时满足“置信度评估”和“结构解析”的需求，这种设计逻辑：

1. contact probs：量化“模型对结构的置信度”
   模型无法保证单次采样的结构“100%符合真实情况”，contact probs（及背后的distogram）是“模型对自身预测的不确定性描述”——例如某token对的contact probs=0.94，说明“模型94%确信它距离≤8Å”，即使单次采样得到的距离是8.2Å（略超阈值），也能通过概率判断“这是小概率偏差，模型仍高度认为它是接触的”（文档4.3节置信头设计逻辑）。
2. 确切距离：提供“可验证的物理结构”
   实验验证（如与PDB结构对比LDDT、计算RMSD）需要“具体的原子坐标”，而非抽象的概率分布。扩散模型的采样过程就是“从概率分布中抽取一个最可能符合物理规律的实例”（文档3.7节，通过加权刚性对齐、键长/键角约束确保合理性），这个实例对应的距离就是mmCIF中的确切值——它是概率分布的“一个具体体现”，而非“与概率矛盾的结果”。

补充：同一seed下5次采样的逻辑

我们前面提到“一个seed采样5次获得5个model”，需额外明确：

• 这里的“5次采样”是同一seed下，扩散模型在去噪过程中加入不同随机噪声的结果（文档3.7.1节，每次采样的噪声ξ_l独立生成），因此5个model会输出5组“contact probs + mmCIF距离”。
• 每组中，contact probs的“概率分布”相对稳定（同一seed下主干网络输出的distogram差异小），但mmCIF的“确切距离”会因采样噪声略有不同（如某token对5次采样的距离可能为7.8Å、8.1Å、7.5Å等）——这恰好体现了“概率分布包含多种可能，单次采样仅对应一种可能”的逻辑，进一步说明两者无矛盾。

这里我做了一个试验，去看看是否真的同一个seed下面不同model获取的contact probs的概率分布是不是稳定的，至少输出的结果上是不是稳定，

我以手头上的一个job为例，

importjson models=[]foriinrange(5):withopen(f"/data2/AlphaFold3-SeqVisToolkit/tests/test_data/1_seed_5_models/fold_p49711_dna_full_data_{i}.json","r")asf:models.append(np.asarray(json.load(f)['contact_probs']))# 比较5个model中任意两个model的contact_probs矩阵的差异foriinrange(5):forjinrange(i+1,5):diff=np.abs(models[i]-models[j])print(f"contact_probs matrix is differnet between model_{i}and model_{j}:{np.any(diff!=0)}")

我的想法很简单，就是contact probs互作概率矩阵，就看看是否这5个model中输出的都一不一样。

结果发现，这些model之间其实contact probs矩阵是一样的：

importjson models=[]foriinrange(5):withopen(f"/data2/AlphaFold3-SeqVisToolkit/tests/test_data/1_seed_5_models/fold_p49711_dna_full_data_{i}.json","r")asf:models.append(np.asarray(json.load(f)['contact_probs']))# 比较5个model中任意两个model的contact_probs矩阵的差异foriinrange(5):forjinrange(i+1,5):diff=np.abs(models[i]-models[j])mean_diff=np.mean(diff)print(f"Mean absolute difference between model_{i}and model_{j}:{mean_diff}")

当然，我没法穷尽所有的生成数据来维持这一个判断，所以这个问题暂时avoid，就是我们尽量避免同一个job跨model比较，最后就是直接选top置信度的即可。

总结

简单来说：

• contact probs是“模型认为该token对距离可能落在什么范围的概率”（像“天气预报说80%概率下雨”）；
• mmCIF中的确切距离是“模型基于概率，采样得到的一个具体距离结果”（像“今天实际下了雨”或“没下雨”，都是80%概率下的可能实例）。
  两者是“概率预测”与“实例结果”的关系，而非矛盾，共同构成了AlphaFold 3从“不确定性量化”到“确定性结构”的完整预测链路。

另外一个问题：计算距离我到底用互作置信数据，还是用唯一采样的距离数据

我的问题本质是“概率置信数据”与“确定结构数据”的应用场景适配，以及跨模型一致性的处理，这个问题其实不是小问题，

1，核心选择原则：根据“后续分析目标”决定用contact prob还是cif距离

两者并非“二选一”，而是服务于不同分析需求，仔细查看原始文献文档的话，文档中明确了它们的功能分工：

关键补充：若后续分析需要“兼顾置信度的硬结构数据”（如筛选“高置信的互作残基对”进行网络分析），可将两者结合——先用contact prob筛选出≥0.7（或自定义阈值，文档中高置信通常≥0.8）的残基对，再用这些对的cif距离进行后续计算，排除低置信的“疑似互作”。

主要问题在于contact probs数据已经用8A进行了定型了，而实际cif结构数据中其实距离并没有限制，

只能说，我们的实际结构cif数据中，对于超过8A的距离，可以使用contact probs的置信处理进行筛选。也就是说，contact probs给我们的是一个距离与8A判断的置信标准，我们要筛选，要么同样用8A，或者用比其更强的条件去筛选。

简单来说，如果要往互作那一边考虑：

• contact prob：以 8Å 为基准，告诉你 “哪些残基对的距离大概率有意义（≤8Å）”；
• cif 距离：提供所有残基对的实际距离，但需通过 contact prob 筛选 ——要么用 “contact prob≥阈值”（对应 8Å 接触的高置信），要么在此基础上加更窄的距离限制（如≤5Å），才能聚焦到有生物学意义的互作数据

2，跨model（不同采样/种子）的置信数据与结构数据一致性处理

同一seed下的不同model（或不同seed的model）理论上会输出不同的contact prob和cif距离，这是模型“概率分布采样”的必然结果（3.7节扩散模型的随机噪声特性），当然我上面的测试只是一个简单的例子，文档5.9节“样本排序”提供了明确的处理方案：

优先选择“高置信的model”作为核心分析对象

文档5.9.3节明确，样本排序需结合pTM（整体结构置信度）、ipTM（界面互作置信度）、pLDDT（局部原子置信度），步骤如下：

1. 对所有model，计算其“整体结构置信度”（如full complex pTM，5.9.1节公式）和“互作相关置信度”（如interface ipTM，5.9.1节公式）；
2. 选择pTM最高、ipTM最高且无明显原子冲突（clashes＜100，5.9.2节定义）的1个model作为“核心model”；
3. 后续分析（无论是用contact prob还是cif距离）均基于该“核心model”，避免跨model的数据混乱。

若需整合多model数据（如提高结果鲁棒性）：用“置信加权”或“多数投票”

若分析需要多model的统计结果（如验证互作的稳定性），可参考文档5.7节“分数聚合”逻辑：

• contact prob整合：对同一残基对，计算所有model的contact prob均值，均值≥0.6（或自定义阈值）视为“稳定高置信互作”；
• cif距离整合：对同一残基对，仅保留“核心model”及其他model中contact prob≥0.7的距离数据，计算中位数（减少极端值影响），再判断是否≤8Å；
• 排除低质量model：若某model的pTM＜0.5（文档中低置信结构阈值），直接排除其数据，避免拉低一致性。

3，总结：流程建议

1. 第一步：筛选核心model
   计算所有model的pTM、ipTM、pLDDT，选择“pTM最高+ipTM最高+无冲突”的model作为核心分析对象（参考文档5.9.3节排序逻辑）。
2. 第二步：按需选择数据类型
   • 若需“互作置信度评估”：用核心model的contact prob；
   • 若需“结构分析（如互作网络、距离统计）”：用核心model的cif距离；
   • 若需“高置信结构分析”：用核心model的contact prob筛选高置信残基对，再用这些对的cif距离。
3. 第三步：跨model验证（可选）
   若需验证结果稳定性，对其他model重复步骤1-2，或用“多model contact prob均值+高置信距离中位数”作为最终结果，确保鲁棒性。

当然，距离数据和分布置信度其实是并行独立的，两者只能说可以用来加强条件并行分析。

简单来说：

pae

1个list，实际上是1个ndarray数据，矩阵，类似于contact probs：

也是计算任意一对token（残基）对之间的指标，token（i，j）就是在i处token对齐，然后j处token预期与真实结构距离差距分布的期望。

因为是矩阵，所以可视化热图很简单：

importnumpyasnp np.asarray(struc_matrix['pae'])importmatplotlib.pyplotasplt plt.imshow(np.asarray(struc_matrix['pae']))plt.colorbar()

token_chain_ids

1个list，数目为结构总token（残基）数目，值为该token所属的chain

有了这个数据，我们就可以知道每一个token也就是每一个残基单位到底是在哪一条链上的了，

那么我们绘制整体的PAE图，既可以全部一起画，也可以在边上标注一个chain的track，或者只绘制某一个chain的PAE，做成类似于下面的这种图：

到时候比如说chain A是蛋白质x，我们就可以直接提取出来单独观察x

token_res_ids

1个list，长度为token数，值为该token所对应的残基编号

3，summary_confidence的解析

这个数据在统计的时候也很重要，同样的，我们导入这个json文件看看

importjson struc_confid_path="/data2/AlphaFold3-SeqVisToolkit/tests/test_data/CTCF_withDNA/fold_flznwitha6csedna_summary_confidences_0.json"withopen(struc_confid_path,'r')asf:confid_data=json.load(f)confid_data

这里的key稍微有点多：

chain_iptm

1个list，长度为chain的数目，值是该chain与其他所有chain的iptm值的均值，也就是本来是chainxchain的iptm方针，现在对行取均值；

数目和前面也对得上，就是14条chain（1蛋白+2DNA+11ion）

然后我们可以用这个iptm均值对特定链的预测结构进行排序

chain_pair_iptm

1个嵌套list，实际上就是ndarray，也就是上面说的方针；

注意这个方针中有None值，

我们来看看是不是每一行的均值就是前面的iptm，

先强制转换为numpy数组，转为float类型可以将None替换为np.nan，然后忽略NA计算每一行均值

importnumpyasnplist(np.nanmean(np.asarray(confid_data['chain_pair_iptm'],dtype=np.float32),axis=1))

两个对比之后，我们可以发现，数据大体上是对应的，就是第一行偏差的有点大，暂时不知道原因，但是逻辑上应该是方针对每一行取均值然后结果为前面的iptm列表

然后，同样的道理，因为是方阵，所以可以简单画一个热图

importmatplotlib.pyplotasplt fig,ax=plt.subplots(figsize=(10,10))ax.imshow(np.asarray(confid_data['chain_pair_iptm'],dtype=np.float64))

但是暂时不知道如何处理中间的这些None值，

importmatplotlib.pyplotasplt fig,ax=plt.subplots(figsize=(10,10))im=ax.imshow(np.asarray(confid_data['chain_pair_iptm'],dtype=np.float64))fig.colorbar(im,ax=ax)

参考：https://www.ebi.ac.uk/training/online/courses/alphafold/inputs-and-outputs/evaluating-alphafolds-predicted-structures-using-confidence-scores/confidence-scores-in-alphafold-multimer/

我们知道pLDDT是score，越高越好，结构越好；PAE是error的期望，越低越好；

pTM和ipTM也都是score，也是越高越好，结构越好；

所以上图，我们大概只能知道chain 0，也就是蛋白质序列，和其他配体的相对位置预测的置信度是比较高的；但是其他配体之间就不好说了。

当然，这个Null的问题其实社区中比较常见：
参考https://github.com/ihmwg/ModelCIF/issues/21

针对iptm以及ptm指标修改的替代方法其实现在也有不少报道，有些打分函数可能会更好，

比如说ipSAE（https://github.com/DunbrackLab/IPSAE）

这一块，只能说是见仁见智吧

chain_pair_pae_min

1个嵌套的list，实际上是ndarray，也就是方针；还是chain数xchain数的方阵，

同样还是有None，可以转float的numpy数组时转换为np.nan；

AlphaFold官方定义是仅限制于链 i 的行和仅限制于链 j 的列中的最低 PAE 值，因为按照前面的说法我们可以知道PAE是token级别的，也就是tokenxtoken（每一对残基），但是现在这个指标chain_pair_pae_min是chain级别的，而每一条chain都有很多个token（残基），所以相当于chain i x chain j有一个token is x token js的子矩阵，然后取这个子矩阵的min值，也就是说取两条链中残基对之间相对位置误差最小的残基对的距离误差期望值。

我们可以很清楚地看到：

14条chainx14条chain，所得到的一个方阵，

我们可以来测试一下到底是不是min值，

我们首先需要token层面的PAE数据，来提取出两个chain集合层面的tokenxtoken的子矩阵，同时得标注出来分别是哪一个chain，然后从这个子矩阵中取min值，检查是否与下面的chain_pair_pae_min对应。

我们此处以A-M/A-N/M-N，也就是蛋白质-双DNA链为对象，

下面进行测试：

importnumpyasnpfromtypingimportDict,Tupledefprocess_token_chains(token_chain_ids):""" 输入1个token_chain_ids列表, 也就是每个残基对应的chain id列表; 返回每一条chain对应的残基索引范围列表. """# convert token_chain_ids to arraytoken_chain_ids=np.asarray(token_chain_ids)# get unique chain ids, and assign index to each chain idunique_chain_ids=np.unique(token_chain_ids)chain_to_index={chain_id:ifori,chain_idinenumerate(unique_chain_ids)}# convert token_chain_ids to index array# original is a 1-dim array, reshape to (1, -1) —— a row vector, (N,) -> (1, N)token_chain_indices=np.asarray([chain_to_index[chain_id]forchain_idintoken_chain_ids]).reshape(1,-1)# get the start and end index for each chain_idchain_to_range:Dict[str,Tuple[int,int]]={}chain_index_to_range:Dict[int,Tuple[int,int]]={}forchain_idinunique_chain_ids:indices=np.where(token_chain_ids==chain_id)[0]chain_to_range[chain_id]=(indices[0],indices[-1])chain_index_to_range[chain_to_index[chain_id]]=(indices[0],indices[-1])returnchain_to_index,chain_to_range,chain_index_to_range,token_chain_indices# 下面开始处理数据importjsonwithopen("/data2/AlphaFold3-SeqVisToolkit/tests/test_data/CTCF_withDNA/fold_flznwitha6csedna_summary_confidences_0.json","r")asf:struc_confid=json.load(f)withopen("/data2/AlphaFold3-SeqVisToolkit/tests/test_data/CTCF_withDNA/fold_flznwitha6csedna_full_data_0.json","r")asf:struc_matrix=json.load(f)chain_to_index,chain_to_range,chain_index_to_range,token_chain_indices=process_token_chains(struc_matrix['token_chain_ids'])# 这个是原始的pae矩阵pae_matrix=np.asarray(struc_matrix['pae'])# 我们从这个矩阵中获取每一条chain与每一条chain之间的最小pae值pae_matrix_chain_min=np.zeros((len(chain_to_index),len(chain_to_index)))fori,(start_i,end_i)inchain_index_to_range.items():forj,(start_j,end_j)inchain_index_to_range.items():sub_pae=pae_matrix[start_i:end_i+1,start_j:end_j+1]sub_min=np.min(sub_pae)pae_matrix_chain_min[i,j]=sub_min# pae_matrix_chain_min是我们从pae中计算出来的, 现在与confidence数据中的chain_pair_pae_min对比一下conf_chain_pair_pae_min=np.asarray(struc_confid['chain_pair_pae_min'])print(pae_matrix_chain_min)print("----")print(conf_chain_pair_pae_min)

对比如下，第一个矩阵是我按照定义从pae子矩阵中自己提取出来的，第二个矩阵是confidence json文件中alphafold自己提供的：

因为confidence数据中有None值，所以系统给出的效果不好，但是从有值的几个地方对比来看，确实基本上数据是符合的。

而且，其实按照定义的话，我们是能够按照自己从PAE矩阵中提取并补完confidence中的chain_pair_pae_min数据的，但是不知道为什么，这么一个小漏洞，alphafold3的输出却没有处理，而是给出了None值？有点小困惑

可视化效果上：下面第一幅图是alphafold3输出，有null值，所以中间空；

第二幅图是我自己依据定义提取出来的，效果完美

importmatplotlib.pyplotasplt plt.imshow(conf_chain_pair_pae_min)plt.colorbar()

plt.imshow(pae_matrix_chain_min)plt.colorbar()

可以看得出来，这些配体（Zn离子）和除蛋白质序列以外的其他任何对象互作、相对位置预测都有很大误差。

较低的值表明，AlphaFold在预测链如何相互作用和相互对齐方面有较高的信心。

chain_ptm

1个list，长度是chain数目，我这里的数据是14条chain，值是该chain的ptm值，我们可以依据ptm score高低对chain进行排序，比如说哪一条chain的score高，可能整条chain预测的结构比较靠谱；

还是有None值，

AlphaFold中出现Null值很频繁，感觉不是个例，我手头上的几个数据的confidence summary json数据中都有Null值，

翻了下官网仓库的issue，有几个问题有关，

https://github.com/google-deepmind/alphafold3/issues/273

1个是用CPU进行推理（没有在GPU上），数值不稳定，然后在PDE输出上出现了NaN这样的null值。

fraction_disordered

1个数值，表示预测结构中有多少比例是无序的，官方定义是通过可及表面积进行计算的

无序比例是预期蛋白质中无序或高度灵活的残基比例的度量。如果pLDDT <50则认为是无序。因此，Fraction_disordered计算如下：（pLDDT < 50的残基数量）/总残基数量。Linker通常被认为是无序的，因为它们具有高迁移率的性质。

has_clash

1个bool值，返回1或0表示，官方定义用于指示结构中是否存在大量冲突原子（超过链的 50%，或具有超过 100 个冲突原子的链）

iptm

1个数值，定义为所有界面预测预测的界面 TM 分数

num_recycles

1个数值，指定AlphaFold为蛋白质序列迭代预测过程的次数

ptm

表示完整结构的预测 TM 分数

ranking_score

1个数值，一个范围在-100 到 1.5 之间的标量，可用于对预测结果进行排序。它将 ptm、iptm、disordered_fraction 和 has_clash 综合为一个数值

4，mmCIF文件的解析

这个文件其实就是我们每次扩散模型采样所生成的一个样本，AlphaFold3所成生的结构数据和我们在PDB或者是其他结构数据库中所存储的其他mmCIF结构文件没有什么区别，

都是按照Structure-Model-Chain-Residue-Atom层级逐层构建的数据结构，

这部分文件的解析，主要是解析其中每一个atom的三维坐标，有了坐标，就有了空间相对远近，就可以进行丰富的作图分析以及绘制了。

这部分工作以及示例，我在此处就不提供代码以及演示了，详情可以参考我在最前面所提到的我开发的那个绘图工具，我将结构解析+空间坐标如何利用+能做什么分析，都集成到了我自己开发的工具中。

当然，空间坐标能够做什么分析，我的工具仅代表我个人观点，仁者见仁智者见智。

另外注意，alphafold3模型给出的结构预测文件，毕竟只是从分布中采出的一个样，只要是采样，其实就有假阳性，至于为什么采这个样，或者说为什么拿这个样作为置信度最好的，置信度的指标并不是万能的。

其实我们仔细看输出结果的cif文件，和其预测的8A置信互作图，就会发现模型输出是存在假阳性情况的。

所以这就是为什么上流的人不看结果，他们更看重模型中间的输出，比如说我就想要距离分布head所输出的那个分布，但是alphafold3并不会在输出数据中提供这个，除非我们自己部署，自己拆分源码，将这个模型输出中的部分数据存下来。

参考

https://biopython.org/docs/latest/Tutorial/chapter_pdb.html（英文原版）

https://github.com/bigwiv/Biopython-cn/blob/095885f51f0f148fa0952bb2b6180693ee39b3be/cn/chr11.rst#L68（中文版开源翻译，但是版本有点老）