Add docstrings to deeprank/features/Edesolv.py. Add Edesolv to docs/d…

…eeprank.features.rst

deeprank/features/Edesolv.py

@@ -3,7 +3,6 @@
 from pdb2sql import pdb2sql, interface
 from deeprank.features import FeatureClass
 from copy import deepcopy
-from os import system
 from io import StringIO
 import numpy as np
 
@@ -16,11 +15,20 @@
 
 
 class Edesolv(FeatureClass):
-#class Edesolv():
 
     # init the class
-    def __init__(self,pdb_data):
-        #super().__init__('Atomic')
+    def __init__(self, pdb_data):
+        """Compute the energy of desolvation feature.
+
+        Biopython 1.79 or later is required for this feature.
+
+        Args:
+            pdb_data (list(byte) or str): pdb data
+
+        Example:
+            >>> edesolv = Edesolv('1AK4.pdb')
+            >>> edesolv.get_feature()
+        """
         self.pdb = pdb_data
         print('PDB DATA: ')
         print(pdb_data)
@@ -29,26 +37,27 @@ def __init__(self,pdb_data):
 
         self.feature_data = {}
         self.feature_data_xyz = {}
+
     # the feature extractor
-
     def get_feature(self, verbose=False):
-    #def __compute_feature__(self, verbose=False):
+        """ Computes the actual Edesolv feature"""
+
         esolcpx = 0.0
         esolfree = 0.0
+
         # create a sql database
         pdb_db = pdb2sql(self.pdb)
         self.db = interface(pdb_db)
+
+        #Create a temporary pdb-like string to calculate complex and free Edesolv
         temp_pdb = StringIO(('\n').join(pdb_db.sql2pdb()))
         temp_pdb.seek(0)
 
-        #Export PDB to parse it in BioPython
-        #pdb_db.exportpdb(fname = '_temp_.pdb')
         # Parse structure in Bio.PDB
         p = PDBParser(QUIET=1) 
-        #struct = p.get_structure('temp', '_temp_.pdb')#(self.pdb.split('/')[-1].split('.')[0], self.pdb) 
-        #system('rm _temp_.pdb')
         struct = p.get_structure('temp', temp_pdb)
 
+        # Identify chains in the Bio.PDB object
         chains = [x.id for x in struct[0] if x.id != ' ']
         chainA = chains[0]
         chainB = chains[1]
@@ -68,6 +77,7 @@ def get_feature(self, verbose=False):
         # get the contact atoms
         indA,indB = list(self.db.get_contact_atoms(chain1=chainA, chain2=chainB).values())
         contact = indA + indB
+
         # extract the atom keys and xyz of the contact CA atoms
         keys = self.db.get('serial,chainID,resName,resSeq,name',rowID=contact)
         xyz = self.db.get('x,y,z',rowID=contact)
@@ -85,106 +95,92 @@ def get_feature(self, verbose=False):
 
 
         # create the dictionary of human readable and xyz-val data
-        hread, xyzval = {},{}
         self.edesolv_data = {}
         self.edesolv_data_xyz = {}
 
         for key, coords in zip(keys, xyz):
             brk_flag = False
             disulf_bond = False
             atom = Atom(key)
-            #if verbose:
-            print('key: ', key)
-            #print('Values for atom n. ', atom.serial)
-            #print('PDB2SQL values:')
-            #print('Name: ', atom.name, 'Resn: ', atom.resn, 'Pos: ', atom.position) # atom.position specifies the position in the residue, if Back Bone or Side Chain
 
+            # Calculate complex sasa
             try:
                 atom.cpx_sasa = struct[0][atom.chainID][atom.resid][atom.name].sasa
             except KeyError: #Handling alternative residues error
                 print('Alternative residue found at:', key)
                 brk_flag = True
-                #raise Exception('KeyError')
                 pass 
+
+            # Calculate unbound (translated) sasa
             try:
                 atom.free_sasa = free_struct[0][atom.chainID][atom.resid][atom.name].sasa
             except KeyError: #Handling alternative residues error
                 print('Alternative residue found at:', key)
                 brk_flag = True
-                #raise Exception('KeyError')
                 pass 
 
             if not brk_flag:
 
-                #if verbose:
-                #print('Complex SASA: ', atom.cpx_sasa)
-                #print('Free    SASA: ', atom.free_sasa)
-
+                # check if the atom is a sulfide involved in a disulfide bond
                 for cys in disulfides_cys:
                     if cys == (atom.chainID, atom.resid):
                         disulf_bond = True
 
-                assign_solv_param(atom, disulf_bond)
-                #if verbose:
-                #print('Solv: ', atom.solv)
+                # Assign solvation parameters
+                Assign_solv_param(atom, disulf_bond)
 
+                #Calculate atom complex edesolv, free edesolv and final edesolv (cpx - free)
                 atom.esolcpx = atom.cpx_sasa * atom.solv
                 atom.esolfree = atom.free_sasa * atom.solv
                 atom.edesolv = atom.esolcpx - atom.esolfree
                 esolcpx += atom.esolcpx
                 esolfree += atom.esolfree
