走向3D:DBC模块

Bio. DBC是一个Biopython模块,专注于处理生物大分子的晶体结构。除其他外,Bio. DBC还包括一个PDBParser类,该类可生成结构对象,该对象可用于以方便的方式访问文件中的原子数据。对解析TSB标头中包含的信息的支持有限。DBC文件格式不再被修改或扩展以支持新内容,PDBx/mminf于2014年成为标准的DBC存档格式。所有全球蛋白质数据库(wwDBC)网站都使用大分子晶体信息文件(mminf)数据字典来描述DBC条目的信息内容。mminf使用灵活且可扩展的关键-值对格式来表示大分子结构数据,并且对单个DBC条目中可以表示的原子、残基或链的数量没有限制(没有拆分条目!)。

读写晶体结构文件

读取mm到岸价格文件

首先创建一个 MMCIFParser 对象:

>>> from Bio.PDB.MMCIFParser import MMCIFParser
>>> parser = MMCIFParser()

然后使用此解析器从mminf文件创建结构对象:

>>> structure = parser.get_structure("1fat", "1fat.cif")

要对mminf文件进行更低级别的访问,您可以使用 MMCIF2Dict 类创建一个Python字典,该字典将mminf文件中的所有mminf标签映射到它们的值。是否有多个值(例如标签的情况 _atom_site.Cartn_y ,它持有 \(y\) 所有原子的坐标)或单个值(例如初始沉积日期),则标签被映射到值列表。该词典是从mminf文件创建的,如下所示:

>>> from Bio.PDB.MMCIF2Dict import MMCIF2Dict
>>> mmcif_dict = MMCIF2Dict("1FAT.cif")

示例:从mminf文件获取溶剂含量:

>>> sc = mmcif_dict["_exptl_crystal.density_percent_sol"]

示例:获取 \(y\) 所有原子的坐标

>>> y_list = mmcif_dict["_atom_site.Cartn_y"]

读取Binaryinf文件

创建一个 BinaryCIFParser 对象:

>>> from Bio.PDB.binary_cif import BinaryCIFParser
>>> parser = BinaryCIFParser()

呼叫 get_structure 带有Binaryinf文件的路径:

>>> parser.get_structure("1GBT", "1gbt.bcif.gz")
<Structure id=1GBT>

读取MMTF格式的文件

您可以使用直接的MMTFParser从文件中读取结构:

>>> from Bio.PDB.mmtf import MMTFParser
>>> structure = MMTFParser.get_structure("PDB/4CUP.mmtf")

或者您可以使用同一个类通过其DBC ID获取结构:

>>> structure = MMTFParser.get_structure_from_url("4CUP")

这为您提供了一个结构对象,就像从TSB或mminf文件中读取一样。

您还可以使用Biopython在内部使用的外部MTF库访问底层数据:

>>> from mmtf import fetch
>>> decoded_data = fetch("4CUP")

例如,您可以仅访问X坐标。

>>> print(decoded_data.x_coord_list)

正在读取DBC文件

首先我们创建一个 PDBParser 对象:

>>> from Bio.PDB.PDBParser import PDBParser
>>> parser = PDBParser(PERMISSIVE=1)

PERMISSIVE 标志指示许多常见问题(请参阅 示例 )将被忽略(但请注意,一些原子和/或残基将丢失)。如果旗帜不存在, PDBConstructionException 如果在解析操作期间检测到任何问题,将生成。

然后,通过让 PDBParser 对象解析DBC文件(在这种情况下称为DBC文件 pdb1fat.ent , 1fat 是用户定义的结构名称):

>>> structure_id = "1fat"
>>> filename = "pdb1fat.ent"
>>> structure = parser.get_structure(structure_id, filename)

You can extract the header and trailer (simple lists of strings) of the PDB file from the PDBParser object with the get_header and get_trailer methods. Note however that many PDB files contain headers with incomplete or erroneous information. Many of the errors have been fixed in the equivalent mmCIF files. Hence, if you are interested in the header information, it is a good idea to extract information from mmCIF files using the ``MMCIF2Dict`` tool described above, instead of parsing the PDB header.

现在澄清了这一点,让我们回到解析DBC标头。结构对象有一个名为 header 这是一个Python字典,将标题记录映射到其值。

例如:

>>> resolution = structure.header["resolution"]
>>> keywords = structure.header["keywords"]

可用的密钥有 name , head , deposition_date , release_date , structure_method , resolution , structure_reference (映射到参考文献列表), journal_reference , author , compound (它映射到包含有关结晶化合物的各种信息的字典), has_missing_residues , missing_residues ,而且 astral (其映射到具有关于域的附加信息(如果存在)的字典)。

has_missing_residues 映射到一个bool,如果至少有一个非空, REMARK 465 找到标题行。在这种情况下,您应该假设实验中使用的分子有一些无法确定ATOM坐标的残基。 missing_residues 映射到包含有关缺失残基信息的字典列表。 The list of missing residues will be empty or incomplete if the PDB header does not follow the template from the PDB specification.

也可以在不创建 Structure 对象,即。直接从DBC文件:

>>> from Bio.PDB import parse_pdb_header
>>> with open(filename, "r") as handle:
...     header_dict = parse_pdb_header(handle)
...

读取PQR文件

为了解析PQR文件,请按照与TSB文件类似的方式继续:

创建一个 PDBParser 对象,使用 is_pqr 旗帜:

>>> from Bio.PDB.PDBParser import PDBParser
>>> pqr_parser = PDBParser(PERMISSIVE=1, is_pqr=True)

is_pqr 标志设置为 True 指示要解析的文件是PQR文件,并且解析器应该读取每个原子条目的原子荷和半径字段。遵循与PQR文件相同的过程,然后生成结构对象,并解析PQR文件。

>>> structure_id = "1fat"
>>> filename = "pdb1fat.ent"
>>> structure = parser.get_structure(structure_id, filename, is_pqr=True)

读取PDBML(PDBML)文件

创建一个 PDBMLParser 对象:

>>> from Bio.PDB.PDBMLParser import PDBMLParser
>>> pdbml_parser = PDBMLParser()

呼叫 get_structure 具有包含PDF格式的DBC结构的文件路径或文件对象:

>>> structure = pdbml_parser.get_structure("1GBT.xml")

正在编写mmCIP文件

MMCIFIO 类可用于将结构写入mminf文件格式:

>>> io = MMCIFIO()
>>> io.set_structure(s)
>>> io.save("out.cif")

Select 类可以以类似的方式使用 PDBIO 下面使用mmCIP词典阅读 MMCIF2Dict 也可以写:

>>> io = MMCIFIO()
>>> io.set_dict(d)
>>> io.save("out.cif")

正在编写DBC文件

使用 PDBIO 为此课程。当然,写出结构的特定部分也很容易。

示例:保存结构

>>> io = PDBIO()
>>> io.set_structure(s)
>>> io.save("out.pdb")

如果您想写出结构的一部分,请使用 Select 类(也在 PDBIO ). Select有四种方法:

  • accept_model(model)

  • accept_chain(chain)

  • accept_residue(residue)

  • accept_atom(atom)

