OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
scoring.py
Go to the documentation of this file.
1 import os
2 from ost import mol
3 from ost import seq
4 from ost import io
5 from ost import conop
6 from ost import settings
7 from ost import geom
8 from ost.mol.alg import lddt
9 from ost.mol.alg import qsscore
10 from ost.mol.alg import chain_mapping
11 from ost.mol.alg import stereochemistry
12 from ost.mol.alg import dockq
13 from ost.mol.alg.lddt import lDDTScorer
14 from ost.mol.alg.qsscore import QSScorer
15 from ost.mol.alg.contact_score import ContactScorer
16 from ost.mol.alg import Molck, MolckSettings
17 from ost import bindings
18 from ost.bindings import cadscore
19 from ost.bindings import tmtools
20 import numpy as np
21 
23  """Scorer specific for a reference/model pair
24 
25  Finds best possible binding site representation of reference in model given
26  lDDT score. Uses :class:`ost.mol.alg.chain_mapping.ChainMapper` to deal with
27  chain mapping.
28 
29  :param reference: Reference structure
30  :type reference: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
31  :param model: Model structure
32  :type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
33  :param residue_number_alignment: Passed to ChainMapper constructor
34  :type residue_number_alignment: :class:`bool`
35  """
36  def __init__(self, reference, model,
37  residue_number_alignment=False):
39  resnum_alignments=residue_number_alignment)
40  self.ref = self.chain_mapper.target
41  self.mdl = model
42 
43  def ScoreBS(self, ligand, radius = 4.0, lddt_radius=10.0):
44  """Computes binding site lDDT score given *ligand*. Best possible
45  binding site representation is selected by lDDT but other scores such as
46  CA based RMSD and GDT are computed too and returned.
47 
48  :param ligand: Defines the scored binding site, i.e. provides positions
49  to perform proximity search
50  :type ligand: r'((Residue)|(Chain)|(Entity))((View)|(Handle))'
51  :param radius: Reference residues with any atom position within *radius*
52  of *ligand* consitute the scored binding site
53  :type radius: :class:`float`
54  :param lddt_radius: Passed as *inclusion_radius* to
55  :class:`ost.mol.alg.lddt.lDDTScorer`
56  :type lddt_radius: :class:`float`
57  :returns: Object of type :class:`ost.mol.alg.chain_mapping.ReprResult`
58  containing all atom lDDT score and mapping information.
59  None if no representation could be found.
60  """
61 
62  # create view of reference binding site
63  ref_residues_hashes = set() # helper to keep track of added residues
64  for ligand_at in ligand.atoms:
65  close_atoms = self.ref.FindWithin(ligand_at.GetPos(), radius)
66  for close_at in close_atoms:
67  ref_res = close_at.GetResidue()
68  h = ref_res.handle.GetHashCode()
69  if h not in ref_residues_hashes:
70  ref_residues_hashes.add(h)
71 
72  # reason for doing that separately is to guarantee same ordering of
73  # residues as in underlying entity. (Reorder by ResNum seems only
74  # available on ChainHandles)
75  ref_bs = self.ref.CreateEmptyView()
76  for ch in self.ref.chains:
77  for r in ch.residues:
78  if r.handle.GetHashCode() in ref_residues_hashes:
79  ref_bs.AddResidue(r, mol.ViewAddFlag.INCLUDE_ALL)
80 
81  # gogogo
82  bs_repr = self.chain_mapper.GetRepr(ref_bs, self.mdl,
83  inclusion_radius = lddt_radius)
84  if len(bs_repr) >= 1:
85  return bs_repr[0]
86  else:
87  return None
88 
89 
90 class Scorer:
91  """ Helper class to access the various scores available from ost.mol.alg
92 
93  Deals with structure cleanup, chain mapping, interface identification etc.
94  Intermediate results are available as attributes.
95 
96  :param model: Model structure - a deep copy is available as :attr:`model`.
97  Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
98  is applied.
99  :type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
100  :param target: Target structure - a deep copy is available as :attr:`target`.
101  Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
102  is applied.
103  :type target: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
104  :param resnum_alignments: Whether alignments between chemically equivalent
105  chains in *model* and *target* can be computed
106  based on residue numbers. This can be assumed in
107  benchmarking setups such as CAMEO/CASP.
108  :type resnum_alignments: :class:`bool`
109  :param molck_settings: Settings used for Molck on *model* and *target*, if
110  set to None, a default object is constructed by
111  setting everything except rm_zero_occ_atoms and
112  colored to True in
113  :class:`ost.mol.alg.MolckSettings` constructor.
114  :type molck_settings: :class:`ost.mol.alg.MolckSettings`
115  :param cad_score_exec: Explicit path to voronota-cadscore executable from
116  voronota installation from
117  https://github.com/kliment-olechnovic/voronota. If
118  not given, voronota-cadscore must be in PATH if any
119  of the CAD score related attributes is requested.
120  :type cad_score_exec: :class:`str`
121  :param custom_mapping: Provide custom chain mapping between *model* and
122  *target*. Dictionary with target chain names as key
123  and model chain names as value.
124  :type custom_mapping: :class:`dict`
125  :param usalign_exec: Explicit path to USalign executable used to compute
126  TM-score. If not given, TM-score will be computed
127  with OpenStructure internal copy of USalign code.
128  :type usalign_exec: :class:`str`
129  :param lddt_no_stereochecks: Whether to compute lDDT without stereochemistry
130  checks
131  :type lddt_no_stereochecks: :class:`bool`
132  :param n_max_naive: Parameter for chain mapping. If the number of possible
133  mappings is <= *n_max_naive*, the full
134  mapping solution space is enumerated to find the
135  the optimum. A heuristic is used otherwise. The default
136  of 40320 corresponds to an octamer (8! = 40320).
137  A structure with stoichiometry A6B2 would be
138  6!*2! = 1440 etc.
139  :type n_max_naive: :class:`int`
140  :param oum: Override USalign Mapping. Inject mapping of :class:`Scorer`
141  object into USalign to compute TM-score. Experimental feature
142  with limitations.
143  :type oum: :class:`bool`
144  :param min_pep_length: Relevant parameter if short peptides are involved in
145  scoring. Minimum peptide length for a chain in the
146  target structure to be considered in chain mapping.
147  The chain mapping algorithm first performs an all vs.
148  all pairwise sequence alignment to identify \"equal\"
149  chains within the target structure. We go for simple
150  sequence identity there. Short sequences can be
151  problematic as they may produce high sequence identity
152  alignments by pure chance.
153  :type min_pep_length: :class:`int`
154  :param min_nuc_length: Relevant parameter if short nucleotides are involved
155  in scoring. Minimum nucleotide length for a chain in
156  the target structure to be considered in chain
157  mapping. The chain mapping algorithm first performs
158  an all vs. all pairwise sequence alignment to
159  identify \"equal\" chains within the target
160  structure. We go for simple sequence identity there.
161  Short sequences can be problematic as they may
162  produce high sequence identity alignments by pure
163  chance.
164  :type min_nuc_length: :class:`int`
165  """
166  def __init__(self, model, target, resnum_alignments=False,
167  molck_settings = None, cad_score_exec = None,
168  custom_mapping=None, usalign_exec = None,
169  lddt_no_stereochecks=False, n_max_naive=40320,
170  oum=False, min_pep_length = 6, min_nuc_length = 4):
171 
172  self._target_orig = target
173  self._model_orig = model
174 
175  if isinstance(self._model_orig, mol.EntityView):
176  self._model = mol.CreateEntityFromView(self._model_orig, False)
177  else:
178  self._model = self._model_orig.Copy()
179 
180  if isinstance(self._target_orig, mol.EntityView):
181  self._target = mol.CreateEntityFromView(self._target_orig, False)
182  else:
183  self._target = self._target_orig.Copy()
184 
185  if molck_settings is None:
186  molck_settings = MolckSettings(rm_unk_atoms=True,
187  rm_non_std=False,
188  rm_hyd_atoms=True,
189  rm_oxt_atoms=True,
190  rm_zero_occ_atoms=False,
191  colored=False,
192  map_nonstd_res=True,
193  assign_elem=True)
194  Molck(self._model, conop.GetDefaultLib(), molck_settings)
195  Molck(self._target, conop.GetDefaultLib(), molck_settings)
196  self._model = self._model.Select("peptide=True or nucleotide=True")
197  self._target = self._target.Select("peptide=True or nucleotide=True")
198 
199  # catch models which have empty chain names
200  for ch in self._model.chains:
201  if ch.GetName().strip() == "":
202  raise RuntimeError("Model chains must have valid chain names")
203  if ch.GetName().strip() == "'" or ch.GetName().strip() == '"':
204  raise RuntimeError("Model chains cannot be named with quote signs (' or \"\")")
205 
206  # catch targets which have empty chain names
207  for ch in self._target.chains:
208  if ch.GetName().strip() == "":
209  raise RuntimeError("Target chains must have valid chain names")
210  if ch.GetName().strip() == "'" or ch.GetName().strip() == '"':
211  raise RuntimeError("Target chains cannot be named with quote signs (' or \"\")")
212 
213  if resnum_alignments:
214  # In case of resnum_alignments, we have some requirements on
215  # residue numbers in the chain mapping: 1) no ins codes 2) strictly
216  # increasing residue numbers.
217  for ch in self._model.chains:
218  ins_codes = [r.GetNumber().GetInsCode() for r in ch.residues]
219  if len(set(ins_codes)) != 1 or ins_codes[0] != '\0':
220  raise RuntimeError("Residue numbers in each model chain "
221  "must not contain insertion codes if "
222  "resnum_alignments are enabled")
223  nums = [r.GetNumber().GetNum() for r in ch.residues]
224  if not all(i < j for i, j in zip(nums, nums[1:])):
225  raise RuntimeError("Residue numbers in each model chain "
226  "must be strictly increasing if "
227  "resnum_alignments are enabled")
228 
229  for ch in self._target.chains:
230  ins_codes = [r.GetNumber().GetInsCode() for r in ch.residues]
231  if len(set(ins_codes)) != 1 or ins_codes[0] != '\0':
232  raise RuntimeError("Residue numbers in each target chain "
233  "must not contain insertion codes if "
234  "resnum_alignments are enabled")
235  nums = [r.GetNumber().GetNum() for r in ch.residues]
236  if not all(i < j for i, j in zip(nums, nums[1:])):
237  raise RuntimeError("Residue numbers in each target chain "
238  "must be strictly increasing if "
239  "resnum_alignments are enabled")
240 
241  if usalign_exec is not None:
242  if not os.path.exists(usalign_exec):
243  raise RuntimeError(f"USalign exec ({usalign_exec}) "
244  f"not found")
245  if not os.access(usalign_exec, os.X_OK):
246  raise RuntimeError(f"USalign exec ({usalign_exec}) "
247  f"is not executable")
248 
249  self.resnum_alignments = resnum_alignments
250  self.cad_score_exec = cad_score_exec
251  self.usalign_exec = usalign_exec
252  self.lddt_no_stereochecks = lddt_no_stereochecks
253  self.n_max_naive = n_max_naive
254  self.oum = oum
255  self.min_pep_length = min_pep_length
256  self.min_nuc_length = min_nuc_length
257 
258  # lazily evaluated attributes
259  self._stereochecked_model = None
260  self._stereochecked_target = None
261  self._model_clashes = None
262  self._model_bad_bonds = None
263  self._model_bad_angles = None
264  self._target_clashes = None
265  self._target_bad_bonds = None
266  self._target_bad_angles = None
267  self._chain_mapper = None
268  self._mapping = None
269  self._model_interface_residues = None
270  self._target_interface_residues = None
271  self._aln = None
272  self._stereochecked_aln = None
273  self._pepnuc_aln = None
274 
275  # lazily constructed scorer objects
276  self._lddt_scorer = None
277  self._bb_lddt_scorer = None
278  self._qs_scorer = None
279  self._contact_scorer = None
280 
281  # lazily computed scores
282  self._lddt = None
283  self._local_lddt = None
284  self._bb_lddt = None
285  self._bb_local_lddt = None
286 
287  self._qs_global = None
288  self._qs_best = None
289  self._qs_target_interfaces = None
290  self._qs_model_interfaces = None
291  self._qs_interfaces = None
292  self._per_interface_qs_global = None
293  self._per_interface_qs_best = None
294 
295  self._contact_target_interfaces = None
296  self._contact_model_interfaces = None
297  self._native_contacts = None
298  self._model_contacts = None
299  self._ics_precision = None
300  self._ics_recall = None
301  self._ics = None
302  self._per_interface_ics_precision = None
303  self._per_interface_ics_recall = None
304  self._per_interface_ics = None
305  self._ips_precision = None
306  self._ips_recall = None
307  self._ips = None
308  self._per_interface_ics_precision = None
309  self._per_interface_ics_recall = None
310  self._per_interface_ics = None
311 
312  self._dockq_target_interfaces = None
313  self._dockq_interfaces = None
314  self._fnat = None
315  self._fnonnat = None
316  self._irmsd = None
317  self._lrmsd = None
318  self._nnat = None
319  self._nmdl = None
320  self._dockq_scores = None
321  self._dockq_ave = None
322  self._dockq_wave = None
323  self._dockq_ave_full = None
324  self._dockq_wave_full = None
325 
326  self._mapped_target_pos = None
327  self._mapped_model_pos = None
329  self._n_target_not_mapped = None
330  self._transform = None
331  self._gdtts = None
332  self._gdtha = None
333  self._rmsd = None
334 
335  self._cad_score = None
336  self._local_cad_score = None
337 
338  self._patch_qs = None
339  self._patch_dockq = None
340 
341  self._tm_score = None
342  self._usalign_mapping = None
343 
344  if custom_mapping is not None:
345  self._set_custom_mapping(custom_mapping)
346 
347  @property
348  def model(self):
349  """ Model with Molck cleanup
350 
351  :type: :class:`ost.mol.EntityHandle`
352  """
353  return self._model
354 
355  @property
356  def model_orig(self):
357  """ The original model passed at object construction
358 
359  :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
360  """
361  return self._model_orig
362 
363  @property
364  def target(self):
365  """ Target with Molck cleanup
366 
367  :type: :class:`ost.mol.EntityHandle`
368  """
369  return self._target
370 
371  @property
372  def target_orig(self):
373  """ The original target passed at object construction
374 
375  :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
376  """
377  return self._target_orig
378 
379  @property
380  def aln(self):
381  """ Alignments of :attr:`model`/:attr:`target` chains
382 
383  Alignments for each pair of chains mapped in :attr:`mapping`.
384  First sequence is target sequence, second sequence the model sequence.
385 
386  :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
387  """
388  if self._aln is None:
389  self._compute_aln()
390  return self._aln
391 
392  @property
393  def stereochecked_aln(self):
394  """ Stereochecked equivalent of :attr:`aln`
395 
396  The alignments may differ, as stereochecks potentially remove residues
397 
398  :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
399  """
400  if self._stereochecked_aln is None:
402  return self._stereochecked_aln
403 
404  @property
405  def pepnuc_aln(self):
406  """ Alignments of :attr:`model_orig`/:attr:`target_orig` chains
407 
408  Selects for peptide and nucleotide residues before sequence
409  extraction. Includes residues that would be removed by molck in
410  structure preprocessing.
411 
412  :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
413  """
414  if self._pepnuc_aln is None:
415  self._compute_pepnuc_aln()
416  return self._pepnuc_aln
417 
418  @property
420  """ View of :attr:`~model` that has stereochemistry checks applied
421 
422  First, a selection for peptide/nucleotide residues is performed,
423  secondly peptide sidechains with stereochemical irregularities are
424  removed (full residue if backbone atoms are involved). Irregularities
425  are clashes or bond lengths/angles more than 12 standard deviations
426  from expected values.
427 
428  :type: :class:`ost.mol.EntityView`
429  """
430  if self._stereochecked_model is None:
431  self._do_stereochecks()
432  return self._stereochecked_model
433 
434  @property
435  def model_clashes(self):
436  """ Clashing model atoms
437 
438  :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
439  """
440  if self._model_clashes is None:
441  self._do_stereochecks()
442  return self._model_clashes
443 
444  @property
445  def model_bad_bonds(self):
446  """ Model bonds with unexpected stereochemistry
447 
448  :type: :class:`list` of
449  :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
450  """
451  if self._model_bad_bonds is None:
452  self._do_stereochecks()
453  return self._model_bad_bonds
454 
455  @property
456  def model_bad_angles(self):
457  """ Model angles with unexpected stereochemistry
458 
459  :type: :class:`list` of
460  :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
461  """
462  if self._model_bad_angles is None:
463  self._do_stereochecks()
464  return self._model_bad_angles
465 
466  @property
468  """ Same as :attr:`~stereochecked_model` for :attr:`~target`
469 
470  :type: :class:`ost.mol.EntityView`
471  """
472  if self._stereochecked_target is None:
473  self._do_stereochecks()
474  return self._stereochecked_target
475 
476  @property
477  def target_clashes(self):
478  """ Clashing target atoms
479 
480  :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
481  """
482  if self._target_clashes is None:
483  self._do_stereochecks()
484  return self._target_clashes
485 
486  @property
487  def target_bad_bonds(self):
488  """ Target bonds with unexpected stereochemistry
489 
490  :type: :class:`list` of
491  :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
492  """
493  if self._target_bad_bonds is None:
494  self._do_stereochecks()
495  return self._target_bad_bonds
496 
497  @property
498  def target_bad_angles(self):
499  """ Target angles with unexpected stereochemistry
500 
501  :type: :class:`list` of
502  :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
503  """
504  if self._target_bad_angles is None:
505  self._do_stereochecks()
506  return self._target_bad_angles
507 
508  @property
509  def chain_mapper(self):
510  """ Chain mapper object for given :attr:`target`
511 
512  :type: :class:`ost.mol.alg.chain_mapping.ChainMapper`
513  """
514  if self._chain_mapper is None:
516  n_max_naive=1e9,
517  resnum_alignments=self.resnum_alignments,
518  min_pep_length=self.min_pep_length,
519  min_nuc_length=self.min_nuc_length)
520  return self._chain_mapper
521 
522  @property
523  def mapping(self):
524  """ Full chain mapping result for :attr:`target`/:attr:`model`
525 
526  :type: :class:`ost.mol.alg.chain_mapping.MappingResult`
527  """
528  if self._mapping is None:
529  self._mapping = \
530  self.chain_mapper.GetMapping(self.model,
531  n_max_naive = self.n_max_naive)
532  return self._mapping
533 
534  @property
536  """ Interface residues in :attr:`~model`
537 
538  Thats all residues having a contact with at least one residue from
539  another chain (CB-CB distance <= 8A, CA in case of Glycine)
540 
541  :type: :class:`dict` with chain names as key and and :class:`list`
542  with residue numbers of the respective interface residues.
543  """
544  if self._model_interface_residues is None:
546  self._get_interface_residues(self.model)
547  return self._model_interface_residues
548 
549  @property
551  """ Same as :attr:`~model_interface_residues` for :attr:`~target`
552 
553  :type: :class:`dict` with chain names as key and and :class:`list`
554  with residue numbers of the respective interface residues.
555  """
556  if self._target_interface_residues is None:
559  return self._target_interface_residues
560 
561  @property
562  def lddt_scorer(self):
563  """ lDDT scorer for :attr:`~stereochecked_target` (default parameters)
564 
565  :type: :class:`ost.mol.alg.lddt.lDDTScorer`
566  """
567  if self._lddt_scorer is None:
568  if self.lddt_no_stereochecks:
569  self._lddt_scorer = lDDTScorer(self.target)
570  else:
572  return self._lddt_scorer
573 
574  @property
575  def bb_lddt_scorer(self):
576  """ Backbone only lDDT scorer for :attr:`~target`
577 
578  No stereochecks applied for bb only lDDT which considers CA atoms
579  for peptides and C3' atoms for nucleotides.
580 
581  :type: :class:`ost.mol.alg.lddt.lDDTScorer`
582  """
583  if self._bb_lddt_scorer is None:
584  self._bb_lddt_scorer = lDDTScorer(self.target, bb_only=True)
585  return self._bb_lddt_scorer
586 
587  @property
588  def qs_scorer(self):
589  """ QS scorer constructed from :attr:`~mapping`
590 
591  The scorer object is constructed with default parameters and relates to
592  :attr:`~model` and :attr:`~target` (no stereochecks).
593 
594  :type: :class:`ost.mol.alg.qsscore.QSScorer`
595  """
596  if self._qs_scorer is None:
597  self._qs_scorer = QSScorer.FromMappingResult(self.mapping)
598  return self._qs_scorer
599 
600  @property
601  def contact_scorer(self):
602  if self._contact_scorer is None:
603  self._contact_scorer = ContactScorer.FromMappingResult(self.mapping)
604  return self._contact_scorer
605 
606 
607  @property
608  def lddt(self):
609  """ Global lDDT score in range [0.0, 1.0]
610 
611  Computed based on :attr:`~stereochecked_model`. In case of oligomers,
612  :attr:`~mapping` is used.
613 
614  :type: :class:`float`
615  """
616  if self._lddt is None:
617  self._compute_lddt()
618  return self._lddt
619 
620  @property
621  def local_lddt(self):
622  """ Per residue lDDT scores in range [0.0, 1.0]
623 
624  Computed based on :attr:`~stereochecked_model` but scores for all
625  residues in :attr:`~model` are reported. If a residue has been removed
626  by stereochemistry checks, the respective score is set to 0.0. If a
627  residue is not covered by the target or is in a chain skipped by the
628  chain mapping procedure (happens for super short chains), the respective
629  score is set to None. In case of oligomers, :attr:`~mapping` is used.
630 
631  :type: :class:`dict`
632  """
633  if self._local_lddt is None:
634  self._compute_lddt()
635  return self._local_lddt
636 
637  @property
638  def bb_lddt(self):
639  """ Backbone only global lDDT score in range [0.0, 1.0]
640 
641  Computed based on :attr:`~model` on backbone atoms only. This is CA for
642  peptides and C3' for nucleotides. No stereochecks are performed. In case
643  of oligomers, :attr:`~mapping` is used.
644 
645  :type: :class:`float`
646  """
647  if self._bb_lddt is None:
648  self._compute_bb_lddt()
649  return self._bb_lddt
650 
651  @property
652  def bb_local_lddt(self):
653  """ Backbone only per residue lDDT scores in range [0.0, 1.0]
654 
655  Computed based on :attr:`~model` on backbone atoms only. This is CA for
656  peptides and C3' for nucleotides. No stereochecks are performed. If a
657  residue is not covered by the target or is in a chain skipped by the
658  chain mapping procedure (happens for super short chains), the respective
659  score is set to None. In case of oligomers, :attr:`~mapping` is used.
660 
661  :type: :class:`dict`
662  """
663  if self._bb_local_lddt is None:
664  self._compute_bb_lddt()
665  return self._bb_local_lddt
666 
667  @property
668  def qs_global(self):
669  """ Global QS-score
670 
671  Computed based on :attr:`model` using :attr:`mapping`
672 
673  :type: :class:`float`
674  """
675  if self._qs_global is None:
676  self._compute_qs()
677  return self._qs_global
678 
679  @property
680  def qs_best(self):
681  """ Global QS-score - only computed on aligned residues
682 
683  Computed based on :attr:`model` using :attr:`mapping`. The QS-score
684  computation only considers contacts between residues with a mapping
685  between target and model. As a result, the score won't be lowered in
686  case of additional chains/residues in any of the structures.
687 
688  :type: :class:`float`
689  """
690  if self._qs_best is None:
691  self._compute_qs()
692  return self._qs_best
693 
694  @property
696  """ Interfaces in :attr:`~target` with non-zero contribution to
697  :attr:`~qs_global`/:attr:`~qs_best`
698 
699  Chain names are lexicographically sorted.
700 
701  :type: :class:`list` of :class:`tuple` with 2 elements each:
702  (trg_ch1, trg_ch2)
703  """
704  if self._qs_target_interfaces is None:
705  self._qs_target_interfaces = self.qs_scorer.qsent1.interacting_chains
706  self._qs_target_interfaces = \
707  [(min(x[0],x[1]), max(x[0],x[1])) for x in self._qs_target_interfaces]
708  return self._qs_target_interfaces
709 
710  @property
712  """ Interfaces in :attr:`~model` with non-zero contribution to
713  :attr:`~qs_global`/:attr:`~qs_best`
714 
715  Chain names are lexicographically sorted.
716 
717  :type: :class:`list` of :class:`tuple` with 2 elements each:
718  (mdl_ch1, mdl_ch2)
719  """
720  if self._qs_model_interfaces is None:
721  self._qs_model_interfaces = self.qs_scorer.qsent2.interacting_chains
722  self._qs_model_interfaces = \
723  [(min(x[0],x[1]), max(x[0],x[1])) for x in self._qs_model_interfaces]
724 
725  return self._qs_model_interfaces
726 
727  @property
728  def qs_interfaces(self):
729  """ Interfaces in :attr:`~qs_target_interfaces` that can be mapped
730  to :attr:`~model`.
731 
732  Target chain names are lexicographically sorted.
733 
734  :type: :class:`list` of :class:`tuple` with 4 elements each:
735  (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
736  """
737  if self._qs_interfaces is None:
738  self._qs_interfaces = list()
739  flat_mapping = self.mapping.GetFlatMapping()
740  for i in self.qs_target_interfaces:
741  if i[0] in flat_mapping and i[1] in flat_mapping:
742  self._qs_interfaces.append((i[0], i[1],
743  flat_mapping[i[0]],
744  flat_mapping[i[1]]))
745  return self._qs_interfaces
746 
747  @property
749  """ QS-score for each interface in :attr:`qs_interfaces`
750 
751  :type: :class:`list` of :class:`float`
752  """
753  if self._per_interface_qs_global is None:
755  return self._per_interface_qs_global
756 
757  @property
759  """ QS-score for each interface in :attr:`qs_interfaces`
760 
761  Only computed on aligned residues
762 
763  :type: :class:`list` of :class:`float`
764  """
765  if self._per_interface_qs_best is None:
767  return self._per_interface_qs_best
768 
769  @property
770  def native_contacts(self):
771  """ Native contacts
772 
773  A contact is a pair or residues from distinct chains that have
774  a minimal heavy atom distance < 5A. Contacts are specified as
775  :class:`tuple` with two strings in format:
776  <cname>.<rnum>.<ins_code>
777 
778  :type: :class:`list` of :class:`tuple`
779  """
780  if self._native_contacts is None:
781  self._native_contacts = self.contact_scorer.cent1.hr_contacts
782  return self._native_contacts
783 
784  @property
785  def model_contacts(self):
786  """ Same for model
787  """
788  if self._model_contacts is None:
789  self._model_contacts = self.contact_scorer.cent2.hr_contacts
790  return self._model_contacts
791 
792  @property
794  """ Interfaces in :class:`target` which have at least one contact
795 
796  Contact as defined in :attr:`~native_contacts`,
797  chain names are lexicographically sorted.
798 
799  :type: :class:`list` of :class:`tuple` with 2 elements each
800  (trg_ch1, trg_ch2)
801  """
802  if self._contact_target_interfaces is None:
803  tmp = self.contact_scorer.cent1.interacting_chains
804  tmp = [(min(x[0],x[1]), max(x[0],x[1])) for x in tmp]
805  self._contact_target_interfaces = tmp
806  return self._contact_target_interfaces
807 
808  @property
810  """ Interfaces in :class:`model` which have at least one contact
811 
812  Contact as defined in :attr:`~native_contacts`,
813  chain names are lexicographically sorted.
814 
815  :type: :class:`list` of :class:`tuple` with 2 elements each
816  (mdl_ch1, mdl_ch2)
817  """
818  if self._contact_model_interfaces is None:
819  tmp = self.contact_scorer.cent2.interacting_chains
820  tmp = [(min(x[0],x[1]), max(x[0],x[1])) for x in tmp]
821  self._contact_model_interfaces = tmp
822  return self._contact_model_interfaces
823 
824  @property
825  def ics_precision(self):
826  """ Fraction of model contacts that are also present in target
827 
828  :type: :class:`float`
829  """
830  if self._ics_precision is None:
831  self._compute_ics_scores()
832  return self._ics_precision
833 
834  @property
835  def ics_recall(self):
836  """ Fraction of target contacts that are correctly reproduced in model
837 
838  :type: :class:`float`
839  """
840  if self._ics_recall is None:
841  self._compute_ics_scores()
842  return self._ics_recall
843 
844  @property
845  def ics(self):
846  """ ICS (Interface Contact Similarity) score
847 
848  Combination of :attr:`~ics_precision` and :attr:`~ics_recall`
849  using the F1-measure
850 
851  :type: :class:`float`
852  """
853  if self._ics is None:
854  self._compute_ics_scores()
855  return self._ics
856 
857  @property
859  """ Per-interface ICS precision
860 
861  :attr:`~ics_precision` for each interface in
862  :attr:`~contact_target_interfaces`
863 
864  :type: :class:`list` of :class:`float`
865  """
866  if self._per_interface_ics_precision is None:
867  self._compute_ics_scores()
868  return self._per_interface_ics_precision
869 
870 
871  @property
873  """ Per-interface ICS recall
874 
875  :attr:`~ics_recall` for each interface in
876  :attr:`~contact_target_interfaces`
877 
878  :type: :class:`list` of :class:`float`
879  """
880  if self._per_interface_ics_recall is None:
881  self._compute_ics_scores()
882  return self._per_interface_ics_recall
883 
884  @property
885  def per_interface_ics(self):
886  """ Per-interface ICS (Interface Contact Similarity) score
887 
888  :attr:`~ics` for each interface in
889  :attr:`~contact_target_interfaces`
890 
891  :type: :class:`float`
892  """
893 
894  if self._per_interface_ics is None:
895  self._compute_ics_scores()
896  return self._per_interface_ics
897 
898 
899  @property
900  def ips_precision(self):
901  """ Fraction of model interface residues that are also interface
902  residues in target
903 
904  :type: :class:`float`
905  """
906  if self._ips_precision is None:
907  self._compute_ips_scores()
908  return self._ips_precision
909 
910  @property
911  def ips_recall(self):
912  """ Fraction of target interface residues that are also interface
913  residues in model
914 
915  :type: :class:`float`
916  """
917  if self._ips_recall is None:
918  self._compute_ips_scores()
919  return self._ips_recall
920 
921  @property
922  def ips(self):
923  """ IPS (Interface Patch Similarity) score
924 
925  Jaccard coefficient of interface residues in target and their mapped
926  counterparts in model
927 
928  :type: :class:`float`
929  """
930  if self._ips is None:
931  self._compute_ips_scores()
932  return self._ips
933 
934  @property
936  """ Per-interface IPS precision
937 
938  :attr:`~ips_precision` for each interface in
939  :attr:`~contact_target_interfaces`
940 
941  :type: :class:`list` of :class:`float`
942  """
943  if self._per_interface_ips_precision is None:
944  self._compute_ips_scores()
945  return self._per_interface_ips_precision
946 
947 
948  @property
950  """ Per-interface IPS recall
951 
952  :attr:`~ips_recall` for each interface in
953  :attr:`~contact_target_interfaces`
954 
955  :type: :class:`list` of :class:`float`
956  """
957  if self._per_interface_ics_recall is None:
958  self._compute_ips_scores()
959  return self._per_interface_ips_recall
960 
961  @property
962  def per_interface_ips(self):
963  """ Per-interface IPS (Interface Patch Similarity) score
964 
965  :attr:`~ips` for each interface in
966  :attr:`~contact_target_interfaces`
967 
968  :type: :class:`list` of :class:`float`
969  """
970 
971  if self._per_interface_ips is None:
972  self._compute_ips_scores()
973  return self._per_interface_ips
974 
975  @property
977  """ Interfaces in :attr:`target` that are relevant for DockQ
978 
979  In principle a subset of :attr:`~contact_target_interfaces` that only
980  contains peptide sequences. Chain names are lexicographically sorted.
981 
982  :type: :class:`list` of :class:`tuple` with 2 elements each:
983  (trg_ch1, trg_ch2)
984  """
985  if self._dockq_target_interfaces is None:
986  self._dockq_target_interfaces = list()
987  pep_seqs = set([s.GetName() for s in self.chain_mapper.polypep_seqs])
988  for interface in self.contact_target_interfaces:
989  if interface[0] in pep_seqs and interface[1] in pep_seqs:
990  self._dockq_target_interfaces.append(interface)
991  return self._dockq_target_interfaces
992 
993  @property
994  def dockq_interfaces(self):
995  """ Interfaces in :attr:`dockq_target_interfaces` that can be mapped
996  to model
997 
998  Target chain names are lexicographically sorted
999 
1000  :type: :class:`list` of :class:`tuple` with 4 elements each:
1001  (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
1002  """
1003  if self._dockq_interfaces is None:
1004  self._dockq_interfaces = list()
1005  flat_mapping = self.mapping.GetFlatMapping()
1006  for i in self.dockq_target_interfaces:
1007  if i[0] in flat_mapping and i[1] in flat_mapping:
1008  self._dockq_interfaces.append((i[0], i[1],
1009  flat_mapping[i[0]],
1010  flat_mapping[i[1]]))
1011  return self._dockq_interfaces
1012 
1013  @property
1014  def dockq_scores(self):
1015  """ DockQ scores for interfaces in :attr:`~dockq_interfaces`
1016 
1017  :class:`list` of :class:`float`
1018  """
1019  if self._dockq_scores is None:
1020  self._compute_dockq_scores()
1021  return self._dockq_scores
1022 
1023  @property
1024  def fnat(self):
1025  """ fnat scores for interfaces in :attr:`~dockq_interfaces`
1026 
1027  fnat: Fraction of native contacts that are also present in model
1028 
1029  :class:`list` of :class:`float`
1030  """
1031  if self._fnat is None:
1032  self._compute_dockq_scores()
1033  return self._fnat
1034 
1035  @property
1036  def nnat(self):
1037  """ N native contacts for interfaces in :attr:`~dockq_interfaces`
1038 
1039  :class:`list` of :class:`int`
1040  """
1041  if self._nnat is None:
1042  self._compute_dockq_scores()
1043  return self._nnat
1044 
1045  @property
1046  def nmdl(self):
1047  """ N model contacts for interfaces in :attr:`~dockq_interfaces`
1048 
1049  :class:`list` of :class:`int`
1050  """
1051  if self._nmdl is None:
1052  self._compute_dockq_scores()
1053  return self._nmdl
1054 
1055  @property
1056  def fnonnat(self):
1057  """ fnonnat scores for interfaces in :attr:`~dockq_interfaces`
1058 
1059  fnat: Fraction of model contacts that are not present in target
1060 
1061  :class:`list` of :class:`float`
1062  """
1063  if self._fnonnat is None:
1064  self._compute_dockq_scores()
1065  return self._fnonnat
1066 
1067  @property
1068  def irmsd(self):
1069  """ irmsd scores for interfaces in :attr:`~dockq_interfaces`
1070 
1071  irmsd: RMSD of interface (RMSD computed on N, CA, C, O atoms) which
1072  consists of each residue that has at least one heavy atom within 10A of
1073  other chain.
1074 
1075  :class:`list` of :class:`float`
1076  """
1077  if self._irmsd is None:
1078  self._compute_dockq_scores()
1079  return self._irmsd
1080 
1081  @property
1082  def lrmsd(self):
1083  """ lrmsd scores for interfaces in :attr:`~dockq_interfaces`
1084 
1085  lrmsd: The interfaces are superposed based on the receptor (rigid
1086  min RMSD superposition) and RMSD for the ligand is reported.
1087  Superposition and RMSD are based on N, CA, C and O positions,
1088  receptor is the chain contributing to the interface with more
1089  residues in total.
1090 
1091  :class:`list` of :class:`float`
1092  """
1093  if self._lrmsd is None:
1094  self._compute_dockq_scores()
1095  return self._lrmsd
1096 
1097  @property
1098  def dockq_ave(self):
1099  """ Average of DockQ scores in :attr:`dockq_scores`
1100 
1101  In its original implementation, DockQ only operates on single
1102  interfaces. Thus the requirement to combine scores for higher order
1103  oligomers.
1104 
1105  :type: :class:`float`
1106  """
1107  if self._dockq_ave is None:
1108  self._compute_dockq_scores()
1109  return self._dockq_ave
1110 
1111  @property
1112  def dockq_wave(self):
1113  """ Same as :attr:`dockq_ave`, weighted by native contacts
1114 
1115  :type: :class:`float`
1116  """
1117  if self._dockq_wave is None:
1118  self._compute_dockq_scores()
1119  return self._dockq_wave
1120 
1121  @property
1122  def dockq_ave_full(self):
1123  """ Same as :attr:`~dockq_ave` but penalizing for missing interfaces
1124 
1125  Interfaces that are not covered in model are added as 0.0
1126  in average computation.
1127 
1128  :type: :class:`float`
1129  """
1130  if self._dockq_ave_full is None:
1131  self._compute_dockq_scores()
1132  return self._dockq_ave_full
1133 
1134  @property
1135  def dockq_wave_full(self):
1136  """ Same as :attr:`~dockq_ave_full`, but weighted
1137 
1138  Interfaces that are not covered in model are added as 0.0 in
1139  average computations and the respective weights are derived from
1140  number of contacts in respective target interface.
1141  """
1142  if self._dockq_wave_full is None:
1143  self._compute_dockq_scores()
1144  return self._dockq_wave_full
1145 
1146  @property
1148  """ Mapped representative positions in target
1149 
1150  Thats CA positions for peptide residues and C3' positions for
1151  nucleotides. Has same length as :attr:`~mapped_model_pos` and mapping
1152  is based on :attr:`~mapping`.
1153 
1154  :type: :class:`ost.geom.Vec3List`
1155  """
1156  if self._mapped_target_pos is None:
1157  self._extract_mapped_pos()
1158  return self._mapped_target_pos
1159 
1160  @property
1161  def mapped_model_pos(self):
1162  """ Mapped representative positions in model
1163 
1164  Thats CA positions for peptide residues and C3' positions for
1165  nucleotides. Has same length as :attr:`~mapped_target_pos` and mapping
1166  is based on :attr:`~mapping`.
1167 
1168  :type: :class:`ost.geom.Vec3List`
1169  """
1170  if self._mapped_model_pos is None:
1171  self._extract_mapped_pos()
1172  return self._mapped_model_pos
1173 
1174  @property
1176  """ :attr:`~mapped_model_pos` with :attr:`~transform` applied
1177 
1178  :type: :class:`ost.geom.Vec3List`
1179  """
1180  if self._transformed_mapped_model_pos is None:
1183  self._transformed_mapped_model_pos.ApplyTransform(self.transform)
1184  return self._transformed_mapped_model_pos
1185 
1186  @property
1188  """ Number of target residues which have no mapping to model
1189 
1190  :type: :class:`int`
1191  """
1192  if self._n_target_not_mapped is None:
1193  self._extract_mapped_pos()
1194  return self._n_target_not_mapped
1195 
1196  @property
1197  def transform(self):
1198  """ Transform: :attr:`~mapped_model_pos` onto :attr:`~mapped_target_pos`
1199 
1200  Computed using Kabsch minimal rmsd algorithm
1201 
1202  :type: :class:`ost.geom.Mat4`
1203  """
1204  if self._transform is None:
1205  try:
1206  res = mol.alg.SuperposeSVD(self.mapped_model_pos,
1207  self.mapped_target_pos)
1208  self._transform = res.transformation
1209  except:
1210  self._transform = geom.Mat4()
1211  return self._transform
1212 
1213  @property
1214  def gdtts(self):
1215  """ GDT with thresholds: 8.0A, 4.0A, 2.0A and 1.0A
1216 
1217  Computed on :attr:`~transformed_mapped_model_pos` and
1218  :attr:`mapped_target_pos`
1219 
1220  :type: :class:`float`
1221  """
1222  if self._gdtts is None:
1223  n = \
1224  self.mapped_target_pos.GetGDTTS(self.transformed_mapped_model_pos,
1225  norm=False)
1226  n_full = 4*len(self.mapped_target_pos) + 4*self.n_target_not_mapped
1227  if n_full > 0:
1228  self._gdtts = float(n) / n_full
1229  else:
1230  self._gdtts = 0.0
1231  return self._gdtts
1232 
1233  @property
1234  def gdtha(self):
1235  """ GDT with thresholds: 4.0A, 2.0A, 1.0A and 0.5A
1236 
1237  Computed on :attr:`~transformed_mapped_model_pos` and
1238  :attr:`mapped_target_pos`
1239 
1240  :type: :class:`float`
1241  """
1242  if self._gdtha is None:
1243  n = \
1244  self.mapped_target_pos.GetGDTHA(self.transformed_mapped_model_pos,
1245  norm=False)
1246  n_full = 4*len(self.mapped_target_pos) + 4*self.n_target_not_mapped
1247  if n_full > 0:
1248  self._gdtha = float(n) / n_full
1249  else:
1250  self._gdtha = 0.0
1251  return self._gdtha
1252 
1253  @property
1254  def rmsd(self):
1255  """ RMSD
1256 
1257  Computed on :attr:`~transformed_mapped_model_pos` and
1258  :attr:`mapped_target_pos`
1259 
1260  :type: :class:`float`
1261  """
1262  if self._rmsd is None:
1263  self._rmsd = \
1264  self.mapped_target_pos.GetRMSD(self.transformed_mapped_model_pos)
1265  return self._rmsd
1266 
1267  @property
1268  def cad_score(self):
1269  """ The global CAD atom-atom (AA) score
1270 
1271  Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
1272  is used.
1273 
1274  :type: :class:`float`
1275  """
1276  if self._cad_score is None:
1277  self._compute_cad_score()
1278  return self._cad_score
1279 
1280  @property
1281  def local_cad_score(self):
1282  """ The per-residue CAD atom-atom (AA) scores
1283 
1284  Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
1285  is used.
1286 
1287  :type: :class:`dict`
1288  """
1289  if self._local_cad_score is None:
1290  self._compute_cad_score()
1291  return self._local_cad_score
1292 
1293  @property
1294  def patch_qs(self):
1295  """ Patch QS-scores for each residue in :attr:`model_interface_residues`
1296 
1297  Representative patches for each residue r in chain c are computed as
1298  follows:
1299 
1300  * mdl_patch_one: All residues in c with CB (CA for GLY) positions within
1301  8A of r and within 12A of residues from any other chain.
1302  * mdl_patch_two: Closest residue x to r in any other chain gets
1303  identified. Patch is then constructed by selecting all residues from
1304  any other chain within 8A of x and within 12A from any residue in c.
1305  * trg_patch_one: Chain name and residue number based mapping from
1306  mdl_patch_one
1307  * trg_patch_two: Chain name and residue number based mapping from
1308  mdl_patch_two
1309 
1310  Results are stored in the same manner as
1311  :attr:`model_interface_residues`, with corresponding scores instead of
1312  residue numbers. Scores for residues which are not
1313  :class:`mol.ChemType.AMINOACIDS` are set to None. Additionally,
1314  interface patches are derived from :attr:`model`. If they contain
1315  residues which are not covered by :attr:`target`, the score is set to
1316  None too.
1317 
1318  :type: :class:`dict` with chain names as key and and :class:`list`
1319  with scores of the respective interface residues.
1320  """
1321  if self._patch_qs is None:
1323  return self._patch_qs
1324 
1325  @property
1326  def patch_dockq(self):
1327  """ Same as :attr:`patch_qs` but for DockQ scores
1328  """
1329  if self._patch_dockq is None:
1331  return self._patch_dockq
1332 
1333  @property
1334  def tm_score(self):
1335  """ TM-score computed with USalign
1336 
1337  USalign executable can be specified with usalign_exec kwarg at Scorer
1338  construction, an OpenStructure internal copy of the USalign code is
1339  used otherwise.
1340 
1341  :type: :class:`float`
1342  """
1343  if self._tm_score is None:
1344  self._compute_tmscore()
1345  return self._tm_score
1346 
1347  @property
1348  def usalign_mapping(self):
1349  """ Mapping computed with USalign
1350 
1351  Dictionary with target chain names as key and model chain names as
1352  values. No guarantee that all chains are mapped. USalign executable
1353  can be specified with usalign_exec kwarg at Scorer construction, an
1354  OpenStructure internal copy of the USalign code is used otherwise.
1355 
1356  :type: :class:`dict`
1357  """
1358  if self._usalign_mapping is None:
1359  self._compute_tmscore()
1360  return self._usalign_mapping
1361 
1362  def _aln_helper(self, target, model):
1363  # perform required alignments - cannot take the alignments from the
1364  # mapping results as we potentially remove stuff there as compared
1365  # to self.model and self.target
1366  trg_seqs = dict()
1367  for ch in target.chains:
1368  cname = ch.GetName()
1369  s = ''.join([r.one_letter_code for r in ch.residues])
1370  s = seq.CreateSequence(ch.GetName(), s)
1371  s.AttachView(target.Select(f"cname='{cname}'"))
1372  trg_seqs[ch.GetName()] = s
1373  mdl_seqs = dict()
1374  for ch in model.chains:
1375  cname = ch.GetName()
1376  s = ''.join([r.one_letter_code for r in ch.residues])
1377  s = seq.CreateSequence(cname, s)
1378  s.AttachView(model.Select(f"cname='{cname}'"))
1379  mdl_seqs[ch.GetName()] = s
1380 
1381  alns = list()
1382  trg_pep_chains = [s.GetName() for s in self.chain_mapper.polypep_seqs]
1383  trg_nuc_chains = [s.GetName() for s in self.chain_mapper.polynuc_seqs]
1384  trg_pep_chains = set(trg_pep_chains)
1385  trg_nuc_chains = set(trg_nuc_chains)
1386  for trg_ch, mdl_ch in self.mapping.GetFlatMapping().items():
1387  if mdl_ch in mdl_seqs and trg_ch in trg_seqs:
1388  if trg_ch in trg_pep_chains:
1389  stype = mol.ChemType.AMINOACIDS
1390  elif trg_ch in trg_nuc_chains:
1391  stype = mol.ChemType.NUCLEOTIDES
1392  else:
1393  raise RuntimeError("Chain name inconsistency... ask "
1394  "Gabriel")
1395  alns.append(self.chain_mapper.Align(trg_seqs[trg_ch],
1396  mdl_seqs[mdl_ch],
1397  stype))
1398  alns[-1].AttachView(0, trg_seqs[trg_ch].GetAttachedView())
1399  alns[-1].AttachView(1, mdl_seqs[mdl_ch].GetAttachedView())
1400  return alns
1401 
1402  def _compute_pepnuc_aln(self):
1403  query = "peptide=true or nucleotide=true"
1404  pep_nuc_target = self.target_orig.Select(query)
1405  pep_nuc_model = self.model_orig.Select(query)
1406  self._pepnuc_aln = self._aln_helper(pep_nuc_target, pep_nuc_model)
1407 
1408  def _compute_aln(self):
1409  self._aln = self._aln_helper(self.target, self.model)
1410 
1411  def _compute_stereochecked_aln(self):
1413  self.stereochecked_model)
1414 
1415  def _compute_lddt(self):
1416  # lDDT requires a flat mapping with mdl_ch as key and trg_ch as value
1417  flat_mapping = self.mapping.GetFlatMapping(mdl_as_key=True)
1418 
1419  # make alignments accessible by mdl seq name
1420  stereochecked_alns = dict()
1421  for aln in self.stereochecked_aln:
1422  mdl_seq = aln.GetSequence(1)
1423  stereochecked_alns[mdl_seq.name] = aln
1424  alns = dict()
1425  for aln in self.aln:
1426  mdl_seq = aln.GetSequence(1)
1427  alns[mdl_seq.name] = aln
1428 
1429  # score variables to be set
1430  lddt_score = None
1431  local_lddt = None
1432 
1433  if self.lddt_no_stereochecks:
1434  lddt_chain_mapping = dict()
1435  for mdl_ch, trg_ch in flat_mapping.items():
1436  if mdl_ch in alns:
1437  lddt_chain_mapping[mdl_ch] = trg_ch
1438  lddt_score = self.lddt_scorer.lDDT(self.model,
1439  chain_mapping = lddt_chain_mapping,
1440  residue_mapping = alns,
1441  check_resnames=False,
1442  local_lddt_prop="lddt")[0]
1443  local_lddt = dict()
1444  for r in self.model.residues:
1445  cname = r.GetChain().GetName()
1446  if cname not in local_lddt:
1447  local_lddt[cname] = dict()
1448  if r.HasProp("lddt"):
1449  score = round(r.GetFloatProp("lddt"), 3)
1450  local_lddt[cname][r.GetNumber()] = score
1451  else:
1452  # not covered by trg or skipped in chain mapping procedure
1453  # the latter happens if its part of a super short chain
1454  local_lddt[cname][r.GetNumber()] = None
1455 
1456  else:
1457  lddt_chain_mapping = dict()
1458  for mdl_ch, trg_ch in flat_mapping.items():
1459  if mdl_ch in stereochecked_alns:
1460  lddt_chain_mapping[mdl_ch] = trg_ch
1461  lddt_score = self.lddt_scorer.lDDT(self.stereochecked_model,
1462  chain_mapping = lddt_chain_mapping,
1463  residue_mapping = stereochecked_alns,
1464  check_resnames=False,
1465  local_lddt_prop="lddt")[0]
1466  local_lddt = dict()
1467  for r in self.model.residues:
1468  cname = r.GetChain().GetName()
1469  if cname not in local_lddt:
1470  local_lddt[cname] = dict()
1471  if r.HasProp("lddt"):
1472  score = round(r.GetFloatProp("lddt"), 3)
1473  local_lddt[cname][r.GetNumber()] = score
1474  else:
1475  # rsc => residue stereo checked...
1476  mdl_res = self.stereochecked_model.FindResidue(cname, r.GetNumber())
1477  if mdl_res.IsValid():
1478  # not covered by trg or skipped in chain mapping procedure
1479  # the latter happens if its part of a super short chain
1480  local_lddt[cname][r.GetNumber()] = None
1481  else:
1482  # opt 1: removed by stereochecks => assign 0.0
1483  # opt 2: removed by stereochecks AND not covered by ref
1484  # => assign None
1485  # fetch trg residue from non-stereochecked aln
1486  trg_r = None
1487  if cname in flat_mapping:
1488  for col in alns[cname]:
1489  if col[0] != '-' and col[1] != '-':
1490  if col.GetResidue(1).number == r.number:
1491  trg_r = col.GetResidue(0)
1492  break
1493  if trg_r is None:
1494  local_lddt[cname][r.GetNumber()] = None
1495  else:
1496  local_lddt[cname][r.GetNumber()] = 0.0
1497 
1498  self._lddt = lddt_score
1499  self._local_lddt = local_lddt
1500 
1501  def _compute_bb_lddt(self):
1502 
1503  # make alignments accessible by mdl seq name
1504  alns = dict()
1505  for aln in self.aln:
1506  mdl_seq = aln.GetSequence(1)
1507  alns[mdl_seq.name] = aln
1508 
1509  # lDDT requires a flat mapping with mdl_ch as key and trg_ch as value
1510  flat_mapping = self.mapping.GetFlatMapping(mdl_as_key=True)
1511  lddt_chain_mapping = dict()
1512  for mdl_ch, trg_ch in flat_mapping.items():
1513  if mdl_ch in alns:
1514  lddt_chain_mapping[mdl_ch] = trg_ch
1515 
1516  lddt_score = self.bb_lddt_scorer.lDDT(self.model,
1517  chain_mapping = lddt_chain_mapping,
1518  residue_mapping = alns,
1519  check_resnames=False,
1520  local_lddt_prop="bb_lddt")[0]
1521  local_lddt = dict()
1522  for r in self.model.residues:
1523  cname = r.GetChain().GetName()
1524  if cname not in local_lddt:
1525  local_lddt[cname] = dict()
1526  if r.HasProp("bb_lddt"):
1527  score = round(r.GetFloatProp("bb_lddt"), 3)
1528  local_lddt[cname][r.GetNumber()] = score
1529  else:
1530  # not covered by trg or skipped in chain mapping procedure
1531  # the latter happens if its part of a super short chain
1532  local_lddt[cname][r.GetNumber()] = None
1533 
1534  self._bb_lddt = lddt_score
1535  self._bb_local_lddt = local_lddt
1536 
1537  def _compute_qs(self):
1538  qs_score_result = self.qs_scorer.Score(self.mapping.mapping)
1539  self._qs_global = qs_score_result.QS_global
1540  self._qs_best = qs_score_result.QS_best
1541 
1542  def _compute_per_interface_qs_scores(self):
1543  self._per_interface_qs_global = list()
1544  self._per_interface_qs_best = list()
1545 
1546  for interface in self.qs_interfaces:
1547  trg_ch1 = interface[0]
1548  trg_ch2 = interface[1]
1549  mdl_ch1 = interface[2]
1550  mdl_ch2 = interface[3]
1551  qs_res = self.qs_scorer.ScoreInterface(trg_ch1, trg_ch2,
1552  mdl_ch1, mdl_ch2)
1553  self._per_interface_qs_best.append(qs_res.QS_best)
1554  self._per_interface_qs_global.append(qs_res.QS_global)
1555 
1556  def _compute_ics_scores(self):
1557  contact_scorer_res = self.contact_scorer.ScoreICS(self.mapping.mapping)
1558  self._ics_precision = contact_scorer_res.precision
1559  self._ics_recall = contact_scorer_res.recall
1560  self._ics = contact_scorer_res.ics
1561  self._per_interface_ics_precision = list()
1562  self._per_interface_ics_recall = list()
1563  self._per_interface_ics = list()
1564  flat_mapping = self.mapping.GetFlatMapping()
1565  for trg_int in self.contact_target_interfaces:
1566  trg_ch1 = trg_int[0]
1567  trg_ch2 = trg_int[1]
1568  if trg_ch1 in flat_mapping and trg_ch2 in flat_mapping:
1569  mdl_ch1 = flat_mapping[trg_ch1]
1570  mdl_ch2 = flat_mapping[trg_ch2]
1571  res = self.contact_scorer.ScoreICSInterface(trg_ch1, trg_ch2,
1572  mdl_ch1, mdl_ch2)
1573  self._per_interface_ics_precision.append(res.precision)
1574  self._per_interface_ics_recall.append(res.recall)
1575  self._per_interface_ics.append(res.ics)
1576  else:
1577  self._per_interface_ics_precision.append(None)
1578  self._per_interface_ics_recall.append(None)
1579  self._per_interface_ics.append(None)
1580 
1581  def _compute_ips_scores(self):
1582  contact_scorer_res = self.contact_scorer.ScoreIPS(self.mapping.mapping)
1583  self._ips_precision = contact_scorer_res.precision
1584  self._ips_recall = contact_scorer_res.recall
1585  self._ips = contact_scorer_res.ips
1586 
1587  self._per_interface_ips_precision = list()
1588  self._per_interface_ips_recall = list()
1589  self._per_interface_ips = list()
1590  flat_mapping = self.mapping.GetFlatMapping()
1591  for trg_int in self.contact_target_interfaces:
1592  trg_ch1 = trg_int[0]
1593  trg_ch2 = trg_int[1]
1594  if trg_ch1 in flat_mapping and trg_ch2 in flat_mapping:
1595  mdl_ch1 = flat_mapping[trg_ch1]
1596  mdl_ch2 = flat_mapping[trg_ch2]
1597  res = self.contact_scorer.ScoreIPSInterface(trg_ch1, trg_ch2,
1598  mdl_ch1, mdl_ch2)
1599  self._per_interface_ips_precision.append(res.precision)
1600  self._per_interface_ips_recall.append(res.recall)
1601  self._per_interface_ips.append(res.ips)
1602  else:
1603  self._per_interface_ips_precision.append(None)
1604  self._per_interface_ips_recall.append(None)
1605  self._per_interface_ips.append(None)
1606 
1607  def _compute_dockq_scores(self):
1608  # lists with values in contact_target_interfaces
1609  self._dockq_scores = list()
1610  self._fnat = list()
1611  self._fnonnat = list()
1612  self._irmsd = list()
1613  self._lrmsd = list()
1614  self._nnat = list()
1615  self._nmdl = list()
1616 
1617  dockq_alns = dict()
1618  for aln in self.aln:
1619  trg_s = aln.GetSequence(0)
1620  mdl_s = aln.GetSequence(1)
1621  dockq_alns[(trg_s.GetName(), mdl_s.GetName())] = aln
1622 
1623  for interface in self.dockq_interfaces:
1624  trg_ch1 = interface[0]
1625  trg_ch2 = interface[1]
1626  mdl_ch1 = interface[2]
1627  mdl_ch2 = interface[3]
1628  aln1 = dockq_alns[(trg_ch1, mdl_ch1)]
1629  aln2 = dockq_alns[(trg_ch2, mdl_ch2)]
1630  res = dockq.DockQ(self.model, self.target, mdl_ch1, mdl_ch2,
1631  trg_ch1, trg_ch2, ch1_aln=aln1,
1632  ch2_aln=aln2)
1633  self._fnat.append(res["fnat"])
1634  self._fnonnat.append(res["fnonnat"])
1635  self._irmsd.append(res["irmsd"])
1636  self._lrmsd.append(res["lrmsd"])
1637  self._dockq_scores.append(res["DockQ"])
1638  self._nnat.append(res["nnat"])
1639  self._nmdl.append(res["nmdl"])
1640 
1641  # keep track of native counts in target interfaces which are
1642  # not covered in model in order to compute
1643  # dockq_ave_full/dockq_wave_full in the end
1644  not_covered_counts = list()
1645  proc_trg_interfaces = set([(x[0], x[1]) for x in self.dockq_interfaces])
1646  for interface in self.dockq_target_interfaces:
1647  if interface not in proc_trg_interfaces:
1648  # let's run DockQ with trg as trg/mdl in order to get the native
1649  # contacts out - no need to pass alns as the residue numbers
1650  # match for sure
1651  trg_ch1 = interface[0]
1652  trg_ch2 = interface[1]
1653  res = dockq.DockQ(self.target, self.target,
1654  trg_ch1, trg_ch2, trg_ch1, trg_ch2)
1655  not_covered_counts.append(res["nnat"])
1656 
1657  # there are 4 types of combined scores
1658  # - simple average
1659  # - average weighted by native_contacts
1660  # - the two above including nonmapped_contact_interfaces => set DockQ to 0.0
1661  scores = np.array(self._dockq_scores)
1662  weights = np.array(self._nnat)
1663  if len(scores) > 0:
1664  self._dockq_ave = np.mean(scores)
1665  else:
1666  self._dockq_ave = 0.0
1667  self._dockq_wave = np.sum(np.multiply(weights/np.sum(weights), scores))
1668  scores = np.append(scores, [0.0]*len(not_covered_counts))
1669  weights = np.append(weights, not_covered_counts)
1670  if len(scores) > 0:
1671  self._dockq_ave_full = np.mean(scores)
1672  else:
1673  self._dockq_ave_full = 0.0
1674  self._dockq_wave_full = np.sum(np.multiply(weights/np.sum(weights),
1675  scores))
1676 
1677  def _extract_mapped_pos(self):
1680  self._n_target_not_mapped = 0
1681  processed_trg_chains = set()
1682  for trg_ch, mdl_ch in self.mapping.GetFlatMapping().items():
1683  processed_trg_chains.add(trg_ch)
1684  aln = self.mapping.alns[(trg_ch, mdl_ch)]
1685  for col in aln:
1686  if col[0] != '-' and col[1] != '-':
1687  trg_res = col.GetResidue(0)
1688  mdl_res = col.GetResidue(1)
1689  trg_at = trg_res.FindAtom("CA")
1690  mdl_at = mdl_res.FindAtom("CA")
1691  if not trg_at.IsValid():
1692  trg_at = trg_res.FindAtom("C3'")
1693  if not mdl_at.IsValid():
1694  mdl_at = mdl_res.FindAtom("C3'")
1695  self._mapped_target_pos.append(trg_at.GetPos())
1696  self._mapped_model_pos.append(mdl_at.GetPos())
1697  elif col[0] != '-':
1698  self._n_target_not_mapped += 1
1699  # count number of trg residues from non-mapped chains
1700  for ch in self.mapping.target.chains:
1701  if ch.GetName() not in processed_trg_chains:
1702  self._n_target_not_mapped += len(ch.residues)
1703 
1704  def _compute_cad_score(self):
1705  if not self.resnum_alignments:
1706  raise RuntimeError("CAD score computations rely on residue numbers "
1707  "that are consistent between target and model "
1708  "chains, i.e. only work if resnum_alignments "
1709  "is True at Scorer construction.")
1710  try:
1711  cad_score_exec = \
1712  settings.Locate("voronota-cadscore",
1713  explicit_file_name=self.cad_score_exec)
1714  except Exception as e:
1715  raise RuntimeError("voronota-cadscore must be in PATH for CAD "
1716  "score scoring") from e
1717  cad_bin_dir = os.path.dirname(cad_score_exec)
1718  m = self.mapping.GetFlatMapping(mdl_as_key=True)
1719  cad_result = cadscore.CADScore(self.model, self.target,
1720  mode = "voronota",
1721  label="localcad",
1722  old_regime=False,
1723  cad_bin_path=cad_bin_dir,
1724  chain_mapping=m)
1725 
1726  local_cad = dict()
1727  for r in self.model.residues:
1728  cname = r.GetChain().GetName()
1729  if cname not in local_cad:
1730  local_cad[cname] = dict()
1731  if r.HasProp("localcad"):
1732  score = round(r.GetFloatProp("localcad"), 3)
1733  local_cad[cname][r.GetNumber()] = score
1734  else:
1735  local_cad[cname][r.GetNumber()] = None
1736 
1737  self._cad_score = cad_result.globalAA
1738  self._local_cad_score = local_cad
1739 
1740  def _get_repr_view(self, ent):
1741  """ Returns view with representative peptide atoms => CB, CA for GLY
1742 
1743  Ensures that each residue has exactly one atom with assertions
1744 
1745  :param ent: Entity for which you want the representative view
1746  :param ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
1747  :returns: :class:`ost.mol.EntityView` derived from ent
1748  """
1749  repr_ent = ent.Select("(aname=\"CB\" or (rname=\"GLY\" and aname=\"CA\"))")
1750  for r in repr_ent.residues:
1751  assert(len(r.atoms) == 1)
1752  return repr_ent
1753 
1754  def _get_interface_residues(self, ent):
1755  """ Get interface residues
1756 
1757  Thats all residues having a contact with at least one residue from another
1758  chain (CB-CB distance <= 8A, CA in case of Glycine)
1759 
1760  :param ent: Model for which to extract interface residues
1761  :type ent: :class:`ost.mol.EntityView`
1762  :returns: :class:`dict` with chain names as key and and :class:`list`
1763  with residue numbers of the respective interface residues.
1764  """
1765  # select for representative positions => CB, CA for GLY
1766  repr_ent = self._get_repr_view(ent)
1767  result = {ch.GetName(): list() for ch in ent.chains}
1768  for ch in ent.chains:
1769  cname = ch.GetName()
1770  sel = repr_ent.Select(f"(cname='{cname}' and 8 <> [cname!='{cname}'])")
1771  result[cname] = [r.GetNumber() for r in sel.residues]
1772  return result
1773 
1774  def _do_stereochecks(self):
1775  """ Perform stereochemistry checks on model and target
1776  """
1777  data = stereochemistry.GetDefaultStereoData()
1778  l_data = stereochemistry.GetDefaultStereoLinkData()
1779 
1780  a, b, c, d = stereochemistry.StereoCheck(self.model, stereo_data = data,
1781  stereo_link_data = l_data)
1782  self._stereochecked_model = a
1783  self._model_clashes = b
1784  self._model_bad_bonds = c
1785  self._model_bad_angles = d
1786 
1787  a, b, c, d = stereochemistry.StereoCheck(self.target, stereo_data = data,
1788  stereo_link_data = l_data)
1789  self._stereochecked_target = a
1790  self._target_clashes = b
1791  self._target_bad_bonds = c
1792  self._target_bad_angles = d
1793 
1794 
1795  def _get_interface_patches(self, mdl_ch, mdl_rnum):
1796  """ Select interface patches representative for specified residue
1797 
1798  The patches for specified residue r in chain c are selected as follows:
1799 
1800  * mdl_patch_one: All residues in c with CB (CA for GLY) positions within 8A
1801  of r and within 12A of residues from any other chain.
1802  * mdl_patch_two: Closest residue x to r in any other chain gets identified.
1803  Patch is then constructed by selecting all residues from any other chain
1804  within 8A of x and within 12A from any residue in c.
1805  * trg_patch_one: Chain name and residue number based mapping from
1806  mdl_patch_one
1807  * trg_patch_two: Chain name and residue number based mapping from
1808  mdl_patch_two
1809 
1810  :param mdl_ch: Name of chain in *self.model* of residue of interest
1811  :type mdl_ch: :class:`str`
1812  :param mdl_rnum: Residue number of residue of interest
1813  :type mdl_rnum: :class:`ost.mol.ResNum`
1814  :returns: Tuple with 5 elements: 1) :class:`bool` flag whether all residues
1815  in *mdl* patches are covered in *trg* 2) mtl_patch_one
1816  3) mdl_patch_two 4) trg_patch_one 5) trg_patch_two
1817  """
1818  # select for representative positions => CB, CA for GLY
1819  repr_mdl = self._get_repr_view(self.model.Select("peptide=true"))
1820 
1821  # get position for specified residue
1822  r = self.model.FindResidue(mdl_ch, mdl_rnum)
1823  if not r.IsValid():
1824  raise RuntimeError(f"Cannot find residue {mdl_rnum} in chain {mdl_ch}")
1825  if r.GetName() == "GLY":
1826  at = r.FindAtom("CA")
1827  else:
1828  at = r.FindAtom("CB")
1829  if not at.IsValid():
1830  raise RuntimeError("Cannot find interface views for res without CB/CA")
1831  r_pos = at.GetPos()
1832 
1833  # mdl_patch_one contains residues from the same chain as r
1834  # => all residues within 8A of r and within 12A of any other chain
1835 
1836  # q1 selects for everything in same chain and within 8A of r_pos
1837  q1 = f"(cname='{mdl_ch}' and 8 <> {{{r_pos[0]},{r_pos[1]},{r_pos[2]}}})"
1838  # q2 selects for everything within 12A of any other chain
1839  q2 = f"(12 <> [cname!={mdl_ch}])"
1840  mdl_patch_one = self.model.CreateEmptyView()
1841  sel = repr_mdl.Select(" and ".join([q1, q2]))
1842  for r in sel.residues:
1843  mdl_r = self.model.FindResidue(r.GetChain().GetName(), r.GetNumber())
1844  mdl_patch_one.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
1845 
1846  # mdl_patch_two contains residues from all other chains. In detail:
1847  # the closest residue to r is identified in any other chain, and the
1848  # patch is filled with residues that are within 8A of that residue and
1849  # within 12A of chain from r
1850  sel = repr_mdl.Select(f"(cname!='{mdl_ch}')")
1851  close_stuff = sel.FindWithin(r_pos, 8)
1852  min_pos = None
1853  min_dist = 42.0
1854  for close_at in close_stuff:
1855  dist = geom.Distance(r_pos, close_at.GetPos())
1856  if dist < min_dist:
1857  min_pos = close_at.GetPos()
1858  min_dist = dist
1859 
1860  # q1 selects for everything not in mdl_ch but within 8A of min_pos
1861  q1 = f"(cname!='{mdl_ch}' and 8 <> {{{min_pos[0]},{min_pos[1]},{min_pos[2]}}})"
1862  # q2 selects for everything within 12A of mdl_ch
1863  q2 = f"(12 <> [cname='{mdl_ch}'])"
1864  mdl_patch_two = self.model.CreateEmptyView()
1865  sel = repr_mdl.Select(" and ".join([q1, q2]))
1866  for r in sel.residues:
1867  mdl_r = self.model.FindResidue(r.GetChain().GetName(), r.GetNumber())
1868  mdl_patch_two.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
1869 
1870  # transfer mdl residues to trg
1871  flat_mapping = self.mapping.GetFlatMapping(mdl_as_key=True)
1872  full_trg_coverage = True
1873  trg_patch_one = self.target.CreateEmptyView()
1874  for r in mdl_patch_one.residues:
1875  trg_r = None
1876  mdl_cname = r.GetChain().GetName()
1877  if mdl_cname in flat_mapping:
1878  aln = self.mapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
1879  for col in aln:
1880  if col[0] != '-' and col[1] != '-':
1881  if col.GetResidue(1).GetNumber() == r.GetNumber():
1882  trg_r = col.GetResidue(0)
1883  break
1884  if trg_r is not None:
1885  trg_patch_one.AddResidue(trg_r.handle,
1886  mol.ViewAddFlag.INCLUDE_ALL)
1887  else:
1888  full_trg_coverage = False
1889 
1890  trg_patch_two = self.target.CreateEmptyView()
1891  for r in mdl_patch_two.residues:
1892  trg_r = None
1893  mdl_cname = r.GetChain().GetName()
1894  if mdl_cname in flat_mapping:
1895  aln = self.mapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
1896  for col in aln:
1897  if col[0] != '-' and col[1] != '-':
1898  if col.GetResidue(1).GetNumber() == r.GetNumber():
1899  trg_r = col.GetResidue(0)
1900  break
1901  if trg_r is not None:
1902  trg_patch_two.AddResidue(trg_r.handle,
1903  mol.ViewAddFlag.INCLUDE_ALL)
1904  else:
1905  full_trg_coverage = False
1906 
1907  return (full_trg_coverage, mdl_patch_one, mdl_patch_two,
1908  trg_patch_one, trg_patch_two)
1909 
1910  def _compute_patchqs_scores(self):
1911  self._patch_qs = dict()
1912  for cname, rnums in self.model_interface_residues.items():
1913  scores = list()
1914  for rnum in rnums:
1915  score = None
1916  r = self.model.FindResidue(cname, rnum)
1917  if r.IsValid() and r.GetChemType() == mol.ChemType.AMINOACIDS:
1918  full_trg_coverage, mdl_patch_one, mdl_patch_two, \
1919  trg_patch_one, trg_patch_two = \
1920  self._get_interface_patches(cname, rnum)
1921  if full_trg_coverage:
1922  score = self._patchqs(mdl_patch_one, mdl_patch_two,
1923  trg_patch_one, trg_patch_two)
1924  scores.append(score)
1925  self._patch_qs[cname] = scores
1926 
1927  def _compute_patchdockq_scores(self):
1928  self._patch_dockq = dict()
1929  for cname, rnums in self.model_interface_residues.items():
1930  scores = list()
1931  for rnum in rnums:
1932  score = None
1933  r = self.model.FindResidue(cname, rnum)
1934  if r.IsValid() and r.GetChemType() == mol.ChemType.AMINOACIDS:
1935  full_trg_coverage, mdl_patch_one, mdl_patch_two, \
1936  trg_patch_one, trg_patch_two = \
1937  self._get_interface_patches(cname, rnum)
1938  if full_trg_coverage:
1939  score = self._patchdockq(mdl_patch_one, mdl_patch_two,
1940  trg_patch_one, trg_patch_two)
1941  scores.append(score)
1942  self._patch_dockq[cname] = scores
1943 
1944  def _patchqs(self, mdl_patch_one, mdl_patch_two, trg_patch_one, trg_patch_two):
1945  """ Score interface residue patches with QS-score
1946 
1947  In detail: Construct two entities with two chains each. First chain
1948  consists of residues from <x>_patch_one and second chain consists of
1949  <x>_patch_two. The returned score is the QS-score between the two
1950  entities
1951 
1952  :param mdl_patch_one: Interface patch representing scored residue
1953  :type mdl_patch_one: :class:`ost.mol.EntityView`
1954  :param mdl_patch_two: Interface patch representing scored residue
1955  :type mdl_patch_two: :class:`ost.mol.EntityView`
1956  :param trg_patch_one: Interface patch representing scored residue
1957  :type trg_patch_one: :class:`ost.mol.EntityView`
1958  :param trg_patch_two: Interface patch representing scored residue
1959  :type trg_patch_two: :class:`ost.mol.EntityView`
1960  :returns: PatchQS score
1961  """
1962  qs_ent_mdl = self._qs_ent_from_patches(mdl_patch_one, mdl_patch_two)
1963  qs_ent_trg = self._qs_ent_from_patches(trg_patch_one, trg_patch_two)
1964 
1965  alnA = seq.CreateAlignment()
1966  s = ''.join([r.one_letter_code for r in mdl_patch_one.residues])
1967  alnA.AddSequence(seq.CreateSequence("A", s))
1968  s = ''.join([r.one_letter_code for r in trg_patch_one.residues])
1969  alnA.AddSequence(seq.CreateSequence("A", s))
1970 
1971  alnB = seq.CreateAlignment()
1972  s = ''.join([r.one_letter_code for r in mdl_patch_two.residues])
1973  alnB.AddSequence(seq.CreateSequence("B", s))
1974  s = ''.join([r.one_letter_code for r in trg_patch_two.residues])
1975  alnB.AddSequence(seq.CreateSequence("B", s))
1976  alns = {("A", "A"): alnA, ("B", "B"): alnB}
1977 
1978  scorer = QSScorer(qs_ent_mdl, [["A"], ["B"]], qs_ent_trg, alns)
1979  score_result = scorer.Score([["A"], ["B"]])
1980 
1981  return score_result.QS_global
1982 
1983  def _patchdockq(self, mdl_patch_one, mdl_patch_two, trg_patch_one,
1984  trg_patch_two):
1985  """ Score interface residue patches with DockQ
1986 
1987  In detail: Construct two entities with two chains each. First chain
1988  consists of residues from <x>_patch_one and second chain consists of
1989  <x>_patch_two. The returned score is the QS-score between the two
1990  entities
1991 
1992  :param mdl_patch_one: Interface patch representing scored residue
1993  :type mdl_patch_one: :class:`ost.mol.EntityView`
1994  :param mdl_patch_two: Interface patch representing scored residue
1995  :type mdl_patch_two: :class:`ost.mol.EntityView`
1996  :param trg_patch_one: Interface patch representing scored residue
1997  :type trg_patch_one: :class:`ost.mol.EntityView`
1998  :param trg_patch_two: Interface patch representing scored residue
1999  :type trg_patch_two: :class:`ost.mol.EntityView`
2000  :returns: DockQ score
2001  """
2002  m = self._qs_ent_from_patches(mdl_patch_one, mdl_patch_two)
2003  t = self._qs_ent_from_patches(trg_patch_one, trg_patch_two)
2004  dockq_result = dockq.DockQ(t, m, "A", "B", "A", "B")
2005  if dockq_result["nnat"] > 0:
2006  return dockq_result["DockQ"]
2007  return 0.0
2008 
2009  def _qs_ent_from_patches(self, patch_one, patch_two):
2010  """ Constructs Entity with two chains named "A" and "B""
2011 
2012  Blindly adds all residues from *patch_one* to chain A and residues from
2013  patch_two to chain B.
2014  """
2015  ent = mol.CreateEntity()
2016  ed = ent.EditXCS()
2017  added_ch = ed.InsertChain("A")
2018  for r in patch_one.residues:
2019  added_r = ed.AppendResidue(added_ch, r.GetName())
2020  added_r.SetChemClass(str(r.GetChemClass()))
2021  for a in r.atoms:
2022  ed.InsertAtom(added_r, a.handle)
2023  added_ch = ed.InsertChain("B")
2024  for r in patch_two.residues:
2025  added_r = ed.AppendResidue(added_ch, r.GetName())
2026  added_r.SetChemClass(str(r.GetChemClass()))
2027  for a in r.atoms:
2028  ed.InsertAtom(added_r, a.handle)
2029  return ent
2030 
2031  def _set_custom_mapping(self, mapping):
2032  """ sets self._mapping with a full blown MappingResult object
2033 
2034  :param mapping: mapping with trg chains as key and mdl ch as values
2035  :type mapping: :class:`dict`
2036  """
2037 
2038  chain_mapper = self.chain_mapper
2039  chem_mapping, chem_group_alns, mdl = \
2040  chain_mapper.GetChemMapping(self.model)
2041 
2042  # now that we have a chem mapping, lets do consistency checks
2043  # - check whether chain names are unique and available in structures
2044  # - check whether the mapped chains actually map to the same chem groups
2045  if len(mapping) != len(set(mapping.keys())):
2046  raise RuntimeError(f"Expect unique trg chain names in mapping. Got "
2047  f"{mapping.keys()}")
2048  if len(mapping) != len(set(mapping.values())):
2049  raise RuntimeError(f"Expect unique mdl chain names in mapping. Got "
2050  f"{mapping.values()}")
2051 
2052  trg_chains = set([ch.GetName() for ch in chain_mapper.target.chains])
2053  mdl_chains = set([ch.GetName() for ch in mdl.chains])
2054  for k,v in mapping.items():
2055  if k not in trg_chains:
2056  raise RuntimeError(f"Target chain \"{k}\" is not available "
2057  f"in target processed for chain mapping "
2058  f"({trg_chains})")
2059  if v not in mdl_chains:
2060  raise RuntimeError(f"Model chain \"{v}\" is not available "
2061  f"in model processed for chain mapping "
2062  f"({mdl_chains})")
2063 
2064  for trg_ch, mdl_ch in mapping.items():
2065  trg_group_idx = None
2066  mdl_group_idx = None
2067  for idx, group in enumerate(chain_mapper.chem_groups):
2068  if trg_ch in group:
2069  trg_group_idx = idx
2070  break
2071  for idx, group in enumerate(chem_mapping):
2072  if mdl_ch in group:
2073  mdl_group_idx = idx
2074  break
2075  if trg_group_idx is None or mdl_group_idx is None:
2076  raise RuntimeError("Could not establish a valid chem grouping "
2077  "of chain names provided in custom mapping.")
2078 
2079  if trg_group_idx != mdl_group_idx:
2080  raise RuntimeError(f"Chem group mismatch in custom mapping: "
2081  f"target chain \"{trg_ch}\" groups with the "
2082  f"following chemically equivalent target "
2083  f"chains: "
2084  f"{chain_mapper.chem_groups[trg_group_idx]} "
2085  f"but model chain \"{mdl_ch}\" maps to the "
2086  f"following target chains: "
2087  f"{chain_mapper.chem_groups[mdl_group_idx]}")
2088 
2089  pairs = set([(trg_ch, mdl_ch) for trg_ch, mdl_ch in mapping.items()])
2090  ref_mdl_alns = \
2091  chain_mapping._GetRefMdlAlns(chain_mapper.chem_groups,
2092  chain_mapper.chem_group_alignments,
2093  chem_mapping,
2094  chem_group_alns,
2095  pairs = pairs)
2096 
2097  # translate mapping format
2098  final_mapping = list()
2099  for ref_chains in chain_mapper.chem_groups:
2100  mapped_mdl_chains = list()
2101  for ref_ch in ref_chains:
2102  if ref_ch in mapping:
2103  mapped_mdl_chains.append(mapping[ref_ch])
2104  else:
2105  mapped_mdl_chains.append(None)
2106  final_mapping.append(mapped_mdl_chains)
2107 
2108  alns = dict()
2109  for ref_group, mdl_group in zip(chain_mapper.chem_groups,
2110  final_mapping):
2111  for ref_ch, mdl_ch in zip(ref_group, mdl_group):
2112  if ref_ch is not None and mdl_ch is not None:
2113  aln = ref_mdl_alns[(ref_ch, mdl_ch)]
2114  trg_view = chain_mapper.target.Select(f"cname='{ref_ch}'")
2115  mdl_view = mdl.Select(f"cname='{mdl_ch}'")
2116  aln.AttachView(0, trg_view)
2117  aln.AttachView(1, mdl_view)
2118  alns[(ref_ch, mdl_ch)] = aln
2119 
2120  self._mapping = chain_mapping.MappingResult(chain_mapper.target, mdl,
2121  chain_mapper.chem_groups,
2122  chem_mapping,
2123  final_mapping, alns)
2124 
2125  def _compute_tmscore(self):
2126  res = None
2127  if self.usalign_exec is None:
2128  if self.oum:
2129  flat_mapping = self.mapping.GetFlatMapping()
2130  res = res = bindings.WrappedMMAlign(self.model, self.target,
2131  mapping=flat_mapping)
2132  else:
2133  res = bindings.WrappedMMAlign(self.model, self.target)
2134  else:
2135  if self.oum:
2136  flat_mapping = self.mapping.GetFlatMapping()
2137  res = tmtools.USAlign(self.model, self.target,
2138  usalign = self.usalign_exec,
2139  custom_chain_mapping = flat_mapping)
2140  else:
2141  res = tmtools.USAlign(self.model, self.target,
2142  usalign = self.usalign_exec)
2143 
2144  self._tm_score = res.tm_score
2145  self._usalign_mapping = dict()
2146  for a,b in zip(res.ent1_mapped_chains, res.ent2_mapped_chains):
2147  self._usalign_mapping[b] = a
void Molck(ost::mol::EntityHandle &ent, ost::conop::CompoundLibPtr lib, const MolckSettings &settings, bool prune=true)
Real Distance(const Vec3 &p1, const Vec3 &p2)
returns the distance between two points
Definition: vecmat3_op.hh:199
definition of EntityView
Definition: entity_view.hh:86