diff --git a/pipeline/pymod/run_engine.py b/pipeline/pymod/run_engine.py index ae0659e16f1d856a552a16f6db22c759fc74f8f5..a8e32b86026e85b3822c19aac6989f8f9d441d3b 100644 --- a/pipeline/pymod/run_engine.py +++ b/pipeline/pymod/run_engine.py @@ -4,6 +4,7 @@ import ost from promod3 import loop +from promod3 import rawmodel def BuildFromRawModel(raw_model, structure_db=None, fragment_db=None, torsion_sampler=None, merge_distance=4): @@ -53,19 +54,19 @@ def BuildFromRawModel(raw_model, structure_db=None, fragment_db=None, # step 1: close small deletions in the raw model ost.LogVerbose("Trying to close small deletions (no. of gaps: %d)." % \ len(raw_model.gaps)) - scorer = raw_model.CloseSmallDeletions(scorer) + scorer = rawmodel.CloseSmallDeletions(raw_model, scorer) # step 2: simple handling of further gaps: merge, then fill from db for distance in range(merge_distance): # step 2a: Combine gaps living close to each other ost.LogVerbose("Trying to merge gaps (%d) with distance %d." % \ (len(raw_model.gaps), distance)) - raw_model.MergeGapsByDistance(distance) + rawmodel.MergeGapsByDistance(raw_model, distance) ost.LogVerbose("%d gap(s) left after merging.\n" % len(raw_model.gaps)+ "Trying to fill loops by database") # step 2b: fit loops into the model - scorer = raw_model.FillLoopsByDatabase(scorer, fragment_db, - structure_db, torsion_sampler, - max_loops_to_search=30) + scorer = rawmodel.FillLoopsByDatabase(raw_model, scorer, fragment_db, + structure_db, torsion_sampler, + max_loops_to_search=30) ost.LogVerbose("%d gap(s) left after database search." % \ len(raw_model.gaps)) if len(raw_model.gaps) == 0: diff --git a/rawmodel/doc/index.rst b/rawmodel/doc/index.rst index 42b7d863d7d759c8370ef4a926eb3f25faeaa587..b36f700adecd7c3c7b588e66a53a7be6ac5e22b9 100644 --- a/rawmodel/doc/index.rst +++ b/rawmodel/doc/index.rst @@ -1,48 +1,28 @@ -:mod:`~promod3.rawmodel` - Coordinate Modelling +:mod:`~promod3.rawmodel` - Protein Modelling ================================================================================ .. module:: promod3.rawmodel - :synopsis: Raw Coordinate Model + :synopsis: Protein Modelling .. currentmodule:: promod3.rawmodel -Functionality to build raw (pseudo) models based on a sequence alignment. -Here is an example of how to build a model from an alignment and a structure. +High-level functionality for protein modelling. +Commonly, your input is a template structure and an alignment of the template to +the desired target sequence. +A protein homology modelling pipeline then has the following main steps: -.. testcode:: rawmodel - :hide: +- Build a raw model from the template (see :func:`BuildRawModel` function) +- Perform loop modelling to close all gaps (see :func:`FillLoopsByDatabase` and TODO function) +- Reconstruct sidechains (see TODO function) +- Minimize energy of final model using molecular mechanics + (see TODO function) - import os - import tempfile - from ost import io - from promod3 import rawmodel - aln = io.LoadAlignment('../tests/rawmodel/data/raw-modeling/seq.fasta') - template_structure = io.LoadPDB('../tests/rawmodel/data/raw-modeling/gly.pdb', - restrict_chains='A') - aln.AttachView(1, template_structure.Select('peptide=true')) - result = rawmodel.BuildRawModel(aln) - (fh, fn) = tempfile.mkstemp(suffix='.pdb') - io.SavePDB(result.model, fn) - os.remove(fn) - -.. doctest:: rawmodel - - from ost import io - from promod3 import rawmodel - - aln = io.LoadAlignment('seq.fasta') - template_structure = io.LoadPDB('gly.pdb', restrict_chains='A') - aln.AttachView(1, template_structure.Select('peptide=true')) - result = rawmodel.BuildRawModel(aln) - io.SavePDB(result.model, 'model.pdb') - - -Raw Coordinate Modelling API +Modelling API -------------------------------------------------------------------------------- -.. function:: BuildRawModel(alignment, calpha_only=False) - BuildRawModel(alignments, calpha_only=False) +.. function:: BuildRawModel(alignment) + BuildRawModel(alignments) Builds a raw (pseudo) model from the alignment. Can either take a single alignment handle or an alignment handle list. Every list item is treated as a @@ -68,12 +48,88 @@ Raw Coordinate Modelling API The returned :class:`ModellingHandle` stores the obtained raw model as well as information about insertions and deletions in the gaps list. - :param calpha_only: If true, only Calpha atoms will be copied. Side chains - and other backbone atoms are completely ignored. + .. testcode:: rawmodel + :hide: + + import os + import tempfile + from ost import io + from promod3 import rawmodel + + aln = io.LoadAlignment('../tests/rawmodel/data/raw-modeling/seq.fasta') + template_structure = io.LoadPDB('../tests/rawmodel/data/raw-modeling/gly.pdb', + restrict_chains='A') + aln.AttachView(1, template_structure.Select('peptide=true')) + result = rawmodel.BuildRawModel(aln) + (fh, fn) = tempfile.mkstemp(suffix='.pdb') + io.SavePDB(result.model, fn) + os.remove(fn) + + .. doctest:: rawmodel + + from ost import io + from promod3 import rawmodel + + aln = io.LoadAlignment('seq.fasta') + template_structure = io.LoadPDB('gly.pdb', restrict_chains='A') + aln.AttachView(1, template_structure.Select('peptide=true')) + result = rawmodel.BuildRawModel(aln) + io.SavePDB(result.model, 'model.pdb') + + :param alignment: Single alignment handle for raw model. + :type alignment: :class:`~ost.seq.AlignmentHandle` + :param alignments: List of alignment handles for raw model with multiple chains. + :type alignments: :class:`~ost.seq.AlignmentList` + + :return: Raw (pseudo) model from the alignment. + :rtype: :class:`ModellingHandle` + :raises: A :exc:`RuntimeError` when the second sequence does not have an - attached structure + attached structure + +.. autofunction:: CloseSmallDeletions + +.. autofunction:: MergeGapsByDistance + +.. autofunction:: FillLoopsByDatabase + +.. function:: ClearGaps(mhandle, gap) + + Removes all gaps from mhandle which are fully enclosed by given gap. + + :param mhandle: Modelling handle on which to apply change. + :type mhandle: :class:`ModellingHandle` + :param gap: Gap defining range in which gaps are to be removed. + :type gap: :class:`StructuralGap` + + :return: Index of next gap in mhandle.gaps after removal. + Returns -1 if last gap was removed. + :rtype: :class:`int` + + :raises: A :exc:`RuntimeError` if any gap in mhandle.gaps is only partially + enclosed by given gap. + +.. function:: MergeGaps(mhandle, index) + + Merges two gaps mhandle.gaps[index] and mhandle.gaps[index+1]. -Modelling Handle + :param mhandle: Modelling handle on which to apply change. + :type mhandle: :class:`ModellingHandle` + :param index: Index of gap to merge with next one. + :type index: :class:`int` + + :raises: A :exc:`RuntimeError` if indices out of range or if trying to merge + N-terminal gap with a C-terminal gap. + +.. function:: RemoveTerminalGaps(mhandle) + + :param mhandle: Modelling handle on which to apply change. + :type mhandle: :class:`ModellingHandle` + + :return: Number of gaps which were removed. + :rtype: :class:`int` + +Modelling Handle class -------------------------------------------------------------------------------- .. class:: ModellingHandle @@ -96,11 +152,188 @@ Modelling Handle :type: :class:`StructuralGapList` - .. automethod:: CloseSmallDeletions + .. attribute:: seqres + + List of sequences with one :class:`~ost.seq.SequenceHandle` for each chain + of target protein. + + :type: :class:`~ost.seq.SequenceList` + + +Gap classes +-------------------------------------------------------------------------------- + +.. class:: StructuralGap(before, after, seq) + + Describes a structural gap, i.e. a loop to be modeled. The gap may either be + terminal or between two defined regions. The gap stores information of the + last residue before and the first residue after the gap as well as the + sequence of gap. Gaps at the N- and C-terminals can be defined by passing + invalid residue handles to `before` or `after`. + + :param before: Fills :attr:`before` + :type before: :class:`ost.mol.ResidueHandle` + :param after: Fills :attr:`after` + :type after: :class:`ost.mol.ResidueHandle` + :param seq: Fills :attr:`seq` + :type seq: :class:`str` + + :raises: A :exc:`RuntimeError` if both residues are invalid or when both + are valid and: + + - residues are from different chains (if both valid) + - `before` is located after `after` + - `seq` has a length which is inconsistent with the gap + + .. method:: GetChainIndex() + + :return: Index of chain, the gap is belonging to + :rtype: :class:`int` + + .. method:: GetChainName() + + :return: Name of chain, the gap is belonging to + :rtype: :class:`str` + + .. method:: GetChain() + + :return: Chain, the gap is belonging to + :rtype: :class:`ost.mol.ChainHandle` + + .. method:: IsNTerminal() + + :return: True, iff gap is N-terminal (i.e. :attr:`before` is invalid + and :attr:`after` is valid) + :rtype: :class:`bool` + + .. method:: IsCTerminal() + + :return: True, iff gap is C-terminal (i.e. :attr:`before` is valid + and :attr:`after` is invalid) + :rtype: :class:`bool` + + .. method:: IsTerminal() + + :return: True, iff gap is N- or C-terminal + :rtype: :class:`bool` + + .. method:: ShiftCTerminal() + + Try to shift gap by one position towards C-terminal. Only possible if new + gap is not terminal and it doesn't try to shift the gap past another gap. + + :return: True, iff shift succeeded (gap is only updated in that case) + :rtype: :class:`bool` + + .. method:: ExtendAtNTerm() + + Try to extend gap at N-terminal end of gap. + Only possible if the gap is not at N-terminal and it doesn't try to + extend the gap past another gap. + + :return: True, iff extend succeeded (gap is only updated in that case) + :rtype: :class:`bool` + + .. method:: ExtendAtCTerm() + + Try to extend gap at C-terminal end of gap. + Only possible if the gap is not at C-terminal and it doesn't try to + extend the gap past another gap. + + :return: True, iff extend succeeded (gap is only updated in that case) + :rtype: :class:`bool` + + .. method:: GetLength() + + :return: Length of the gap. + :rtype: :class:`int` + + .. method:: Copy() + + :return: Copy of the gap. + :rtype: :class:`StructuralGap` + + .. attribute:: length + + Alias for :meth:`GetLength()` (read-only, :class:`int`) + + .. attribute:: seq + + Sequence string for the gap (read-only, :class:`str`) + + .. attribute:: before + + Residue before the gap (read-only, :class:`ost.mol.ResidueHandle`) + + .. attribute:: after + + Residue after the gap (read-only, :class:`ost.mol.ResidueHandle`) + + .. attribute:: full_seq + + Full sequence, including stem residues (read-only) + +.. class:: StructuralGapList + + Represents a :class:`list` of :class:`StructuralGap`. + + +Gap Extender classes +-------------------------------------------------------------------------------- + +The extender classes work on a given :class:`StructuralGap` and provide an +Extend() function to propose new gaps for loop modelling. The function returns +False if no new extension possible. + +.. class:: GapExtender(gap) + + The extender cycles through the following steps: + + .. code-block:: none + + - + -- + -- + --- + --- + --- + ---- + ---- + ---- + ---- + + :param gap: The gap which will be extended by :meth:`Extend`. + :type gap: :class:`StructuralGap` + + .. method:: Extend() + + Tries to extend `gap`. + + :return: False, iff the gap cannot be extended any further. + :rtype: :class:`bool` + +.. class:: ScoringGapExtender(gap, extension_penalty, penalties) + + The extender scores possible gap extensions and returns them in order of + their score when :meth:`Extend` is called. + The score is penalized according to length and according to certain (well + conserved) regions in the structure. + score = length * `extension_penalty` + sum( `penalties` [i] ) + (i = res.num. of residues in extension) + + :param gap: The gap which will be extended by :meth:`Extend`. + :type gap: :class:`StructuralGap` + :param extension_penalty: Penalty for length of gap. + :type extension_penalty: :class:`float` + :param penalties: Penalty for each residue added to gap. + :type penalties: :class:`list` of :class:`float` + + .. method:: Extend() - .. automethod:: MergeGapsByDistance + Tries to extend `gap`. - .. automethod:: FillLoopsByDatabase + :return: False, iff the gap cannot be extended any further. + :rtype: :class:`bool` .. LocalWords: currentmodule promod aln AttachView BuildRawModel pdb calpha .. LocalWords: ModellingHandle StructuralGapList rawmodel Modelling os ost diff --git a/rawmodel/pymod/__init__.py b/rawmodel/pymod/__init__.py index dc608ebbc8b44f72112c974c7443ad72b66242db..78a1f850948f6a3fd1aac61bc49d104a20fa4d4f 100644 --- a/rawmodel/pymod/__init__.py +++ b/rawmodel/pymod/__init__.py @@ -2,12 +2,4 @@ Initialise the rawmodel module. """ from _rawmodel import * -from _closegaps import _CloseSmallDeletions, _MergeGapsByDistance, \ - _FillLoopsByDatabase - -# this should attach _CloseSmallDeletions to ModellingHandle as class method. -setattr(ModellingHandle, 'CloseSmallDeletions', _CloseSmallDeletions) -# this should attach _MergeGapsByDistance to ModellingHandle as class method. -setattr(ModellingHandle, 'MergeGapsByDistance', _MergeGapsByDistance) -# this should attach _FillGapsByDatabase to ModellingHandle as class method. -setattr(ModellingHandle, 'FillLoopsByDatabase', _FillLoopsByDatabase) +from _closegaps import * diff --git a/rawmodel/pymod/_closegaps.py b/rawmodel/pymod/_closegaps.py index f15999da1f4c9f38f35b2ca41b6e1d0dcb63bd0d..6792787e65fc88a6c898d9b0d5e3781b4c1d4129 100644 --- a/rawmodel/pymod/_closegaps.py +++ b/rawmodel/pymod/_closegaps.py @@ -1,6 +1,6 @@ -'''Functions to be 'injected' into the RawModellingResult class in the -__init__.py file. Do not use directly, they belong to the RawModellingResult -class and are to be called as class method, for instances of this class. +'''Functions added as high-level functionality to modelling module in the +__init__.py file. To be used directly by passing a ModellingHandle instance +as argument. ''' import ost @@ -8,8 +8,8 @@ import ost from . import _rawmodel as rawmodel from promod3 import loop -def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, - e_thresh=200): +def CloseSmallDeletions(mhandle, scorer, extension_steps=9, clash_thresh=1.0, + e_thresh=200): '''Close small deletions by relaxing neighbouring residues. Small deletions in the template from the target-template alignment have a @@ -18,7 +18,7 @@ def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, already in the raw-model. After closure some checks are done to see if the solution is stereochemically sensible. - Closed gaps are removed from :attr:`self.gaps`. + Closed gaps are removed from :attr:`mhandle.gaps`. .. testcode:: closesmalldel :hide: @@ -30,11 +30,11 @@ def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, aln = ost.seq.CreateAlignment(ost.seq.CreateSequence('trg', 'GGG-GGG'), ost.seq.CreateSequence('tpl', 'GGGAGGG')) aln.AttachView(1, tpl.CreateFullView()) - rmodel = rawmodel.BuildRawModel(aln) - assert len(rmodel.gaps) == 1 - scorer = loop.SetupBackboneScorer(rmodel) - rmodel.CloseSmallDeletions(scorer) - assert len(rmodel.gaps) == 0 + mhandle = rawmodel.BuildRawModel(aln) + assert len(mhandle.gaps) == 1 + scorer = loop.SetupBackboneScorer(mhandle) + rawmodel.CloseSmallDeletions(mhandle, scorer) + assert len(mhandle.gaps) == 0 .. doctest:: closesmalldel @@ -45,10 +45,12 @@ def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, tpl = ost.io.LoadPDB('gly.pdb') aln = ost.io.LoadAlignment('seq.fasta') aln.AttachView(1, tpl.CreateFullView()) - rmodel = rawmodel.BuildRawModel(aln) - scorer = loop.SetupBackboneScorer(rmodel) - rmodel.CloseSmallDeletions(scorer) + mhandle = rawmodel.BuildRawModel(aln) + scorer = loop.SetupBackboneScorer(mhandle) + rawmodel.CloseSmallDeletions(mhandle, scorer) + :param mhandle: Modelling handle on which to apply change. + :type mhandle: :class:`ModellingHandle` :param scorer: A scorer dedicated to this raw model. :type scorer: :class:`~promod3.loop.BackboneLoopScorer` @@ -72,11 +74,11 @@ def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, ''' current_gap_index = 0 - # iterating self.gaps. The number of gaps may change during the process, + # iterating mhandle.gaps. The number of gaps may change during the process, # hence we run by 'while', comparing to the updated list of gaps. If a gap - # gets closed, it is deleted from self.gaps, otherwise, current_gap_index + # gets closed, it is deleted from mhandle.gaps, otherwise, current_gap_index # as a counter is increased. - while current_gap_index < len(self.gaps): + while current_gap_index < len(mhandle.gaps): # in the end this determines if a gap was closed or not success = False @@ -84,9 +86,9 @@ def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, # transform the template sequence into the target sequence, aa's vanish, # so the target sequence has gap caracters, the template not. # If we are not looking at a deletion, do nothing. - if len(self.gaps[current_gap_index].seq) == 0: + if len(mhandle.gaps[current_gap_index].seq) == 0: # work on a copy of the gap, if not closed in the end, no harm done - current_gap = self.gaps[current_gap_index].Copy() + current_gap = mhandle.gaps[current_gap_index].Copy() current_chain = current_gap.GetChain() current_chain_index = current_gap.GetChainIndex() @@ -132,14 +134,14 @@ def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, potential_e < e_thresh and \ scorer.TransOmegaTorsions(bb_list): ost.LogVerbose("Closed: %s by relaxing %s" % \ - (self.gaps[current_gap_index], current_gap)) + (mhandle.gaps[current_gap_index], current_gap)) chain = current_gap.before.GetChain() bb_list.InsertInto(chain, current_gap.before.GetNumber().GetNum(), current_gap.after.GetNumber().GetNum()) scorer.SetEnvironment(\ bb_list, current_gap.before.GetNumber().GetNum()) - self.ClearGaps(current_gap) + rawmodel.ClearGaps(mhandle, current_gap) success = True break @@ -151,7 +153,7 @@ def _CloseSmallDeletions(self, scorer, extension_steps=9, clash_thresh=1.0, return scorer -def _MergeGapsByDistance(self, distance): +def MergeGapsByDistance(mhandle, distance): '''Merge 2 neighbouring gaps by deleting residues in-between. Check if two neighbouring gaps are at max. *distance* residues apart from @@ -170,10 +172,10 @@ def _MergeGapsByDistance(self, distance): ost.seq.CreateSequence('tpl', 'NN----A----LF')) aln.AttachView(1, tpl.CreateFullView()) - rmodel = rawmodel.BuildRawModel(aln) - assert len(rmodel.gaps) == 2 - rmodel.MergeGapsByDistance(0) - assert len(rmodel.gaps) == 1 + mhandle = rawmodel.BuildRawModel(aln) + assert len(mhandle.gaps) == 2 + rawmodel.MergeGapsByDistance(mhandle, 0) + assert len(mhandle.gaps) == 1 .. doctest:: mergegapsbydist @@ -183,9 +185,11 @@ def _MergeGapsByDistance(self, distance): tpl = ost.io.LoadPDB('1mcg.pdb') aln = ost.io.LoadAlignment('1mcg_aln.fasta') aln.AttachView(1, tpl.CreateFullView()) - rmodel = rawmodel.BuildRawModel(aln) - rmodel.MergeGapsByDistance(0) + mhandle = rawmodel.BuildRawModel(aln) + rawmodel.MergeGapsByDistance(mhandle, 0) + :param mhandle: Modelling handle on which to apply change. + :type mhandle: :class:`ModellingHandle` :param distance: The max. no. of residues between two gaps up to which merge happens. :type distance: :class:`int` @@ -205,9 +209,9 @@ def _MergeGapsByDistance(self, distance): while try_again: try_again = False # iterate all but the last gap, since we are always looking ahead - for i in range(len(self.gaps) - 1): - current_gap = self.gaps[i] - next_gap = self.gaps[i+1] + for i in range(len(mhandle.gaps) - 1): + current_gap = mhandle.gaps[i] + next_gap = mhandle.gaps[i+1] # check that we are on the same chain if current_gap.GetChain() != next_gap.GetChain(): continue @@ -219,15 +223,15 @@ def _MergeGapsByDistance(self, distance): - current_gap.after.GetNumber().GetNum() if dist <= distance: # gaps are close enough, combine! combine! - self.MergeGaps(i) + rawmodel.MergeGaps(mhandle, i) ost.LogVerbose("Merged gap %s and %s into %s" % \ - (current_gap, next_gap, self.gaps[i])) + (current_gap, next_gap, mhandle.gaps[i])) try_again = True break -def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, - torsion_sampler, max_loops_to_search=30, - max_db_loop_len=12): +def FillLoopsByDatabase(mhandle, scorer, fragment_db, structure_db, + torsion_sampler, max_loops_to_search=30, + max_db_loop_len=12): '''Try to fill up loops from a structural database. Usually this will extend the gaps a bit to match candidates from the @@ -247,13 +251,13 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, aln.AttachView(1, tpl.CreateFullView()) - rmodel = rawmodel.BuildRawModel(aln) - assert len(rmodel.gaps) == 1 - scorer = loop.SetupBackboneScorer(rmodel) - rmodel.FillLoopsByDatabase(scorer, loop.LoadFragDB(), - loop.LoadStructureDB(), - loop.LoadTorsionSamplerCoil()) - assert len(rmodel.gaps) == 0 + mhandle = rawmodel.BuildRawModel(aln) + assert len(mhandle.gaps) == 1 + scorer = loop.SetupBackboneScorer(mhandle) + rawmodel.FillLoopsByDatabase(mhandle, scorer, loop.LoadFragDB(), + loop.LoadStructureDB(), + loop.LoadTorsionSamplerCoil()) + assert len(mhandle.gaps) == 0 .. doctest:: fillloopsbydb @@ -264,11 +268,14 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, aln = ost.io.LoadAlignment('2dbs.fasta') aln.AttachView(1, tpl.CreateFullView()) - rmodel = rawmodel.BuildRawModel(aln) - scorer = loop.SetupBackboneScorer(rmodel) - scorer = rmodel.FillLoopsByDatabase(scorer, loop.LoadFragDB(), - loop.LoadStructureDB(), - loop.LoadTorsionSamplerCoil()) + mhandle = rawmodel.BuildRawModel(aln) + scorer = loop.SetupBackboneScorer(mhandle) + scorer = rawmodel.FillLoopsByDatabase(mhandle, scorer, loop.LoadFragDB(), + loop.LoadStructureDB(), + loop.LoadTorsionSamplerCoil()) + + :param mhandle: Modelling handle on which to apply change. + :type mhandle: :class:`ModellingHandle` :param scorer: A scorer dedicated to this raw model. :type scorer: :class:`~promod3.loop.BackboneLoopScorer` @@ -335,18 +342,18 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, # Iterate all gaps. Since the number of gaps may change, always compare to # an updated list. gap_idx is only increased when necessary, e.g. current # gap could not be removed. - while gap_idx < len(self.gaps): + while gap_idx < len(mhandle.gaps): # count how many loops have been found/ searched so far, needed as a # stop criterion at some point to start investigating found loops. found_loops = 0 # terminal gaps are no loops, database needs 2 stems to work on - if self.gaps[gap_idx].IsNTerminal() or self.gaps[gap_idx].IsCTerminal(): + if mhandle.gaps[gap_idx].IsNTerminal() or mhandle.gaps[gap_idx].IsCTerminal(): gap_idx += 1 continue # get info on current gap - actual_gap = self.gaps[gap_idx].Copy() + actual_gap = mhandle.gaps[gap_idx].Copy() actual_length = actual_gap.length actual_extender = rawmodel.GapExtender(actual_gap) min_before_resnum = 100000000000 @@ -354,7 +361,7 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, actual_candidates = list() actual_extended_gaps = list() actual_chain_idx = actual_gap.GetChainIndex() - actual_chain = self.model.chains[actual_chain_idx] + actual_chain = mhandle.model.chains[actual_chain_idx] # Find loop candidates. while actual_length <= max_db_loop_len: @@ -407,7 +414,7 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, # all loop candidates will be scored along the full max. extension ever # reached in the search before, so we build an overall frame, where we # insert the loops - frame_seq = self.seqres[actual_chain_idx]\ + frame_seq = mhandle.seqres[actual_chain_idx]\ [min_before_resnum-1:max_after_resnum] frame_backbone_list = loop.BackboneList(frame_seq) actual_res_num = ost.mol.ResNum(min_before_resnum) @@ -466,7 +473,7 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, if optimal_candidate != -1: idx_a = back_mapper[optimal_candidate][0] idx_b = back_mapper[optimal_candidate][1] - actual_candidates[idx_a].InsertInto(self.model, idx_b) + actual_candidates[idx_a].InsertInto(mhandle.model, idx_b) bb_list = actual_candidates[idx_a][idx_b].bb_list scorer.SetEnvironment(\ bb_list, @@ -474,7 +481,7 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, actual_chain_idx) ost.LogVerbose("Resolved %s by filling %s" % \ (str(actual_gap), str(actual_extended_gaps[idx_a]))) - next_gap = self.ClearGaps(actual_extended_gaps[idx_a]) + next_gap = rawmodel.ClearGaps(mhandle, actual_extended_gaps[idx_a]) if next_gap == -1: break else: @@ -484,12 +491,12 @@ def _FillLoopsByDatabase(self, scorer, fragment_db, structure_db, gap_idx += 1 return scorer -__all__ = ['_CloseSmallDeletions', '_MergeGapsByDistance', - '_FillLoopsByDatabase'] +__all__ = ['CloseSmallDeletions', 'MergeGapsByDistance', + 'FillLoopsByDatabase'] -# LocalWords: modeling stereochemically param idx RawModellingResult init +# LocalWords: modeling stereochemically param idx init # LocalWords: py ost pylint rawmodel promod CloseSmallDeletions testcode -# LocalWords: fillloopsbydb tpl aln trg TLNGFTVPAGNTLV LNPDKGATVTMA rmodel +# LocalWords: fillloopsbydb tpl aln trg TLNGFTVPAGNTLV LNPDKGATVTMA mhandle # LocalWords: NGGTLLIPNGTYHFLGIQMKSNVHIRVE AttachView CreateFullView len # LocalWords: BuildRawModel SetupBackboneScorer FillLoopsByDatabase doctest # LocalWords: LoadFragDB LoadStructureDB LoadTorsionSamplerCoil dbs fasta diff --git a/rawmodel/pymod/export_gap.cc b/rawmodel/pymod/export_gap.cc index caf2e22346f616036843a719eca1fd9239969351..27192da45454ebfe32df110bf01a2c0a721d78fd 100644 --- a/rawmodel/pymod/export_gap.cc +++ b/rawmodel/pymod/export_gap.cc @@ -29,7 +29,6 @@ void export_gap() .def("ExtendAtNTerm", &StructuralGap::ExtendAtNTerm) .def("ExtendAtCTerm", &StructuralGap::ExtendAtCTerm) .def("GetLength", &StructuralGap::GetLength) - .def("Transfer", &StructuralGap::Transfer) .def("Copy", &StructuralGap::Copy) .add_property("length", &StructuralGap::GetLength) .def_readonly("seq", &StructuralGap::sequence) diff --git a/rawmodel/pymod/export_model.cc b/rawmodel/pymod/export_model.cc index a153ba8d45dca7ea6e1dd503053139fb81593f46..7cce177126a55d556e016bc0f19d319fc9db5be3 100644 --- a/rawmodel/pymod/export_model.cc +++ b/rawmodel/pymod/export_model.cc @@ -23,10 +23,10 @@ void export_model() .def_readwrite("model", &ModellingHandle::model) .def_readwrite("gaps", &ModellingHandle::gaps) .def_readwrite("seqres", &ModellingHandle::seqres) - .def("ClearGaps",&ModellingHandle::ClearGaps,(arg("gap"))) - .def("MergeGaps",&ModellingHandle::MergeGaps,(arg("index"))) - .def("RemoveTerminalGaps",&ModellingHandle::RemoveTerminalGaps) ; + def("ClearGaps",ClearGaps,(arg("gap"))); + def("MergeGaps",MergeGaps,(arg("index"))); + def("RemoveTerminalGaps",RemoveTerminalGaps); def("BuildRawModel", BuildRawModelHandle, (arg("aln"), arg("include_ligands")=false, diff --git a/rawmodel/src/gap.cc b/rawmodel/src/gap.cc index 3713839ddc3bcb22d6c33cc48e92507acf614d29..af3d103d3f860f8bf85479ab0280b084231648f7 100644 --- a/rawmodel/src/gap.cc +++ b/rawmodel/src/gap.cc @@ -62,6 +62,11 @@ bool StructuralGap::ExtendAtCTerm() if (!next.IsValid()) { return false; } + ResNum next_num=next.GetNumber(); + ResNum exp_num=after.GetNumber()+1; + if (next_num!=exp_num) { + return false; + } sequence+=after.GetOneLetterCode(); after=next; return true; @@ -77,6 +82,11 @@ bool StructuralGap::ExtendAtNTerm() if (!prev.IsValid()) { return false; } + ResNum next_num=prev.GetNumber(); + ResNum exp_num=before.GetNumber()-1; + if (next_num!=exp_num) { + return false; + } sequence=String(1, before.GetOneLetterCode())+sequence; before=prev; return true; @@ -119,26 +129,4 @@ String StructuralGap::AsString() const return ss.str(); } -/// \brief returns a copy of the gap, with residues pointing to ent, instead -/// of the current entity -StructuralGap StructuralGap::Transfer(const EntityHandle& ent) const -{ - ChainHandle chain=ent.GetChainList()[0]; - ResidueHandle new_before; - if (before) { - new_before=chain.FindResidue(before.GetNumber()); - if (!new_before.IsValid()) { - LOG_WARNING("residue " << before << " doesn't exist"); - } - } - ResidueHandle new_after; - if (after) { - new_after=chain.FindResidue(after.GetNumber()); - if (!new_after.IsValid()) { - LOG_WARNING("residue " << after << " doesn't exist"); - } - } - return StructuralGap(new_before, new_after, sequence); -} - }} diff --git a/rawmodel/src/gap.hh b/rawmodel/src/gap.hh index 684a1c7ae5e0b987cee5b4ae506581bbb348b856..86f852366b1b98e97dc4e0569b46a5767fdb0c33 100644 --- a/rawmodel/src/gap.hh +++ b/rawmodel/src/gap.hh @@ -75,15 +75,28 @@ struct StructuralGap { /// \brief get length of the gap int GetLength() const { return sequence.length(); } - /// \brief extend the gap at the C-terminus. + /// \brief Try to extend gap at C-terminal end of gap. /// - /// returns true if the gap could be extended, false if not. + /// Only possible if the gap is not at C-terminal and it doesn't try to + /// extend the gap past another gap. + /// + /// \returns True, iff extend succeeded (gap is only updated in that case) bool ExtendAtCTerm(); - /// \ brief Extend gap at N-terminal end of gap. Only possible if the - /// C-terminal end points to a valid residue. + /// \brief Try to extend gap at N-terminal end of gap. + /// + /// Only possible if the gap is not at N-terminal and it doesn't try to + /// extend the gap past another gap. + /// + /// \returns True, iff extend succeeded (gap is only updated in that case) bool ExtendAtNTerm(); + /// \brief Try to shift gap by one position towards C-terminal. + /// + /// Only possible if new gap is not terminal and it doesn't try to shift the + /// gap past another gap. + /// + /// \returns True, iff shift succeeded (gap is only updated in that case) bool ShiftCTerminal(); String AsString() const; @@ -114,10 +127,6 @@ struct StructuralGap { if (after) { return after.GetEntity(); } return ost::mol::EntityHandle(); } - - /// \brief returns a copy of the gap, with residues pointing to ent, instead - /// of the current entity - StructuralGap Transfer(const ost::mol::EntityHandle& ent) const; StructuralGap Copy() { return StructuralGap(before, after, sequence); } diff --git a/rawmodel/src/gap_extender.hh b/rawmodel/src/gap_extender.hh index 33523af270a98a33892e1333f75aafeb0dd23567..c8656355f6e6b0de408fd01444a106667843cc63 100644 --- a/rawmodel/src/gap_extender.hh +++ b/rawmodel/src/gap_extender.hh @@ -33,8 +33,6 @@ typedef boost::shared_ptr<ScoringGapExtender> ScoringGapExtenderPtr; /// /// When the gap can not be extended any further, GapExtender::Extend() /// returns false. -/// -/// Note that Extend may turn an internal gap into a terminal gap. class GapExtender { public: GapExtender(StructuralGap& g): gap(g), initial_length(g.GetLength()), diff --git a/rawmodel/src/model.cc b/rawmodel/src/model.cc index 3d01b00ec11ac346731d0a88c96d9e31e2014132..6c6729bced09a00a0f58266be42999857a37769e 100644 --- a/rawmodel/src/model.cc +++ b/rawmodel/src/model.cc @@ -64,13 +64,13 @@ bool CheckCalphaAtom(ResidueView res) } -int ModellingHandle::ClearGaps(const StructuralGap& gap){ +int ClearGaps(ModellingHandle& mhandle, const StructuralGap& gap) { String chain_name = gap.GetChainName(); int next_gap = -1; if(gap.IsNTerminal()){ - for(StructuralGapList::iterator i = gaps.begin(); i != gaps.end();){ + for(StructuralGapList::iterator i = mhandle.gaps.begin(); i != mhandle.gaps.end();){ if(i->GetChainName() != chain_name){ ++i; continue; @@ -83,14 +83,14 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ ss << " is only overlapping not fully enclosing!"; throw promod3::Error(ss.str()); } - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } }else if(i->IsNTerminal()){ if(i->after.GetNumber() <= gap.after.GetNumber()){ - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } }else{ @@ -102,8 +102,8 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ throw promod3::Error(ss.str()); } } - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } ++i; @@ -112,7 +112,7 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ } if(gap.IsCTerminal()){ - for(StructuralGapList::iterator i = gaps.begin(); i != gaps.end();){ + for(StructuralGapList::iterator i = mhandle.gaps.begin(); i != mhandle.gaps.end();){ if(i->GetChainName() != chain_name){ ++i; continue; @@ -125,14 +125,14 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ ss << " is only overlapping not fully enclosing!"; throw promod3::Error(ss.str()); } - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } }else if(i->IsCTerminal()){ if(i->before.GetNumber() >= gap.before.GetNumber()){ - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } }else{ @@ -144,8 +144,8 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ throw promod3::Error(ss.str()); } } - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } ++i; @@ -153,7 +153,7 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ return next_gap; } - for(StructuralGapList::iterator i = gaps.begin(); i != gaps.end();){ + for(StructuralGapList::iterator i = mhandle.gaps.begin(); i != mhandle.gaps.end();){ if(i->GetChainName() != chain_name){ ++i; continue; @@ -172,8 +172,8 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ ss << " is only overlapping not fully enclosing!"; throw promod3::Error(ss.str()); } - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } } @@ -190,15 +190,15 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ ss << " is only overlapping not fully enclosing!"; throw promod3::Error(ss.str()); } - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } }else{ if(i->before.GetNumber() >= gap.before.GetNumber() && i->after.GetNumber() <= gap.after.GetNumber()){ - i = gaps.erase(i); - next_gap = i-gaps.begin(); + i = mhandle.gaps.erase(i); + next_gap = i-mhandle.gaps.begin(); continue; } if(i->before.GetNumber() > gap.before.GetNumber() && @@ -222,38 +222,38 @@ int ModellingHandle::ClearGaps(const StructuralGap& gap){ return next_gap; } -void ModellingHandle::MergeGaps(uint index){ +void MergeGaps(ModellingHandle& mhandle, uint index) { - if(index >= gaps.size()-1){ + if(index >= mhandle.gaps.size()-1){ throw promod3::Error("Invalid index provided when merging gaps!"); } - if(gaps[index].GetChain() != gaps[index+1].GetChain()){ + if(mhandle.gaps[index].GetChain() != mhandle.gaps[index+1].GetChain()){ throw promod3::Error("Subsequent gap must be of the same chain to be merged in!"); } - if(gaps[index].IsNTerminal() && gaps[index+1].IsCTerminal()){ + if(mhandle.gaps[index].IsNTerminal() && mhandle.gaps[index+1].IsCTerminal()){ std::stringstream ss; - ss << "You don't want to merge an NTerminal gap with a CTerminal Gap... "; + ss << "You don't want to merge an N-terminal gap with a C-terminal gap... "; ss << "You might want to use the RemoveTerminalGaps function to get rid of them!"; throw promod3::Error(ss.str()); } //assumes, that the gaps are sequential! - String full_seq = seqres[gaps[index].GetChainIndex()].GetGaplessString(); + String full_seq = mhandle.seqres[mhandle.gaps[index].GetChainIndex()].GetGaplessString(); int before_number, after_number; - if(gaps[index].IsNTerminal()) before_number = 0; - else before_number = gaps[index].before.GetNumber().GetNum(); + if(mhandle.gaps[index].IsNTerminal()) before_number = 0; + else before_number = mhandle.gaps[index].before.GetNumber().GetNum(); - if(gaps[index+1].IsCTerminal()) after_number = full_seq.size()+1; - else after_number = gaps[index+1].after.GetNumber().GetNum(); + if(mhandle.gaps[index+1].IsCTerminal()) after_number = full_seq.size()+1; + else after_number = mhandle.gaps[index+1].after.GetNumber().GetNum(); String gap_seq = full_seq.substr(before_number,after_number-before_number-1); //kill all residues in between the two new stems - ost::mol::XCSEditor ed = model.EditXCS(ost::mol::BUFFERED_EDIT); - ost::mol::ChainHandle chain = gaps[index].GetChain(); + ost::mol::XCSEditor ed = mhandle.model.EditXCS(ost::mol::BUFFERED_EDIT); + ost::mol::ChainHandle chain = mhandle.gaps[index].GetChain(); ost::mol::ResNum i(before_number+1); ost::mol::ResNum e(after_number); for(; i < e; ++i){ @@ -261,20 +261,20 @@ void ModellingHandle::MergeGaps(uint index){ if(res.IsValid()) ed.DeleteResidue(res); } - StructuralGap new_gap(gaps[index].before, - gaps[index+1].after, + StructuralGap new_gap(mhandle.gaps[index].before, + mhandle.gaps[index+1].after, gap_seq); - gaps[index] = new_gap; + mhandle.gaps[index] = new_gap; - gaps.erase(gaps.begin()+index+1); + mhandle.gaps.erase(mhandle.gaps.begin()+index+1); } -int ModellingHandle::RemoveTerminalGaps(){ +int RemoveTerminalGaps(ModellingHandle& mhandle) { int removed_gaps = 0; - for(StructuralGapList::iterator i = gaps.begin(); i != gaps.end(); ){ + for(StructuralGapList::iterator i = mhandle.gaps.begin(); i != mhandle.gaps.end(); ){ if(i->IsNTerminal() || i->IsCTerminal()){ - i = gaps.erase(i); + i = mhandle.gaps.erase(i); ++removed_gaps; } else ++i; diff --git a/rawmodel/src/model.hh b/rawmodel/src/model.hh index 74eb2b6b90825e5b3a9bb5a238a17494f84f13cf..fa42133843b77b8abd557a5cac70a77f051915a1 100644 --- a/rawmodel/src/model.hh +++ b/rawmodel/src/model.hh @@ -27,17 +27,31 @@ struct ModellingHandle { return !this->operator==(rhs); } - int ClearGaps(const StructuralGap& gap); - - void MergeGaps(uint index); - - int RemoveTerminalGaps(); - ost::mol::EntityHandle model; StructuralGapList gaps; ost::seq::SequenceList seqres; }; +/// \brief Removes all gaps from mhandle which are fully enclosed by given gap. +/// +/// \returns Index of next gap in mhandle.gaps after removal. +/// Returns -1 if last gap was removed. +/// +/// \throws promod3::Error if any gap in mhandle.gaps is only partially +/// enclosed by given gap. +int ClearGaps(ModellingHandle& mhandle, const StructuralGap& gap); + +/// \brief Merges two gaps mhandle.gaps[index] and mhandle.gaps[index+1]. +/// +/// \throws promod3::Error if indices out of range or if trying to merge +/// N-terminal gap with a C-terminal gap. +void MergeGaps(ModellingHandle& mhandle, uint index); + +/// \brief Removes all gaps in mhandle which are at N-terminal or C-terminal. +/// +/// \returns Number of gaps which were removed. +int RemoveTerminalGaps(ModellingHandle& mhandle); + /// \brief copies all atom of src_res to dst_res /// \param has_cbeta will be set to true if the src_res has a cbeta and the // dst_residue is not a glycine diff --git a/rawmodel/tests/test_close_gaps.py b/rawmodel/tests/test_close_gaps.py index 643be41ca14895be6df43d7093fa16f56f13a311..c4b0155f7a11f607d57b27703530e997d647c225 100644 --- a/rawmodel/tests/test_close_gaps.py +++ b/rawmodel/tests/test_close_gaps.py @@ -36,7 +36,7 @@ class CloseGapsTests(unittest.TestCase): self.assertEqual(len(rmodel.gaps), 1) # obtain the scorer scorer = loop.SetupBackboneScorer(rmodel) - rmodel.CloseSmallDeletions(scorer) + rawmodel.CloseSmallDeletions(rmodel, scorer) self.assertEqual(len(rmodel.gaps), 0) self.assertEqual(self.log.messages['VERBOSE'], ['Assigning MOL_IDs', @@ -55,7 +55,7 @@ class CloseGapsTests(unittest.TestCase): self.assertEqual(len(rmodel.gaps), 2) self.assertEqual(str(rmodel.gaps[0]), 'A.ASP2-(FAGD)-A.THR7') self.assertEqual(str(rmodel.gaps[1]), 'A.THR7-(KNLG)-A.HIS12') - rmodel.MergeGapsByDistance(0) + rawmodel.MergeGapsByDistance(rmodel,0) self.assertEqual(len(rmodel.gaps), 1) self.assertEqual(str(rmodel.gaps[0]), 'A.ASP2-(FAGDTKNLG)-A.HIS12') self.assertEqual(self.log.messages['VERBOSE'], @@ -74,7 +74,7 @@ class CloseGapsTests(unittest.TestCase): aln.AttachView(1, tpl.CreateFullView()) rmodel = rawmodel.BuildRawModel(aln) self.assertEqual(len(rmodel.gaps), 2) - rmodel.MergeGapsByDistance(4) + rawmodel.MergeGapsByDistance(rmodel,4) self.assertEqual(len(rmodel.gaps), 2) def testMergeGapsByDistanceOneTerminal(self): @@ -85,7 +85,7 @@ class CloseGapsTests(unittest.TestCase): aln.AttachView(1, tpl.CreateFullView()) rmodel = rawmodel.BuildRawModel(aln) self.assertEqual(len(rmodel.gaps), 2) - rmodel.MergeGapsByDistance(2) + rawmodel.MergeGapsByDistance(rmodel,2) self.assertEqual(len(rmodel.gaps), 1) def testFillGapsByDatabase(self): @@ -99,9 +99,9 @@ class CloseGapsTests(unittest.TestCase): rmodel = rawmodel.BuildRawModel(aln) self.assertEqual(len(rmodel.gaps), 1) scorer = loop.SetupBackboneScorer(rmodel) - rmodel.FillLoopsByDatabase(scorer, loop.LoadFragDB(), - loop.LoadStructureDB(), - loop.LoadTorsionSamplerCoil()) + rawmodel.FillLoopsByDatabase(rmodel, scorer, loop.LoadFragDB(), + loop.LoadStructureDB(), + loop.LoadTorsionSamplerCoil()) self.assertEqual(len(rmodel.gaps), 0) self.assertEqual(self.log.messages['VERBOSE'], ['Assigning MOL_IDs', diff --git a/sidechain/doc/index.rst b/sidechain/doc/index.rst index d49a3f43d6685d3cc269b45e0158a52a58265582..1004293800d3edc4a0d399770c25618a21fe24c0 100644 --- a/sidechain/doc/index.rst +++ b/sidechain/doc/index.rst @@ -33,8 +33,6 @@ The simplest usage of the module is provided by the :func:`Reconstruct` function .. method:: Reconstruct(prot,[,keep_sidechains=False,build_disulfids,rotamer_model="frm",consider_hbonds=True,rotamer_library=None,add_polar_hydrogens=False]) - .. currentmodule:: promod3.sidechain - The function takes a structure and reconstructs its sidechains given the input parameters.