默认情况下,每个方法都返回1(这意味着模型/链/残基/原子包含在输出中)。通过子类化 Select 并在适当时返回0,您可以从输出中排除模型、链等。也许很烦人,但非常强大。下面的代码仅写出了甘氨残基:

>>> class GlySelect(Select):
...     def accept_residue(self, residue):
...         if residue.get_name() == "GLY":
...             return True
...         else:
...             return False
...
>>> io = PDBIO()
>>> io.set_structure(s)
>>> io.save("gly_only.pdb", GlySelect())

如果这对您来说太复杂了, Dice 模块包含一个方便的 extract 写出开始残基和结束残基之间链中所有残基的函数。

编写PQR文件

使用 PDBIO 类,就像对PDB文件一样,使用标志 is_pqr=True . PDBOP方法也可以用于PQR文件。

示例:编写PQR文件

>>> io = PDBIO(is_pqr=True)
>>> io.set_structure(s)
>>> io.save("out.pdb")

编写MTF文件

要将结构写入MTF文件格式:

>>> from Bio.PDB.mmtf import MMTFIO
>>> io = MMTFIO()
>>> io.set_structure(s)
>>> io.save("out.mmtf")

Select 类可以如上所述使用。请注意,标准MMTF文件中包含的绑定信息、二级结构分配和一些其他信息不会写出来,因为不容易从结构对象中确定。此外,在标准MMTF文件中分组为同一实体的分子将被视为单独的实体, MMTFIO .

结构表示

A的总体布局 Structure 对象遵循所谓的SMCRA(结构/模型/链/残差/原子)架构:

  • 结构由模型组成

  • 模型由链条组成

  • 链由残基组成

  • 残基由原子组成

这是许多结构生物学家/生物信息学家思考结构的方式,并提供了一种简单但有效的处理结构的方法。额外的东西本质上是在需要时添加的。的统一建模图 Structure 对象(忘记 Disordered 目前的班级)显示在 图 3 .这样的数据结构不一定最适合结构的大分子含量的表示,但对于良好解释描述该结构的文件(通常是TSB或MMCIF文件)中存在的数据来说,它绝对是必要的。如果此层次结构无法表示结构文件的内容,则相当肯定该文件包含错误,或者至少没有明确描述该结构。如果无法生成SMURA数据结构,则有理由怀疑存在问题。因此,解析DBC文件可用于检测可能的问题。我们将在部分中给出几个例子 示例 .

SMURA体系结构的URL图  ``Structure`` 类

图 3 SMURA体系结构的URL图 Structure

这用于代表大分子结构。带钻石的整线表示聚合,带箭头的整线表示引用,带三角形的整线表示继承,带三角形的虚线表示界面实现。

结构、模型、链和剩余都是实体基本类的子集。Atom类仅(部分)实现Entity接口(因为Atom没有子级)。

对于每个Entity Subclass,您可以通过使用该子类的唯一id作为键来提取该子类(例如,您可以通过使用原子名称字符串作为键来从Residue对象中提取Atom对象,您可以通过使用其链标识符作为键来从Model对象中提取Chain对象)。

无序的原子和残基由DisorderedAtom和DisorderedResidue类别表示,这两个类别都是DisorderedEntityWrapper基本类别的子集。它们隐藏了与无序相关的复杂性,并与Atom和Residue对象完全相同。

一般来说,子实体对象(即Atom、Residue、Chain、Model)可以通过使用id作为密钥从其父对象(即分别为Residue、Chain、Model、Architecture)中提取。

>>> child_entity = parent_entity[child_id]

您还可以获取父实体对象的所有子实体的列表。请注意,此列表以特定方式排序(例如,根据Model对象中Chain对象的链标识符)。

>>> child_list = parent_entity.get_list()

您还可以从孩子那里获取父母:

>>> parent_entity = child_entity.get_parent()

在SMURA层次结构的所有级别,您还可以提取 full id .完整id是一个元组,包含从顶部对象(Structure)到当前对象的所有id。一个Residue对象的完整id例如是这样的:

>>> full_id = residue.get_full_id()
>>> print(full_id)
("1abc", 0, "A", ("", 10, "A"))

这相当于:

  • 带有id的结构 "1abc"

  • 有id的模特 0

  • 带有id的链条 "A"

  • 带有id的剩余物 ("", 10, "A")

残基id表明该残基不是异源残基(也不是水),因为它有一个空白的异源字段,它的序列标识符是10,它的插入代码是 "A" .

要获取实体的id,请使用 get_id 方法:

>>> entity.get_id()

您可以使用 has_id 方法:

>>> entity.has_id(entity_id)

实体的长度等于其子女数量:

>>> nr_children = len(entity)

可以从父实体中删除、重命名、添加等子实体,但这不包括任何健全性检查(例如,可以将具有相同id的两个残基添加到一条链中)。这确实应该通过一个包含完整性检查的漂亮Decorator类来完成,但如果您想使用原始接口,您可以查看代码(recty.py)。

结构

结构对象位于层次结构的顶部。它的id是用户给定的字符串。该结构包含许多模型孩子。大多数晶体结构(但不是全部)包含单个模型,而核磁共振结构通常由多个模型组成。大部分分子的晶体结构混乱也可能导致几种模型。

模型

Model对象的id是一个integer,它源自模型在分析的文件中的位置(它们从0开始自动编号)。晶体结构通常只有一个模型(id为0),而核磁共振文件通常有多个模型。然而许多DBC解析器假设只有一个模型,即 Structure 阶级 Bio.PDB 其设计使得它可以轻松处理具有多个型号的DBC文件。

例如,要从结构对象获取第一个模型,请使用

>>> first_model = structure[0]

Model对象存储Chain子项的列表。

Chain对象的id源自TSB/mminf文件中的链标识符,并且是单个字符(通常是字母)。Model对象中的每个Chain都有一个唯一的id。例如,要从Model对象获取具有标识符“A”的Chain对象,请使用

>>> chain_A = model["A"]

Chain对象存储剩余子项列表。

残基

residue id是一个包含三个元素的数组:

  • hetero-field (hetfield):这是

    • 'W' 对于水分子;

    • 'H_' 随后是其它杂残基的残基名称(例如, 'H_GLC' 在葡萄糖分子的情况下);

    • 标准氨基酸和核酸为空白。

    采用该方案的原因见第节 相关问题 .

  • sequence identifier (resseq),描述链中残基位置的一个积分(例如,100);

  • insertion code (icode);字符串,例如“A”。插入代码有时用于保留某种理想的残基编号方案。Se 80插入突变体(例如插入在Thr 80和Asn 81残基之间)可以例如具有如下序列标识符和插入代码:Thr 80 A、Se 80 B、Asn 81。通过这种方式,残基编号方案与野生型结构保持一致。

因此,上述葡萄糖残基的id将是 (’H_GLC’, 100, ’A’) .如果异源标志和插入代码为空,则可以单独使用序列标识符:

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

