获取一组残差的所有邻居

Question

我有一个保存在centerResidueList = [100, 140, 170, 53] 中的残基编号列表，我正在尝试从这组残基中获取所有相邻的残基。
目前我正在使用下面的脚本，我是否处理整个 PDB 文件并生成距离为 10.0 的原子对列表，然后遍历列表并检查 all_neighbors 列表中的残基编号是否对应于残基编号在centerResidueList。

from Bio.PDB import *

centerResidueList = [100, 140, 170, 53]
neighbours_resi_number = []
structure = PDBParser().get_structure('X', "1xxx.pdb") 
atom_list = Selection.unfold_entities(structure, 'A') 
ns = NeighborSearch(atom_list)
all_neighbors = ns.search_all(10.0, "R") 
for residuepair in all_neighbors:
    resi_number = residuepair[0].id[1]
    if resi_number in centerResidueList:
        resi_number_partner = residuepair[1].id[1]
        neighbours_resi_number.append(resi_number_partner)

首先，我怎样才能只使用 CA 原子创建 atom_list？

其次，residuepair[0].id[1] 是生成残数的正确方法吗（它有效，但有什么方法可以得到）？

最后，有没有更好的解决方案来实现这一点？

Answer 1

使用NeighborSearch 绝对是正确的想法——它构造了一个k-d tree，它可以非常快速地查找最近的邻居。

如果您只有几个残基要搜索，我会在这些残基的原子上使用search() 方法（可能只是为了速度而使用它们的 CA 原子）。这将比使用search_all() 然后过滤更有效。我会回答你的两个问题，然后在底部提供完整的解决方案。

---

如何仅使用 CA 原子创建 atom_list？

您可以使用filter 或列表推导（我认为列表推导更具可读性）：

atom_list = [atom for atom in structure.get_atoms() if atom.name == 'CA']

---

其次，residuepair[0].id[1] 是生成残数的正确方法吗（它有效，但有什么方法可以得到）？

这绝对是正确的方法。 然而（这是一个重要的警告），请注意这不会处理带有insertion codes 的残基。为什么不自己处理Residue 对象呢？

---

我的代码：

from Bio.PDB import NeighborSearch, PDBParser, Selection


structure = PDBParser().get_structure('X', "1xxx.pdb")

chain = structure[0]['A']  # Supply chain name for "center residues"
center_residues = [chain[resi] for resi in [100, 140, 170, 53]]
center_atoms = Selection.unfold_entities(center_residues, 'A')

atom_list = [atom for atom in structure.get_atoms() if atom.name == 'CA']
ns = NeighborSearch(atom_list)

# Set comprehension (Python 2.7+, use `set()` on a generator/list for < 2.7)
nearby_residues = {res for center_atom in center_atoms
                   for res in ns.search(center_atom.coord, 10, 'R')}

# Print just the residue number (WARNING: does not account for icodes)
print sorted(res.id[1] for res in nearby_residues)