Allow to control the property string to extract solvent accessibilities when adding to StructureDB

doc/tests/scripts/loop_structure_db.py

@@ -31,8 +31,9 @@ for i in ids:
     prof_path = os.path.join(prof_dir, pdb_id + chain_id + ".hhm")
 
     # let's load the structure
-    structure = io.LoadPDB(structure_path)\
-                  .Select("cname=" + chain_id + " and peptide=True")
+    structure = io.LoadPDB(structure_path).Select("peptide=True")
+    single_chain_structure = structure.Select("cname=" + chain_id)
+    
     # and the according profile in hhm format
     prof = io.LoadSequenceProfile(prof_path)
 
@@ -43,9 +44,12 @@ for i in ids:
     # as a final step we need the surface as calculated by msms
     surf = msms.CalculateSurface(structure)[0]
 
-    # let's add it
+    # let's add it! Please note, that we calculated solvent 
+    # accessiblity / secondary structure and surface on the full
+    # structure but only add one of the chains.
     structure_db.AddCoordinates(pdb_id, chain_id,
-                                structure.chains[0], surf, prof)
+                                single_chain_structure.chains[0], 
+                                surf, prof)
 
 # we now have two structures in the database...
 

loop/doc/structure_db.rst

@@ -6,15 +6,14 @@ Structural Data
 The structural database serves as a container for structural backbone
 and profile data. It can be filled with chains of pdb structures with their
 corresponding profiles as they are produced by the HHSuite tools [soding2005]_. 
-Some of the extracted features expect the secondary structure, as well as the solvent
-accessibility to be assigned using dssp [kabsch1983]_.
+Structural and profile data get complemented by with additional information.
 Following features get stored on a per residue basis:
 
 * The amino acid one letter code
 * The coordinates of the backbone atoms (N,CA,C,O)
 * The phi/psi dihedral angles
 * The secondary structure state as defined by dssp
-* The solvent accessibility in square Angstrom as calculated by dssp
+* The solvent accessibility in square Angstrom 
 * The residue depth defined as the average distance from all atoms of a
   residue to the closest surface vertex as calculated by msms [sanner1996]_.
   This is a simplified version of the residue depth as discussed in
@@ -132,18 +131,29 @@ database, you might want to consider two things:
               db has been loaded with load_frequencies=False (enforces only
               complete databases to be saved down)
 
-  .. method:: AddCoordinates(pdb_id, chain_name, chain, surf, prof)
+  .. method:: AddCoordinates(pdb_id, chain_name, chain, surf, prof, \
+                             solvent_accessibility_string="solvent_accessibility")
 
     This method takes a structural chain and searches the longest stretch of
     connected residues containing all necessary backbone atoms. This stretch
     gets then added to the database. Due to technical reasons, The maximal
     extent along one of the base axis cannot exceed 650 A.
+    Following features are expected to be set on a per residue level in *chain*:
+    The secondary structure (e.g. call :meth:`ost.mol.alg.AssignSecStruct` on the
+    full :class:`ost.mol.EntityView` the *chain* belongs to) and
+    the solvent accessibility assigned using a float property with name 
+    *solvent_accessibility_string* as name on a per residue level (e.g. call 
+    :meth:`ost.mol.alg.Accessibility` on the full :class:`ost.mol.EntityView`
+    the chain belongs to and take care to pass the appropriate
+    *solvent_accessibility_string*).
 
     :param pdb_id:      4-letter code of the structure the chain belongs to
     :param chain_name:  Name of the chain consisting of one letter
     :param chain:       The actual chain 
     :param surf:        A surface describing the solvent accessible surface
-    :param prof:        Profile information for this structure.
+                        (we advise you to provide a surface estimated with the
+                        full :class:`ost.mol.EntityView` the *chain* belongs to)
+    :param prof:        Profile information for this *chain*.
 
     :type pdb_id:       :class:`str`
     :type chain_name:   :class:`str`

loop/pymod/export_structure_db.cc

@@ -87,8 +87,9 @@ void export_StructureDB(){
     .staticmethod("LoadPortable")
     .def("SavePortable", &StructureDB::SavePortable, (arg("filename")))
     .def("AddCoordinates", &StructureDB::AddCoordinates,
-         (arg("pdb_id"), arg("chain_name"), arg("chain"), arg("surf"),
-          arg("prof")))
+         (arg("pdb_id"), arg("chain_name"), 
+          arg("chain"), arg("surf"), arg("prof"), 
+          arg("solvent_accessibility_string")="solvent_accessibility"))
     .def("GetCoordIndex", &StructureDB::GetCoordIndex,
          (arg("pdb_id"), arg("chain_name")))
     .def("GetBackboneList", &wrap_get_bb_list_one,

loop/src/structure_db.cc

@@ -307,7 +307,8 @@ int StructureDB::AddCoordinates(const String& pdb_id,
                                 const String& chain_name,
                                 const ost::mol::ChainView& chain,
                                 const ost::mol::SurfaceHandle& surf,
-                                ost::seq::ProfileHandlePtr prof) {
+                                ost::seq::ProfileHandlePtr prof,
+                                const String& solvent_accessibility_string) {
 
   if(!frequencies_present_){
     std::stringstream ss;
@@ -416,14 +417,15 @@ int StructureDB::AddCoordinates(const String& pdb_id,
       dssp_state = 'C';
     }
     dssp_states.push_back(dssp_state);
-    if(i->HasProp("solvent_accessibility")){
-      solvent_accessibilities.push_back(i->GetFloatProp("solvent_accessibility"));
+    if(i->HasProp(solvent_accessibility_string)){
+      Real acc = i->GetFloatProp(solvent_accessibility_string);
+      solvent_accessibilities.push_back(acc);
     }
     else{
-      std::stringstream ss;
-      ss << "Residue \"" << i->GetQualifiedName() <<"\" has no property \"solvent_accessibility\"";
-      ss << "adding 0.0 instead!";
-      std::cerr<<ss.str()<<std::endl;
+      String err = "Residue \"" + i->GetQualifiedName() + "\" has no ";
+      err += ("property " + solvent_accessibility_string + " describing ");
+      err += "its solvent accessibility. Adding 0.0 instead.";
+      std::cerr<<err<<std::endl;
       solvent_accessibilities.push_back(0.0);
     }
   }

loop/src/structure_db.hh

@@ -260,7 +260,9 @@ public:
   int AddCoordinates(const String& pdb_id, const String& chain_name,
                      const ost::mol::ChainView& chain,
                      const ost::mol::SurfaceHandle& surf,
-                     ost::seq::ProfileHandlePtr prof);
+                     ost::seq::ProfileHandlePtr prof,
+                     const String& solvent_accessibilty_string = 
+                     "solvent_accessibility");
 
   int GetCoordIndex(const String& pdb_id, const String& chain_name) const;