异源标记的原因是,许多很多DBC文件对氨基酸和异源残基或水使用相同的序列标识符,如果不使用异源标记,这将产生明显的问题。

毫不奇怪,Residue对象存储了一组Atom子项。它还包含一个字符串,指定残基名称(例如“RST”)和残基的段标识符(X-PLOR用户众所周知,但不用于SMURA数据结构的构建)。

让我们看看一些例子。插入代码为空的Asn 10将具有剩余id (’ ’, 10, ’) .水10将有残留id (’W’, 10, ’) .具有序列标识符10的葡萄糖分子(残基名称为GLC的杂残基)将具有残基id (’H_GLC’, 10, ’) .这样,三个残基(具有相同的插入代码和序列标识符)可以成为同一链的一部分,因为它们的残基id是不同的。

在大多数情况下,hettag和插入代码字段将为空,例如 (’ ’, 10, ’) .在这些情况下,序列标识符可以用作完整id的快捷方式:

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

Chain对象中的每个Residue对象都应该有一个唯一的id。但是,无序的残基是以特殊的方式处理的,如部分所述 点突变 .

Residue对象具有许多附加方法:

>>> 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) 测试残留物是否是氨基酸。

原子

Atom对象存储与原子关联的数据,并且没有子项。原子的id是其原子名称(例如“OG”代表Se残基的侧链氧)。Atom id在剩余物中必须是唯一的。同样,无序原子也是例外,如部分所述 无序原子 .

原子id只是原子名称(例如。 ’CA’ ).在实践中,原子名称是通过删除TSB文件中原子名称中的所有空白来创建的。

然而,在DBC文件中,空间可以是原子名称的一部分。通常,钙原子被称为 ’CA..’ 为了将它们与C区分开来 \(\alpha\) 原子(称为 ’.CA.’ ).如果剥离空间会产生问题(即。两个原子被称为 ’CA’ 在同一残留物中)保持空间。

在DBC文件中,原子名称由4个字符组成,通常带有前置和尾部空白。通常可以删除这些空间以方便使用(例如,氨基酸C \(\alpha\) 原子被标记为“. CA”。在DBC文件中,其中点代表空间)。为了生成原子名称(从而生成原子id),将删除空白,除非这会导致残留中的名称冲突(即两个原子对象具有相同的原子名称和id)。在后一种情况下,尝试包含空间的原子名称。例如,当一个残基包含名称为“. CA”的原子时,可能会发生这种情况。和“CA..”,尽管这不太可能。

存储的原子数据包括原子名称、原子坐标(包括标准差,如果存在)、B因子(包括各向异性B因子和标准差,如果存在)、altloc说明符和完整的原子名称(包括空间)。不存储较少使用的项目,例如有时在TSB文件中指定的原子元素号或原子荷。

要操作原子坐标,请使用 transform 方法 Atom object.使用 set_coord 方法直接指定原子坐标。

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.")

为了表示原子坐标,使用siguaj、各向异性B因子和sigATM Numpy阵列。

get_vector 方法返回一个 Vector 的坐标的对象表示 Atom 对象,允许您对原子坐标进行载体操作。 Vector 实现全套3D向量运算、矩阵相乘(左和右)以及一些高级的旋转相关运算。

As an example of the capabilities of Bio.PDB’s Vector module, suppose that you would like to find the position of a Gly residue’s C\(\beta\) atom, if it had one. Rotating the N atom of the Gly residue along the C\(\alpha\)-C bond over -120 degrees roughly puts it in the position of a virtual C\(\beta\) atom. Here’s how to do it, making use of the rotaxis method (which can be used to construct a rotation around a certain axis) of the Vector module:

# 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

这个例子表明,可以对原子数据进行一些非常重要的载体操作,这可能非常有用。除了所有常见的载体操作(交叉(使用 * * )和点(使用 * )产品、角度、规范等)和上述 rotaxis 函数 Vector 模块还有旋转的方法 (rotmat )或反映 (refmat )一个载体在另一个载体之上。

提取特定 Atom/Residue/Chain/Model 从结构

以下是一些例子:

>>> model = structure[0]
>>> chain = model["A"]
>>> residue = chain[100]
>>> atom = residue["CA"]

请注意,您可以使用快捷方式:

>>> atom = structure[0]["A"][100]["CA"]

障碍

Bio. DBC可以处理无序原子和点突变(即同一位置的Gly和ALA残基)。

一般方法

应该从两个角度来处理混乱:原子和残余的角度。总的来说,我们试图概括由无序引起的所有复杂性。如果您只是想循环所有C \(\alpha\) 原子,你并不关心某些残基是否具有无序的侧链。另一方面,也应该可以在数据结构中完全表示混乱。因此,无序的原子或残基存储在特殊物体中,这些物体的行为就好像没有无序一样。这是通过仅代表无序原子或残基的子集来实现的。选择哪个子集(例如,使用Se残基的两个无序OG侧链原子位置中的哪一个)可以由用户指定。

无序原子

无序原子由普通原子代表 Atom 对象,但全部 Atom 代表同一物理原子的对象存储在 DisorderedAtom 对象(参见 图 3 ).每个 Atom 中对象 DisorderedAtom 可以使用其altloc说明符对对象进行唯一索引。的 DisorderedAtom 对象将所有未捕获的方法调用转发到所选Atom对象,默认情况下是代表占用率最高的原子的对象。用户当然可以更改选择的 Atom 对象,利用其altloc说明符。这样,原子无序就可以正确地表示,而没有太多额外的复杂性。换句话说,如果你对原子无序不感兴趣,你就不会被它困扰。

每个无序原子都有一个特征altloc标识符。您可以指定 DisorderedAtom 对象的行为应该像 Atom 与特定altloc标识符关联的对象:

>>> atom.disordered_select("A")  # select altloc A atom
>>> print(atom.get_altloc())
"A"
>>> atom.disordered_select("B")  # select altloc B atom
>>> print(atom.get_altloc())
"B"

无序残基

常见情况

最常见的情况是包含一个或多个无序原子的残基。这显然可以通过使用DisorderedAtom对象来表示无序原子,并像普通Atom对象一样将DisorderedAtom对象存储在Residue对象中来解决。DisorderedAtom的行为与普通原子(实际上是占用率最高的原子)完全相同,方法是将所有未捕获的方法调用转发到其包含的Atom对象之一(所选Atom对象)。

点突变

当疾病是由于点突变引起的时,即当晶体中存在两个或更多个肽点突变体时,就会出现特殊情况。这方面的一个例子可以在TSB结构1 EN 2中找到。

由于这些残基属于不同的残基类型(例如,假设Se 60和Cys 60),因此它们不应储存在单个残基中 Residue 对象与常见情况一样。在这种情况下,每个残基由一个代表 Residue 对象,以及两者 Residue 对象存储在单个 DisorderedResidue 对象(参见 图 3 ).

DisorderedResidue 对象将所有未捕获的方法转发给所选的 Residue 对象(默认为最后一个 Residue 添加的对象),因此表现得像普通残留物。每个 Residue 中对象 DisorderedResidue 对象可以通过其残基名称唯一识别。在上面的例子中,残基Se 60在 DisorderedResidue 对象,而残基Cys 60将具有id“CYS”。用户可以选择活动的 Residue 中对象 DisorderedResidue 通过此id对象。

示例:假设一条链在第10位有一个点突变,由Se和Cys残基组成。确保该链的残基10表现为Cys残基。

>>> residue = chain[10]
>>> residue.disordered_select("CYS")

此外,您还可以获得所有列表 Atom 对象(即。所有 DisorderedAtom 对象被“拆开”到其个人 Atom 对象)使用 get_unpacked_list 的方法 (Disordered)Residue object.

乙醇残留物

相关问题

异源残基的一个常见问题是,同一链中存在的几个异源和非异源残基共享相同的序列标识符(和插入代码)。因此,为了为每个异源残基生成唯一的id,水和其他异源残基以不同的方式处理。

请记住,Residue对象的id是tuple(hetfield,resseq,icode)。hetfield是空的(“”)代表氨基酸和核酸,而字符串代表水和其他异源残基。hetfield的内容解释如下。

水残留物

残渣的hetfield字符串由字母“W”组成。因此,水的典型残留id是(“W”,1,“”)。

其他异源残基

其他异源残基的hetfield字符串以“H_”开头,后面是残基名称。残基名称为“GLC”的葡萄糖分子将具有hetfield“H_GLC”。例如,其残基id可以是(“H_GLC”,1,“”)。

分析结构

测量距离

原子的负操作符已过载,以返回两个原子之间的距离。

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

测量角度

使用原子坐标的向量表示, calc_angle 职能从 Vector 模块:

>>> vector1 = atom1.get_vector()
>>> vector2 = atom2.get_vector()
>>> vector3 = atom3.get_vector()
>>> angle = calc_angle(vector1, vector2, vector3)

测量扭转角

使用原子坐标的向量表示, calc_dihedral 职能从 Vector 模块:

>>> vector1 = atom1.get_vector()
>>> vector2 = atom2.get_vector()
>>> vector3 = atom3.get_vector()
>>> vector4 = atom4.get_vector()
>>> angle = calc_dihedral(vector1, vector2, vector3, vector4)

内部坐标-距离、角度、扭转角、距离图等

蛋白质结构通常以相对于固定原点的3D XYZ坐标形式提供,如TSB或mminf文件中。的 internal_coords 模块有助于将该系统转换为键长、角和二面角。除了支持标准 psi, phi, chi 等计算蛋白质结构时,这种表示对于平移和旋转是不变的,并且该实现对于结构分析具有多个益处。

首先在这里加载一些模块,以获取后面的示例:

>>> from Bio.PDB.PDBParser import PDBParser
>>> from Bio.PDB.Chain import Chain
>>> from Bio.PDB.internal_coords import *
>>> from Bio.PDB.PICIO import write_PIC, read_PIC, read_PIC_seq
>>> from Bio.PDB.ic_rebuild import write_PDB, IC_duplicate, structure_rebuild_test
>>> from Bio.PDB.SCADIO import write_SCAD
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> from Bio.PDB.PDBIO import PDBIO
>>> import numpy as np

平行二面面、角度和键长

我们从计算结构内部坐标的简单情况开始:

>>> # load a structure as normal, get first chain
>>> parser = PDBParser()
>>> myProtein = parser.get_structure("1a8o", "1A8O.pdb")
>>> myChain = myProtein[0]["A"]
>>> # compute bond lengths, angles, dihedral angles
>>> myChain.atom_to_internal_coordinates(verbose=True)
chain break at THR  186  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at THR  216  due to MaxPeptideBond (1.4 angstroms) exceeded

通过删除 verbose=True 上面的选项。为了避免创建中断并允许不切实际的长N-C键,请重写类变量 MaxPeptideBond ,例如:

>>> IC_Chain.MaxPeptideBond = 4.0
>>> myChain.internal_coord = None  # force re-loading structure data with new cutoff
>>> myChain.atom_to_internal_coordinates(verbose=True)

At this point the values are available at both the chain and residue level. The first residue of 1A8O is HETATM MSE (selenomethionine), so we investigate residue 2 below using either canonical names or atom specifiers. Here we obtain the chi1 dihedral and tau angles by name and by atom sequence, and the C\(\alpha\)-C\(\beta\) distance by specifying the atom pair:

>>> r2 = myChain.child_list[1]
>>> r2
<Residue ASP het=  resseq=152 icode= >
>>> r2ic = r2.internal_coord
>>> print(r2ic, ":", r2ic.pretty_str(), ":", r2ic.rbase, ":", r2ic.lc)
('1a8o', 0, 'A', (' ', 152, ' ')) : ASP  152  : (152, None, 'D') : D
>>> r2chi1 = r2ic.get_angle("chi1")
>>> print(round(r2chi1, 2))
-144.86
>>> r2ic.get_angle("chi1") == r2ic.get_angle("N:CA:CB:CG")
True
>>> print(round(r2ic.get_angle("tau"), 2))
113.45
>>> r2ic.get_angle("tau") == r2ic.get_angle("N:CA:C")
True
>>> print(round(r2ic.get_length("CA:CB"), 2))
1.53

Chain.internal_coord 对象包含面体(3个键合原子)和二面体(4个键合原子)对象的数组和字典。词典按二元组编制索引 AtomKey 对象; AtomKey 对象捕获残基位置、插入代码、1个或3个字符的残基名称、原子名称、altloc和占用。

下面我们得到了同样的 chi1tau 通过索引 Chain 直接使用数组 AtomKey s索引 Chain 数组:

>>> myCic = myChain.internal_coord

>>> r2chi1_object = r2ic.pick_angle("chi1")
>>> # or same thing (as for get_angle() above):
>>> r2chi1_object == r2ic.pick_angle("N:CA:CB:CG")
True
>>> r2chi1_key = r2chi1_object.atomkeys
>>> r2chi1_key  # r2chi1_key is tuple of AtomKeys
(152_D_N, 152_D_CA, 152_D_CB, 152_D_CG)

>>> r2chi1_index = myCic.dihedraNdx[r2chi1_key]
>>> # or same thing:
>>> r2chi1_index == r2chi1_object.ndx
True
>>> print(round(myCic.dihedraAngle[r2chi1_index], 2))
-144.86
>>> # also:
>>> r2chi1_object == myCic.dihedra[r2chi1_key]
True

>>> # hedra angles are similar:
>>> r2tau = r2ic.pick_angle("tau")
>>> print(round(myCic.hedraAngle[r2tau.ndx], 2))
113.45

Chain 级别更复杂(不建议)。如图所示,多个面体将在不同位置共享单个键:

>>> r2CaCb = r2ic.pick_length("CA:CB")  # returns list of hedra containing bond
>>> r2CaCb[0][0].atomkeys
(152_D_CB, 152_D_CA, 152_D_C)
>>> print(round(myCic.hedraL12[r2CaCb[0][0].ndx], 2))  # position 1-2
1.53
>>> r2CaCb[0][1].atomkeys
(152_D_N, 152_D_CA, 152_D_CB)
>>> print(round(myCic.hedraL23[r2CaCb[0][1].ndx], 2))  # position 2-3
1.53
>>> r2CaCb[0][2].atomkeys
(152_D_CA, 152_D_CB, 152_D_CG)
>>> print(round(myCic.hedraL12[r2CaCb[0][2].ndx], 2))  # position 1-2
1.53

