6 from ost
import settings
8 from ost
import LogScript, LogInfo, LogDebug
19 from ost
import bindings
25 """Scorer specific for a reference/model pair
27 Finds best possible binding site representation of reference in model given
28 lDDT score. Uses :class:`ost.mol.alg.chain_mapping.ChainMapper` to deal with
31 :param reference: Reference structure
32 :type reference: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
33 :param model: Model structure
34 :type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
35 :param residue_number_alignment: Passed to ChainMapper constructor
36 :type residue_number_alignment: :class:`bool`
39 residue_number_alignment=False):
41 resnum_alignments=residue_number_alignment)
45 def ScoreBS(self, ligand, radius = 4.0, lddt_radius=10.0):
46 """Computes binding site lDDT score given *ligand*. Best possible
47 binding site representation is selected by lDDT but other scores such as
48 CA based RMSD and GDT are computed too and returned.
50 :param ligand: Defines the scored binding site, i.e. provides positions
51 to perform proximity search
52 :type ligand: r'((Residue)|(Chain)|(Entity))((View)|(Handle))'
53 :param radius: Reference residues with any atom position within *radius*
54 of *ligand* consitute the scored binding site
55 :type radius: :class:`float`
56 :param lddt_radius: Passed as *inclusion_radius* to
57 :class:`ost.mol.alg.lddt.lDDTScorer`
58 :type lddt_radius: :class:`float`
59 :returns: Object of type :class:`ost.mol.alg.chain_mapping.ReprResult`
60 containing all atom lDDT score and mapping information.
61 None if no representation could be found.
65 ref_residues_hashes = set()
66 for ligand_at
in ligand.atoms:
67 close_atoms = self.
refref.FindWithin(ligand_at.GetPos(), radius)
68 for close_at
in close_atoms:
69 ref_res = close_at.GetResidue()
70 h = ref_res.handle.GetHashCode()
71 if h
not in ref_residues_hashes:
72 ref_residues_hashes.add(h)
77 ref_bs = self.
refref.CreateEmptyView()
78 for ch
in self.
refref.chains:
80 if r.handle.GetHashCode()
in ref_residues_hashes:
81 ref_bs.AddResidue(r, mol.ViewAddFlag.INCLUDE_ALL)
85 inclusion_radius = lddt_radius)
93 """ Helper class to access the various scores available from ost.mol.alg
95 Deals with structure cleanup, chain mapping, interface identification etc.
96 Intermediate results are available as attributes.
98 :param model: Model structure - a deep copy is available as :attr:`model`.
99 Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
101 :type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
102 :param target: Target structure - a deep copy is available as :attr:`target`.
103 Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
105 :type target: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
106 :param resnum_alignments: Whether alignments between chemically equivalent
107 chains in *model* and *target* can be computed
108 based on residue numbers. This can be assumed in
109 benchmarking setups such as CAMEO/CASP.
110 :type resnum_alignments: :class:`bool`
111 :param molck_settings: Settings used for Molck on *model* and *target*, if
112 set to None, a default object is constructed by
113 setting everything except rm_zero_occ_atoms and
115 :class:`ost.mol.alg.MolckSettings` constructor.
116 :type molck_settings: :class:`ost.mol.alg.MolckSettings`
117 :param cad_score_exec: Explicit path to voronota-cadscore executable from
118 voronota installation from
119 https://github.com/kliment-olechnovic/voronota. If
120 not given, voronota-cadscore must be in PATH if any
121 of the CAD score related attributes is requested.
122 :type cad_score_exec: :class:`str`
123 :param custom_mapping: Provide custom chain mapping between *model* and
124 *target*. Dictionary with target chain names as key
125 and model chain names as value.
126 :type custom_mapping: :class:`dict`
127 :param usalign_exec: Explicit path to USalign executable used to compute
128 TM-score. If not given, TM-score will be computed
129 with OpenStructure internal copy of USalign code.
130 :type usalign_exec: :class:`str`
131 :param lddt_no_stereochecks: Whether to compute lDDT without stereochemistry
133 :type lddt_no_stereochecks: :class:`bool`
134 :param n_max_naive: Parameter for chain mapping. If the number of possible
135 mappings is <= *n_max_naive*, the full
136 mapping solution space is enumerated to find the
137 the optimum. A heuristic is used otherwise. The default
138 of 40320 corresponds to an octamer (8! = 40320).
139 A structure with stoichiometry A6B2 would be
141 :type n_max_naive: :class:`int`
142 :param oum: Override USalign Mapping. Inject mapping of :class:`Scorer`
143 object into USalign to compute TM-score. Experimental feature
145 :type oum: :class:`bool`
146 :param min_pep_length: Relevant parameter if short peptides are involved in
147 scoring. Minimum peptide length for a chain in the
148 target structure to be considered in chain mapping.
149 The chain mapping algorithm first performs an all vs.
150 all pairwise sequence alignment to identify \"equal\"
151 chains within the target structure. We go for simple
152 sequence identity there. Short sequences can be
153 problematic as they may produce high sequence identity
154 alignments by pure chance.
155 :type min_pep_length: :class:`int`
156 :param min_nuc_length: Relevant parameter if short nucleotides are involved
157 in scoring. Minimum nucleotide length for a chain in
158 the target structure to be considered in chain
159 mapping. The chain mapping algorithm first performs
160 an all vs. all pairwise sequence alignment to
161 identify \"equal\" chains within the target
162 structure. We go for simple sequence identity there.
163 Short sequences can be problematic as they may
164 produce high sequence identity alignments by pure
166 :type min_nuc_length: :class:`int`
167 :param lddt_add_mdl_contacts: lDDT specific flag. Only using contacts in
168 lDDT that are within a certain distance
169 threshold in the target does not penalize
170 for added model contacts. If set to True, this
171 flag will also consider target contacts
172 that are within the specified distance
173 threshold in the model but not necessarily in
174 the target. No contact will be added if the
175 respective atom pair is not resolved in the
177 :type lddt_add_mdl_contacts: :class:`bool`
179 def __init__(self, model, target, resnum_alignments=False,
180 molck_settings = None, cad_score_exec = None,
181 custom_mapping=None, usalign_exec = None,
182 lddt_no_stereochecks=False, n_max_naive=40320,
183 oum=False, min_pep_length = 6, min_nuc_length = 4,
184 lddt_add_mdl_contacts=False):
199 if molck_settings
is None:
204 rm_zero_occ_atoms=
False,
208 LogScript(
"Cleaning up input structures")
209 Molck(self.
_model_model, conop.GetDefaultLib(), molck_settings)
210 Molck(self.
_target_target, conop.GetDefaultLib(), molck_settings)
211 self.
_model_model = self.
_model_model.Select(
"peptide=True or nucleotide=True")
212 self.
_target_target = self.
_target_target.Select(
"peptide=True or nucleotide=True")
215 for ch
in self.
_model_model.chains:
216 if ch.GetName().strip() ==
"":
217 raise RuntimeError(
"Model chains must have valid chain names")
218 if ch.GetName().strip() ==
"'" or ch.GetName().strip() ==
'"':
219 raise RuntimeError(
"Model chains cannot be named with quote signs (' or \"\")")
222 for ch
in self.
_target_target.chains:
223 if ch.GetName().strip() ==
"":
224 raise RuntimeError(
"Target chains must have valid chain names")
225 if ch.GetName().strip() ==
"'" or ch.GetName().strip() ==
'"':
226 raise RuntimeError(
"Target chains cannot be named with quote signs (' or \"\")")
228 if resnum_alignments:
232 for ch
in self.
_model_model.chains:
233 ins_codes = [r.GetNumber().GetInsCode()
for r
in ch.residues]
234 if len(set(ins_codes)) != 1
or ins_codes[0] !=
'\0':
235 raise RuntimeError(
"Residue numbers in each model chain "
236 "must not contain insertion codes if "
237 "resnum_alignments are enabled")
238 nums = [r.GetNumber().GetNum()
for r
in ch.residues]
239 if not all(i < j
for i, j
in zip(nums, nums[1:])):
240 raise RuntimeError(
"Residue numbers in each model chain "
241 "must be strictly increasing if "
242 "resnum_alignments are enabled")
244 for ch
in self.
_target_target.chains:
245 ins_codes = [r.GetNumber().GetInsCode()
for r
in ch.residues]
246 if len(set(ins_codes)) != 1
or ins_codes[0] !=
'\0':
247 raise RuntimeError(
"Residue numbers in each target chain "
248 "must not contain insertion codes if "
249 "resnum_alignments are enabled")
250 nums = [r.GetNumber().GetNum()
for r
in ch.residues]
251 if not all(i < j
for i, j
in zip(nums, nums[1:])):
252 raise RuntimeError(
"Residue numbers in each target chain "
253 "must be strictly increasing if "
254 "resnum_alignments are enabled")
256 if usalign_exec
is not None:
257 if not os.path.exists(usalign_exec):
258 raise RuntimeError(f
"USalign exec ({usalign_exec}) "
260 if not os.access(usalign_exec, os.X_OK):
261 raise RuntimeError(f
"USalign exec ({usalign_exec}) "
262 f
"is not executable")
299 self.
_lddt_lddt =
None
331 self.
_fnat_fnat =
None
335 self.
_nnat_nnat =
None
336 self.
_nmdl_nmdl =
None
362 self.
_rmsd_rmsd =
None
373 if custom_mapping
is not None:
375 LogDebug(
"Scorer sucessfully initialized")
379 """ Model with Molck cleanup
381 :type: :class:`ost.mol.EntityHandle`
387 """ The original model passed at object construction
389 :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
395 """ Target with Molck cleanup
397 :type: :class:`ost.mol.EntityHandle`
403 """ The original target passed at object construction
405 :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
411 """ Alignments of :attr:`model`/:attr:`target` chains
413 Alignments for each pair of chains mapped in :attr:`mapping`.
414 First sequence is target sequence, second sequence the model sequence.
416 :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
418 if self.
_aln_aln
is None:
424 """ Stereochecked equivalent of :attr:`aln`
426 The alignments may differ, as stereochecks potentially remove residues
428 :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
436 """ Alignments of :attr:`model_orig`/:attr:`target_orig` chains
438 Selects for peptide and nucleotide residues before sequence
439 extraction. Includes residues that would be removed by molck in
440 structure preprocessing.
442 :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
450 """ View of :attr:`~model` that has stereochemistry checks applied
452 First, a selection for peptide/nucleotide residues is performed,
453 secondly peptide sidechains with stereochemical irregularities are
454 removed (full residue if backbone atoms are involved). Irregularities
455 are clashes or bond lengths/angles more than 12 standard deviations
456 from expected values.
458 :type: :class:`ost.mol.EntityView`
466 """ Clashing model atoms
468 :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
476 """ Model bonds with unexpected stereochemistry
478 :type: :class:`list` of
479 :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
487 """ Model angles with unexpected stereochemistry
489 :type: :class:`list` of
490 :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
498 """ Same as :attr:`~stereochecked_model` for :attr:`~target`
500 :type: :class:`ost.mol.EntityView`
508 """ Clashing target atoms
510 :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
518 """ Target bonds with unexpected stereochemistry
520 :type: :class:`list` of
521 :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
529 """ Target angles with unexpected stereochemistry
531 :type: :class:`list` of
532 :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
540 """ Chain mapper object for given :attr:`target`
542 :type: :class:`ost.mol.alg.chain_mapping.ChainMapper`
554 """ Full chain mapping result for :attr:`target`/:attr:`model`
556 Computed with :func:`ost.mol.alg.ChainMapper.GetMapping`
558 :type: :class:`ost.mol.alg.chain_mapping.MappingResult`
561 LogScript(
"Computing chain mapping")
569 """ Full chain mapping result for :attr:`target`/:attr:`model`
571 Computed with :func:`ost.mol.alg.ChainMapper.GetRMSDMapping`
573 :type: :class:`ost.mol.alg.chain_mapping.MappingResult`
576 LogScript(
"Computing rigid chain mapping")
583 """ Interface residues in :attr:`~model`
585 Thats all residues having a contact with at least one residue from
586 another chain (CB-CB distance <= 8A, CA in case of Glycine)
588 :type: :class:`dict` with chain names as key and and :class:`list`
589 with residue numbers of the respective interface residues.
598 """ Same as :attr:`~model_interface_residues` for :attr:`~target`
600 :type: :class:`dict` with chain names as key and and :class:`list`
601 with residue numbers of the respective interface residues.
610 """ lDDT scorer for :attr:`~stereochecked_target` (default parameters)
612 :type: :class:`ost.mol.alg.lddt.lDDTScorer`
623 """ Backbone only lDDT scorer for :attr:`~target`
625 No stereochecks applied for bb only lDDT which considers CA atoms
626 for peptides and C3' atoms for nucleotides.
628 :type: :class:`ost.mol.alg.lddt.lDDTScorer`
636 """ QS scorer constructed from :attr:`~mapping`
638 The scorer object is constructed with default parameters and relates to
639 :attr:`~model` and :attr:`~target` (no stereochecks).
641 :type: :class:`ost.mol.alg.qsscore.QSScorer`
656 """ Global lDDT score in range [0.0, 1.0]
658 Computed based on :attr:`~stereochecked_model`. In case of oligomers,
659 :attr:`~mapping` is used.
661 :type: :class:`float`
663 if self.
_lddt_lddt
is None:
665 return self.
_lddt_lddt
669 """ Per residue lDDT scores in range [0.0, 1.0]
671 Computed based on :attr:`~stereochecked_model` but scores for all
672 residues in :attr:`~model` are reported. If a residue has been removed
673 by stereochemistry checks, the respective score is set to 0.0. If a
674 residue is not covered by the target or is in a chain skipped by the
675 chain mapping procedure (happens for super short chains), the respective
676 score is set to None. In case of oligomers, :attr:`~mapping` is used.
686 """ Backbone only global lDDT score in range [0.0, 1.0]
688 Computed based on :attr:`~model` on backbone atoms only. This is CA for
689 peptides and C3' for nucleotides. No stereochecks are performed. In case
690 of oligomers, :attr:`~mapping` is used.
692 :type: :class:`float`
700 """ Backbone only per residue lDDT scores in range [0.0, 1.0]
702 Computed based on :attr:`~model` on backbone atoms only. This is CA for
703 peptides and C3' for nucleotides. No stereochecks are performed. If a
704 residue is not covered by the target or is in a chain skipped by the
705 chain mapping procedure (happens for super short chains), the respective
706 score is set to None. In case of oligomers, :attr:`~mapping` is used.
718 Computed based on :attr:`model` using :attr:`mapping`
720 :type: :class:`float`
728 """ Global QS-score - only computed on aligned residues
730 Computed based on :attr:`model` using :attr:`mapping`. The QS-score
731 computation only considers contacts between residues with a mapping
732 between target and model. As a result, the score won't be lowered in
733 case of additional chains/residues in any of the structures.
735 :type: :class:`float`
743 """ Interfaces in :attr:`~target` with non-zero contribution to
744 :attr:`~qs_global`/:attr:`~qs_best`
746 Chain names are lexicographically sorted.
748 :type: :class:`list` of :class:`tuple` with 2 elements each:
759 """ Interfaces in :attr:`~model` with non-zero contribution to
760 :attr:`~qs_global`/:attr:`~qs_best`
762 Chain names are lexicographically sorted.
764 :type: :class:`list` of :class:`tuple` with 2 elements each:
776 """ Interfaces in :attr:`~qs_target_interfaces` that can be mapped
779 Target chain names are lexicographically sorted.
781 :type: :class:`list` of :class:`tuple` with 4 elements each:
782 (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
786 flat_mapping = self.
mappingmapping.GetFlatMapping()
788 if i[0]
in flat_mapping
and i[1]
in flat_mapping:
796 """ QS-score for each interface in :attr:`qs_interfaces`
798 :type: :class:`list` of :class:`float`
806 """ QS-score for each interface in :attr:`qs_interfaces`
808 Only computed on aligned residues
810 :type: :class:`list` of :class:`float`
820 A contact is a pair or residues from distinct chains that have
821 a minimal heavy atom distance < 5A. Contacts are specified as
822 :class:`tuple` with two strings in format:
823 <cname>.<rnum>.<ins_code>
825 :type: :class:`list` of :class:`tuple`
841 """ Interfaces in :class:`target` which have at least one contact
843 Contact as defined in :attr:`~native_contacts`,
844 chain names are lexicographically sorted.
846 :type: :class:`list` of :class:`tuple` with 2 elements each
851 tmp = [(min(x[0],x[1]), max(x[0],x[1]))
for x
in tmp]
857 """ Interfaces in :class:`model` which have at least one contact
859 Contact as defined in :attr:`~native_contacts`,
860 chain names are lexicographically sorted.
862 :type: :class:`list` of :class:`tuple` with 2 elements each
867 tmp = [(min(x[0],x[1]), max(x[0],x[1]))
for x
in tmp]
873 """ Fraction of model contacts that are also present in target
875 :type: :class:`float`
883 """ Fraction of target contacts that are correctly reproduced in model
885 :type: :class:`float`
893 """ ICS (Interface Contact Similarity) score
895 Combination of :attr:`~ics_precision` and :attr:`~ics_recall`
898 :type: :class:`float`
900 if self.
_ics_ics
is None:
906 """ Per-interface ICS precision
908 :attr:`~ics_precision` for each interface in
909 :attr:`~contact_target_interfaces`
911 :type: :class:`list` of :class:`float`
920 """ Per-interface ICS recall
922 :attr:`~ics_recall` for each interface in
923 :attr:`~contact_target_interfaces`
925 :type: :class:`list` of :class:`float`
933 """ Per-interface ICS (Interface Contact Similarity) score
935 :attr:`~ics` for each interface in
936 :attr:`~contact_target_interfaces`
938 :type: :class:`float`
948 """ Fraction of model interface residues that are also interface
951 :type: :class:`float`
959 """ Fraction of target interface residues that are also interface
962 :type: :class:`float`
970 """ IPS (Interface Patch Similarity) score
972 Jaccard coefficient of interface residues in target and their mapped
973 counterparts in model
975 :type: :class:`float`
977 if self.
_ips_ips
is None:
983 """ Per-interface IPS precision
985 :attr:`~ips_precision` for each interface in
986 :attr:`~contact_target_interfaces`
988 :type: :class:`list` of :class:`float`
997 """ Per-interface IPS recall
999 :attr:`~ips_recall` for each interface in
1000 :attr:`~contact_target_interfaces`
1002 :type: :class:`list` of :class:`float`
1010 """ Per-interface IPS (Interface Patch Similarity) score
1012 :attr:`~ips` for each interface in
1013 :attr:`~contact_target_interfaces`
1015 :type: :class:`list` of :class:`float`
1024 """ Interfaces in :attr:`target` that are relevant for DockQ
1026 In principle a subset of :attr:`~contact_target_interfaces` that only
1027 contains peptide sequences. Chain names are lexicographically sorted.
1029 :type: :class:`list` of :class:`tuple` with 2 elements each:
1034 pep_seqs = set([s.GetName()
for s
in self.
chain_mapperchain_mapper.polypep_seqs])
1036 if interface[0]
in pep_seqs
and interface[1]
in pep_seqs:
1042 """ Interfaces in :attr:`dockq_target_interfaces` that can be mapped
1045 Target chain names are lexicographically sorted
1047 :type: :class:`list` of :class:`tuple` with 4 elements each:
1048 (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
1052 flat_mapping = self.
mappingmapping.GetFlatMapping()
1054 if i[0]
in flat_mapping
and i[1]
in flat_mapping:
1057 flat_mapping[i[1]]))
1062 """ DockQ scores for interfaces in :attr:`~dockq_interfaces`
1064 :class:`list` of :class:`float`
1072 """ fnat scores for interfaces in :attr:`~dockq_interfaces`
1074 fnat: Fraction of native contacts that are also present in model
1076 :class:`list` of :class:`float`
1078 if self.
_fnat_fnat
is None:
1080 return self.
_fnat_fnat
1084 """ N native contacts for interfaces in :attr:`~dockq_interfaces`
1086 :class:`list` of :class:`int`
1088 if self.
_nnat_nnat
is None:
1090 return self.
_nnat_nnat
1094 """ N model contacts for interfaces in :attr:`~dockq_interfaces`
1096 :class:`list` of :class:`int`
1098 if self.
_nmdl_nmdl
is None:
1100 return self.
_nmdl_nmdl
1104 """ fnonnat scores for interfaces in :attr:`~dockq_interfaces`
1106 fnat: Fraction of model contacts that are not present in target
1108 :class:`list` of :class:`float`
1116 """ irmsd scores for interfaces in :attr:`~dockq_interfaces`
1118 irmsd: RMSD of interface (RMSD computed on N, CA, C, O atoms) which
1119 consists of each residue that has at least one heavy atom within 10A of
1122 :class:`list` of :class:`float`
1124 if self.
_irmsd_irmsd
is None:
1130 """ lrmsd scores for interfaces in :attr:`~dockq_interfaces`
1132 lrmsd: The interfaces are superposed based on the receptor (rigid
1133 min RMSD superposition) and RMSD for the ligand is reported.
1134 Superposition and RMSD are based on N, CA, C and O positions,
1135 receptor is the chain contributing to the interface with more
1138 :class:`list` of :class:`float`
1140 if self.
_lrmsd_lrmsd
is None:
1146 """ Average of DockQ scores in :attr:`dockq_scores`
1148 In its original implementation, DockQ only operates on single
1149 interfaces. Thus the requirement to combine scores for higher order
1152 :type: :class:`float`
1160 """ Same as :attr:`dockq_ave`, weighted by native contacts
1162 :type: :class:`float`
1170 """ Same as :attr:`~dockq_ave` but penalizing for missing interfaces
1172 Interfaces that are not covered in model are added as 0.0
1173 in average computation.
1175 :type: :class:`float`
1183 """ Same as :attr:`~dockq_ave_full`, but weighted
1185 Interfaces that are not covered in model are added as 0.0 in
1186 average computations and the respective weights are derived from
1187 number of contacts in respective target interface.
1195 """ Mapped representative positions in target
1197 Thats CA positions for peptide residues and C3' positions for
1198 nucleotides. Has same length as :attr:`~mapped_model_pos` and mapping
1199 is based on :attr:`~mapping`.
1201 :type: :class:`ost.geom.Vec3List`
1209 """ Mapped representative positions in model
1211 Thats CA positions for peptide residues and C3' positions for
1212 nucleotides. Has same length as :attr:`~mapped_target_pos` and mapping
1213 is based on :attr:`~mapping`.
1215 :type: :class:`ost.geom.Vec3List`
1223 """ :attr:`~mapped_model_pos` with :attr:`~transform` applied
1225 :type: :class:`ost.geom.Vec3List`
1235 """ Number of target residues which have no mapping to model
1245 """ Transform: :attr:`~mapped_model_pos` onto :attr:`~mapped_target_pos`
1247 Computed using Kabsch minimal rmsd algorithm
1249 :type: :class:`ost.geom.Mat4`
1255 self.
_transform_transform = res.transformation
1262 """ Mapped representative positions in target
1264 Thats CA positions for peptide residues and C3' positions for
1265 nucleotides. Has same length as :attr:`~rigid_mapped_model_pos` and mapping
1266 is based on :attr:`~rigid_mapping`.
1268 :type: :class:`ost.geom.Vec3List`
1276 """ Mapped representative positions in model
1278 Thats CA positions for peptide residues and C3' positions for
1279 nucleotides. Has same length as :attr:`~mapped_target_pos` and mapping
1280 is based on :attr:`~rigid_mapping`.
1282 :type: :class:`ost.geom.Vec3List`
1290 """ :attr:`~rigid_mapped_model_pos` with :attr:`~rigid_transform` applied
1292 :type: :class:`ost.geom.Vec3List`
1302 """ Number of target residues which have no rigid mapping to model
1312 """ Transform: :attr:`~rigid_mapped_model_pos` onto :attr:`~rigid_mapped_target_pos`
1314 Computed using Kabsch minimal rmsd algorithm
1316 :type: :class:`ost.geom.Mat4`
1329 """ Fraction CA (C3' for nucleotides) that can be superposed within 0.5A
1331 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1332 Similar iterative algorithm as LGA tool
1334 :type: :class:`float`
1336 if self.
_gdt_05_gdt_05
is None:
1340 self.
_gdt_05_gdt_05 = float(n) / n_full
1347 """ Fraction CA (C3' for nucleotides) that can be superposed within 1.0A
1349 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1350 Similar iterative algorithm as LGA tool
1352 :type: :class:`float`
1354 if self.
_gdt_1_gdt_1
is None:
1358 self.
_gdt_1_gdt_1 = float(n) / n_full
1365 """ Fraction CA (C3' for nucleotides) that can be superposed within 2.0A
1367 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1368 Similar iterative algorithm as LGA tool
1371 :type: :class:`float`
1373 if self.
_gdt_2_gdt_2
is None:
1377 self.
_gdt_2_gdt_2 = float(n) / n_full
1384 """ Fraction CA (C3' for nucleotides) that can be superposed within 4.0A
1386 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1387 Similar iterative algorithm as LGA tool
1389 :type: :class:`float`
1391 if self.
_gdt_4_gdt_4
is None:
1395 self.
_gdt_4_gdt_4 = float(n) / n_full
1402 """ Fraction CA (C3' for nucleotides) that can be superposed within 8.0A
1404 Similar iterative algorithm as LGA tool
1406 :type: :class:`float`
1408 if self.
_gdt_8_gdt_8
is None:
1412 self.
_gdt_8_gdt_8 = float(n) / n_full
1420 """ avg GDT with thresholds: 8.0A, 4.0A, 2.0A and 1.0A
1422 :type: :class:`float`
1424 if self.
_gdtts_gdtts
is None:
1425 LogScript(
"Computing GDT-TS score")
1431 """ avg GDT with thresholds: 4.0A, 2.0A, 1.0A and 0.5A
1433 :type: :class:`float`
1435 if self.
_gdtha_gdtha
is None:
1436 LogScript(
"Computing GDT-HA score")
1444 Computed on :attr:`~rigid_transformed_mapped_model_pos` and
1445 :attr:`rigid_mapped_target_pos`
1447 :type: :class:`float`
1449 if self.
_rmsd_rmsd
is None:
1450 LogScript(
"Computing RMSD")
1453 return self.
_rmsd_rmsd
1457 """ The global CAD atom-atom (AA) score
1459 Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
1462 :type: :class:`float`
1470 """ The per-residue CAD atom-atom (AA) scores
1472 Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
1475 :type: :class:`dict`
1483 """ Patch QS-scores for each residue in :attr:`model_interface_residues`
1485 Representative patches for each residue r in chain c are computed as
1488 * mdl_patch_one: All residues in c with CB (CA for GLY) positions within
1489 8A of r and within 12A of residues from any other chain.
1490 * mdl_patch_two: Closest residue x to r in any other chain gets
1491 identified. Patch is then constructed by selecting all residues from
1492 any other chain within 8A of x and within 12A from any residue in c.
1493 * trg_patch_one: Chain name and residue number based mapping from
1495 * trg_patch_two: Chain name and residue number based mapping from
1498 Results are stored in the same manner as
1499 :attr:`model_interface_residues`, with corresponding scores instead of
1500 residue numbers. Scores for residues which are not
1501 :class:`mol.ChemType.AMINOACIDS` are set to None. Additionally,
1502 interface patches are derived from :attr:`model`. If they contain
1503 residues which are not covered by :attr:`target`, the score is set to
1506 :type: :class:`dict` with chain names as key and and :class:`list`
1507 with scores of the respective interface residues.
1515 """ Same as :attr:`patch_qs` but for DockQ scores
1523 """ TM-score computed with USalign
1525 USalign executable can be specified with usalign_exec kwarg at Scorer
1526 construction, an OpenStructure internal copy of the USalign code is
1529 :type: :class:`float`
1537 """ Mapping computed with USalign
1539 Dictionary with target chain names as key and model chain names as
1540 values. No guarantee that all chains are mapped. USalign executable
1541 can be specified with usalign_exec kwarg at Scorer construction, an
1542 OpenStructure internal copy of the USalign code is used otherwise.
1544 :type: :class:`dict`
1550 def _aln_helper(self, target, model):
1555 for ch
in target.chains:
1556 cname = ch.GetName()
1557 s =
''.join([r.one_letter_code
for r
in ch.residues])
1558 s = seq.CreateSequence(ch.GetName(), s)
1559 s.AttachView(target.Select(f
"cname={mol.QueryQuoteName(cname)}"))
1560 trg_seqs[ch.GetName()] = s
1562 for ch
in model.chains:
1563 cname = ch.GetName()
1564 s =
''.join([r.one_letter_code
for r
in ch.residues])
1565 s = seq.CreateSequence(cname, s)
1566 s.AttachView(model.Select(f
"cname={mol.QueryQuoteName(cname)}"))
1567 mdl_seqs[ch.GetName()] = s
1570 trg_pep_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polypep_seqs]
1571 trg_nuc_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polynuc_seqs]
1572 trg_pep_chains = set(trg_pep_chains)
1573 trg_nuc_chains = set(trg_nuc_chains)
1574 for trg_ch, mdl_ch
in self.
mappingmapping.GetFlatMapping().items():
1575 if mdl_ch
in mdl_seqs
and trg_ch
in trg_seqs:
1576 if trg_ch
in trg_pep_chains:
1577 stype = mol.ChemType.AMINOACIDS
1578 elif trg_ch
in trg_nuc_chains:
1579 stype = mol.ChemType.NUCLEOTIDES
1581 raise RuntimeError(
"Chain name inconsistency... ask "
1583 alns.append(self.
chain_mapperchain_mapper.Align(trg_seqs[trg_ch],
1586 alns[-1].AttachView(0, trg_seqs[trg_ch].GetAttachedView())
1587 alns[-1].AttachView(1, mdl_seqs[mdl_ch].GetAttachedView())
1590 def _compute_pepnuc_aln(self):
1591 query =
"peptide=true or nucleotide=true"
1592 pep_nuc_target = self.
target_origtarget_orig.Select(query)
1593 pep_nuc_model = self.
model_origmodel_orig.Select(query)
1596 def _compute_aln(self):
1599 def _compute_stereochecked_aln(self):
1603 def _compute_lddt(self):
1604 LogScript(
"Computing all-atom lDDT")
1606 flat_mapping = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
1609 stereochecked_alns = dict()
1611 mdl_seq = aln.GetSequence(1)
1612 stereochecked_alns[mdl_seq.name] = aln
1614 for aln
in self.
alnaln:
1615 mdl_seq = aln.GetSequence(1)
1616 alns[mdl_seq.name] = aln
1623 lddt_chain_mapping = dict()
1624 for mdl_ch, trg_ch
in flat_mapping.items():
1626 lddt_chain_mapping[mdl_ch] = trg_ch
1628 chain_mapping = lddt_chain_mapping,
1629 residue_mapping = alns,
1630 check_resnames=
False,
1631 local_lddt_prop=
"lddt",
1634 for r
in self.
modelmodel.residues:
1635 cname = r.GetChain().GetName()
1636 if cname
not in local_lddt:
1637 local_lddt[cname] = dict()
1638 if r.HasProp(
"lddt"):
1639 score = round(r.GetFloatProp(
"lddt"), 3)
1640 local_lddt[cname][r.GetNumber()] = score
1644 local_lddt[cname][r.GetNumber()] =
None
1647 lddt_chain_mapping = dict()
1648 for mdl_ch, trg_ch
in flat_mapping.items():
1649 if mdl_ch
in stereochecked_alns:
1650 lddt_chain_mapping[mdl_ch] = trg_ch
1652 chain_mapping = lddt_chain_mapping,
1653 residue_mapping = stereochecked_alns,
1654 check_resnames=
False,
1655 local_lddt_prop=
"lddt",
1658 for r
in self.
modelmodel.residues:
1659 cname = r.GetChain().GetName()
1660 if cname
not in local_lddt:
1661 local_lddt[cname] = dict()
1662 if r.HasProp(
"lddt"):
1663 score = round(r.GetFloatProp(
"lddt"), 3)
1664 local_lddt[cname][r.GetNumber()] = score
1668 if mdl_res.IsValid():
1671 local_lddt[cname][r.GetNumber()] =
None
1678 if cname
in flat_mapping:
1679 for col
in alns[cname]:
1680 if col[0] !=
'-' and col[1] !=
'-':
1681 if col.GetResidue(1).number == r.number:
1682 trg_r = col.GetResidue(0)
1685 local_lddt[cname][r.GetNumber()] =
None
1687 local_lddt[cname][r.GetNumber()] = 0.0
1689 self.
_lddt_lddt = lddt_score
1692 def _compute_bb_lddt(self):
1693 LogScript(
"Computing backbone lDDT")
1696 for aln
in self.
alnaln:
1697 mdl_seq = aln.GetSequence(1)
1698 alns[mdl_seq.name] = aln
1701 flat_mapping = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
1702 lddt_chain_mapping = dict()
1703 for mdl_ch, trg_ch
in flat_mapping.items():
1705 lddt_chain_mapping[mdl_ch] = trg_ch
1708 chain_mapping = lddt_chain_mapping,
1709 residue_mapping = alns,
1710 check_resnames=
False,
1711 local_lddt_prop=
"bb_lddt",
1714 for r
in self.
modelmodel.residues:
1715 cname = r.GetChain().GetName()
1716 if cname
not in local_lddt:
1717 local_lddt[cname] = dict()
1718 if r.HasProp(
"bb_lddt"):
1719 score = round(r.GetFloatProp(
"bb_lddt"), 3)
1720 local_lddt[cname][r.GetNumber()] = score
1724 local_lddt[cname][r.GetNumber()] =
None
1729 def _compute_qs(self):
1730 LogScript(
"Computing global QS-score")
1731 qs_score_result = self.
qs_scorerqs_scorer.Score(self.
mappingmapping.mapping)
1732 self.
_qs_global_qs_global = qs_score_result.QS_global
1733 self.
_qs_best_qs_best = qs_score_result.QS_best
1735 def _compute_per_interface_qs_scores(self):
1736 LogScript(
"Computing per-interface QS-score")
1741 trg_ch1 = interface[0]
1742 trg_ch2 = interface[1]
1743 mdl_ch1 = interface[2]
1744 mdl_ch2 = interface[3]
1745 qs_res = self.
qs_scorerqs_scorer.ScoreInterface(trg_ch1, trg_ch2,
1750 def _compute_ics_scores(self):
1751 LogScript(
"Computing ICS scores")
1754 self.
_ics_recall_ics_recall = contact_scorer_res.recall
1755 self.
_ics_ics = contact_scorer_res.ics
1759 flat_mapping = self.
mappingmapping.GetFlatMapping()
1761 trg_ch1 = trg_int[0]
1762 trg_ch2 = trg_int[1]
1763 if trg_ch1
in flat_mapping
and trg_ch2
in flat_mapping:
1764 mdl_ch1 = flat_mapping[trg_ch1]
1765 mdl_ch2 = flat_mapping[trg_ch2]
1766 res = self.
contact_scorercontact_scorer.ScoreICSInterface(trg_ch1, trg_ch2,
1776 def _compute_ips_scores(self):
1777 LogScript(
"Computing IPS scores")
1780 self.
_ips_recall_ips_recall = contact_scorer_res.recall
1781 self.
_ips_ips = contact_scorer_res.ips
1786 flat_mapping = self.
mappingmapping.GetFlatMapping()
1788 trg_ch1 = trg_int[0]
1789 trg_ch2 = trg_int[1]
1790 if trg_ch1
in flat_mapping
and trg_ch2
in flat_mapping:
1791 mdl_ch1 = flat_mapping[trg_ch1]
1792 mdl_ch2 = flat_mapping[trg_ch2]
1793 res = self.
contact_scorercontact_scorer.ScoreIPSInterface(trg_ch1, trg_ch2,
1803 def _compute_dockq_scores(self):
1804 LogScript(
"Computing DockQ")
1807 self.
_fnat_fnat = list()
1809 self.
_irmsd_irmsd = list()
1810 self.
_lrmsd_lrmsd = list()
1811 self.
_nnat_nnat = list()
1812 self.
_nmdl_nmdl = list()
1815 for aln
in self.
alnaln:
1816 trg_s = aln.GetSequence(0)
1817 mdl_s = aln.GetSequence(1)
1818 dockq_alns[(trg_s.GetName(), mdl_s.GetName())] = aln
1821 trg_ch1 = interface[0]
1822 trg_ch2 = interface[1]
1823 mdl_ch1 = interface[2]
1824 mdl_ch2 = interface[3]
1825 aln1 = dockq_alns[(trg_ch1, mdl_ch1)]
1826 aln2 = dockq_alns[(trg_ch2, mdl_ch2)]
1827 res = dockq.DockQ(self.
modelmodel, self.
targettarget, mdl_ch1, mdl_ch2,
1828 trg_ch1, trg_ch2, ch1_aln=aln1,
1830 self.
_fnat_fnat.append(res[
"fnat"])
1831 self.
_fnonnat_fnonnat.append(res[
"fnonnat"])
1832 self.
_irmsd_irmsd.append(res[
"irmsd"])
1833 self.
_lrmsd_lrmsd.append(res[
"lrmsd"])
1835 self.
_nnat_nnat.append(res[
"nnat"])
1836 self.
_nmdl_nmdl.append(res[
"nmdl"])
1841 not_covered_counts = list()
1842 proc_trg_interfaces = set([(x[0], x[1])
for x
in self.
dockq_interfacesdockq_interfaces])
1844 if interface
not in proc_trg_interfaces:
1848 trg_ch1 = interface[0]
1849 trg_ch2 = interface[1]
1850 res = dockq.DockQ(self.
targettarget, self.
targettarget,
1851 trg_ch1, trg_ch2, trg_ch1, trg_ch2)
1852 not_covered_counts.append(res[
"nnat"])
1859 weights = np.array(self.
_nnat_nnat)
1864 self.
_dockq_wave_dockq_wave = np.sum(np.multiply(weights/np.sum(weights), scores))
1865 scores = np.append(scores, [0.0]*len(not_covered_counts))
1866 weights = np.append(weights, not_covered_counts)
1871 self.
_dockq_wave_full_dockq_wave_full = np.sum(np.multiply(weights/np.sum(weights),
1874 def _extract_mapped_pos(self):
1878 processed_trg_chains = set()
1879 for trg_ch, mdl_ch
in self.
mappingmapping.GetFlatMapping().items():
1880 processed_trg_chains.add(trg_ch)
1881 aln = self.
mappingmapping.alns[(trg_ch, mdl_ch)]
1883 if col[0] !=
'-' and col[1] !=
'-':
1884 trg_res = col.GetResidue(0)
1885 mdl_res = col.GetResidue(1)
1886 trg_at = trg_res.FindAtom(
"CA")
1887 mdl_at = mdl_res.FindAtom(
"CA")
1888 if not trg_at.IsValid():
1889 trg_at = trg_res.FindAtom(
"C3'")
1890 if not mdl_at.IsValid():
1891 mdl_at = mdl_res.FindAtom(
"C3'")
1897 for ch
in self.
mappingmapping.target.chains:
1898 if ch.GetName()
not in processed_trg_chains:
1902 def _extract_rigid_mapped_pos(self):
1906 processed_trg_chains = set()
1907 for trg_ch, mdl_ch
in self.
rigid_mappingrigid_mapping.GetFlatMapping().items():
1908 processed_trg_chains.add(trg_ch)
1909 aln = self.
rigid_mappingrigid_mapping.alns[(trg_ch, mdl_ch)]
1911 if col[0] !=
'-' and col[1] !=
'-':
1912 trg_res = col.GetResidue(0)
1913 mdl_res = col.GetResidue(1)
1914 trg_at = trg_res.FindAtom(
"CA")
1915 mdl_at = mdl_res.FindAtom(
"CA")
1916 if not trg_at.IsValid():
1917 trg_at = trg_res.FindAtom(
"C3'")
1918 if not mdl_at.IsValid():
1919 mdl_at = mdl_res.FindAtom(
"C3'")
1926 if ch.GetName()
not in processed_trg_chains:
1930 def _compute_cad_score(self):
1932 raise RuntimeError(
"CAD score computations rely on residue numbers "
1933 "that are consistent between target and model "
1934 "chains, i.e. only work if resnum_alignments "
1935 "is True at Scorer construction.")
1937 LogScript(
"Computing CAD score")
1939 settings.Locate(
"voronota-cadscore",
1941 except Exception
as e:
1942 raise RuntimeError(
"voronota-cadscore must be in PATH for CAD "
1943 "score scoring")
from e
1944 cad_bin_dir = os.path.dirname(cad_score_exec)
1945 m = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
1946 cad_result = cadscore.CADScore(self.
modelmodel, self.
targettarget,
1950 cad_bin_path=cad_bin_dir,
1954 for r
in self.
modelmodel.residues:
1955 cname = r.GetChain().GetName()
1956 if cname
not in local_cad:
1957 local_cad[cname] = dict()
1958 if r.HasProp(
"localcad"):
1959 score = round(r.GetFloatProp(
"localcad"), 3)
1960 local_cad[cname][r.GetNumber()] = score
1962 local_cad[cname][r.GetNumber()] =
None
1964 self.
_cad_score_cad_score = cad_result.globalAA
1967 def _get_repr_view(self, ent):
1968 """ Returns view with representative peptide atoms => CB, CA for GLY
1970 Ensures that each residue has exactly one atom with assertions
1972 :param ent: Entity for which you want the representative view
1973 :param ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
1974 :returns: :class:`ost.mol.EntityView` derived from ent
1976 repr_ent = ent.Select(
"(aname=\"CB\" or (rname=\"GLY\" and aname=\"CA\"))")
1977 for r
in repr_ent.residues:
1978 assert(len(r.atoms) == 1)
1981 def _get_interface_residues(self, ent):
1982 """ Get interface residues
1984 Thats all residues having a contact with at least one residue from another
1985 chain (CB-CB distance <= 8A, CA in case of Glycine)
1987 :param ent: Model for which to extract interface residues
1988 :type ent: :class:`ost.mol.EntityView`
1989 :returns: :class:`dict` with chain names as key and and :class:`list`
1990 with residue numbers of the respective interface residues.
1994 result = {ch.GetName(): list()
for ch
in ent.chains}
1995 for ch
in ent.chains:
1996 cname = ch.GetName()
1997 sel = repr_ent.Select(f
"(cname={mol.QueryQuoteName(cname)} and 8 <> [cname!={mol.QueryQuoteName(cname)}])")
1998 result[cname] = [r.GetNumber()
for r
in sel.residues]
2001 def _do_stereochecks(self):
2002 """ Perform stereochemistry checks on model and target
2004 LogInfo(
"Performing stereochemistry checks on model and target")
2005 data = stereochemistry.GetDefaultStereoData()
2006 l_data = stereochemistry.GetDefaultStereoLinkData()
2008 a, b, c, d = stereochemistry.StereoCheck(self.
modelmodel, stereo_data = data,
2009 stereo_link_data = l_data)
2015 a, b, c, d = stereochemistry.StereoCheck(self.
targettarget, stereo_data = data,
2016 stereo_link_data = l_data)
2023 def _get_interface_patches(self, mdl_ch, mdl_rnum):
2024 """ Select interface patches representative for specified residue
2026 The patches for specified residue r in chain c are selected as follows:
2028 * mdl_patch_one: All residues in c with CB (CA for GLY) positions within 8A
2029 of r and within 12A of residues from any other chain.
2030 * mdl_patch_two: Closest residue x to r in any other chain gets identified.
2031 Patch is then constructed by selecting all residues from any other chain
2032 within 8A of x and within 12A from any residue in c.
2033 * trg_patch_one: Chain name and residue number based mapping from
2035 * trg_patch_two: Chain name and residue number based mapping from
2038 :param mdl_ch: Name of chain in *self.model* of residue of interest
2039 :type mdl_ch: :class:`str`
2040 :param mdl_rnum: Residue number of residue of interest
2041 :type mdl_rnum: :class:`ost.mol.ResNum`
2042 :returns: Tuple with 5 elements: 1) :class:`bool` flag whether all residues
2043 in *mdl* patches are covered in *trg* 2) mtl_patch_one
2044 3) mdl_patch_two 4) trg_patch_one 5) trg_patch_two
2050 r = self.
modelmodel.FindResidue(mdl_ch, mdl_rnum)
2052 raise RuntimeError(f
"Cannot find residue {mdl_rnum} in chain {mdl_ch}")
2053 if r.GetName() ==
"GLY":
2054 at = r.FindAtom(
"CA")
2056 at = r.FindAtom(
"CB")
2057 if not at.IsValid():
2058 raise RuntimeError(
"Cannot find interface views for res without CB/CA")
2065 q1 = f
"(cname={mol.QueryQuoteName(mdl_ch)} and 8 <> {{{r_pos[0]},{r_pos[1]},{r_pos[2]}}})"
2067 q2 = f
"(12 <> [cname!={mol.QueryQuoteName(mdl_ch)}])"
2068 mdl_patch_one = self.
modelmodel.CreateEmptyView()
2069 sel = repr_mdl.Select(
" and ".join([q1, q2]))
2070 for r
in sel.residues:
2071 mdl_r = self.
modelmodel.FindResidue(r.GetChain().GetName(), r.GetNumber())
2072 mdl_patch_one.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
2078 sel = repr_mdl.Select(f
"(cname!={mol.QueryQuoteName(mdl_ch)})")
2079 close_stuff = sel.FindWithin(r_pos, 8)
2082 for close_at
in close_stuff:
2085 min_pos = close_at.GetPos()
2089 q1 = f
"(cname!={mol.QueryQuoteName(mdl_ch)} and 8 <> {{{min_pos[0]},{min_pos[1]},{min_pos[2]}}})"
2091 q2 = f
"(12 <> [cname={mol.QueryQuoteName(mdl_ch)}])"
2092 mdl_patch_two = self.
modelmodel.CreateEmptyView()
2093 sel = repr_mdl.Select(
" and ".join([q1, q2]))
2094 for r
in sel.residues:
2095 mdl_r = self.
modelmodel.FindResidue(r.GetChain().GetName(), r.GetNumber())
2096 mdl_patch_two.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
2099 flat_mapping = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
2100 full_trg_coverage =
True
2101 trg_patch_one = self.
targettarget.CreateEmptyView()
2102 for r
in mdl_patch_one.residues:
2104 mdl_cname = r.GetChain().GetName()
2105 if mdl_cname
in flat_mapping:
2106 aln = self.
mappingmapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
2108 if col[0] !=
'-' and col[1] !=
'-':
2109 if col.GetResidue(1).GetNumber() == r.GetNumber():
2110 trg_r = col.GetResidue(0)
2112 if trg_r
is not None:
2113 trg_patch_one.AddResidue(trg_r.handle,
2114 mol.ViewAddFlag.INCLUDE_ALL)
2116 full_trg_coverage =
False
2118 trg_patch_two = self.
targettarget.CreateEmptyView()
2119 for r
in mdl_patch_two.residues:
2121 mdl_cname = r.GetChain().GetName()
2122 if mdl_cname
in flat_mapping:
2123 aln = self.
mappingmapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
2125 if col[0] !=
'-' and col[1] !=
'-':
2126 if col.GetResidue(1).GetNumber() == r.GetNumber():
2127 trg_r = col.GetResidue(0)
2129 if trg_r
is not None:
2130 trg_patch_two.AddResidue(trg_r.handle,
2131 mol.ViewAddFlag.INCLUDE_ALL)
2133 full_trg_coverage =
False
2135 return (full_trg_coverage, mdl_patch_one, mdl_patch_two,
2136 trg_patch_one, trg_patch_two)
2138 def _compute_patchqs_scores(self):
2139 LogScript(
"Computing patch QS-scores")
2145 r = self.
modelmodel.FindResidue(cname, rnum)
2146 if r.IsValid()
and r.GetChemType() == mol.ChemType.AMINOACIDS:
2147 full_trg_coverage, mdl_patch_one, mdl_patch_two, \
2148 trg_patch_one, trg_patch_two = \
2150 if full_trg_coverage:
2151 score = self.
_patchqs_patchqs(mdl_patch_one, mdl_patch_two,
2152 trg_patch_one, trg_patch_two)
2153 scores.append(score)
2156 def _compute_patchdockq_scores(self):
2157 LogScript(
"Computing patch DockQ scores")
2163 r = self.
modelmodel.FindResidue(cname, rnum)
2164 if r.IsValid()
and r.GetChemType() == mol.ChemType.AMINOACIDS:
2165 full_trg_coverage, mdl_patch_one, mdl_patch_two, \
2166 trg_patch_one, trg_patch_two = \
2168 if full_trg_coverage:
2169 score = self.
_patchdockq_patchdockq(mdl_patch_one, mdl_patch_two,
2170 trg_patch_one, trg_patch_two)
2171 scores.append(score)
2174 def _patchqs(self, mdl_patch_one, mdl_patch_two, trg_patch_one, trg_patch_two):
2175 """ Score interface residue patches with QS-score
2177 In detail: Construct two entities with two chains each. First chain
2178 consists of residues from <x>_patch_one and second chain consists of
2179 <x>_patch_two. The returned score is the QS-score between the two
2182 :param mdl_patch_one: Interface patch representing scored residue
2183 :type mdl_patch_one: :class:`ost.mol.EntityView`
2184 :param mdl_patch_two: Interface patch representing scored residue
2185 :type mdl_patch_two: :class:`ost.mol.EntityView`
2186 :param trg_patch_one: Interface patch representing scored residue
2187 :type trg_patch_one: :class:`ost.mol.EntityView`
2188 :param trg_patch_two: Interface patch representing scored residue
2189 :type trg_patch_two: :class:`ost.mol.EntityView`
2190 :returns: PatchQS score
2195 alnA = seq.CreateAlignment()
2196 s =
''.join([r.one_letter_code
for r
in mdl_patch_one.residues])
2197 alnA.AddSequence(seq.CreateSequence(
"A", s))
2198 s =
''.join([r.one_letter_code
for r
in trg_patch_one.residues])
2199 alnA.AddSequence(seq.CreateSequence(
"A", s))
2201 alnB = seq.CreateAlignment()
2202 s =
''.join([r.one_letter_code
for r
in mdl_patch_two.residues])
2203 alnB.AddSequence(seq.CreateSequence(
"B", s))
2204 s =
''.join([r.one_letter_code
for r
in trg_patch_two.residues])
2205 alnB.AddSequence(seq.CreateSequence(
"B", s))
2206 alns = {(
"A",
"A"): alnA, (
"B",
"B"): alnB}
2208 scorer =
QSScorer(qs_ent_mdl, [[
"A"], [
"B"]], qs_ent_trg, alns)
2209 score_result = scorer.Score([[
"A"], [
"B"]])
2211 return score_result.QS_global
2213 def _patchdockq(self, mdl_patch_one, mdl_patch_two, trg_patch_one,
2215 """ Score interface residue patches with DockQ
2217 In detail: Construct two entities with two chains each. First chain
2218 consists of residues from <x>_patch_one and second chain consists of
2219 <x>_patch_two. The returned score is the QS-score between the two
2222 :param mdl_patch_one: Interface patch representing scored residue
2223 :type mdl_patch_one: :class:`ost.mol.EntityView`
2224 :param mdl_patch_two: Interface patch representing scored residue
2225 :type mdl_patch_two: :class:`ost.mol.EntityView`
2226 :param trg_patch_one: Interface patch representing scored residue
2227 :type trg_patch_one: :class:`ost.mol.EntityView`
2228 :param trg_patch_two: Interface patch representing scored residue
2229 :type trg_patch_two: :class:`ost.mol.EntityView`
2230 :returns: DockQ score
2234 dockq_result = dockq.DockQ(t, m,
"A",
"B",
"A",
"B")
2235 if dockq_result[
"nnat"] > 0:
2236 return dockq_result[
"DockQ"]
2239 def _qs_ent_from_patches(self, patch_one, patch_two):
2240 """ Constructs Entity with two chains named "A" and "B""
2242 Blindly adds all residues from *patch_one* to chain A and residues from
2243 patch_two to chain B.
2245 ent = mol.CreateEntity()
2247 added_ch = ed.InsertChain(
"A")
2248 for r
in patch_one.residues:
2249 added_r = ed.AppendResidue(added_ch, r.GetName())
2250 added_r.SetChemClass(str(r.GetChemClass()))
2252 ed.InsertAtom(added_r, a.handle)
2253 added_ch = ed.InsertChain(
"B")
2254 for r
in patch_two.residues:
2255 added_r = ed.AppendResidue(added_ch, r.GetName())
2256 added_r.SetChemClass(str(r.GetChemClass()))
2258 ed.InsertAtom(added_r, a.handle)
2261 def _set_custom_mapping(self, mapping):
2262 """ sets self._mapping with a full blown MappingResult object
2264 :param mapping: mapping with trg chains as key and mdl ch as values
2265 :type mapping: :class:`dict`
2267 LogInfo(
"Setting custom chain mapping")
2270 chem_mapping, chem_group_alns, mdl = \
2271 chain_mapper.GetChemMapping(self.
modelmodel)
2276 if len(mapping) != len(set(mapping.keys())):
2277 raise RuntimeError(f
"Expect unique trg chain names in mapping. Got "
2278 f
"{mapping.keys()}")
2279 if len(mapping) != len(set(mapping.values())):
2280 raise RuntimeError(f
"Expect unique mdl chain names in mapping. Got "
2281 f
"{mapping.values()}")
2283 trg_chains = set([ch.GetName()
for ch
in chain_mapper.target.chains])
2284 mdl_chains = set([ch.GetName()
for ch
in mdl.chains])
2285 for k,v
in mapping.items():
2286 if k
not in trg_chains:
2287 raise RuntimeError(f
"Target chain \"{k}\" is not available "
2288 f
"in target processed for chain mapping "
2290 if v
not in mdl_chains:
2291 raise RuntimeError(f
"Model chain \"{v}\" is not available "
2292 f
"in model processed for chain mapping "
2295 for trg_ch, mdl_ch
in mapping.items():
2296 trg_group_idx =
None
2297 mdl_group_idx =
None
2298 for idx, group
in enumerate(chain_mapper.chem_groups):
2302 for idx, group
in enumerate(chem_mapping):
2306 if trg_group_idx
is None or mdl_group_idx
is None:
2307 raise RuntimeError(
"Could not establish a valid chem grouping "
2308 "of chain names provided in custom mapping.")
2310 if trg_group_idx != mdl_group_idx:
2311 raise RuntimeError(f
"Chem group mismatch in custom mapping: "
2312 f
"target chain \"{trg_ch}\" groups with the "
2313 f
"following chemically equivalent target "
2315 f
"{chain_mapper.chem_groups[trg_group_idx]} "
2316 f
"but model chain \"{mdl_ch}\" maps to the "
2317 f
"following target chains: "
2318 f
"{chain_mapper.chem_groups[mdl_group_idx]}")
2320 pairs = set([(trg_ch, mdl_ch)
for trg_ch, mdl_ch
in mapping.items()])
2322 chain_mapping._GetRefMdlAlns(chain_mapper.chem_groups,
2323 chain_mapper.chem_group_alignments,
2329 final_mapping = list()
2330 for ref_chains
in chain_mapper.chem_groups:
2331 mapped_mdl_chains = list()
2332 for ref_ch
in ref_chains:
2333 if ref_ch
in mapping:
2334 mapped_mdl_chains.append(mapping[ref_ch])
2336 mapped_mdl_chains.append(
None)
2337 final_mapping.append(mapped_mdl_chains)
2340 for ref_group, mdl_group
in zip(chain_mapper.chem_groups,
2342 for ref_ch, mdl_ch
in zip(ref_group, mdl_group):
2343 if ref_ch
is not None and mdl_ch
is not None:
2344 aln = ref_mdl_alns[(ref_ch, mdl_ch)]
2345 trg_view = chain_mapper.target.Select(f
"cname={mol.QueryQuoteName(ref_ch)}")
2346 mdl_view = mdl.Select(f
"cname={mol.QueryQuoteName(mdl_ch)}")
2347 aln.AttachView(0, trg_view)
2348 aln.AttachView(1, mdl_view)
2349 alns[(ref_ch, mdl_ch)] = aln
2352 chain_mapper.chem_groups,
2354 final_mapping, alns)
2356 def _compute_tmscore(self):
2359 LogScript(
"Computing patch TM-score with USalign exectuable")
2361 flat_mapping = self.
mappingmapping.GetFlatMapping()
2362 LogInfo(
"Overriding TM-score chain mapping")
2363 res = res = bindings.WrappedMMAlign(self.
modelmodel, self.
targettarget,
2364 mapping=flat_mapping)
2366 res = bindings.WrappedMMAlign(self.
modelmodel, self.
targettarget)
2368 LogScript(
"Computing patch TM-score with built-in USalign")
2370 LogInfo(
"Overriding TM-score chain mapping")
2371 flat_mapping = self.
mappingmapping.GetFlatMapping()
2372 res = tmtools.USAlign(self.
modelmodel, self.
targettarget,
2374 custom_chain_mapping = flat_mapping)
2376 res = tmtools.USAlign(self.
modelmodel, self.
targettarget,
2381 for a,b
in zip(res.ent1_mapped_chains, res.ent2_mapped_chains):
def target_interface_residues(self)
def _patchqs(self, mdl_patch_one, mdl_patch_two, trg_patch_one, trg_patch_two)
_per_interface_ips_recall
def per_interface_ics(self)
def target_bad_angles(self)
def _get_repr_view(self, ent)
def qs_target_interfaces(self)
def rigid_mapped_model_pos(self)
def _do_stereochecks(self)
def _compute_per_interface_qs_scores(self)
def model_bad_angles(self)
_per_interface_ics_recall
_per_interface_ips_precision
def model_interface_residues(self)
_contact_target_interfaces
def stereochecked_target(self)
def per_interface_ips_precision(self)
def dockq_wave_full(self)
def mapped_model_pos(self)
def _compute_dockq_scores(self)
def _aln_helper(self, target, model)
def model_bad_bonds(self)
def per_interface_ics_recall(self)
def _compute_bb_lddt(self)
def transformed_mapped_model_pos(self)
def usalign_mapping(self)
def native_contacts(self)
def local_cad_score(self)
def _compute_patchqs_scores(self)
def per_interface_ips(self)
def _compute_stereochecked_aln(self)
_rigid_transformed_mapped_model_pos
def __init__(self, model, target, resnum_alignments=False, molck_settings=None, cad_score_exec=None, custom_mapping=None, usalign_exec=None, lddt_no_stereochecks=False, n_max_naive=40320, oum=False, min_pep_length=6, min_nuc_length=4, lddt_add_mdl_contacts=False)
def rigid_transformed_mapped_model_pos(self)
def n_target_not_mapped(self)
def contact_model_interfaces(self)
_transformed_mapped_model_pos
def _patchdockq(self, mdl_patch_one, mdl_patch_two, trg_patch_one, trg_patch_two)
def dockq_interfaces(self)
def per_interface_ips_recall(self)
def per_interface_ics_precision(self)
def _compute_ics_scores(self)
_rigid_n_target_not_mapped
def _get_interface_residues(self, ent)
def rigid_transform(self)
def _compute_ips_scores(self)
def _extract_mapped_pos(self)
def rigid_mapped_target_pos(self)
def _compute_tmscore(self)
def dockq_target_interfaces(self)
def stereochecked_aln(self)
_per_interface_ics_precision
def _compute_patchdockq_scores(self)
def _get_interface_patches(self, mdl_ch, mdl_rnum)
def _extract_rigid_mapped_pos(self)
_contact_model_interfaces
_target_interface_residues
def contact_target_interfaces(self)
def per_interface_qs_global(self)
def _qs_ent_from_patches(self, patch_one, patch_two)
def stereochecked_model(self)
_model_interface_residues
def _compute_pepnuc_aln(self)
def per_interface_qs_best(self)
def rigid_n_target_not_mapped(self)
def _set_custom_mapping(self, mapping)
def mapped_target_pos(self)
def target_bad_bonds(self)
def qs_model_interfaces(self)
def _compute_cad_score(self)
def ScoreBS(self, ligand, radius=4.0, lddt_radius=10.0)
def __init__(self, reference, model, residue_number_alignment=False)
Real DLLEXPORT_OST_GEOM Distance(const Line2 &l, const Vec2 &v)
void Molck(ost::mol::EntityHandle &ent, ost::conop::CompoundLibPtr lib, const MolckSettings &settings, bool prune=true)
void GDT(const geom::Vec3List &mdl_pos, const geom::Vec3List &ref_pos, int window_size, int max_windows, Real distance_thresh, int &n_superposed, geom::Mat4 &transform)