Skip to content
Snippets Groups Projects
Commit 79eb5516 authored by Studer Gabriel's avatar Studer Gabriel
Browse files

lDDT: add flag to comput CA-lDDT

parent 93edac50
No related branches found
No related tags found
No related merge requests found
......@@ -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
"""
......
......@@ -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():
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment