Using NeighborSearch
is definitely the right idea - it constructs a k-d tree, which performs very fast lookups on nearest neighbors.
If you have just a few residues to search around, I would use the search()
method on those residues' atoms (perhaps just their CA atoms for speed). This will be more efficient than using search_all()
then filtering. I'll answer your two questions, then provide a complete solution at the bottom.
How can I just create the atom_list using CA atoms only?
You can either use filter
, or a list comprehension (I think list comprehensions are more readable):
atom_list = [atom for atom in structure.get_atoms() if atom.name == 'CA']
Second, is
residuepair[0].id[1]
the correct way of generating the residue numbers (it works but is there a method to get this)?
This is definitely the right approach. However (and this is a significant caveat), note that this will not handle residues with insertion codes. Why not deal with the Residue
objects themselves?
My code:
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)