dssp style secondary structure assignment

modules/mol/alg/doc/molalg.rst

@@ -920,6 +920,15 @@ Algorithms on Structures
   :return: The summed solvent accessibilty of each atom in *ent*.
 
 
+.. method:: AssignSecStruct(ent)
+
+  Assigns secondary structures to all residues based on hydrogen bond patterns
+  as described by DSSP.
+
+  :param ent:           Entity on which to assign secondary structures
+  :type ent:            :class:`~ost.mol.EntityView`/
+                        :class:`~ost.mol.EntityHandle`
+
 
 .. _traj-analysis:
 

modules/mol/alg/pymod/CMakeLists.txt

@@ -6,6 +6,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES
   export_structure_analysis.cc
   export_contact_overlap.cc
   export_accessibility.cc
+  export_sec_structure.cc
 )
 
 set(OST_MOL_ALG_PYMOD_MODULES

modules/mol/alg/pymod/export_sec_structure.cc

+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2017 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+
+
+#include <boost/python.hpp>
+
+#include <ost/mol/alg/sec_struct.hh>
+
+using namespace boost::python;
+
+namespace{
+
+void AssignSecStructView(ost::mol::EntityView& view) {
+  ost::mol::alg::AssignSecStruct(view);
+}
+
+void AssignSecStructHandle(ost::mol::EntityHandle& handle) {
+  ost::mol::alg::AssignSecStruct(handle);
+}
+
+} // ns
+
+void export_sec_struct() {
+  def("AssignSecStruct", &AssignSecStructView, (arg("ent")));
+  def("AssignSecStruct", &AssignSecStructHandle, (arg("ent")));
+}
+

modules/mol/alg/pymod/wrap_mol_alg.cc

@@ -41,6 +41,7 @@ void export_StructureAnalysis();
 void export_Clash();
 void export_contact_overlap();
 void export_accessibility();
+void export_sec_struct();
 #if OST_IMG_ENABLED
 void export_entity_to_density();
 #endif
@@ -111,6 +112,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
   export_Clash();
   export_contact_overlap();
   export_accessibility();
+  export_sec_struct();
   #if OST_IMG_ENABLED
   export_entity_to_density();
   #endif

modules/mol/alg/src/CMakeLists.txt

@@ -19,6 +19,7 @@ set(OST_MOL_ALG_HEADERS
   domain_find.hh
   similarity_matrix.hh
   accessibility.hh
+  sec_struct.hh
 )
 
 set(OST_MOL_ALG_SOURCES
@@ -41,6 +42,7 @@ set(OST_MOL_ALG_SOURCES
   domain_find.cc
   similarity_matrix.cc
   accessibility.cc
+  sec_struct.cc
 )
 
 set(MOL_ALG_DEPS ost_mol ost_seq)

modules/mol/alg/src/sec_struct.hh

+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2017 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+#ifndef OST_SEC_STRUCT_HH
+#define OST_SEC_STRUCT_HH
+
+#include <vector>
+
+#include <ost/base.hh>
+#include <ost/geom/geom.hh>
+#include <ost/message.hh>
+#include <ost/mol/mol.hh>
+
+namespace ost { namespace mol{ namespace alg{
+
+inline Real DSSPHBondEnergy(const geom::Vec3& h_pos, const geom::Vec3& n_pos,
+                            const geom::Vec3& c_pos, const geom::Vec3& o_pos) {
+  Real on = 1.0/geom::Distance(o_pos,n_pos);
+  Real ch = 1.0/geom::Distance(c_pos,h_pos);
+  Real oh = 1.0/geom::Distance(o_pos,h_pos);
+  Real cn = 1.0/geom::Distance(c_pos,n_pos);
+  return 27.888 * (on+ch-oh-cn);
+}
+
+// Raw estimation of secondary structure
+//
+// This function is not intended for Python export, since the input has to be 
+// prepared carefully. It basically estimates the secondary structure of a 
+// stretch of amino acids based on the hydrogen bond pattern as described for 
+// the dssp tool. 
+//
+// To define the hydrogen bonds you need to provide the vectors donor_for_one
+// and donor_for two. The index of an acceptor residue appears in one of the
+// two vectors if the corresponding hbond energy is < -0.5 
+// (be aware of prolines that can't be donors!).
+// There are two of those vectors, because for every residue we store the two
+// lowest energy acceptors. If there is no acceptor available, the value
+// at this position must be -1.
+// The connected_to_next contains zeros and ones that defines, whether there
+// is a peptide bond towards the next residue and the ca_positions vector
+// is self explaining.
+// As an additional feature you can also provide the according data for the
+// full structure but only estimate the secondary structure for a small
+// stretch, that gets defined by start_idx and size.
+// For an example usage have a look at the AssignSecStruct functions.
+
+String RawEstimateSS(const std::vector<geom::Vec3>& ca_positions, 
+                     int start_idx, int size,
+	                   const std::vector<int>& donor_for_one,
+                     const std::vector<int>& donor_for_two, 
+                     const std::vector<int>& connected_to_next);
+
+void PrepareSSData(const ost::mol::ResidueViewList& res_list,
+                 std::vector<int>& res_indices,
+                 std::vector<geom::Vec3>& ca_positions,
+                 std::vector<int>& donor_for_one,
+                 std::vector<int>& donor_for_two,
+                 std::vector<int>& connected_to_next);
+
+void AssignSecStruct(ost::mol::EntityView& ent);
+
+void AssignSecStruct(ost::mol::EntityHandle& ent);
+
+}}} //ns
+
+#endif

modules/mol/alg/tests/CMakeLists.txt

@@ -2,15 +2,16 @@ set(OST_MOL_ALG_UNIT_TESTS
   test_superposition.cc
   tests.cc
   test_consistency_checks.cc
+  test_partial_sec_struct_assignment.cc
   test_pdbize.py
   test_convenient_superpose.py
   test_hbond.py
   test_accessibility.py
+  test_sec_struct.py
 )
 
 if (COMPOUND_LIB)
   list(APPEND OST_MOL_ALG_UNIT_TESTS test_qsscoring.py)
 endif()
 
-ost_unittest(MODULE mol_alg SOURCES "${OST_MOL_ALG_UNIT_TESTS}")
-
+ost_unittest(MODULE mol_alg SOURCES "${OST_MOL_ALG_UNIT_TESTS}" LINK ost_io)

modules/mol/alg/tests/test_partial_sec_struct_assignment.cc

+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2017 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+
+
+#include <ost/mol/alg/sec_struct.hh>
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+#include <boost/test/auto_unit_test.hpp>
+#include <ost/mol/mol.hh>
+#include <ost/io/mol/pdb_reader.hh>
+#include <ost/mol/builder.hh>
+#include <ost/message.hh>
+
+
+
+BOOST_AUTO_TEST_SUITE( mol_alg );
+
+BOOST_AUTO_TEST_CASE(partial_sec_struct_assignment) {
+
+  // load a test structure
+  String filepath = "testfiles/1a0s.pdb";
+  ost::io::PDBReader reader(filepath, ost::io::IOProfile());
+  ost::mol::EntityHandle test_ent = ost::mol::CreateEntity();
+  reader.Import(test_ent);
+
+  ost::mol::EntityView view = test_ent.CreateFullView();
+  ost::mol::ResidueViewList res_list = view.GetResidueList();
+
+  std::vector<int> res_indices;
+  std::vector<geom::Vec3> ca_positions;
+  std::vector<int> donor_for_one;
+  std::vector<int> donor_for_two;
+  std::vector<int> connected_to_next;
+
+  ost::mol::alg::PrepareSSData(res_list, res_indices, ca_positions,
+                               donor_for_one, donor_for_two,
+                               connected_to_next);
+
+  String full_pred = ost::mol::alg::RawEstimateSS(ca_positions,
+                                                  0, ca_positions.size(),
+                                                  donor_for_one, donor_for_two,
+                                                  connected_to_next);
+
+  // select some start indices to test
+  int max_num = ca_positions.size();
+  std::vector<int> test_locations;
+  test_locations.push_back(max_num/10);
+  test_locations.push_back((max_num/10)*2);
+  test_locations.push_back((max_num/10)*3);
+  test_locations.push_back((max_num/10)*4);
+  test_locations.push_back((max_num/10)*5);
+  test_locations.push_back((max_num/10)*6);
+  test_locations.push_back((max_num/10)*7);
+  test_locations.push_back((max_num/10)*8);
+  test_locations.push_back((max_num/10)*9);
+  int stretch_length = std::min(12, max_num - test_locations.back());
+
+  for(uint i = 0; i < test_locations.size(); ++i) {
+    String expected_ss = full_pred.substr(test_locations[i], stretch_length);
+    String calculated_ss = ost::mol::alg::RawEstimateSS(ca_positions,
+                                                        test_locations[i],
+                                                        stretch_length,
+                                                        donor_for_one,
+                                                        donor_for_two,
+                                                        connected_to_next);
+
+    BOOST_CHECK(expected_ss == calculated_ss);
+  }
+
+}
+
+BOOST_AUTO_TEST_SUITE_END();

modules/mol/alg/tests/test_sec_struct.py

+from ost import io, mol, settings
+import unittest
+import os
+from ost.bindings import dssp
+
+
+class TestSecStruct(unittest.TestCase):
+  
+  def testSecStruct(self):
+
+    # unit test only makes sense, when a dssp binary is around
+    dssp_path = None
+    try:
+      dssp_path = settings.locate("dssp")
+    except:
+      try:
+        dssp_path = settings.locate("mkdssp")
+      except:
+        print "Could not find dssp, could not compare sec struct assignment..."
+
+    dssp_ent = io.LoadPDB(os.path.join("testfiles", "1a0s.pdb"))
+    ost_ent = io.LoadPDB(os.path.join("testfiles", "1a0s.pdb"))
+
+    dssp.AssignDSSP(dssp_ent, dssp_bin = dssp_path)
+    mol.alg.AssignSecStruct(ost_ent)
+
+    for a, b in zip(dssp_ent.residues, ost_ent.residues):
+      self.assertTrue(str(a.GetSecStructure()) == str(b.GetSecStructure()))
+      
+
+if __name__ == "__main__":
+
+  from ost import testutils
+  testutils.RunTests()