diff --git a/modules/mol/alg/pymod/lddt.py b/modules/mol/alg/pymod/lddt.py index 58e423491429c502836118cdfa68bd71eb9b89a7..c14ea43a20fea23d1c9fe1801ed774125fe788c8 100644 --- a/modules/mol/alg/pymod/lddt.py +++ b/modules/mol/alg/pymod/lddt.py @@ -121,6 +121,11 @@ class lDDTScorer: original residue numbers were. :type seqres_mapping: :class:`dict` (key: :class:`str`, value: :class:`ost.seq.AlignmentHandle`) + :param calpha: Only consider atoms with name "CA". Technically this sets + the expected atom names for each residue name to ["CA"], thus + invalidating *compound_lib*. No check whether the target + residues are actually amino acids! + :type calpha: :class:`bool` :raises: :class:`RuntimeError` if *target* contains compound which is not in *compound_lib*, :class:`RuntimeError` if *symmetry_settings* specifies symmetric atoms that are not present in the according @@ -132,7 +137,6 @@ class lDDTScorer: (seq1 contains gaps, mismatch in seq1/seq2, seq2 does not match residues in corresponding chains). """ - def __init__( self, target, @@ -141,6 +145,7 @@ class lDDTScorer: sequence_separation=0, symmetry_settings=None, seqres_mapping=dict(), + calpha=False ): self.target = target @@ -152,6 +157,10 @@ class lDDTScorer: else: self.symmetry_settings = symmetry_settings + # whether to only consider atoms with name "CA", invalidates + # *compound_lib* + self.calpha=calpha + # names of heavy atoms of each unique compound present in *target* as # extracted from *compound_lib*, e.g. # self.compound_anames["GLY"] = ["N", "CA", "C", "O"] @@ -202,7 +211,7 @@ class lDDTScorer: # setup members defined above self._SetupEnv(self.compound_lib, self.symmetry_settings, - seqres_mapping) + seqres_mapping, self.calpha) # distance related members are lazily computed as they're affected # by different flavours of lDDT (e.g. lDDT including inter-chain @@ -560,7 +569,8 @@ class lDDTScorer: dummy_scorer = lDDTScorer(model.Select("cname="+ch_name), self.compound_lib, symmetry_settings = sm, - inclusion_radius = self.inclusion_radius) + inclusion_radius = self.inclusion_radius, + calpha = self.calpha) penalty += dummy_scorer.n_distances return penalty @@ -622,7 +632,8 @@ class lDDTScorer: return rnums - def _SetupEnv(self, compound_lib, symmetry_settings, seqres_mapping): + def _SetupEnv(self, compound_lib, symmetry_settings, seqres_mapping, + calpha): """Sets target related lDDTScorer members defined in constructor No distance related members - see _SetupDistances @@ -682,31 +693,10 @@ class lDDTScorer: self.chain_res_start_indices.append(len(self.compound_names)) for r, rnum in zip(chain.residues, residue_numbers[ch_name]): if r.name not in self.compound_anames: - compound = compound_lib.FindCompound(r.name) - if compound is None: - raise RuntimeError(f"compound_lib has no entry for {r}") - atom_names = list() - for atom_spec in compound.GetAtomSpecs(): - if atom_spec.element not in ["H", "D"]: - atom_names.append(atom_spec.name) - self.compound_anames[r.name] = atom_names - self.compound_symmetric_atoms[r.name] = list() - if r.name in symmetry_settings.symmetric_compounds: - for pair in symmetry_settings.symmetric_compounds[ - r.name - ]: - try: - a = self.compound_anames[r.name].index(pair[0]) - b = self.compound_anames[r.name].index(pair[1]) - except: - msg = f"Could not find symmetric atoms " - msg += f"({pair[0]}, {pair[1]}) for {r.name} " - msg += f"as specified in SymmetrySettings in " - msg += f"compound from component dictionary. " - msg += f"Atoms in compound: " - msg += f"{self.compound_anames[r.name]}" - raise RuntimeError(msg) - self.compound_symmetric_atoms[r.name].append((a, b)) + # sets compound info in self.compound_anames and + # self.compound_symmetric_atoms + self._SetupCompound(r, compound_lib, symmetry_settings, + calpha) self.res_start_indices.append(current_idx) self.res_mapper[(ch_name, rnum)] = len(self.compound_names) @@ -729,6 +719,39 @@ class lDDTScorer: ) self.n_atoms = current_idx + def _SetupCompound(self, r, compound_lib, symmetry_settings, calpha): + """fill self.compound_anames/self.compound_symmetric_atoms + """ + if calpha: + # throw away compound_lib info + self.compound_anames[r.name] = ["CA"] + self.compound_symmetric_atoms[r.name] = list() + else: + atom_names = list() + symmetric_atoms = list() + compound = compound_lib.FindCompound(r.name) + if compound is None: + raise RuntimeError(f"no entry for {r} in compound_lib") + for atom_spec in compound.GetAtomSpecs(): + if atom_spec.element not in ["H", "D"]: + atom_names.append(atom_spec.name) + self.compound_anames[r.name] = atom_names + if r.name in symmetry_settings.symmetric_compounds: + for pair in symmetry_settings.symmetric_compounds[r.name]: + try: + a = self.compound_anames[r.name].index(pair[0]) + b = self.compound_anames[r.name].index(pair[1]) + except: + msg = f"Could not find symmetric atoms " + msg += f"({pair[0]}, {pair[1]}) for {r.name} " + msg += f"as specified in SymmetrySettings in " + msg += f"compound from component dictionary. " + msg += f"Atoms in compound: " + msg += f"{self.compound_anames[r.name]}" + raise RuntimeError(msg) + symmetric_atoms.append((a, b)) + self.compound_symmetric_atoms[r.name] = symmetric_atoms + def _SetupDistances(self): """Compute distance related members of lDDTScorer """ diff --git a/modules/mol/alg/tests/test_lddt.py b/modules/mol/alg/tests/test_lddt.py index 0264b55f0a968f84d8bdee9698bb69814379d72e..0a0e365b1d2dff552ccfb4d45adc2ecae3f40776 100644 --- a/modules/mol/alg/tests/test_lddt.py +++ b/modules/mol/alg/tests/test_lddt.py @@ -148,6 +148,30 @@ class TestlDDT(unittest.TestCase): scorer = lDDTScorer(target, conop.GetDefaultLib(), sequence_separation=0) + def test_calpha(self): + model = _LoadFile("7SGN_C_model.pdb") + target = _LoadFile("7SGN_C_target.pdb") + + # do scoring and select aname=CA + scorer = lDDTScorer(target.Select("aname=CA"), + conop.GetDefaultLib()) + score_one, per_res_scores_one = scorer.lDDT(model) + score_two, per_res_scores_two = scorer.lDDT(model.Select("aname=CA")) + + # no selection, just setting calpha flag should give the same + scorer = lDDTScorer(target, + conop.GetDefaultLib(), + calpha=True) + score_three, per_res_scores_three = scorer.lDDT(model) + + # check + self.assertAlmostEqual(score_one, score_two, places=5) + self.assertAlmostEqual(score_one, score_three, places=5) + for a,b in zip(per_res_scores_one, per_res_scores_two): + self.assertAlmostEqual(a, b, places=5) + for a,b in zip(per_res_scores_one, per_res_scores_three): + self.assertAlmostEqual(a, b, places=5) + if __name__ == "__main__": from ost import testutils if testutils.SetDefaultCompoundLib():