请使用 Residue 水平 set_length :math:'函数改为。

测试结构的完整性

缺失原子和其他问题可能会在重建结构时造成问题。使用 structure_rebuild_test :math:`` to determine quickly if a structure has sufficient data for a clean rebuild. Add ` ' verbose=True '和/或检查结果字典以获取更多详细信息:

>>> # check myChain makes sense (can get angles and rebuild same structure)
>>> resultDict = structure_rebuild_test(myChain)
>>> resultDict["pass"]
True

改造和重建结构

It’s preferable to use the residue level set_angle:math:`` and set_length:math:`` facilities for modifying internal coordinates rather than directly accessing the Chain structures. While directly modifying hedra angles is safe, bond lengths appear in multiple overlapping hedra as noted above, and this is handled by set_length:math:. When applied to a dihedral angle, ``set_angle:math:`` will wrap the result to +/-180 and rotate adjacent dihedra as well (such as both bonds for an isoleucine chi1 angle - which is probably what you want).

>>> # rotate residue 2 chi1 angle by -120 degrees
>>> r2ic.set_angle("chi1", r2chi1 - 120.0)
>>> print(round(r2ic.get_angle("chi1"), 2))
95.14
>>> r2ic.set_length("CA:CB", 1.49)
>>> print(round(myCic.hedraL12[r2CaCb[0][0].ndx], 2))  # Cb-Ca-C position 1-2
1.49

从内部坐标重建结构是一个简单的调用 internal_to_atom_coordinates() :

>>> myChain.internal_to_atom_coordinates()

>>> # just for proof:
>>> myChain.internal_coord = None  # all internal_coord data removed, only atoms left
>>> myChain.atom_to_internal_coordinates()  # re-generate internal coordinates
>>> r2ic = myChain.child_list[1].internal_coord
>>> print(round(r2ic.get_angle("chi1"), 2))  # show measured values match what was set above
95.14
>>> print(round(myCic.hedraL23[r2CaCb[0][1].ndx], 2))  # N-Ca-Cb position 2-3
1.49

生成的结构可以像正常情况一样用PD贝洛编写:

write_PDB(myProtein, "myChain.pdb")
# or just the ATOM records without headers:
io = PDBIO()
io.set_structure(myProtein)
io.save("myChain2.pdb")

蛋白质内部坐标(.pic)文件和默认值

文件格式在 PICIO 将蛋白质链描述为相对于初始坐标的面体和两面体的模块。除了残基序列信息之外的文件的所有部分(例如, (’1A8O’, 0, ’A’, (’ ’, 153, ’)) ILE )是可选的,如果未指定,将填写默认值,并且 read_PIC :math:`` is called with the ` ' defaults=True '选项。默认值从2019年9月Dunbrack cullpdb_pc20_res2.2_R1.0开始计算。

在这里,我们将“myChain”写成 .pic 内部坐标规范文件,然后将其读回为“myProtein 2”。

# write chain as 'protein internal coordinates' (.pic) file
write_PIC(myProtein, "myChain.pic")
# read .pic file
myProtein2 = read_PIC("myChain.pic")

由于所有内部坐标值都可以用默认值替换, PICIO.read_PIC_seq :math:'作为实用函数提供,用于从输入序列创建有效(主要是螺旋形)默认结构:

# create default structure for random sequence by reading as .pic file
myProtein3 = read_PIC_seq(
    SeqRecord(
        Seq("GAVLIMFPSTCNQYWDEHKR"),
        id="1RND",
        description="my random sequence",
    )
)
myProtein3.internal_to_atom_coordinates()
write_PDB(myProtein3, "myRandom.pdb")

探索例如所需的准确性可能会很有趣 omega 角(180.0)、多面体角和/或键长。picFlags选项用于 write_PIC :math:'启用此功能,允许将数据选择写入.pic文件,而不是保留未指定以获取默认值。

可以进行各种组合,并且提供一些预设,例如 classic 只会写 psi, phi, tau 、Pro omega 和侧链 chi 指向.pic文件:

write_PIC(myProtein, "myChain.pic", picFlags=IC_Residue.pic_flags.classic)
myProtein2 = read_PIC("myChain.pic", defaults=True)

初始化全原子AtomArray

Biopython中的所有3D XYZ坐标 Atom 对象被移动到 Chain 类,并在早期步骤中由Numpy“视图”到该数组中取代 atom_to_internal_coordinates :math:``. Software accessing Biopython ` “Atom”坐标不受影响,但新阵列可能会为未来的工作提供效率。

不像 Atom XYZ坐标, AtomArray 坐标是齐性的,这意味着它们是像这样的数组 [ x y z 1.0] 第四个元素是1.0。这促进了在整个 internal_coords module.有相应的 AtomArrayIndex 字典、地图 AtomKeys 他们的坐标。

在这里,我们演示了特定C \(\beta\) 从数组中删除原子,然后表明修改数组值会修改 Atom 对象同时:

>>> # access the array of all atoms for the chain, e.g. r2 above is residue 152 C-beta
>>> r2_cBeta_index = myChain.internal_coord.atomArrayIndex[AtomKey("152_D_CB")]
>>> r2_cBeta_coords = myChain.internal_coord.atomArray[r2_cBeta_index]
>>> print(np.round(r2_cBeta_coords, 2))
[-0.75 -1.18 -0.51  1.  ]

>>> # the Biopython Atom coord array is now a view into atomArray, so
>>> assert r2_cBeta_coords[1] == r2["CB"].coord[1]
>>> r2_cBeta_coords[1] += 1.0  # change the Y coord 1 angstrom
>>> assert r2_cBeta_coords[1] == r2["CB"].coord[1]
>>> # they are always the same (they share the same memory)
>>> r2_cBeta_coords[1] -= 1.0  # restore

请注意,很容易“破坏”Atom coord数组和链atomArray之间的视图链接。当直接修改Atom坐标时,请使用逐个元素复制的语法来避免这种情况:

# use these:
myAtom1.coord[:] = myAtom2.coord
myAtom1.coord[...] = myAtom2.coord
myAtom1.coord[:] = [1, 2, 3]
for i in range(3):
    myAtom1.coord[i] = myAtom2.coord[i]

# do not use:
myAtom1.coord = myAtom2.coord
myAtom1.coord = [1, 2, 3]

使用 atomArrayIndex 和知识 AtomKey 类使我们能够创建Numpy“选择器”,如下所示,以仅提取C的数组 \(\alpha\) 原子坐标:

>>> # create a selector to filter just the C-alpha atoms from the all atom array
>>> atmNameNdx = AtomKey.fields.atm
>>> aaI = myChain.internal_coord.atomArrayIndex
>>> CaSelect = [aaI.get(k) for k in aaI.keys() if k.akl[atmNameNdx] == "CA"]
>>> # now the ordered array of C-alpha atom coordinates is:
>>> CA_coords = myChain.internal_coord.atomArray[CaSelect]
>>> # note this uses Numpy fancy indexing, so CA_coords is a new copy
>>> # (if you modify it, the original atomArray is unaffected)

距离图

的益处 atomArray 从它生成距离图是一行 Numpy 代码:

np.linalg.norm(atomArray[:, None, :] - atomArray[None, :, :], axis=-1)

尽管它很复杂,但这个习语很难记住,并且在上面的形式中生成全原子距离,而不是经典的C \(\alpha\) 可以根据需要绘制。的 distance_plot :math:`` method wraps the line above and accepts an optional selector like ` “Caselect”在上一节中定义。看到 图 4 .

