diff --git a/modules/mol/base/pymod/export_editors.cc b/modules/mol/base/pymod/export_editors.cc
index aae903f7d276721545a60af51cd1ed0dff227263..68f48cb910663590c2112408f4203891bedbaf19 100644
--- a/modules/mol/base/pymod/export_editors.cc
+++ b/modules/mol/base/pymod/export_editors.cc
@@ -248,7 +248,7 @@ void export_Editors()
.def("AddTorsion", &EditorBase::AddTorsion)
.def("ReorderResidues",&EditorBase::ReorderResidues)
.def("ReorderAllResidues",&EditorBase::ReorderAllResidues)
- .def("RenumberAllResidues",&EditorBase::RenumberAllResidues)
+ .def("RenumberChain",&EditorBase::RenumberChain)
;
void (XCSEditor::*apply_transform1)(const geom::Mat4&) = &XCSEditor::ApplyTransform;
diff --git a/modules/mol/base/pymod/export_residue.cc b/modules/mol/base/pymod/export_residue.cc
index 284c2f0bc294825e5751a3f441e377d1d1666517..8f0e169944f856f6f47a0e6e37cf2068ddad9177 100644
--- a/modules/mol/base/pymod/export_residue.cc
+++ b/modules/mol/base/pymod/export_residue.cc
@@ -122,6 +122,10 @@ void export_Residue()
.def(self+int())
.def(self-int())
;
+ class_<ResNumList>("ResNumList", init<>())
+ .def(vector_indexing_suite<ResNumList>())
+ .def(geom::VectorAdditions<ResNumList>())
+ ;
implicitly_convertible<int, ResNum>();
scope().attr("PEPTIDE_LINKING")=char(ChemClass::PEPTIDE_LINKING);
diff --git a/modules/mol/base/src/editor_base.cc b/modules/mol/base/src/editor_base.cc
index be1b3f6d18b4ca84ad603a9a29867de6b32065a6..4d4111c8014b9934fa3321f3984cff606465daf1 100644
--- a/modules/mol/base/src/editor_base.cc
+++ b/modules/mol/base/src/editor_base.cc
@@ -245,6 +245,12 @@ BondHandle EditorBase::Connect(const AtomHandle& first,
return this->Connect(first,second,0.0, 0.0, 0.0, 1);
}
+void EditorBase::RenumberChain(ChainHandle chain, const ResNumList& new_numbers)
+{
+ CheckHandleValidity(chain);
+ chain.Impl()->RenumberAllResidues(new_numbers);
+}
+
BondHandle EditorBase::Connect(const AtomHandle& first,
const AtomHandle& second,
Real len, Real theta, Real phi,
diff --git a/modules/mol/base/src/editor_base.hh b/modules/mol/base/src/editor_base.hh
index d1674a68a9d2cde585b04e07c8ac2c5e6d58b1f8..5e1fbcdef596675758332fcc0477dfc69f85c6e7 100644
--- a/modules/mol/base/src/editor_base.hh
+++ b/modules/mol/base/src/editor_base.hh
@@ -216,7 +216,7 @@ public:
void RenameResidue(ResidueHandle res, const String& new_name);
void SetResidueNumber(ResidueHandle res, const ResNum& num);
-
+
void RenameChain(ChainHandle chain, const String& new_name);
/// \brief Assign type of chain according to ChainType.
@@ -281,6 +281,9 @@ public:
/// If set to false, residues will continously be renumbered ongoing from start.
/// Otherwise the spacings between the residues are kept.
void RenumberAllResidues(int start, bool keep_spacing);
+
+
+ void RenumberChain(ChainHandle chain, const ResNumList& new_numbers);
/// \brief Get edit mode of editor
EditMode GetMode() const {return mode_;}
diff --git a/modules/mol/base/src/impl/chain_impl.cc b/modules/mol/base/src/impl/chain_impl.cc
index 8b3cbe271eef0e860d78ea5039ac7f00a35e191e..2ec5b8694668d91c207679a99de9c0b2a80157f9 100644
--- a/modules/mol/base/src/impl/chain_impl.cc
+++ b/modules/mol/base/src/impl/chain_impl.cc
@@ -516,6 +516,17 @@ void ChainImpl::RenumberAllResidues(int start, bool keep_spacing)
UpdateShifts();
}
+void ChainImpl::RenumberAllResidues(const ResNumList& new_numbers)
+{
+ if (new_numbers.size() != residue_list_.size()) {
+ throw Error("number of residues and residue numbers must match");
+ }
+ for (size_t i = 0; i<new_numbers.size(); ++i) {
+ residue_list_[i]->SetNumber(new_numbers[i]);
+ }
+ this->UpdateShifts();
+}
+
void ChainImpl::SetInSequence(const int index)
{
ResNum num=residue_list_[index]->GetNumber();
diff --git a/modules/mol/base/src/impl/chain_impl.hh b/modules/mol/base/src/impl/chain_impl.hh
index a2c607ec4e2df3cca3783a8a14532e795536593e..82334c6fbb054eb200db5dacaec0587c9d5511d1 100644
--- a/modules/mol/base/src/impl/chain_impl.hh
+++ b/modules/mol/base/src/impl/chain_impl.hh
@@ -175,6 +175,8 @@ public:
void RenumberAllResidues(int start, bool keep_spacing);
+ void RenumberAllResidues(const ResNumList& new_numbers);
+
int GetIndex(const ResidueImplPtr& res) const;
void AssignSecondaryStructure(SecStructure ss,
const ResNum& start,
diff --git a/modules/mol/base/src/residue_prop.hh b/modules/mol/base/src/residue_prop.hh
index 4a411cca6e27250fdd23d4b538ffea6064b8d106..d4d79fd688ed96cf6f197f13e4116270e84ff212 100644
--- a/modules/mol/base/src/residue_prop.hh
+++ b/modules/mol/base/src/residue_prop.hh
@@ -19,6 +19,7 @@
#ifndef OST_RESIDUE_PROP_HH
#define OST_RESIDUE_PROP_HH
+#include <vector>
#include <boost/operators.hpp>
#include <ost/mol/module_config.hh>
@@ -122,6 +123,7 @@ private:
};
typedef String ResidueKey;
+typedef std::vector<ResNum> ResNumList;
inline std::ostream& operator<<(std::ostream& os, const ResNum& n)
{
diff --git a/modules/seq/alg/pymod/renumber.py b/modules/seq/alg/pymod/renumber.py
index 092850f2a356c4b240f3afe9f42ebab730c60b70..22873111e877c66da19ee2af80f49ce169ccd4ef 100644
--- a/modules/seq/alg/pymod/renumber.py
+++ b/modules/seq/alg/pymod/renumber.py
@@ -1,6 +1,43 @@
from ost import io, seq, mol, conop
from ost import *
+
+def _RenumberSeq(seq_handle):
+ if not seq_handle.HasAttachedView():
+ raise RuntimeError("Sequence Handle has no attached view")
+ ev = seq_handle.attached_view.CreateEmptyView()
+ new_numbers = mol.ResNumList()
+ for pos in range(len(seq_handle)):
+ if seq_handle[pos]!='-':
+ r=seq_handle.GetResidue(pos)
+ if r.IsValid():
+ #print seq_handle[pos],r.number.num,pos+1
+ ev.AddResidue(r, mol.INCLUDE_ALL)
+ new_numbers.append(pos+1)
+ else:
+ raise RuntimeError('Error: renumbering failed at position %s' %pos)
+ return ev, new_numbers
+
+def _RenumberAln(aln, seq_index):
+ if not aln.sequences[seq_index].HasAttachedView():
+ raise RuntimeError("Sequence Handle has no attached view")
+ counter=0
+ ev = aln.sequences[seq_index].attached_view.CreateEmptyView()
+ new_numbers = mol.ResNumList()
+ for col in aln:
+ if col[0]!='-' and col[seq_index]!='-':
+ if col[0]!=col[seq_index]:
+ raise RuntimeError("residue mismatch at position %d (%s vs %s) (renumbering failed)"%(counter, col[0],col[1]))
+ rnum=aln.GetSequence(seq_index).GetResidueIndex(counter)
+ r=aln.GetSequence(seq_index).GetResidue(counter)
+ if not r.IsValid():
+ raise RuntimeError("invalid residue at postion %s (renumbering failed)" %(counter))
+ ev.AddResidue(r, mol.INCLUDE_ALL)
+ new_numbers.append(counter+1)
+ counter+=1
+ return ev, new_numbers
+
+
def Renumber(seq_handle, sequence_number_with_attached_view=1):
"""
Function to renumber an entity according to an alignment between the model sequence
@@ -21,68 +58,11 @@ def Renumber(seq_handle, sequence_number_with_attached_view=1):
"""
if isinstance(seq_handle, seq.SequenceHandle):
- if seq_handle.HasAttachedView()==False:
- raise RuntimeError, "Sequence Handle has no attached view"
- changed_residue_count=0
- renumberingFlag = False
- ent_n=mol.CreateEntity()
- ed=ent_n.EditXCS()
- c=ed.InsertChain(" ")
- for pos in range(len(seq_handle)):
- if seq_handle[pos]!='-':
- r=seq_handle.GetResidue(pos)
- if r.IsValid():
- #print seq_handle[pos],r.number.num,pos+1
- if r.number.num!=pos+1:
- changed_residue_count+=1
- renumberingFlag = True
- r_n=ed.AppendResidue(c,r.name, mol.ResNum(pos+1))
- for atom in r.atoms:
- ed.InsertAtom(r_n,atom.name,atom.pos,element=atom.element)
- else:
- err='Error: renumbering failed at position %s' %pos
- raise RuntimeError, err
- if renumberingFlag == True:
- err = 'Warning: %s residues have been renumbered!' %changed_residue_count
- LogInfo(err)
-
- # FIXME: BZDNG-416
- #conop.ConnectAll(ent_n)
- return ent_n
-
+ ev, new_numbers = _RenumberSeq(seq_handle)
elif isinstance(seq_handle, seq.AlignmentHandle):
- if seq_handle.GetSequence(sequence_number_with_attached_view).HasAttachedView()==False:
- raise RuntimeError, "Sequence Handle has no attached view"
- dir(seq_handle)
- counter=0
- changed_residue_count=0
- renumberingFlag = False
- ent_n=mol.CreateEntity()
- ed=ent_n.EditXCS()
- c=ed.InsertChain(seq_handle.GetSequence(sequence_number_with_attached_view).GetAttachedView().chains[0].name)
- for col in seq_handle:
- if col[0]!='-' and col[1]!='-':
- if col[0]==col[1]:
- rnum=seq_handle.GetSequence(sequence_number_with_attached_view).GetResidueIndex(counter)
- r=seq_handle.GetSequence(sequence_number_with_attached_view).GetResidue(counter)
- if r.IsValid():
- if r.number.num!=counter+1:
- changed_residue_count+=1
- renumberingFlag = True
- r_n=ed.AppendResidue(c,r.name, mol.ResNum(counter+1))
- for atom in r.atoms:
- ed.InsertAtom(r_n,atom.name,atom.pos,atom.element, atom.b_factor,
- atom.occupancy, atom.is_hetatom)
-
- else:
- raise RuntimeError("invalide residue at postion %s (renumbering failed)" %(counter))
- else:
- raise RuntimeError("residue mismatch at position %d (%s vs %s) (renumbering failed)"%(counter, col[0],col[1]))
- counter+=1
- if renumberingFlag == True:
- err = 'Warning: %s residues have been renumbered!' %changed_residue_count
- LogInfo(err)
- # FIXME: BZDNG-416
- #conop.ConnectAll(ent_n)
- return ent_n
+ ev, new_numbers = _RenumberAln(seq_handle, sequence_number_with_attached_view)
+ ev.AddAllInclusiveBonds()
+ new_ent = mol.CreateEntityFromView(ev,False);
+ new_ent.EditXCS().RenumberChain(new_ent.chains[0], new_numbers)
+ return new_ent
diff --git a/modules/seq/alg/tests/test_renumber.py b/modules/seq/alg/tests/test_renumber.py
index 51e4727d3cfb7b3f87fea3517f9664fe7c896b16..4a39420af9a938c8ab3a3e44d47bcca76403711a 100644
--- a/modules/seq/alg/tests/test_renumber.py
+++ b/modules/seq/alg/tests/test_renumber.py
@@ -5,6 +5,11 @@ from ost import seq
from ost.bindings.clustalw import *
from ost.seq.alg import renumber
+try:
+ clustalw_path=settings.Locate(('clustalw', 'clustalw2'))
+except(settings.FileNotFound):
+ clustalw_path=None
+
class TestRenumber(unittest.TestCase):
def setUp(self):
@@ -16,14 +21,53 @@ class TestRenumber(unittest.TestCase):
self.peptide_del_4 = io.LoadEntity("testfiles/peptide_del_4.pdb")
self.peptide_mutation_3 = io.LoadEntity("testfiles/peptide_mutation_3.pdb")
+ def testRenumbersChainsBasedOnSequence(self):
+ aln_seq = seq.CreateSequence('A', 'MP-T---NA')
+ aln_seq.AttachView(self.peptide_original.Select(''))
+
+ renumbered = renumber.Renumber(aln_seq)
+
+ res = renumbered.residues
+ self.assertEqual(len(res), 5)
+ self.assertEqual(res[0].number, 1)
+ self.assertEqual(res[1].number, 2)
+ self.assertEqual(res[2].number, 4)
+ self.assertEqual(res[3].number, 8)
+ self.assertEqual(res[4].number, 9)
+
+ def testRenumbersChainsBasedOnAlignment(self):
+ aln_seq = seq.CreateSequence('A', 'MP-T---NA')
+ aln_seq.AttachView(self.peptide_original.Select(''))
+ aln = seq.CreateAlignment(seq.CreateSequence('A', 'MP-T-XXNA'), aln_seq)
+
+ renumbered = renumber.Renumber(aln)
+
+ res = renumbered.residues
+ self.assertEqual(len(res), 5)
+ self.assertEqual(res[0].number, 1)
+ self.assertEqual(res[1].number, 2)
+ self.assertEqual(res[2].number, 4)
+ self.assertEqual(res[3].number, 8)
+ self.assertEqual(res[4].number, 9)
+
+ def testRenumberPreservesBonds(self):
+ aln_seq = seq.CreateSequence('A', 'MP-T---NA')
+ aln_seq.AttachView(self.peptide_original.Select(''))
+
+ renumbered = renumber.Renumber(aln_seq)
+
+ self.assertTrue(mol.BondExists(renumbered.chains[0].FindAtom(1, 'N'),
+ renumbered.chains[0].FindAtom(1, 'CA')))
- def testPeptidePlusFive(self):
+ def testClustalWPeptidePlusFive(self):
"""
All residue numbers shifted by 5.
Check whether internal atom order changes while renumbering
(a new entity is generated in the edit_mode)
TODO: add more basic tests: are all properties preserved?
"""
+ if not clustalw_path:
+ return
model_seq=seq.SequenceFromChain(" ", self.peptide_plus_5.chains[0])
model_seq.name="model"
aln=ClustalW(self.target_seq,model_seq)
@@ -38,10 +82,12 @@ class TestRenumber(unittest.TestCase):
"Renumbering failed on atom level: restoring from ResNum+5"
- def testPeptideRandom(self):
+ def testClustalWPeptideRandom(self):
"""
Change residue names in random order
"""
+ if not clustalw_path:
+ return
model_seq=seq.SequenceFromChain(" ", self.peptide_random.chains[0])
model_seq.name="model"
aln=ClustalW(self.target_seq,model_seq)
@@ -56,10 +102,12 @@ class TestRenumber(unittest.TestCase):
"Renumbering failed on atom level: restoring from random residue numbers"
- def testPeptideDel_1_2(self):
+ def testClustalWPeptideDel_1_2(self):
"""
First two residues were removed
"""
+ if not clustalw_path:
+ return
model_seq=seq.SequenceFromChain(" ", self.peptide_del_1_2.chains[0])
model_seq.name="model"
aln=ClustalW(self.target_seq,model_seq)
@@ -76,10 +124,12 @@ class TestRenumber(unittest.TestCase):
"Renumbering failed on atom level: restoring from random residue numbers"
- def testPeptideDel_4(self):
+ def testClustalWPeptideDel_4(self):
"""
Residues in the middle (position 4) was removed
"""
+ if not clustalw_path:
+ return
model_seq=seq.SequenceFromChain(" ", self.peptide_del_4.chains[0])
model_seq.name="model"
aln=ClustalW(self.target_seq,model_seq)
@@ -98,10 +148,12 @@ class TestRenumber(unittest.TestCase):
"Renumbering failed on atom level: restoring from random residue numbers"
- def testPeptideMutation_3(self):
+ def testClustalWPeptideMutation_3(self):
"""
Mutation to GLY at postion 3
"""
+ if not clustalw_path:
+ return
model_seq=seq.SequenceFromChain(" ", self.peptide_mutation_3.chains[0])
model_seq.name="model"
aln=ClustalW(self.target_seq,model_seq)
@@ -113,10 +165,7 @@ class TestRenumber(unittest.TestCase):
if __name__ == "__main__":
# test renumbering
# test if clustalw package is available on system, otherwise ignore tests
- try:
- clustalw_path=settings.Locate(('clustalw', 'clustalw2'))
- except(settings.FileNotFound):
- print "Could not find clustalw executable: ignoring unit tests"
- sys.exit(0)
+ if not clustalw_path:
+ print "Could not find clustalw executable: ignoring some renumber unit tests"
from ost import testutils
testutils.RunTests()
diff --git a/modules/seq/base/src/sequence_handle.hh b/modules/seq/base/src/sequence_handle.hh
index 937387c797bc42f28b0fe902f1e9af256d6fa416..d7247cbd10abca25bdc9fc3c00b60dac7c6e9915 100644
--- a/modules/seq/base/src/sequence_handle.hh
+++ b/modules/seq/base/src/sequence_handle.hh
@@ -284,7 +284,7 @@ public:
///
/// The sequence is mapped onto the chain with given name
void AttachView(const mol::EntityView& view, const String& chain_name);
-
+
/// \internal
SequenceHandle(const impl::SequenceImplPtr& impl);