企业官网建设流程全解析

app开发软件排行榜,seo排名点击工具,17858833595做网站,旅游网站网页设计代码文章目录先放结论背景二#xff0c;走向3D#xff1a;BioPython的PDB模块1#xff0c;坐标文件读取2#xff0c;文件内容访问#xff08;主要是字典数据#xff09;3#xff0c;写入PDB/mmCIF文件4#xff0c;结构表示4.1 结构4.2 模型4.3 链4.4 残基4.5 原子4.6 从结构…文章目录先放结论背景二走向3DBioPython的PDB模块1坐标文件读取2文件内容访问主要是字典数据3写入PDB/mmCIF文件4结构表示4.1 结构4.2 模型4.3 链4.4 残基4.5 原子4.6 从结构中提取指定的Atom/Residue/Chain/Model5无序Disorder5.1 一般性方法5.2 无序原子5.3 无序残基6非标准残基 Hetero residues7, 浏览Structure对象解析mmCIF文件提取一些Model、Chain、Residue和Atom对象迭代遍历一个结构中的所有原子遍历模型中的所有残基从链中提取异质残基如resseq 10的葡萄糖GLC部分打印链中所有异质残基输出一个结构分子中所有B因子大于50的CA原子的坐标输出所有含无序原子的残基遍历所有无序原子并选取所有具有altloc A的原子如果有的话从 Structure 对象中提取多肽获取结构的序列8分析结构8.1 度量距离8.2 度量角度8.3 度量扭转角三多链复合物结构分析我的一个例子四以AlphaFold的预测结果文件为例1Alphafold3输出文件组织2Full_data的解析atom_chain_idsatom_plddtscontact_probscontact prob的准确定义与计算基于文档4.4节、5.7节为什么用8Å作为“接触”的硬性阈值基于文档2.5.1节、6.3节对照我的疑问为什么输出是小数如0.94而非0/1基于文档3.7节、4.3节另外一个自然而然的问题分布概率与确定距离本质两者属于模型输出的“不同维度”无矛盾逻辑上是从“概率分布”到“确定结构”对照我的疑问为什么同一model中会同时存在“概率”和“确定距离”补充同一seed下5次采样的逻辑总结另外一个问题计算距离我到底用互作置信数据还是用唯一采样的距离数据1核心选择原则根据“后续分析目标”决定用contact prob还是cif距离2跨model不同采样/种子的置信数据与结构数据一致性处理3总结流程建议paetoken_chain_idstoken_res_ids3summary_confidence的解析chain_iptmchain_pair_iptmchain_pair_pae_minchain_ptmfraction_disorderedhas_clashiptmnum_recyclesptmranking_score4mmCIF文件的解析参考先放结论本篇博客所讨论的结构数据处理、以及可视化内容我按照自己的想法已经实现、整合成了1个工具并做了一定的开发。详细使用以及文档可以参考https://github.com/MaybeBio/AlphaFold3-SeqVisToolkit欢迎PR和Issue背景我有了alphafold3预测的结构数据这个数据大多是PDB格式或者是mmCIF格式也就是生物大分子晶体结构数据。现在我们想要解析、挖掘这个结构预测文件中的空间坐标信息希望能够得到一些有用的信息但是怎么访问、处理、解析以及可视化这个文件怎么整合一维的序列特征信息整合到三维的结构可视化中二走向3DBioPython的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(PERMISSIVE1) # 接着通过 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.html2文件内容访问主要是字典数据其中比较重要的一点是结构对象的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_noDNAfrom 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 对象的整体布局遵循称为SMCRAStructure/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 方法来检查这个实体是否有子类具有给定identity.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。这样一来残基编号体制保持与野生型结构一致。也就是三元组hetfieldresseqicode。比如说上面的葡萄酸残基id可以是(’H_GLC’, 100, ’A’)如果异质标签和插入码为空那么可以只使用序列标识符这其实就是一个正常的残基了# Full id residuechain[( , 100, )] # Shortcut id residuechain[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, ’ ’)可以作为一个keyvalue我们可以填写该残基真实的字符串也就是说按照正常逻辑来讲我们可以下面这样一个残基对象存储着一个子原子集它还包含一个表示残基名称的字符串如 “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()) # SEGIDSEGID暂时不清楚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) # 结构IDresidue.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(fAtom Name原子名字: {atom.get_name()}) print(fAtom Id原子编号: {atom.get_id()}) print(fAtom Coordinates原子坐标: {atom.get_coord()}) print(fAtom Vector原子向量: {atom.get_vector()}) print(fAtom B-factor原子B因子: {atom.get_bfactor()}) print(fAtom Occupancy原子占有率: {atom.get_occupancy()}) print(fAtom Altloc原子替代位置: {atom.get_altloc()}) print(fAtom Sigatm原子参数标准偏差: {atom.get_sigatm()}) print(fAtom Siguij各向异性B因子标准偏差: {atom.get_siguij()}) print(fAtom Anisou各向异性B因子: {atom.get_anisou()}) print(fAtom 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_originca此处我还是以我的例子进行解释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_originca这个例子展示了在原子数据上能进行一些相当不平凡的向量运算这些运算会很有用。除了所有常用向量运算叉积用 ** 点积用 * 角度 取范数等和上述提到的 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] atom5无序DisorderBio.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-mutations6非标准残基 Hetero residues7, 浏览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明显的是 Aatom, Rresidue, Cchain, Mmodel, Sstructure 。而且我们可以用这种标记返回层次中的上层如从一个 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(fnumber of CA with B-factor 50.0: {ca_number}) print(fratio 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() ?ppbfrom Bio.PDB.Polypeptide import PPBuilder, CaPPBuilder # Using C-N ppbPPBuilder() for pp in ppb.build_peptides(structure): print(pp.get_sequence()) # Using CA-CA ppbCaPPBuilder() 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) angle8.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 parserMMCIFParser()CTCF_2_DNAparser.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]structureCTCF_2_DNA modelstructure[0]model.child_dict可以看到正好26条chain对应上了262蛋白质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查看一下objectmodel[A]forattrindir(object):ifattr.startswith(_):continueprint(f{attr} —— {getattr(object, attr)})我们先看一下这个链的所有残基同样的res.id残基id是三元组的存储数据格式我们先来检查一下数据forchaininmodel:forresinchain: print(res.id[0])我们可以看到输出有Zn锌原子杂原子我们回过头来看一下我们的这几条链chainforchaininmodel: 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互作结构数据中的展示所以这里我就全过滤掉了。总之现在我们的结构数据中的各种chainAB是蛋白质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,standardTrue))elifchain.idin[Y,Z]:forresinchain: print(res, Polypeptide.is_aa(res,standardTrue))效果如下当然我们能够用这个函数确定一个残基是否是氨基酸以及是否是标准的但是核酸呢我们之所以在这里需要确定这些残基的身份类型其中一大原因在于我们需要依据残基的身份类型来确定相应的操作逻辑比如说我要计算残基之间的距离对氨基酸我就可以选择α-C但是核酸呢一般是选择编号为1的C位其他类型的残基就更不一样了所以为了后续计算的需求我们需要先确定当前残基的身份。from Bio.PDBimportPolypeptideforresinmodel[A]or model[B]:ifPolypeptide.is_aa(res,standardTrue): print(res.get_resname(), res.child_dict)首先是蛋白质我们以氨基酸的CA也就是α-C作为其表征原子比如说距离计算时可以用该原子表征该氨基酸我们可以先看一下这些表征原子的三维坐标from Bio.PDBimportPolypeptideforresinmodel[A]or model[B]:ifPolypeptide.is_aa(res,standardTrue): 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]: chainmodel[chain_id]forresinchain: print(res.id, res.get_resname())我们不调用get_atom方法我们直接for循环从原子自身层面上看看有没有什么特征forchain_idin[Y,Z]: chainmodel[chain_id]forresinchain:foratominres: print(chain_id, res.id, atom.name)关键的1号位C如果我们用1号位的C去限制DNAforchain_idin[Y,Z]:chainmodel[chain_id]forresinchain:ifC1inres:print(chain_id,res.id)我们其实可以发现正好是将42位长度的残基都整理了出来正好是DNA双链长度或者其实实在没办法我们可以先定义1个核酸规范残基名称的列表或者是字典然后再用resname判断逻辑# 定义常见的核酸残基名称nucleic_acids{DA,DT,DC,DG,DI,A,C,G,U,I}res_nameres.resname.strip()# 判断是否可能是核酸is_nucleicres_nameinnucleic_acids或者稍微复杂一点其实我们可以仅仅从残基的resname这个类属性上去解析一切残基分子的身份只要我们为每一个残基定义的名称列表足够例如def obtain_resname(res):ifres.resname[:2]CA:resnameCAelifres.resname[:2]FE:resnameFEelifres.resname[:2]CU:resnameCUelse: resnameres.resname.strip()ifresnameinMETAL:returnMelse:returnresname事实上很多方法都是直接从残基名字中进行推断的# Only count RNA residuesnum_residues0forresinc:if(res.resname.strip()Aor res.resname.strip()Cor res.resname.strip()Gor res.resname.strip()U): num_residues1总之除了蛋白质可以用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_resself.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 sequencereturnproteineliffirst_resinself.dna:#if the first residue is in the dna dictionary, we have a dna sequencereturndnaelse:returnrna我的判断逻辑代码如下elifC1inres: rep_atomres[C1]rnameres.resname.strip()ifrnamein{DA,DC,DG,DT}: rtypeDNAelifrnamein{A,C,G,U}: rtypeRNAelse: rtypeNucleic# Fallback for modified bases至于其余的在我们当前这个结构中其实就是配体了然后配体的话我这里统一取第1个原子forchaininmodel:ifchain.idin[A,B,Y,Z]:continueforresinchain: print(res, list(res.get_atoms())[0])四以AlphaFold的预测结果文件为例1Alphafold3输出文件组织以我手头上的一个预测结果文件为例我们下载下来的数据一般组织如下一般是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}}]}]2Full_data的解析此处我还是以上面例子的full_data为例进行解析importjson struc_file_path/data2/AlphaFold3-SeqVisToolkit/tests/test_data/CTCF_withDNA/fold_flznwitha6csedna_full_data_0.jsonwithopen(struc_file_path,r)asf:struc_matrixjson.load(f)struc_matrix我们先看一下这里有哪些数据atom_chain_ids1个list长度为结构中原子数值为每一个原子对应的链ID比如说是A chain或者是B chain等可以通过集合简单查看其中有多少个chainset(struc_matrix[atom_chain_ids]),len(set(struc_matrix[atom_chain_ids]))比如说我这个例子里一共有14条chain符合我们前面的数据1条蛋白质序列2条DNA单链序列11个Zn离子121114并且我们可以看得出来DNA是单链单独建模的离子则是单独一个建模的另外用字典简单统一一下这个结构中每一条chain有多少个原子fromcollectionsimportCounter Counter(struc_matrix[atom_chain_ids])最高的那个是蛋白质CTCF一共有5784个氨基酸原子然后是两条DNA单链最后是1个Zn离子单独一条链单独1个原子atom_plddts1个list长度为结构中原子数目值为每一个原子具体的pLDDT置信值我们可以依据这个指标将每一个原子划分一下置信类别然后再上色就可以绘制出AlphaFold一样的效果contact_probs1个list实际应该视为array数据其实就是一个矩阵但是注意这里是token数xtoken数的矩阵不是原子数。importnumpyasnp np.asarray(struc_matrix[contact_probs])我们可以看到822确实不是原子数这里指的是token数原子数目的话我们其实很好理解就是一个原子一个单位那这里的token应该如何理解呢其实就是一个残基、碱基或者是一个重原子那么这里的822727CTCF蛋白有727个残基11Zn原子42x2DNA双链共42个碱基对正好对应上所以这里的token我们应该理解为残基/生物分子基本单位级别的一般高于atom原子级别。那么这里的contact prob其实就是每一个残基之间距离关系的一个映射描述的就是残基之间的距离关系。值都是0-1之间描述任意两个残基之间的互作概率。我们可以简单地将其进行可视化importnumpyasnp np.asarray(struc_matrix[contact_probs])importmatplotlib.pyplotasplt plt.imshow(np.asarray(struc_matrix[contact_probs]),cmaphot,interpolationnearest)plt.colorbar()可以发现互作概率绝大多数残基之间都是很低的注意我这里说的是残基并不局限于蛋白质序列残基因为全长822包括727的蛋白质序列当然我这里可以将727的蛋白质序列单独画1个热图出来所以我们很自然地想要看一下AlphaFold对于contact的定义是什么结构方面的contact基本上都是threshold就是划1个阈值比如说这里的8A当然常规的想法肯定是硬阈值8以内是18以外是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的计算其核心细节如下依赖token代表原子文档4.4节明确计算接触概率时需先确定每个token的“代表原子”如标准氨基酸用Cα/Cβ、标准核苷酸用C1’、配体用单原子token自身仅基于这两个代表原子的距离判断是否“接触”。从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.30.50.140.94。核心作用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团队基于生物分子互作的物理特性模型训练目标确定的核心原因有两点匹配生物分子“有效互作”的实际距离范围文档2.5.1节定义“界面interface”时已明确“链间重原子最小距离5Å”为真实互作但在contact prob计算中放宽到8Å是为了覆盖“弱互作/间接接触”的场景如配体与蛋白的疏水作用、核酸与蛋白的远程静电作用避免遗漏潜在的功能相关token对。同时文档6.3节提到核酸相关LDDT计算用30Å半径、蛋白用15Å可佐证“不同分子类型的互作距离阈值需适配大小/作用方式”8Å是兼顾“特异性”与“完整性”的折中值。模型训练的一致性与鲁棒性若阈值过窄如5Å可能因模型对“距离预测的微小偏差”误判为“不接触”导致权重丢失若过宽如10Å则会引入大量无关token对稀释有效信息。8Å是团队通过训练验证的“最优阈值”——既能稳定区分“真实接触”与“随机距离”又能适配不同分子蛋白、核酸、配体的接触特性保证跨数据集PDB、蒸馏集训练的一致性。当然这些只是我推断的毕竟原始文献附录中对于距离数据的使用比较多8A并不是一个唯一的数据标准但就像超参数调参一样这些数据很难一句话解释为什么这么设置。对照我的疑问为什么输出是小数如0.94而非0/1基于文档3.7节、4.3节这是由AlphaFold 3的概率预测机制和模型不确定性决定的完全符合文档中的模型设计逻辑distogram本质是概率分布而非确定距离模型并非直接预测“token对(i,j)的距离是7.2Å”而是输出“距离落在各个bin的概率”如文档4.4节所述distogram head的输出是binned概率而非确定值。contact prob是这些bin的概率总和自然是0-1之间的小数——例如0.94代表“该token对距离≤8Å的总概率为94%”而非“100%接触”。体现模型的不确定性小数结果是模型对“接触与否”的置信度体现。例如0.94意味着“模型高度确信接触但仍有6%的概率距离8Å”可能源于MSA信息不足、模板缺失等若为0.5则代表模型对“是否接触”判断模糊。这种概率化输出比“0/1硬判断”更能反映真实生物结构的不确定性如部分动态互作的构象波动也更利于后续gPDE评分的精准加权。总而言之其实我们拿到mmCIF数据真正进行CIF parse解析之后对于1个model我们可以获取任意一对token之间的预测距离当我们seed设置的越多获取的数据越多其实类似于构建了一个分布数据之后我们也同样可以得到一个tokenij之间距离分布的数据我们再用硬性阈值就会得到一个是否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对距离”的两种互补输出服务于不同目标不存在逻辑冲突输出类型核心定义基于文档来源模块作用contact probs接触概率从distogram距离分布直方图中提取的“token对距离≤8Å的概率总和”文档4.4节、5.7节公式16注释是0-1的概率值Distogram Head距离分布头作为gPDE评分的权重文档5.7节量化“模型对该token对是否接触的置信度”mmCIF中的确切距离从模型最终输出的原子坐标{x_pred_l}中计算的“token代表原子间的实际距离”如蛋白Cα、核酸C1’是具体数值如7.2ÅDiffusion Module扩散模块提供“可解析的物理结构”用于与实验结构对比如LDDT评分、RMSD计算逻辑上是从“概率分布”到“确定结构”同一个seed下的1个model对应1次采样中这两个输出是“先后生成、前者指导后者”的关系具体流程完全匹配文档3.7节扩散模块和4.4节距离分布预测的设计第一步模型先预测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对距离范围”的概率判断反映“不确定性”。第二步基于概率分布通过扩散模型采样得到确定结构模型的Diffusion Module扩散模块会以“主干网络输出的单表示/对表示”“distogram隐含的距离概率”为条件通过200步去噪过程文档3.7.1节从纯噪声中采样生成唯一的原子坐标集合{x_pred_l}文档18算法SampleDiffusion。这个坐标集合就是mmCIF结构的来源——将其按token代表原子如Cα、C1’计算距离得到的就是“确切距离”本质是从概率分布中采样出的“一个具体实例”反映“确定性结果”。对照我的疑问为什么同一model中会同时存在“概率”和“确定距离”这是AlphaFold 3“概率建模生成式预测”框架的必然结果核心目的是同时满足“置信度评估”和“结构解析”的需求这种设计逻辑contact probs量化“模型对结构的置信度”模型无法保证单次采样的结构“100%符合真实情况”contact probs及背后的distogram是“模型对自身预测的不确定性描述”——例如某token对的contact probs0.94说明“模型94%确信它距离≤8Å”即使单次采样得到的距离是8.2Å略超阈值也能通过概率判断“这是小概率偏差模型仍高度认为它是接触的”文档4.3节置信头设计逻辑。确切距离提供“可验证的物理结构”实验验证如与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(i1,5):diffnp.abs(models[i]-models[j])print(fcontact_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(i1,5):diffnp.abs(models[i]-models[j])mean_diffnp.mean(diff)print(fMean 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 prob1. contact prob是从distogram距离分布直方图中提取的“距离≤8Å的概率总和”4.4节、5.7节公式16直接反映模型对“残基是否互作”的置信度 2. 文档5.7节用contact prob作为gPDE评分权重正是因为它能区分“高置信互作”和“低置信互作”。进行“基于结构的硬分析”如计算互作网络、结合口袋距离、与实验结构对比cif文件计算的实际距离1. cif文件的原子坐标是扩散模型采样生成的“确定物理结构”3.7节SampleDiffusion算法其距离是可直接用于结构分析的具体数值 2. 文档6.3节评估时用cif距离计算LDDT、RMSD等指标正是因为需要“确定结构”作为分析基础。关键补充若后续分析需要“兼顾置信度的硬结构数据”如筛选“高置信的互作残基对”进行网络分析可将两者结合——先用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局部原子置信度步骤如下对所有model计算其“整体结构置信度”如full complex pTM5.9.1节公式和“互作相关置信度”如interface ipTM5.9.1节公式选择pTM最高、ipTM最高且无明显原子冲突clashes1005.9.2节定义的1个model作为“核心model”后续分析无论是用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的pTM0.5文档中低置信结构阈值直接排除其数据避免拉低一致性。3总结流程建议第一步筛选核心model计算所有model的pTM、ipTM、pLDDT选择“pTM最高ipTM最高无冲突”的model作为核心分析对象参考文档5.9.3节排序逻辑。第二步按需选择数据类型若需“互作置信度评估”用核心model的contact prob若需“结构分析如互作网络、距离统计”用核心model的cif距离若需“高置信结构分析”用核心model的contact prob筛选高置信残基对再用这些对的cif距离。第三步跨model验证可选若需验证结果稳定性对其他model重复步骤1-2或用“多model contact prob均值高置信距离中位数”作为最终结果确保鲁棒性。当然距离数据和分布置信度其实是并行独立的两者只能说可以用来加强条件并行分析。简单来说pae1个list实际上是1个ndarray数据矩阵类似于contact probs也是计算任意一对token残基对之间的指标tokenij就是在i处token对齐然后j处token预期与真实结构距离差距分布的期望。因为是矩阵所以可视化热图很简单importnumpyasnp np.asarray(struc_matrix[pae])importmatplotlib.pyplotasplt plt.imshow(np.asarray(struc_matrix[pae]))plt.colorbar()token_chain_ids1个list数目为结构总token残基数目值为该token所属的chain有了这个数据我们就可以知道每一个token也就是每一个残基单位到底是在哪一条链上的了那么我们绘制整体的PAE图既可以全部一起画也可以在边上标注一个chain的track或者只绘制某一个chain的PAE做成类似于下面的这种图到时候比如说chain A是蛋白质x我们就可以直接提取出来单独观察xtoken_res_ids1个list长度为token数值为该token所对应的残基编号3summary_confidence的解析这个数据在统计的时候也很重要同样的我们导入这个json文件看看importjson struc_confid_path/data2/AlphaFold3-SeqVisToolkit/tests/test_data/CTCF_withDNA/fold_flznwitha6csedna_summary_confidences_0.jsonwithopen(struc_confid_path,r)asf:confid_datajson.load(f)confid_data这里的key稍微有点多chain_iptm1个list长度为chain的数目值是该chain与其他所有chain的iptm值的均值也就是本来是chainxchain的iptm方针现在对行取均值数目和前面也对得上就是14条chain1蛋白2DNA11ion然后我们可以用这个iptm均值对特定链的预测结构进行排序chain_pair_iptm1个嵌套list实际上就是ndarray也就是上面说的方针注意这个方针中有None值我们来看看是不是每一行的均值就是前面的iptm先强制转换为numpy数组转为float类型可以将None替换为np.nan然后忽略NA计算每一行均值importnumpyasnplist(np.nanmean(np.asarray(confid_data[chain_pair_iptm],dtypenp.float32),axis1))两个对比之后我们可以发现数据大体上是对应的就是第一行偏差的有点大暂时不知道原因但是逻辑上应该是方针对每一行取均值然后结果为前面的iptm列表然后同样的道理因为是方阵所以可以简单画一个热图importmatplotlib.pyplotasplt fig,axplt.subplots(figsize(10,10))ax.imshow(np.asarray(confid_data[chain_pair_iptm],dtypenp.float64))但是暂时不知道如何处理中间的这些None值importmatplotlib.pyplotasplt fig,axplt.subplots(figsize(10,10))imax.imshow(np.asarray(confid_data[chain_pair_iptm],dtypenp.float64))fig.colorbar(im,axax)参考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指标修改的替代方法其实现在也有不少报道有些打分函数可能会更好比如说ipSAEhttps://github.com/DunbrackLab/IPSAE这一块只能说是见仁见智吧chain_pair_pae_min1个嵌套的list实际上是ndarray也就是方针还是chain数xchain数的方阵同样还是有None可以转float的numpy数组时转换为np.nanAlphaFold官方定义是仅限制于链 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_idsnp.asarray(token_chain_ids)# get unique chain ids, and assign index to each chain idunique_chain_idsnp.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_indicesnp.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:indicesnp.where(token_chain_idschain_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_confidjson.load(f)withopen(/data2/AlphaFold3-SeqVisToolkit/tests/test_data/CTCF_withDNA/fold_flznwitha6csedna_full_data_0.json,r)asf:struc_matrixjson.load(f)chain_to_index,chain_to_range,chain_index_to_range,token_chain_indicesprocess_token_chains(struc_matrix[token_chain_ids])# 这个是原始的pae矩阵pae_matrixnp.asarray(struc_matrix[pae])# 我们从这个矩阵中获取每一条chain与每一条chain之间的最小pae值pae_matrix_chain_minnp.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_paepae_matrix[start_i:end_i1,start_j:end_j1]sub_minnp.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_minnp.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_ptm1个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/2731个是用CPU进行推理没有在GPU上数值不稳定然后在PDE输出上出现了NaN这样的null值。fraction_disordered1个数值表示预测结构中有多少比例是无序的官方定义是通过可及表面积进行计算的无序比例是预期蛋白质中无序或高度灵活的残基比例的度量。如果pLDDT 50则认为是无序。因此Fraction_disordered计算如下pLDDT 50的残基数量/总残基数量。Linker通常被认为是无序的因为它们具有高迁移率的性质。has_clash1个bool值返回1或0表示官方定义用于指示结构中是否存在大量冲突原子超过链的 50%或具有超过 100 个冲突原子的链iptm1个数值定义为所有界面预测预测的界面 TM 分数num_recycles1个数值指定AlphaFold为蛋白质序列迭代预测过程的次数ptm表示完整结构的预测 TM 分数ranking_score1个数值一个范围在-100 到 1.5 之间的标量可用于对预测结果进行排序。它将 ptm、iptm、disordered_fraction 和 has_clash 综合为一个数值4mmCIF文件的解析这个文件其实就是我们每次扩散模型采样所生成的一个样本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中文版开源翻译但是版本有点老