GSOC2010 Joao

作者 & 導師

導師

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 的儲存庫、錯誤追蹤系統等）
收集科學文獻並討論一些待實作的方法。

第 1 週 [5 月 31 日 - 6 月 6 日]

重新編號結構的殘基

讀取 SEQRES 記錄以考慮間隙
或者讀取 ATOM 記錄。

探測結構中的二硫鍵

透過 NeighbourSearch 類別
也在標頭中使用 SSBOND

提取生物單元

REMARK350 包含旋轉和平移資訊
如果沒有 REMARK，則不執行任何操作。

第 2 週 [6 月 7 日 – 6 月 13 日]

結構氫化

透過與 WHATIF 伺服器的介面新增所有/極性氫。
可選定義一組 pH 值
- pKa 值演算法

氫化報告

產生結構中極性氫原子的簡短清單。
- 鏈 | 殘基 [編號] | 原子

第 3-5 週 [6 月 14 日 - 7 月 4 日]

移除紊亂的原子

在 Biopython wiki 中提出的解決方案

殘基名稱正規化

建立不同命名法之間的轉換表（在社群連結期間研究它們）
撰寫函式，使給定結構符合給定軟體命名法
- Amber
- CNS/HADDOCK
- GROMACS

粗略晶粒結構

實作降低結構複雜度的函式
- 1pt*c-alpha
- 2pt*c-alpha / c-beta
- 3pt*c-alpha / c-beta / 側鏈偽質心或側鏈質心

第 6 週（期中）[7 月 5 日 - 7 月 11 日]

徹底測試和鞏固功能。
為每個功能撰寫文件和範例，以包含在 Biopython 的 Wiki 和 Bio.PDB 的 FAQ 中。
期中評估。與導師討論專案目前狀態，並調整後續時程表以符合專案需求。

第 7 週 [7 月 12 日 - 7 月 19 日]

新增對 Biopython 的 MODELLER PIR 格式支援

格式描述
SeqIO
AlignIO

允許將結構物件轉換為序列物件

基於 Bio.PDB.Polypeptide 函式

第 8-10 週 [7 月 20 日 - 8 月 9 日]

新增序列/結構同源性函式

建立對 Biopython 的 BLAST 介面的呼叫
- 允許直接從結構物件進行 blast（例如 protein.find_homoseq() ）
- 返回具有 E 值 * 字典（名稱、比對長度等）的元組清單
建立與結構同源性網路服務的介面
- 例如 Dali 伺服器
- 返回具有 Z 分數 * 字典（名稱等）的元組清單

實作基本結構驗證檢查

透過 NeighbourSearch 類別
- 相同電荷接觸
- 原子衝突
透過 ResidueDepth 類別
- 埋藏電荷
介面 WHATIF PDBReport 網路服務
- 剖析 WARNING 和 ERROR 訊息

第 11 週 [8 月 10 日 - 8 月 17 日]

檢閱文件、程式碼，為新函式撰寫測試。

專案程式碼

託管於此 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.'

PDB 檔案的氫化

在導師和我之間的討論之後，我們決定最好不僅為此目的包含一個網路伺服器，還包括一個本機演算法。這不會在使用者缺乏網際網路連線時限制使用者。

已實作 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 物件功能。基本上，結構的所有原子現在都會在建立結構時定義其質量（請查看 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 追蹤 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 物件會從殘基中移除，並新增一個對應於使用者選擇位置的單一 Atom 物件（keep_loc 引數），預設為 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 格式

MODELLER PIR 格式支援已加入 SeqIO，並以 ‘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')