6 from ost
import settings
8 from ost
import LogScript, LogInfo, LogDebug
20 from ost
import bindings
26 """Scorer specific for a reference/model pair
28 Finds best possible binding site representation of reference in model given
29 LDDT score. Uses :class:`ost.mol.alg.chain_mapping.ChainMapper` to deal with
32 :param reference: Reference structure
33 :type reference: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
34 :param model: Model structure
35 :type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
36 :param residue_number_alignment: Passed to ChainMapper constructor
37 :type residue_number_alignment: :class:`bool`
40 residue_number_alignment=False):
42 resnum_alignments=residue_number_alignment)
46 def ScoreBS(self, ligand, radius = 4.0, lddt_radius=10.0):
47 """Computes binding site LDDT score given *ligand*. Best possible
48 binding site representation is selected by LDDT but other scores such as
49 CA based RMSD and GDT are computed too and returned.
51 :param ligand: Defines the scored binding site, i.e. provides positions
52 to perform proximity search
53 :type ligand: r'((Residue)|(Chain)|(Entity))((View)|(Handle))'
54 :param radius: Reference residues with any atom position within *radius*
55 of *ligand* consitute the scored binding site
56 :type radius: :class:`float`
57 :param lddt_radius: Passed as *inclusion_radius* to
58 :class:`ost.mol.alg.lddt.lDDTScorer`
59 :type lddt_radius: :class:`float`
60 :returns: Object of type :class:`ost.mol.alg.chain_mapping.ReprResult`
61 containing all atom LDDT score and mapping information.
62 None if no representation could be found.
66 ref_residues_hashes = set()
67 for ligand_at
in ligand.atoms:
68 close_atoms = self.
refref.FindWithin(ligand_at.GetPos(), radius)
69 for close_at
in close_atoms:
70 ref_res = close_at.GetResidue()
71 h = ref_res.handle.GetHashCode()
72 if h
not in ref_residues_hashes:
73 ref_residues_hashes.add(h)
78 ref_bs = self.
refref.CreateEmptyView()
79 for ch
in self.
refref.chains:
81 if r.handle.GetHashCode()
in ref_residues_hashes:
82 ref_bs.AddResidue(r, mol.ViewAddFlag.INCLUDE_ALL)
86 inclusion_radius = lddt_radius)
94 """ Helper class to access the various scores available from ost.mol.alg
96 Deals with structure cleanup, chain mapping, interface identification etc.
97 Intermediate results are available as attributes.
99 :param model: Model structure - a deep copy is available as :attr:`~model`.
100 Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
102 :type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
103 :param target: Target structure - a deep copy is available as :attr:`~target`.
104 Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
106 :type target: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
107 :param resnum_alignments: Whether alignments between chemically equivalent
108 chains in *model* and *target* can be computed
109 based on residue numbers. This can be assumed in
110 benchmarking setups such as CAMEO/CASP.
111 :type resnum_alignments: :class:`bool`
112 :param molck_settings: Settings used for Molck on *model* and *target*, if
113 set to None, a default object is constructed by
114 setting everything except rm_zero_occ_atoms and
116 :class:`ost.mol.alg.MolckSettings` constructor.
117 :type molck_settings: :class:`ost.mol.alg.MolckSettings`
118 :param cad_score_exec: Explicit path to voronota-cadscore executable from
119 voronota installation from
120 https://github.com/kliment-olechnovic/voronota. If
121 not given, voronota-cadscore must be in PATH if any
122 of the CAD score related attributes is requested.
123 :type cad_score_exec: :class:`str`
124 :param custom_mapping: Provide custom chain mapping between *model* and
125 *target*. Dictionary with target chain names as key
126 and model chain names as value.
127 :attr:`~mapping` is constructed from this.
128 :type custom_mapping: :class:`dict`
129 :param custom_rigid_mapping: Provide custom chain mapping between *model*
130 and *target*. Dictionary with target chain
131 names as key and model chain names as value.
132 :attr:`~rigid_mapping` is constructed from
134 :type custom_rigid_mapping: :class:`dict`
135 :param usalign_exec: Explicit path to USalign executable used to compute
136 TM-score. If not given, TM-score will be computed
137 with OpenStructure internal copy of USalign code.
138 :type usalign_exec: :class:`str`
139 :param lddt_no_stereochecks: Whether to compute LDDT without stereochemistry
141 :type lddt_no_stereochecks: :class:`bool`
142 :param n_max_naive: Parameter for chain mapping. If the number of possible
143 mappings is <= *n_max_naive*, the full
144 mapping solution space is enumerated to find the
145 the optimum. A heuristic is used otherwise. The default
146 of 40320 corresponds to an octamer (8! = 40320).
147 A structure with stoichiometry A6B2 would be
149 :type n_max_naive: :class:`int`
150 :param oum: Override USalign Mapping. Inject rigid_mapping of
151 :class:`Scorer` object into USalign to compute TM-score.
152 Experimental feature with limitations.
153 :type oum: :class:`bool`
154 :param min_pep_length: Relevant parameter if short peptides are involved in
155 scoring. Minimum peptide length for a chain in the
156 target structure to be considered in chain mapping.
157 The chain mapping algorithm first performs an all vs.
158 all pairwise sequence alignment to identify \"equal\"
159 chains within the target structure. We go for simple
160 sequence identity there. Short sequences can be
161 problematic as they may produce high sequence identity
162 alignments by pure chance.
163 :type min_pep_length: :class:`int`
164 :param min_nuc_length: Relevant parameter if short nucleotides are involved
165 in scoring. Minimum nucleotide length for a chain in
166 the target structure to be considered in chain
167 mapping. The chain mapping algorithm first performs
168 an all vs. all pairwise sequence alignment to
169 identify \"equal\" chains within the target
170 structure. We go for simple sequence identity there.
171 Short sequences can be problematic as they may
172 produce high sequence identity alignments by pure
174 :type min_nuc_length: :class:`int`
175 :param lddt_add_mdl_contacts: LDDT specific flag. Only using contacts in
176 LDDT that are within a certain distance
177 threshold in the target does not penalize
178 for added model contacts. If set to True, this
179 flag will also consider target contacts
180 that are within the specified distance
181 threshold in the model but not necessarily in
182 the target. No contact will be added if the
183 respective atom pair is not resolved in the
185 :type lddt_add_mdl_contacts: :class:`bool`
186 :param dockq_capri_peptide: Flag that changes two things in the way DockQ
187 and its underlying scores are computed which is
188 proposed by the CAPRI community when scoring
189 peptides (PMID: 31886916).
190 ONE: Two residues are considered in contact if
191 any of their atoms is within 5A. This is
192 relevant for fnat and fnonat scores. CAPRI
193 suggests to lower this threshold to 4A for
194 protein-peptide interactions.
195 TWO: irmsd is computed on interface residues. A
196 residue is defined as interface residue if any
197 of its atoms is within 10A of another chain.
198 CAPRI suggests to lower the default of 10A to
199 8A in combination with only considering CB atoms
200 for protein-peptide interactions.
201 This flag has no influence on patch_dockq
203 :type dockq_capri_peptide: :class:`bool`
204 :param lddt_symmetry_settings: Passed as *symmetry_settings* parameter to
205 LDDT scorer. Default: None
206 :type lddt_symmetry_settings: :class:`ost.mol.alg.lddt.SymmetrySettings`
207 :param lddt_inclusion_radius: LDDT inclusion radius.
208 :param pep_seqid_thr: Parameter that affects identification of identical
209 chains in target - see
210 :class:`ost.mol.alg.chain_mapping.ChainMapper`
211 :type pep_seqid_thr: :class:`float`
212 :param nuc_seqid_thr: Parameter that affects identification of identical
213 chains in target - see
214 :class:`ost.mol.alg.chain_mapping.ChainMapper`
215 :type nuc_seqid_thr: :class:`float`
216 :param mdl_map_pep_seqid_thr: Parameter that affects mapping of model chains
217 to target chains - see
218 :class:`ost.mol.alg.chain_mapping.ChainMapper`
219 :type mdl_map_pep_seqid_thr: :class:`float`
220 :param mdl_map_nuc_seqid_thr: Parameter that affects mapping of model chains
221 to target chains - see
222 :class:`ost.mol.alg.chain_mapping.ChainMapper`
223 :type mdl_map_nuc_seqid_thr: :class:`float`
225 def __init__(self, model, target, resnum_alignments=False,
226 molck_settings = None, cad_score_exec = None,
227 custom_mapping=None, custom_rigid_mapping=None,
228 usalign_exec = None, lddt_no_stereochecks=False,
229 n_max_naive=40320, oum=False, min_pep_length = 6,
230 min_nuc_length = 4, lddt_add_mdl_contacts=False,
231 dockq_capri_peptide=False, lddt_symmetry_settings = None,
232 lddt_inclusion_radius = 15.0,
233 pep_seqid_thr = 95., nuc_seqid_thr = 95.,
234 mdl_map_pep_seqid_thr = 0.,
235 mdl_map_nuc_seqid_thr = 0.):
250 if molck_settings
is None:
255 rm_zero_occ_atoms=
False,
259 LogScript(
"Cleaning up input structures")
260 Molck(self.
_model_model, conop.GetDefaultLib(), molck_settings)
261 Molck(self.
_target_target, conop.GetDefaultLib(), molck_settings)
263 if resnum_alignments:
276 self.
_model_model = self.
_model_model.Select(
"peptide=True or nucleotide=True")
277 self.
_target_target = self.
_target_target.Select(
"peptide=True or nucleotide=True")
280 for ch
in self.
_model_model.chains:
281 if ch.GetName().strip() ==
"":
282 raise RuntimeError(
"Model chains must have valid chain names")
283 if ch.GetName().strip() ==
"'" or ch.GetName().strip() ==
'"':
284 raise RuntimeError(
"Model chains cannot be named with quote signs (' or \"\")")
287 for ch
in self.
_target_target.chains:
288 if ch.GetName().strip() ==
"":
289 raise RuntimeError(
"Target chains must have valid chain names")
290 if ch.GetName().strip() ==
"'" or ch.GetName().strip() ==
'"':
291 raise RuntimeError(
"Target chains cannot be named with quote signs (' or \"\")")
293 if resnum_alignments:
297 for ch
in self.
_model_model.chains:
298 ins_codes = [r.GetNumber().GetInsCode()
for r
in ch.residues]
299 if len(set(ins_codes)) != 1
or ins_codes[0] !=
'\0':
300 raise RuntimeError(
"Residue numbers in each model chain "
301 "must not contain insertion codes if "
302 "resnum_alignments are enabled")
303 nums = [r.GetNumber().GetNum()
for r
in ch.residues]
304 if not all(i < j
for i, j
in zip(nums, nums[1:])):
305 raise RuntimeError(
"Residue numbers in each model chain "
306 "must be strictly increasing if "
307 "resnum_alignments are enabled")
309 for ch
in self.
_target_target.chains:
310 ins_codes = [r.GetNumber().GetInsCode()
for r
in ch.residues]
311 if len(set(ins_codes)) != 1
or ins_codes[0] !=
'\0':
312 raise RuntimeError(
"Residue numbers in each target chain "
313 "must not contain insertion codes if "
314 "resnum_alignments are enabled")
315 nums = [r.GetNumber().GetNum()
for r
in ch.residues]
316 if not all(i < j
for i, j
in zip(nums, nums[1:])):
317 raise RuntimeError(
"Residue numbers in each target chain "
318 "must be strictly increasing if "
319 "resnum_alignments are enabled")
321 if usalign_exec
is not None:
322 if not os.path.exists(usalign_exec):
323 raise RuntimeError(f
"USalign exec ({usalign_exec}) "
325 if not os.access(usalign_exec, os.X_OK):
326 raise RuntimeError(f
"USalign exec ({usalign_exec}) "
327 f
"is not executable")
374 self.
_lddt_lddt =
None
425 self.
_fnat_fnat =
None
429 self.
_nnat_nnat =
None
430 self.
_nmdl_nmdl =
None
461 self.
_rmsd_rmsd =
None
472 if custom_mapping
is not None:
475 if custom_rigid_mapping
is not None:
478 LogDebug(
"Scorer sucessfully initialized")
482 """ Model with Molck cleanup
484 :type: :class:`ost.mol.EntityHandle`
490 """ The original model passed at object construction
492 :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
498 """ Target with Molck cleanup
500 :type: :class:`ost.mol.EntityHandle`
506 """ The original target passed at object construction
508 :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
514 """ Alignments of :attr:`~model`/:attr:`~target` chains
516 Alignments for each pair of chains mapped in :attr:`~mapping`.
517 First sequence is target sequence, second sequence the model sequence.
519 :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
521 if self.
_aln_aln
is None:
527 """ Stereochecked equivalent of :attr:`~aln`
529 The alignments may differ, as stereochecks potentially remove residues
531 :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
539 """ Alignments of :attr:`~model_orig`/:attr:`~target_orig` chains
541 Selects for peptide and nucleotide residues before sequence
542 extraction. Includes residues that would be removed by molck in
543 structure preprocessing.
545 :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
553 """ Alignments of :attr:`~trimmed_model`/:attr:`~target` chains
555 Alignments for each pair of chains mapped in :attr:`~mapping`.
556 First sequence is target sequence, second sequence the model sequence.
558 :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
566 """ View of :attr:`~model` that has stereochemistry checks applied
568 First, a selection for peptide/nucleotide residues is performed,
569 secondly peptide sidechains with stereochemical irregularities are
570 removed (full residue if backbone atoms are involved). Irregularities
571 are clashes or bond lengths/angles more than 12 standard deviations
572 from expected values.
574 :type: :class:`ost.mol.EntityView`
582 """ Clashing model atoms
584 :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
592 """ Model bonds with unexpected stereochemistry
594 :type: :class:`list` of
595 :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
603 """ Model angles with unexpected stereochemistry
605 :type: :class:`list` of
606 :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
614 """ Same as :attr:`~stereochecked_model` for :attr:`~target`
616 :type: :class:`ost.mol.EntityView`
624 """ Clashing target atoms
626 :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
634 """ Target bonds with unexpected stereochemistry
636 :type: :class:`list` of
637 :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
645 """ Target angles with unexpected stereochemistry
647 :type: :class:`list` of
648 :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
656 """ :attr:`~model` trimmed to target
658 Removes residues that are not covered by :class:`target` given
659 :attr:`~mapping`. In other words: no model residues without experimental
660 evidence from :class:`target`.
662 :type: :class:`ost.mol.EntityView`
670 """ Chain mapper object for given :attr:`~target`
672 :type: :class:`ost.mol.alg.chain_mapping.ChainMapper`
688 """ Full chain mapping result for :attr:`~target`/:attr:`~model`
690 Computed with :func:`ost.mol.alg.ChainMapper.GetMapping`
692 :type: :class:`ost.mol.alg.chain_mapping.MappingResult`
695 LogScript(
"Computing chain mapping")
703 """ Full chain mapping result for :attr:`~target`/:attr:`~model`
705 Computed with :func:`ost.mol.alg.ChainMapper.GetRMSDMapping`
707 :type: :class:`ost.mol.alg.chain_mapping.MappingResult`
710 LogScript(
"Computing rigid chain mapping")
717 """ Interface residues in :attr:`~model`
719 Thats all residues having a contact with at least one residue from
720 another chain (CB-CB distance <= 8A, CA in case of Glycine)
722 :type: :class:`dict` with chain names as key and and :class:`list`
723 with residue numbers of the respective interface residues.
732 """ Same as :attr:`~model_interface_residues` for :attr:`~target`
734 :type: :class:`dict` with chain names as key and and :class:`list`
735 with residue numbers of the respective interface residues.
744 """ LDDT scorer for :attr:`~target`/:attr:`~stereochecked_target`
746 Depending on :attr:`~lddt_no_stereocheck` and
747 :attr:`~lddt_symmetry_settings`.
749 :type: :class:`ost.mol.alg.lddt.lDDTScorer`
764 """ LDDT scorer for :attr:`~target`, restricted to representative
767 No stereochecks applied for bb only LDDT which considers CA atoms
768 for peptides and C3' atoms for nucleotides.
770 :type: :class:`ost.mol.alg.lddt.lDDTScorer`
780 """ QS scorer constructed from :attr:`~mapping`
782 The scorer object is constructed with default parameters and relates to
783 :attr:`~model` and :attr:`~target` (no stereochecks).
785 :type: :class:`ost.mol.alg.qsscore.QSScorer`
801 self.
mappingmapping.chem_groups,
808 """ Global LDDT score in range [0.0, 1.0]
810 Computed based on :attr:`~stereochecked_model`. In case of oligomers,
811 :attr:`~mapping` is used.
813 :type: :class:`float`
815 if self.
_lddt_lddt
is None:
817 return self.
_lddt_lddt
821 """ Per residue LDDT scores in range [0.0, 1.0]
823 Computed based on :attr:`~stereochecked_model` but scores for all
824 residues in :attr:`~model` are reported. If a residue has been removed
825 by stereochemistry checks, the respective score is set to 0.0. If a
826 residue is not covered by the target or is in a chain skipped by the
827 chain mapping procedure (happens for super short chains), the respective
828 score is set to None. In case of oligomers, :attr:`~mapping` is used.
838 """ Per atom LDDT scores in range [0.0, 1.0]
840 Computed based on :attr:`~stereochecked_model` but scores for all
841 atoms in :attr:`~model` are reported. If an atom has been removed
842 by stereochemistry checks, the respective score is set to 0.0. If an
843 atom is not covered by the target or is in a chain skipped by the
844 chain mapping procedure (happens for super short chains), the respective
845 score is set to None. In case of oligomers, :attr:`~mapping` is used.
855 """ Global LDDT score restricted to representative backbone atoms in
858 Computed based on :attr:`~model` on representative backbone atoms only.
859 This is CA for peptides and C3' for nucleotides. No stereochecks are
860 performed. In case of oligomers, :attr:`~mapping` is used.
862 :type: :class:`float`
870 """ Per residue LDDT scores restricted to representative backbone atoms
873 Computed based on :attr:`~model` on representative backbone atoms only.
874 This is CA for peptides and C3' for nucleotides. No stereochecks are
875 performed. If a residue is not covered by the target or is in a chain
876 skipped by the chain mapping procedure (happens for super short
877 chains), the respective score is set to None. In case of oligomers,
878 :attr:`~mapping` is used.
888 """ Global interface LDDT score in range [0.0, 1.0]
890 This is LDDT only based on inter-chain contacts. Value is None if no
891 such contacts are present. For example if we're dealing with a monomer.
892 Computed based on :attr:`~stereochecked_model` and :attr:`~mapping` for
895 :type: :class:`float`
897 if self.
_ilddt_ilddt
is None:
908 Computed based on :attr:`~model` using :attr:`~mapping`
910 :type: :class:`float`
918 """ Global QS-score - only computed on aligned residues
920 Computed based on :attr:`~model` using :attr:`~mapping`. The QS-score
921 computation only considers contacts between residues with a mapping
922 between target and model. As a result, the score won't be lowered in
923 case of additional chains/residues in any of the structures.
925 :type: :class:`float`
933 """ Interfaces in :attr:`~target` with non-zero contribution to
934 :attr:`~qs_global`/:attr:`~qs_best`
936 Chain names are lexicographically sorted.
938 :type: :class:`list` of :class:`tuple` with 2 elements each:
949 """ Interfaces in :attr:`~model` with non-zero contribution to
950 :attr:`~qs_global`/:attr:`~qs_best`
952 Chain names are lexicographically sorted.
954 :type: :class:`list` of :class:`tuple` with 2 elements each:
966 """ Interfaces in :attr:`~qs_target_interfaces` that can be mapped
969 Target chain names are lexicographically sorted.
971 :type: :class:`list` of :class:`tuple` with 4 elements each:
972 (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
976 flat_mapping = self.
mappingmapping.GetFlatMapping()
978 if i[0]
in flat_mapping
and i[1]
in flat_mapping:
986 """ QS-score for each interface in :attr:`~qs_interfaces`
988 :type: :class:`list` of :class:`float`
996 """ QS-score for each interface in :attr:`~qs_interfaces`
998 Only computed on aligned residues
1000 :type: :class:`list` of :class:`float`
1010 A contact is a pair or residues from distinct chains that have
1011 a minimal heavy atom distance < 5A. Contacts are specified as
1012 :class:`tuple` with two strings in format:
1013 <cname>.<rnum>.<ins_code>
1015 :type: :class:`list` of :class:`tuple`
1023 """ Same for :attr:`~model`
1031 """ Same for :attr:`~trimmed_model`
1039 """ Interfaces in :class:`target` which have at least one contact
1041 Contact as defined in :attr:`~native_contacts`,
1042 chain names are lexicographically sorted.
1044 :type: :class:`list` of :class:`tuple` with 2 elements each
1049 tmp = [(min(x[0],x[1]), max(x[0],x[1]))
for x
in tmp]
1055 """ Interfaces in :class:`model` which have at least one contact
1057 Contact as defined in :attr:`~native_contacts`,
1058 chain names are lexicographically sorted.
1060 :type: :class:`list` of :class:`tuple` with 2 elements each
1065 tmp = [(min(x[0],x[1]), max(x[0],x[1]))
for x
in tmp]
1071 """ Fraction of model contacts that are also present in target
1073 :type: :class:`float`
1081 """ Fraction of target contacts that are correctly reproduced in model
1083 :type: :class:`float`
1091 """ ICS (Interface Contact Similarity) score
1093 Combination of :attr:`~ics_precision` and :attr:`~ics_recall`
1094 using the F1-measure
1096 :type: :class:`float`
1098 if self.
_ics_ics
is None:
1100 return self.
_ics_ics
1104 """ Per-interface ICS precision
1106 :attr:`~ics_precision` for each interface in
1107 :attr:`~contact_target_interfaces`
1109 :type: :class:`list` of :class:`float`
1118 """ Per-interface ICS recall
1120 :attr:`~ics_recall` for each interface in
1121 :attr:`~contact_target_interfaces`
1123 :type: :class:`list` of :class:`float`
1131 """ Per-interface ICS (Interface Contact Similarity) score
1133 :attr:`~ics` for each interface in
1134 :attr:`~contact_target_interfaces`
1136 :type: :class:`float`
1145 """ Fraction of model interface residues that are also interface
1148 :type: :class:`float`
1156 """ Fraction of target interface residues that are also interface
1159 :type: :class:`float`
1167 """ IPS (Interface Patch Similarity) score
1169 Jaccard coefficient of interface residues in target and their mapped
1170 counterparts in model
1172 :type: :class:`float`
1174 if self.
_ips_ips
is None:
1176 return self.
_ips_ips
1180 """ Same as :attr:`~ics` but with trimmed model
1182 Model is trimmed to residues which can me mapped to target in order
1183 to not penalize contacts in the model for which we have no experimental
1186 :type: :class:`float`
1194 """ Same as :attr:`~ics_precision` but with trimmed model
1196 Model is trimmed to residues which can me mapped to target in order
1197 to not penalize contacts in the model for which we have no experimental
1200 :type: :class:`float`
1208 """ Same as :attr:`~ics_recall` but with trimmed model
1210 Model is trimmed to residues which can me mapped to target in order
1211 to not penalize contacts in the model for which we have no experimental
1214 :type: :class:`float`
1222 """ Same as :attr:`~per_interface_ics_precision` but with :attr:`~trimmed_model`
1224 :attr:`~ics_precision_trimmed` for each interface in
1225 :attr:`~contact_target_interfaces`
1227 :type: :class:`list` of :class:`float`
1236 """ Same as :attr:`~per_interface_ics_recall` but with :attr:`~trimmed_model`
1238 :attr:`~ics_recall_trimmed` for each interface in
1239 :attr:`~contact_target_interfaces`
1241 :type: :class:`list` of :class:`float`
1249 """ Same as :attr:`~per_interface_ics` but with :attr:`~trimmed_model`
1251 :attr:`~ics` for each interface in
1252 :attr:`~contact_target_interfaces`
1254 :type: :class:`float`
1263 """ Same as :attr:`~ips` but with trimmed model
1265 Model is trimmed to residues which can me mapped to target in order
1266 to not penalize contacts in the model for which we have no experimental
1269 :type: :class:`float`
1277 """ Same as :attr:`~ips_precision` but with trimmed model
1279 Model is trimmed to residues which can me mapped to target in order
1280 to not penalize contacts in the model for which we have no experimental
1283 :type: :class:`float`
1291 """ Same as :attr:`~ips_recall` but with trimmed model
1293 Model is trimmed to residues which can me mapped to target in order
1294 to not penalize contacts in the model for which we have no experimental
1297 :type: :class:`float`
1305 """ Per-interface IPS precision
1307 :attr:`~ips_precision` for each interface in
1308 :attr:`~contact_target_interfaces`
1310 :type: :class:`list` of :class:`float`
1318 """ Per-interface IPS recall
1320 :attr:`~ips_recall` for each interface in
1321 :attr:`~contact_target_interfaces`
1323 :type: :class:`list` of :class:`float`
1331 """ Per-interface IPS (Interface Patch Similarity) score
1333 :attr:`~ips` for each interface in
1334 :attr:`~contact_target_interfaces`
1336 :type: :class:`list` of :class:`float`
1345 """ Same as :attr:`~per_interface_ips_precision` but with :attr:`~trimmed_model`
1347 :attr:`~ips_precision_trimmed` for each interface in
1348 :attr:`~contact_target_interfaces`
1350 :type: :class:`list` of :class:`float`
1359 """ Same as :attr:`~per_interface_ips_recall` but with :attr:`~trimmed_model`
1361 :attr:`~ics_recall_trimmed` for each interface in
1362 :attr:`~contact_target_interfaces`
1364 :type: :class:`list` of :class:`float`
1372 """ Same as :attr:`~per_interface_ips` but with :attr:`~trimmed_model`
1374 :attr:`~ics` for each interface in
1375 :attr:`~contact_target_interfaces`
1377 :type: :class:`float`
1386 """ Interfaces in :attr:`~target` that are relevant for DockQ
1388 All interfaces in :attr:`~target` with non-zero contacts that are
1389 relevant for DockQ. Includes protein-protein, protein-nucleotide and
1390 nucleotide-nucleotide interfaces. Chain names for each interface are
1391 lexicographically sorted.
1393 :type: :class:`list` of :class:`tuple` with 2 elements each:
1403 contact_d = contact_d)
1406 interfaces = cent.interacting_chains
1407 interfaces = [(min(x[0],x[1]), max(x[0],x[1]))
for x
in interfaces]
1409 pep_seqs = set([s.GetName()
for s
in self.
chain_mapperchain_mapper.polypep_seqs])
1410 nuc_seqs = set([s.GetName()
for s
in self.
chain_mapperchain_mapper.polynuc_seqs])
1412 seqs = pep_seqs.union(nuc_seqs)
1414 for interface
in interfaces:
1415 if interface[0]
in seqs
and interface[1]
in seqs:
1422 """ Interfaces in :attr:`~dockq_target_interfaces` that can be mapped
1425 Target chain names are lexicographically sorted
1427 :type: :class:`list` of :class:`tuple` with 4 elements each:
1428 (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
1432 flat_mapping = self.
mappingmapping.GetFlatMapping()
1434 if i[0]
in flat_mapping
and i[1]
in flat_mapping:
1437 flat_mapping[i[1]]))
1442 """ DockQ scores for interfaces in :attr:`~dockq_interfaces`
1444 :class:`list` of :class:`float`
1452 """ fnat scores for interfaces in :attr:`~dockq_interfaces`
1454 fnat: Fraction of native contacts that are also present in model
1456 :class:`list` of :class:`float`
1458 if self.
_fnat_fnat
is None:
1460 return self.
_fnat_fnat
1464 """ N native contacts for interfaces in :attr:`~dockq_interfaces`
1466 :class:`list` of :class:`int`
1468 if self.
_nnat_nnat
is None:
1470 return self.
_nnat_nnat
1474 """ N model contacts for interfaces in :attr:`~dockq_interfaces`
1476 :class:`list` of :class:`int`
1478 if self.
_nmdl_nmdl
is None:
1480 return self.
_nmdl_nmdl
1484 """ fnonnat scores for interfaces in :attr:`~dockq_interfaces`
1486 fnat: Fraction of model contacts that are not present in target
1488 :class:`list` of :class:`float`
1496 """ irmsd scores for interfaces in :attr:`~dockq_interfaces`
1498 irmsd: RMSD of interface (RMSD computed on backbone atoms) which
1499 consists of each residue that has at least one heavy atom within 10A of
1500 other chain. Backbone atoms for proteins: "CA","C","N","O", for
1501 nucleotides: "P", "OP1", "OP2", "O2'", "O3'", "O4'", "O5'", "C1'",
1502 "C2'", "C3'", "C4'", "C5'".
1504 :class:`list` of :class:`float`
1506 if self.
_irmsd_irmsd
is None:
1512 """ lrmsd scores for interfaces in :attr:`~dockq_interfaces`
1514 lrmsd: The two chains involved in the interface are superposed based on
1515 the receptor (rigid min RMSD superposition) and the ligand RMSD is
1516 reported. Receptor is the chain with more residues. Superposition and
1517 RMSD is computed on same backbone atoms as :attr:`~irmsd`.
1519 :class:`list` of :class:`float`
1521 if self.
_lrmsd_lrmsd
is None:
1527 """ Average of DockQ scores in :attr:`~dockq_scores`
1529 In its original implementation, DockQ only operates on single
1530 interfaces. Thus the requirement to combine scores for higher order
1533 :type: :class:`float`
1541 """ Same as :attr:`~dockq_ave`, weighted by native contacts
1543 :type: :class:`float`
1551 """ Same as :attr:`~dockq_ave` but penalizing for missing interfaces
1553 Interfaces that are not covered in model are added as 0.0
1554 in average computation.
1556 :type: :class:`float`
1564 """ Same as :attr:`~dockq_ave_full`, but weighted
1566 Interfaces that are not covered in model are added as 0.0 in
1567 average computations and the respective weights are derived from
1568 number of contacts in respective target interface.
1576 """ Mapped representative positions in target
1578 Thats CA positions for peptide residues and C3' positions for
1579 nucleotides. Has same length as :attr:`~mapped_model_pos` and mapping
1580 is based on :attr:`~mapping`.
1582 :type: :class:`ost.geom.Vec3List`
1590 """ Mapped representative positions in target
1592 Thats the equivalent of :attr:`~mapped_target_pos` but containing more
1593 backbone atoms (N, CA, C for peptide residues and O5', C5', C4', C3', O3
1594 for nucleotide residues). mapping is based on :attr:`~mapping`.
1596 :type: :class:`ost.geom.Vec3List`
1604 """ Mapped representative positions in model
1606 Thats CA positions for peptide residues and C3' positions for
1607 nucleotides. Has same length as :attr:`~mapped_target_pos` and mapping
1608 is based on :attr:`~mapping`.
1610 :type: :class:`ost.geom.Vec3List`
1618 """ Mapped representative positions in model
1620 Thats the equivalent of :attr:`~mapped_model_pos` but containing more
1621 backbone atoms (N, CA, C for peptide residues and O5', C5', C4', C3', O3
1622 for nucleotide residues). mapping is based on :attr:`~mapping`.
1624 :type: :class:`ost.geom.Vec3List`
1632 """ :attr:`~mapped_model_pos` with :attr:`~transform` applied
1634 :type: :class:`ost.geom.Vec3List`
1644 """ Number of target residues which have no mapping to model
1654 """ Transform: :attr:`~mapped_model_pos` onto :attr:`~mapped_target_pos`
1656 Computed using Kabsch minimal rmsd algorithm. If number of positions
1657 is too small (< 3), :attr:`~mapped_model_pos_full_bb` and
1658 :attr:`~mapped_target_pos_full_bb` are used.
1660 :type: :class:`ost.geom.Mat4`
1667 self.
_transform_transform = res.transformation
1674 self.
_transform_transform = res.transformation
1679 """ Mapped representative positions in target
1681 Thats CA positions for peptide residues and C3' positions for
1682 nucleotides. Has same length as :attr:`~rigid_mapped_model_pos` and mapping
1683 is based on :attr:`~rigid_mapping`.
1685 :type: :class:`ost.geom.Vec3List`
1693 """ Mapped representative positions in target
1695 Thats the equivalent of :attr:`~rigid_mapped_target_pos` but containing
1696 more backbone atoms (N, CA, C for peptide residues and O5', C5', C4',
1697 C3', O3 for nucleotide residues). mapping is based on :attr:`~mapping`.
1699 :type: :class:`ost.geom.Vec3List`
1707 """ Mapped representative positions in model
1709 Thats CA positions for peptide residues and C3' positions for
1710 nucleotides. Has same length as :attr:`~mapped_target_pos` and mapping
1711 is based on :attr:`~rigid_mapping`.
1713 :type: :class:`ost.geom.Vec3List`
1721 """ Mapped representative positions in model
1723 Thats the equivalent of :attr:`~rigid_mapped_model_pos` but containing
1724 more backbone atoms (N, CA, C for peptide residues and O5', C5', C4',
1725 C3', O3 for nucleotide residues). mapping is based on :attr:`~mapping`.
1727 :type: :class:`ost.geom.Vec3List`
1735 """ :attr:`~rigid_mapped_model_pos` with :attr:`~rigid_transform` applied
1737 :type: :class:`ost.geom.Vec3List`
1747 """ Number of target residues which have no rigid mapping to model
1757 """ Transform: :attr:`~rigid_mapped_model_pos` onto :attr:`~rigid_mapped_target_pos`
1759 Computed using Kabsch minimal rmsd algorithm. If number of positions
1760 is too small (< 3), :attr:`~rigid_mapped_model_pos_full_bb` and
1761 :attr:`~rigid_mapped_target_pos_full_bb` are used.
1763 :type: :class:`ost.geom.Mat4`
1782 """ Fraction CA (C3' for nucleotides) that can be superposed within 0.5A
1784 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1785 Similar iterative algorithm as LGA tool
1787 :type: :class:`float`
1789 if self.
_gdt_05_gdt_05
is None:
1792 for window_size
in wsizes:
1795 window_size, 1000, 0.5)[0]
1800 self.
_gdt_05_gdt_05 = float(n) / n_full
1807 """ Fraction CA (C3' for nucleotides) that can be superposed within 1.0A
1809 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1810 Similar iterative algorithm as LGA tool
1812 :type: :class:`float`
1814 if self.
_gdt_1_gdt_1
is None:
1817 for window_size
in wsizes:
1820 window_size, 1000, 1.0)[0]
1825 self.
_gdt_1_gdt_1 = float(n) / n_full
1832 """ Fraction CA (C3' for nucleotides) that can be superposed within 2.0A
1834 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1835 Similar iterative algorithm as LGA tool
1838 :type: :class:`float`
1840 if self.
_gdt_2_gdt_2
is None:
1843 for window_size
in wsizes:
1846 window_size, 1000, 2.0)[0]
1851 self.
_gdt_2_gdt_2 = float(n) / n_full
1858 """ Fraction CA (C3' for nucleotides) that can be superposed within 4.0A
1860 Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
1861 Similar iterative algorithm as LGA tool
1863 :type: :class:`float`
1865 if self.
_gdt_4_gdt_4
is None:
1868 for window_size
in wsizes:
1871 window_size, 1000, 4.0)[0]
1876 self.
_gdt_4_gdt_4 = float(n) / n_full
1883 """ Fraction CA (C3' for nucleotides) that can be superposed within 8.0A
1885 Similar iterative algorithm as LGA tool
1887 :type: :class:`float`
1889 if self.
_gdt_8_gdt_8
is None:
1892 for window_size
in wsizes:
1895 window_size, 1000, 8.0)[0]
1900 self.
_gdt_8_gdt_8 = float(n) / n_full
1908 """ avg GDT with thresholds: 8.0A, 4.0A, 2.0A and 1.0A
1910 :type: :class:`float`
1912 if self.
_gdtts_gdtts
is None:
1913 LogScript(
"Computing GDT-TS score")
1919 """ avg GDT with thresholds: 4.0A, 2.0A, 1.0A and 0.5A
1921 :type: :class:`float`
1923 if self.
_gdtha_gdtha
is None:
1924 LogScript(
"Computing GDT-HA score")
1932 Computed on :attr:`~rigid_transformed_mapped_model_pos` and
1933 :attr:`~rigid_mapped_target_pos`
1935 :type: :class:`float`
1937 if self.
_rmsd_rmsd
is None:
1938 LogScript(
"Computing RMSD")
1941 return self.
_rmsd_rmsd
1945 """ The global CAD atom-atom (AA) score
1947 Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
1950 :type: :class:`float`
1958 """ The per-residue CAD atom-atom (AA) scores
1960 Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
1963 :type: :class:`dict`
1971 """ Patch QS-scores for each residue in :attr:`~model_interface_residues`
1973 Representative patches for each residue r in chain c are computed as
1976 * mdl_patch_one: All residues in c with CB (CA for GLY) positions within
1977 8A of r and within 12A of residues from any other chain.
1978 * mdl_patch_two: Closest residue x to r in any other chain gets
1979 identified. Patch is then constructed by selecting all residues from
1980 any other chain within 8A of x and within 12A from any residue in c.
1981 * trg_patch_one: Chain name and residue number based mapping from
1983 * trg_patch_two: Chain name and residue number based mapping from
1986 Results are stored in the same manner as
1987 :attr:`~model_interface_residues`, with corresponding scores instead of
1988 residue numbers. Scores for residues which are not
1989 :class:`mol.ChemType.AMINOACIDS` are set to None. Additionally,
1990 interface patches are derived from :attr:`~model`. If they contain
1991 residues which are not covered by :attr:`~target`, the score is set to
1994 :type: :class:`dict` with chain names as key and and :class:`list`
1995 with scores of the respective interface residues.
2003 """ Same as :attr:`~patch_qs` but for DockQ scores
2011 """ TM-score computed with USalign
2013 USalign executable can be specified with usalign_exec kwarg at Scorer
2014 construction, an OpenStructure internal copy of the USalign code is
2017 :type: :class:`float`
2025 """ Mapping computed with USalign
2027 Dictionary with target chain names as key and model chain names as
2028 values. No guarantee that all chains are mapped. USalign executable
2029 can be specified with usalign_exec kwarg at Scorer construction, an
2030 OpenStructure internal copy of the USalign code is used otherwise.
2032 :type: :class:`dict`
2038 def _aln_helper(self, target, model):
2043 for ch
in target.chains:
2044 cname = ch.GetName()
2045 s =
''.join([r.one_letter_code
for r
in ch.residues])
2046 s = seq.CreateSequence(ch.GetName(), s)
2047 s.AttachView(target.Select(f
"cname={mol.QueryQuoteName(cname)}"))
2048 trg_seqs[ch.GetName()] = s
2050 for ch
in model.chains:
2051 cname = ch.GetName()
2052 s =
''.join([r.one_letter_code
for r
in ch.residues])
2053 s = seq.CreateSequence(cname, s)
2054 s.AttachView(model.Select(f
"cname={mol.QueryQuoteName(cname)}"))
2055 mdl_seqs[ch.GetName()] = s
2058 trg_pep_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polypep_seqs]
2059 trg_nuc_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polynuc_seqs]
2060 trg_pep_chains = set(trg_pep_chains)
2061 trg_nuc_chains = set(trg_nuc_chains)
2062 for trg_ch, mdl_ch
in self.
mappingmapping.GetFlatMapping().items():
2063 if mdl_ch
in mdl_seqs
and trg_ch
in trg_seqs:
2064 if trg_ch
in trg_pep_chains:
2065 stype = mol.ChemType.AMINOACIDS
2066 elif trg_ch
in trg_nuc_chains:
2067 stype = mol.ChemType.NUCLEOTIDES
2069 raise RuntimeError(
"Chain name inconsistency... ask "
2071 alns.append(self.
chain_mapperchain_mapper.Align(trg_seqs[trg_ch],
2074 alns[-1].AttachView(0, trg_seqs[trg_ch].GetAttachedView())
2075 alns[-1].AttachView(1, mdl_seqs[mdl_ch].GetAttachedView())
2078 def _compute_pepnuc_aln(self):
2079 query =
"peptide=true or nucleotide=true"
2080 pep_nuc_target = self.
target_origtarget_orig.Select(query)
2081 pep_nuc_model = self.
model_origmodel_orig.Select(query)
2084 def _compute_aln(self):
2087 def _compute_stereochecked_aln(self):
2091 def _compute_lddt(self):
2092 LogScript(
"Computing all-atom LDDT")
2094 flat_mapping = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
2097 stereochecked_alns = dict()
2099 mdl_seq = aln.GetSequence(1)
2100 stereochecked_alns[mdl_seq.name] = aln
2102 for aln
in self.
alnaln:
2103 mdl_seq = aln.GetSequence(1)
2104 alns[mdl_seq.name] = aln
2109 aa_local_lddt =
None
2112 lddt_chain_mapping = dict()
2113 for mdl_ch, trg_ch
in flat_mapping.items():
2115 lddt_chain_mapping[mdl_ch] = trg_ch
2117 chain_mapping = lddt_chain_mapping,
2118 residue_mapping = alns,
2119 check_resnames=
False,
2120 local_lddt_prop=
"lddt",
2122 set_atom_props=
True)[0]
2124 aa_local_lddt = dict()
2125 for r
in self.
modelmodel.residues:
2127 cname = r.GetChain().GetName()
2128 if cname
not in local_lddt:
2129 local_lddt[cname] = dict()
2130 aa_local_lddt[cname] = dict()
2132 rnum = r.GetNumber()
2133 if rnum
not in aa_local_lddt[cname]:
2134 aa_local_lddt[cname][rnum] = dict()
2136 if r.HasProp(
"lddt"):
2137 score = round(r.GetFloatProp(
"lddt"), 3)
2138 local_lddt[cname][rnum] = score
2142 local_lddt[cname][rnum] =
None
2145 if a.HasProp(
"lddt"):
2146 score = round(a.GetFloatProp(
"lddt"), 3)
2147 aa_local_lddt[cname][rnum][a.GetName()] = score
2152 aa_local_lddt[cname][rnum][a.GetName()] =
None
2156 lddt_chain_mapping = dict()
2157 for mdl_ch, trg_ch
in flat_mapping.items():
2158 if mdl_ch
in stereochecked_alns:
2159 lddt_chain_mapping[mdl_ch] = trg_ch
2161 chain_mapping = lddt_chain_mapping,
2162 residue_mapping = stereochecked_alns,
2163 check_resnames=
False,
2164 local_lddt_prop=
"lddt",
2166 set_atom_props=
True)[0]
2168 aa_local_lddt = dict()
2169 for r
in self.
modelmodel.residues:
2170 cname = r.GetChain().GetName()
2171 if cname
not in local_lddt:
2172 local_lddt[cname] = dict()
2173 aa_local_lddt[cname] = dict()
2174 rnum = r.GetNumber()
2175 if rnum
not in aa_local_lddt[cname]:
2176 aa_local_lddt[cname][rnum] = dict()
2178 if r.HasProp(
"lddt"):
2179 score = round(r.GetFloatProp(
"lddt"), 3)
2180 local_lddt[cname][rnum] = score
2186 if a.HasProp(
"lddt"):
2187 score = round(a.GetFloatProp(
"lddt"), 3)
2188 aa_local_lddt[cname][rnum][a.GetName()] = score
2198 if cname
in flat_mapping:
2199 for col
in alns[cname]:
2200 if col[0] !=
'-' and col[1] !=
'-':
2201 if col.GetResidue(1).number == r.number:
2202 trg_r = col.GetResidue(0)
2204 if trg_r
is not None:
2205 trg_cname = trg_r.GetChain().GetName()
2206 trg_rnum = trg_r.GetNumber()
2217 if trg_r
is not None and not trg_r.FindAtom(a.GetName()).IsValid():
2219 aa_local_lddt[cname][rnum][a.GetName()] =
None
2220 elif trg_r
is not None and trg_r.FindAtom(a.GetName()).IsValid()
and \
2221 mdl_r
is not None and not mdl_r.FindAtom(a.GetName()).IsValid():
2223 aa_local_lddt[cname][rnum][a.GetName()] = 0.0
2226 aa_local_lddt[cname][rnum][a.GetName()] =
None
2230 if mdl_res.IsValid():
2233 local_lddt[cname][rnum] =
None
2235 aa_local_lddt[cname][rnum][a.GetName()] =
None
2242 if cname
in flat_mapping:
2243 for col
in alns[cname]:
2244 if col[0] !=
'-' and col[1] !=
'-':
2245 if col.GetResidue(1).number == r.number:
2246 trg_r = col.GetResidue(0)
2248 if trg_r
is not None:
2249 trg_cname = trg_r.GetChain().GetName()
2250 trg_rnum = trg_r.GetNumber()
2257 local_lddt[cname][rnum] =
None
2259 aa_local_lddt[cname][rnum][a.GetName()] =
None
2261 local_lddt[cname][rnum] = 0.0
2263 if trg_r.FindAtom(a.GetName()).IsValid():
2264 aa_local_lddt[cname][rnum][a.GetName()] = 0.0
2266 aa_local_lddt[cname][rnum][a.GetName()] =
None
2268 self.
_lddt_lddt = lddt_score
2272 def _compute_bb_lddt(self):
2273 LogScript(
"Computing backbone LDDT")
2276 for aln
in self.
alnaln:
2277 mdl_seq = aln.GetSequence(1)
2278 alns[mdl_seq.name] = aln
2281 flat_mapping = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
2282 lddt_chain_mapping = dict()
2283 for mdl_ch, trg_ch
in flat_mapping.items():
2285 lddt_chain_mapping[mdl_ch] = trg_ch
2288 chain_mapping = lddt_chain_mapping,
2289 residue_mapping = alns,
2290 check_resnames=
False,
2291 local_lddt_prop=
"bb_lddt",
2294 for r
in self.
modelmodel.residues:
2295 cname = r.GetChain().GetName()
2296 if cname
not in local_lddt:
2297 local_lddt[cname] = dict()
2298 if r.HasProp(
"bb_lddt"):
2299 score = round(r.GetFloatProp(
"bb_lddt"), 3)
2300 local_lddt[cname][r.GetNumber()] = score
2304 local_lddt[cname][r.GetNumber()] =
None
2309 def _compute_ilddt(self):
2310 LogScript(
"Computing all-atom iLDDT")
2312 flat_mapping = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
2316 for aln
in self.
alnaln:
2317 mdl_seq = aln.GetSequence(1)
2318 alns[mdl_seq.name] = aln
2319 lddt_chain_mapping = dict()
2320 for mdl_ch, trg_ch
in flat_mapping.items():
2322 lddt_chain_mapping[mdl_ch] = trg_ch
2324 chain_mapping = lddt_chain_mapping,
2325 residue_mapping = alns,
2326 check_resnames=
False,
2327 local_lddt_prop=
"lddt",
2329 no_intrachain=
True)[0]
2333 mdl_seq = aln.GetSequence(1)
2334 alns[mdl_seq.name] = aln
2335 lddt_chain_mapping = dict()
2336 for mdl_ch, trg_ch
in flat_mapping.items():
2338 lddt_chain_mapping[mdl_ch] = trg_ch
2340 chain_mapping = lddt_chain_mapping,
2341 residue_mapping = alns,
2342 check_resnames=
False,
2343 local_lddt_prop=
"lddt",
2345 no_intrachain=
True)[0]
2348 def _compute_qs(self):
2349 LogScript(
"Computing global QS-score")
2350 qs_score_result = self.
qs_scorerqs_scorer.Score(self.
mappingmapping.mapping)
2351 self.
_qs_global_qs_global = qs_score_result.QS_global
2352 self.
_qs_best_qs_best = qs_score_result.QS_best
2354 def _compute_per_interface_qs_scores(self):
2355 LogScript(
"Computing per-interface QS-score")
2360 trg_ch1 = interface[0]
2361 trg_ch2 = interface[1]
2362 mdl_ch1 = interface[2]
2363 mdl_ch2 = interface[3]
2364 qs_res = self.
qs_scorerqs_scorer.ScoreInterface(trg_ch1, trg_ch2,
2369 def _compute_ics_scores(self):
2370 LogScript(
"Computing ICS scores")
2373 self.
_ics_recall_ics_recall = contact_scorer_res.recall
2374 self.
_ics_ics = contact_scorer_res.ics
2378 flat_mapping = self.
mappingmapping.GetFlatMapping()
2380 trg_ch1 = trg_int[0]
2381 trg_ch2 = trg_int[1]
2382 if trg_ch1
in flat_mapping
and trg_ch2
in flat_mapping:
2383 mdl_ch1 = flat_mapping[trg_ch1]
2384 mdl_ch2 = flat_mapping[trg_ch2]
2385 res = self.
contact_scorercontact_scorer.ScoreICSInterface(trg_ch1, trg_ch2,
2395 def _trim_model(self):
2396 trimmed_mdl = mol.CreateEntityFromView(self.
mappingmapping.model,
True)
2397 trimmed_aln = dict()
2399 for trg_cname, mdl_cname
in self.
mappingmapping.GetFlatMapping().items():
2400 mdl_ch = trimmed_mdl.FindChain(mdl_cname)
2401 aln = self.
mappingmapping.alns[(trg_cname, mdl_cname)]
2404 assert(mdl_ch.GetResidueCount() == \
2405 len(aln.GetSequence(1).GetGaplessString()))
2407 mdl_residues = mdl_ch.residues
2409 aligned_mdl_seq = [
'-'] * aln.GetLength()
2410 for col_idx, col
in enumerate(aln):
2411 if col[0] !=
'-' and col[1] !=
'-':
2412 mdl_res = mdl_residues[mdl_res_idx]
2413 mdl_res.SetIntProp(
"aligned", 1)
2414 aligned_mdl_seq[col_idx] = col[1]
2417 aligned_mdl_seq =
''.join(aligned_mdl_seq)
2418 aligned_mdl_seq = seq.CreateSequence(mdl_cname, aligned_mdl_seq)
2420 new_aln = seq.CreateAlignment()
2421 new_aln.AddSequence(aln.GetSequence(0))
2422 new_aln.AddSequence(aligned_mdl_seq)
2423 trimmed_aln[(trg_cname, mdl_cname)] = new_aln
2425 self.
_trimmed_model_trimmed_model = trimmed_mdl.Select(
"graligned:0=1")
2428 def _compute_ics_scores_trimmed(self):
2429 LogScript(
"Computing ICS scores trimmed")
2443 flat_mapping = self.
mappingmapping.GetFlatMapping()
2445 trg_ch1 = trg_int[0]
2446 trg_ch2 = trg_int[1]
2447 if trg_ch1
in flat_mapping
and trg_ch2
in flat_mapping:
2448 mdl_ch1 = flat_mapping[trg_ch1]
2449 mdl_ch2 = flat_mapping[trg_ch2]
2460 def _compute_ips_scores(self):
2461 LogScript(
"Computing IPS scores")
2464 self.
_ips_recall_ips_recall = contact_scorer_res.recall
2465 self.
_ips_ips = contact_scorer_res.ips
2470 flat_mapping = self.
mappingmapping.GetFlatMapping()
2472 trg_ch1 = trg_int[0]
2473 trg_ch2 = trg_int[1]
2474 if trg_ch1
in flat_mapping
and trg_ch2
in flat_mapping:
2475 mdl_ch1 = flat_mapping[trg_ch1]
2476 mdl_ch2 = flat_mapping[trg_ch2]
2477 res = self.
contact_scorercontact_scorer.ScoreIPSInterface(trg_ch1, trg_ch2,
2487 def _compute_ips_scores_trimmed(self):
2488 LogScript(
"Computing IPS scores trimmed")
2501 flat_mapping = self.
mappingmapping.GetFlatMapping()
2503 trg_ch1 = trg_int[0]
2504 trg_ch2 = trg_int[1]
2505 if trg_ch1
in flat_mapping
and trg_ch2
in flat_mapping:
2506 mdl_ch1 = flat_mapping[trg_ch1]
2507 mdl_ch2 = flat_mapping[trg_ch2]
2518 def _compute_dockq_scores(self):
2519 LogScript(
"Computing DockQ")
2522 raise RuntimeError(
"Cannot compute DockQ for reference structures "
2523 "with nucleotide chains if dockq_capri_peptide "
2528 self.
_fnat_fnat = list()
2530 self.
_irmsd_irmsd = list()
2531 self.
_lrmsd_lrmsd = list()
2532 self.
_nnat_nnat = list()
2533 self.
_nmdl_nmdl = list()
2536 for aln
in self.
alnaln:
2537 trg_s = aln.GetSequence(0)
2538 mdl_s = aln.GetSequence(1)
2539 dockq_alns[(trg_s.GetName(), mdl_s.GetName())] = aln
2542 trg_ch1 = interface[0]
2543 trg_ch2 = interface[1]
2544 mdl_ch1 = interface[2]
2545 mdl_ch2 = interface[3]
2546 aln1 = dockq_alns[(trg_ch1, mdl_ch1)]
2547 aln2 = dockq_alns[(trg_ch2, mdl_ch2)]
2549 res = dockq.DockQ(self.
modelmodel, self.
targettarget, mdl_ch1, mdl_ch2,
2550 trg_ch1, trg_ch2, ch1_aln=aln1,
2551 ch2_aln=aln2, contact_dist_thresh = 4.0,
2552 interface_dist_thresh=8.0,
2553 interface_cb =
True)
2555 res = dockq.DockQ(self.
modelmodel, self.
targettarget, mdl_ch1, mdl_ch2,
2556 trg_ch1, trg_ch2, ch1_aln=aln1,
2559 self.
_fnat_fnat.append(res[
"fnat"])
2560 self.
_fnonnat_fnonnat.append(res[
"fnonnat"])
2561 self.
_irmsd_irmsd.append(res[
"irmsd"])
2562 self.
_lrmsd_lrmsd.append(res[
"lrmsd"])
2564 self.
_nnat_nnat.append(res[
"nnat"])
2565 self.
_nmdl_nmdl.append(res[
"nmdl"])
2570 not_covered_counts = list()
2571 proc_trg_interfaces = set([(x[0], x[1])
for x
in self.
dockq_interfacesdockq_interfaces])
2573 if interface
not in proc_trg_interfaces:
2577 trg_ch1 = interface[0]
2578 trg_ch2 = interface[1]
2581 res = dockq.DockQ(self.
targettarget, self.
targettarget,
2582 trg_ch1, trg_ch2, trg_ch1, trg_ch2,
2583 contact_dist_thresh = 4.0,
2584 interface_dist_thresh=8.0,
2585 interface_cb =
True)
2587 res = dockq.DockQ(self.
targettarget, self.
targettarget,
2588 trg_ch1, trg_ch2, trg_ch1, trg_ch2)
2590 not_covered_counts.append(res[
"nnat"])
2597 weights = np.array(self.
_nnat_nnat)
2602 self.
_dockq_wave_dockq_wave = np.sum(np.multiply(weights/np.sum(weights), scores))
2603 scores = np.append(scores, [0.0]*len(not_covered_counts))
2604 weights = np.append(weights, not_covered_counts)
2609 self.
_dockq_wave_full_dockq_wave_full = np.sum(np.multiply(weights/np.sum(weights),
2612 def _extract_mapped_pos(self):
2616 processed_trg_chains = set()
2617 for trg_ch, mdl_ch
in self.
mappingmapping.GetFlatMapping().items():
2618 processed_trg_chains.add(trg_ch)
2619 aln = self.
mappingmapping.alns[(trg_ch, mdl_ch)]
2621 if col[0] !=
'-' and col[1] !=
'-':
2622 trg_res = col.GetResidue(0)
2623 mdl_res = col.GetResidue(1)
2624 trg_at = trg_res.FindAtom(
"CA")
2625 mdl_at = mdl_res.FindAtom(
"CA")
2626 if not trg_at.IsValid():
2627 trg_at = trg_res.FindAtom(
"C3'")
2628 if not mdl_at.IsValid():
2629 mdl_at = mdl_res.FindAtom(
"C3'")
2635 for ch
in self.
mappingmapping.target.chains:
2636 if ch.GetName()
not in processed_trg_chains:
2639 def _extract_mapped_pos_full_bb(self):
2642 exp_pep_atoms = [
"N",
"CA",
"C"]
2643 exp_nuc_atoms = [
"\"O5'\"",
"\"C5'\"",
"\"C4'\"",
"\"C3'\"",
"\"O3'\""]
2644 trg_pep_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polypep_seqs]
2645 trg_nuc_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polynuc_seqs]
2646 for trg_ch, mdl_ch
in self.
mappingmapping.GetFlatMapping().items():
2647 aln = self.
mappingmapping.alns[(trg_ch, mdl_ch)]
2648 trg_ch = aln.GetSequence(0).GetName()
2649 if trg_ch
in trg_pep_chains:
2650 exp_atoms = exp_pep_atoms
2651 elif trg_ch
in trg_nuc_chains:
2652 exp_atoms = exp_nuc_atoms
2655 raise RuntimeError(
"Unexpected error - contact OST developer")
2657 if col[0] !=
'-' and col[1] !=
'-':
2658 trg_res = col.GetResidue(0)
2659 mdl_res = col.GetResidue(1)
2660 for aname
in exp_atoms:
2661 trg_at = trg_res.FindAtom(aname)
2662 mdl_at = mdl_res.FindAtom(aname)
2663 if not (trg_at.IsValid()
and mdl_at.IsValid()):
2665 raise RuntimeError(
"Unexpected error - contact OST "
2671 def _extract_rigid_mapped_pos(self):
2675 processed_trg_chains = set()
2676 for trg_ch, mdl_ch
in self.
rigid_mappingrigid_mapping.GetFlatMapping().items():
2677 processed_trg_chains.add(trg_ch)
2678 aln = self.
rigid_mappingrigid_mapping.alns[(trg_ch, mdl_ch)]
2680 if col[0] !=
'-' and col[1] !=
'-':
2681 trg_res = col.GetResidue(0)
2682 mdl_res = col.GetResidue(1)
2683 trg_at = trg_res.FindAtom(
"CA")
2684 mdl_at = mdl_res.FindAtom(
"CA")
2685 if not trg_at.IsValid():
2686 trg_at = trg_res.FindAtom(
"C3'")
2687 if not mdl_at.IsValid():
2688 mdl_at = mdl_res.FindAtom(
"C3'")
2695 if ch.GetName()
not in processed_trg_chains:
2698 def _extract_rigid_mapped_pos_full_bb(self):
2701 exp_pep_atoms = [
"N",
"CA",
"C"]
2702 exp_nuc_atoms = [
"\"O5'\"",
"\"C5'\"",
"\"C4'\"",
"\"C3'\"",
"\"O3'\""]
2703 trg_pep_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polypep_seqs]
2704 trg_nuc_chains = [s.GetName()
for s
in self.
chain_mapperchain_mapper.polynuc_seqs]
2705 for trg_ch, mdl_ch
in self.
rigid_mappingrigid_mapping.GetFlatMapping().items():
2706 aln = self.
mappingmapping.alns[(trg_ch, mdl_ch)]
2707 trg_ch = aln.GetSequence(0).GetName()
2708 if trg_ch
in trg_pep_chains:
2709 exp_atoms = exp_pep_atoms
2710 elif trg_ch
in trg_nuc_chains:
2711 exp_atoms = exp_nuc_atoms
2714 raise RuntimeError(
"Unexpected error - contact OST developer")
2716 if col[0] !=
'-' and col[1] !=
'-':
2717 trg_res = col.GetResidue(0)
2718 mdl_res = col.GetResidue(1)
2719 for aname
in exp_atoms:
2720 trg_at = trg_res.FindAtom(aname)
2721 mdl_at = mdl_res.FindAtom(aname)
2722 if not (trg_at.IsValid()
and mdl_at.IsValid()):
2724 raise RuntimeError(
"Unexpected error - contact OST "
2729 def _compute_cad_score(self):
2731 raise RuntimeError(
"CAD score computations rely on residue numbers "
2732 "that are consistent between target and model "
2733 "chains, i.e. only work if resnum_alignments "
2734 "is True at Scorer construction.")
2736 LogScript(
"Computing CAD score")
2738 settings.Locate(
"voronota-cadscore",
2740 except Exception
as e:
2741 raise RuntimeError(
"voronota-cadscore must be in PATH for CAD "
2742 "score scoring")
from e
2743 cad_bin_dir = os.path.dirname(cad_score_exec)
2744 m = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
2745 cad_result = cadscore.CADScore(self.
modelmodel, self.
targettarget,
2749 cad_bin_path=cad_bin_dir,
2753 for r
in self.
modelmodel.residues:
2754 cname = r.GetChain().GetName()
2755 if cname
not in local_cad:
2756 local_cad[cname] = dict()
2757 if r.HasProp(
"localcad"):
2758 score = round(r.GetFloatProp(
"localcad"), 3)
2759 local_cad[cname][r.GetNumber()] = score
2761 local_cad[cname][r.GetNumber()] =
None
2763 self.
_cad_score_cad_score = cad_result.globalAA
2766 def _get_repr_view(self, ent):
2767 """ Returns view with representative peptide atoms => CB, CA for GLY
2769 Ensures that each residue has exactly one atom with assertions
2771 :param ent: Entity for which you want the representative view
2772 :param ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
2773 :returns: :class:`ost.mol.EntityView` derived from ent
2775 repr_ent = ent.Select(
"(aname=\"CB\" or (rname=\"GLY\" and aname=\"CA\"))")
2776 for r
in repr_ent.residues:
2777 assert(len(r.atoms) == 1)
2780 def _get_interface_residues(self, ent):
2781 """ Get interface residues
2783 Thats all residues having a contact with at least one residue from another
2784 chain (CB-CB distance <= 8A, CA in case of Glycine)
2786 :param ent: Model for which to extract interface residues
2787 :type ent: :class:`ost.mol.EntityView`
2788 :returns: :class:`dict` with chain names as key and and :class:`list`
2789 with residue numbers of the respective interface residues.
2793 result = {ch.GetName(): list()
for ch
in ent.chains}
2794 for ch
in ent.chains:
2795 cname = ch.GetName()
2796 sel = repr_ent.Select(f
"(cname={mol.QueryQuoteName(cname)} and 8 <> [cname!={mol.QueryQuoteName(cname)}])")
2797 result[cname] = [r.GetNumber()
for r
in sel.residues]
2800 def _do_stereochecks(self):
2801 """ Perform stereochemistry checks on model and target
2803 LogInfo(
"Performing stereochemistry checks on model and target")
2804 data = stereochemistry.GetDefaultStereoData()
2805 l_data = stereochemistry.GetDefaultStereoLinkData()
2807 a, b, c, d = stereochemistry.StereoCheck(self.
modelmodel, stereo_data = data,
2808 stereo_link_data = l_data)
2814 a, b, c, d = stereochemistry.StereoCheck(self.
targettarget, stereo_data = data,
2815 stereo_link_data = l_data)
2822 def _get_interface_patches(self, mdl_ch, mdl_rnum):
2823 """ Select interface patches representative for specified residue
2825 The patches for specified residue r in chain c are selected as follows:
2827 * mdl_patch_one: All residues in c with CB (CA for GLY) positions within 8A
2828 of r and within 12A of residues from any other chain.
2829 * mdl_patch_two: Closest residue x to r in any other chain gets identified.
2830 Patch is then constructed by selecting all residues from any other chain
2831 within 8A of x and within 12A from any residue in c.
2832 * trg_patch_one: Chain name and residue number based mapping from
2834 * trg_patch_two: Chain name and residue number based mapping from
2837 :param mdl_ch: Name of chain in *self.model* of residue of interest
2838 :type mdl_ch: :class:`str`
2839 :param mdl_rnum: Residue number of residue of interest
2840 :type mdl_rnum: :class:`ost.mol.ResNum`
2841 :returns: Tuple with 5 elements: 1) :class:`bool` flag whether all residues
2842 in *mdl* patches are covered in *trg* 2) mtl_patch_one
2843 3) mdl_patch_two 4) trg_patch_one 5) trg_patch_two
2849 r = self.
modelmodel.FindResidue(mdl_ch, mdl_rnum)
2851 raise RuntimeError(f
"Cannot find residue {mdl_rnum} in chain {mdl_ch}")
2852 if r.GetName() ==
"GLY":
2853 at = r.FindAtom(
"CA")
2855 at = r.FindAtom(
"CB")
2856 if not at.IsValid():
2857 raise RuntimeError(
"Cannot find interface views for res without CB/CA")
2864 q1 = f
"(cname={mol.QueryQuoteName(mdl_ch)} and 8 <> {{{r_pos[0]},{r_pos[1]},{r_pos[2]}}})"
2866 q2 = f
"(12 <> [cname!={mol.QueryQuoteName(mdl_ch)}])"
2867 mdl_patch_one = self.
modelmodel.CreateEmptyView()
2868 sel = repr_mdl.Select(
" and ".join([q1, q2]))
2869 for r
in sel.residues:
2870 mdl_r = self.
modelmodel.FindResidue(r.GetChain().GetName(), r.GetNumber())
2871 mdl_patch_one.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
2877 sel = repr_mdl.Select(f
"(cname!={mol.QueryQuoteName(mdl_ch)})")
2878 close_stuff = sel.FindWithin(r_pos, 8)
2881 for close_at
in close_stuff:
2884 min_pos = close_at.GetPos()
2888 q1 = f
"(cname!={mol.QueryQuoteName(mdl_ch)} and 8 <> {{{min_pos[0]},{min_pos[1]},{min_pos[2]}}})"
2890 q2 = f
"(12 <> [cname={mol.QueryQuoteName(mdl_ch)}])"
2891 mdl_patch_two = self.
modelmodel.CreateEmptyView()
2892 sel = repr_mdl.Select(
" and ".join([q1, q2]))
2893 for r
in sel.residues:
2894 mdl_r = self.
modelmodel.FindResidue(r.GetChain().GetName(), r.GetNumber())
2895 mdl_patch_two.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
2898 flat_mapping = self.
mappingmapping.GetFlatMapping(mdl_as_key=
True)
2899 full_trg_coverage =
True
2900 trg_patch_one = self.
targettarget.CreateEmptyView()
2901 for r
in mdl_patch_one.residues:
2903 mdl_cname = r.GetChain().GetName()
2904 if mdl_cname
in flat_mapping:
2905 aln = self.
mappingmapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
2907 if col[0] !=
'-' and col[1] !=
'-':
2908 if col.GetResidue(1).GetNumber() == r.GetNumber():
2909 trg_r = col.GetResidue(0)
2911 if trg_r
is not None:
2912 trg_patch_one.AddResidue(trg_r.handle,
2913 mol.ViewAddFlag.INCLUDE_ALL)
2915 full_trg_coverage =
False
2917 trg_patch_two = self.
targettarget.CreateEmptyView()
2918 for r
in mdl_patch_two.residues:
2920 mdl_cname = r.GetChain().GetName()
2921 if mdl_cname
in flat_mapping:
2922 aln = self.
mappingmapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
2924 if col[0] !=
'-' and col[1] !=
'-':
2925 if col.GetResidue(1).GetNumber() == r.GetNumber():
2926 trg_r = col.GetResidue(0)
2928 if trg_r
is not None:
2929 trg_patch_two.AddResidue(trg_r.handle,
2930 mol.ViewAddFlag.INCLUDE_ALL)
2932 full_trg_coverage =
False
2934 return (full_trg_coverage, mdl_patch_one, mdl_patch_two,
2935 trg_patch_one, trg_patch_two)
2937 def _compute_patchqs_scores(self):
2938 LogScript(
"Computing patch QS-scores")
2944 r = self.
modelmodel.FindResidue(cname, rnum)
2945 if r.IsValid()
and r.GetChemType() == mol.ChemType.AMINOACIDS:
2946 full_trg_coverage, mdl_patch_one, mdl_patch_two, \
2947 trg_patch_one, trg_patch_two = \
2949 if full_trg_coverage:
2950 score = self.
_patchqs_patchqs(mdl_patch_one, mdl_patch_two,
2951 trg_patch_one, trg_patch_two)
2952 scores.append(score)
2955 def _compute_patchdockq_scores(self):
2956 LogScript(
"Computing patch DockQ scores")
2962 r = self.
modelmodel.FindResidue(cname, rnum)
2963 if r.IsValid()
and r.GetChemType() == mol.ChemType.AMINOACIDS:
2964 full_trg_coverage, mdl_patch_one, mdl_patch_two, \
2965 trg_patch_one, trg_patch_two = \
2967 if full_trg_coverage:
2968 score = self.
_patchdockq_patchdockq(mdl_patch_one, mdl_patch_two,
2969 trg_patch_one, trg_patch_two)
2970 scores.append(score)
2973 def _patchqs(self, mdl_patch_one, mdl_patch_two, trg_patch_one, trg_patch_two):
2974 """ Score interface residue patches with QS-score
2976 In detail: Construct two entities with two chains each. First chain
2977 consists of residues from <x>_patch_one and second chain consists of
2978 <x>_patch_two. The returned score is the QS-score between the two
2981 :param mdl_patch_one: Interface patch representing scored residue
2982 :type mdl_patch_one: :class:`ost.mol.EntityView`
2983 :param mdl_patch_two: Interface patch representing scored residue
2984 :type mdl_patch_two: :class:`ost.mol.EntityView`
2985 :param trg_patch_one: Interface patch representing scored residue
2986 :type trg_patch_one: :class:`ost.mol.EntityView`
2987 :param trg_patch_two: Interface patch representing scored residue
2988 :type trg_patch_two: :class:`ost.mol.EntityView`
2989 :returns: PatchQS score
2994 alnA = seq.CreateAlignment()
2995 s =
''.join([r.one_letter_code
for r
in mdl_patch_one.residues])
2996 alnA.AddSequence(seq.CreateSequence(
"A", s))
2997 s =
''.join([r.one_letter_code
for r
in trg_patch_one.residues])
2998 alnA.AddSequence(seq.CreateSequence(
"A", s))
3000 alnB = seq.CreateAlignment()
3001 s =
''.join([r.one_letter_code
for r
in mdl_patch_two.residues])
3002 alnB.AddSequence(seq.CreateSequence(
"B", s))
3003 s =
''.join([r.one_letter_code
for r
in trg_patch_two.residues])
3004 alnB.AddSequence(seq.CreateSequence(
"B", s))
3005 alns = {(
"A",
"A"): alnA, (
"B",
"B"): alnB}
3007 scorer =
QSScorer(qs_ent_mdl, [[
"A"], [
"B"]], qs_ent_trg, alns)
3008 score_result = scorer.Score([[
"A"], [
"B"]])
3010 return score_result.QS_global
3012 def _patchdockq(self, mdl_patch_one, mdl_patch_two, trg_patch_one,
3014 """ Score interface residue patches with DockQ
3016 In detail: Construct two entities with two chains each. First chain
3017 consists of residues from <x>_patch_one and second chain consists of
3018 <x>_patch_two. The returned score is the QS-score between the two
3021 :param mdl_patch_one: Interface patch representing scored residue
3022 :type mdl_patch_one: :class:`ost.mol.EntityView`
3023 :param mdl_patch_two: Interface patch representing scored residue
3024 :type mdl_patch_two: :class:`ost.mol.EntityView`
3025 :param trg_patch_one: Interface patch representing scored residue
3026 :type trg_patch_one: :class:`ost.mol.EntityView`
3027 :param trg_patch_two: Interface patch representing scored residue
3028 :type trg_patch_two: :class:`ost.mol.EntityView`
3029 :returns: DockQ score
3033 dockq_result = dockq.DockQ(t, m,
"A",
"B",
"A",
"B")
3034 if dockq_result[
"nnat"] > 0:
3035 return dockq_result[
"DockQ"]
3038 def _qs_ent_from_patches(self, patch_one, patch_two):
3039 """ Constructs Entity with two chains named "A" and "B""
3041 Blindly adds all residues from *patch_one* to chain A and residues from
3042 patch_two to chain B.
3044 ent = mol.CreateEntity()
3046 added_ch = ed.InsertChain(
"A")
3047 for r
in patch_one.residues:
3048 added_r = ed.AppendResidue(added_ch, r.GetName())
3049 added_r.SetChemClass(str(r.GetChemClass()))
3051 ed.InsertAtom(added_r, a.handle)
3052 added_ch = ed.InsertChain(
"B")
3053 for r
in patch_two.residues:
3054 added_r = ed.AppendResidue(added_ch, r.GetName())
3055 added_r.SetChemClass(str(r.GetChemClass()))
3057 ed.InsertAtom(added_r, a.handle)
3060 def _construct_custom_mapping(self, mapping):
3061 """ constructs a full blown MappingResult object from simple dict
3063 :param mapping: mapping with trg chains as key and mdl ch as values
3064 :type mapping: :class:`dict`
3066 LogInfo(
"Setting custom chain mapping")
3069 chem_mapping, chem_group_alns, unmapped_mdl_chains, mdl = \
3070 chain_mapper.GetChemMapping(self.
modelmodel)
3075 if len(mapping) != len(set(mapping.keys())):
3076 raise RuntimeError(f
"Expect unique trg chain names in mapping. Got "
3077 f
"{mapping.keys()}")
3078 if len(mapping) != len(set(mapping.values())):
3079 raise RuntimeError(f
"Expect unique mdl chain names in mapping. Got "
3080 f
"{mapping.values()}")
3082 trg_chains = set([ch.GetName()
for ch
in chain_mapper.target.chains])
3083 mdl_chains = set([ch.GetName()
for ch
in mdl.chains])
3084 for k,v
in mapping.items():
3085 if k
not in trg_chains:
3086 raise RuntimeError(f
"Target chain \"{k}\" is not available "
3087 f
"in target processed for chain mapping "
3089 if v
not in mdl_chains:
3090 raise RuntimeError(f
"Model chain \"{v}\" is not available "
3091 f
"in model processed for chain mapping "
3094 for trg_ch, mdl_ch
in mapping.items():
3095 trg_group_idx =
None
3096 mdl_group_idx =
None
3097 for idx, group
in enumerate(chain_mapper.chem_groups):
3101 for idx, group
in enumerate(chem_mapping):
3105 if trg_group_idx
is None or mdl_group_idx
is None:
3106 raise RuntimeError(
"Could not establish a valid chem grouping "
3107 "of chain names provided in custom mapping.")
3109 if trg_group_idx != mdl_group_idx:
3110 raise RuntimeError(f
"Chem group mismatch in custom mapping: "
3111 f
"target chain \"{trg_ch}\" groups with the "
3112 f
"following chemically equivalent target "
3114 f
"{chain_mapper.chem_groups[trg_group_idx]} "
3115 f
"but model chain \"{mdl_ch}\" maps to the "
3116 f
"following target chains: "
3117 f
"{chain_mapper.chem_groups[mdl_group_idx]}")
3119 pairs = set([(trg_ch, mdl_ch)
for trg_ch, mdl_ch
in mapping.items()])
3121 chain_mapping._GetRefMdlAlns(chain_mapper.chem_groups,
3122 chain_mapper.chem_group_alignments,
3128 final_mapping = list()
3129 for ref_chains
in chain_mapper.chem_groups:
3130 mapped_mdl_chains = list()
3131 for ref_ch
in ref_chains:
3132 if ref_ch
in mapping:
3133 mapped_mdl_chains.append(mapping[ref_ch])
3135 mapped_mdl_chains.append(
None)
3136 final_mapping.append(mapped_mdl_chains)
3139 for ref_group, mdl_group
in zip(chain_mapper.chem_groups,
3141 for ref_ch, mdl_ch
in zip(ref_group, mdl_group):
3142 if ref_ch
is not None and mdl_ch
is not None:
3143 aln = ref_mdl_alns[(ref_ch, mdl_ch)]
3144 trg_view = chain_mapper.target.Select(f
"cname={mol.QueryQuoteName(ref_ch)}")
3145 mdl_view = mdl.Select(f
"cname={mol.QueryQuoteName(mdl_ch)}")
3146 aln.AttachView(0, trg_view)
3147 aln.AttachView(1, mdl_view)
3148 alns[(ref_ch, mdl_ch)] = aln
3151 chain_mapper.chem_groups,
3153 unmapped_mdl_chains,
3154 final_mapping, alns)
3156 def _compute_tmscore(self):
3159 LogScript(
"Computing TM-score with built-in USalign")
3161 flat_mapping = self.
rigid_mappingrigid_mapping.GetFlatMapping()
3162 LogInfo(
"Overriding TM-score chain mapping")
3163 res = res = bindings.WrappedMMAlign(self.
modelmodel, self.
targettarget,
3164 mapping=flat_mapping)
3166 res = bindings.WrappedMMAlign(self.
modelmodel, self.
targettarget)
3168 LogScript(
"Computing TM-score with USalign exectuable")
3170 LogInfo(
"Overriding TM-score chain mapping")
3171 flat_mapping = self.
rigid_mappingrigid_mapping.GetFlatMapping()
3172 res = tmtools.USAlign(self.
modelmodel, self.
targettarget,
3174 custom_chain_mapping = flat_mapping)
3176 res = tmtools.USAlign(self.
modelmodel, self.
targettarget,
3181 for a,b
in zip(res.ent1_mapped_chains, res.ent2_mapped_chains):
3184 def _resnum_connect(self, ent):
3186 for ch
in ent.chains:
3187 res_list = ch.residues
3188 for i
in range(len(res_list) - 1):
3191 if ra.GetNumber().GetNum() + 1 == rb.GetNumber().GetNum():
3192 if ra.IsPeptideLinking()
and rb.IsPeptideLinking():
3193 c = ra.FindAtom(
"C")
3194 n = rb.FindAtom(
"N")
3195 if c.IsValid()
and n.IsValid()
and not mol.BondExists(c, n):
3197 ed = ent.EditXCS(mol.BUFFERED_EDIT)
3199 elif ra.IsNucleotideLinking()
and rb.IsNucleotideLinking():
3200 o = ra.FindAtom(
"O3'")
3201 p = rb.FindAtom(
"P")
3202 if o.IsValid()
and p.IsValid()
and not mol.BondExists(o, p):
3204 ed = ent.EditXCS(mol.BUFFERED_EDIT)
3209 __all__ = (
'lDDTBSScorer',
'Scorer',)
def target_interface_residues(self)
_per_interface_ips_trimmed
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)
_rigid_mapped_target_pos_full_bb
_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)
_mapped_model_pos_full_bb
def dockq_wave_full(self)
def mapped_model_pos(self)
def rigid_mapped_target_pos_full_bb(self)
_per_interface_ics_precision_trimmed
def _compute_dockq_scores(self)
def _aln_helper(self, target, model)
def model_bad_bonds(self)
def per_interface_ics_recall(self)
def ics_precision_trimmed(self)
def _compute_bb_lddt(self)
def transformed_mapped_model_pos(self)
_per_interface_ips_recall_trimmed
_per_interface_ics_recall_trimmed
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 rigid_transformed_mapped_model_pos(self)
def trimmed_model_contacts(self)
def _compute_ips_scores_trimmed(self)
def n_target_not_mapped(self)
def per_interface_ips_trimmed(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)
_rigid_mapped_model_pos_full_bb
def dockq_interfaces(self)
def per_interface_ips_recall(self)
def per_interface_ics_precision(self)
def _compute_ics_scores(self)
def _resnum_connect(self, ent)
_rigid_n_target_not_mapped
def _get_interface_residues(self, ent)
def ips_precision_trimmed(self)
def rigid_transform(self)
_mapped_target_pos_full_bb
def ics_recall_trimmed(self)
def mapped_target_pos_full_bb(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)
def trimmed_contact_scorer(self)
def per_interface_ics_recall_trimmed(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 rigid_mapped_model_pos_full_bb(self)
def stereochecked_model(self)
_model_interface_residues
def _compute_pepnuc_aln(self)
def per_interface_qs_best(self)
def _extract_rigid_mapped_pos_full_bb(self)
def per_interface_ips_precision_trimmed(self)
def _compute_ics_scores_trimmed(self)
_per_interface_ips_precision_trimmed
def rigid_n_target_not_mapped(self)
def per_interface_ics_precision_trimmed(self)
def __init__(self, model, target, resnum_alignments=False, molck_settings=None, cad_score_exec=None, custom_mapping=None, custom_rigid_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, dockq_capri_peptide=False, lddt_symmetry_settings=None, lddt_inclusion_radius=15.0, pep_seqid_thr=95., nuc_seqid_thr=95., mdl_map_pep_seqid_thr=0., mdl_map_nuc_seqid_thr=0.)
def per_interface_ips_recall_trimmed(self)
def mapped_target_pos(self)
def mapped_model_pos_full_bb(self)
def target_bad_bonds(self)
def qs_model_interfaces(self)
def _compute_cad_score(self)
def ips_recall_trimmed(self)
def per_interface_ics_trimmed(self)
_per_interface_ics_trimmed
def _construct_custom_mapping(self, mapping)
def _extract_mapped_pos_full_bb(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)