From 184652ac7d52e31f47f9d8435eec29cd16dddeb1 Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Tue, 25 Feb 2025 22:23:37 +0100
Subject: [PATCH] chain_mapping/scoring refactoring/bugfixes
- avoids requirement of having views attached to each and every sequence
that comes out of chain_mapper. Original motivation of this was that
we're more flexible when it comes to residue number alignments, especially
when we want to provide custom SEQRES. As an additional treat, we get some
performance increase. In all-vs-all chain mapping scenarios where all
structures are already in memory, the creation of all these views requires
significant resources.
- unit test of Scorer objects
- bugfix in IPS: The reported per-interface ips precision/recall values
were wrong. The actual per-interface ips score as well as the global
ips score were not affected.
- bugfix in lDDT: The scorer doesn't use the alignments that come out
of the chain mapper but create their own. One set of alignments for
the full model/target and one set of alignments for stereochecked
model/target. These alignments can be inconsistent if residues in
the stereochemistry checks are removed at unfortunate locations
(beginning/end) or simply if a lot of residues are removed. The
stereochecked alns are not independently computed anymore, but
start now from the full alignments and just introduce gaps for
residues that are completely removed. This does not affect scoring
runs with residue number alignments enabled.
---
actions/ost-compare-structures | 94 +++-
modules/mol/alg/pymod/chain_mapping.py | 514 ++++++++------------
modules/mol/alg/pymod/contact_score.py | 24 +-
modules/mol/alg/pymod/scoring.py | 367 +++++++++-----
modules/mol/alg/tests/CMakeLists.txt | 3 +-
modules/mol/alg/tests/test_chain_mapping.py | 49 +-
modules/mol/alg/tests/test_scoring.py | 242 +++++++++
7 files changed, 819 insertions(+), 474 deletions(-)
create mode 100644 modules/mol/alg/tests/test_scoring.py
diff --git a/actions/ost-compare-structures b/actions/ost-compare-structures
index cd2fc76aa..649522d99 100644
--- a/actions/ost-compare-structures
+++ b/actions/ost-compare-structures
@@ -781,22 +781,48 @@ def _AlnToFastaStr(aln):
s2 = aln.GetSequence(1)
return f">reference:{s1.name}\n{str(s1)}\n>model:{s2.name}\n{str(s2)}"
-def _GetInconsistentResidues(alns):
+def _GetInconsistentResidues(alns, ref, mdl):
lst = list()
for aln in alns:
+
+ ref_ch = ref.FindChain(aln.GetSequence(0).GetName())
+ mdl_ch = mdl.FindChain(aln.GetSequence(1).GetName())
+
+ if not ref_ch.IsValid():
+ raise RuntimeError("ref lacks requested chain in _GetInconsistentResidues")
+
+ if not mdl_ch.IsValid():
+ raise RuntimeError("mdl lacks requested chain in _GetInconsistentResidues")
+
+ if len(aln.GetSequence(0).GetGaplessString()) != ref_ch.GetResidueCount():
+ raise RuntimeError("aln/chain mismatch in _GetAlignedResidues")
+ if len(aln.GetSequence(1).GetGaplessString()) != mdl_ch.GetResidueCount():
+ raise RuntimeError("aln/chain mismatch in _GetAlignedResidues")
+
+ ref_res = ref_ch.residues
+ mdl_res = mdl_ch.residues
+
+ ref_res_idx = 0
+ mdl_res_idx = 0
+
for col in aln:
- r1 = col.GetResidue(0)
- r2 = col.GetResidue(1)
- if r1.IsValid() and r2.IsValid() and r1.GetName() != r2.GetName():
- ch_1 = r1.GetChain().name
- num_1 = r1.number.num
- ins_code_1 = r1.number.ins_code.strip("\u0000")
- id_1 = f"{ch_1}.{num_1}.{ins_code_1}"
- ch_2 = r2.GetChain().name
- num_2 = r2.number.num
- ins_code_2 = r2.number.ins_code.strip("\u0000")
- id_2 = f"{ch_2}.{num_2}.{ins_code_2}"
- lst.append(f"{id_1}-{id_2}")
+ if col[0] != '-' and col[1] != '-':
+ r1 = ref_res[ref_res_idx]
+ r2 = mdl_res[mdl_res_idx]
+ if r1.IsValid() and r2.IsValid() and r1.GetName() != r2.GetName():
+ ch_1 = r1.GetChain().name
+ num_1 = r1.number.num
+ ins_code_1 = r1.number.ins_code.strip("\u0000")
+ id_1 = f"{ch_1}.{num_1}.{ins_code_1}"
+ ch_2 = r2.GetChain().name
+ num_2 = r2.number.num
+ ins_code_2 = r2.number.ins_code.strip("\u0000")
+ id_2 = f"{ch_2}.{num_2}.{ins_code_2}"
+ lst.append(f"{id_1}-{id_2}")
+ if col[0] != '-':
+ ref_res_idx += 1
+ if col[1] != '-':
+ mdl_res_idx += 1
return lst
def _LocalScoresToJSONDict(score_dict):
@@ -842,29 +868,52 @@ def _PatchScoresToJSONList(interface_dict, score_dict):
json_list.append(_RoundOrNone(item))
return json_list
-def _GetAlignedResidues(aln):
+def _GetAlignedResidues(aln, ref, mdl):
aligned_residues = list()
for a in aln:
mdl_lst = list()
ref_lst = list()
+
+ ref_ch = ref.FindChain(a.GetSequence(0).GetName())
+ mdl_ch = mdl.FindChain(a.GetSequence(1).GetName())
+
+ if not ref_ch.IsValid():
+ raise RuntimeError("ref lacks requested chain in _GetAlignedResidues")
+
+ if not mdl_ch.IsValid():
+ raise RuntimeError("mdl lacks requested chain in _GetAlignedResidues")
+
+ if len(a.GetSequence(0).GetGaplessString()) != ref_ch.GetResidueCount():
+ raise RuntimeError("aln/chain mismatch in _GetAlignedResidues")
+ if len(a.GetSequence(1).GetGaplessString()) != mdl_ch.GetResidueCount():
+ raise RuntimeError("aln/chain mismatch in _GetAlignedResidues")
+
+ ref_res = ref_ch.residues
+ mdl_res = mdl_ch.residues
+
+ ref_res_idx = 0
+ mdl_res_idx = 0
+
for c in a:
- mdl_r = c.GetResidue(1)
- ref_r = c.GetResidue(0)
- if mdl_r.IsValid():
+ if c[1] != '-':
+ mdl_r = mdl_res[mdl_res_idx]
olc = mdl_r.one_letter_code
num = mdl_r.GetNumber().num
ins_code = mdl_r.GetNumber().ins_code.strip("\u0000")
mdl_lst.append({"olc": olc,
"num": f"{num}.{ins_code}"})
+ mdl_res_idx += 1
else:
mdl_lst.append(None)
- if ref_r.IsValid():
+ if c[0] != '-':
+ ref_r = ref_res[ref_res_idx]
olc = ref_r.one_letter_code
num = ref_r.GetNumber().num
ins_code = ref_r.GetNumber().ins_code.strip("\u0000")
ref_lst.append({"olc": olc,
"num": f"{num}.{ins_code}"})
+ ref_res_idx += 1
else:
ref_lst.append(None)
@@ -901,7 +950,7 @@ def _Process(model, reference, args, model_format, reference_format):
mdl_map_pep_seqid_thr = args.chem_map_seqid_thresh,
mdl_map_nuc_seqid_thr = args.chem_map_seqid_thresh)
- ir = _GetInconsistentResidues(scorer.aln)
+ ir = _GetInconsistentResidues(scorer.aln, scorer.target, scorer.model)
if len(ir) > 0 and args.enforce_consistency:
raise RuntimeError(f"Inconsistent residues observed: {' '.join(ir)}")
@@ -916,13 +965,16 @@ def _Process(model, reference, args, model_format, reference_format):
out["inconsistent_residues"] = ir
if args.dump_aligned_residues:
- out["aligned_residues"] = _GetAlignedResidues(scorer.aln)
+ out["aligned_residues"] = _GetAlignedResidues(scorer.aln, scorer.target,
+ scorer.model)
if args.dump_pepnuc_alns:
out["pepnuc_aln"] = [_AlnToFastaStr(aln) for aln in scorer.pepnuc_aln]
if args.dump_pepnuc_aligned_residues:
- out["pepnuc_aligned_residues"] = _GetAlignedResidues(scorer.pepnuc_aln)
+ out["pepnuc_aligned_residues"] = _GetAlignedResidues(scorer.pepnuc_aln,
+ scorer.pepnuc_target,
+ scorer.pepnuc_model)
if args.lddt:
out["lddt"] = _RoundOrNone(scorer.lddt)
diff --git a/modules/mol/alg/pymod/chain_mapping.py b/modules/mol/alg/pymod/chain_mapping.py
index 9347d3161..04d278745 100644
--- a/modules/mol/alg/pymod/chain_mapping.py
+++ b/modules/mol/alg/pymod/chain_mapping.py
@@ -113,8 +113,7 @@ class MappingResult:
Each alignment is accessible with ``alns[(t_chain,m_chain)]``. First
sequence is the sequence of :attr:`target` chain, second sequence the
- one from :attr:`~model`. The respective :class:`ost.mol.EntityView` are
- attached with :func:`ost.seq.ConstSequenceHandle.AttachView`.
+ one from :attr:`~model`.
:type: :class:`dict` with key: :class:`tuple` of :class:`str`, value:
:class:`ost.seq.AlignmentHandle`
@@ -588,6 +587,12 @@ class ChainMapper:
self.n_max_naive = n_max_naive
self.mdl_map_pep_seqid_thr = mdl_map_pep_seqid_thr
self.mdl_map_nuc_seqid_thr = mdl_map_nuc_seqid_thr
+ self.pep_subst_mat = pep_subst_mat
+ self.pep_gap_open = pep_gap_open
+ self.pep_gap_ext = pep_gap_ext
+ self.nuc_subst_mat = nuc_subst_mat
+ self.nuc_gap_open = nuc_gap_open
+ self.nuc_gap_ext = nuc_gap_ext
# lazy computed attributes
self._chem_groups = None
@@ -595,15 +600,6 @@ class ChainMapper:
self._chem_group_ref_seqs = None
self._chem_group_types = None
- # helper class to generate pairwise alignments
- self.aligner = _Aligner(resnum_aln = resnum_alignments,
- pep_subst_mat = pep_subst_mat,
- pep_gap_open = pep_gap_open,
- pep_gap_ext = pep_gap_ext,
- nuc_subst_mat = nuc_subst_mat,
- nuc_gap_open = nuc_gap_open,
- nuc_gap_ext = nuc_gap_ext)
-
# target structure preprocessing
self._target, self._polypep_seqs, self._polynuc_seqs = \
self.ProcessStructure(target)
@@ -669,14 +665,7 @@ class ChainMapper:
:type: :class:`ost.seq.AlignmentList`
"""
if self._chem_group_alignments is None:
- self._chem_group_alignments, self._chem_group_types = \
- _GetChemGroupAlignments(self.polypep_seqs, self.polynuc_seqs,
- self.aligner,
- pep_seqid_thr=self.pep_seqid_thr,
- min_pep_length=self.min_pep_length,
- nuc_seqid_thr=self.nuc_seqid_thr,
- min_nuc_length=self.min_nuc_length)
-
+ self._SetChemGroupAlignments()
return self._chem_group_alignments
@property
@@ -694,7 +683,6 @@ class ChainMapper:
for a in self.chem_group_alignments:
s = seq.CreateSequence(a.GetSequence(0).GetName(),
a.GetSequence(0).GetGaplessString())
- s.AttachView(a.GetSequence(0).GetAttachedView())
self._chem_group_ref_seqs.AddSequence(s)
return self._chem_group_ref_seqs
@@ -710,14 +698,7 @@ class ChainMapper:
:type: :class:`list` of :class:`ost.mol.ChemType`
"""
if self._chem_group_types is None:
- self._chem_group_alignments, self._chem_group_types = \
- _GetChemGroupAlignments(self.polypep_seqs, self.polynuc_seqs,
- self.aligner,
- pep_seqid_thr=self.pep_seqid_thr,
- min_pep_length=self.min_pep_length,
- nuc_seqid_thr=self.nuc_seqid_thr,
- min_nuc_length=self.min_nuc_length)
-
+ self._SetChemGroupAlignments()
return self._chem_group_types
def GetChemMapping(self, model):
@@ -744,12 +725,11 @@ class ChainMapper:
alns = [seq.AlignmentList() for x in self.chem_groups]
for s in mdl_pep_seqs:
- idx, aln = _MapSequence(self.chem_group_ref_seqs,
- self.chem_group_types,
- s, mol.ChemType.AMINOACIDS,
- self.aligner,
- seq_id_thr = self.mdl_map_pep_seqid_thr,
- min_aln_length = self.min_pep_length)
+ idx, aln = self._MapSequence(self.chem_group_ref_seqs,
+ self.chem_group_types,
+ s, mol.ChemType.AMINOACIDS, mdl,
+ seq_id_thr = self.mdl_map_pep_seqid_thr,
+ min_aln_length = self.min_pep_length)
if idx is None:
unmapped_mdl_chains.append(s.GetName())
else:
@@ -757,12 +737,11 @@ class ChainMapper:
alns[idx].append(aln)
for s in mdl_nuc_seqs:
- idx, aln = _MapSequence(self.chem_group_ref_seqs,
- self.chem_group_types,
- s, mol.ChemType.NUCLEOTIDES,
- self.aligner,
- seq_id_thr = self.mdl_map_nuc_seqid_thr,
- min_aln_length = self.min_nuc_length)
+ idx, aln = self._MapSequence(self.chem_group_ref_seqs,
+ self.chem_group_types,
+ s, mol.ChemType.NUCLEOTIDES, mdl,
+ seq_id_thr = self.mdl_map_nuc_seqid_thr,
+ min_aln_length = self.min_nuc_length)
if idx is None:
unmapped_mdl_chains.append(s.GetName())
else:
@@ -878,8 +857,6 @@ class ChainMapper:
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
- aln.AttachView(0, _CSel(self.target, [ref_ch]))
- aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
unmapped_mdl_chains, one_to_one, alns)
@@ -921,8 +898,6 @@ class ChainMapper:
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
- aln.AttachView(0, _CSel(self.target, [ref_ch]))
- aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
@@ -1018,8 +993,6 @@ class ChainMapper:
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
- aln.AttachView(0, _CSel(self.target, [ref_ch]))
- aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
unmapped_mdl_chains, one_to_one, alns)
@@ -1062,8 +1035,6 @@ class ChainMapper:
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
- aln.AttachView(0, _CSel(self.target, [ref_ch]))
- aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
@@ -1138,8 +1109,6 @@ class ChainMapper:
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
- aln.AttachView(0, _CSel(self.target, [ref_ch]))
- aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
unmapped_mdl_chains, one_to_one, alns)
@@ -1210,8 +1179,6 @@ class ChainMapper:
for ref_ch, mdl_ch in zip(ref_group, mdl_group):
if ref_ch is not None and mdl_ch is not None:
aln = ref_mdl_alns[(ref_ch, mdl_ch)]
- aln.AttachView(0, _CSel(self.target, [ref_ch]))
- aln.AttachView(1, _CSel(mdl, [mdl_ch]))
alns[(ref_ch, mdl_ch)] = aln
return MappingResult(self.target, mdl, self.chem_groups, chem_mapping,
@@ -1600,7 +1567,6 @@ class ChainMapper:
s = ''.join([r.one_letter_code for r in ch.residues])
s = seq.CreateSequence(ch.GetName(), s)
- s.AttachView(_CSel(view, [ch.GetName()]))
if n_pep == n_res:
polypep_seqs.AddSequence(s)
elif n_nuc == n_res:
@@ -1616,16 +1582,17 @@ class ChainMapper:
f"- mapping failed")
# select for chains for which we actually extracted the sequence
- chain_names = [s.GetAttachedView().chains[0].name for s in polypep_seqs]
- chain_names += [s.GetAttachedView().chains[0].name for s in polynuc_seqs]
+ chain_names = [s.name for s in polypep_seqs]
+ chain_names += [s.name for s in polynuc_seqs]
view = _CSel(view, chain_names)
return (view, polypep_seqs, polynuc_seqs)
- def Align(self, s1, s2, stype):
+ def NWAlign(self, s1, s2, stype):
""" Access to internal sequence alignment functionality
- Alignment parameterization is setup at ChainMapper construction
+ Performs Needleman-Wunsch alignment with parameterization
+ setup at ChainMapper construction
:param s1: First sequence to align - must have view attached in case
of resnum_alignments
@@ -1638,59 +1605,11 @@ class ChainMapper:
:class:`ost.mol.ChemType.NUCLEOTIDES`]
:returns: Pairwise alignment of s1 and s2
"""
- if stype not in [mol.ChemType.AMINOACIDS, mol.ChemType.NUCLEOTIDES]:
- raise RuntimeError("stype must be ost.mol.ChemType.AMINOACIDS or "
- "ost.mol.ChemType.NUCLEOTIDES")
- return self.aligner.Align(s1, s2, chem_type = stype)
-
-
-# INTERNAL HELPERS
-##################
-class _Aligner:
- def __init__(self, pep_subst_mat = seq.alg.BLOSUM62, pep_gap_open = -5,
- pep_gap_ext = -2, nuc_subst_mat = seq.alg.NUC44,
- nuc_gap_open = -4, nuc_gap_ext = -4, resnum_aln = False):
- """ Helper class to compute alignments
-
- Sets default values for substitution matrix, gap open and gap extension
- penalties. They are only used in default mode (Needleman-Wunsch aln).
- If *resnum_aln* is True, only residue numbers of views that are attached
- to input sequences are considered.
- """
- self.pep_subst_mat = pep_subst_mat
- self.pep_gap_open = pep_gap_open
- self.pep_gap_ext = pep_gap_ext
- self.nuc_subst_mat = nuc_subst_mat
- self.nuc_gap_open = nuc_gap_open
- self.nuc_gap_ext = nuc_gap_ext
- self.resnum_aln = resnum_aln
-
- def Align(self, s1, s2, chem_type=None):
- if self.resnum_aln:
- return self.ResNumAlign(s1, s2)
- else:
- if chem_type is None:
- raise RuntimeError("Must specify chem_type for NW alignment")
- return self.NWAlign(s1, s2, chem_type)
-
- def NWAlign(self, s1, s2, chem_type):
- """ Returns pairwise alignment using Needleman-Wunsch algorithm
-
- :param s1: First sequence to align
- :type s1: :class:`ost.seq.SequenceHandle`
- :param s2: Second sequence to align
- :type s2: :class:`ost.seq.SequenceHandle`
- :param chem_type: Must be in [:class:`ost.mol.ChemType.AMINOACIDS`,
- :class:`ost.mol.ChemType.NUCLEOTIDES`], determines
- substitution matrix and gap open/extension penalties
- :type chem_type: :class:`ost.mol.ChemType`
- :returns: Alignment with s1 as first and s2 as second sequence
- """
- if chem_type == mol.ChemType.AMINOACIDS:
+ if stype == mol.ChemType.AMINOACIDS:
return seq.alg.SemiGlobalAlign(s1, s2, self.pep_subst_mat,
gap_open=self.pep_gap_open,
gap_ext=self.pep_gap_ext)[0]
- elif chem_type == mol.ChemType.NUCLEOTIDES:
+ elif stype == mol.ChemType.NUCLEOTIDES:
return seq.alg.SemiGlobalAlign(s1, s2, self.nuc_subst_mat,
gap_open=self.nuc_gap_open,
gap_ext=self.nuc_gap_ext)[0]
@@ -1698,45 +1617,214 @@ class _Aligner:
raise RuntimeError("Invalid ChemType")
return aln
- def ResNumAlign(self, s1, s2):
- """ Returns pairwise alignment using residue numbers of attached views
-
- Assumes that there are no insertion codes (alignment only on numerical
- component) and that resnums are strictly increasing (fast min/max
- identification). These requirements are assured if a structure has been
- processed by :class:`ChainMapper`.
+ def ResNumAlign(self, s1, s2, s1_ent, s2_ent):
+ """ Access to internal sequence alignment functionality
+
+ Performs residue number based alignment. Residue numbers are extracted
+ from *s1_ent*/*s2_ent*.
- :param s1: First sequence to align, must have :class:`ost.mol.EntityView`
- attached
+ :param s1: First sequence to align
:type s1: :class:`ost.seq.SequenceHandle`
- :param s2: Second sequence to align, must have :class:`ost.mol.EntityView`
- attached
+ :param s2: Second sequence to align
:type s2: :class:`ost.seq.SequenceHandle`
+ :param s1_ent: Structure as source for residue numbers in *s1*. Must
+ contain a chain named after sequence name in *s1*. This
+ chain must have the exact same number of residues as
+ characters in s1.
+ :type s1_ent: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
+ :param s2_ent: Same for *s2*.
+ :type s2_ent: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
+
+ :returns: Pairwise alignment of s1 and s2
"""
- assert(s1.HasAttachedView())
- assert(s2.HasAttachedView())
- v1 = s1.GetAttachedView()
- rnums1 = [r.GetNumber().GetNum() for r in v1.residues]
- v2 = s2.GetAttachedView()
- rnums2 = [r.GetNumber().GetNum() for r in v2.residues]
+ ch1 = s1_ent.FindChain(s1.name)
+ ch2 = s2_ent.FindChain(s2.name)
+
+ if not ch1.IsValid():
+ raise RuntimeError("s1_ent lacks requested chain in ResNumAlign")
+
+ if not ch2.IsValid():
+ raise RuntimeError("s2_ent lacks requested chain in ResNumAlign")
+
+ if len(s1) != ch1.GetResidueCount():
+ raise RuntimeError("Sequence/structure mismatch in ResNumAlign")
+
+ if len(s2) != ch2.GetResidueCount():
+ raise RuntimeError("Sequence/structure mismatch in ResNumAlign")
+
+ rnums1 = [r.GetNumber().GetNum() for r in ch1.residues]
+ rnums2 = [r.GetNumber().GetNum() for r in ch2.residues]
min_num = min(rnums1[0], rnums2[0])
max_num = max(rnums1[-1], rnums2[-1])
aln_length = max_num - min_num + 1
aln_s1 = ['-'] * aln_length
- for r, rnum in zip(v1.residues, rnums1):
- aln_s1[rnum-min_num] = r.one_letter_code
+ for olc, rnum in zip(s1, rnums1):
+ aln_s1[rnum-min_num] = olc
aln_s2 = ['-'] * aln_length
- for r, rnum in zip(v2.residues, rnums2):
- aln_s2[rnum-min_num] = r.one_letter_code
+ for olc, rnum in zip(s2, rnums2):
+ aln_s2[rnum-min_num] = olc
aln = seq.CreateAlignment()
aln.AddSequence(seq.CreateSequence(s1.GetName(), ''.join(aln_s1)))
aln.AddSequence(seq.CreateSequence(s2.GetName(), ''.join(aln_s2)))
return aln
+ def _SetChemGroupAlignments(self):
+ """Sets self._chem_group_alignments and self._chem_group_types
+ """
+ pep_groups = self._GroupSequences(self.polypep_seqs, self.pep_seqid_thr,
+ self.min_pep_length,
+ mol.ChemType.AMINOACIDS)
+ nuc_groups = self._GroupSequences(self.polynuc_seqs, self.nuc_seqid_thr,
+ self.min_nuc_length,
+ mol.ChemType.NUCLEOTIDES)
+ group_types = [mol.ChemType.AMINOACIDS] * len(pep_groups)
+ group_types += [mol.ChemType.NUCLEOTIDES] * len(nuc_groups)
+ groups = pep_groups
+ groups.extend(nuc_groups)
+ self._chem_group_alignments = groups
+ self._chem_group_types = group_types
+
+ def _GroupSequences(self, seqs, seqid_thr, min_length, chem_type):
+ """Get list of alignments representing groups of equivalent sequences
+
+ :param seqid_thr: Threshold used to decide when two chains are identical.
+ :type seqid_thr: :class:`float`
+ :param gap_thr: Additional threshold to avoid gappy alignments with high
+ seqid. Number of aligned columns must be at least this
+ number.
+ :type gap_thr: :class:`int`
+ :param aligner: Helper class to generate pairwise alignments
+ :type aligner: :class:`_Aligner`
+ :param chem_type: ChemType of seqs, must be in
+ [:class:`ost.mol.ChemType.AMINOACIDS`,
+ :class:`ost.mol.ChemType.NUCLEOTIDES`]
+ :type chem_type: :class:`ost.mol.ChemType`
+ :returns: A list of alignments, one alignment for each group
+ with longest sequence (reference) as first sequence.
+ :rtype: :class:`ost.seq.AlignmentList`
+ """
+ groups = list()
+ for s_idx in range(len(seqs)):
+ matching_group = None
+ for g_idx in range(len(groups)):
+ for g_s_idx in range(len(groups[g_idx])):
+ if self.resnum_alignments:
+ aln = self.ResNumAlign(seqs[s_idx],
+ seqs[groups[g_idx][g_s_idx]],
+ self.target, self.target)
+ else:
+ aln = self.NWAlign(seqs[s_idx],
+ seqs[groups[g_idx][g_s_idx]],
+ chem_type)
+ sid, n_aligned = _GetAlnPropsOne(aln)
+ if sid >= seqid_thr and n_aligned >= min_length:
+ matching_group = g_idx
+ break
+ if matching_group is not None:
+ break
+
+ if matching_group is None:
+ groups.append([s_idx])
+ else:
+ groups[matching_group].append(s_idx)
+
+ # sort based on sequence length
+ sorted_groups = list()
+ for g in groups:
+ if len(g) > 1:
+ tmp = sorted([[len(seqs[i]), i] for i in g], reverse=True)
+ sorted_groups.append([x[1] for x in tmp])
+ else:
+ sorted_groups.append(g)
+
+ # translate from indices back to sequences and directly generate alignments
+ # of the groups with the longest (first) sequence as reference
+ aln_list = seq.AlignmentList()
+ for g in sorted_groups:
+ if len(g) == 1:
+ # aln with one single sequence
+ aln_list.append(seq.CreateAlignment(seqs[g[0]]))
+ else:
+ # obtain pairwise aln of first sequence (reference) to all others
+ alns = seq.AlignmentList()
+ i = g[0]
+ for j in g[1:]:
+ if self.resnum_alignments:
+ aln = self.ResNumAlign(seqs[i], seqs[j],
+ self.target, self.target)
+ else:
+ aln = self.NWAlign(seqs[i], seqs[j], chem_type)
+ alns.append(aln)
+ # and merge
+ aln_list.append(seq.alg.MergePairwiseAlignments(alns, seqs[i]))
+
+ return aln_list
+
+ def _MapSequence(self, ref_seqs, ref_types, s, s_type, s_ent,
+ seq_id_thr=0.0, min_aln_length=0):
+ """Tries top map *s* onto any of the sequences in *ref_seqs*
+
+ Computes alignments of *s* to each of the reference sequences of equal type
+ and sorts them by seqid*fraction_covered (seqid: sequence identity of
+ aligned columns in alignment, fraction_covered: Fraction of non-gap
+ characters in reference sequence that are covered by non-gap characters in
+ *s*). Best scoring mapping is returned. Optionally, *seq_id*/
+ *min_aln_length* thresholds can be enabled to avoid non-sensical mappings.
+ However, *min_aln_length* only has an effect if *seq_id_thr* > 0!!!
+
+ :param ref_seqs: Reference sequences
+ :type ref_seqs: :class:`ost.seq.SequenceList`
+ :param ref_types: Types of reference sequences, e.g.
+ ost.mol.ChemType.AminoAcids
+ :type ref_types: :class:`list` of :class:`ost.mol.ChemType`
+ :param s: Sequence to map
+ :type s: :class:`ost.seq.SequenceHandle`
+ :param s_type: Type of *s*, only try mapping to sequences in *ref_seqs*
+ with equal type as defined in *ref_types*
+ :param s_ent: Entity which represents *s*. Only relevant in case of
+ residue number alignments.
+ :type s_ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
+ :param seq_id_thr: Minimum sequence identity to be considered as match
+ :type seq_id_thr: :class:`float`
+ :param min_aln_length: Minimum number of aligned columns to be considered
+ as match. Only has an effect if *seq_id_thr* > 0!
+ :type min_aln_length: :class:`int`
+ :returns: Tuple with two elements. 1) index of sequence in *ref_seqs* to
+ which *s* can be mapped 2) Pairwise sequence alignment with
+ sequence from *ref_seqs* as first sequence. Both elements are
+ None if no mapping can be found or if thresholds are not
+ fulfilled for any alignment.
+ :raises: :class:`RuntimeError` if mapping is ambiguous, i.e. *s*
+ successfully maps to more than one sequence in *ref_seqs*
+ """
+ scored_alns = list()
+ for ref_idx, ref_seq in enumerate(ref_seqs):
+ if ref_types[ref_idx] == s_type:
+ if self.resnum_alignments:
+ aln = self.ResNumAlign(ref_seq, s, self.target, s_ent)
+ else:
+ aln = self.NWAlign(ref_seq, s, s_type)
+ seqid, n_tot, n_aligned = _GetAlnPropsTwo(aln)
+ if seq_id_thr > 0:
+ if seqid >= seq_id_thr and n_aligned >= min_aln_length:
+ fraction_covered = float(n_aligned)/n_tot
+ score = seqid * fraction_covered
+ scored_alns.append((score, ref_idx, aln))
+ else:
+ fraction_covered = float(n_aligned)/n_tot
+ score = seqid * fraction_covered
+ scored_alns.append((score, ref_idx, aln))
+
+ if len(scored_alns) == 0:
+ return (None, None) # no mapping possible...
+
+ scored_alns = sorted(scored_alns, key=lambda x: x[0], reverse=True)
+ return (scored_alns[0][1], scored_alns[0][2])
+
def _GetAlnPropsTwo(aln):
"""Returns basic properties of *aln* version two...
@@ -1765,176 +1853,6 @@ def _GetAlnPropsOne(aln):
n_aligned = sum([1 for col in aln if (col[0] != '-' and col[1] != '-')])
return (seqid, n_aligned)
-def _GetChemGroupAlignments(pep_seqs, nuc_seqs, aligner, pep_seqid_thr=95.,
- min_pep_length=6, nuc_seqid_thr=95.,
- min_nuc_length=4):
- """Returns alignments with groups of chemically equivalent chains
-
- :param pep_seqs: List of polypeptide sequences
- :type pep_seqs: :class:`seq.SequenceList`
- :param nuc_seqs: List of polynucleotide sequences
- :type nuc_seqs: :class:`seq.SequenceList`
- :param aligner: Helper class to generate pairwise alignments
- :type aligner: :class:`_Aligner`
- :param pep_seqid_thr: Threshold used to decide when two peptide chains are
- identical. 95 percent tolerates the few mutations
- crystallographers like to do.
- :type pep_seqid_thr: :class:`float`
- :param min_pep_length: Additional threshold to avoid gappy alignments with high
- seqid. Number of aligned columns must be at least this
- number.
- :type min_pep_length: :class:`int`
- :param nuc_seqid_thr: Nucleotide equivalent of *pep_seqid_thr*
- :type nuc_seqid_thr: :class:`float`
- :param min_nuc_length: Nucleotide equivalent of *min_pep_length*
- :type min_nuc_length: :class:`int`
- :returns: Tuple with first element being an AlignmentList. Each alignment
- represents a group of chemically equivalent chains and the first
- sequence is the longest. Second element is a list of equivalent
- length specifying the types of the groups. List elements are in
- [:class:`ost.ChemType.AMINOACIDS`,
- :class:`ost.ChemType.NUCLEOTIDES`]
- """
- pep_groups = _GroupSequences(pep_seqs, pep_seqid_thr, min_pep_length, aligner,
- mol.ChemType.AMINOACIDS)
- nuc_groups = _GroupSequences(nuc_seqs, nuc_seqid_thr, min_nuc_length, aligner,
- mol.ChemType.NUCLEOTIDES)
- group_types = [mol.ChemType.AMINOACIDS] * len(pep_groups)
- group_types += [mol.ChemType.NUCLEOTIDES] * len(nuc_groups)
- groups = pep_groups
- groups.extend(nuc_groups)
- return (groups, group_types)
-
-def _GroupSequences(seqs, seqid_thr, min_length, aligner, chem_type):
- """Get list of alignments representing groups of equivalent sequences
-
- :param seqid_thr: Threshold used to decide when two chains are identical.
- :type seqid_thr: :class:`float`
- :param gap_thr: Additional threshold to avoid gappy alignments with high
- seqid. Number of aligned columns must be at least this
- number.
- :type gap_thr: :class:`int`
- :param aligner: Helper class to generate pairwise alignments
- :type aligner: :class:`_Aligner`
- :param chem_type: ChemType of seqs which is passed to *aligner*, must be in
- [:class:`ost.mol.ChemType.AMINOACIDS`,
- :class:`ost.mol.ChemType.NUCLEOTIDES`]
- :type chem_type: :class:`ost.mol.ChemType`
- :returns: A list of alignments, one alignment for each group
- with longest sequence (reference) as first sequence.
- :rtype: :class:`ost.seq.AlignmentList`
- """
- groups = list()
- for s_idx in range(len(seqs)):
- matching_group = None
- for g_idx in range(len(groups)):
- for g_s_idx in range(len(groups[g_idx])):
- aln = aligner.Align(seqs[s_idx], seqs[groups[g_idx][g_s_idx]],
- chem_type)
- sid, n_aligned = _GetAlnPropsOne(aln)
- if sid >= seqid_thr and n_aligned >= min_length:
- matching_group = g_idx
- break
- if matching_group is not None:
- break
-
- if matching_group is None:
- groups.append([s_idx])
- else:
- groups[matching_group].append(s_idx)
-
- # sort based on sequence length
- sorted_groups = list()
- for g in groups:
- if len(g) > 1:
- tmp = sorted([[len(seqs[i]), i] for i in g], reverse=True)
- sorted_groups.append([x[1] for x in tmp])
- else:
- sorted_groups.append(g)
-
- # translate from indices back to sequences and directly generate alignments
- # of the groups with the longest (first) sequence as reference
- aln_list = seq.AlignmentList()
- for g in sorted_groups:
- if len(g) == 1:
- # aln with one single sequence
- aln_list.append(seq.CreateAlignment(seqs[g[0]]))
- else:
- # obtain pairwise aln of first sequence (reference) to all others
- alns = seq.AlignmentList()
- i = g[0]
- for j in g[1:]:
- alns.append(aligner.Align(seqs[i], seqs[j], chem_type))
- # and merge
- aln_list.append(seq.alg.MergePairwiseAlignments(alns, seqs[i]))
-
- # transfer attached views
- seq_dict = {s.GetName(): s for s in seqs}
- for aln_idx in range(len(aln_list)):
- for aln_s_idx in range(aln_list[aln_idx].GetCount()):
- s_name = aln_list[aln_idx].GetSequence(aln_s_idx).GetName()
- s = seq_dict[s_name]
- aln_list[aln_idx].AttachView(aln_s_idx, s.GetAttachedView())
-
- return aln_list
-
-def _MapSequence(ref_seqs, ref_types, s, s_type, aligner,
- seq_id_thr=0.0, min_aln_length=0):
- """Tries top map *s* onto any of the sequences in *ref_seqs*
-
- Computes alignments of *s* to each of the reference sequences of equal type
- and sorts them by seqid*fraction_covered (seqid: sequence identity of
- aligned columns in alignment, fraction_covered: Fraction of non-gap
- characters in reference sequence that are covered by non-gap characters in
- *s*). Best scoring mapping is returned. Optionally, *seq_id*/
- *min_aln_length* thresholds can be enabled to avoid non-sensical mappings.
- However, *min_aln_length* only has an effect if *seq_id_thr* > 0!!!
-
- :param ref_seqs: Reference sequences
- :type ref_seqs: :class:`ost.seq.SequenceList`
- :param ref_types: Types of reference sequences, e.g.
- ost.mol.ChemType.AminoAcids
- :type ref_types: :class:`list` of :class:`ost.mol.ChemType`
- :param s: Sequence to map
- :type s: :class:`ost.seq.SequenceHandle`
- :param s_type: Type of *s*, only try mapping to sequences in *ref_seqs*
- with equal type as defined in *ref_types*
- :param aligner: Helper class to generate pairwise alignments
- :type aligner: :class:`_Aligner`
- :param seq_id_thr: Minimum sequence identity to be considered as match
- :type seq_id_thr: :class:`float`
- :param min_aln_length: Minimum number of aligned columns to be considered
- as match. Only has an effect if *seq_id_thr* > 0!
- :type min_aln_length: :class:`int`
- :returns: Tuple with two elements. 1) index of sequence in *ref_seqs* to
- which *s* can be mapped 2) Pairwise sequence alignment with
- sequence from *ref_seqs* as first sequence. Both elements are
- None if no mapping can be found or if thresholds are not
- fulfilled for any alignment.
- :raises: :class:`RuntimeError` if mapping is ambiguous, i.e. *s*
- successfully maps to more than one sequence in *ref_seqs*
- """
- scored_alns = list()
- for ref_idx, ref_seq in enumerate(ref_seqs):
- if ref_types[ref_idx] == s_type:
- aln = aligner.Align(ref_seq, s, s_type)
- seqid, n_tot, n_aligned = _GetAlnPropsTwo(aln)
- if seq_id_thr > 0:
- if seqid >= seq_id_thr and n_aligned >= min_aln_length:
- fraction_covered = float(n_aligned)/n_tot
- score = seqid * fraction_covered
- scored_alns.append((score, ref_idx, aln))
- else:
- fraction_covered = float(n_aligned)/n_tot
- score = seqid * fraction_covered
- scored_alns.append((score, ref_idx, aln))
-
- if len(scored_alns) == 0:
- return (None, None) # no mapping possible...
-
- scored_alns = sorted(scored_alns, key=lambda x: x[0], reverse=True)
- return (scored_alns[0][1], scored_alns[0][2])
-
def _GetRefMdlAlns(ref_chem_groups, ref_chem_group_msas, mdl_chem_groups,
mdl_chem_group_alns, pairs=None):
""" Get all possible ref/mdl chain alignments given chem group mapping
diff --git a/modules/mol/alg/pymod/contact_score.py b/modules/mol/alg/pymod/contact_score.py
index 845f620e0..0af2d467a 100644
--- a/modules/mol/alg/pymod/contact_score.py
+++ b/modules/mol/alg/pymod/contact_score.py
@@ -384,7 +384,7 @@ class ContactScorerResultIPS:
:type: :class:`int`
"""
- return self._n_trg_contacts
+ return self._n_trg_int_res
@property
def n_mdl_int_res(self):
@@ -738,19 +738,29 @@ class ContactScorer:
trg_int_r = (trg_ch2, trg_ch1)
mdl_int_r = (mdl_ch2, mdl_ch1)
+ trg_contacts = None
if trg_int in self.cent1.contacts:
- n_trg = len(self.cent1.contacts[trg_int])
+ trg_contacts = self.cent1.contacts[trg_int]
elif trg_int_r in self.cent1.contacts:
- n_trg = len(self.cent1.contacts[trg_int_r])
- else:
+ trg_contacts = self.cent1.contacts[trg_int_r]
+
+ if trg_contacts is None:
n_trg = 0
+ else:
+ n_trg = len(set([x[0] for x in trg_contacts]))
+ n_trg += len(set([x[1] for x in trg_contacts]))
+ mdl_contacts = None
if mdl_int in self.cent2.contacts:
- n_mdl = len(self.cent2.contacts[mdl_int])
+ mdl_contacts = self.cent2.contacts[mdl_int]
elif mdl_int_r in self.cent2.contacts:
- n_mdl = len(self.cent2.contacts[mdl_int_r])
- else:
+ mdl_contacts = self.cent2.contacts[mdl_int_r]
+
+ if mdl_contacts is None:
n_mdl = 0
+ else:
+ n_mdl = len(set([x[0] for x in mdl_contacts]))
+ n_mdl += len(set([x[1] for x in mdl_contacts]))
_, _, n_union, n_intersection = self._MappedInterfaceScores(trg_int, mdl_int)
return ContactScorerResultIPS(n_trg, n_mdl, n_union, n_intersection)
diff --git a/modules/mol/alg/pymod/scoring.py b/modules/mol/alg/pymod/scoring.py
index cea7db318..f1d164a8e 100644
--- a/modules/mol/alg/pymod/scoring.py
+++ b/modules/mol/alg/pymod/scoring.py
@@ -22,6 +22,48 @@ from ost.bindings import cadscore
from ost.bindings import tmtools
import numpy as np
+def _GetAlignedResidues(aln, s1_ent, s2_ent):
+ """ Yields aligned residues
+
+ :param aln: The alignment with 2 sequences defining a residue-by-residue
+ relationship.
+ :type aln: :class:`ost.seq.AlignmentHandle`
+ :param s1_ent: Structure representing first sequence in *aln*.
+ One chain must be named after the first sequence and the
+ number of residues must match the number of non-gap
+ characters.
+ :type s1_ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
+ :param s2_ent: Same for second sequence in *aln*.
+ :type s2_ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
+ """
+ s1_ch = s1_ent.FindChain(aln.GetSequence(0).GetName())
+ s2_ch = s2_ent.FindChain(aln.GetSequence(1).GetName())
+
+ if not s1_ch.IsValid():
+ raise RuntimeError("s1_ent lacks required chain in _GetAlignedResidues")
+
+ if not s2_ch.IsValid():
+ raise RuntimeError("s2_ent lacks required chain in _GetAlignedResidues")
+
+ if len(aln.GetSequence(0).GetGaplessString()) != s1_ch.GetResidueCount():
+ raise RuntimeError("aln/chain mismatch in _GetAlignedResidues")
+ if len(aln.GetSequence(1).GetGaplessString()) != s2_ch.GetResidueCount():
+ raise RuntimeError("aln/chain mismatch in _GetAlignedResidues")
+
+ s1_res = s1_ch.residues
+ s2_res = s2_ch.residues
+
+ s1_res_idx = 0
+ s2_res_idx = 0
+
+ for col in aln:
+ if col[0] != '-' and col[1] != '-':
+ yield (s1_res[s1_res_idx], s2_res[s2_res_idx])
+ if col[0] != '-':
+ s1_res_idx += 1
+ if col[1] != '-':
+ s2_res_idx += 1
+
class lDDTBSScorer:
"""Scorer specific for a reference/model pair
@@ -237,6 +279,10 @@ class Scorer:
self._target_orig = target
self._model_orig = model
+ # lazily computed versions of target_orig and model_orig
+ self._pepnuc_target = None
+ self._pepnuc_model = None
+
if isinstance(self._model_orig, mol.EntityView):
self._model = mol.CreateEntityFromView(self._model_orig, False)
else:
@@ -493,6 +539,19 @@ class Scorer:
"""
return self._model_orig
+ @property
+ def pepnuc_model(self):
+ """ A selection of :attr:`~model_orig`
+
+ Only contains peptide and nucleotide residues
+
+ :type: :class:`ost.mol.EntityView`
+ """
+ if self._pepnuc_model is None:
+ query = "peptide=true or nucleotide=true"
+ self._pepnuc_model = self.model_orig.Select(query)
+ return self._pepnuc_model
+
@property
def target(self):
""" Target with Molck cleanup
@@ -509,6 +568,19 @@ class Scorer:
"""
return self._target_orig
+ @property
+ def pepnuc_target(self):
+ """ A selection of :attr:`~target_orig`
+
+ Only contains peptide and nucleotide residues
+
+ :type: :class:`ost.mol.EntityView`
+ """
+ if self._pepnuc_target is None:
+ query = "peptide=true or nucleotide=true"
+ self._pepnuc_target = self.target_orig.Select(query)
+ return self._pepnuc_target
+
@property
def aln(self):
""" Alignments of :attr:`~model`/:attr:`~target` chains
@@ -2044,14 +2116,12 @@ class Scorer:
cname = ch.GetName()
s = ''.join([r.one_letter_code for r in ch.residues])
s = seq.CreateSequence(ch.GetName(), s)
- s.AttachView(target.Select(f"cname={mol.QueryQuoteName(cname)}"))
trg_seqs[ch.GetName()] = s
mdl_seqs = dict()
for ch in model.chains:
cname = ch.GetName()
s = ''.join([r.one_letter_code for r in ch.residues])
s = seq.CreateSequence(cname, s)
- s.AttachView(model.Select(f"cname={mol.QueryQuoteName(cname)}"))
mdl_seqs[ch.GetName()] = s
alns = list()
@@ -2068,25 +2138,73 @@ class Scorer:
else:
raise RuntimeError("Chain name inconsistency... ask "
"Gabriel")
- alns.append(self.chain_mapper.Align(trg_seqs[trg_ch],
+ if self.resnum_alignments:
+ aln = self.chain_mapper.ResNumAlign(trg_seqs[trg_ch],
+ mdl_seqs[mdl_ch],
+ target, model)
+ else:
+ aln = self.chain_mapper.NWAlign(trg_seqs[trg_ch],
mdl_seqs[mdl_ch],
- stype))
- alns[-1].AttachView(0, trg_seqs[trg_ch].GetAttachedView())
- alns[-1].AttachView(1, mdl_seqs[mdl_ch].GetAttachedView())
- return alns
+ stype)
- def _compute_pepnuc_aln(self):
- query = "peptide=true or nucleotide=true"
- pep_nuc_target = self.target_orig.Select(query)
- pep_nuc_model = self.model_orig.Select(query)
- self._pepnuc_aln = self._aln_helper(pep_nuc_target, pep_nuc_model)
+ alns.append(aln)
+ return alns
def _compute_aln(self):
self._aln = self._aln_helper(self.target, self.model)
def _compute_stereochecked_aln(self):
- self._stereochecked_aln = self._aln_helper(self.stereochecked_target,
- self.stereochecked_model)
+ # lets not redo the alignment and derive it from self.aln
+ alns = list()
+ for a in self.aln:
+ trg_s = a.GetSequence(0)
+ mdl_s = a.GetSequence(1)
+ trg_ch = self.target.FindChain(trg_s.name)
+ mdl_ch = self.model.FindChain(mdl_s.name)
+
+ sc_trg_olc = ['-'] * len(trg_s)
+ sc_mdl_olc = ['-'] * len(mdl_s)
+
+ sc_trg_ch = self.stereochecked_target.FindChain(trg_s.name)
+ if sc_trg_ch.IsValid():
+ # there is the theoretical possibility that the full chain
+ # has been removed in stereochemistry checks...
+ trg_residues = trg_ch.residues
+ res_idx = 0
+ for olc_idx, olc in enumerate(trg_s):
+ if olc != '-':
+ r = trg_residues[res_idx]
+ sc_r = sc_trg_ch.FindResidue(r.GetNumber())
+ if sc_r.IsValid():
+ sc_trg_olc[olc_idx] = sc_r.one_letter_code
+ res_idx += 1
+
+ sc_mdl_ch = self.stereochecked_model.FindChain(mdl_s.name)
+ if sc_mdl_ch.IsValid():
+ # there is the theoretical possibility that the full chain
+ # has been removed in stereochemistry checks...
+ mdl_residues = mdl_ch.residues
+ res_idx = 0
+ for olc_idx, olc in enumerate(mdl_s):
+ if olc != '-':
+ r = mdl_residues[res_idx]
+ sc_r = sc_mdl_ch.FindResidue(r.GetNumber())
+ if sc_r.IsValid():
+ sc_mdl_olc[olc_idx] = sc_r.one_letter_code
+ res_idx += 1
+
+ sc_trg_s = seq.CreateSequence(trg_s.name, ''.join(sc_trg_olc))
+ sc_mdl_s = seq.CreateSequence(mdl_s.name, ''.join(sc_mdl_olc))
+ new_a = seq.CreateAlignment()
+ new_a.AddSequence(sc_trg_s)
+ new_a.AddSequence(sc_mdl_s)
+ alns.append(new_a)
+
+ self._stereochecked_aln = alns
+
+ def _compute_pepnuc_aln(self):
+ self._pepnuc_aln = self._aln_helper(self.pepnuc_target,
+ self.pepnuc_model)
def _compute_lddt(self):
LogScript("Computing all-atom LDDT")
@@ -2167,10 +2285,12 @@ class Scorer:
local_lddt = dict()
aa_local_lddt = dict()
for r in self.model.residues:
+
cname = r.GetChain().GetName()
if cname not in local_lddt:
local_lddt[cname] = dict()
aa_local_lddt[cname] = dict()
+
rnum = r.GetNumber()
if rnum not in aa_local_lddt[cname]:
aa_local_lddt[cname][rnum] = dict()
@@ -2179,8 +2299,8 @@ class Scorer:
score = round(r.GetFloatProp("lddt"), 3)
local_lddt[cname][rnum] = score
- trg_r = None
- mdl_r = None
+ trg_r = None # represents stereochecked trg res
+ mdl_r = None # represents stereochecked mdl res
for a in r.atoms:
if a.HasProp("lddt"):
@@ -2195,18 +2315,24 @@ class Scorer:
# stereochecks but is there in stereochecked
# target => 0.0
if trg_r is None:
+ # let's first see if we find that target residue
+ # in the non-stereochecked target
+ tmp = None
if cname in flat_mapping:
- for col in alns[cname]:
- if col[0] != '-' and col[1] != '-':
- if col.GetResidue(1).number == r.number:
- trg_r = col.GetResidue(0)
- break
- if trg_r is not None:
- trg_cname = trg_r.GetChain().GetName()
- trg_rnum = trg_r.GetNumber()
- tmp = self.stereochecked_target.FindResidue(trg_cname,
- trg_rnum)
+ for x, y in _GetAlignedResidues(alns[cname],
+ self.target,
+ self.model):
+ if y.number == r.number:
+ tmp = x
+ break
+ if tmp is not None:
+ # we have it in the non-stereochecked target!
+ tmp_cname = tmp.GetChain().GetName()
+ tmp_rnum = tmp.GetNumber()
+ tmp = self.stereochecked_target.FindResidue(tmp_cname,
+ tmp_rnum)
if tmp.IsValid():
+ # And it's there in the stereochecked target too!
trg_r = tmp
if mdl_r is None:
@@ -2214,12 +2340,23 @@ class Scorer:
if tmp.IsValid():
mdl_r = tmp
- if trg_r is not None and not trg_r.FindAtom(a.GetName()).IsValid():
- # opt 1
+ if trg_r is None:
+ # opt 1 - the whole target residue is not there
+ # this is actually an impossibility, as we have
+ # a score for the full mdl residue set
+ aa_local_lddt[cname][rnum][a.GetName()] = None
+ elif trg_r is not None and not trg_r.FindAtom(a.GetName()).IsValid():
+ # opt 1 - the target residue is there but not the atom
aa_local_lddt[cname][rnum][a.GetName()] = None
+ elif trg_r is not None and trg_r.FindAtom(a.GetName()).IsValid() and \
+ mdl_r is None:
+ # opt 2 - trg atom is there but full model residue is removed
+ # this is actuall an impossibility, as we have
+ # a score for the full mdl residue set
+ aa_local_lddt[cname][rnum][a.GetName()] = 0.0
elif trg_r is not None and trg_r.FindAtom(a.GetName()).IsValid() and \
mdl_r is not None and not mdl_r.FindAtom(a.GetName()).IsValid():
- # opt 2
+ # opt 2 - trg atom is there but model atom is removed
aa_local_lddt[cname][rnum][a.GetName()] = 0.0
else:
# unknown issue
@@ -2237,20 +2374,25 @@ class Scorer:
# opt 1: removed by stereochecks => assign 0.0
# opt 2: removed by stereochecks AND not covered by ref
# => assign None
+
# fetch trg residue from non-stereochecked aln
trg_r = None
if cname in flat_mapping:
- for col in alns[cname]:
- if col[0] != '-' and col[1] != '-':
- if col.GetResidue(1).number == r.number:
- trg_r = col.GetResidue(0)
- break
- if trg_r is not None:
- trg_cname = trg_r.GetChain().GetName()
- trg_rnum = trg_r.GetNumber()
- tmp = self.stereochecked_target.FindResidue(trg_cname,
- trg_rnum)
+ tmp = None
+ for x, y in _GetAlignedResidues(alns[cname],
+ self.target,
+ self.model):
+ if y.number == r.number:
+ tmp = x
+ break
+ if tmp is not None:
+ # we have it in the non-stereochecked target!
+ tmp_cname = tmp.GetChain().GetName()
+ tmp_rnum = tmp.GetNumber()
+ tmp = self.stereochecked_target.FindResidue(tmp_cname,
+ tmp_rnum)
if tmp.IsValid():
+ # And it's there in the stereochecked target too!
trg_r = tmp
if trg_r is None:
@@ -2617,30 +2759,30 @@ class Scorer:
for trg_ch, mdl_ch in self.mapping.GetFlatMapping().items():
processed_trg_chains.add(trg_ch)
aln = self.mapping.alns[(trg_ch, mdl_ch)]
- for col in aln:
- if col[0] != '-' and col[1] != '-':
- trg_res = col.GetResidue(0)
- mdl_res = col.GetResidue(1)
- trg_at = trg_res.FindAtom("CA")
- mdl_at = mdl_res.FindAtom("CA")
- if not trg_at.IsValid():
- trg_at = trg_res.FindAtom("C3'")
- if not mdl_at.IsValid():
- mdl_at = mdl_res.FindAtom("C3'")
- self._mapped_target_pos.append(trg_at.GetPos())
- self._mapped_model_pos.append(mdl_at.GetPos())
- elif col[0] != '-':
- self._n_target_not_mapped += 1
+ n_mapped = 0
+ for trg_res, mdl_res in _GetAlignedResidues(aln,
+ self.mapping.target,
+ self.mapping.model):
+ trg_at = trg_res.FindAtom("CA")
+ mdl_at = mdl_res.FindAtom("CA")
+ if not trg_at.IsValid():
+ trg_at = trg_res.FindAtom("C3'")
+ if not mdl_at.IsValid():
+ mdl_at = mdl_res.FindAtom("C3'")
+ self._mapped_target_pos.append(trg_at.GetPos())
+ self._mapped_model_pos.append(mdl_at.GetPos())
+ n_mapped += 1
+ self._n_target_not_mapped += (len(aln.GetSequence(0).GetGaplessString())-n_mapped)
# count number of trg residues from non-mapped chains
for ch in self.mapping.target.chains:
if ch.GetName() not in processed_trg_chains:
- self._n_target_not_mapped += len(ch.residues)
+ self._n_target_not_mapped += ch.GetResidueCount()
def _extract_mapped_pos_full_bb(self):
self._mapped_target_pos_full_bb = geom.Vec3List()
self._mapped_model_pos_full_bb = geom.Vec3List()
exp_pep_atoms = ["N", "CA", "C"]
- exp_nuc_atoms = ["\"O5'\"", "\"C5'\"", "\"C4'\"", "\"C3'\"", "\"O3'\""]
+ exp_nuc_atoms = ["O5'", "C5'", "C4'", "C3'", "O3'"]
trg_pep_chains = [s.GetName() for s in self.chain_mapper.polypep_seqs]
trg_nuc_chains = [s.GetName() for s in self.chain_mapper.polynuc_seqs]
for trg_ch, mdl_ch in self.mapping.GetFlatMapping().items():
@@ -2653,19 +2795,18 @@ class Scorer:
else:
# this should be guaranteed by the chain mapper
raise RuntimeError("Unexpected error - contact OST developer")
- for col in aln:
- if col[0] != '-' and col[1] != '-':
- trg_res = col.GetResidue(0)
- mdl_res = col.GetResidue(1)
- for aname in exp_atoms:
- trg_at = trg_res.FindAtom(aname)
- mdl_at = mdl_res.FindAtom(aname)
- if not (trg_at.IsValid() and mdl_at.IsValid()):
- # this should be guaranteed by the chain mapper
- raise RuntimeError("Unexpected error - contact OST "
- "developer")
- self._mapped_target_pos_full_bb.append(trg_at.GetPos())
- self._mapped_model_pos_full_bb.append(mdl_at.GetPos())
+ for trg_res, mdl_res in _GetAlignedResidues(aln,
+ self.mapping.target,
+ self.mapping.model):
+ for aname in exp_atoms:
+ trg_at = trg_res.FindAtom(aname)
+ mdl_at = mdl_res.FindAtom(aname)
+ if not (trg_at.IsValid() and mdl_at.IsValid()):
+ # this should be guaranteed by the chain mapper
+ raise RuntimeError("Unexpected error - contact OST "
+ "developer")
+ self._mapped_target_pos_full_bb.append(trg_at.GetPos())
+ self._mapped_model_pos_full_bb.append(mdl_at.GetPos())
def _extract_rigid_mapped_pos(self):
@@ -2676,34 +2817,35 @@ class Scorer:
for trg_ch, mdl_ch in self.rigid_mapping.GetFlatMapping().items():
processed_trg_chains.add(trg_ch)
aln = self.rigid_mapping.alns[(trg_ch, mdl_ch)]
- for col in aln:
- if col[0] != '-' and col[1] != '-':
- trg_res = col.GetResidue(0)
- mdl_res = col.GetResidue(1)
- trg_at = trg_res.FindAtom("CA")
- mdl_at = mdl_res.FindAtom("CA")
- if not trg_at.IsValid():
- trg_at = trg_res.FindAtom("C3'")
- if not mdl_at.IsValid():
- mdl_at = mdl_res.FindAtom("C3'")
- self._rigid_mapped_target_pos.append(trg_at.GetPos())
- self._rigid_mapped_model_pos.append(mdl_at.GetPos())
- elif col[0] != '-':
- self._rigid_n_target_not_mapped += 1
+ n_mapped = 0
+ for trg_res, mdl_res in _GetAlignedResidues(aln,
+ self.rigid_mapping.target,
+ self.rigid_mapping.model):
+ trg_at = trg_res.FindAtom("CA")
+ mdl_at = mdl_res.FindAtom("CA")
+ if not trg_at.IsValid():
+ trg_at = trg_res.FindAtom("C3'")
+ if not mdl_at.IsValid():
+ mdl_at = mdl_res.FindAtom("C3'")
+ self._rigid_mapped_target_pos.append(trg_at.GetPos())
+ self._rigid_mapped_model_pos.append(mdl_at.GetPos())
+ n_mapped += 1
+
+ self._rigid_n_target_not_mapped += (len(aln.GetSequence(0).GetGaplessString())-n_mapped)
# count number of trg residues from non-mapped chains
for ch in self.rigid_mapping.target.chains:
if ch.GetName() not in processed_trg_chains:
- self._rigid_n_target_not_mapped += len(ch.residues)
+ self._rigid_n_target_not_mapped += ch.GetResidueCount()
def _extract_rigid_mapped_pos_full_bb(self):
self._rigid_mapped_target_pos_full_bb = geom.Vec3List()
self._rigid_mapped_model_pos_full_bb = geom.Vec3List()
exp_pep_atoms = ["N", "CA", "C"]
- exp_nuc_atoms = ["\"O5'\"", "\"C5'\"", "\"C4'\"", "\"C3'\"", "\"O3'\""]
+ exp_nuc_atoms = ["O5'", "C5'", "C4'", "C3'", "O3'"]
trg_pep_chains = [s.GetName() for s in self.chain_mapper.polypep_seqs]
trg_nuc_chains = [s.GetName() for s in self.chain_mapper.polynuc_seqs]
for trg_ch, mdl_ch in self.rigid_mapping.GetFlatMapping().items():
- aln = self.mapping.alns[(trg_ch, mdl_ch)]
+ aln = self.rigid_mapping.alns[(trg_ch, mdl_ch)]
trg_ch = aln.GetSequence(0).GetName()
if trg_ch in trg_pep_chains:
exp_atoms = exp_pep_atoms
@@ -2712,19 +2854,18 @@ class Scorer:
else:
# this should be guaranteed by the chain mapper
raise RuntimeError("Unexpected error - contact OST developer")
- for col in aln:
- if col[0] != '-' and col[1] != '-':
- trg_res = col.GetResidue(0)
- mdl_res = col.GetResidue(1)
- for aname in exp_atoms:
- trg_at = trg_res.FindAtom(aname)
- mdl_at = mdl_res.FindAtom(aname)
- if not (trg_at.IsValid() and mdl_at.IsValid()):
- # this should be guaranteed by the chain mapper
- raise RuntimeError("Unexpected error - contact OST "
- "developer")
- self._rigid_mapped_target_pos_full_bb.append(trg_at.GetPos())
- self._rigid_mapped_model_pos_full_bb.append(mdl_at.GetPos())
+ for trg_res, mdl_res in _GetAlignedResidues(aln,
+ self.rigid_mapping.target,
+ self.rigid_mapping.model):
+ for aname in exp_atoms:
+ trg_at = trg_res.FindAtom(aname)
+ mdl_at = mdl_res.FindAtom(aname)
+ if not (trg_at.IsValid() and mdl_at.IsValid()):
+ # this should be guaranteed by the chain mapper
+ raise RuntimeError("Unexpected error - contact OST "
+ "developer")
+ self._rigid_mapped_target_pos_full_bb.append(trg_at.GetPos())
+ self._rigid_mapped_model_pos_full_bb.append(mdl_at.GetPos())
def _compute_cad_score(self):
if not self.resnum_alignments:
@@ -2813,12 +2954,12 @@ class Scorer:
a, b, c, d = stereochemistry.StereoCheck(self.target, stereo_data = data,
stereo_link_data = l_data)
+
self._stereochecked_target = a
self._target_clashes = b
self._target_bad_bonds = c
self._target_bad_angles = d
-
def _get_interface_patches(self, mdl_ch, mdl_rnum):
""" Select interface patches representative for specified residue
@@ -2897,34 +3038,38 @@ class Scorer:
# transfer mdl residues to trg
flat_mapping = self.mapping.GetFlatMapping(mdl_as_key=True)
full_trg_coverage = True
- trg_patch_one = self.target.CreateEmptyView()
+ trg_patch_one = self.mapping.target.CreateEmptyView()
for r in mdl_patch_one.residues:
trg_r = None
mdl_cname = r.GetChain().GetName()
if mdl_cname in flat_mapping:
aln = self.mapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
- for col in aln:
- if col[0] != '-' and col[1] != '-':
- if col.GetResidue(1).GetNumber() == r.GetNumber():
- trg_r = col.GetResidue(0)
- break
+ for x,y in _GetAlignedResidues(aln,
+ self.mapping.target,
+ self.mapping.model):
+ if y.GetNumber() == r.GetNumber():
+ trg_r = x
+ break
+
if trg_r is not None:
trg_patch_one.AddResidue(trg_r.handle,
mol.ViewAddFlag.INCLUDE_ALL)
else:
full_trg_coverage = False
- trg_patch_two = self.target.CreateEmptyView()
+ trg_patch_two = self.mapping.target.CreateEmptyView()
for r in mdl_patch_two.residues:
trg_r = None
mdl_cname = r.GetChain().GetName()
if mdl_cname in flat_mapping:
aln = self.mapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
- for col in aln:
- if col[0] != '-' and col[1] != '-':
- if col.GetResidue(1).GetNumber() == r.GetNumber():
- trg_r = col.GetResidue(0)
- break
+ for x,y in _GetAlignedResidues(aln,
+ self.mapping.target,
+ self.mapping.model):
+ if y.GetNumber() == r.GetNumber():
+ trg_r = x
+ break
+
if trg_r is not None:
trg_patch_two.AddResidue(trg_r.handle,
mol.ViewAddFlag.INCLUDE_ALL)
diff --git a/modules/mol/alg/tests/CMakeLists.txt b/modules/mol/alg/tests/CMakeLists.txt
index e21f07468..1c030fec3 100644
--- a/modules/mol/alg/tests/CMakeLists.txt
+++ b/modules/mol/alg/tests/CMakeLists.txt
@@ -21,7 +21,8 @@ if (COMPOUND_LIB)
list(APPEND OST_MOL_ALG_UNIT_TESTS test_qsscoring.py
test_nonstandard.py
test_chain_mapping.py
- test_ligand_scoring.py)
+ test_ligand_scoring.py
+ test_scoring.py)
endif()
ost_unittest(MODULE mol_alg SOURCES "${OST_MOL_ALG_UNIT_TESTS}" LINK ost_io)
diff --git a/modules/mol/alg/tests/test_chain_mapping.py b/modules/mol/alg/tests/test_chain_mapping.py
index d7ec18c40..3183ed7c8 100644
--- a/modules/mol/alg/tests/test_chain_mapping.py
+++ b/modules/mol/alg/tests/test_chain_mapping.py
@@ -63,15 +63,6 @@ class TestChainMapper(unittest.TestCase):
self.assertEqual(str(mapper.polynuc_seqs[0]), str(nuc_s_one))
self.assertEqual(str(mapper.polynuc_seqs[1]), str(nuc_s_two))
- for s in mapper.polypep_seqs:
- self.assertTrue(s.HasAttachedView())
- for s in mapper.polynuc_seqs:
- self.assertTrue(s.HasAttachedView())
- self.assertTrue(_CompareViews(mapper.polypep_seqs[0].GetAttachedView(), pep_view_one))
- self.assertTrue(_CompareViews(mapper.polypep_seqs[1].GetAttachedView(), pep_view_two))
- self.assertTrue(_CompareViews(mapper.polynuc_seqs[0].GetAttachedView(), nuc_view_one))
- self.assertTrue(_CompareViews(mapper.polynuc_seqs[1].GetAttachedView(), nuc_view_two))
-
# peptide sequences should be in the same group, the nucleotides not
self.assertEqual(len(mapper.chem_group_alignments), 3)
self.assertEqual(len(mapper.chem_groups), 3)
@@ -89,11 +80,6 @@ class TestChainMapper(unittest.TestCase):
self.assertEqual(str(mapper.chem_group_ref_seqs[0]), str(pep_s_one))
self.assertEqual(str(mapper.chem_group_ref_seqs[1]), str(nuc_s_one))
self.assertEqual(str(mapper.chem_group_ref_seqs[2]), str(nuc_s_two))
- for s in mapper.chem_group_ref_seqs:
- self.assertTrue(s.HasAttachedView())
- self.assertTrue(_CompareViews(mapper.chem_group_ref_seqs[0].GetAttachedView(), pep_view_one))
- self.assertTrue(_CompareViews(mapper.chem_group_ref_seqs[1].GetAttachedView(), nuc_view_one))
- self.assertTrue(_CompareViews(mapper.chem_group_ref_seqs[2].GetAttachedView(), nuc_view_two))
# check chem_group_alignments attribute
self.assertEqual(len(mapper.chem_group_alignments), 3)
@@ -108,14 +94,6 @@ class TestChainMapper(unittest.TestCase):
self.assertEqual(s0.GetGaplessString(), str(nuc_s_one))
s0 = mapper.chem_group_alignments[2].GetSequence(0)
self.assertEqual(s0.GetGaplessString(), str(nuc_s_two))
- self.assertTrue(mapper.chem_group_alignments[0].GetSequence(0).HasAttachedView())
- self.assertTrue(mapper.chem_group_alignments[0].GetSequence(1).HasAttachedView())
- self.assertTrue(mapper.chem_group_alignments[1].GetSequence(0).HasAttachedView())
- self.assertTrue(mapper.chem_group_alignments[2].GetSequence(0).HasAttachedView())
- self.assertTrue(_CompareViews(mapper.chem_group_alignments[0].GetSequence(0).GetAttachedView(), pep_view_one))
- self.assertTrue(_CompareViews(mapper.chem_group_alignments[0].GetSequence(1).GetAttachedView(), pep_view_two))
- self.assertTrue(_CompareViews(mapper.chem_group_alignments[1].GetSequence(0).GetAttachedView(), nuc_view_one))
- self.assertTrue(_CompareViews(mapper.chem_group_alignments[2].GetSequence(0).GetAttachedView(), nuc_view_two))
# ensure that error is triggered if there are insertion codes
# and resnum_alignments are enabled
@@ -218,19 +196,6 @@ class TestChainMapper(unittest.TestCase):
self.assertEqual(alns[2][0].GetSequence(0).GetGaplessString(), str(ref_nuc_s_two))
self.assertEqual(alns[2][0].GetSequence(1).GetGaplessString(), str(mdl_nuc_s_one))
- self.assertTrue(alns[0][0].GetSequence(0).HasAttachedView())
- self.assertTrue(alns[0][0].GetSequence(1).HasAttachedView())
- self.assertTrue(alns[0][1].GetSequence(0).HasAttachedView())
- self.assertTrue(alns[0][1].GetSequence(1).HasAttachedView())
- self.assertTrue(alns[2][0].GetSequence(0).HasAttachedView())
- self.assertTrue(alns[2][0].GetSequence(1).HasAttachedView())
- self.assertTrue(_CompareViews(alns[0][0].GetSequence(0).GetAttachedView(),ref_pep_view_one))
- self.assertTrue(_CompareViews(alns[0][0].GetSequence(1).GetAttachedView(),mdl_pep_view_one))
- self.assertTrue(_CompareViews(alns[0][1].GetSequence(0).GetAttachedView(),ref_pep_view_one))
- self.assertTrue(_CompareViews(alns[0][1].GetSequence(1).GetAttachedView(),mdl_pep_view_two))
- self.assertTrue(_CompareViews(alns[2][0].GetSequence(0).GetAttachedView(),ref_nuc_view_two))
- self.assertTrue(_CompareViews(alns[2][0].GetSequence(1).GetAttachedView(),mdl_nuc_view_one))
-
# test for unmapped mdl chains, i.e. chains in the mdl that are not
# present in ref
mdl = _LoadFile("mdl_different_chain_mdl.pdb")
@@ -355,7 +320,7 @@ class TestChainMapper(unittest.TestCase):
self.assertTrue(mdl_s is not None)
self.assertEqual(ref_s_type, mdl_s_type)
- aln = mapper.Align(ref_s, mdl_s, ref_s_type)
+ aln = mapper.NWAlign(ref_s, mdl_s, ref_s_type)
self.assertEqual(ref_aln.GetSequence(0).GetName(),
aln.GetSequence(0).GetName())
self.assertEqual(ref_aln.GetSequence(1).GetName(),
@@ -415,6 +380,18 @@ class TestChainMapper(unittest.TestCase):
for ref_ch, mdl_ch in repr_mapping.items():
self.assertEqual(mdl_ch, flat_mapping[ref_ch])
+ def test_resnum_aln_mapping(self):
+ ref = _LoadFile("3l1p.1.pdb")
+ mdl = _LoadFile("3l1p.1_model.pdb")
+ mapper = ChainMapper(ref, resnum_alignments=True)
+
+ # Again, no in depth testing. Its only about checking if residue number
+ # alignments run through
+
+ # lDDT based chain mappings
+ naive_lddt_res = mapper.GetlDDTMapping(mdl, strategy="naive")
+ self.assertEqual(naive_lddt_res.mapping, [['X', 'Y'],[None],['Z']])
+
def test_misc(self):
# check for triggered error when no chain fulfills length threshold
diff --git a/modules/mol/alg/tests/test_scoring.py b/modules/mol/alg/tests/test_scoring.py
new file mode 100644
index 000000000..162808afc
--- /dev/null
+++ b/modules/mol/alg/tests/test_scoring.py
@@ -0,0 +1,242 @@
+import unittest, os, sys
+import ost
+from ost import io, mol, geom
+# check if we can import: fails if numpy or scipy not available
+try:
+ from ost.mol.alg.scoring import *
+except ImportError:
+ print("Failed to import scoring.py. Happens when numpy or scipy "\
+ "missing. Ignoring test_scoring.py tests.")
+ sys.exit(0)
+
+def _LoadFile(file_name):
+ """Helper to avoid repeating input path over and over."""
+ return io.LoadPDB(os.path.join('testfiles', file_name))
+
+class TestScorer(unittest.TestCase):
+
+ # compare to hardcoded values - no in depth testing
+ # this should be sufficient to flag issues when changes are introduced
+
+ def test_scorer_lddt(self):
+
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_ref.pdb")
+
+ scorer = Scorer(mdl, trg)
+
+ # check global lDDT values
+ self.assertAlmostEqual(scorer.lddt, 0.539, 3)
+ self.assertAlmostEqual(scorer.bb_lddt, 0.622, 3)
+ self.assertAlmostEqual(scorer.ilddt, 0.282, 3)
+
+ for ch in scorer.model.chains:
+ self.assertTrue(ch.name in scorer.local_lddt)
+ self.assertEqual(len(scorer.local_lddt[ch.name]), ch.GetResidueCount())
+ self.assertTrue(ch.name in scorer.bb_local_lddt)
+ self.assertEqual(len(scorer.bb_local_lddt[ch.name]), ch.GetResidueCount())
+
+ for ch in scorer.model.chains:
+ self.assertTrue(ch.name in scorer.aa_local_lddt)
+ self.assertEqual(len(scorer.aa_local_lddt[ch.name]), ch.GetResidueCount())
+
+ # check some random per-residue/per-atom scores
+ self.assertEqual(scorer.local_lddt["B"][mol.ResNum(42)], 0.659, 3)
+ self.assertEqual(scorer.local_lddt["A"][mol.ResNum(142)], 0.849, 3)
+ self.assertEqual(scorer.bb_local_lddt["B"][mol.ResNum(42)], 0.782, 3)
+ self.assertEqual(scorer.bb_local_lddt["A"][mol.ResNum(142)], 0.910, 3)
+ self.assertEqual(scorer.aa_local_lddt["B"][mol.ResNum(42)]["CA"], 0.718, 3)
+ self.assertEqual(scorer.aa_local_lddt["A"][mol.ResNum(142)]["CB"], 0.837, 3)
+
+ # test stereochemistry checks related behaviour
+
+ # stereochemistry issue in mdl sidechain
+ bad_mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ ed = bad_mdl.EditXCS()
+ at = bad_mdl.FindResidue("B", mol.ResNum(42)).FindAtom("CD")
+ pos = at.GetPos()
+ new_pos = geom.Vec3(pos[0], pos[1], pos[2] + 1.0)
+ ed.SetAtomPos(at, new_pos)
+ scorer = Scorer(bad_mdl, trg)
+ # original score without stereochemistry issues: 0.659
+ # now it should be much worse but above zero
+ penalized_score = scorer.local_lddt["B"][mol.ResNum(42)]
+ self.assertTrue(penalized_score < 0.5 and penalized_score > 0.1)
+
+ # let's make it really bad, i.e. involve a backbone atom
+ at = bad_mdl.FindResidue("B", mol.ResNum(42)).FindAtom("CA")
+ pos = at.GetPos()
+ new_pos = geom.Vec3(pos[0], pos[1], pos[2] + 1.0)
+ ed.SetAtomPos(at, new_pos)
+ scorer = Scorer(bad_mdl, trg)
+ # original score without stereochemistry issues: 0.659
+ # now it should be 0.0
+ penalized_score = scorer.local_lddt["B"][mol.ResNum(42)]
+ self.assertEqual(penalized_score, 0.0)
+
+ # let's fiddle around in trg stereochemistry
+ bad_trg = _LoadFile("1eud_ref.pdb")
+ ed = bad_trg.EditXCS()
+ # thats the atom that should map to B.42 in mdl
+ at = bad_trg.FindResidue("B", mol.ResNum(5)).FindAtom("CD")
+ pos = at.GetPos()
+ new_pos = geom.Vec3(pos[0], pos[1], pos[2] + 1.0)
+ ed.SetAtomPos(at, new_pos)
+ scorer = Scorer(mdl, bad_trg)
+ # there is no reference info anymore on the whole sidechain
+ # The scores of the sidechain atoms should thus be None
+ self.assertTrue(scorer.aa_local_lddt["B"][mol.ResNum(42)]["CD"] is None)
+ # but not the sidechain atom scores from backbone atoms!
+ self.assertFalse(scorer.aa_local_lddt["B"][mol.ResNum(42)]["CA"] is None)
+ # also the full per-residue score is still a valid number
+ self.assertFalse(scorer.local_lddt["B"][mol.ResNum(42)] is None)
+
+ bad_trg = _LoadFile("1eud_ref.pdb")
+ ed = bad_trg.EditXCS()
+ at = bad_trg.FindResidue("B", mol.ResNum(5)).FindAtom("CA")
+ pos = at.GetPos()
+ new_pos = geom.Vec3(pos[0], pos[1], pos[2] + 1.0)
+
+ ed.SetAtomPos(at, new_pos)
+ scorer = Scorer(mdl, bad_trg)
+
+ # all scores should be None now for this residue
+ self.assertTrue(scorer.aa_local_lddt["B"][mol.ResNum(42)]["CD"] is None)
+ self.assertTrue(scorer.aa_local_lddt["B"][mol.ResNum(42)]["CA"] is None)
+ self.assertTrue(scorer.local_lddt["B"][mol.ResNum(42)] is None)
+
+
+ def test_scorer_qsscore(self):
+
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_ref.pdb")
+
+ scorer = Scorer(mdl, trg)
+
+ # check qs-score related values
+ self.assertAlmostEqual(scorer.qs_global, 0.321, 3)
+ self.assertAlmostEqual(scorer.qs_best, 0.932, 3)
+ self.assertEqual(len(scorer.qs_interfaces), 1)
+ self.assertEqual(scorer.qs_interfaces[0], ("A", "B", "A", "B"))
+
+ # should be equal global scores since we're only dealing with
+ # a single interface
+ self.assertAlmostEqual(scorer.per_interface_qs_global[0], 0.321, 3)
+ self.assertAlmostEqual(scorer.per_interface_qs_best[0], 0.932, 3)
+
+ def test_scorer_rigid_scores(self):
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_ref.pdb")
+ scorer = Scorer(mdl, trg)
+ self.assertAlmostEqual(scorer.gdtts, 0.616, 3)
+ self.assertAlmostEqual(scorer.gdtha, 0.473, 3)
+ self.assertAlmostEqual(scorer.rmsd, 2.944, 3)
+
+ def test_scorer_contacts(self):
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_ref.pdb")
+ scorer = Scorer(mdl, trg)
+
+ self.assertEqual(len(scorer.model_contacts), 48)
+ self.assertEqual(len(scorer.native_contacts), 140)
+
+ self.assertAlmostEqual(scorer.ics, 0.415, 3)
+ self.assertAlmostEqual(scorer.ics_precision, 0.812, 3)
+ self.assertAlmostEqual(scorer.ics_recall, 0.279, 3)
+ self.assertAlmostEqual(scorer.ips, 0.342, 3)
+ self.assertAlmostEqual(scorer.ips_precision, 0.891, 3)
+ self.assertAlmostEqual(scorer.ips_recall, 0.357, 3)
+
+
+ # per interface scores should be equal since we're only dealing with one
+ # interface
+ self.assertEqual(len(scorer.per_interface_ics), 1)
+ self.assertAlmostEqual(scorer.per_interface_ics[0], 0.415, 3)
+ self.assertAlmostEqual(scorer.per_interface_ics_precision[0], 0.812, 3)
+ self.assertAlmostEqual(scorer.per_interface_ics_recall[0], 0.279, 3)
+
+ self.assertEqual(len(scorer.per_interface_ips), 1)
+ self.assertAlmostEqual(scorer.per_interface_ips[0], 0.342, 3)
+ self.assertAlmostEqual(scorer.per_interface_ips_precision[0], 0.891, 3)
+ self.assertAlmostEqual(scorer.per_interface_ips_recall[0], 0.357, 3)
+
+ def test_scorer_dockq(self):
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_ref.pdb")
+ scorer = Scorer(mdl, trg)
+ self.assertEqual(scorer.dockq_interfaces, [("A", "B", "A", "B")])
+ self.assertAlmostEqual(scorer.dockq_scores[0], 0.559, 3)
+ self.assertAlmostEqual(scorer.fnat[0], 0.279, 3)
+ self.assertAlmostEqual(scorer.fnonnat[0], 0.188, 2)
+ self.assertAlmostEqual(scorer.irmsd[0], 0.988, 3)
+ self.assertAlmostEqual(scorer.lrmsd[0], 5.533, 3)
+ self.assertEqual(scorer.nnat[0], 140)
+ self.assertEqual(scorer.nmdl[0], 48)
+
+ # ave and wave values are the same as scorer.dockq_scores[0]
+ self.assertAlmostEqual(scorer.dockq_wave, scorer.dockq_scores[0], 7)
+ self.assertAlmostEqual(scorer.dockq_ave, scorer.dockq_scores[0], 7)
+
+ def test_scorer_patch_scores(self):
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_ref.pdb")
+ scorer = Scorer(mdl, trg)
+
+ patch_qs = scorer.patch_qs
+ patch_dockq = scorer.patch_dockq
+
+ # check some random values
+ self.assertAlmostEqual(patch_qs["B"][5], 0.925, 3)
+ self.assertAlmostEqual(patch_qs["B"][17], 0.966, 3)
+ self.assertAlmostEqual(patch_qs["A"][4], 0.973, 3)
+
+ self.assertAlmostEqual(patch_dockq["B"][5], 0.858, 3)
+ self.assertAlmostEqual(patch_dockq["B"][17], 0.965, 3)
+ self.assertAlmostEqual(patch_dockq["A"][4], 0.973, 3)
+
+ def test_scorer_trimmed_contacts(self):
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_mdl_partial-dimer.pdb")
+
+ scorer = Scorer(mdl, trg)
+
+ # mdl and trg are the same
+ self.assertEqual(scorer.ics, 1.0)
+ self.assertEqual(scorer.ips, 1.0)
+
+ # let's cut a critical interface loop in the target
+ trg = trg.Select("cname=A or (cname=B and rnum!=157:168)")
+ scorer = Scorer(mdl, trg)
+
+ # mdl is now longer than trg which lowers ICS/IPS
+ self.assertTrue(scorer.ics < 0.75)
+ self.assertTrue(scorer.ips < 0.75)
+
+ # but if we use the trimmed versions, it should go up to 1.0
+ # again
+ self.assertEqual(scorer.ics_trimmed, 1.0)
+ self.assertEqual(scorer.ips_trimmed, 1.0)
+
+ # lets see if the trimmed model has the right
+ # residues missing
+ for r in scorer.model.residues:
+ cname = r.GetChain().GetName()
+ rnum = r.GetNumber()
+ trimmed_r = scorer.trimmed_model.FindResidue(cname, rnum)
+ if cname == "B" and (rnum.num >= 157 and rnum.num <= 168):
+ self.assertFalse(trimmed_r.IsValid())
+ else:
+ self.assertTrue(trimmed_r.IsValid())
+
+ def test_scorer_tmscore(self):
+ mdl = _LoadFile("1eud_mdl_partial-dimer.pdb")
+ trg = _LoadFile("1eud_ref.pdb")
+ scorer = Scorer(mdl, trg)
+ self.assertAlmostEqual(scorer.tm_score, 0.711, 3)
+
+if __name__ == "__main__":
+ from ost import testutils
+ if testutils.DefaultCompoundLibIsSet():
+ testutils.RunTests()
+ else:
+ print('No compound lib available. Ignoring test_scoring.py tests.')
--
GitLab