# create a C-alpha distance plot
caDistances = myChain.internal_coord.distance_plot(CaSelect)
# display with e.g. MatPlotLib:
import matplotlib.pyplot as plt

plt.imshow(caDistances, cmap="hot", interpolation="nearest")
plt.show()
DBC文件1A8 O的C阿尔法距离图

图 4 DBC文件1A 8 O(HIV病毒外壳C末端结构域)的C-阿尔法距离图

从远处建造结构

全原子距离图是蛋白质结构的另一种表示,也不受平移和旋转的影响,但缺乏手征信息(镜像结构将生成相同的距离图)。通过将距离矩阵与每个二面角的符号相结合,可以重新生成内部坐标。

这项工作使用了由Hedronmeter Blue开发的方程,在https://math.stackexchange.com/a/49340/409和http://daylateanddollarshort.com/mathdocs/Heron-like-Results-for-Tetrahedral-Volume.pdf上进行了讨论。

首先,我们从“myChain”中提取距离和手征值:

>>> ## create the all-atom distance plot
>>> distances = myCic.distance_plot()
>>> ## get the signs of the dihedral angles
>>> chirality = myCic.dihedral_signs()

我们需要匹配“myChain”的有效数据结构才能正确重建它;使用 read_PIC_seq :math:'上面在一般情况下适用,但这里使用的1A 8 O示例具有一定的ALTSYS复杂性,这是序列本身不会产生的。对于演示来说,最简单的方法是简单地复制“myChain”结构,但我们将所有原子和内部坐标链数组设置为0(仅用于演示),只是为了确保没有数据来自原始结构:

>>> ## get new, empty data structure : copy data structure from myChain
>>> myChain2 = IC_duplicate(myChain)[0]["A"]
>>> cic2 = myChain2.internal_coord

>>> ## clear the new atomArray and di/hedra value arrays, just for proof
>>> cic2.atomArray = np.zeros((cic2.AAsiz, 4), dtype=np.float64)
>>> cic2.dihedraAngle[:] = 0.0
>>> cic2.hedraAngle[:] = 0.0
>>> cic2.hedraL12[:] = 0.0
>>> cic2.hedraL23[:] = 0.0

该方法是根据距离图数据重新生成内部坐标,然后根据内部坐标生成原子坐标,如上所示。为了将最终生成的结构放置在与初始结构相同的坐标空间中,我们仅复制前三个N-C的坐标 \(\alpha\) - 从“myChain”的链起点到“myChain 2”结构的C原子(这仅需要证明结尾的等效性):

>>> ## copy just the first N-Ca-C coords so structures will superimpose:
>>> cic2.copy_initNCaCs(myChain.internal_coord)

distance_to_internal_coordinates :math:`` routine needs arrays of the six inter-atom distances for each dihedron for the target structure. The convenience routine ` distplot_to_dh_arrays` :math:`` extracts these values from the previously generated distance matrix as needed, and may be replaced by a user method to write these data to the arrays in the ` `Chain.internal_coords``对象

>>> ## copy distances to chain arrays:
>>> cic2.distplot_to_dh_arrays(distances, chirality)
>>> ## compute angles and dihedral angles from distances:
>>> cic2.distance_to_internal_coordinates()

以下步骤从新生成的“myChain 2”内部坐标生成原子坐标,然后使用Numpy allclose :math:'用于确认所有值是否匹配优于TSB文件分辨率的例行程序:

>>> ## generate XYZ coordinates from internal coordinates:
>>> myChain2.internal_to_atom_coordinates()
>>> ## confirm result atomArray matches original structure:
>>> np.allclose(cic2.atomArray, myCic.atomArray)
True

请注意,此过程不使用整个距离矩阵,而仅使用每个二面角的四个原子之间的六个局部距离。

叠加残留物及其邻近区域

internal_coords 模块依赖于不同坐标空间之间的原子坐标转换,以计算扭转角和重建结构。每个二面体都有一个坐标空间变换,将其第一个原子放置在XZ平面上,第二个原子放置在原点,第三个原子放置在+Z轴上,以及相应的反向变换,将其返回到原始结构中的坐标。这些变换矩阵可供使用,如下所示。通过明智地选择参考二面体,可以在多个蛋白质结构中研究和可视化成对和更高级的残基相互反应,例如 图 5 .

DBC文件3PBL中的邻近苯色素侧链

图 5 DBC文件3PBL中的邻近苯色素侧链(人多巴胺D3受体)

This example superimposes each PHE residue in a chain on its N-C\(\alpha\)-C\(\beta\) atoms, and presents all PHEs in the chain in the respective coordinate space as a simple demonstration. A more realistic exploration of pairwise sidechain interactions would examine a dataset of structures and filter for interaction classes as discussed in the relevant literature.

# superimpose all phe-phe pairs - quick hack just to demonstrate concept
# for analyzing pairwise residue interactions.  Generates PDB ATOM records
# placing each PHE at origin and showing all other PHEs in environment

## shorthand for key variables:
cic = myChain.internal_coord
resNameNdx = AtomKey.fields.resname
aaNdx = cic.atomArrayIndex

## select just PHE atoms:
pheAtomSelect = [aaNdx.get(k) for k in aaNdx.keys() if k.akl[resNameNdx] == "F"]
aaF = cic.atomArray[pheAtomSelect]  # numpy fancy indexing makes COPY not view

for ric in cic.ordered_aa_ic_list:  # internal_coords version of get_residues()
    if ric.lc == "F":  # if PHE, get transform matrices for chi1 dihedral
        chi1 = ric.pick_angle("chi1")  # N:CA:CB:CG space has C-alpha at origin
        cst = np.transpose(chi1.cst)  # transform TO chi1 space
        # rcst = np.transpose(chi1.rcst)  # transform FROM chi1 space (not needed here)
        cic.atomArray[pheAtomSelect] = aaF.dot(cst)  # transform just the PHEs
        for res in myChain.get_residues():  # print PHEs in new coordinate space
            if res.resname in ["PHE"]:
                print(res.internal_coord.pdb_residue_string())
        cic.atomArray[pheAtomSelect] = aaF  # restore coordinate space from copy

3D打印蛋白质结构

