-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathprotein_interface.py
198 lines (156 loc) · 8.23 KB
/
protein_interface.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
import mdtraj as md
import numpy as np
import itertools
from topology_objects import Residue, Atom, store_residue
from utils import resid_from_aidx, ang, index_pairs, ion_index_pairs, aromatic_cation_index_pairs
from hbonds import find_hbond_triplets
from aromatics import aromatic_centroid_orth
class Interface:
'''
Holds information about the protein-protein (or protein-ligand) interface.
Attributes
----------
t : md.Trajectory
Input trajectory.
sel_receptor : str
Selection string for receptor.
sel_ligand : str
Selection string for ligand.
'''
def __init__(self, t, sel_receptor, sel_ligand):
self.t = t
print(f'Analyzing interactions in trajectory with {t.n_frames} frames and {t.n_atoms} atoms.')
# atom indices of all receptor and ligand atoms
self.idx_receptor = t.top.select(sel_receptor)
self.idx_ligand = t.top.select(sel_ligand)
#residue indices of all receptor and ligand residues
self.resid_receptor = resid_from_aidx(self.t, self.idx_receptor)
self.resid_ligand = resid_from_aidx(self.t, self.idx_ligand)
print(f'Receptor has {len(self.idx_receptor)} atoms and {len(self.resid_receptor)} residues.')
print(f'Ligand has {len(self.idx_ligand)} atoms and {len(self.resid_ligand)} residues.')
# residue objects for the interface
interface_resid_receptor, interface_resid_ligand = self.get_interface_residue_idx(cutoff=0.5)
self.interface_residues_receptor = [store_residue(self.t, x) for x in interface_resid_receptor]
self.interface_residues_ligand = [store_residue(self.t, x) for x in interface_resid_ligand]
# atom objects for atoms in interface residues
interface_atoms_receptor = []
interface_atoms_ligand = []
for r in self.interface_residues_receptor:
for a in r.atoms:
interface_atoms_receptor.append(a)
for r in self.interface_residues_ligand:
for a in r.atoms:
interface_atoms_ligand.append(a)
polar_interface_atoms_receptor = [a for a in interface_atoms_receptor if a.is_polar]
polar_interface_atoms_ligand = [a for a in interface_atoms_ligand if a.is_polar]
acceptor_interface_atoms_receptor = [a for a in interface_atoms_receptor if a.is_hbond_acceptor]
acceptor_interface_atoms_ligand = [a for a in interface_atoms_ligand if a.is_hbond_acceptor]
donor_interface_atoms_receptor = [a for a in interface_atoms_receptor if a.is_hbond_donor]
donor_interface_atoms_ligand = [a for a in interface_atoms_ligand if a.is_hbond_donor]
cation_interface_atoms_receptor = [a for a in interface_atoms_receptor if a.is_cation]
cation_interface_atoms_ligand = [a for a in interface_atoms_ligand if a.is_cation]
anion_interface_atoms_receptor = [a for a in interface_atoms_receptor if a.is_anion]
anion_interface_atoms_ligand = [a for a in interface_atoms_ligand if a.is_anion]
print('Interface:')
print(f'Receptor: {len(self.interface_residues_receptor)} residues, {len(interface_atoms_receptor)} atoms.')
print(f'Ligand: {len(self.interface_residues_ligand)} residues, {len(interface_atoms_ligand)} atoms.')
# receptor-ligand pairs (and triplets for hbonds)
self.residue_pairs = index_pairs(self.interface_residues_receptor, self.interface_residues_ligand)
self.aromatic_residue_pairs = index_pairs(self.interface_residues_receptor, self.interface_residues_ligand, aromatic_residues_only=True)
self.polar_atom_pairs = index_pairs(polar_interface_atoms_receptor, polar_interface_atoms_ligand)
self.ion_atom_pairs = ion_index_pairs(cat_rec=cation_interface_atoms_receptor, an_rec=anion_interface_atoms_receptor, cat_lig=cation_interface_atoms_ligand, an_lig=anion_interface_atoms_ligand)
self.hbond_triplets = find_hbond_triplets(rec_donors=donor_interface_atoms_receptor, rec_acceptors=acceptor_interface_atoms_receptor, lig_donors=donor_interface_atoms_ligand, lig_acceptors=acceptor_interface_atoms_ligand)
self.aromatic_cation_pairs = aromatic_cation_index_pairs(rec=self.interface_residues_receptor, cat_rec=cation_interface_atoms_receptor, lig=self.interface_residues_ligand, cat_lig=cation_interface_atoms_ligand)
# contacts
self.residue_contacts = None
self.get_residue_contacts()
self.polar_atom_contacts = self.get_atom_contacts(self.polar_atom_pairs, cutoff=0.37)
self.ionic_contacts = self.get_atom_contacts(self.ion_atom_pairs, cutoff=0.37)
# h-bonds
self.hbonds = self.get_hbonds(self.hbond_triplets)
# aromatic interactions
self.aromatic_cation = self.get_cation_pi(self.aromatic_cation_pairs, cutoff=0.45)
def get_interface_residue_idx(self, cutoff=0.5, frame=0):
'''
Get the residue indices of the interface residues based on contact distance between ligand and receptor residues.
Parameters
----------
cutoff : float
Cutoff distance.
frame : int
Only one frame of the trajectory is considered. Default 0.
Returns
-------
interface_resid_receptor : list of int
Residue indices of receptor residues in the interface.
interface_resid_ligand : list of int
Residue indices of ligand indices in the interface.
'''
interface_resid_receptor = []
interface_resid_ligand = []
# compute residue contacts of all ligand residue pairs
residue_pairs = np.array([x for x in itertools.product(self.resid_ligand, self.resid_receptor)])
contacts, _ = md.compute_contacts(self.t[frame], residue_pairs, scheme='closest-heavy')
interface_pairs = np.where(contacts[0] < cutoff)[0]
for i in interface_pairs:
rp = residue_pairs[i]
if rp[0] not in interface_resid_ligand:
interface_resid_ligand.append(rp[0])
if rp[1] not in interface_resid_receptor:
interface_resid_receptor.append(rp[1])
return interface_resid_receptor, interface_resid_ligand
def get_residue_contacts(self, cutoff=0.35):
'''
Calculate residue contacts, where the closest heavy atoms of residues of ligand and receptor come closer than the cutoff distance.
Parameters
----------
cutoff : float
Distance cutoff in nm.
'''
residue_pairs = self.residue_pairs
contacts, _ = md.compute_contacts(self.t, residue_pairs, scheme='closest-heavy')
self.residue_contacts = contacts < cutoff
def get_atom_contacts(self, pairs, cutoff=0.37):
'''
Calculate contacts between atoms.
Parameters
----------
pairs : list of tuple of int
Atom pairs.
cutoff : float
Distance cutoff in nm.
'''
dists = md.compute_distances(self.t, pairs)
return dists < cutoff
def get_hbonds(self, triplets):
'''
Calculate H bonds based on D-H-A triplets and the Baker Hubbard criterion.
'''
HA_cutoff = 0.25
DHA_cutoff = 2.09
# HA pairs for distance calculation
ha_pairs = [x[1:3] for x in triplets]
# distance and angle calculations
dists = md.compute_distances(self.t, ha_pairs)
angs = md.compute_angles(self.t, triplets)
# boolean masks from distance and angles
mask_dist = dists < HA_cutoff
mask_ang = angs < DHA_cutoff
mask = mask_dist * mask_ang
return mask
def get_cation_pi(self, pairs, cutoff=0.45):
'''
Calculate cation-pi interactions.
'''
contacts = np.zeros((len(self.t), len(pairs)), dtype=bool)
for i, p in enumerate(pairs):
c_ar, o_ar = aromatic_centroid_orth(self.t, p[0])
c_cat = self.t.xyz[:, p[1]]
d = np.linalg.norm(c_ar - c_cat, axis=1)
contacts[:, i] = d < cutoff
return contacts
if __name__ == '__main__':
t = md.load('test/test.xtc', top='test/test.pdb')
sel_receptor = 'resid 0 to 222'
sel_ligand = 'resid 223 to 280'
I = Interface(t, sel_receptor=sel_receptor, sel_ligand=sel_ligand)