OpenStructure
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
ligand_scoring_base.py
Go to the documentation of this file.
1 from contextlib import contextmanager
2 import numpy as np
3 import networkx
4 
5 import ost
6 from ost import io
7 from ost import mol
8 from ost import conop
9 from ost import LogWarning, LogScript, LogInfo, LogVerbose, LogDebug, GetVerbosityLevel, PushVerbosityLevel, PopVerbosityLevel
10 from ost.mol.alg import chain_mapping
11 
12 
13 @contextmanager
14 def _SinkVerbosityLevel(n=1):
15  """ Context manager to temporarily reduce the verbosity level by n.
16 
17  Example usage:
18  with _SinkVerbosityLevel(2):
19  LogVerbose("Test")
20  Will only write "Test" in script level (2 above)
21  """
22  PushVerbosityLevel(GetVerbosityLevel() - n)
23  try:
24  yield
25  except:
26  raise
27  finally:
29 
30 
31 def _QualifiedAtomNotation(a):
32  """Return a parsable string of the atom in the format:
33  ChainName.ResidueNumber.InsertionCode.AtomName."""
34  resnum = a.residue.number
35  return "{cname}.{rnum}.{ins_code}.{aname}".format(
36  cname=a.chain.name,
37  rnum=resnum.num,
38  ins_code=resnum.ins_code.strip("\u0000"),
39  aname=a.name,
40  )
41 
42 
43 def _QualifiedResidueNotation(r):
44  """Return a parsable string of the residue in the format:
45  ChainName.ResidueNumber.InsertionCode."""
46  resnum = r.number
47  return "{cname}.{rnum}.{ins_code}".format(
48  cname=r.chain.name,
49  rnum=resnum.num,
50  ins_code=resnum.ins_code.strip("\u0000"),
51  )
52 
53 
54 def CleanHydrogens(ent, clib):
55  """ Ligand scoring helper - Returns copy of *ent* without hydrogens
56 
57  Non-standard hydrogen naming can cause trouble in residue property
58  assignment which is done by the :class:`ost.conop.RuleBasedProcessor` when
59  loading. In fact, residue property assignment is not done for every residue
60  that has unknown atoms according to the chemical component dictionary. This
61  function therefore re-processes the entity after removing hydrogens.
62 
63  :param ent: Entity to clean
64  :type ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
65  :param clib: Compound library to perform re-processing after hydrogen
66  removal.
67  :type clib: :class:`ost.conop.CompoundLib`
68  :returns: Cleaned and re-processed ent
69  """
70  cleaned_ent = mol.CreateEntityFromView(ent.Select(
71  "ele != H and ele != D"), include_exlusive_atoms=False)
72  # process again to set missing residue properties due to non standard
73  # hydrogens
74  processor = conop.RuleBasedProcessor(clib)
75  processor.Process(cleaned_ent)
76  return cleaned_ent
77 
78 
79 def MMCIFPrep(mmcif_path, biounit=None, extract_nonpoly=False,
80  fault_tolerant=False):
81  """ Ligand scoring helper - Prepares :class:`LigandScorer` input from mmCIF
82 
83  Only performs gentle cleanup of hydrogen atoms. Further cleanup is delegated
84  to :class:`LigandScorer`.
85 
86  :param mmcif_path: Path to mmCIF file that contains polymer and optionally
87  non-polymer entities
88  :type mmcif_path: :class:`str`
89  :param biounit: If given, construct specified biounit from mmCIF AU
90  :type biounit: :class:`str`
91  :param extract_nonpoly: Additionally returns a list of
92  :class:`ost.mol.EntityHandle`
93  objects representing all non-polymer (ligand)
94  entities.
95  :type extract_nonpoly: :class:`bool`
96  :param fault_tolerant: Passed as parameter to :func:`ost.io.LoadMMCIF`
97  :type fault_tolerant: :class:`bool`
98  :returns: :class:`ost.mol.EntityHandle` which only contains polymer
99  entities representing the receptor structure. If *extract_nonpoly*
100  is True, a tuple is returned which additionally contains a
101  :class:`list` of :class:`ost.mol.EntityHandle`, where each
102  :class:`ost.mol.EntityHandle` represents a non-polymer (ligand).
103  """
104  clib = conop.GetDefaultLib()
105  if not clib:
106  ost.LogError("A compound library is required. "
107  "Please refer to the OpenStructure website: "
108  "https://openstructure.org/docs/conop/compoundlib/.")
109  raise RuntimeError("No compound library found")
110 
111  mmcif_entity, mmcif_info = io.LoadMMCIF(mmcif_path, info=True,
112  fault_tolerant=fault_tolerant)
113  mmcif_entity = CleanHydrogens(mmcif_entity, clib)
114 
115  # get AU chain names representing polymer entities
116  polymer_entity_ids = mmcif_info.GetEntityIdsOfType("polymer")
117  polymer_chain_names = list()
118  for ch in mmcif_entity.chains:
119  if mmcif_info.GetMMCifEntityIdTr(ch.name) in polymer_entity_ids:
120  polymer_chain_names.append(ch.name)
121 
122  # get AU chain names representing non-polymer entities
123  non_polymer_entity_ids = mmcif_info.GetEntityIdsOfType("non-polymer")
124  non_polymer_chain_names = list()
125  for ch in mmcif_entity.chains:
126  if mmcif_info.GetMMCifEntityIdTr(ch.name) in non_polymer_entity_ids:
127  non_polymer_chain_names.append(ch.name)
128 
129  # construct biounit if necessary
130  if biounit is not None:
131  biounit_found = False
132  for bu in mmcif_info.biounits:
133  if bu.id == biounit:
134  mmcif_entity = mol.alg.CreateBU(mmcif_entity, bu)
135  biounit_found = True
136  break
137  if not biounit_found:
138  raise RuntimeError(f"Specified biounit '{biounit}' not in "
139  f"{mmcif_path}")
140 
141  # assign generic properties for selection later on
142  non_poly_id = 0
143  for ch in mmcif_entity.chains:
144  cname = None
145  if biounit is not None:
146  # if a biounit is constructed, you get chain names like: 1.YOLO
147  # we cannot simply split by '.' since '.' is an allowed character
148  # in chain names. => split by first occurence
149  dot_index = ch.name.find('.')
150  if dot_index == -1:
151  cname = ch.name
152  else:
153  cname = ch.name[dot_index+1:]
154  else:
155  cname = ch.name
156 
157  if cname in polymer_chain_names:
158  ch.SetIntProp("poly", 1)
159  if cname in non_polymer_chain_names:
160  ch.SetIntProp("nonpolyid", non_poly_id)
161  non_poly_id += 1
162 
163  poly_sel = mmcif_entity.Select("gcpoly:0=1")
164  poly_ent = mol.CreateEntityFromView(poly_sel, True)
165 
166  if extract_nonpoly == False:
167  return poly_ent
168 
169  non_poly_sel = mmcif_entity.Select("gcnonpoly:0=1")
170  non_poly_entities = list()
171  for i in range(non_poly_id):
172  view = mmcif_entity.Select(f"gcnonpolyid:{non_poly_id}={i}")
173  if view.GetResidueCount() != 1:
174  raise RuntimeError(f"Expect non-polymer entities in "
175  f"{mmcif_path} to contain exactly 1 "
176  f"residue. Got {ch.GetResidueCount()} "
177  f"in chain {ch.name}")
178  compound = clib.FindCompound(view.residues[0].name)
179  if compound is None:
180  raise RuntimeError(f"Can only extract non-polymer entities if "
181  f"respective residues are available in PDB "
182  f"component dictionary. Can't find "
183  f"\"{view.residues[0].name}\"")
184 
185  non_poly_entities.append(mol.CreateEntityFromView(view, True))
186 
187  return (poly_ent, non_poly_entities)
188 
189 
190 def PDBPrep(pdb_path, fault_tolerant=False):
191  """ Ligand scoring helper - Prepares :class:`LigandScorer` input from PDB
192 
193  Only performs gentle cleanup of hydrogen atoms. Further cleanup is delegated
194  to :class:`LigandScorer`. There is no logic to extract ligands from PDB
195  files. Ligands must be provided separately as SDF files in these cases.
196 
197  :param pdb_path: Path to PDB file that contains polymer entities
198  :type pdb_path: :class:`str`
199  :param fault_tolerant: Passed as parameter to :func:`ost.io.LoadPDB`
200  :type fault_tolerant: :class:`bool`
201  :returns: :class:`EntityHandle` from loaded file.
202  """
203  clib = conop.GetDefaultLib()
204  if not clib:
205  ost.LogError("A compound library is required. "
206  "Please refer to the OpenStructure website: "
207  "https://openstructure.org/docs/conop/compoundlib/.")
208  raise RuntimeError("No compound library found")
209 
210  pdb_entity = io.LoadPDB(pdb_path, fault_tolerant=fault_tolerant)
211  pdb_entity = CleanHydrogens(pdb_entity, clib)
212 
213  return pdb_entity
214 
215 
217  """ Scorer to compute various small molecule ligand (non polymer) scores.
218 
219  :class:`LigandScorer` is an abstract base class dealing with all the setup,
220  data storage, enumerating ligand symmetries and target/model ligand
221  matching/assignment. But actual score computation is delegated to child
222  classes.
223 
224  At the moment, two such classes are available:
225 
226  * :class:`ost.mol.alg.ligand_scoring_lddtpli.LDDTPLIScorer`
227  that assesses the conservation of protein-ligand
228  contacts (LDDT-PLI);
229  * :class:`ost.mol.alg.ligand_scoring_scrmsd.SCRMSDScorer`
230  that computes a binding-site superposed, symmetry-corrected RMSD
231  (BiSyRMSD) and ligand pocket LDDT (LDDT-LP).
232 
233  All versus all scores are available through the lazily computed
234  :attr:`score_matrix`. However, many things can go wrong... be it even
235  something as simple as two ligands not matching. Error states therefore
236  encode scoring issues. An Issue for a particular ligand is indicated by a
237  non-zero state in :attr:`model_ligand_states`/:attr:`target_ligand_states`.
238  This invalidates pairwise scores of such a ligand with all other ligands.
239  This and other issues in pairwise score computation are reported in
240  :attr:`state_matrix` which has the same size as :attr:`score_matrix`.
241  Only if the respective location is 0, a valid pairwise score can be
242  expected. The states and their meaning can be explored with code::
243 
244  for state_code, (short_desc, desc) in scorer_obj.state_decoding.items():
245  print(state_code)
246  print(short_desc)
247  print(desc)
248 
249  A common use case is to derive a one-to-one mapping between ligands in
250  the model and the target for which :class:`LigandScorer` provides an
251  automated :attr:`assignment` procedure.
252  By default, only exact matches between target and model ligands are
253  considered. This is a problem when the target only contains a subset
254  of the expected atoms (for instance if atoms are missing in an
255  experimental structure, which often happens in the PDB). With
256  `substructure_match=True`, complete model ligands can be scored against
257  partial target ligands. One problem with this approach is that it is
258  very easy to find good matches to small, irrelevant ligands like EDO, CO2
259  or GOL. The assignment algorithm therefore considers the coverage,
260  expressed as the fraction of atoms of the model ligand atoms covered in the
261  target. Higher coverage matches are prioritized, but a match with a better
262  score will be preferred if it falls within a window of `coverage_delta`
263  (by default 0.2) of a worse-scoring match. As a result, for instance,
264  with a delta of 0.2, a low-score match with coverage 0.96 would be
265  preferred over a high-score match with coverage 0.70.
266 
267  Assumptions:
268 
269  Unlike most of OpenStructure, this class does not assume that the ligands
270  (either for the model or the target) are part of the PDB component
271  dictionary. They may have arbitrary residue names. Residue names do not
272  have to match between the model and the target. Matching is based on
273  the calculation of isomorphisms which depend on the atom element name and
274  atom connectivity (bond order is ignored).
275  It is up to the caller to ensure that the connectivity of atoms is properly
276  set before passing any ligands to this class. Ligands with improper
277  connectivity will lead to bogus results.
278 
279  This only applies to the ligand. The rest of the model and target
280  structures (protein, nucleic acids) must still follow the usual rules and
281  contain only residues from the compound library. Structures are cleaned up
282  according to constructor documentation. We advise to
283  use the :func:`MMCIFPrep` and :func:`PDBPrep` for loading which already
284  clean hydrogens and, in the case of MMCIF, optionally extract ligands ready
285  to be used by the :class:`LigandScorer` based on "non-polymer" entity types.
286  In case of PDB file format, ligands must be loaded separately as SDF files.
287 
288  Only polymers (protein and nucleic acids) of model and target are considered
289  for ligand binding sites. The
290  :class:`ost.mol.alg.chain_mapping.ChainMapper` is used to enumerate possible
291  mappings of these chains. In short: identical chains in the target are
292  grouped based on pairwise sequence identity
293  (see pep_seqid_thr/nuc_seqid_thr param). Each model chain is assigned to
294  one of these groups (see mdl_map_pep_seqid_thr/mdl_map_nuc_seqid_thr param).
295  To avoid spurious matches, only polymers of a certain length are considered
296  in this matching procedure (see min_pep_length/min_nuc_length param).
297  Shorter polymers are never mapped and do not contribute to scoring.
298 
299  Here is an example of how to setup a scorer::
300 
301  from ost.mol.alg.ligand_scoring_scrmsd import SCRMSDScorer
302 
303  # Load data
304  # Structure model in PDB format, containing the receptor only
305  model = PDBPrep("path_to_model.pdb")
306  # Ligand model as SDF file
307  model_ligand = io.LoadEntity("path_to_ligand.sdf", format="sdf")
308  # Target loaded from mmCIF, containing the ligand
309  target, target_ligands = MMCIFPrep("path_to_target.cif",
310  extract_nonpoly=True)
311 
312  # Setup scorer object and compute SCRMSD
313  model_ligands = [model_ligand.Select("ele != H")]
314  sc = SCRMSDScorer(model, target, model_ligands, target_ligands)
315 
316  # Perform assignment and read respective scores
317  for lig_pair in sc.assignment:
318  trg_lig = sc.target_ligands[lig_pair[0]]
319  mdl_lig = sc.model_ligands[lig_pair[1]]
320  score = sc.score_matrix[lig_pair[0], lig_pair[1]]
321  print(f"Score for {trg_lig} and {mdl_lig}: {score}")
322 
323  # check cleanup in model and target structure:
324  print("model cleanup:", sc.model_cleanup_log)
325  print("target cleanup:", sc.target_cleanup_log)
326 
327 
328  :param model: Model structure - a deep copy is available as :attr:`model`.
329  The model undergoes the following cleanup steps which are
330  dependent on :class:`ost.conop.CompoundLib` returned by
331  :func:`ost.conop.GetDefaultLib`: 1) removal
332  of hydrogens, 2) removal of residues for which there is no
333  entry in :class:`ost.conop.CompoundLib`, 3) removal of
334  residues that are not peptide linking or nucleotide linking
335  according to :class:`ost.conop.CompoundLib` 4) removal of
336  atoms that are not defined for respective residues in
337  :class:`ost.conop.CompoundLib`. Except step 1), every cleanup
338  is logged with :class:`ost.LogLevel` Warning and a report is
339  available as :attr:`model_cleanup_log`.
340  :type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
341  :param target: Target structure - same processing as *model*.
342  :type target: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
343  :param model_ligands: Model ligands, as a list of
344  :class:`ost.mol.ResidueHandle`/
345  :class:`ost.mol.ResidueView`/
346  :class:`ost.mol.EntityHandle`/
347  :class:`ost.mol.EntityView`. For
348  :class:`ost.mol.EntityHandle`/
349  :class:`ost.mol.EntityView`, each residue is
350  considered to be an individual ligand.
351  All ligands are copied into a separate
352  :class:`ost.mol.EntityHandle` available as
353  :attr:`model_ligand_ent` and the respective
354  list of ligands is available as :attr:`model_ligands`.
355  :type model_ligands: :class:`list`
356  :param target_ligands: Target ligands, same processing as model ligands.
357  :type target_ligands: :class:`list`
358  :param resnum_alignments: Whether alignments between chemically equivalent
359  chains in *model* and *target* can be computed
360  based on residue numbers. This can be assumed in
361  benchmarking setups such as CAMEO/CASP.
362  :type resnum_alignments: :class:`bool`
363  :param substructure_match: Set this to True to allow incomplete (i.e.
364  partially resolved) target ligands.
365  :type substructure_match: :class:`bool`
366  :param coverage_delta: the coverage delta for partial ligand assignment.
367  :type coverage_delta: :class:`float`
368  :param max_symmetries: If more than that many isomorphisms exist for
369  a target-ligand pair, it will be ignored and reported
370  as unassigned.
371  :type max_symmetries: :class:`int`
372  :param min_pep_length: Relevant parameter if short peptides are involved in
373  the polymer binding site. Minimum peptide length for
374  a chain to be considered in chain mapping.
375  The chain mapping algorithm first performs an all vs.
376  all pairwise sequence alignment to identify \"equal\"
377  chains within the target structure. We go for simple
378  sequence identity there. Short sequences can be
379  problematic as they may produce high sequence identity
380  alignments by pure chance.
381  :type min_pep_length: :class:`int`
382  :param min_nuc_length: Same for nucleotides
383  :type min_nuc_length: :class:`int`
384  :param pep_seqid_thr: Parameter that affects identification of identical
385  chains in target - see
386  :class:`ost.mol.alg.chain_mapping.ChainMapper`
387  :type pep_seqid_thr: :class:`float`
388  :param nuc_seqid_thr: Parameter that affects identification of identical
389  chains in target - see
390  :class:`ost.mol.alg.chain_mapping.ChainMapper`
391  :type nuc_seqid_thr: :class:`float`
392  :param mdl_map_pep_seqid_thr: Parameter that affects mapping of model chains
393  to target chains - see
394  :class:`ost.mol.alg.chain_mapping.ChainMapper`
395  :type mdl_map_pep_seqid_thr: :class:`float`
396  :param mdl_map_nuc_seqid_thr: Parameter that affects mapping of model chains
397  to target chains - see
398  :class:`ost.mol.alg.chain_mapping.ChainMapper`
399  :type mdl_map_nuc_seqid_thr: :class:`float`
400  """
401 
402  def __init__(self, model, target, model_ligands, target_ligands,
403  resnum_alignments=False, substructure_match=False,
404  coverage_delta=0.2, max_symmetries=1e5,
405  rename_ligand_chain=False, min_pep_length = 6,
406  min_nuc_length = 4, pep_seqid_thr = 95.,
407  nuc_seqid_thr = 95.,
408  mdl_map_pep_seqid_thr = 0.,
409  mdl_map_nuc_seqid_thr = 0.):
410 
411  if isinstance(model, mol.EntityView):
412  self._model_model = mol.CreateEntityFromView(model, False)
413  elif isinstance(model, mol.EntityHandle):
414  self._model_model = model.Copy()
415  else:
416  raise RuntimeError("model must be of type EntityView/EntityHandle")
417 
418  if isinstance(target, mol.EntityView):
419  self._target_target = mol.CreateEntityFromView(target, False)
420  elif isinstance(target, mol.EntityHandle):
421  self._target_target = target.Copy()
422  else:
423  raise RuntimeError("target must be of type EntityView/EntityHandle")
424 
425  clib = conop.GetDefaultLib()
426  if not clib:
427  ost.LogError("A compound library is required. "
428  "Please refer to the OpenStructure website: "
429  "https://openstructure.org/docs/conop/compoundlib/.")
430  raise RuntimeError("No compound library found")
431  self._target_target, self._target_cleanup_log_target_cleanup_log = \
432  self._cleanup_polymer_ent_cleanup_polymer_ent(self._target_target, clib)
433  self._model_model, self._model_cleanup_log_model_cleanup_log = \
434  self._cleanup_polymer_ent_cleanup_polymer_ent(self._model_model, clib)
435 
436  # keep ligands separate from polymer entities
437  self._target_ligand_ent_target_ligand_ent = mol.CreateEntity()
438  self._model_ligand_ent_model_ligand_ent = mol.CreateEntity()
439  target_ligand_ent_ed = self._target_ligand_ent_target_ligand_ent.EditXCS(mol.BUFFERED_EDIT)
440  model_ligand_ent_ed = self._model_ligand_ent_model_ligand_ent.EditXCS(mol.BUFFERED_EDIT)
441 
442  self._target_ligands_target_ligands = list()
443  for l in target_ligands:
444  if isinstance(l, mol.EntityView) or isinstance(l, mol.EntityHandle):
445  for r in l.residues:
446  self._target_ligands_target_ligands.append(self._copy_ligand_copy_ligand(r, self._target_ligand_ent_target_ligand_ent,
447  target_ligand_ent_ed,
448  rename_ligand_chain))
449  elif isinstance(l, mol.ResidueHandle) or isinstance(l, mol.ResidueView):
450  self._target_ligands_target_ligands.append(self._copy_ligand_copy_ligand(l, self._target_ligand_ent_target_ligand_ent,
451  target_ligand_ent_ed,
452  rename_ligand_chain))
453  else:
454  raise RuntimeError("ligands must be of type EntityView/"
455  "EntityHandle/ResidueView/ResidueHandle")
456 
457  self._model_ligands_model_ligands = list()
458  for l in model_ligands:
459  if isinstance(l, mol.EntityView) or isinstance(l, mol.EntityHandle):
460  for r in l.residues:
461  self._model_ligands_model_ligands.append(self._copy_ligand_copy_ligand(r, self._model_ligand_ent_model_ligand_ent,
462  model_ligand_ent_ed,
463  rename_ligand_chain))
464  elif isinstance(l, mol.ResidueHandle) or isinstance(l, mol.ResidueView):
465  self._model_ligands_model_ligands.append(self._copy_ligand_copy_ligand(l, self._model_ligand_ent_model_ligand_ent,
466  model_ligand_ent_ed,
467  rename_ligand_chain))
468  else:
469  raise RuntimeError("ligands must be of type EntityView/"
470  "EntityHandle/ResidueView/ResidueHandle")
471 
472 
473  if len(self.model_ligandsmodel_ligands) == 0:
474  LogWarning("No ligands in the model")
475  if len(self.target_ligandstarget_ligands) == 0:
476  raise ValueError("No ligand in the model and in the target")
477 
478  self._resnum_alignments_resnum_alignments = resnum_alignments
479  self._substructure_match_substructure_match = substructure_match
480  self._coverage_delta_coverage_delta = coverage_delta
481  self._max_symmetries_max_symmetries = max_symmetries
482  self._min_pep_length_min_pep_length = min_pep_length
483  self._min_nuc_length_min_nuc_length = min_nuc_length
484  self._pep_seqid_thr_pep_seqid_thr = pep_seqid_thr
485  self._nuc_seqid_thr_nuc_seqid_thr = nuc_seqid_thr
486  self._mdl_map_pep_seqid_thr_mdl_map_pep_seqid_thr = mdl_map_pep_seqid_thr
487  self._mdl_map_nuc_seqid_thr_mdl_map_nuc_seqid_thr = mdl_map_nuc_seqid_thr
488 
489  # lazily computed attributes
490  self.__chain_mapper__chain_mapper = None
491  self.__chem_mapping__chem_mapping = None
492  self.__chem_group_alns__chem_group_alns = None
493  self.__ref_mdl_alns__ref_mdl_alns = None
494  self.__unmapped_mdl_chains__unmapped_mdl_chains = None
495  self.__chain_mapping_mdl__chain_mapping_mdl = None
496 
497  # keep track of states
498  # simple integers instead of enums - encoding available in
499  # self.state_decoding
500  self._state_matrix_state_matrix = None
501  self._model_ligand_states_model_ligand_states = None
502  self._target_ligand_states_target_ligand_states = None
503 
504  # score matrices
505  self._score_matrix_score_matrix = None
506  self._coverage_matrix_coverage_matrix = None
507  self._aux_matrix_aux_matrix = None
508 
509  # assignment and derived data
510  self._assignment_assignment = None
511  self._score_dict_score_dict = None
512  self._aux_dict_aux_dict = None
513 
514  # human readable description of states - child class must extend with
515  # with child class specific states
516  # each state code comes with a tuple of two elements:
517  # 1) short description 2) human readable description
518  # The actual states are set in _compute_scores in :class:`LigandScorer`
519  # or _compute_score of the child class.
520  if self.substructure_matchsubstructure_match:
521  iso = "subgraph isomorphism"
522  else:
523  iso = "full graph isomorphism"
524 
525  self.state_decodingstate_decoding = \
526  {0: ("OK", "OK"),
527  1: ("identity", f"Ligands could not be matched (by {iso})"),
528  2: ("symmetries", "Too many symmetries between ligand atoms were "
529  "found - increasing max_symmetries might help"),
530  3: ("no_iso", "No full isomorphic match could be found - enabling "
531  "substructure_match might allow a match"),
532  4: ("disconnected", "Ligand graph is disconnected"),
533  5: ("stoichiometry", "Ligand was already assigned to another ligand "
534  "(different stoichiometry)"),
535  6: ("single_ligand_issue", "Cannot compute valid pairwise score as "
536  "either model or target ligand have non-zero state."),
537  9: ("unknown", "An unknown error occured in LigandScorer")}
538 
539 
540  @property
541  def model(self):
542  """ Model receptor structure
543 
544  Processed according to docs in :class:`LigandScorer` constructor
545  """
546  return self._model_model
547 
548  @property
549  def target(self):
550  """ Target receptor structure
551 
552  Processed according to docs in :class:`LigandScorer` constructor
553  """
554  return self._target_target
555 
556  @property
557  def model_cleanup_log(self):
558  """ Reports residues/atoms that were removed in model during cleanup
559 
560  Residues and atoms are described as :class:`str` in format
561  <chain_name>.<resnum>.<ins_code> (residue) and
562  <chain_name>.<resnum>.<ins_code>.<aname> (atom).
563 
564  :class:`dict` with keys:
565 
566  * 'cleaned_residues': another :class:`dict` with keys:
567 
568  * 'no_clib': residues that have been removed because no entry could be
569  found in :class:`ost.conop.CompoundLib`
570  * 'not_linking': residues that have been removed because they're not
571  peptide or nucleotide linking according to
572  :class:`ost.conop.CompoundLib`
573 
574  * 'cleaned_atoms': another :class:`dict` with keys:
575 
576  * 'unknown_atoms': atoms that have been removed as they're not part
577  of their respective residue according to
578  :class:`ost.conop.CompoundLib`
579  """
580  return self._model_cleanup_log_model_cleanup_log
581 
582  @property
584  """ Same for target
585  """
586  return self._target_cleanup_log_target_cleanup_log
587 
588  @property
589  def model_ligands(self):
590  """ Residues representing model ligands
591 
592  :class:`list` of :class:`ost.mol.ResidueHandle`
593  """
594  return self._model_ligands_model_ligands
595 
596  @property
597  def target_ligands(self):
598  """ Residues representing target ligands
599 
600  :class:`list` of :class:`ost.mol.ResidueHandle`
601  """
602  return self._target_ligands_target_ligands
603 
604  @property
605  def resnum_alignments(self):
606  """ Given at :class:`LigandScorer` construction
607  """
608  return self._resnum_alignments_resnum_alignments
609 
610  @property
611  def min_pep_length(self):
612  """ Given at :class:`LigandScorer` construction
613  """
614  return self._min_pep_length_min_pep_length
615 
616  @property
617  def min_nuc_length(self):
618  """ Given at :class:`LigandScorer` construction
619  """
620  return self._min_nuc_length_min_nuc_length
621 
622  @property
623  def pep_seqid_thr(self):
624  """ Given at :class:`LigandScorer` construction
625  """
626  return self._pep_seqid_thr_pep_seqid_thr
627 
628  @property
629  def nuc_seqid_thr(self):
630  """ Given at :class:`LigandScorer` construction
631  """
632  return self._nuc_seqid_thr_nuc_seqid_thr
633 
634  @property
636  """ Given at :class:`LigandScorer` construction
637  """
638  return self._mdl_map_pep_seqid_thr_mdl_map_pep_seqid_thr
639 
640  @property
642  """ Given at :class:`LigandScorer` construction
643  """
644  return self._mdl_map_nuc_seqid_thr_mdl_map_nuc_seqid_thr
645 
646  @property
648  """ Given at :class:`LigandScorer` construction
649  """
650  return self._substructure_match_substructure_match
651 
652  @property
653  def coverage_delta(self):
654  """ Given at :class:`LigandScorer` construction
655  """
656  return self._coverage_delta_coverage_delta
657 
658  @property
659  def max_symmetries(self):
660  """ Given at :class:`LigandScorer` construction
661  """
662  return self._max_symmetries_max_symmetries
663 
664  @property
665  def state_matrix(self):
666  """ Encodes states of ligand pairs
667 
668  Ligand pairs can be matched and a valid score can be expected if
669  respective location in this matrix is 0.
670  Target ligands are in rows, model ligands in columns. States are encoded
671  as integers <= 9. Larger numbers encode errors for child classes.
672  Use something like ``self.state_decoding[3]`` to get a decscription.
673 
674  :rtype: :class:`~numpy.ndarray`
675  """
676  if self._state_matrix_state_matrix is None:
677  self._compute_scores_compute_scores()
678  return self._state_matrix_state_matrix
679 
680  @property
682  """ Encodes states of model ligands
683 
684  Non-zero state in any of the model ligands invalidates the full
685  respective column in :attr:`~state_matrix`.
686 
687  :rtype: :class:`~numpy.ndarray`
688  """
689  if self._model_ligand_states_model_ligand_states is None:
690  self._compute_scores_compute_scores()
691  return self._model_ligand_states_model_ligand_states
692 
693  @property
695  """ Encodes states of target ligands
696 
697  Non-zero state in any of the target ligands invalidates the full
698  respective row in :attr:`~state_matrix`.
699 
700  :rtype: :class:`~numpy.ndarray`
701  """
702  if self._target_ligand_states_target_ligand_states is None:
703  self._compute_scores_compute_scores()
704  return self._target_ligand_states_target_ligand_states
705 
706  @property
707  def score_matrix(self):
708  """ Get the matrix of scores.
709 
710  Target ligands are in rows, model ligands in columns.
711 
712  NaN values indicate that no value could be computed (i.e. different
713  ligands). In other words: values are only valid if the respective
714  location in :attr:`~state_matrix` is 0.
715 
716  :rtype: :class:`~numpy.ndarray`
717  """
718  if self._score_matrix_score_matrix is None:
719  self._compute_scores_compute_scores()
720  return self._score_matrix_score_matrix
721 
722  @property
723  def coverage_matrix(self):
724  """ Get the matrix of model ligand atom coverage in the target.
725 
726  Target ligands are in rows, model ligands in columns.
727 
728  NaN values indicate that no value could be computed (i.e. different
729  ligands). In other words: values are only valid if the respective
730  location in :attr:`~state_matrix` is 0. If `substructure_match=False`,
731  only full match isomorphisms are considered, and therefore only values
732  of 1.0 can be observed.
733 
734  :rtype: :class:`~numpy.ndarray`
735  """
736  if self._coverage_matrix_coverage_matrix is None:
737  self._compute_scores_compute_scores()
738  return self._coverage_matrix_coverage_matrix
739 
740  @property
741  def aux_matrix(self):
742  """ Get the matrix of scorer specific auxiliary data.
743 
744  Target ligands are in rows, model ligands in columns.
745 
746  Auxiliary data consists of arbitrary data dicts which allow a child
747  class to provide additional information for a scored ligand pair.
748  empty dictionaries indicate that the child class simply didn't return
749  anything or that no value could be computed (e.g. different ligands).
750  In other words: values are only valid if respective location in the
751  :attr:`~state_matrix` is 0.
752 
753  :rtype: :class:`~numpy.ndarray`
754  """
755  if self._aux_matrix_aux_matrix is None:
756  self._compute_scores_compute_scores()
757  return self._aux_matrix_aux_matrix
758 
759  @property
760  def assignment(self):
761  """ Ligand assignment based on computed scores
762 
763  Implements a greedy algorithm to assign target and model ligands
764  with each other. Starts from each valid ligand pair as indicated
765  by a state of 0 in :attr:`state_matrix`. Each iteration first selects
766  high coverage pairs. Given max_coverage defined as the highest
767  coverage observed in the available pairs, all pairs with coverage
768  in [max_coverage-*coverage_delta*, max_coverage] are selected.
769  The best scoring pair among those is added to the assignment
770  and the whole process is repeated until there are no ligands to
771  assign anymore.
772 
773  :rtype: :class:`list` of :class:`tuple` (trg_lig_idx, mdl_lig_idx)
774  """
775  if self._assignment_assignment is None:
776  self._assignment_assignment = list()
777  # Build working array that contains tuples for all mdl/trg ligand
778  # pairs with valid score as indicated by a state of 0:
779  # (score, coverage, trg_ligand_idx, mdl_ligand_idx)
780  tmp = list()
781  for trg_idx in range(self.score_matrixscore_matrix.shape[0]):
782  for mdl_idx in range(self.score_matrixscore_matrix.shape[1]):
783  if self.state_matrixstate_matrix[trg_idx, mdl_idx] == 0:
784  tmp.append((self.score_matrixscore_matrix[trg_idx, mdl_idx],
785  self.coverage_matrixcoverage_matrix[trg_idx, mdl_idx],
786  trg_idx, mdl_idx))
787 
788  # sort by score, such that best scoring item is in front
789  if self._score_dir_score_dir() == '+':
790  tmp.sort(reverse=True)
791  elif self._score_dir_score_dir() == '-':
792  tmp.sort()
793  else:
794  raise RuntimeError("LigandScorer._score_dir must return one in "
795  "['+', '-']")
796 
797  LogScript("Computing ligand assignment")
798  while len(tmp) > 0:
799  # select high coverage ligand pairs in working array
800  coverage_thresh = max([x[1] for x in tmp]) - self.coverage_deltacoverage_delta
801  top_coverage = [x for x in tmp if x[1] >= coverage_thresh]
802 
803  # working array is sorted by score => just pick first one
804  a = top_coverage[0][2] # selected trg_ligand_idx
805  b = top_coverage[0][3] # selected mdl_ligand_idx
806  self._assignment_assignment.append((a, b))
807 
808  # kick out remaining pairs involving these ligands
809  tmp = [x for x in tmp if (x[2] != a and x[3] != b)]
810 
811  return self._assignment_assignment
812 
813  @property
814  def score(self):
815  """ Get a dictionary of score values, keyed by model ligand
816 
817  Extract score with something like:
818  ``scorer.score[lig.GetChain().GetName()][lig.GetNumber()]``.
819  The returned scores are based on :attr:`~assignment`.
820 
821  :rtype: :class:`dict`
822  """
823  if self._score_dict_score_dict is None:
824  self._score_dict_score_dict = dict()
825  for (trg_lig_idx, mdl_lig_idx) in self.assignmentassignment:
826  mdl_lig = self.model_ligandsmodel_ligands[mdl_lig_idx]
827  cname = mdl_lig.GetChain().GetName()
828  rnum = mdl_lig.GetNumber()
829  if cname not in self._score_dict_score_dict:
830  self._score_dict_score_dict[cname] = dict()
831  score = self.score_matrixscore_matrix[trg_lig_idx, mdl_lig_idx]
832  self._score_dict_score_dict[cname][rnum] = score
833  return self._score_dict_score_dict
834 
835  @property
836  def aux(self):
837  """ Get a dictionary of score details, keyed by model ligand
838 
839  Extract dict with something like:
840  ``scorer.score[lig.GetChain().GetName()][lig.GetNumber()]``.
841  The returned info dicts are based on :attr:`~assignment`. The content is
842  documented in the respective child class.
843 
844  :rtype: :class:`dict`
845  """
846  if self._aux_dict_aux_dict is None:
847  self._aux_dict_aux_dict = dict()
848  for (trg_lig_idx, mdl_lig_idx) in self.assignmentassignment:
849  mdl_lig = self.model_ligandsmodel_ligands[mdl_lig_idx]
850  cname = mdl_lig.GetChain().GetName()
851  rnum = mdl_lig.GetNumber()
852  if cname not in self._aux_dict_aux_dict:
853  self._aux_dict_aux_dict[cname] = dict()
854  d = self.aux_matrixaux_matrix[trg_lig_idx, mdl_lig_idx]
855  self._aux_dict_aux_dict[cname][rnum] = d
856  return self._aux_dict_aux_dict
857 
858 
859  @property
861  """ Get indices of target ligands which are not assigned
862 
863  :rtype: :class:`list` of :class:`int`
864  """
865  # compute on-the-fly, no need for caching
866  assigned = set([x[0] for x in self.assignmentassignment])
867  return [x for x in range(len(self.target_ligandstarget_ligands)) if x not in assigned]
868 
869  @property
871  """ Get indices of model ligands which are not assigned
872 
873  :rtype: :class:`list` of :class:`int`
874  """
875  # compute on-the-fly, no need for caching
876  assigned = set([x[1] for x in self.assignmentassignment])
877  return [x for x in range(len(self.model_ligandsmodel_ligands)) if x not in assigned]
878 
879  def get_target_ligand_state_report(self, trg_lig_idx):
880  """ Get summary of states observed with respect to all model ligands
881 
882  Mainly for debug purposes
883 
884  :param trg_lig_idx: Index of target ligand for which report should be
885  generated
886  :type trg_lig_idx: :class:`int`
887  """
888  return self._get_report_get_report(self.target_ligand_statestarget_ligand_states[trg_lig_idx],
889  self.state_matrixstate_matrix[trg_lig_idx,:])
890 
891  def get_model_ligand_state_report(self, mdl_lig_idx):
892  """ Get summary of states observed with respect to all target ligands
893 
894  Mainly for debug purposes
895 
896  :param mdl_lig_idx: Index of model ligand for which report should be
897  generated
898  :type mdl_lig_idx: :class:`int`
899  """
900  return self._get_report_get_report(self.model_ligand_statesmodel_ligand_states[mdl_lig_idx],
901  self.state_matrixstate_matrix[:, mdl_lig_idx])
902 
903  def _get_report(self, ligand_state, pair_states):
904  """ Helper
905  """
906  pair_report = list()
907  for s in np.unique(pair_states):
908  desc = self.state_decodingstate_decoding[s]
909  indices = np.flatnonzero(pair_states == s).tolist()
910  pair_report.append({"state": s,
911  "short desc": desc[0],
912  "desc": desc[1],
913  "indices": indices})
914 
915  desc = self.state_decodingstate_decoding[ligand_state]
916  ligand_report = {"state": ligand_state,
917  "short desc": desc[0],
918  "desc": desc[1]}
919 
920  return (ligand_report, pair_report)
921 
923  """ Makes an educated guess why target ligand is not assigned
924 
925  This either returns actual error states or custom states that are
926  derived from them. Currently, the following reasons are reported:
927 
928  * `no_ligand`: there was no ligand in the model.
929  * `disconnected`: the ligand graph was disconnected.
930  * `identity`: the ligand was not found in the model (by graph
931  isomorphism). Check your ligand connectivity.
932  * `no_iso`: no full isomorphic match could be found. Try enabling
933  `substructure_match=True` if the target ligand is incomplete.
934  * `symmetries`: too many symmetries were found (by graph isomorphisms).
935  Try to increase `max_symmetries`.
936  * `stoichiometry`: there was a possible assignment in the model, but
937  the model ligand was already assigned to a different target ligand.
938  This indicates different stoichiometries.
939  * `no_contact` (LDDT-PLI only): There were no LDDT contacts between
940  the binding site and the ligand, and LDDT-PLI is undefined.
941  * `target_binding_site` (SCRMSD only): no polymer residues were in
942  proximity of the target ligand.
943  * `model_binding_site` (SCRMSD only): the binding site was not found
944  in the model. Either the binding site was not modeled or the model
945  ligand was positioned too far in combination with
946  `full_bs_search=False`.
947 
948  :param trg_lig_idx: Index of target ligand
949  :type trg_lig_idx: :class:`int`
950  :returns: :class:`tuple` with two elements: 1) keyword 2) human readable
951  sentence describing the issue, (\"unknown\",\"unknown\") if
952  nothing obvious can be found.
953  :raises: :class:`RuntimeError` if specified target ligand is assigned
954  """
955  if trg_lig_idx not in self.unassigned_target_ligandsunassigned_target_ligands:
956  raise RuntimeError("Specified target ligand is not unassigned")
957 
958  # hardcoded tuple if there is simply nothing we can assign it to
959  if len(self.model_ligandsmodel_ligands) == 0:
960  return ("no_ligand", "No ligand in the model")
961 
962  # if something with the ligand itself is wrong, we can be pretty sure
963  # thats why the ligand is unassigned
964  if self.target_ligand_statestarget_ligand_states[trg_lig_idx] != 0:
965  return self.state_decodingstate_decoding[self.target_ligand_statestarget_ligand_states[trg_lig_idx]]
966 
967  # The next best guess comes from looking at pair states
968  tmp = np.unique(self.state_matrixstate_matrix[trg_lig_idx,:])
969 
970  # In case of any 0, it could have been assigned so it's probably
971  # just not selected due to different stoichiometry - this is no
972  # defined state, we just return a hardcoded tuple in this case
973  if 0 in tmp:
974  return ("stoichiometry",
975  "Ligand was already assigned to an other "
976  "model ligand (different stoichiometry)")
977 
978  # maybe it's a symmetry issue?
979  if 2 in tmp:
980  return self.state_decodingstate_decoding[2]
981 
982  # if the state is 6 (single_ligand_issue), there is an issue with its
983  # target counterpart.
984  if 6 in tmp:
985  mdl_idx = np.where(self.state_matrixstate_matrix[trg_lig_idx,:]==6)[0]
986  for i in mdl_idx:
987  if self.model_ligand_statesmodel_ligand_states[i] == 0:
988  raise RuntimeError("This should never happen")
989  if self.model_ligand_statesmodel_ligand_states[i] != 4 or len(tmp) == 1:
990  # Don't report disconnected if only 1 model ligand is
991  # disconnected, unless that's the only reason
992  return self.state_decodingstate_decoding[self.model_ligand_statesmodel_ligand_states[i]]
993 
994  # get rid of remaining single ligand issues (only disconnected error)
995  if 6 in tmp and len(tmp) > 1:
996  tmp = tmp[tmp!=6]
997 
998  # prefer everything over identity state
999  if 1 in tmp and len(tmp) > 1:
1000  tmp = tmp[tmp!=1]
1001 
1002  # just return whatever is left
1003  return self.state_decodingstate_decoding[tmp[0]]
1004 
1005 
1007  """ Makes an educated guess why model ligand is not assigned
1008 
1009  This either returns actual error states or custom states that are
1010  derived from them. Currently, the following reasons are reported:
1011 
1012  * `no_ligand`: there was no ligand in the target.
1013  * `disconnected`: the ligand graph is disconnected.
1014  * `identity`: the ligand was not found in the target (by graph or
1015  subgraph isomorphism). Check your ligand connectivity.
1016  * `no_iso`: no full isomorphic match could be found. Try enabling
1017  `substructure_match=True` if the target ligand is incomplete.
1018  * `symmetries`: too many symmetries were found (by graph isomorphisms).
1019  Try to increase `max_symmetries`.
1020  * `stoichiometry`: there was a possible assignment in the target, but
1021  the model target was already assigned to a different model ligand.
1022  This indicates different stoichiometries.
1023  * `no_contact` (LDDT-PLI only): There were no LDDT contacts between
1024  the binding site and the ligand, and LDDT-PLI is undefined.
1025  * `target_binding_site` (SCRMSD only): a potential assignment was found
1026  in the target, but there were no polymer residues in proximity of the
1027  ligand in the target.
1028  * `model_binding_site` (SCRMSD only): a potential assignment was
1029  found in the target, but no binding site was found in the model.
1030  Either the binding site was not modeled or the model ligand was
1031  positioned too far in combination with `full_bs_search=False`.
1032 
1033  :param mdl_lig_idx: Index of model ligand
1034  :type mdl_lig_idx: :class:`int`
1035  :returns: :class:`tuple` with two elements: 1) keyword 2) human readable
1036  sentence describing the issue, (\"unknown\",\"unknown\") if
1037  nothing obvious can be found.
1038  :raises: :class:`RuntimeError` if specified model ligand is assigned
1039  """
1040  if mdl_lig_idx not in self.unassigned_model_ligandsunassigned_model_ligands:
1041  raise RuntimeError("Specified model ligand is not unassigned")
1042 
1043  # hardcoded tuple if there is simply nothing we can assign it to
1044  if len(self.target_ligandstarget_ligands) == 0:
1045  return ("no_ligand", "No ligand in the target")
1046 
1047  # if something with the ligand itself is wrong, we can be pretty sure
1048  # thats why the ligand is unassigned
1049  if self.model_ligand_statesmodel_ligand_states[mdl_lig_idx] != 0:
1050  return self.state_decodingstate_decoding[self.model_ligand_statesmodel_ligand_states[mdl_lig_idx]]
1051 
1052  # The next best guess comes from looking at pair states
1053  tmp = np.unique(self.state_matrixstate_matrix[:,mdl_lig_idx])
1054 
1055  # In case of any 0, it could have been assigned so it's probably
1056  # just not selected due to different stoichiometry - this is no
1057  # defined state, we just return a hardcoded tuple in this case
1058  if 0 in tmp:
1059  return ("stoichiometry",
1060  "Ligand was already assigned to an other "
1061  "model ligand (different stoichiometry)")
1062 
1063  # maybe its a symmetry issue?
1064  if 2 in tmp:
1065  return self.state_decodingstate_decoding[2]
1066 
1067  # if the state is 6 (single_ligand_issue), there is an issue with its
1068  # target counterpart.
1069  if 6 in tmp:
1070  trg_idx = np.where(self.state_matrixstate_matrix[:,mdl_lig_idx]==6)[0]
1071  for i in trg_idx:
1072  if self.target_ligand_statestarget_ligand_states[i] == 0:
1073  raise RuntimeError("This should never happen")
1074  if self.target_ligand_statestarget_ligand_states[i] != 4 or len(tmp) == 1:
1075  # Don't report disconnected if only 1 target ligand is
1076  # disconnected, unless that's the only reason
1077  return self.state_decodingstate_decoding[self.target_ligand_statestarget_ligand_states[i]]
1078 
1079  # get rid of remaining single ligand issues (only disconnected error)
1080  if 6 in tmp and len(tmp) > 1:
1081  tmp = tmp[tmp!=6]
1082 
1083  # prefer everything over identity state
1084  if 1 in tmp and len(tmp) > 1:
1085  tmp = tmp[tmp!=1]
1086 
1087  # just return whatever is left
1088  return self.state_decodingstate_decoding[tmp[0]]
1089 
1090  @property
1092  return_dict = dict()
1093  for i in self.unassigned_model_ligandsunassigned_model_ligands:
1094  lig = self.model_ligandsmodel_ligands[i]
1095  cname = lig.GetChain().GetName()
1096  rnum = lig.GetNumber()
1097  if cname not in return_dict:
1098  return_dict[cname] = dict()
1099  return_dict[cname][rnum] = \
1100  self.guess_model_ligand_unassigned_reasonguess_model_ligand_unassigned_reason(i)[0]
1101  return return_dict
1102 
1103  @property
1105  return_dict = dict()
1106  for i in self.unassigned_target_ligandsunassigned_target_ligands:
1107  lig = self.target_ligandstarget_ligands[i]
1108  cname = lig.GetChain().GetName()
1109  rnum = lig.GetNumber()
1110  if cname not in return_dict:
1111  return_dict[cname] = dict()
1112  return_dict[cname][rnum] = \
1113  self.guess_target_ligand_unassigned_reasonguess_target_ligand_unassigned_reason(i)[0]
1114  return return_dict
1115 
1116  @property
1117  def _chain_mapper(self):
1118  """ Chain mapper object for the given :attr:`target`.
1119 
1120  Can be used by child classes if needed, constructed with
1121  *resnum_alignments* flag
1122 
1123  :type: :class:`ost.mol.alg.chain_mapping.ChainMapper`
1124  """
1125  if self.__chain_mapper__chain_mapper is None:
1126  with _SinkVerbosityLevel():
1127  self.__chain_mapper__chain_mapper = \
1128  chain_mapping.ChainMapper(self.targettarget,
1129  n_max_naive=1e9,
1130  resnum_alignments=self.resnum_alignmentsresnum_alignments,
1131  min_pep_length=self.min_pep_lengthmin_pep_length,
1132  min_nuc_length=self.min_nuc_lengthmin_nuc_length,
1133  pep_seqid_thr=self.pep_seqid_thrpep_seqid_thr,
1134  nuc_seqid_thr=self.nuc_seqid_thrnuc_seqid_thr,
1135  mdl_map_pep_seqid_thr=self.mdl_map_pep_seqid_thrmdl_map_pep_seqid_thr,
1136  mdl_map_nuc_seqid_thr=self.mdl_map_nuc_seqid_thrmdl_map_nuc_seqid_thr)
1137  return self.__chain_mapper__chain_mapper
1138 
1139  @property
1140  def _chem_mapping(self):
1141  if self.__chem_mapping__chem_mapping is None:
1142  self.__chem_mapping__chem_mapping, self.__chem_group_alns__chem_group_alns, \
1143  self.__unmapped_mdl_chains__unmapped_mdl_chains, self.__chain_mapping_mdl__chain_mapping_mdl = \
1144  self._chain_mapper_chain_mapper.GetChemMapping(self.modelmodel)
1145  return self.__chem_mapping__chem_mapping
1146 
1147  @property
1148  def _chem_group_alns(self):
1149  if self.__chem_group_alns__chem_group_alns is None:
1150  self.__chem_mapping__chem_mapping, self.__chem_group_alns__chem_group_alns, \
1151  self.__unmapped_mdl_chains__unmapped_mdl_chains, self.__chain_mapping_mdl__chain_mapping_mdl = \
1152  self._chain_mapper_chain_mapper.GetChemMapping(self.modelmodel)
1153  return self.__chem_group_alns__chem_group_alns
1154 
1155  @property
1156  def _ref_mdl_alns(self):
1157  if self.__ref_mdl_alns__ref_mdl_alns is None:
1158  self.__ref_mdl_alns__ref_mdl_alns = \
1159  chain_mapping._GetRefMdlAlns(self._chain_mapper_chain_mapper.chem_groups,
1160  self._chain_mapper_chain_mapper.chem_group_alignments,
1161  self._chem_mapping_chem_mapping,
1162  self._chem_group_alns_chem_group_alns)
1163  return self.__ref_mdl_alns__ref_mdl_alns
1164 
1165  @property
1166  def _chain_mapping_mdl(self):
1167  if self.__chain_mapping_mdl__chain_mapping_mdl is None:
1168  with _SinkVerbosityLevel():
1169  self.__chem_mapping__chem_mapping, self.__chem_group_alns__chem_group_alns, \
1170  self.__unmapped_mdl_chains__unmapped_mdl_chains, self.__chain_mapping_mdl__chain_mapping_mdl = \
1171  self._chain_mapper_chain_mapper.GetChemMapping(self.modelmodel)
1172  return self.__chain_mapping_mdl__chain_mapping_mdl
1173 
1174  @property
1175  def _unmapped_mdl_chains(self):
1176  if self.__unmapped_mdl_chains__unmapped_mdl_chains is None:
1177  self.__chem_mapping__chem_mapping, self.__chem_group_alns__chem_group_alns, \
1178  self.__unmapped_mdl_chains__unmapped_mdl_chains, self.__chain_mapping_mdl__chain_mapping_mdl = \
1179  self._chain_mapper_chain_mapper.GetChemMapping(self.modelmodel)
1180  return self.__unmapped_mdl_chains__unmapped_mdl_chains
1181 
1182  def _compute_scores(self):
1183  """
1184  Compute score for every possible target-model ligand pair and store the
1185  result in internal matrices.
1186  """
1187 
1190  shape = (len(self.target_ligandstarget_ligands), len(self.model_ligandsmodel_ligands))
1191  self._score_matrix_score_matrix = np.full(shape, np.nan, dtype=np.float32)
1192  self._coverage_matrix_coverage_matrix = np.full(shape, np.nan, dtype=np.float32)
1193  self._state_matrix_state_matrix = np.full(shape, 0, dtype=np.int32)
1194  self._model_ligand_states_model_ligand_states = np.zeros(len(self.model_ligandsmodel_ligands))
1195  self._target_ligand_states_target_ligand_states = np.zeros(len(self.target_ligandstarget_ligands))
1196  self._aux_matrix_aux_matrix = np.empty(shape, dtype=dict)
1197 
1198  # precompute ligand graphs
1199  target_graphs = \
1200  [_ResidueToGraph(l, by_atom_index=True) for l in self.target_ligandstarget_ligands]
1201  model_graphs = \
1202  [_ResidueToGraph(l, by_atom_index=True) for l in self.model_ligandsmodel_ligands]
1203 
1204  for g_idx, g in enumerate(target_graphs):
1205  if not networkx.is_connected(g):
1206  self._target_ligand_states_target_ligand_states[g_idx] = 4
1207  self._state_matrix_state_matrix[g_idx,:] = 6
1208  msg = "Disconnected graph observed for target ligand "
1209  msg += str(self.target_ligandstarget_ligands[g_idx])
1210  LogWarning(msg)
1211 
1212  for g_idx, g in enumerate(model_graphs):
1213  if not networkx.is_connected(g):
1214  self._model_ligand_states_model_ligand_states[g_idx] = 4
1215  self._state_matrix_state_matrix[:,g_idx] = 6
1216  msg = "Disconnected graph observed for model ligand "
1217  msg += str(self.model_ligandsmodel_ligands[g_idx])
1218  LogWarning(msg)
1219 
1220 
1221  LogScript("Computing pairwise scores for all %s x %s ligands" % shape)
1222  for target_id, target_ligand in enumerate(self.target_ligandstarget_ligands):
1223  LogInfo("Analyzing target ligand %s" % target_ligand)
1224 
1225  if self._target_ligand_states_target_ligand_states[target_id] == 4:
1226  # Disconnected graph - already updated states and reported
1227  # to LogVerbose
1228  continue
1229 
1230  for model_id, model_ligand in enumerate(self.model_ligandsmodel_ligands):
1231  LogInfo("Comparing to model ligand %s" % model_ligand)
1232 
1233 
1236 
1237  if self._model_ligand_states_model_ligand_states[model_id] == 4:
1238  # Disconnected graph - already updated states and reported
1239  # to LogVerbose
1240  continue
1241 
1242  try:
1243  symmetries = ComputeSymmetries(
1244  model_ligand, target_ligand,
1245  substructure_match=self.substructure_matchsubstructure_match,
1246  by_atom_index=True,
1247  max_symmetries=self.max_symmetriesmax_symmetries,
1248  model_graph = model_graphs[model_id],
1249  target_graph = target_graphs[target_id])
1250  LogInfo("Ligands %s and %s symmetry match" % (
1251  str(model_ligand), str(target_ligand)))
1252  except NoSymmetryError:
1253  # Ligands are different - skip
1254  LogInfo("No symmetry between %s and %s" % (
1255  str(model_ligand), str(target_ligand)))
1256  self._state_matrix_state_matrix[target_id, model_id] = 1
1257  continue
1258  except TooManySymmetriesError:
1259  # Ligands are too symmetrical - skip
1260  LogWarning("Too many symmetries between %s and %s" % (
1261  str(model_ligand), str(target_ligand)))
1262  self._state_matrix_state_matrix[target_id, model_id] = 2
1263  continue
1264  except NoIsomorphicSymmetryError:
1265  # Ligands are different - skip
1266  LogInfo("No isomorphic symmetry between %s and %s" % (
1267  str(model_ligand), str(target_ligand)))
1268  self._state_matrix_state_matrix[target_id, model_id] = 3
1269  continue
1270  except DisconnectedGraphError:
1271  # this should never happen as we guard against
1272  # DisconnectedGraphError when precomputing the graph
1273  LogError("Disconnected graph observed for %s and %s" % (
1274  str(model_ligand), str(target_ligand)))
1275  # just set both ligand states to 4
1276  self._model_ligand_states_model_ligand_states[model_id] = 4
1277  self._model_ligand_states_model_ligand_states[target_id] = 4
1278  self._state_matrix_state_matrix[target_id, model_id] = 6
1279  continue
1280 
1281 
1284  score, pair_state, target_ligand_state, model_ligand_state, aux\
1285  = self._compute_compute(symmetries, target_ligand, model_ligand)
1286 
1287 
1290 
1291  # Ensure that returned states are associated with a
1292  # description. This is a requirement when subclassing
1293  # LigandScorer => state_decoding dict from base class must
1294  # be modified in subclass constructor
1295  if pair_state not in self.state_decodingstate_decoding or \
1296  target_ligand_state not in self.state_decodingstate_decoding or \
1297  model_ligand_state not in self.state_decodingstate_decoding:
1298  raise RuntimeError(f"Subclass returned state "
1299  f"\"{state}\" for which no "
1300  f"description is available. Point "
1301  f"the developer of the used scorer "
1302  f"to this error message.")
1303 
1304  # if any of the ligand states is non-zero, this must be marked
1305  # by the scorer subclass by specifying a specific pair state
1306  if target_ligand_state != 0 and pair_state != 6:
1307  raise RuntimeError("Observed non-zero target-ligand "
1308  "state which must trigger certain "
1309  "pair state. Point the developer of "
1310  "the used scorer to this error message")
1311 
1312  if model_ligand_state != 0 and pair_state != 6:
1313  raise RuntimeError("Observed non-zero model-ligand "
1314  "state which must trigger certain "
1315  "pair state. Point the developer of "
1316  "the used scorer to this error message")
1317 
1318  self._state_matrix_state_matrix[target_id, model_id] = pair_state
1319  self._target_ligand_states_target_ligand_states[target_id] = target_ligand_state
1320  self._model_ligand_states_model_ligand_states[model_id] = model_ligand_state
1321  if pair_state == 0:
1322  if score is None or np.isnan(score):
1323  raise RuntimeError("LigandScorer returned invalid "
1324  "score despite 0 error state")
1325  # it's a valid score!
1326  self._score_matrix_score_matrix[target_id, model_id] = score
1327  cvg = len(symmetries[0][0]) / len(model_ligand.atoms)
1328  self._coverage_matrix_coverage_matrix[target_id, model_id] = cvg
1329  self._aux_matrix_aux_matrix[target_id, model_id] = aux
1330 
1331  def _compute(self, symmetries, target_ligand, model_ligand):
1332  """ Compute score for specified ligand pair - defined by child class
1333 
1334  Raises :class:`NotImplementedError` if not implemented by child class.
1335 
1336  :param symmetries: Defines symmetries between *target_ligand* and
1337  *model_ligand*. Return value of
1338  :func:`ComputeSymmetries`
1339  :type symmetries: :class:`list` of :class:`tuple` with two elements
1340  each: 1) :class:`list` of atom indices in
1341  *target_ligand* 2) :class:`list` of respective atom
1342  indices in *model_ligand*
1343  :param target_ligand: The target ligand
1344  :type target_ligand: :class:`ost.mol.ResidueHandle` or
1345  :class:`ost.mol.ResidueView`
1346  :param model_ligand: The model ligand
1347  :type model_ligand: :class:`ost.mol.ResidueHandle` or
1348  :class:`ost.mol.ResidueView`
1349 
1350  :returns: A :class:`tuple` with three elements: 1) a score
1351  (:class:`float`) 2) state (:class:`int`).
1352  3) auxiliary data for this ligand pair (:class:`dict`).
1353  If state is 0, the score and auxiliary data will be
1354  added to :attr:`~score_matrix` and :attr:`~aux_matrix` as well
1355  as the respective value in :attr:`~coverage_matrix`.
1356  Returned score must be valid in this case (not None/NaN).
1357  Child specific non-zero states must be >= 10.
1358  """
1359  raise NotImplementedError("_compute must be implemented by child "
1360  "class")
1361 
1362  def _score_dir(self):
1363  """ Return direction of score - defined by child class
1364 
1365  Relevant for ligand assignment. Must return a string in ['+', '-'].
1366  '+' for ascending scores, i.e. higher is better (lddt etc.)
1367  '-' for descending scores, i.e. lower is better (rmsd etc.)
1368  """
1369  raise NotImplementedError("_score_dir must be implemented by child "
1370  "class")
1371 
1372  def _copy_ligand(self, l, ent, ed, rename_ligand_chain):
1373  """ Copies ligand into entity and returns residue handle
1374  """
1375  new_chain = ent.FindChain(l.chain.name)
1376  if not new_chain.IsValid():
1377  new_chain = ed.InsertChain(l.chain.name)
1378  ed.SetChainType(new_chain, mol.ChainType.CHAINTYPE_NON_POLY)
1379  else:
1380  # Does a residue with the same name already exist?
1381  already_exists = new_chain.FindResidue(l.number).IsValid()
1382  if already_exists:
1383  if rename_ligand_chain:
1384  chain_ext = 2 # Extend the chain name by this
1385  while True:
1386  new_chain_name = \
1387  l.chain.name + "_" + str(chain_ext)
1388  new_chain = ent.FindChain(new_chain_name)
1389  if new_chain.IsValid():
1390  chain_ext += 1
1391  continue
1392  else:
1393  new_chain = \
1394  ed.InsertChain(new_chain_name)
1395  break
1396  LogInfo("Moved ligand residue %s to new chain %s" % (
1397  l.qualified_name, new_chain.name))
1398  else:
1399  msg = \
1400  "A residue number %s already exists in chain %s" % (
1401  l.number, l.chain.name)
1402  raise RuntimeError(msg)
1403 
1404  # Add the residue with its original residue number
1405  new_res = ed.AppendResidue(new_chain, l.name, l.number)
1406  # Add atoms
1407  for old_atom in l.atoms:
1408  ed.InsertAtom(new_res, old_atom.name, old_atom.pos,
1409  element=old_atom.element, occupancy=old_atom.occupancy,
1410  b_factor=old_atom.b_factor, is_hetatm=old_atom.is_hetatom)
1411  # Add bonds
1412  for old_atom in l.atoms:
1413  for old_bond in old_atom.bonds:
1414  new_first = new_res.FindAtom(old_bond.first.name)
1415  new_second = new_res.FindAtom(old_bond.second.name)
1416  ed.Connect(new_first, new_second, old_bond.bond_order)
1417  new_res.SetIsLigand(True)
1418  return new_res
1419 
1420  def _cleanup_polymer_ent(self, ent, clib):
1421  """ In principle molck light but logs LigandScorer specific warnings
1422 
1423  Only to be applied to polymer entity
1424 
1425  1) removes atoms with elements set to H or D (not logged as there is no
1426  effect on scoring)
1427  2) removes residues with no entry in component dictionary
1428  3) removes all residues that are not peptide_liking or
1429  nucleotide_linking according component dictionary
1430  4) removes unknown atoms according to component dictionary
1431  5) reruns processor
1432  """
1433 
1434  cleanup_log = {"cleaned_residues": {"no_clib": [],
1435  "not_linking": []},
1436  "cleaned_atoms": {"unknown_atoms": []}}
1437 
1438  # 1)
1439  hydrogen_sel = ent.Select("ele == H or ele == D")
1440  if hydrogen_sel.GetAtomCount() > 0:
1441  # just do, no logging as it has no effect on scoring
1442  ent = ost.mol.CreateEntityFromView(ent.Select(
1443  "ele != H and ele != D"), include_exlusive_atoms=False)
1444 
1445  # extract residues/atoms for 2), 3) and 4)
1446  not_in_clib = list()
1447  not_linking = list()
1448  unknown_atom = list()
1449 
1450  for r in ent.residues:
1451  comp = clib.FindCompound(r.GetName())
1452  if comp is None:
1453  not_in_clib.append(r)
1454  continue
1455  if not (comp.IsPeptideLinking() or comp.IsNucleotideLinking()):
1456  not_linking.append(r)
1457  continue
1458  comp_anames = set([a.name for a in comp.atom_specs])
1459  for a in r.atoms:
1460  if a.name not in comp_anames:
1461  unknown_atom.append(a)
1462 
1463  # 2)
1464  if len(not_in_clib) > 0:
1465  cleanup_log["cleaned_residues"]["no_clib"] = \
1466  [_QualifiedResidueNotation(r) for r in not_in_clib]
1467  msg = ("Expect all residues in receptor structure to be defined in "
1468  "compound library but this is not the case. "
1469  "Please refer to the OpenStructure website if an updated "
1470  "library is sufficient to solve the problem: "
1471  "https://openstructure.org/docs/conop/compoundlib/ "
1472  "These residues will not be considered for scoring (which "
1473  "may also affect pre-processing steps such as alignments): ")
1474  msg += ','.join(cleanup_log["cleaned_residues"]["no_clib"])
1475  ost.LogWarning(msg)
1476  for r in not_in_clib:
1477  r.SetIntProp("clib", 1)
1478  ent = ost.mol.CreateEntityFromView(ent.Select("grclib:0!=1"),
1479  include_exlusive_atoms=False)
1480 
1481  # 3)
1482  if len(not_linking) > 0:
1483  cleanup_log["cleaned_residues"]["not_linking"] = \
1484  [_QualifiedResidueNotation(r) for r in not_linking]
1485  msg = ("Expect all residues in receptor structure to be peptide "
1486  "linking or nucleotide linking according to the compound "
1487  "library but this is not the case. "
1488  "Please refer to the OpenStructure website if an updated "
1489  "library is sufficient to solve the problem: "
1490  "https://openstructure.org/docs/conop/compoundlib/ "
1491  "These residues will not be considered for scoring (which "
1492  "may also affect pre-processing steps such as alignments): ")
1493  msg += ','.join(cleanup_log["cleaned_residues"]["not_linking"])
1494  ost.LogWarning(msg)
1495  for r in not_linking:
1496  r.SetIntProp("linking", 1)
1497  ent = ost.mol.CreateEntityFromView(ent.Select("grlinking:0!=1"),
1498  include_exlusive_atoms=False)
1499 
1500  # 4)
1501  if len(unknown_atom) > 0:
1502  cleanup_log["cleaned_atoms"]["unknown_atoms"] = \
1503  [_QualifiedAtomNotation(a) for a in unknown_atom]
1504  msg = ("Expect atom names according to the compound library but "
1505  "this is not the case."
1506  "Please refer to the OpenStructure website if an updated "
1507  "library is sufficient to solve the problem: "
1508  "https://openstructure.org/docs/conop/compoundlib/ "
1509  "These atoms will not be considered for scoring: ")
1510  msg += ','.join(cleanup_log["cleaned_atoms"]["unknown_atoms"])
1511  ost.LogWarning(msg)
1512  for a in unknown_atom:
1513  a.SetIntProp("unknown", 1)
1514  ent = ost.mol.CreateEntityFromView(ent.Select("gaunknown:0!=1"),
1515  include_exlusive_atoms=False)
1516 
1517  # 5)
1518  processor = conop.RuleBasedProcessor(clib)
1519  processor.Process(ent)
1520 
1521  return (ent, cleanup_log)
1522 
1523 
1524 def _ResidueToGraph(residue, by_atom_index=False):
1525  """Return a NetworkX graph representation of the residue.
1526 
1527  :param residue: the residue from which to derive the graph
1528  :type residue: :class:`ost.mol.ResidueHandle` or
1529  :class:`ost.mol.ResidueView`
1530  :param by_atom_index: Set this parameter to True if you need the nodes to
1531  be labeled by atom index (within the residue).
1532  Otherwise, if False, the nodes will be labeled by
1533  atom names.
1534  :type by_atom_index: :class:`bool`
1535  :rtype: :class:`~networkx.classes.graph.Graph`
1536 
1537  Nodes are labeled with the Atom's uppercase
1538  :attr:`~ost.mol.AtomHandle.element`.
1539  """
1540  nxg = networkx.Graph()
1541 
1542  for atom in residue.atoms:
1543  nxg.add_node(atom.name, element=atom.element.upper())
1544 
1545  # This will list all edges twice - once for every atom of the pair.
1546  # But as of NetworkX 3.0 adding the same edge twice has no effect,
1547  # so we're good.
1548  nxg.add_edges_from([(
1549  b.first.name,
1550  b.second.name) for a in residue.atoms for b in a.GetBondList()])
1551 
1552  if by_atom_index:
1553  nxg = networkx.relabel_nodes(nxg,
1554  {a: b for a, b in zip(
1555  [a.name for a in residue.atoms],
1556  range(len(residue.atoms)))},
1557  True)
1558  return nxg
1559 
1560 def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False,
1561  by_atom_index=False, return_symmetries=True,
1562  max_symmetries=1e6, model_graph = None,
1563  target_graph = None):
1564  """Return a list of symmetries (isomorphisms) of the model onto the target
1565  residues.
1566 
1567  :param model_ligand: The model ligand
1568  :type model_ligand: :class:`ost.mol.ResidueHandle` or
1569  :class:`ost.mol.ResidueView`
1570  :param target_ligand: The target ligand
1571  :type target_ligand: :class:`ost.mol.ResidueHandle` or
1572  :class:`ost.mol.ResidueView`
1573  :param substructure_match: Set this to True to allow partial ligands
1574  in the reference.
1575  :type substructure_match: :class:`bool`
1576  :param by_atom_index: Set this parameter to True if you need the symmetries
1577  to refer to atom index (within the residue).
1578  Otherwise, if False, the symmetries refer to atom
1579  names.
1580  :type by_atom_index: :class:`bool`
1581  :type return_symmetries: If Truthy, return the mappings, otherwise simply
1582  return True if a mapping is found (and raise if
1583  no mapping is found). This is useful to quickly
1584  find out if a mapping exist without the expensive
1585  step to find all the mappings.
1586  :type return_symmetries: :class:`bool`
1587  :param max_symmetries: If more than that many isomorphisms exist, raise
1588  a :class:`TooManySymmetriesError`. This can only be assessed by
1589  generating at least that many isomorphisms and can take some time.
1590  :type max_symmetries: :class:`int`
1591  :raises: :class:`NoSymmetryError` when no symmetry can be found;
1592  :class:`NoIsomorphicSymmetryError` in case of isomorphic
1593  subgraph but *substructure_match* is False;
1594  :class:`TooManySymmetriesError` when more than `max_symmetries`
1595  isomorphisms are found; :class:`DisconnectedGraphError` if
1596  graph for *model_ligand*/*target_ligand* is disconnected.
1597  """
1598 
1599  # Get the Graphs of the ligands
1600  if model_graph is None:
1601  model_graph = _ResidueToGraph(model_ligand,
1602  by_atom_index=by_atom_index)
1603 
1604  if not networkx.is_connected(model_graph):
1605  msg = "Disconnected graph for model ligand %s" % model_ligand
1606  raise DisconnectedGraphError(msg)
1607 
1608 
1609  if target_graph is None:
1610  target_graph = _ResidueToGraph(target_ligand,
1611  by_atom_index=by_atom_index)
1612 
1613  if not networkx.is_connected(target_graph):
1614  msg = "Disconnected graph for target ligand %s" % target_ligand
1615  raise DisconnectedGraphError(msg)
1616 
1617  # Note the argument order (model, target) which differs from spyrmsd.
1618  # This is because a subgraph of model is isomorphic to target - but not the
1619  # opposite as we only consider partial ligands in the reference.
1620  # Make sure to generate the symmetries correctly in the end
1621  gm = networkx.algorithms.isomorphism.GraphMatcher(
1622  model_graph, target_graph, node_match=lambda x, y:
1623  x["element"] == y["element"])
1624  if gm.is_isomorphic():
1625  if not return_symmetries:
1626  return True
1627  symmetries = []
1628  for i, isomorphism in enumerate(gm.isomorphisms_iter()):
1629  if i >= max_symmetries:
1630  raise TooManySymmetriesError(
1631  "Too many symmetries between %s and %s" % (
1632  str(model_ligand), str(target_ligand)))
1633  symmetries.append((list(isomorphism.values()),
1634  list(isomorphism.keys())))
1635  assert len(symmetries) > 0
1636  LogVerbose("Found %s isomorphic mappings (symmetries)" % len(symmetries))
1637  elif gm.subgraph_is_isomorphic() and substructure_match:
1638  if not return_symmetries:
1639  return True
1640  symmetries = []
1641  for i, isomorphism in enumerate(gm.subgraph_isomorphisms_iter()):
1642  if i >= max_symmetries:
1643  raise TooManySymmetriesError(
1644  "Too many symmetries between %s and %s" % (
1645  str(model_ligand), str(target_ligand)))
1646  symmetries.append((list(isomorphism.values()),
1647  list(isomorphism.keys())))
1648  assert len(symmetries) > 0
1649  # Assert that all the atoms in the target are part of the substructure
1650  assert len(symmetries[0][0]) == len(target_ligand.atoms)
1651  n_sym = len(symmetries)
1652  LogVerbose("Found %s subgraph isomorphisms (symmetries)" % n_sym)
1653  elif gm.subgraph_is_isomorphic():
1654  LogVerbose("Found subgraph isomorphisms (symmetries), but"
1655  " ignoring because substructure_match=False")
1656  raise NoIsomorphicSymmetryError("No symmetry between %s and %s" % (
1657  str(model_ligand), str(target_ligand)))
1658  else:
1659  LogVerbose("Found no isomorphic mappings (symmetries)")
1660  raise NoSymmetryError("No symmetry between %s and %s" % (
1661  str(model_ligand), str(target_ligand)))
1662 
1663  return symmetries
1664 
1665 class NoSymmetryError(ValueError):
1666  """ Exception raised when no symmetry can be found.
1667  """
1668  pass
1669 
1670 class NoIsomorphicSymmetryError(ValueError):
1671  """ Exception raised when no isomorphic symmetry can be found
1672 
1673  There would be isomorphic subgraphs for which symmetries can
1674  be found, but substructure_match is disabled
1675  """
1676  pass
1677 
1678 class TooManySymmetriesError(ValueError):
1679  """ Exception raised when too many symmetries are found.
1680  """
1681  pass
1682 
1683 class DisconnectedGraphError(Exception):
1684  """ Exception raised when the ligand graph is disconnected.
1685  """
1686  pass
1687 
1688 # specify public interface
1689 __all__ = ('CleanHydrogens', 'MMCIFPrep', 'PDBPrep',
1690  'LigandScorer', 'ComputeSymmetries', 'NoSymmetryError',
1691  'NoIsomorphicSymmetryError', 'TooManySymmetriesError',
1692  'DisconnectedGraphError')
Protein or molecule.
definition of EntityView
Definition: entity_view.hh:86
def _copy_ligand(self, l, ent, ed, rename_ligand_chain)
def __init__(self, model, target, model_ligands, target_ligands, resnum_alignments=False, substructure_match=False, coverage_delta=0.2, max_symmetries=1e5, rename_ligand_chain=False, min_pep_length=6, min_nuc_length=4, pep_seqid_thr=95., nuc_seqid_thr=95., mdl_map_pep_seqid_thr=0., mdl_map_nuc_seqid_thr=0.)
def _compute(self, symmetries, target_ligand, model_ligand)
def _get_report(self, ligand_state, pair_states)
def PopVerbosityLevel()
Definition: __init__.py:232
def PushVerbosityLevel(value)
Definition: __init__.py:229
def MMCIFPrep(mmcif_path, biounit=None, extract_nonpoly=False, fault_tolerant=False)
def PDBPrep(pdb_path, fault_tolerant=False)
def ComputeSymmetries(model_ligand, target_ligand, substructure_match=False, by_atom_index=False, return_symmetries=True, max_symmetries=1e6, model_graph=None, target_graph=None)
EntityHandle DLLEXPORT_OST_MOL CreateEntityFromView(const EntityView &view, bool include_exlusive_atoms, EntityHandle handle=EntityHandle())
create new entity handle from entity view