OpenSCAD(https:openscad.org)是一种用于创建实体3D CAD对象的语言。OpenSCAD中提供了根据内部坐标构建蛋白质结构的算法,并提供了描述结构的数据,以便可以生成适合3D打印的模型。虽然其他软件可以生成RTL数据作为3D打印的渲染选项(例如Chimera,https://www.cgl.ucsf.edu/chimera/),但这种方法生成球体和圆柱体作为输出,因此更适合与3D打印蛋白质结构相关的修改。可以在OpenSCAD代码中选择单个残基和键进行特殊处理,例如通过尺寸突出显示或在特定位置添加可旋转键(例如,请参阅https://www.thingiverse.com/thing:3957471)。

# write OpenSCAD program of spheres and cylinders to 3d print myChain backbone
## set atom load filter to accept backbone only:
IC_Residue.accept_atoms = IC_Residue.accept_backbone
## set chain break cutoff very high to bridge missing residues with long bonds
IC_Chain.MaxPeptideBond = 4.0
## delete existing data to force re-read of all atoms with attributes set above:
myChain.internal_coord = None
write_SCAD(myChain, "myChain.scad", scale=10.0)

internal_coords 控制属性

中提供了一些控制属性 internal_coords 类用于在计算内部坐标时修改或过滤数据。这些都列在表中  Bio.PDB.internal_coords中的控制属性。 :

表 3 Bio.PDB.internal_coords中的控制属性。

属性

默认

效果

AtomKey

D2h

Convert D atoms to H if True

IC_Chain

MaxPeptideBond

1.4

不存在链断裂的最大C-N长度;增大以连接3D模型的缺失残基

IC_Residue

accept_atoms

主线,氢原子

覆盖以删除部分或所有侧链、H、D

IC_Residue

accept_resnames

CYG、YCM、不详

3-要处理的HEATM的字母名称,仅限于主干,除非添加到ic_data.py

IC_Residue

gly_Cbeta

覆盖以生成Gly C \(\beta\) 基于数据库平均值的原子

确定原子与原子的接触

使用 NeighborSearch 以执行邻居查找。邻居查找是使用用C编写的KD树模块完成的(请参阅 KDTree 模块中的类 Bio.PDB.kdtrees ),使其速度非常快。它还包括一种快速方法来查找彼此一定距离内的所有点对。

两个结构叠加

使用 Superimposer 对象初始化两个坐标集。该对象计算旋转和平移矩阵,该矩阵将两个原子列表相互旋转,以使其RMSD最小化。当然,这两个列表需要包含相同数量的原子。的 Superimposer 对象还可以将旋转/平移应用于原子列表。旋转和平移作为一个数组存储在 rotran 属性 Superimposer 对象(注意旋转是右乘!)。RMSD存储在 rmsd 属性

使用的算法 Superimposer 来自Golub & Van Loan [Golub1989] 并利用奇异值分解(这是在一般 Bio.SVDSuperimposer 模块)。

例如:

>>> sup = Superimposer()
# Specify the atom lists
# 'fixed' and 'moving' are lists of Atom objects
# The moving atoms will be put on the fixed atoms
>>> sup.set_atoms(fixed, moving)
# Print rotation/translation/rmsd
>>> print(sup.rotran)
>>> print(sup.rms)
# Apply rotation/translation to the moving atoms
>>> sup.apply(moving)

要根据其活性位点对两种结构进行量化,请使用活性位点原子来计算旋转/平移矩阵(如上所述),并将其应用于整个分子。

计算半球曝光

半球暴露(HSE)是一种新的2D溶剂暴露测量方法 [Hamelryck2005]. 基本上,它统计了C的数量 \(\alpha\) 残基周围的原子沿其侧链方向和相反方向(半径在13 GHz内)。尽管它很简单,但它的表现优于许多其他溶剂暴露指标。

HSE comes in two flavors: HSE\(\alpha\) and HSE\(\beta\). The former only uses the C\(\alpha\) atom positions, while the latter uses the C\(\alpha\) and C\(\beta\) atom positions. The HSE measure is calculated by the HSExposure class, which can also calculate the contact number. The latter class has methods which return dictionaries that map a Residue object to its corresponding HSE\(\alpha\), HSE\(\beta\) and contact number values.

例如:

>>> model = structure[0]
>>> hse = HSExposure()
# Calculate HSEalpha
>>> exp_ca = hse.calc_hs_exposure(model, option="CA3")
# Calculate HSEbeta
>>> exp_cb = hse.calc_hs_exposure(model, option="CB")
# Calculate classical coordination number
>>> exp_fs = hse.calc_fs_exposure(model)
# Print HSEalpha for a residue
>>> print(exp_ca[some_residue])

确定二级结构

对于此功能,您需要安装DS SP(并获得其许可证-免费用于学术用途,请参阅https://swift.cmbi.umpn.nl/gv/dssp/)。然后使用 DSSP 类,它映射 Residue 对象其二级结构(和可及表面积)。表中列出了DS SP代码  Bio.PDB中的DS SP代码。 .注意DSSP(程序,因此类)不能处理多个模型!

表 4 Bio.PDB中的DS SP代码。

代码

二级结构

H

\(\alpha\) - 螺旋

B

分离 \(\beta\) - 桥残基

E

G

3-10螺旋

I

\(\Pi\) - 螺旋

T

反过来

S

弯曲

其他

DSSP 类还可用于计算残留物的可及表面积。但另请参阅部分 计算残留深度 .

计算残留深度

Residue depth is the average distance of a residue’s atoms from the solvent accessible surface. It’s a fairly new and very powerful parameterization of solvent accessibility. For this functionality, you need to install Michel Sanner’s MSMS program (https://www.scripps.edu/sanner/html/msms_home.html). Then use the ResidueDepth class. This class behaves as a dictionary which maps Residue objects to corresponding (residue depth, C\(\alpha\) depth) tuples. The C\(\alpha\) depth is the distance of a residue’s C\(\alpha\) atom to the solvent accessible surface.

例如:

>>> model = structure[0]
>>> rd = ResidueDepth(model, pdb_file)
>>> residue_depth, ca_depth = rd[some_residue]

您还可以访问分子表面本身(通过 get_surface 函数),采用具有表面点的数字Python数组的形式。

DBC文件中的常见问题

众所周知,许多DBC文件包含语义错误(不是结构本身,而是它们在DBC文件中的表示)。Bio. ZB尝试通过两种方式处理这个问题。PDBParser对象可以以两种方式运行:限制性方式和许可性方式(默认方式)。

例如:

# Permissive parser
>>> parser = PDBParser(PERMISSIVE=1)
>>> parser = PDBParser()  # The same (default)
# Strict parser
>>> strict_parser = PDBParser(PERMISSIVE=0)

在许可状态(DID)下,明显包含错误的DBC文件将被“更正”(即省略了一些残基或原子)。这些错误包括:

  • 具有相同标识符的多个残基

  • 具有相同标识符的多个原子(考虑altloc标识符)

这些错误表明DBC文件中存在真正的问题(详细信息,请参阅Hamelryck和Manderick,2003 [Hamelryck2003A]) .在限制状态下,出现错误的DBC文件会导致发生异常。这对于查找TSB文件中的错误很有用。

然而,有些错误会自动纠正。通常,每个无序原子都应该有一个非空白的altloc标识符。然而,有许多结构不遵循这一惯例,并且对于同一原子的两个无序位置有空白和非空白标识符。这会自动以正确的方式解释。

有时,一个结构包含一系列属于链A的残基,然后是属于链B的残基,然后是属于链A的残基,即链被“断裂”。这也是正确的解释。

示例

PDBParser/Architecture类别在大约800个结构上进行了测试(每个结构属于一个独特的SCOP超家族)。这大约需要20分钟,即每个结构平均1.5秒。在1000 MHz PC上解析包含约64000个原子的大核糖体亚基(1FKK)的结构需要10秒。

在无法构建明确的数据结构的情况下,会生成三个异常。在这三种情况下,可能的原因是DBC文件中的错误,应该更正。在这些情况下生成异常比在数据结构中错误描述结构的机会要好得多。

重复残基

一种结构在一条链中包含两个氨基酸残基,具有相同的序列标识符(resseq 3)和icode。检查后发现该链含有残留物Thr A3、..、Gly A202、Leu A3、Glu A204。显然,Leu A3应该是Leu A203。结构1FFK存在一些类似的情况(例如,其含有Gly B64、Met B65、Glu B65、Thr B67,即残基Glu B65应该是Glu B66)。

重复原子

结构1 EJG在A链第22位含有Se/Pro点突变。反过来,第22号序列包含一些无序原子。正如预期的那样,属于22号序列的所有原子都有一个非空白altloc指定符(B或C)。Pro 22的所有原子都具有altloc A,除了N原子具有空白altloc。这会产生一个例外,因为属于点突变处两个残基的所有原子都应该具有非空白altloc。事实证明,这个原子可能由Se和Pro 22共享,因为Se 22缺少N原子。这再次指出了文件中的一个问题:N原子应该同时存在于Se和Pro残基中,在这两种情况下都与合适的altloc标识符相关。

自动校正

有些错误非常常见,并且可以很容易地纠正,而没有太大做出错误解释的风险。这些案例如下。

无序原子的空白altloc

通常,每个无序原子都应该有一个非空白的altloc标识符。然而,有许多结构不遵循这一惯例,并且对于同一原子的两个无序位置有空白和非空白标识符。这会自动以正确的方式解释。

断裂链

有时,一个结构包含一系列属于链A的残基,然后是属于链B的残基,然后是属于链A的残基,即链被“断裂”。这是正确的解释。

致命错误

有时DBC文件无法明确解释。不会猜测并冒出错的风险,而是会生成异常,并且用户需要更正DBC文件。这些案例如下。

重复残基

链中的所有残基都应该有一个唯一的id。该id是基于:

  • 序列标识符(resseq)。

  • 插入代码(icode)。

  • hetfield字符串(“W”代表waters,“H_”后面跟着残基名称代表其他异源残基)

  • 点突变情况下残基的残基名称(将Residue对象存储在DisorderedResidue对象中)。

如果这并没有导致一个唯一的id,那么很可能是错误的,并生成一个异常。

重复原子

残基中的所有原子都应该有一个唯一的id。该id是基于:

  • 原子名称(不带空位,或者如果出现问题则带空位)。

  • altloc说明符。

如果这并没有导致一个唯一的id,那么很可能是错误的,并生成一个异常。

初始化蛋白质数据库

从蛋白质数据库下载结构

可以使用 retrieve_pdb_file 法在 PDBList object.此方法的参数是结构的DBC标识符。

>>> pdbl = PDBList()
>>> pdbl.retrieve_pdb_file("1FAT")

PDBList 类还可以用作命令行工具:

python PDBList.py 1fat

下载的文件将被调用 pdb1fat.ent 并存储在当前工作目录中。注意到 retrieve_pdb_file 方法还有一个可选参数 pdir 它指定存储下载的DBC文件的特定目录。

retrieve_pdb_file 方法还有一些选项来指定用于下载的压缩格式以及用于本地解压缩的程序(默认 .Z 格式和 gunzip ).此外,还可以在创建 PDBList object.默认情况下,使用全球蛋白质数据库(ftp://ftp.wwpdb.org/pub/pdb/data/structures/divided/pdb/)的服务器。有关更多详细信息,请参阅API文档。再次感谢克里斯蒂安·罗瑟捐赠这个模块。

下载整个DBC

以下命令将所有DBC文件存储在 /data/pdb 目录:

python PDBList.py all /data/pdb

python PDBList.py all /data/pdb -d

为此,调用API方法 download_entire_pdb .添加 -d 选项将所有文件存储在同一目录中。否则,它们将根据其DBC ID分类到DBC风格的子目录中。根据流量不同,完整下载需要2-4天。

保持TSB的本地副本最新

这也可以使用 PDBList object.人们只是创建一个 PDBList 对象(指定存在DBC本地副本的目录)并调用 update_pdb 方法:

>>> pl = PDBList(pdb="/data/pdb")
>>> pl.update_pdb()

当然,人们可以每周 cronjob 以保持本地副本自动更新。也可以指定PDB ftp站点(参见API文档)。

PDBList 有一些可以使用的额外方法。的 get_all_obsolete 方法可用于获取所有过时的DBC条目的列表。的 changed_this_week 方法可用于获取本周添加、修改或废弃的条目。有关以下可能性的更多信息 PDBList 请参阅API文档。

一般问题

Bio. DBC的测试结果如何?

实际上,相当不错。Bio. DBC已在来自TSB的近5500个结构上进行了广泛测试-所有结构似乎都得到了正确解析。更多详细信息请参阅Bio. DBC Bioinformatics文章。Bio. DBC已经/正在许多研究项目中用作可靠的工具。事实上,我几乎每天都使用Bio. DBC进行研究,并继续致力于改进它并添加新功能。

它有多快?

PDBParser 在大约800个结构上进行了性能测试(每个结构属于一个独特的SCOP超家族)。这大约需要20分钟,即每个结构平均1.5秒。在1000 MHz PC上解析包含约64000个原子的大核糖体亚基(1FKK)的结构需要10秒。简而言之:对于许多应用程序来说,它已经足够快了。

是否支持分子图形?

不是直接的,主要是因为已经有相当多的基于Python/Python感知的解决方案,可以与Bio. DBC一起使用。我的选择是Pyol,顺便说一句(我已经成功地将其与Bio. DBC一起使用,并且很快/有一天Bio. DBC中可能会有特定的PyMol模块)。基于Python/感知的分子图形解决方案包括:

谁在使用Bio. DBC?

PDB被用于DISEMBL的构建,DISEMBL是一个预测蛋白质中无序区域的网络服务器(http://dis.embl.de/)。PDB还被用于在PDB中进行蛋白质结构之间的活性位点相似性的大规模搜索 [Hamelryck2003B], 并开发识别线性二级结构元素的新算法 [Majumdar2005].

从对功能和信息的请求来看,Bio. DBC也被多家LCC(大型制药公司:-)使用。