-                #if verbose:
-                #print('atom_esolcpx: ', atom.esolcpx)
-                #print('atom_esolfree: ', atom.esolfree)
-                #print('atom_Edesolv: ', atom.edesolv)
-                #print('\n')
-                #raise Exception('OK.')
-
-                # human readable
-                # { (chainID,resName,resSeq,name) : [val] }
-                hread[tuple(key)] = [atom.edesolv] #Put Edesolv here!
 
                 # xyz-val
                 # { (0|1,x,y,z) : [val] }
                 chain = [{chainA:0,chainB:1}[key[1]]]
-                #k = tuple(chain + xyz)
                 k = tuple(chain + coords)
-                #print(k, type(k))
-                #xyzval[k] = [atom.edesolv]
 
+                #Human readable
                 self.edesolv_data[(key[1], key[3], key[2], key[4])] = [atom.edesolv]
+                #XYZ data
                 self.edesolv_data_xyz[k] = [atom.edesolv]
 
-
-        #self.feature_data['Edesolv'] = hread
-        #self.feature_data_xyz['Edesolv'] = xyzval
-        print('edesolv_data: ', self.edesolv_data)
-        print('edesolv_data_xyz: ', self.edesolv_data_xyz)
         self.feature_data['Edesolv'] = self.edesolv_data
         self.feature_data_xyz['Edesolv'] = self.edesolv_data_xyz
 
-        #Edesolv = esolcpx - esolfree
-        #return Edesolv
-
 ##################################################################
 
 ##################################################################
 
-def __compute_feature__(pdb_data, featgrp, featgrp_raw, chainA, chainB):
+def __compute_feature__(pdb_data, featgrp, featgrp_raw):
+    """Main function called in deeprank for the feature calculations.
 
-        edesolv_feat = Edesolv(pdb_data)
-        edesolv_feat.get_feature()
+    Args:
+        pdb_data (list(bytes)): pdb information
+        featgrp (str): name of the group where to save xyz-val data
+        featgrp_raw (str): name of the group where to save human readable data
+    """
+
+    edesolv_feat = Edesolv(pdb_data)
+    edesolv_feat.get_feature()
 
-        # export in the hdf5 file
-        edesolv_feat.export_dataxyz_hdf5(featgrp)
-        edesolv_feat.export_data_hdf5(featgrp_raw)
+    # export in the hdf5 file
+    edesolv_feat.export_dataxyz_hdf5(featgrp)
+    edesolv_feat.export_data_hdf5(featgrp_raw)
 
-        # close
-        edesolv_feat.db._close()        
+    # close
+    edesolv_feat.db._close()        
 
 class Atom():
 
-    def __init__(self,atom_specs): # [498, 'A', 'HIS', 54, 'CB']
+    def __init__(self,atom_specs):
+        """Atom class necessary to calculate atomic-level surface area,
+        solvation parameters and desolvation energy.
+        
+        Args:
+            atom_specs("""
+
         serial = atom_specs[0]
         chainID = atom_specs[1]
         resn = atom_specs[2]
@@ -202,7 +198,18 @@ def __init__(self,atom_specs): # [498, 'A', 'HIS', 54, 'CB']
             self.position = 'SC'
 
 
-def assign_solv_param(atom, disulf_bond=False):
+def Assign_solv_param(atom, disulf_bond=False):
+    """Assigns solvation parameter to an Atom object
+
+    Args:
+        atom (Atom): Edesolv.Atom object.
+        disulf_bond (bool, optional): Disulfide bond flag. If true, the atom is a sulfide involved in a disulfie bond. Defaults to False.
+
+    Returns:
+        None.
+
+    """
+
     arofac = 6.26
     alifac = 1.27
     polfac = 2.30
@@ -236,7 +243,7 @@ def assign_solv_param(atom, disulf_bond=False):
         atom.solv = 0.0000
 
     '''
-    # Corarse-Grained model
+    # Corarse-Grained model from HADDOCK code for eventual future implementation
     if atom.position == "BB" and atom.resn=="ALA":
         atom.solv = -0.0107
     elif atom.position == "BB" and atom.resn=="GLY":
@@ -367,8 +374,14 @@ def get_disulfide_bonds(pdb_data):
     return disulfide_cys
 
 def get_dihedral(p):
-    """Praxeolitic formula
-    1 sqrt, 1 cross product"""
+    """Get dihedral angle from four points (atoms) 3D coordinates.
+    
+    Args:
+        p (list): list of four points 3D coordinates. Each element must be a list of floats.
+        
+    Returns:
+        dihedral (float): Dihedral angle.
+        """
     p0 = p[0]
     p1 = p[1]
     p2 = p[2]
@@ -394,4 +407,5 @@ def get_dihedral(p):
     # v and w may not be normalized but that's fine since tan is y/x
     x = np.dot(v, w)
     y = np.dot(np.cross(b1, v), w)
-    return np.degrees(np.arctan2(y, x))
+    dihedral = np.degrees(np.arctan2(y, x))
+    return dihedral

docs/deeprank.features.rst

@@ -35,6 +35,16 @@ Buried Surface Area
 
 .. autofunction:: deeprank.features.BSA.__compute_feature__
 
+Energy of desolvation
+-------------------
+
+.. automodule:: deeprank.features.Edesolv
+    :members:
+    :undoc-members:
+    :private-members:
+
+.. autofunction:: deeprank.features.Edesolv.__compute_feature__
+
 FullPSSM
 --------