diff --git a/modules/mol/alg/doc/molalg.rst b/modules/mol/alg/doc/molalg.rst index 7758f3ff290de4812b9635514cf0f7a2f51e2f0a..e64bdcfd798f882b6e70ae12432d37f1ed0b9699 100644 --- a/modules/mol/alg/doc/molalg.rst +++ b/modules/mol/alg/doc/molalg.rst @@ -41,8 +41,14 @@ to -1, the each frame is superposed to the previous frame. :returns: A newly created coord group containing the superposed frames. - - + +.. autofunction:: ParseAtomNames + +.. autofunction:: MatchResidueByNum + +.. autofunction:: Superpose + + Steric Clashes -------------------------------------------------------------------------------- diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt index a1880895f2096619147137d66cfa94c0480ce7f2..2fbf108cbe9133fcab98cd55b30fb99c74dabd2a 100644 --- a/modules/mol/alg/pymod/CMakeLists.txt +++ b/modules/mol/alg/pymod/CMakeLists.txt @@ -7,6 +7,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES set(OST_MOL_ALG_PYMOD_MODULES "__init__.py" views.py + superpose.py ) if (ENABLE_IMG) diff --git a/modules/mol/alg/pymod/__init__.py b/modules/mol/alg/pymod/__init__.py index 2e6198f819ed18fe8d559866e65c6229b1eb336d..0a4a26fa99bd6fc44c915b0ae544e99386a8a86a 100644 --- a/modules/mol/alg/pymod/__init__.py +++ b/modules/mol/alg/pymod/__init__.py @@ -1 +1,2 @@ from _mol_alg import * +from ost.mol.alg.superpose import * diff --git a/modules/mol/alg/pymod/superpose.py b/modules/mol/alg/pymod/superpose.py new file mode 100644 index 0000000000000000000000000000000000000000..44c337d8788d91e9daf33c32b6b14dd976addfe3 --- /dev/null +++ b/modules/mol/alg/pymod/superpose.py @@ -0,0 +1,157 @@ +""" +Superposition of structures made simple. + +Authors: Stefan Bienert +""" + +import math +import ost.mol.alg + + +def ParseAtomNames(atoms): + """ + Parses different representations of a list of atom names and returns a + :class:`set`, understandable by :func:`~ost.mol.alg.MatchResidueByNum`. In + essence, this function translates + + * None to ``None`` + + * 'all' to ``None`` + + * 'backbone' to ``set(['N', 'CA', 'C', 'O'])`` + + * 'aname1, aname2' to ``set(['aname1', 'aname2'])`` + + * ``['aname1', 'aname2']`` to ``set(['aname1', 'aname2'])`` + + :param atoms: Identifier or list of atoms + :type atoms: :class:`str`, :class:`list`, :class:`set` + """ + if atoms==None: + return None + if isinstance(atoms, str): + if atoms.upper()=='ALL': + return None + if atoms.upper()=='BACKBONE': + return set(['N', 'CA', 'C', 'O']) + ## if no recognised expression, split at ',' + return set([a.strip() for a in atoms.split(',')]) + return set(atoms) + + +def _EmptyView(view): + """ + for internal use, only + """ + if isinstance(view, ost.mol.EntityHandle): + return view.CreateEmptyView() + return view.handle.CreateEmptyView() + + +def _fetch_atoms(r_a, r_b, result_a, result_b, atmset): + """ + for internal use, only + """ + ## compare atoms of residues + for a_a in r_a.GetAtomList(): + if atmset==None or a_a.name in atmset: + a_b = r_b.FindAtom(a_a.name) + if a_b.IsValid(): + result_a.AddAtom(a_a) + result_b.AddAtom(a_b) + return result_a, result_b + + +def MatchResidueByNum(ent_a, ent_b, atoms='all'): + """ + Returns a tuple of views containing exactly the same number of atoms. + Residues are matched by residue number. A subset of atoms to be included in + the views can be specified in the **atoms** argument. Regardless of what the + list of **atoms** says, only those present in two matched residues will be + included in the views. Chains are processed in the order they occur in the + entities. If **ent_a** and **ent_b** contain a different number of chains, + processing stops with the lower count. + + :param ent_a: The first entity + :type ent_a: :class:`~ost.mol.EntityView` or :class:`~ost.mol.EntityHandle` + :param ent_b: The second entity + :type ent_b: :class:`~ost.mol.EntityView` or :class:`~ost.mol.EntityHandle` + :param atoms: The subset of atoms to be included in the two views. + :type atoms: :class:`str`, :class:`list`, :class:`set` + """ + ## init. final views + result_a=_EmptyView(ent_a) + result_b=_EmptyView(ent_b) + ## get lower no. of chains + if ent_a.chain_count < ent_b.chain_count: + n_chains=ent_a.chain_count + else: + n_chains=ent_b.chain_count + ## get a set of atoms or None + atmset=ParseAtomNames(atoms) + ## iterate chains + for i in range(0, n_chains): + chain_a=ent_a.chains[i] + chain_b=ent_b.chains[i] + ## residues of chains need to be consecutively numbered (sorted) + if chain_a.InSequence() and chain_b.InSequence(): + residues_a=iter(chain_a.residues) + residues_b=iter(chain_b.residues) + ## check residues & copy to views + try: + while True: + r_a=residues_a.next() + r_b=residues_b.next() + while r_a.number<r_b.number: + r_a=residues_a.next() + while r_b.number<r_a.number: + r_b=residues_b.next() + assert r_a.number==r_b.number + result_a,result_b=_fetch_atoms(r_a, r_b, result_a, result_b, atmset) + except StopIteration: + pass + else: + ## iterate one list of residues, search in other list + residues_a=iter(chain_a.residues) + try: + while True: + r_a=residues_a.next() + r_b=chain_b.FindResidue(r_a.number) + if r_b.IsValid(): + result_a,result_b=_fetch_atoms(r_a, r_b, result_a, result_b, atmset) + except StopIteration: + pass + result_a.AddAllInclusiveBonds() + result_b.AddAllInclusiveBonds() + return result_a, result_b + + +def Superpose(ent_a, ent_b, match='number', atoms='all'): + """ + Superposes the first entity onto the second. To do so, two views are created, + returned with the result. **atoms** describes what goes in to these views and + **match** the selection method. For superposition, + :func:`~ost.mol.alg.SuperposeSVD` is called. For matching, following methods + are recognised: + + * ``number`` - select residues by residue number, includes **atoms**, calls + :func:`~ost.mol.alg.MatchResidueByNum` + + :param ent_a: The first entity + :type ent_a: :class:`~ost.mol.EntityView` or :class:`~ost.mol.EntityHandle` + :param ent_b: The second entity + :type ent_b: :class:`~ost.mol.EntityView` or :class:`~ost.mol.EntityHandle` + :param match: Method to gather residues/ atoms + :type match: :class:`str` + :param atoms: The subset of atoms to be used in the superposition. + :type atoms: :class:`str`, :class:`list`, :class:`set` + """ + not_supported="Superpose called with unsupported matching request." + ## create views to superpose + if match.upper()=='NUMBER': + view_a, view_b=MatchResidueByNum(ent_a, ent_b, atoms) + else: + raise ValueError(not_supported) + ## action + res=ost.mol.alg.SuperposeSVD(view_a, view_b) + return res diff --git a/modules/mol/alg/src/svd_superpose.hh b/modules/mol/alg/src/svd_superpose.hh index 3580055cae58f4e1d4e5c5f4dce72fc2c1547e94..dad42504c8c026343b9164b2052519f1ad3927ff 100644 --- a/modules/mol/alg/src/svd_superpose.hh +++ b/modules/mol/alg/src/svd_superpose.hh @@ -51,7 +51,7 @@ struct DLLEXPORT_OST_MOL_ALG SuperpositionResult { class SuperposerSVDImpl; -/// \brief effiently superpose a bunch of models with the same number of atoms +/// \brief efficiently superpose a bunch of models with the same number of atoms /// Choose either two EntityViews or two AtomViewLists. class DLLEXPORT_OST_MOL_ALG SuperposerSVD { public: diff --git a/modules/mol/alg/tests/test_convenient_superpose.py b/modules/mol/alg/tests/test_convenient_superpose.py index bc051b8e1cc59d1079c06055b35af88c6b60a9bc..0b8e5bc4f0613f03f3af27fb92866c35c22a7e90 100644 --- a/modules/mol/alg/tests/test_convenient_superpose.py +++ b/modules/mol/alg/tests/test_convenient_superpose.py @@ -9,6 +9,7 @@ class TestConvenientSuperpose(unittest.TestCase): def runAtomOrdering(self, view1, view2): # call atom ordering function here + view1, view2 = ost.mol.alg.MatchResidueByNum(view1, view2) return view1, view2 def assertEqualAtomOrder(self, view1, view2): @@ -77,7 +78,7 @@ class TestConvenientSuperpose(unittest.TestCase): ent_view_wrong.AddChain(c) for r in c.residues: ent_view_wrong.AddResidue(r) - atoms = r.atoms + atoms = list(r.atoms) random.shuffle(atoms) for a in atoms: ent_view_wrong.AddAtom(a) @@ -86,17 +87,12 @@ class TestConvenientSuperpose(unittest.TestCase): view1, view2 = self.runAtomOrdering(ent_view_wrong, ent_view) self.assertEqualAtomOrder(view1, view2) - def testWrongResidueOrder(self): - ''' - there is a bug in ost (BZDNG-260) that gives weird atom counts here, thus - this test is skipped for the moment - - + def testWrongResidueOrder(self): ent_view = io.LoadEntity(os.path.join("testfiles","1aho.pdb")).Select("") ent_view_wrong = ent_view.CreateEmptyView() for c in ent_view.chains: ent_view_wrong.AddChain(c) - residues = c.residues + residues = list(c.residues) random.shuffle(residues) for r in residues: ent_view_wrong.AddResidue(r) @@ -106,8 +102,6 @@ class TestConvenientSuperpose(unittest.TestCase): self.assertEqualAtomOrder(view1, view2) view1, view2 = self.runAtomOrdering(ent_view_wrong, ent_view) self.assertEqualAtomOrder(view1, view2) - ''' - pass if __name__ == "__main__": - unittest.main() \ No newline at end of file + unittest.main()