diff --git a/actions/ost-compare-structures b/actions/ost-compare-structures index a16c76f06ce60137bff4fc2f9ac9de443e672585..8364be451edcec48e10512166b97d8e03374e502 100644 --- a/actions/ost-compare-structures +++ b/actions/ost-compare-structures @@ -743,12 +743,11 @@ def _Main(): ost.LogInfo(str(lddt_settings).rstrip()) ost.LogInfo("===") oligo_lddt_scorer = qs_scorer.GetOligoLDDTScorer(lddt_settings) - for scorer_index, lddt_scorer in enumerate( - oligo_lddt_scorer.sc_lddt_scorers): - # Get chains and renumber according to alignment (for lDDT) - model_chain = lddt_scorer.model.chains[0].GetName() - reference_chain = \ - lddt_scorer.references[0].chains[0].GetName() + for mapped_lddt_scorer in oligo_lddt_scorer.mapped_lddt_scorers: + # Get data + lddt_scorer = mapped_lddt_scorer.lddt_scorer + model_chain = mapped_lddt_scorer.model_chain_name + reference_chain = mapped_lddt_scorer.reference_chain_name if skip_score: ost.LogInfo( " --> Skipping single chain lDDT because " @@ -785,8 +784,7 @@ def _Main(): "total_contacts": lddt_scorer.total_contacts} if opts.save_per_residue_scores: per_residue_sc = \ - oligo_lddt_scorer.GetPerResidueScores( - scorer_index) + mapped_lddt_scorer.GetPerResidueScores() ost.LogInfo("Per residue local lDDT (reference):") ost.LogInfo("Chain\tResidue Number\tResidue Name" "\tlDDT\tConserved Contacts\tTotal " diff --git a/modules/mol/alg/pymod/qsscoring.py b/modules/mol/alg/pymod/qsscoring.py index 8686ff05c3f4b2d7cc11949f401db30226c5b495..aa55203a930c5ee20372df794a38a2bfc2a9af55 100644 --- a/modules/mol/alg/pymod/qsscoring.py +++ b/modules/mol/alg/pymod/qsscoring.py @@ -355,8 +355,9 @@ class QSscorer: There will be one alignment for each mapped chain and they are ordered by their chain names in :attr:`qs_ent_1`. - The sequences of the alignments have views attached into - :attr:`QSscoreEntity.ent` of :attr:`qs_ent_1` and :attr:`qs_ent_2`. + The sequences of the alignments are named according to the chain name and + have views attached into :attr:`QSscoreEntity.ent` of :attr:`qs_ent_1` and + :attr:`qs_ent_2`. :getter: Computed on first use (cached) :type: :class:`list` of :class:`~ost.seq.AlignmentHandle` @@ -832,6 +833,9 @@ def GetContacts(entity, calpha_only, dist_thr=12.0): # DONE return contacts +############################################################################### +# Oligo-lDDT scores +############################################################################### class OligoLDDTScorer(object): """A simple class to calculate oligomeric lDDT score.""" @@ -855,7 +859,7 @@ class OligoLDDTScorer(object): self._lddt_ref = None self._lddt_mdl = None self._oligo_lddt_scorer = None - self._sc_lddt_scorers = None + self._mapped_lddt_scorers = None @property def lddt_ref(self): @@ -890,48 +894,136 @@ class OligoLDDTScorer(object): return self._oligo_lddt_scorer @property - def sc_lddt_scorers(self): - if self._sc_lddt_scorers is None: - self._sc_lddt_scorers = list() + def mapped_lddt_scorers(self): + if self._mapped_lddt_scorers is None: + self._mapped_lddt_scorers = list() for aln in self.alignments: - # Get chains and renumber according to alignment (for lDDT) - reference = Renumber( - aln.GetSequence(0), - old_number_label=self._old_number_label).CreateFullView() - refseq = seq.CreateSequence( - "reference_renumbered", - aln.GetSequence(0).GetString()) - refseq.AttachView(reference) - aln.AddSequence(refseq) - model = Renumber( - aln.GetSequence(1), - old_number_label=self._old_number_label).CreateFullView() - modelseq = seq.CreateSequence( - "model_renumbered", - aln.GetSequence(1).GetString()) - modelseq.AttachView(model) - aln.AddSequence(modelseq) - # Filter to CA-only if desired (done after AttachView to not mess it up) - if self.calpha_only: - lddt_scorer = lDDTScorer( - references=[reference.Select('aname=CA')], - model=model.Select('aname=CA'), - settings=self.settings) - else: - lddt_scorer = lDDTScorer( - references=[reference], - model=model, - settings=self.settings) - lddt_scorer.alignment = aln # a bit of a hack - self._sc_lddt_scorers.append(lddt_scorer) - return self._sc_lddt_scorers - - def GetPerResidueScores(self, scorer_index): + mapped_lddt_scorer = MappedLDDTScorer(aln, self.calpha_only, + self.settings) + self._mapped_lddt_scorers.append(mapped_lddt_scorer) + return self._mapped_lddt_scorers + + @property + def sc_lddt_scorers(self): + return [mls.lddt_scorer for mls in self.mapped_lddt_scorers] + + @property + def sc_lddt(self): + if self._sc_lddt is None: + self._sc_lddt = list() + for lddt_scorer in self.sc_lddt_scorers: + try: + self._sc_lddt.append(lddt_scorer.global_score) + except Exception as ex: + LogError('Single chain lDDT failed:', str(ex)) + self._sc_lddt.append(0.0) + return self._sc_lddt + + @property + def weighted_lddt(self): + if self._weighted_lddt is None: + scores = [s.global_score for s in self.sc_lddt_scorers] + weights = [s.total_contacts for s in self.sc_lddt_scorers] + nominator = sum([s * w for s, w in zip(scores, weights)]) + denominator = sum(weights) + if denominator > 0: + self._weighted_lddt = nominator / float(denominator) + else: + self._weighted_lddt = 0.0 + return self._weighted_lddt + + ############################################################################## + # Class internal helpers (anything that doesnt easily work without this class) + ############################################################################## + + def _PrepareOligoEntities(self): + # simple wrapper to avoid code duplication + self._lddt_ref, self._lddt_mdl = _MergeAlignedChains(self.alignments, + self.ref, + self.mdl, + self.calpha_only) + + +class MappedLDDTScorer(object): + """A simple class to calculate a single-chain lDDT score on a given chain to + chain mapping as extracted from :class:`OligoLDDTScorer`. + + :param alignment: Sets :attr:`alignment` + :param calpha_only: Sets :attr:`calpha_only` + :param settings: Sets :attr:`settings` + + .. attribute:: alignment + + Alignment with two sequences named according to the mapped chains and with + views attached to both sequences (e.g. one of the items of + :attr:`QSScorer.alignments`). + + The first sequence is assumed to be the reference and the second one the + model. Since the lDDT score is not symmetric (extra residues in model are + ignored), the order is important. + + :type: :class:`~ost.seq.AlignmentHandle` + + .. attribute:: calpha_only + + If true, restricts lDDT score to CA only. + + :type: :class:`bool` + + .. attribute:: settings + + Settings to use for lDDT scoring. + + :type: :class:`~mol.alg.lDDTSettings` + + .. attribute:: lddt_scorer + + lDDT Scorer object for the given chains. + + :type: :class:`~mol.alg.lDDTScorer` + + .. attribute:: reference_chain_name + + Chain name of the reference. + + :type: :class:`str` + + .. attribute:: model_chain_name + + Chain name of the model. + + :type: :class:`str` + """ + def __init__(self, alignment, calpha_only, settings): + # prepare fields + self.alignment = alignment + self.calpha_only = calpha_only + self.settings = settings + self.lddt_scorer = None # set in _InitScorer + self.reference_chain_name = alignment.sequences[0].name + self.model_chain_name = alignment.sequences[1].name + self._old_number_label = "old_num" + self._extended_alignment = None # set in _InitScorer + # initialize lDDT scorer + self._InitScorer() + + def GetPerResidueScores(self): + """ + :return: Scores for each residue + :rtype: :class:`list` of :class:`dict` with one item for each residue + existing in model and reference: + + - "residue_number": Residue number in reference chain + - "residue_name": Residue number in reference chain + - "lddt": local lDDT + - "conserved_contacts": number of conserved contacts + - "total_contacts": total number of contacts + """ scores = list() assigned_residues = list() # Make sure the score is calculated - self.sc_lddt_scorers[scorer_index].global_score - for col in self.sc_lddt_scorers[scorer_index].alignment: + self.lddt_scorer.global_score + for col in self._extended_alignment: if col[0] != "-" and col.GetResidue(3).IsValid(): ref_res = col.GetResidue(0) mdl_res = col.GetResidue(1) @@ -968,42 +1060,43 @@ class OligoLDDTScorer(object): "total_contacts": mdl_res_renum.GetFloatProp(self.settings.label + "_total")}) return scores - @property - def sc_lddt(self): - if self._sc_lddt is None: - self._sc_lddt = list() - for lddt_scorer in self.sc_lddt_scorers: - try: - self._sc_lddt.append(lddt_scorer.global_score) - except Exception as ex: - LogError('Single chain lDDT failed:', str(ex)) - self._sc_lddt.append(0.0) - return self._sc_lddt - - @property - def weighted_lddt(self): - if self._weighted_lddt is None: - scores = [s.global_score for s in self.sc_lddt_scorers] - weights = [s.total_contacts for s in self.sc_lddt_scorers] - nominator = sum([s * w for s, w in zip(scores, weights)]) - denominator = sum(weights) - if denominator > 0: - self._weighted_lddt = nominator / float(denominator) - else: - self._weighted_lddt = 0.0 - return self._weighted_lddt - ############################################################################## # Class internal helpers (anything that doesnt easily work without this class) ############################################################################## - def _PrepareOligoEntities(self): - # simple wrapper to avoid code duplication - self._lddt_ref, self._lddt_mdl = _MergeAlignedChains(self.alignments, - self.ref, - self.mdl, - self.calpha_only) - + def _InitScorer(self): + # Use copy of alignment (extended by 2 extra sequences for renumbering) + aln = self.alignment.Copy() + # Get chains and renumber according to alignment (for lDDT) + reference = Renumber( + aln.GetSequence(0), + old_number_label=self._old_number_label).CreateFullView() + refseq = seq.CreateSequence( + "reference_renumbered", + aln.GetSequence(0).GetString()) + refseq.AttachView(reference) + aln.AddSequence(refseq) + model = Renumber( + aln.GetSequence(1), + old_number_label=self._old_number_label).CreateFullView() + modelseq = seq.CreateSequence( + "model_renumbered", + aln.GetSequence(1).GetString()) + modelseq.AttachView(model) + aln.AddSequence(modelseq) + # Filter to CA-only if desired (done after AttachView to not mess it up) + if self.calpha_only: + self.lddt_scorer = lDDTScorer( + references=[reference.Select('aname=CA')], + model=model.Select('aname=CA'), + settings=self.settings) + else: + self.lddt_scorer = lDDTScorer( + references=[reference], + model=model, + settings=self.settings) + # Store alignment for later + self._extended_alignment = aln ############################################################################### # HELPERS @@ -2528,4 +2621,4 @@ def _ComputeLDDTScore(ref, mdl): # specify public interface __all__ = ('QSscoreError', 'QSscorer', 'QSscoreEntity', 'FilterContacts', - 'GetContacts', 'OligoLDDTScorer') + 'GetContacts', 'OligoLDDTScorer', 'MappedLDDTScorer')