Finish docstrings, update readme

README.md

@@ -24,11 +24,60 @@ Load your environnement
 conda activate dssp
 ```
 
+### Running
+
 You should be ready to run the program, by calling, for example
 ```
 python3 src/dssp.py data/1est.pdb
 ```
 
+And if you want the extended output, use the -v flag, or --verbose
+```
+python3 src/dssp.py data/1est.pdb -v
+```
+
+To get usage help, run
+```
+python3 src/dssp.py data/1est.pdb -h
+```
+
+## Output
+
+To understand the output, here is a short description of the DSSP output, as followed from the official website : https://swift.cmbi.umcn.nl/gv/dssp/
+
+*   RESIDUE
+
+two columns of residue numbers. First column is DSSP's sequential residue number, starting at the first residue actually in the data set and including chain breaks; this number is used to refer to residues throughout. Second column gives crystallographers' 'residue sequence number','insertion code' and 'chain identifier' (see protein data bank file record format manual), given for reference only. 
+
+*   AA
+
+one letter amino acid code
+
+
+*   BP1 
+
+residue number of first bridge partner
+
+*   TCO
+
+cosine of angle between C=O of residue I and C=O of residue I-1. For α-helices, TCO is near +1, for β-sheets TCO is near -1. Not used for structure definition.
+
+*   KAPPA
+
+virtual bond angle (bend angle) defined by the three Cα atoms of residues I-2,I,I+2. Used to define bend (structure code 'S').
+
+*   ALPHA
+
+virtual torsion angle (dihedral angle) defined by the four Cα atoms of residues I-1,I,I+1,I+2.Used to define chirality (structure code '+' or '-').
+
+*   PHI PSI
+
+IUPAC peptide backbone torsion angles
+
+*   X-CA Y-CA Z-CA
+
+echo of Cα atom coordinates 
+
 ## Author
 
 * **FOREST Thomas** - *M2BI* - [Univ-Paris-Diderot.fr](https://www.univ-paris-diderot.fr/) - thomas.forest@etu.univ-paris-diderot.fr

doc/notes.org

@@ -1,37 +0,0 @@
-Dépôt output DSSP:
-http://ftp.cbi.pku.edu.cn/pub/database/DSSP/latest_release/
-
-
-structure basée sur liaisons H
-
-- n-turns : liaison H entre CO et NH du res i+n, ou n = [3,4,5]
-- bridges : LH avec res non proches
-
-* hélices : 4-turns répétés
-* brins beta : repetition bridges
-
-common imperfections : helical kinks + beta-bulges
-
-propriétés déduites geométriquement : bends, chirality, SS bonds, solvent exposure
-
-
-Hbond(i,j) = [E<-0.5kcal/mole]
-
-electrostatic interaction energy between twho H-bonding groups by placing partial charges on the C,O(+q1,-q1) and N,H(-q2,+q2) atoms 
-
-E = q1q2(1/r(ON) + 1/r(CH) - 1/r(OH) - 1/r(CN))*f
-avec
-q1 = 0.42e
-q2 = 0.20e
-e la charge élémentaire d'un electro (e = 1,602 176 634 ×10−19 C)
-r en angström
-f (facteur de dimmensionnalité) = 332
-E en kcal/mol
-
-Bonne LH si E <= -3 kcal/mol
-
-Cutoff choisi permet une distance N-O jusqu'à 2.2 Angström
-N-H bound distance : 1.030 A
-https://pubs.acs.org/doi/pdf/10.1021/ja00313a047
-
-Angles : https://www.sciencedirect.com/topics/biochemistry-genetics-and-molecular-biology/peptide-bond

src/chem.py

@@ -2,8 +2,30 @@ import math
 from geom import *
 
 class Atom:
+    """Defines Atom objects and their associated methods.
 
+    This function is used to instenciate Atom objects with
+    convenient fields, and a method used to compute
+    distance between atoms.
+    """
     def dist_atoms(self, atom2):
+        """Compute distance between two atoms.
+
+        Parameters
+        ----------
+        self : Atom object
+            First atom involved to compute the distance.
+        atom2 : Atom object
+            Second atom involved to compute the distance.
+
+        This function is used to get the distance between
+        two atoms in 3D space.
+
+        Returns
+        -------
+        float
+            Distance, in angströms.
+        """
         return(math.sqrt((self.coord_x-atom2.coord_x)**2 +
                          (self.coord_y-atom2.coord_y)**2 +
                          (self.coord_z-atom2.coord_z)**2))
@@ -11,6 +33,25 @@ class Atom:
 
     def __init__(self, atom_id, atom_name, res_name, chain_id,
                  res_seq_nb, insertion_code, coordinates):
+        """Constructor of Atom class object.
+
+        Parameters
+        ----------
+        self : Atom object
+            The new atom being created.
+        atom_id : int
+            Sequential PDB atom number.
+        res_name : str
+            3-letter code of residue.
+        chain_id : str
+            Unique chain ID.
+        res_seq_nb : int
+            Actual amino acid sequence position.
+        insertion_code : str
+            Insertion letter when added in case of ambiguity.
+        coordinates : list[float, ...]
+            X,Y,Z coordinates of atom.
+        """
         self.atom_id = atom_id
         self.atom_name = atom_name
         self.res_name = res_name
@@ -53,6 +94,17 @@ class Helix(Structure):
         
 class Residue:
     def __init__(self, atoms_list, indice):
+        """Constructor of Residue class object.
+
+        Parameters
+        ----------
+        self : Residue object
+            The new residue being created.
+        atoms_list : list[Atom, ...]
+            List of Atom objects contained in this residue.
+        indice : int
+            Sequential Residue number.
+        """
         self.atoms = {}
         for atom in atoms_list:
             self.atoms[atom.atom_name] = atom
@@ -64,6 +116,24 @@ class Residue:
             self.indice = indice
             
     def get_amino_letter(self, res_name):
+        """Get 1-letter code of amino acid 3-letter code.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        res_name : str
+            3-letter code of the residue.
+
+        This function is used to convert the 3-letter code
+        format of the residue name to its 1-letter code
+        correspondance.
+
+        Returns
+        -------
+        str
+            3-letter amino acid code.
+        """
         code3 = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY',
                  'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER',
                  'THR', 'TRP', 'TYR', 'VAL']
@@ -75,6 +145,26 @@ class Residue:
         return code1[code3.index(res_name)]
 
     def h_bond(self, res2):
+        """Compute the H-Bond interaction between
+        two residues, in kcal/mole.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        res2 : Residue object
+            Second residue to evaluate the bond
+            energy.
+
+        This function computes the bond energy between 
+        two residues based on N,H,C and O atoms of each
+        residue.
+
+        Returns
+        -------
+        str
+            3-letter amino acid code.
+        """
         if("H" not in res2.atoms.keys()):
             return(False)
         # dimensionnal factor, in kcal/mole
@@ -97,6 +187,23 @@ class Residue:
     def get_turns(self, residues):
         """
         Get all the turns from a specific residue.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        residues : list[Residue, ...]
+            List of all the Residues in the protein.
+
+        This function determines if there is a n-turn
+        at a specific residue, and its type. 3,4,5-turn...
+        
+        Returns
+        -------
+        Turn object | Boolean 
+            If there is a n-turn, a Turn object is returned, 
+            if not, False is returned.
+        
         """
         turns = {}
         i = residues.index(self)
@@ -111,6 +218,25 @@ class Residue:
         return False
 
     def get_bends(self, residues):
+        """Get the bend from a specific residue.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        residues : list[Residue, ...]
+            List of all the Residues in the protein.
+
+        This function determines if there is a bend
+        at a specific residue. Used to determine Kappa
+        angle.
+        
+        Returns
+        -------
+        list[float, str] 
+            Kappa angle value, and S to signify bend presence,
+            if Kappa>70°.
+        """
         i = residues.index(self)
         if i >=2 and i <len(residues)-2:
             angle = math.degrees(vector_angles(
@@ -125,9 +251,24 @@ class Residue:
             return [360.0, '']
 
     def get_helix(self, residues):
-        """
+        """Return if there is an helix at a given residue,
-        Return if there is an helix at a given residue,
         as well as its type.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        residues : list[Residue, ...]
+            List of all the Residues in the protein.
+
+        This function determines if there is an helix
+        at a specific residue. 
+
+        Returns
+        -------
+        list[str, int]
+            The type of helix G(=3),H(=4),I(=5), and the sequential
+            position of residue where the helix starts.
         """
         i = residues.index(self)
         # if there are no turns or it is the first residue, skip
@@ -139,6 +280,21 @@ class Residue:
         return(False)
 
     def get_ladder(self, residues):
+        """Return if there is a ladder at a given residue.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        residues : list[Residue, ...]
+            List of all the Residues in the protein.
+
+        Returns
+        -------
+        dict | Boolean
+            If there is a ladder, return it as a dictionnary, else returns False.
+        
+        """
         i = residues.index(self)
         if i != 0:
             if self.get_bridges(residues):
@@ -153,6 +309,23 @@ class Residue:
         return False   
 
     def get_tco(self, residues):
+        """Cosine of angle between C=O of residue I and C=O of residue I-1. 
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        residues : list[Residue, ...]
+            List of all the Residues in the protein.
+
+        For α-helices, TCO is near +1, for β-sheets TCO is near -1. 
+        Not used for structure definition. 
+
+        Returns
+        -------
+        float
+            Cosine of angle.
+        """
         i = residues.index(self)
         if(i!=0):
             res2 = residues[i-1]
@@ -166,6 +339,24 @@ class Residue:
         return(math.cos(angle))
 
     def get_chirality(self, residues):
+        """Return the chirality at a given residue.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        residues : list[Residue, ...]
+            List of all the Residues in the protein.
+
+        Virtual torsion angle (dihedral angle) defined by the four 
+        Cα atoms of residues I-1,I,I+1,I+2.Used to define chirality 
+        (structure code '+' or '-'). 
+
+        Returns
+        -------
+        list[float, str]
+            Dihedral angle and chirality sign. 
+        """
         i = residues.index(self)
         if (i >=1 and i < len(residues)-2):
             chirality = {}
@@ -187,6 +378,22 @@ class Residue:
         return [angle, sign]
     
     def get_phi_psi(self, residues):
+        """Compute Phi and Psi angles of protein
+        backbone around the peptidic bound at each
+        residue.
+
+        Parameters
+        ----------
+        self : Residue object
+            Reference to Residue instance.
+        residues : list[Residue, ...]
+            List of all the Residues in the protein.
+
+        Returns
+        -------
+        tuple(float, float)
+        Phi and Psi angles.
+        """
         i = residues.index(self)  
         if(i==0):
             phi = 360.0
@@ -205,6 +412,23 @@ class Residue:
 
         return((phi, psi))
 def get_bridges(residues):
+    """Construct all bridges between residues in the 
+    protein.
+    
+    Parameters
+    ----------
+    residues : list[Residue, ...]
+        List of all the Residues in the protein.
+    
+    Constructs a list of Bridges objects, characterized by
+    their involved residues positions, sequential positions,
+    type (parallel or antiparallel).
+    
+    Returns
+    -------
+    list(Bridge)
+    Returns the list of bridges.
+    """
     bridges = {}
     bridge = {}
     strongest_bridge = {}
@@ -277,6 +501,23 @@ def get_bridges(residues):
         return(False)
 
 def get_ladders(bridges, residues):
+    """Construct all ladders in the protein.
+    
+    Parameters
+    ----------
+    residues : list[Residue, ...]
+        List of all the Residues in the protein.
+    
+    Constructs a dictionnary of ladders, characterized by
+    their involved residues positions, sequential positions,
+    involved bridges. Ladders are made of consecutive bridges
+    of the same type.
+    
+    Returns
+    -------
+    dict
+        Returns the dictionnary of ladders.
+    """
     ladders = {}
     i = 1
     while i < len(residues):
@@ -299,6 +540,9 @@ def get_ladders(bridges, residues):
     return ladders
 
 def connected_ladders(ladd_1, ladd_2):
+    """Check if two ladders, ladd_1 and ladd_2, 
+    are linked by one shared bridge.
+    """
     links = []
     for bridge in ladd_1['bridges']:
         if bridge.res_partner in res_list(ladd_2):
@@ -335,6 +579,9 @@ def get_sheets(ladders):
     return sheets
    
 def res_list(ladder):
+    """Iterate over ladders, from start to end,
+    to get the list of residues inside it.
+    """
     l=[]
     for i in range(ladder['i'], ladder['j']):
         l.append(i)
@@ -357,6 +604,9 @@ def build_turns_patterns(residues):
     return[turns_3, turns_4, turns_5]
 
 def build_helix_patterns(residues):
+    """Builds series of letters (G, H  or I) corresponding
+    to helix types (3, 4 or 5). Used for the final output.
+    """
     helix_3 = {}
     helix_4 = {}
     helix_5 = {}
@@ -375,6 +625,10 @@ def build_helix_patterns(residues):
     return[helix_3, helix_4, helix_5]
 
 def print_helix_pattern(residues, res, helix):
+    """Used to print the Helix type at a specific
+    residue position if there is one. Used for the final
+    output.
+    """
     i = residues.index(res)+1
     if i in helix.keys():
         return (helix[i])
@@ -382,6 +636,11 @@ def print_helix_pattern(residues, res, helix):
         return(' ')
 
 def print_turn_pattern(residues, res, turns):
+    """Builds series of symbols (>, <  or n) corresponding
+    to turn types (n=3, 4 or 5). Used for the final output.
+    '>' represent start of the turn series.
+    '<' represent end of the turn.
+    """
     i = residues.index(res)+1
     if i in turns.keys() and not i-1 in turns.keys():
         return(">")
@@ -392,7 +651,15 @@ def print_turn_pattern(residues, res, turns):
     else:
         return(' ')
     
-def raw_ladders_print(ladders, ):
+def raw_ladders_print(ladders):
+    """Function to print roughly the connected ladders by bridges in 
+    the final output.
+    This has been use to illustrate sheets principles.
+    Only shown if the -v or --verbose flag is used when running
+    dssp.py.
+
+    """
+    print("### CONNECTED LADDERS ###")
     for ladd1 in ladders:
         ladd_1 = ladders[ladd1]
         for ladd2 in ladders:

src/geom.py

@@ -1,5 +1,10 @@
 import math
 import numpy as np
+
+# A couple of obvious functions for vector manipulation.
+# Lack of time caused this not to be OOP-implemented, or furthermore
+# documented.
+
 def vector_from_pos(a, b):
     xAB = b[0]-a[0]
     yAB = b[1]-a[1]

src/pdb.py

@@ -112,11 +112,30 @@ class PDBFile:
             #Need to add hydrogens 
             self.add_hydrogens(filename)
 
-    def check_hydrogens(self, atoms):
-        print("ENTER CHECK HYDROGEN")
-        return True
-
     def add_hydrogens(self, filename, output_pdb=None):
+        """Add hydrogens to N atoms of the protein backbone.
+
+        Parameters
+        ----------
+        self : PDBFile
+            The instance of PDB File to get the header from.
+
+        filename : str
+            The name of the file to be exploited by PyMOL lib.
+        output_pdb : str
+            If defined, the name of the file to be saved py PyMOL,
+            with hydrogens added.
+
+        This function is used to add the missing H atoms on the
+        Nitrogens atoms of the protein backbone, when H are missing
+        from PDB file. H are needed to perform all the bond interactions-
+        based calculations.
+
+        Returns
+        -------
+        dict
+            A dictionnary of all Hydrogens added to the backbone.
+        """
         pymol.finish_launching(['pymol', '-qc'])
         pymol.cmd.load(filename)
         pymol.cmd.select("nitrogens",'name n')  
@@ -129,7 +148,7 @@ class PDBFile:
                                 'else ' \
                                 'stored.pos[resi].append([name,resi,x,y,z])')
         if(output_pdb!=None):
-            pymol.cmd.save(output_file)
+            pymol.cmd.save(output_pdb)
         for resi in self.residues:
             if(str(resi.resid) in pymol.stored.pos.keys()):
                 # iterate over residue H atoms