diff --git a/README.md b/README.md index 7814c0e..854386c 100644 --- a/README.md +++ b/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 diff --git a/doc/Notice.pdf b/doc/Notice.pdf deleted file mode 100644 index d8c54fa..0000000 Binary files a/doc/Notice.pdf and /dev/null differ diff --git a/doc/Sujets_Projets_Courts_2019_2020.pdf b/doc/Sujets_Projets_Courts_2019_2020.pdf deleted file mode 100644 index 74bdb31..0000000 Binary files a/doc/Sujets_Projets_Courts_2019_2020.pdf and /dev/null differ diff --git a/doc/notes.org b/doc/notes.org deleted file mode 100644 index 9049b87..0000000 --- a/doc/notes.org +++ /dev/null @@ -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 diff --git a/doc/pdb_format_33_rot90.png b/doc/pdb_format_33_rot90.png deleted file mode 100644 index 1908be7..0000000 Binary files a/doc/pdb_format_33_rot90.png and /dev/null differ diff --git a/src/chem.py b/src/chem.py index 79797ae..d6ca229 100644 --- a/src/chem.py +++ b/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 =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: diff --git a/src/geom.py b/src/geom.py index e4e5b9c..28e27a1 100644 --- a/src/geom.py +++ b/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] diff --git a/src/pdb.py b/src/pdb.py index f655648..6f3674f 100644 --- a/src/pdb.py +++ b/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