João Rodrigues [email protected]
导师
Eric Talevich
Diana Jaunzeikare
Peter Cock
Biopython 是生物信息学和计算生物学中非常流行的库。它的 Bio.PDB
模块最初由 Thomas Hamelryck 开发,是结构生物学家的一个简单而强大的工具。尽管它提供可靠的 PDB 解析器功能,并允许对大分子进行多种计算(最近邻搜索、RMS),但它仍然缺少许多研究人员日常工作中需要用到的功能。在模块的巧妙结构和良好的整体组织下,可以将探测结构中的二硫键并相应添加极性氢原子等功能整合到 Bio.PDB
中。链移除和残基重命名(为了适应不同的现有命名法)以及重新编号等美观操作也将会受到社区的热烈欢迎。
可以改进 Bio.PDB
的另一个方面是为大分子模拟中的重量级软件(例如 MODELLER、GROMACS、AutoDock、HADDOCK)提供一个平滑的集成/交互层。有人可能会认为最简单的解决方案是为这些软件包的功能和例程编写挂钩。但是,我认为,像最近开发的 edPDB 或更完整的 Biskit 库 这样的项目使得这些接口工作变得多余。相反,我认为将这些软件的输入/输出格式包含在 Biopython 的 SeqIO
和 AlignIO
模块中更有利。这与为模型验证/结构检查服务/软件创建接口一起,将允许 Biopython 用作模拟前和模拟后的工具。最终,它将为其包含在结构建模、分子动力学和对接模拟的管道和工作流程中铺平道路。
以下进度表被安排得灵活,这意味着一些功能可能会提前完成。此外,周次包含了对功能的文档和单元测试工作,并在项目的两个点(中途,最后一周)对这些工作进行更长时间的审查。
熟悉开发环境(GitHub 帐户、Git、Biopython 的代码库、错误跟踪系统等)
收集科学文献,并讨论一些待实现的方法。
NeighbourSearch
类彻底测试和巩固功能。
为每个功能编写文档和示例,以包含在 Biopython 的 Wiki 和 Bio.PDB
的 FAQ 中。
中期评估。与导师讨论项目的当前状态,并调整后续计划以符合项目的需要。
SeqIO
AlignIO
Bio.PDB.Polypeptide
函数protein.find_homoseq()
)托管在 此 GitHub 分支 上
由于我添加了一些仅对蛋白质有用的/逻辑的方法,因此将它们公开在 Structure.py
中以适用于每种分子可能会产生误导。我们决定添加一个 as_protein()
方法,允许访问特定于蛋白质的方法。以下示例演示了此调用的工作方式。请注意,search_ss_bonds
方法在 dir(s)
中不存在,但在 dir(prot)
中存在。
>>> from Bio.PDB import PDBParser
>>> p = PDBParser()
>>> s = p.get_structure("example", "4PTI.pdb")
>>> dir(s) # Cut for viewing purposes
['__doc__', ... , 'renumber_residues', 'set_parent', 'xtra']
>>> prot = s.as_protein()
>>> dir(prot)
['__doc__', ... , 'renumber_residues', 'search_ss_bonds', 'set_parent', 'xtra']
由于 parse_pdb_header
远非最佳,并且很可能在将来发生变化,因此我选择放弃读取 SEQREQ 记录来说明间隙。但是,忽略此信息并根据 ATOM 记录进行重新编号会导致我们丢失关于间隙的信息。我选择从所有残基中减去第一个残基编号-1,从而使编号从 1 开始,并且仍然保留间隙。我还添加了一个参数(start),允许用户设置计数的起始编号。
示例
from Bio.PDB import PDBParser
p = PDBParser()
s = p.get_structure("example", "1IHM.pdb")
print(list(s.get_residues())[0])
# <Residue ASP het= resseq=1029 icode= >
s.renumber_residues()
print(list(s.get_residues())[0])
# <Residue ASP het= resseq=1 icode= >
SEQRES 中的相同原理也适用于排除查找 SSBOND。此外,我并没有使用 NeighborSearch
来查找键距内的半胱氨酸对,而是使用了减号运算符,因为它已被重载以返回两个原子之间的距离(FAQ 的第 10 页)。文献中引用的平均距离是 2.05A,但其他软件包和我的测试将 3.0A 设置为一个良好的阈值。但是,用户可以手动设置自己的阈值。
该函数返回一个包含残基对元组的迭代器。
from Bio.PDB import PDBParser
p = PDBParser()
s = p.get_structure("example", "4PTI.pdb")
prot = s.as_protein()
for bond in prot.search_ss_bonds():
print(bond)
# (<Residue CYS het= resseq=5 icode= >, <Residue CYS het= resseq=55 icode= >)
# (<Residue CYS het= resseq=14 icode= >, <Residue CYS het= resseq=38 icode= >)
# (<Residue CYS het= resseq=30 icode= >, <Residue CYS het= resseq=51 icode= >)
在 parse_pdb_header
中添加了对 REMARK350 的解析,因为已经为另一个 REMARK 部分编写了一些代码。这将提取标题中的变换矩阵和平移向量,然后将其馈送到 Structure
函数中。每个新的旋转结构都作为新的 MODEL 创建。我选择这样做是因为晶体结构很少有多个 MODEL 实例,而且因为 NMR 模型没有 REMARK 350(至少据我所知)。
from Bio.PDB import PDBParser
p = PDBParser()
s1 = p.get_structure("a", "4PTI.pdb")
s1.build_biological_unit()
# 'Processed 0 transformations on the structure.' # Identity matrix is ignored.
s2 = p.get_structure("b", "homol_1bd8.pdb") # A homology model
s2.build_biological_unit()
# 'PDB File lacks appropriate REMARK 350 entries to build Biological Unit.'
s3 = p.get_structure("c", "1IHM.pdb")
s3.build_biological_unit()
# 'Processed 59 transformations on the structure.'
在与导师的讨论后,我们决定,也许最好不仅包括一个用于此目的的 Web 服务器,而且还包括一个本地算法。这样一来,当用户没有互联网连接时,不会限制用户。
WHATIF 质子化服务的接口已经实现,但目前应视为高度实验性。与该服务器的接口包括为类似 PDBXML 的格式编写一个小的解析器,该解析器预计其初始版本中存在严重错误。我运行了一些简单的测试,它可以工作。它尚不支持水分子,也不支持除蛋白质以外的任何其他分子。这些问题有望在以后解决。
对于那些有勇气想要测试它(并帮助我调试它)的人,这里是一个示例用法。
from Bio.Struct.WWW import WHATIF
from Bio import Struct
server = (
WHATIF.WHATIF()
) # Performs a sort of PING to the server. Gracefully exits if the servers are down.
# Get the protein structure
structure = Struct.read("4PTI.pdb")
protein = structure.as_protein() # This excludes water molecules
# Upload the structure to the WHATIF server
# This should convert the structure from a Structure object to a string via tempfile and PDBIO
# I was having some issues uploading structures...
id = server.UploadPDB(protein)
# Protonate
# Returns a Structure Object / WARNING! Bug prone for now.
protein_h = server.PDBasXMLwithSymwithPolarH(id)
关于本地实现,经过大量的阅读,我决定使用 PyMol 的算法。它似乎允许对任何结构进行质子化,无论其性质如何(蛋白质、DNA 等)。它的矢量和矩阵操作很可能可以使用 Numpy 和 Biopython 的 Vector
模块进行优化。此第一个实现仅适用于蛋白质。我将在以后添加通用分子支持。
from Bio import Struct
from Bio.Struct import Hydrogenate as H
s = Struct.read("1ctf.pdb")
p = s.as_protein()
prot = H.Hydrogenate_Protein()
prot.add_hydrogens(p)
质心函数最初是作为新模块 Bio.Struct.Geometry
的一部分开发的。它允许计算几何中心(所有质量相等)和质心(考虑原子的元素质量)。质量是原子对象的新功能,源自 此列表 和 PyMol。本质上,在创建结构时,结构中的所有原子都会定义其质量(请检查 Atom.py
和 此线程 以了解详细信息)。这显然是实验性的。
要计算任何实体(结构、模型、链、残基)或原子列表的质心
>>> from Bio.Struct.Geometry import center_of_mass
>>> from Bio import Struct
>>> s = Struct.read('4PTI.pdb')
>>> print(center_of_mass.__doc__)
Returns gravitic or geometric center of mass of an Entity.
Geometric assumes all masses are equal (geometric=True)
Defaults to Gravitic.
>>> print(center_of_mass(s))
[14.833301303933874, 21.431581746366263, 4.1218478418007134]
>>> print(center_of_mass(s, geometric=True))
[14.805324902127458, 21.365571977563405, 4.1108949403803985]
截至目前,支持 3 种 CG 模型。
1) CA-Trace 2) ENCAD 3 点模型 (CA、O、侧链珠) 3) MARTINI 蛋白质模型 (BB、侧链点 [S1 到 S4])
一个示例,从上面的 s 结构中获取
>>> p = s.as_protein() # To expose the CG method
>>> ca_trace = p.coarse_grain()
>>> # One atom per residue
>>> print(len(list(p.get_residues())) == len(list(ca_trace.get_atoms())))
True
>>> cg_encad = p.coarse_grain('ENCAD_3P')
>>> for residue in cg_encad.get_residues():
... print(residue.resname, residue.child_list)
...
ARG [<Atom CA>, <Atom O>, <Atom CMA>]
PRO [<Atom CA>, <Atom O>, <Atom CMA>]
ASP [<Atom CA>, <Atom O>, <Atom CMA>]
PHE [<Atom CA>, <Atom O>, <Atom CMA>]
CYS [<Atom CA>, <Atom O>, <Atom CMA>]
LEU [<Atom CA>, <Atom O>, <Atom CMA>]
GLU [<Atom CA>, <Atom O>, <Atom CMA>]
PRO [<Atom CA>, <Atom O>, <Atom CMA>]
PRO [<Atom CA>, <Atom O>, <Atom CMA>]
TYR [<Atom CA>, <Atom O>, <Atom CMA>]
...
CYS [<Atom CA>, <Atom O>, <Atom CMA>]
GLY [<Atom CA>, <Atom O>]
GLY [<Atom CA>, <Atom O>]
ALA [<Atom CA>, <Atom O>, <Atom CMA>]
>>> cg_martini = p.coarse_grain('MARTINI')
>>> for residue in cg_martini.get_residues():
... print(residue.resname, residue.child_list)
...
ARG [<Atom BB>, <Atom S1>, <Atom S2>]
PRO [<Atom BB>, <Atom S1>]
ASP [<Atom BB>, <Atom S1>]
PHE [<Atom BB>, <Atom S1>, <Atom S2>, <Atom S3>]
CYS [<Atom BB>, <Atom S1>]
LEU [<Atom BB>, <Atom S1>]
GLU [<Atom BB>, <Atom S1>]
PRO [<Atom BB>, <Atom S1>]
PRO [<Atom BB>, <Atom S1>]
TYR [<Atom BB>, <Atom S1>, <Atom S2>, <Atom S3>]
......
CYS [<Atom BB>, <Atom S1>]
GLY [<Atom BB>]
GLY [<Atom BB>]
ALA [<Atom BB>]
作为 Structure.py
的一部分实现,并且松散地基于 Ramon Crehuet 的贡献。DisorderedAtom
对象从残基中移除,并且添加一个与用户选择的位置(keep_loc
参数)相对应的单个 Atom
对象,默认值为 A。
一个示例,仍然保留上面的 s
>>> s = s.remove_disordered_atoms(verbose=True)
0 residues were modified
>>> # Now if we load a structure with disordered atoms
>>> ds = Struct.read('1MC2.pdb')
>>> ds.remove_disordered_atoms(verbose=True)
Residue TRP:1010 has 8 disordered atoms: CD1/CD2/NE1/CE2/CE3/CZ2/CZ3/CH2
Residue VAL:1018 has 3 disordered atoms: CB/CG1/CG2
Residue LEU:1024 has 4 disordered atoms: CB/CG/CD1/CD2
Residue ARG:1043 has 7 disordered atoms: CB/CG/CD/NE/CZ/NH1/NH2
Residue MET:1092 has 4 disordered atoms: CB/CG/SD/CE
Residue ARG:1107 has 7 disordered atoms: CB/CG/CD/NE/CZ/NH1/NH2
Residue GLU:1108 has 4 disordered atoms: CG/CD/OE1/OE2
Residue ASP:1111 has 4 disordered atoms: CB/CG/OD1/OD2
Residue SER:1116 has 1 disordered atoms: OG
Residue SER:1131 has 1 disordered atoms: O
10 residues were modified
Biopython 支持 BLAST(通过 NCBI 服务器进行本地和远程搜索)。我们桥接了 Bio.PDB
和 Bio.Blast
模块,以便更容易地搜索序列同源物。目前,它支持通过 Bio.Blast.NCBIWWW
进行远程 BLAST,并且作为一个黑盒工作 - 即用户无法更改任何搜索参数。如果用户想完全使用 BLAST,则应使用常规的 BLAST
模块。这只是一个便利函数。
该功能仅对 Protein
对象开放。它使用结构对象的序列查询 NCBI BLAST 服务器的 PDB 子集数据库,并针对短序列(少于 15 个残基)自动调整参数。
它返回一个按期望值排序的列表,其中包含一些信息值(e 值、同一性、正值、间隙)、匹配的 PDB 代码和比对结果。
>>> from Bio import Struct
>>> s = Struct.read('1A8O.pdb')
>>> p = s.as_protein()
>>> seq_homologues = p.find_seq_homologues()
>>> for homologues in seq_homologues:
... print(homologues[0], homologues[1])
... print(homologues[-1])
... print()
...
2BUO 1.82482e-31
DIRQGPKEPFRDYVDRFYKTLRAEQASQEVKNW-TETLLVQNANPDCKTILKALGPGATLEE--TACQG
DIRQGPKEPFRDYVDRFYKTLRAEQASQEVKNW TETLLVQNANPDCKTILKALGPGATLEE TACQG
DIRQGPKEPFRDYVDRFYKTLRAEQASQEVKNWMTETLLVQNANPDCKTILKALGPGATLEEMMTACQG
1AUM 1.82482e-31
DIRQGPKEPFRDYVDRFYKTLRAEQASQEVKNW-TETLLVQNANPDCKTILKALGPGATLEE--TACQG
DIRQGPKEPFRDYVDRFYKTLRAEQASQEVKNW TETLLVQNANPDCKTILKALGPGATLEE TACQG
DIRQGPKEPFRDYVDRFYKTLRAEQASQEVKNWMTETLLVQNANPDCKTILKALGPGATLEEMMTACQG
...
对 SeqIO
添加了 MODELLER PIR 格式支持,称为 'pir-modeller'。目前,该格式可以读取但不能写入。以下是一个格式示例以及解析器使用示例。
>P1;5fd1
structureX:5fd1:1 :A:106 :A:ferredoxin:Azotobacter vinelandii: 1.90: 0.19
AFVVTDNCIKCKYTDCVEVCPVDCFYEGPNFLVIHPDECIDCALCEPECPAQAIFSEDEVPEDMQEFIQLNAELA
EVWPNITEKKDPLPDAEDWDGVKGKLQHLER*
>>> from Bio import SeqIO
>>> for i in SeqIO.parse("test_pir.txt", "pir-modeller"):
... print(i)
...
ID: 5fd1
Name: 5fd1
Description: ferredoxin
Number of features: 0
/r_factor= 0.19
/end_residue=106
/initial_chain=a
/end_chain=a
/record_type=X-Ray Structure
/initial_residue=1
/resolution= 1.90
/source_organism=Azotobacter vinelandii
Seq('AFVVTDNCIKCKYTDCVEVCPVDCFYEGPNFLVIHPDECIDCALCEPECPAQAI...LER')