OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
contact_score.py
Go to the documentation of this file.
1 import itertools
2 import numpy as np
3 
4 import time
5 from ost import mol
6 from ost import geom
7 from ost import io
8 
10  """ Helper object for Contact-score computation
11  """
12  def __init__(self, ent, contact_d = 5.0, contact_mode="aa"):
13 
14  if contact_mode not in ["aa", "repr"]:
15  raise RuntimeError("contact_mode must be in [\"aa\", \"repr\"]")
16 
17  if contact_mode == "repr":
18  for r in ent.residues:
19  repr_at = None
20  if r.IsPeptideLinking():
21  cb = r.FindAtom("CB")
22  if cb.IsValid():
23  repr_at = cb
24  elif r.GetName() == "GLY":
25  ca = r.FindAtom("CA")
26  if ca.IsValid():
27  repr_at = ca
28  elif r.IsNucleotideLinking():
29  c3 = r.FindAtom("C3'")
30  if c3.IsValid():
31  repr_at = c3
32  else:
33  raise RuntimeError(f"Only support peptide and nucleotide "
34  f"residues in \"repr\" contact mode. "
35  f"Problematic residue: {r}")
36  if repr_at is None:
37  raise RuntimeError(f"Residue {r} has no required "
38  f"representative atom (CB for peptide "
39  f"residues (CA for GLY) C3' for "
40  f"nucleotide residues.")
41 
42  self._contact_mode = contact_mode
43 
44  if self.contact_mode == "aa":
45  self._view = ent.CreateFullView()
46  elif self.contact_mode == "repr":
47  pep_query = "(peptide=true and (aname=\"CB\" or (rname=\"GLY\" and aname=\"CA\")))"
48  nuc_query = "(nucleotide=True and aname=\"C3'\")"
49  self._view = ent.Select(" or ".join([pep_query, nuc_query]))
50  self._contact_d = contact_d
51 
52  # the following attributes will be lazily evaluated
53  self._chain_names = None
54  self._interacting_chains = None
55  self._sequence = dict()
56  self._contacts = None
57  self._hr_contacts = None
58  self._interface_residues = None
59  self._hr_interface_residues = None
60 
61  @property
62  def view(self):
63  """ The structure depending on *contact_mode*
64 
65  Full view in case of "aa", view that only contains representative
66  atoms in case of "repr".
67 
68  :type: :class:`ost.mol.EntityView`
69  """
70  return self._view
71 
72  @property
73  def contact_mode(self):
74  """ The contact mode
75 
76  Can either be "aa", meaning that all atoms are considered to identify
77  contacts, or "repr" which only considers distances between
78  representative atoms. For peptides thats CB (CA for GLY), for
79  nucleotides thats C3'.
80 
81  :type: :class:`str`
82  """
83  return self._contact_mode
84 
85  @property
86  def contact_d(self):
87  """ Pairwise distance of residues to be considered as contacts
88 
89  Given at :class:`ContactScorer` construction
90 
91  :type: :class:`float`
92  """
93  return self._contact_d
94 
95  @property
96  def chain_names(self):
97  """ Chain names in :attr:`~view`
98 
99  Names are sorted
100 
101  :type: :class:`list` of :class:`str`
102  """
103  if self._chain_names is None:
104  self._chain_names = sorted([ch.name for ch in self.view.chains])
105  return self._chain_names
106 
107  @property
109  """ Pairs of chains in :attr:`~view` with at least one contact
110 
111  :type: :class:`list` of :class:`tuples`
112  """
113  if self._interacting_chains is None:
114  self._interacting_chains = list(self.contacts.keys())
115  return self._interacting_chains
116 
117  @property
118  def contacts(self):
119  """ Interchain contacts
120 
121  Organized as :class:`dict` with key (cname1, cname2) and values being
122  a set of tuples with the respective residue indices.
123  cname1 < cname2 evaluates to True.
124  """
125  if self._contacts is None:
126  self._SetupContacts()
127  return self._contacts
128 
129  @property
130  def hr_contacts(self):
131  """ Human readable interchain contacts
132 
133  Human readable version of :attr:`~contacts`. Simple list with tuples
134  containing two strings specifying the residues in contact. Format:
135  <cname>.<rnum>.<ins_code>
136  """
137  if self._hr_contacts is None:
138  self._SetupContacts()
139  return self._hr_contacts
140 
141  @property
143  """ Interface residues
144 
145  Residues in each chain that are in contact with any other chain.
146  Organized as :class:`dict` with key cname and values the respective
147  residue indices in a :class:`set`.
148  """
149  if self._interface_residues is None:
151  return self._interface_residues
152 
153  @property
155  """ Human readable interface residues
156 
157  Human readable version of :attr:`interface_residues`. :class:`list` of
158  strings specifying the interface residues in format:
159  <cname>.<rnum>.<ins_code>
160  """
161  if self._interface_residues is None:
163  return self._hr_interface_residues
164 
165  def GetChain(self, chain_name):
166  """ Get chain by name
167 
168  :param chain_name: Chain in :attr:`~view`
169  :type chain_name: :class:`str`
170  """
171  chain = self.view.FindChain(chain_name)
172  if not chain.IsValid():
173  raise RuntimeError(f"view has no chain named \"{chain_name}\"")
174  return chain
175 
176  def GetSequence(self, chain_name):
177  """ Get sequence of chain
178 
179  Returns sequence of specified chain as raw :class:`str`
180 
181  :param chain_name: Chain in :attr:`~view`
182  :type chain_name: :class:`str`
183  """
184  if chain_name not in self._sequence:
185  ch = self.GetChain(chain_name)
186  s = ''.join([r.one_letter_code for r in ch.residues])
187  self._sequence[chain_name] = s
188  return self._sequence[chain_name]
189 
190  def _SetupContacts(self):
191  # this function is incredibly inefficient... if performance is an issue,
192  # go ahead and optimize
193  self._contacts = dict()
194  self._hr_contacts = list()
195 
196  # set indices relative to full view
197  for ch in self.view.chains:
198  for r_idx, r in enumerate(ch.residues):
199  r.SetIntProp("contact_idx", r_idx)
200 
201  for cname in self.chain_names:
202  # q1 selects stuff in current chain that is close to any other chain
203  q1 = f"cname={cname} and {self.contact_d} <> [cname!={cname}]"
204  # q2 selects stuff in other chains that is close to current chain
205  q2 = f"cname!={cname} and {self.contact_d} <> [cname={cname}]"
206  v1 = self.view.Select(q1)
207  v2 = self.view.Select(q2)
208  v1_p = [geom.Vec3List([a.pos for a in r.atoms]) for r in v1.residues]
209  for r1, p1 in zip(v1.residues, v1_p):
210  for ch2 in v2.chains:
211  cname2 = ch2.GetName()
212  if cname2 > cname:
213  v2_p = [geom.Vec3List([a.pos for a in r.atoms]) for r in ch2.residues]
214  for r2, p2 in zip(ch2.residues, v2_p):
215  if p1.IsWithin(p2, self.contact_d):
216  cname_key = (cname, cname2)
217  if cname_key not in self._contacts:
218  self._contacts[cname_key] = set()
219  self._contacts[cname_key].add((r1.GetIntProp("contact_idx"),
220  r2.GetIntProp("contact_idx")))
221  rnum1 = r1.GetNumber()
222  hr1 = f"{cname}.{rnum1.num}.{rnum1.ins_code}"
223  rnum2 = r2.GetNumber()
224  hr2 = f"{cname2}.{rnum2.num}.{rnum2.ins_code}"
225  self._hr_contacts.append((hr1.strip("\u0000"),
226  hr2.strip("\u0000")))
227 
228  def _SetupInterfaceResidues(self):
229  self._interface_residues = {cname: set() for cname in self.chain_names}
230  for k,v in self.contacts.items():
231  for item in v:
232  self._interface_residues[k[0]].add(item[0])
233  self._interface_residues[k[1]].add(item[1])
234 
235  def _SetupHRInterfaceResidues(self):
236  interface_residues = set()
237  for item in self.hr_contacts:
238  interface_residues.add(item[0])
239  interface_residues.add(item[1])
240  self._hr_interface_residues = list(interface_residues)
241 
242 
244  """
245  Holds data relevant to compute ics
246  """
247  def __init__(self, n_trg_contacts, n_mdl_contacts, n_union, n_intersection):
248  self._n_trg_contacts = n_trg_contacts
249  self._n_mdl_contacts = n_mdl_contacts
250  self._n_union = n_union
251  self._n_intersection = n_intersection
252 
253  @property
254  def n_trg_contacts(self):
255  """ Number of contacts in target
256 
257  :type: :class:`int`
258  """
259  return self._n_trg_contacts
260 
261  @property
262  def n_mdl_contacts(self):
263  """ Number of contacts in model
264 
265  :type: :class:`int`
266  """
267  return self._n_mdl_contacts
268 
269  @property
270  def precision(self):
271  """ Precision of model contacts
272 
273  The fraction of model contacts that are also present in target
274 
275  :type: :class:`int`
276  """
277  if self._n_mdl_contacts != 0:
278  return self._n_intersection / self._n_mdl_contacts
279  else:
280  return 0.0
281 
282  @property
283  def recall(self):
284  """ Recall of model contacts
285 
286  The fraction of target contacts that are also present in model
287 
288  :type: :class:`int`
289  """
290  if self._n_trg_contacts != 0:
291  return self._n_intersection / self._n_trg_contacts
292  else:
293  return 0.0
294 
295  @property
296  def ics(self):
297  """ The Interface Contact Similarity score (ICS)
298 
299  Combination of :attr:`precision` and :attr:`recall` using the F1-measure
300 
301  :type: :class:`float`
302  """
303  p = self.precision
304  r = self.recall
305  nominator = p*r
306  denominator = p + r
307  if denominator != 0.0:
308  return 2*nominator/denominator
309  else:
310  return 0.0
311 
313  """
314  Holds data relevant to compute ips
315  """
316  def __init__(self, n_trg_int_res, n_mdl_int_res, n_union, n_intersection):
317  self._n_trg_int_res = n_trg_int_res
318  self._n_mdl_int_res = n_mdl_int_res
319  self._n_union = n_union
320  self._n_intersection = n_intersection
321 
322  @property
323  def n_trg_int_res(self):
324  """ Number of interface residues in target
325 
326  :type: :class:`int`
327  """
328  return self._n_trg_contacts
329 
330  @property
331  def n_mdl_int_res(self):
332  """ Number of interface residues in model
333 
334  :type: :class:`int`
335  """
336  return self._n_mdl_int_res
337 
338  @property
339  def precision(self):
340  """ Precision of model interface residues
341 
342  The fraction of model interface residues that are also interface
343  residues in target
344 
345  :type: :class:`int`
346  """
347  if self._n_mdl_int_res != 0:
348  return self._n_intersection / self._n_mdl_int_res
349  else:
350  return 0.0
351 
352  @property
353  def recall(self):
354  """ Recall of model interface residues
355 
356  The fraction of target interface residues that are also interface
357  residues in model
358 
359  :type: :class:`int`
360  """
361  if self._n_trg_int_res != 0:
362  return self._n_intersection / self._n_trg_int_res
363  else:
364  return 0.0
365 
366  @property
367  def ips(self):
368  """ The Interface Patch Similarity score (IPS)
369 
370  Jaccard coefficient of interface residues in model/target.
371  Technically thats :attr:`intersection`/:attr:`union`
372 
373  :type: :class:`float`
374  """
375  if(self._n_union > 0):
376  return self._n_intersection/self._n_union
377  return 0.0
378 
380  """ Helper object to compute Contact scores
381 
382  Tightly integrated into the mechanisms from the chain_mapping module.
383  The prefered way to derive an object of type :class:`ContactScorer` is
384  through the static constructor: :func:`~FromMappingResult`.
385 
386  Usage is the same as for :class:`ost.mol.alg.QSScorer`
387  """
388 
389  def __init__(self, target, chem_groups, model, alns,
390  contact_mode="aa", contact_d=5.0):
391  self._cent1 = ContactEntity(target, contact_mode = contact_mode,
392  contact_d = contact_d)
393  # ensure that target chain names match the ones in chem_groups
394  chem_group_ch_names = list(itertools.chain.from_iterable(chem_groups))
395  if self._cent1.chain_names != sorted(chem_group_ch_names):
396  raise RuntimeError(f"Expect exact same chain names in chem_groups "
397  f"and in target (which is processed to only "
398  f"contain peptides/nucleotides). target: "
399  f"{self._cent1.chain_names}, chem_groups: "
400  f"{chem_group_ch_names}")
401 
402  self._chem_groups = chem_groups
403  self._cent2 = ContactEntity(model, contact_mode = contact_mode,
404  contact_d = contact_d)
405  self._alns = alns
406 
407  # cache for mapped interface scores
408  # key: tuple of tuple ((qsent1_ch1, qsent1_ch2),
409  # ((qsent2_ch1, qsent2_ch2))
410  # value: tuple with four numbers required for computation of
411  # per-interface scores.
412  # The first two are relevant for ICS, the others for per
413  # interface IPS.
414  # 1: n_union_contacts
415  # 2: n_intersection_contacts
416  # 3: n_union_interface_residues
417  # 4: n_intersection_interface_residues
418  self._mapped_cache_interface = dict()
419 
420  # cache for mapped single chain scores
421  # for interface residues of single chains
422  # key: tuple: (qsent1_ch, qsent2_ch)
423  # value: tuple with two numbers required for computation of IPS
424  # 1: n_union
425  # 2: n_intersection
426  self._mapped_cache_sc = dict()
427 
428  @staticmethod
429  def FromMappingResult(mapping_result, contact_mode="aa", contact_d = 5.0):
430  """ The preferred way to get a :class:`ContactScorer`
431 
432  Static constructor that derives an object of type :class:`ContactScorer`
433  using a :class:`ost.mol.alg.chain_mapping.MappingResult`
434 
435  :param mapping_result: Data source
436  :type mapping_result: :class:`ost.mol.alg.chain_mapping.MappingResult`
437  """
438  contact_scorer = ContactScorer(mapping_result.target,
439  mapping_result.chem_groups,
440  mapping_result.model,
441  mapping_result.alns,
442  contact_mode = contact_mode,
443  contact_d = contact_d)
444  return contact_scorer
445 
446  @property
447  def cent1(self):
448  """ Represents *target*
449 
450  :type: :class:`ContactEntity`
451  """
452  return self._cent1
453 
454  @property
455  def chem_groups(self):
456  """ Groups of chemically equivalent chains in *target*
457 
458  Provided at object construction
459 
460  :type: :class:`list` of :class:`list` of :class:`str`
461  """
462  return self._chem_groups
463 
464  @property
465  def cent2(self):
466  """ Represents *model*
467 
468  :type: :class:`ContactEntity`
469  """
470  return self._cent2
471 
472  @property
473  def alns(self):
474  """ Alignments between chains in :attr:`~cent1` and :attr:`~cent2`
475 
476  Provided at object construction. Each alignment is accessible with
477  ``alns[(t_chain,m_chain)]``. First sequence is the sequence of the
478  respective chain in :attr:`~cent1`, second sequence the one from
479  :attr:`~cent2`.
480 
481  :type: :class:`dict` with key: :class:`tuple` of :class:`str`, value:
482  :class:`ost.seq.AlignmentHandle`
483  """
484  return self._alns
485 
486  def ScoreICS(self, mapping, check=True):
487  """ Computes ICS given chain mapping
488 
489  Again, the preferred way is to get *mapping* is from an object
490  of type :class:`ost.mol.alg.chain_mapping.MappingResult`.
491 
492  :param mapping: see
493  :attr:`ost.mol.alg.chain_mapping.MappingResult.mapping`
494  :type mapping: :class:`list` of :class:`list` of :class:`str`
495  :param check: Perform input checks, can be disabled for speed purposes
496  if you know what you're doing.
497  :type check: :class:`bool`
498  :returns: Result object of type :class:`ContactScorerResultICS`
499  """
500 
501  if check:
502  # ensure that dimensionality of mapping matches self.chem_groups
503  if len(self.chem_groups) != len(mapping):
504  raise RuntimeError("Dimensions of self.chem_groups and mapping "
505  "must match")
506  for a,b in zip(self.chem_groups, mapping):
507  if len(a) != len(b):
508  raise RuntimeError("Dimensions of self.chem_groups and "
509  "mapping must match")
510  # ensure that chain names in mapping are all present in cent2
511  for name in itertools.chain.from_iterable(mapping):
512  if name is not None and name not in self.cent2.chain_names:
513  raise RuntimeError(f"Each chain in mapping must be present "
514  f"in self.cent2. No match for "
515  f"\"{name}\"")
516 
517  flat_mapping = dict()
518  for a, b in zip(self.chem_groups, mapping):
519  flat_mapping.update({x: y for x, y in zip(a, b) if y is not None})
520 
521  return self.ICSFromFlatMapping(flat_mapping)
522 
523  def ScoreICSInterface(self, trg_ch1, trg_ch2, mdl_ch1, mdl_ch2):
524  """ Computes ICS scores only considering one interface
525 
526  This only works for interfaces that are computed in :func:`Score`, i.e.
527  interfaces for which the alignments are set up correctly.
528 
529  :param trg_ch1: Name of first interface chain in target
530  :type trg_ch1: :class:`str`
531  :param trg_ch2: Name of second interface chain in target
532  :type trg_ch2: :class:`str`
533  :param mdl_ch1: Name of first interface chain in model
534  :type mdl_ch1: :class:`str`
535  :param mdl_ch2: Name of second interface chain in model
536  :type mdl_ch2: :class:`str`
537  :returns: Result object of type :class:`ContactScorerResultICS`
538  :raises: :class:`RuntimeError` if no aln for trg_ch1/mdl_ch1 or
539  trg_ch2/mdl_ch2 is available.
540  """
541  if (trg_ch1, mdl_ch1) not in self.alns:
542  raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
543  f"mdl_ch1 ({mdl_ch1}) available. Did you "
544  f"construct the QSScorer object from a "
545  f"MappingResult and are trg_ch1 and mdl_ch1 "
546  f"mapped to each other?")
547  if (trg_ch2, mdl_ch2) not in self.alns:
548  raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
549  f"mdl_ch1 ({mdl_ch1}) available. Did you "
550  f"construct the QSScorer object from a "
551  f"MappingResult and are trg_ch1 and mdl_ch1 "
552  f"mapped to each other?")
553  trg_int = (trg_ch1, trg_ch2)
554  mdl_int = (mdl_ch1, mdl_ch2)
555  trg_int_r = (trg_ch2, trg_ch1)
556  mdl_int_r = (mdl_ch2, mdl_ch1)
557 
558  if trg_int in self.cent1.contacts:
559  n_trg = len(self.cent1.contacts[trg_int])
560  elif trg_int_r in self.cent1.contacts:
561  n_trg = len(self.cent1.contacts[trg_int_r])
562  else:
563  n_trg = 0
564 
565  if mdl_int in self.cent2.contacts:
566  n_mdl = len(self.cent2.contacts[mdl_int])
567  elif mdl_int_r in self.cent2.contacts:
568  n_mdl = len(self.cent2.contacts[mdl_int_r])
569  else:
570  n_mdl = 0
571 
572  n_union, n_intersection, _, _ = self._MappedInterfaceScores(trg_int, mdl_int)
573  return ContactScorerResultICS(n_trg, n_mdl, n_union, n_intersection)
574 
575  def ICSFromFlatMapping(self, flat_mapping):
576  """ Same as :func:`ScoreICS` but with flat mapping
577 
578  :param flat_mapping: Dictionary with target chain names as keys and
579  the mapped model chain names as value
580  :type flat_mapping: :class:`dict` with :class:`str` as key and value
581  :returns: Result object of type :class:`ContactScorerResultICS`
582  """
583  n_trg = sum([len(x) for x in self.cent1.contacts.values()])
584  n_mdl = sum([len(x) for x in self.cent2.contacts.values()])
585  n_union = 0
586  n_intersection = 0
587 
588  processed_cent2_interfaces = set()
589  for int1 in self.cent1.interacting_chains:
590  if int1[0] in flat_mapping and int1[1] in flat_mapping:
591  int2 = (flat_mapping[int1[0]], flat_mapping[int1[1]])
592  a, b, _, _ = self._MappedInterfaceScores(int1, int2)
593  n_union += a
594  n_intersection += b
595  processed_cent2_interfaces.add((min(int2), max(int2)))
596 
597  # process interfaces that only exist in qsent2
598  r_flat_mapping = {v:k for k,v in flat_mapping.items()} # reverse mapping
599  for int2 in self.cent2.interacting_chains:
600  if int2 not in processed_cent2_interfaces:
601  if int2[0] in r_flat_mapping and int2[1] in r_flat_mapping:
602  int1 = (r_flat_mapping[int2[0]], r_flat_mapping[int2[1]])
603  a, b, _, _ = self._MappedInterfaceScores(int1, int2)
604  n_union += a
605  n_intersection += b
606 
607  return ContactScorerResultICS(n_trg, n_mdl,
608  n_union, n_intersection)
609 
610  def ScoreIPS(self, mapping, check=True):
611  """ Computes IPS given chain mapping
612 
613  Again, the preferred way is to get *mapping* is from an object
614  of type :class:`ost.mol.alg.chain_mapping.MappingResult`.
615 
616  :param mapping: see
617  :attr:`ost.mol.alg.chain_mapping.MappingResult.mapping`
618  :type mapping: :class:`list` of :class:`list` of :class:`str`
619  :param check: Perform input checks, can be disabled for speed purposes
620  if you know what you're doing.
621  :type check: :class:`bool`
622  :returns: Result object of type :class:`ContactScorerResultIPS`
623  """
624 
625  if check:
626  # ensure that dimensionality of mapping matches self.chem_groups
627  if len(self.chem_groups) != len(mapping):
628  raise RuntimeError("Dimensions of self.chem_groups and mapping "
629  "must match")
630  for a,b in zip(self.chem_groups, mapping):
631  if len(a) != len(b):
632  raise RuntimeError("Dimensions of self.chem_groups and "
633  "mapping must match")
634  # ensure that chain names in mapping are all present in cent2
635  for name in itertools.chain.from_iterable(mapping):
636  if name is not None and name not in self.cent2.chain_names:
637  raise RuntimeError(f"Each chain in mapping must be present "
638  f"in self.cent2. No match for "
639  f"\"{name}\"")
640 
641  flat_mapping = dict()
642  for a, b in zip(self.chem_groups, mapping):
643  flat_mapping.update({x: y for x, y in zip(a, b) if y is not None})
644 
645  return self.IPSFromFlatMapping(flat_mapping)
646 
647  def ScoreIPSInterface(self, trg_ch1, trg_ch2, mdl_ch1, mdl_ch2):
648  """ Computes IPS scores only considering one interface
649 
650  This only works for interfaces that are computed in :func:`Score`, i.e.
651  interfaces for which the alignments are set up correctly.
652 
653  :param trg_ch1: Name of first interface chain in target
654  :type trg_ch1: :class:`str`
655  :param trg_ch2: Name of second interface chain in target
656  :type trg_ch2: :class:`str`
657  :param mdl_ch1: Name of first interface chain in model
658  :type mdl_ch1: :class:`str`
659  :param mdl_ch2: Name of second interface chain in model
660  :type mdl_ch2: :class:`str`
661  :returns: Result object of type :class:`ContactScorerResultIPS`
662  :raises: :class:`RuntimeError` if no aln for trg_ch1/mdl_ch1 or
663  trg_ch2/mdl_ch2 is available.
664  """
665  if (trg_ch1, mdl_ch1) not in self.alns:
666  raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
667  f"mdl_ch1 ({mdl_ch1}) available. Did you "
668  f"construct the QSScorer object from a "
669  f"MappingResult and are trg_ch1 and mdl_ch1 "
670  f"mapped to each other?")
671  if (trg_ch2, mdl_ch2) not in self.alns:
672  raise RuntimeError(f"No aln between trg_ch1 ({trg_ch1}) and "
673  f"mdl_ch1 ({mdl_ch1}) available. Did you "
674  f"construct the QSScorer object from a "
675  f"MappingResult and are trg_ch1 and mdl_ch1 "
676  f"mapped to each other?")
677  trg_int = (trg_ch1, trg_ch2)
678  mdl_int = (mdl_ch1, mdl_ch2)
679  trg_int_r = (trg_ch2, trg_ch1)
680  mdl_int_r = (mdl_ch2, mdl_ch1)
681 
682  if trg_int in self.cent1.contacts:
683  n_trg = len(self.cent1.contacts[trg_int])
684  elif trg_int_r in self.cent1.contacts:
685  n_trg = len(self.cent1.contacts[trg_int_r])
686  else:
687  n_trg = 0
688 
689  if mdl_int in self.cent2.contacts:
690  n_mdl = len(self.cent2.contacts[mdl_int])
691  elif mdl_int_r in self.cent2.contacts:
692  n_mdl = len(self.cent2.contacts[mdl_int_r])
693  else:
694  n_mdl = 0
695 
696  _, _, n_union, n_intersection = self._MappedInterfaceScores(trg_int, mdl_int)
697  return ContactScorerResultIPS(n_trg, n_mdl, n_union, n_intersection)
698 
699 
700  def IPSFromFlatMapping(self, flat_mapping):
701  """ Same as :func:`ScoreIPS` but with flat mapping
702 
703  :param flat_mapping: Dictionary with target chain names as keys and
704  the mapped model chain names as value
705  :type flat_mapping: :class:`dict` with :class:`str` as key and value
706  :returns: Result object of type :class:`ContactScorerResultIPS`
707  """
708  n_trg = sum([len(x) for x in self.cent1.interface_residues.values()])
709  n_mdl = sum([len(x) for x in self.cent2.interface_residues.values()])
710  n_union = 0
711  n_intersection = 0
712 
713  processed_cent2_chains = set()
714  for trg_ch in self.cent1.chain_names:
715  if trg_ch in flat_mapping:
716  a, b = self._MappedSCScores(trg_ch, flat_mapping[trg_ch])
717  n_union += a
718  n_intersection += b
719  processed_cent2_chains.add(flat_mapping[trg_ch])
720  else:
721  n_union += len(self.cent1.interface_residues[trg_ch])
722 
723  for mdl_ch in self._cent2.chain_names:
724  if mdl_ch not in processed_cent2_chains:
725  n_union += len(self.cent2.interface_residues[mdl_ch])
726 
727  return ContactScorerResultIPS(n_trg, n_mdl,
728  n_union, n_intersection)
729 
730 
731  def _MappedInterfaceScores(self, int1, int2):
732  key_one = (int1, int2)
733  if key_one in self._mapped_cache_interface:
734  return self._mapped_cache_interface[key_one]
735  key_two = ((int1[1], int1[0]), (int2[1], int2[0]))
736  if key_two in self._mapped_cache_interface:
737  return self._mapped_cache_interface[key_two]
738 
739  a, b, c, d = self._InterfaceScores(int1, int2)
740  self._mapped_cache_interface[key_one] = (a, b, c, d)
741  return (a, b, c, d)
742 
743  def _InterfaceScores(self, int1, int2):
744  if int1 in self.cent1.contacts:
745  ref_contacts = self.cent1.contacts[int1]
746  elif (int1[1], int1[0]) in self.cent1.contacts:
747  ref_contacts = self.cent1.contacts[(int1[1], int1[0])]
748  # need to reverse contacts
749  ref_contacts = set([(x[1], x[0]) for x in ref_contacts])
750  else:
751  ref_contacts = set() # no contacts at all
752 
753  if int2 in self.cent2.contacts:
754  mdl_contacts = self.cent2.contacts[int2]
755  elif (int2[1], int2[0]) in self.cent2.contacts:
756  mdl_contacts = self.cent2.contacts[(int2[1], int2[0])]
757  # need to reverse contacts
758  mdl_contacts = set([(x[1], x[0]) for x in mdl_contacts])
759  else:
760  mdl_contacts = set() # no contacts at all
761 
762  # indices in contacts lists are specific to the respective
763  # structures, need manual mapping from alignments
764  ch1_aln = self.alns[(int1[0], int2[0])]
765  ch2_aln = self.alns[(int1[1], int2[1])]
766  mapped_ref_contacts = set()
767  mapped_mdl_contacts = set()
768  for c in ref_contacts:
769  mapped_c = (ch1_aln.GetPos(0, c[0]), ch2_aln.GetPos(0, c[1]))
770  mapped_ref_contacts.add(mapped_c)
771  for c in mdl_contacts:
772  mapped_c = (ch1_aln.GetPos(1, c[0]), ch2_aln.GetPos(1, c[1]))
773  mapped_mdl_contacts.add(mapped_c)
774 
775  contact_union = len(mapped_ref_contacts.union(mapped_mdl_contacts))
776  contact_intersection = len(mapped_ref_contacts.intersection(mapped_mdl_contacts))
777 
778  # above, we computed the union and intersection on actual
779  # contacts. Here, we do the same on interface residues
780 
781  # process interface residues of chain one in interface
782  tmp_ref = set([x[0] for x in mapped_ref_contacts])
783  tmp_mdl = set([x[0] for x in mapped_mdl_contacts])
784  intres_union = len(tmp_ref.union(tmp_mdl))
785  intres_intersection = len(tmp_ref.intersection(tmp_mdl))
786 
787  # process interface residues of chain two in interface
788  tmp_ref = set([x[1] for x in mapped_ref_contacts])
789  tmp_mdl = set([x[1] for x in mapped_mdl_contacts])
790  intres_union += len(tmp_ref.union(tmp_mdl))
791  intres_intersection += len(tmp_ref.intersection(tmp_mdl))
792 
793  return (contact_union, contact_intersection,
794  intres_union, intres_intersection)
795 
796  def _MappedSCScores(self, ref_ch, mdl_ch):
797  if (ref_ch, mdl_ch) in self._mapped_cache_sc:
798  return self._mapped_cache_sc[(ref_ch, mdl_ch)]
799  n_union, n_intersection = self._SCScores(ref_ch, mdl_ch)
800  self._mapped_cache_sc[(ref_ch, mdl_ch)] = (n_union, n_intersection)
801  return (n_union, n_intersection)
802 
803  def _SCScores(self, ch1, ch2):
804  ref_int_res = self.cent1.interface_residues[ch1]
805  mdl_int_res = self.cent2.interface_residues[ch2]
806  aln = self.alns[(ch1, ch2)]
807  mapped_ref_int_res = set()
808  mapped_mdl_int_res = set()
809  for r_idx in ref_int_res:
810  mapped_ref_int_res.add(aln.GetPos(0, r_idx))
811  for r_idx in mdl_int_res:
812  mapped_mdl_int_res.add(aln.GetPos(1, r_idx))
813  return(len(mapped_ref_int_res.union(mapped_mdl_int_res)),
814  len(mapped_ref_int_res.intersection(mapped_mdl_int_res)))
815 
816 # specify public interface
817 __all__ = ('ContactEntity', 'ContactScorerResultICS', 'ContactScorerResultIPS', 'ContactScorer')