cleaner output

src/chem.py

@@ -0,0 +1,507 @@
+import math
+#from structure import *
+from geom import *
+
+class Atom:
+
+    def dist_atoms(self, atom2):
+        return(math.sqrt((self.coord_x-atom2.coord_x)**2 +
+                         (self.coord_y-atom2.coord_y)**2 +
+                         (self.coord_z-atom2.coord_z)**2))
+
+
+    def __init__(self, atom_id, atom_name, res_name, chain_id,
+                 res_seq_nb, insertion_code, coordinates):
+        self.atom_id = atom_id
+        self.atom_name = atom_name
+        self.res_name = res_name
+        self.chain_id = chain_id
+        self.res_seq_nb = res_seq_nb
+        self.insertion_code = insertion_code
+        self.coord_x = coordinates[0]
+        self.coord_y = coordinates[1]
+        self.coord_z = coordinates[2]
+        self.coords = coordinates
+
+class Structure:
+    def __init__(self, res):
+        self.res = res
+
+class Turn(Structure):
+
+    def __init__(self, turn_type, res):
+        self.turn_type = turn_type
+        Structure.res = res
+                        
+class Bridge(Structure):
+
+    def __init__(self, bridge_type, res_num, res_partner, indices):
+        self.bridge_type = bridge_type
+        self.res_num = res_num
+        self.res_partner = res_partner
+        Structure.res = res_num
+        self.i = indices[0]
+        self.j = indices[1]
+
+class Helix(Structure):
+
+    def __init__(self, residues, res_num, helix_type):
+        self.residues = residues
+        self.res_num = res_num
+        Structure.res = res_num
+        self.helix_type = helix_type
+
+        
+class Residue:
+    def __init__(self, atoms_list, indice):
+        self.atoms = {}
+        for atom in atoms_list:
+            self.atoms[atom.atom_name] = atom
+            self.resid = atom.res_seq_nb
+            self.res_name = atom.res_name
+            self.res_letter = self.get_amino_letter(atom.res_name)
+            self.chain_id = atom.chain_id
+            self.insertion_code = atom.insertion_code
+            self.indice = indice
+            
+    def get_amino_letter(self, res_name):
+        code3 = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY',
+                 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER',
+                 'THR', 'TRP', 'TYR', 'VAL']
+
+
+        code1 = ['A','R','N','D','C','E','Q','G','H','I','L','K','M','F','P',
+                 'S','T','W','Y','V']
+
+        return code1[code3.index(res_name)]
+
+    def h_bond(self, res2):
+        if("H" not in res2.atoms.keys()):
+            return(False)
+        # dimensionnal factor, in kcal/mole
+        f = 332
+        # partial charges
+        q1 = 0.42
+        q2 = 0.20
+        # distance between O-N atoms, in angströms
+        r_ON = self.atoms["O"].dist_atoms(res2.atoms["N"])
+        # distance between C-H atoms, in angströms
+        r_CH = self.atoms["C"].dist_atoms(res2.atoms["H"])
+        # distance between O-H atoms, in angströms
+        r_OH = self.atoms["O"].dist_atoms(res2.atoms["H"])
+        # distance between C-N atoms, in angströms
+        r_CN = self.atoms["C"].dist_atoms(res2.atoms["N"])
+        # electrostatic interaction energy, in kcal/mole
+        E = q1*q2*(1/r_ON + 1/r_CH - 1/r_OH - 1/r_CN)*f
+        return(E)
+
+    def get_turns(self, residues):
+        """
+        Get all the turns from a specific residue.
+        """
+        turns = {}
+        i = residues.index(self)
+        k = 0
+        for j in range(3,6):
+                if(i+j<len(residues)):
+                    if(self.h_bond(residues[i+j])<-0.5):
+                        k = j
+        if k != 0:    
+            #print(k,"TURN", residues[i].resid, residues[i+k].resid)
+            return Turn(k,residues[i].resid)
+
+        return False
+
+    def get_bends(self, residues):
+        i = residues.index(self)
+        if i >=2 and i <len(residues)-2:
+            angle = math.degrees(vector_angles(vectors_substr(position_vector(residues[i].atoms["CA"].coords),
+                                                              position_vector(residues[i-2].atoms["CA"].coords)),
+                                               vectors_substr(position_vector(residues[i+2].atoms["CA"].coords),
+                                                              position_vector(residues[i].atoms["CA"].coords))))
+            if(angle>70):
+                return [angle, 'S']
+            return [angle, '']
+        else:
+            return [360.0, '']
+
+    def get_bridges(self, residues):
+        bridges = {}
+        bridge = {}
+        strongest_bridge = {}
+        i = residues.index(self)
+        if(i >= 1 and i < len(residues)-4):
+            E_min = 0
+            for j in range(i+2,len(residues)-1):
+                # select triplet with the minimal energy 
+
+                if(residues[i-1].h_bond(residues[j])<-0.5
+                   and residues[j].h_bond(residues[i+1])<-0.5):
+                    bridge = {'res1':residues[i-1].h_bond(residues[j]),
+                              'res2':residues[j].h_bond(residues[i+1]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"para"}
+
+                if(residues[j-1].h_bond(residues[i])<-0.5
+                   and residues[i].h_bond(residues[j+1])<-0.5):
+                    bridge = {'res1':residues[j-1].h_bond(residues[i]),
+                              'res2':residues[i].h_bond(residues[j+1]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"para"}
+
+                if(residues[i].h_bond(residues[j])<-0.5
+                   and residues[j].h_bond(residues[i])<-0.5):
+                    bridge = {'res1':residues[i].h_bond(residues[j]),
+                              'res2':residues[j].h_bond(residues[i]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"anti"}
+
+                if(residues[i-1].h_bond(residues[j+1])<-0.5
+                   and residues[j-1].h_bond(residues[i+1])<-0.5):
+                    bridge = {'res1':residues[i-1].h_bond(residues[j+1]),
+                              'res2':residues[j-1].h_bond(residues[i+1]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"anti"}
+
+                if(bridge):
+                    if(bridge['res1']+bridge['res2']<E_min):
+                        E_min = bridge['res1']+bridge['res2']
+                        strongest_bridge = bridge
+                        bridge = {}
+                        coord_bridge = [i,j]
+            # finally add the strongest bridge at i and j pos
+            if(strongest_bridge):
+                bridges[coord_bridge[0]] = (Bridge(strongest_bridge['btype'],
+                                                   strongest_bridge['ipos'],
+                                                   strongest_bridge['jpos']))
+                
+                bridges[coord_bridge[1]] = (Bridge(strongest_bridge['btype'],
+                                                   strongest_bridge['jpos'],
+                                                   strongest_bridge['ipos']))
+        if(len(bridges)>0):
+            return(bridges[coord_bridge[0]])
+        else:
+            return(False)
+
+    def get_helix(self, residues):
+        """
+        Return if there is an helix at a given residue,
+        as well as its type.
+        """
+        i = residues.index(self)
+        # if there are no turns or it is the first residue, skip
+        if i == 0:
+            return False
+        
+        if(self.get_turns(residues) and residues[i-1].get_turns(residues)):
+            #print(self.get_turns(residues).turn_type,"- HELIX at", residues[i].indice)
+            return(self.get_turns(residues).turn_type, residues[i].indice)
+        return(False)
+
+    def get_ladder(self, residues):
+        #ladders = {}
+        i = residues.index(self)
+        if i != 0:
+            if self.get_bridges(residues):
+                if (residues[i-1].get_bridges(residues)):
+                    local_bridge = self.get_bridges(residues)
+                    consec_bridge = residues[i-1].get_bridges(residues)
+                    if local_bridge.bridge_type == consec_bridge.bridge_type:
+                        #print("ladder", consec_bridge.res_num, local_bridge.res_num)
+                        ladder = {'start':consec_bridge.res_num,
+                                  'end':local_bridge.res_num,
+                                  'bridges':[consec_bridge, local_bridge]}
+                        return ladder
+        return False   
+
+    def get_tco(self, residues):
+        i = residues.index(self)
+        if(i!=0):
+            res2 = residues[i-1]
+            CO_res1 = vector_from_pos(self.atoms["C"].coords,
+                                      self.atoms["O"].coords)
+            CO_res2 = vector_from_pos(res2.atoms["C"].coords,
+                                      res2.atoms["O"].coords)
+            angle = vector_angles(CO_res1, CO_res2)
+        else:
+            angle = math.pi/2
+        return(math.cos(angle))
+
+    def get_chirality(self, residues):
+        i = residues.index(self)
+        if (i >=1 and i < len(residues)-2):
+            chirality = {}
+            angle = calc_dihedral(residues[i-1].atoms["CA"].coords,
+                                  residues[i].atoms["CA"].coords,
+                                  residues[i+1].atoms["CA"].coords,
+                                  residues[i+2].atoms["CA"].coords)
+
+            if(angle>0 and angle<=180):
+                sign="+"
+
+            if(angle<=0 and angle > -180):
+                sign="-"
+
+        else:
+            angle = 360.0
+            sign = ''
+
+        return [angle, sign]
+    
+    def get_phi_psi(self, residues):
+        i = residues.index(self)  
+        if(i==0):
+            phi = 360.0
+        else:
+            phi = calc_dihedral(residues[i-1].atoms["C"].coords,
+                                residues[i].atoms["N"].coords,
+                                residues[i].atoms["CA"].coords,
+                                residues[i].atoms["C"].coords)
+        if(i==len(residues)-1):
+            psi = 360.0
+        else:
+            psi = calc_dihedral(residues[i].atoms["N"].coords,
+                                residues[i].atoms["CA"].coords,
+                                residues[i].atoms["C"].coords,
+                                residues[i+1].atoms["N"].coords)
+
+        return((phi, psi))
+def get_bridges(residues):
+    bridges = {}
+    bridge = {}
+    strongest_bridge = {}
+    for i in range(1,len(residues)-4):
+        E_min = 0
+        for j in range(i+2,len(residues)-1):
+            # select triplet with the minimal energy 
+            
+            if(residues[i-1].h_bond(residues[j])<-0.5
+               and residues[j].h_bond(residues[i+1])<-0.5):
+                bridge = {'res1':residues[i-1].h_bond(residues[j]),
+                          'res2':residues[j].h_bond(residues[i+1]),
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"para"}
+                   
+            if(residues[j-1].h_bond(residues[i])<-0.5
+               and residues[i].h_bond(residues[j+1])<-0.5):
+                bridge = {'res1':residues[j-1].h_bond(residues[i]),
+                          'res2':residues[i].h_bond(residues[j+1]),
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"para"}
+
+            if(residues[i].h_bond(residues[j])<-0.5
+               and residues[j].h_bond(residues[i])<-0.5):
+                bridge = {'res1':residues[i].h_bond(residues[j]),
+                          'res2':residues[j].h_bond(residues[i]),
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"anti"}
+                    
+            if(residues[i-1].h_bond(residues[j+1])<-0.5
+               and residues[j-1].h_bond(residues[i+1])<-0.5):
+                bridge = {'res1':residues[i-1].h_bond(residues[j+1]),
+                          'res2':residues[j-1].h_bond(residues[i+1]),
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"anti"}
+
+            if(bridge):
+                if(bridge['res1']+bridge['res2']<E_min):
+                    E_min = bridge['res1']+bridge['res2']
+                    strongest_bridge = bridge
+                    coord_bridge = [i,j]
+                    bridge = {}
+
+        # finally add the strongest bridge at i and j pos
+        if(strongest_bridge):
+
+            bridges[strongest_bridge['i']] = (Bridge(strongest_bridge['btype'],
+                                                     strongest_bridge['ipos'],
+                                                     strongest_bridge['jpos'],
+                                                     [strongest_bridge['i'],
+                                                      strongest_bridge['j']]))
+                                              
+            bridges[strongest_bridge['j']] = (Bridge(strongest_bridge['btype'],
+                                                     strongest_bridge['ipos'],
+                                                     strongest_bridge['jpos'],
+                                                     [strongest_bridge['i'],
+                                                      strongest_bridge['j']]))
+    if(len(bridges)>0):
+        return(bridges)
+    else:
+        return(False)
+
+def get_ladders(bridges, residues):
+    ladders = {}
+    i = 1
+    while i < len(residues):
+        k = 1
+        if i in bridges.keys():
+            temp_bridges = [bridges[i]]
+            while ((i+k in bridges.keys()) and
+                   (bridges[i].bridge_type == bridges[i+k].bridge_type)):
+                temp_bridges.append(bridges[i+k])
+                k+=1
+        if k>1:
+            #print("ladder", bridges[i].res_num, bridges[i+k-1].res_num)
+            ladders[i] = {'start':bridges[i].res_num,
+                          'end':bridges[i+k-1].res_num,
+                          'bridges':temp_bridges,
+                          'i':i,
+                          'j':i+k-1}
+            i+=k-1
+        else:
+            i+=1
+    return ladders
+
+def connected_ladders(ladd_1, ladd_2):
+    links = []
+    for bridge in ladd_1['bridges']:
+        if bridge.res_partner in res_list(ladd_2):
+            return ladd_2
+
+    return False
+
+def connected_ladders2(ladd_1, ladd_2):
+    links = []
+    for bridge in ladd_1['bridges']:
+        if bridge.res_partner in res_list(ladd_2):
+            return([ladd_1['i'], ladd_1['j'], bridge.i, bridge.j,
+                    ladd_2['i'], ladd_2['j']])
+            #return ladd_2
+
+    return False
+
+
+def get_sheets(ladders):
+    """
+    Bridges between ladders.
+    Check if 1 bridge between one ladder and one or more other ladders.
+    Iterate over all residues of one ladder and check if bridge with other residues
+    of the other ladders.
+    """
+    ladds = [ elem for elem in ladders.values() ]
+    sheets = {}
+
+    corresp = {}
+    for ladd1 in ladds:
+        for ladd2 in ladds:
+            if connected_ladders(ladd1, ladd2)!=False:
+                corresp_list = [ elem for elem in corresp.keys() ]
+                if ladd1['i'] not in corresp_list and ladd2['i'] not in corresp_list:
+                    ind = len(sheets.keys())
+                    sheets[ind] = []
+                    sheets[ind].append(ladd1)
+                    sheets[ind].append(ladd2)
+                    corresp[ladd1['i']] = ind
+                    corresp[ladd2['i']] = ind
+                elif ladd2 not in corresp_list and ladd1 in corresp_list:
+                    sheets[corresp[ladd1['i']]].append(ladd2)
+                    corresp[ladd2['i']] = corresp[ladd1['i']]
+                elif ladd1 not in corresp_list and ladd2 in corresp_list:
+                    sheets[corresp[ladd2['i']]].append(ladd1)
+                    corresp[ladd1['i']] = corresp[ladd2['i']]
+                    
+    return sheets
+
+def get_sheets2(ladders):
+    """
+    Bridges between ladders.
+    Check if 1 bridge between one ladder and one or more other ladders.
+    Iterate over all residues of one ladder and check if bridge with other residues
+    of the other ladders.
+    """
+    sheets = {}    
+    for ladder in ladders:
+        for ladd2 in ladders:
+            if connected_ladders(ladders[ladder], ladders[ladd2]):
+                bridge_i = connected_ladders2(ladders[ladder], ladders[ladd2])[2]
+                bridge_j = connected_ladders2(ladders[ladder], ladders[ladd2])[3]
+                print("ladder",ladders[ladder]['i'], ladders[ladder]['j'],"bridge",bridge_i, bridge_j,
+                      "ladder 2",ladders[ladd2]['i'], ladders[ladd2]['j'])
+
+               
+def res_list(ladder):
+    # TODO : method in ladder class
+    l=[]
+    for i in range(ladder['i'], ladder['j']):
+        l.append(i)
+    return(l)
+
+def build_turns_patterns(residues):
+    turns_3 = {}
+    turns_4 = {}
+    turns_5 = {}
+    for i,res in enumerate(residues):
+        turn = residues[i].get_turns(residues)
+        if(turn):
+            for k in range(turn.turn_type):
+                if turn.turn_type == 3:
+                    turns_3[i+1+k] = turn.turn_type
+                if turn.turn_type == 4:
+                    turns_4[i+1+k] = turn.turn_type
+                if turn.turn_type == 5:
+                    turns_5[i+1+k] = turn.turn_type
+    return[turns_3, turns_4, turns_5]
+
+def build_helix_patterns(residues):
+    helix_3 = {}
+    helix_4 = {}
+    helix_5 = {}
+    for i,res in enumerate(residues):
+        helix = residues[i].get_helix(residues)
+        if(helix):
+            helix_type = residues[i].get_helix(residues)[0]
+            helix_pos = residues[i].get_helix(residues)[1]
+            #print("TYPE", helix_type)
+            for k in range(helix_type):
+                if helix_type == 3:
+                    helix_3[i+1+k] = "G"
+                if helix_type == 4:
+                    helix_4[i+1+k] = "H"
+                if helix_type == 5:
+                    helix_5[i+1+k] = "I"
+        #print(helix_3)
+    return[helix_3, helix_4, helix_5]
+
+def print_helix_pattern(residues, res, helix):
+    i = residues.index(res)+1
+    if i in helix.keys():
+        return (helix[i])
+    else:
+        return(' ')
+
+def print_turn_pattern(residues, res, turns):
+    i = residues.index(res)+1
+    if i in turns.keys() and not i-1 in turns.keys():
+        return(">")
+    elif i in turns.keys() and i-1 in turns.keys():
+        return(turns[i])
+    elif i not in turns.keys() and i-1 in turns.keys():
+        return("<")
+    else:
+        return(' ')
+    
+def raw_ladders_print(ladders, ):
+    for ladd1 in ladders:
+        ladd_1 = ladders[ladd1]
+        for ladd2 in ladders:
+            ladd_2 = ladders[ladd2]
+            for bridge in ladd_1['bridges']:
+                if bridge.j in res_list(ladd_2):  
+                    print("ladder", ladd_1['i'],"-",ladd_1['j'], "|", bridge.i, "...", bridge.j, "| ladder", ladd_2['i'], ladd_2['j'])
+

src/chem.py~

@@ -1,5 +1,7 @@
 import math
-from structure import *
+#from structure import *
+from geom import *
+
 class Atom:
 
     def dist_atoms(self, atom2):
@@ -21,6 +23,35 @@ class Atom:
         self.coord_z = coordinates[2]
         self.coords = coordinates
 
+class Structure:
+    def __init__(self, res):
+        self.res = res
+
+class Turn(Structure):
+
+    def __init__(self, turn_type, res):
+        self.turn_type = turn_type
+        Structure.res = res
+                        
+class Bridge(Structure):
+
+    def __init__(self, bridge_type, res_num, res_partner, indices):
+        self.bridge_type = bridge_type
+        self.res_num = res_num
+        self.res_partner = res_partner
+        Structure.res = res_num
+        self.i = indices[0]
+        self.j = indices[1]
+
+class Helix(Structure):
+
+    def __init__(self, residues, res_num, helix_type):
+        self.residues = residues
+        self.res_num = res_num
+        Structure.res = res_num
+        self.helix_type = helix_type
+
+        
 class Residue:
     def __init__(self, atoms_list, indice):
         self.atoms = {}

src/dssp.py

@@ -1,104 +1,86 @@
 import pdb
-import atom
 import sys
 from structure import *
-from atom import *
+from chem import *
+import argparse
 
-def print_dssp():
+from datetime import date
 
-    empt_line=" " * 145
+def print_dssp(pdb_file, verbose=False, raw_ladd=False):
+    # gets export date
+    export_date = date.today().strftime("%d/%m/%Y")
       
+    ################### OUTPUT ####################
+    # print DSSP-style formatted header of PDB file
+    print("################### DSSP OUTPUT ####################\n####",
+          "Exported on",export_date,
+          "\n####################################################")
+    for elem in pdb_file.get_header() : 
+        print(pdb_file.get_header()[elem], end="")
+
+    # Get preliminary data for print loop
+    residues = pdb_file.residues
+    bridges = get_bridges(residues)
+    ladders = get_ladders(bridges, residues)
+    sheets = get_sheets(ladders)  
+    helix = build_helix_patterns(residues)
+    # iterating over residues
+    # print header of table
+    print('   # RESIDUE AA  STRUCTURE     BP1   TCO   KAPPA  ALPHA   ' \
+    'PHI    PSI    X-CA   Y-CA   Z-CA')
+    for i,res in enumerate(residues):
+        #res.get_turns(residues, turns)
+        kappa = res.get_bends(residues)[0]
+        bend_symbol = res.get_bends(residues)[1]
+        t_co = res.get_tco(residues)
+        alpha = res.get_chirality(residues)[0]
+        chirality = res.get_chirality(residues)[1]
+        phi = res.get_phi_psi(residues)[0]
+        psi = res.get_phi_psi(residues)[1]
+        x_ca = res.atoms["CA"].coord_x
+        y_ca = res.atoms["CA"].coord_y
+        z_ca = res.atoms["CA"].coord_z
+        turns = build_turns_patterns(residues)
+        turn_3 = print_turn_pattern(residues, res, turns[0])
+        turn_4 = print_turn_pattern(residues, res, turns[1])
+        turn_5 = print_turn_pattern(residues, res, turns[2])
+        helix_3 = print_helix_pattern(residues, res, helix[0])
+        helix_4 = print_helix_pattern(residues, res, helix[1])
+        helix_5 = print_helix_pattern(residues, res, helix[2])
+        if i in bridges.keys():
+            bp1 = bridges[i].j
+        else:
+            bp1 = 0
+       
+        print('{0: >5}{1: >5}{2: >2}{3: >2}{4: >2}{5: >2}{6: >2}{7: >2}'\
+              '{8: >2}{9: >2}{10: >2}{11: >1}{12: >5}{13: >7}{14: >7}' \
+              '{15: >7}{16: >7}{17: >7}{18: >7}{19: >7}'\
+              '{20: >7}'.format(i+1, res.resid,
+                                res.chain_id,
+                                res.res_letter,
+                                helix_3, helix_4,
+                                helix_5, turn_3,
+                                turn_4, turn_5, bend_symbol,
+                                chirality, bp1,round(t_co, 3),
+                                round(kappa,1), round(alpha, 1),
+                                round(phi, 1), round(psi,1),
+                                round(x_ca,1), round(y_ca,1),
+                                round(z_ca,1)))
+    if(raw_ladd):
+        raw_ladders_print(ladders)
+
+parser = argparse.ArgumentParser()
 if(len(sys.argv)<2):
     print("Not enough arguments! Run with --help to learn more about proper"
           "call structure and parameters.")
 else:
+    parser.add_argument("file", type=str, help="input file in PDB format")
+    parser.add_argument("-v","--verbose", help="add ladders verbose output",
+                        action="store_true")
+    args = parser.parse_args()
     pdb_file = pdb.PDBFile(sys.argv[1])
-  
+    if(args.verbose):
-    #turns = get_turns(pdb_file.residues)
+        print_dssp(pdb_file, raw_ladd=True)
-    #get_helix(pdb_file.residues, turns)
-    #get_bends(pdb_file.residues)
-    # get_chirality(pdb_file.residues)
-
-
-    # print("NBRIDGES",len(bridges))
-
-    # bridges = get_bridges(pdb_file.residues)
-
-    # get_bonds(pdb_file.residues)
-
-
-    # e_min = 0
-    # for res in residues:
-    #     ene = residues[1].h_bond(res)
-    #     if ene <= e_min:
-    #         e_min = ene
-    #         best_res = res
-    #print(residues[1].resid, best_res.resid,e_min)
-
-    #get_phi_psi(residues)
-    #print(residues[2].atoms, residues[0].resid, residues[1].resid)
-    #print(get_TCO(residues[2],residues[3]))
-
-
-# for i,ladder in enumerate(ladders.values()):
-#     print(chr(65+i))
-    
-# for ind, sheet in sheets.items():
-#     print(ind, sheet)
-#     for ladder in sheet:
-#         print(ladder['start'], ladder['end'])
-
-
-################### OUTPUT ####################
-# print DSSP-style formatted header of PDB file
-for elem in pdb_file.get_header() : 
-    print(pdb_file.get_header()[elem], end="")
-
-# Get preliminary data for print loop
-residues = pdb_file.residues
-bridges = get_bridges(residues)
-ladders = get_ladders(bridges, residues)
-sheets = get_sheets(ladders)  
-helix = build_helix_patterns(residues)
-# iterating over residues
-
-
-for i,res in enumerate(residues):
-    #res.get_turns(residues, turns)
-    kappa = res.get_bends(residues)[0]
-    bend_symbol = res.get_bends(residues)[1]
-    t_co = res.get_tco(residues)
-    alpha = res.get_chirality(residues)[0]
-    chirality = res.get_chirality(residues)[1]
-    phi = res.get_phi_psi(residues)[0]
-    psi = res.get_phi_psi(residues)[1]
-    x_ca = res.atoms["CA"].coord_x
-    y_ca = res.atoms["CA"].coord_y
-    z_ca = res.atoms["CA"].coord_z
-    turns = build_turns_patterns(residues)
-    turn_3 = print_turn_pattern(residues, res, turns[0])
-    turn_4 = print_turn_pattern(residues, res, turns[1])
-    turn_5 = print_turn_pattern(residues, res, turns[2])
-    helix_3 = print_helix_pattern(residues, res, helix[0])
-    helix_4 = print_helix_pattern(residues, res, helix[1])
-    helix_5 = print_helix_pattern(residues, res, helix[2])
-    if i in bridges.keys():
-        bp1 = bridges[i].j
     else:
-        bp1 = 0
+        print_dssp(pdb_file)
-    print(i+1, res.resid, res.chain_id, res.res_letter, helix_3, helix_4, helix_5, turn_3, turn_4, turn_5, bend_symbol, chirality, bp1,round(t_co, 3), round(kappa,1),
-           round(alpha, 1), round(phi, 1), round(psi,1), round(x_ca,1), round(y_ca,1), round(z_ca,1))
     
-
-for ladd1 in ladders:
-    ladd_1 = ladders[ladd1]
-    for ladd2 in ladders:
-        ladd_2 = ladders[ladd2]
-        for bridge in ladd_1['bridges']:
-            if bridge.j in res_list(ladd_2):  
-                print("ladder", ladd_1['i'],"-",ladd_1['j'], "|", bridge.i, "...", bridge.j, "| ladder", ladd_2['i'], ladd_2['j'])
-
-# for ind, sheet in sheets.items():
-#     print(ind, sheet)
-#     for ladder in sheet:
-#         print(ladder['start'], ladder['end'])

src/geom.py~

@@ -1,66 +0,0 @@
-
-def vector_from_pos(a, b):
-    xAB = b[0]-a[0]
-    yAB = b[1]-a[1]
-    zAB = b[2]-a[2]
-    coord_AB = [xAB,yAB,zAB]
-    return coord_AB
-
-def vector_norm(v):
-    norm = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
-    return norm
-
-def dot_product(v1,v2):
-    dot_product = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
-    return(dot_product)
-
-def position_vector(c):
-    vector = vector_from_pos([0,0,0],c)
-    return vector
-
-def vectors_substr(v1, v2):
-    return ([v1[0]-v2[0],
-             v1[1]-v2[1],
-             v1[2]-v2[2]])
-
-def vector_angles(v1,v2):
-    dot_prod = dot_product(v1,v2)
-    norm_v1 = vector_norm(v1)
-    norm_v2 = vector_norm(v2)
-    term = dot_prod/(abs(norm_v1)*abs(norm_v2))
-    rad_angle = math.acos(term)
-    return rad_angle
-
-def calc_dihedral(u1, u2, u3, u4):
-    """ Calculate dihedral angle method. From bioPython.PDB
-    (adapted to np.array)
-    Calculate the dihedral angle between 4 vectors
-    representing 4 connected points. The angle is in
-    [-pi, pi].
-
-    Adapted function of dihedral_angle_from_coordinates.py 
-    by Eric Alcaide.
-    Source : https://gist.github.com/EricAlcaide
-    URL : https://gist.github.com/EricAlcaide/30d6bfdb8358d3a57d010c9a501fda56
-    """
-    #convert coords to numpy arrays
-    u1 = np.array(u1)
-    u2 = np.array(u2)
-    u3 = np.array(u3)
-    u4 = np.array(u4)
-    
-    a1 = u2 - u1
-    a2 = u3 - u2
-    a3 = u4 - u3
-
-    v1 = np.cross(a1, a2)
-    v1 = v1 / (v1 * v1).sum(-1)**0.5
-    v2 = np.cross(a2, a3)
-    v2 = v2 / (v2 * v2).sum(-1)**0.5
-    porm = np.sign((v1 * a3).sum(-1))
-    rad = np.arccos((v1*v2).sum(-1) / ((v1**2).sum(-1) * (v2**2).sum(-1))**0.5)
-    if not porm == 0:
-        rad = rad * porm
-
-    deg_angle = rad*(180/math.pi)
-    return deg_angle

src/pdb.py

@@ -1,5 +1,5 @@
 # custom imports
-from atom import *
+from chem import *
 # import libs
 import sys
 import pymol
@@ -67,9 +67,7 @@ class PDBFile:
         # hydrogens should represent in average 50% of total atoms...
         # We use 30% threshold...
         if(count_h/len(temp_atoms)<0.30):
-            #if(output_pdb==None):
+            #Need to add hydrogens 
-            print("Need to add hydrogens ! If you want the modified PDB file, " 
-                  "please use the -o output.pdb argument")
             self.add_hydrogens(filename)
 
     def check_hydrogens(self, atoms):