如何计算具有多个 models/conformations 的蛋白质的平均结构
How to calculate the average structure of a protein with multiple models/conformations
我有一个 PDB 文件 '1abz' (https://files.rcsb.org/view/1ABZ.pdb),其中包含具有 23 种不同模型(编号为 MODEL 1-23)的蛋白质结构的坐标。请忽略 header 注释,有趣的信息从第 276 行开始 'MODEL 1'。
我想计算蛋白质的平均结构。蛋白质的 PDB 文件包含多个 conformations/models,我想计算每个残基的单个原子的平均坐标,这样我最终得到一个 conformation/model。
我不知道如何使用 Biopython 做到这一点,所以我尝试使用 Pandas 来计算平均坐标。我想我已经设法计算出平均值,但现在的问题是我有一个不再是 PDB 格式的 csv 文件,所以我无法将此文件加载到 PyMol 中。
我的问题是,如何将我的 csv 文件转换为 PDB 格式。更好的是,如何在不影响原始 pdb 文件格式的情况下获得 Biopython 或 Python 中的平均坐标?
这是我用来计算 pandas 中的平均坐标的代码。
#I first converted my pdb file to a csv file
import pandas as pd
import re
pdbfile = '1abz.pdb'
df = pd.DataFrame(columns=['Model','Residue','Seq','Atom','x','y','z']) #make dataframe object
i = 0 #counter
b = re.compile("MODEL\s+(\d+)")
regex1 = "([A-Z]+)\s+(\d+)\s+([^\s]+)\s+([A-Z]+)[+-]?\s+([A-Z]|)"
regex2 = "\s+(\d+)\s+([+-]?\d+\.\d+\s+[+-]?\d+\.\d+\s+[+-]?\d+\.\d+)"
reg = re.compile(regex1+regex2)
with open(pdbfile) as f:
columns = ('label', 'ident', 'atomName', 'residue', 'chain', 'sequence', 'x', 'y', 'z', 'occ', 'temp', 'element')
data = []
for line in f:
n = b.match(line)
if n:
modelNum = n.group(1)
m = reg.match(line)
if m:
d = dict(zip(columns, line.split()))
d['model'] = int(modelNum)
data.append(d)
df = pd.DataFrame(data)
df.to_csv(pdbfile[:-3]+'csv', header=True, sep='\t', mode='w')
#Then I calculated the average coordinates
df = pd.read_csv('1abz.csv', delim_whitespace = True, usecols = [0,5,7,8,10,11,12])
df1 = df.groupby(['atomName','residue','sequence'],as_index=False)['x','y','z'].mean()
df1.to_csv('avg_coord.csv', header=True, sep='\t', mode='w')
这在biopython中当然是可行的。让我用一个忽略 pdb 文件中的 HETRES 的例子来帮助你:
首先用你所有的模型解析 pdb 文件:
import Bio.PDB
import numpy as np
parser = Bio.PDB.PDBParser(QUIET=True) # Don't show me warnings
structure = parser.get_structure('1abz', '1abz.pdb') # id of pdb file and location
好的,现在我们有了文件的内容,假设你的所有模型中都有相同的原子,得到一个列表,每个原子都有一个唯一的标识符(例如:chain + residue pos + atom name) :
atoms = [a.parent.parent.id + '-' + str(a.parent.id[1]) + '-' + a.name for a in structure[0].get_atoms() if a.parent.id[0] == ' '] # obtained from model '0'
请注意,我忽略了带有 a.parent.id[0] == ' '
的杂残基。现在让我们得到每个原子的平均值:
atom_avgs = {}
for atom in atoms:
atom_avgs[atom] = []
for model in structure:
atom_ = atom.split('-')
coor = model[atom_[0]][int(atom_[1])][atom_[2]].coord
atom_avgs[atom].append(coor)
atom_avgs[atom] = sum(atom_avgs[atom]) / len(atom_avgs[atom]) # average
现在让我们使用结构的一个模型创建新的 pdb:
ns = Bio.PDB.StructureBuilder.Structure('id=1baz') # new structure
ns.add(structure[0]) # add model 0
for atom in ns[0].get_atoms():
chain = atom.parent.parent
res = atom.parent
if res.id[0] != ' ':
chain.detach_child(res) # detach hetres
else:
coor = atom_avgs[chain.id + '-' + str(res.id[1]) + '-' + atom.name]
atom.coord = coor
现在让我们编写 pdb
io = Bio.PDB.PDBIO()
io.set_structure(ns)
io.save('new_1abz.pdb')
我有一个 PDB 文件 '1abz' (https://files.rcsb.org/view/1ABZ.pdb),其中包含具有 23 种不同模型(编号为 MODEL 1-23)的蛋白质结构的坐标。请忽略 header 注释,有趣的信息从第 276 行开始 'MODEL 1'。
我想计算蛋白质的平均结构。蛋白质的 PDB 文件包含多个 conformations/models,我想计算每个残基的单个原子的平均坐标,这样我最终得到一个 conformation/model。
我不知道如何使用 Biopython 做到这一点,所以我尝试使用 Pandas 来计算平均坐标。我想我已经设法计算出平均值,但现在的问题是我有一个不再是 PDB 格式的 csv 文件,所以我无法将此文件加载到 PyMol 中。
我的问题是,如何将我的 csv 文件转换为 PDB 格式。更好的是,如何在不影响原始 pdb 文件格式的情况下获得 Biopython 或 Python 中的平均坐标?
这是我用来计算 pandas 中的平均坐标的代码。
#I first converted my pdb file to a csv file
import pandas as pd
import re
pdbfile = '1abz.pdb'
df = pd.DataFrame(columns=['Model','Residue','Seq','Atom','x','y','z']) #make dataframe object
i = 0 #counter
b = re.compile("MODEL\s+(\d+)")
regex1 = "([A-Z]+)\s+(\d+)\s+([^\s]+)\s+([A-Z]+)[+-]?\s+([A-Z]|)"
regex2 = "\s+(\d+)\s+([+-]?\d+\.\d+\s+[+-]?\d+\.\d+\s+[+-]?\d+\.\d+)"
reg = re.compile(regex1+regex2)
with open(pdbfile) as f:
columns = ('label', 'ident', 'atomName', 'residue', 'chain', 'sequence', 'x', 'y', 'z', 'occ', 'temp', 'element')
data = []
for line in f:
n = b.match(line)
if n:
modelNum = n.group(1)
m = reg.match(line)
if m:
d = dict(zip(columns, line.split()))
d['model'] = int(modelNum)
data.append(d)
df = pd.DataFrame(data)
df.to_csv(pdbfile[:-3]+'csv', header=True, sep='\t', mode='w')
#Then I calculated the average coordinates
df = pd.read_csv('1abz.csv', delim_whitespace = True, usecols = [0,5,7,8,10,11,12])
df1 = df.groupby(['atomName','residue','sequence'],as_index=False)['x','y','z'].mean()
df1.to_csv('avg_coord.csv', header=True, sep='\t', mode='w')
这在biopython中当然是可行的。让我用一个忽略 pdb 文件中的 HETRES 的例子来帮助你:
首先用你所有的模型解析 pdb 文件:
import Bio.PDB
import numpy as np
parser = Bio.PDB.PDBParser(QUIET=True) # Don't show me warnings
structure = parser.get_structure('1abz', '1abz.pdb') # id of pdb file and location
好的,现在我们有了文件的内容,假设你的所有模型中都有相同的原子,得到一个列表,每个原子都有一个唯一的标识符(例如:chain + residue pos + atom name) :
atoms = [a.parent.parent.id + '-' + str(a.parent.id[1]) + '-' + a.name for a in structure[0].get_atoms() if a.parent.id[0] == ' '] # obtained from model '0'
请注意,我忽略了带有 a.parent.id[0] == ' '
的杂残基。现在让我们得到每个原子的平均值:
atom_avgs = {}
for atom in atoms:
atom_avgs[atom] = []
for model in structure:
atom_ = atom.split('-')
coor = model[atom_[0]][int(atom_[1])][atom_[2]].coord
atom_avgs[atom].append(coor)
atom_avgs[atom] = sum(atom_avgs[atom]) / len(atom_avgs[atom]) # average
现在让我们使用结构的一个模型创建新的 pdb:
ns = Bio.PDB.StructureBuilder.Structure('id=1baz') # new structure
ns.add(structure[0]) # add model 0
for atom in ns[0].get_atoms():
chain = atom.parent.parent
res = atom.parent
if res.id[0] != ' ':
chain.detach_child(res) # detach hetres
else:
coor = atom_avgs[chain.id + '-' + str(res.id[1]) + '-' + atom.name]
atom.coord = coor
现在让我们编写 pdb
io = Bio.PDB.PDBIO()
io.set_structure(ns)
io.save('new_1abz.pdb')