hdiff output

r29514/chirality.py 2017-01-21 10:36:52.274250000 +0000 r29513/chirality.py 2017-01-21 10:36:52.502250000 +0000
  1: #!/usr/bin/python  1: #!/usr/bin/python
  2:   2: 
  3: import read_amber as ra  3: import read_amber as ra
  4: import networkx as nx  4: import networkx as nx
  5: import sys  5: import sys
  6: from collections import defaultdict 
  7:   6: 
  8: class GhostAtom(ra.Atom):  7: class GhostAtom(ra.Atom):
  9:   pass  8:   pass
 10:   9: 
 11: # Set the default valence to 0, if the element isn't present in the defined_valences dictionary. 10: valences = { "H": 1,
 12: valences = defaultdict(int) 11:          "C": 4,
 13: defined_valences = { "H": 1, 12:          "N": 3,
 14:                      "C": 4, 13:          "O": 2,
 15:                      "N": 3, 14:          "S": 2,
 16:                      "O": 2, 15:          "Se": 2,
 17:                      "S": 2, 16:          "P": 5 }
 18:                      "Se": 2, 
 19:                      "P": 5 } 
 20:  
 21: for k, v in defined_valences.items(): 
 22:     valences[k] = v 
 23:  17: 
 24: def chiral_candidates(atoms): 18: def chiral_candidates(atoms):
 25:   """ 
 26:   Returns those atoms which have 4 neighbours in the atom graph (and thus might be chiral). 
 27:  
 28:   :param atoms: atom graph 
 29:   :return: list of read_amber.Atom objects which have 4 neighbours in the atoms graph 
 30:   """ 
 31:   candidates = [atom for atom in atoms.nodes() if len(nx.edges(atoms, atom)) == 4] 19:   candidates = [atom for atom in atoms.nodes() if len(nx.edges(atoms, atom)) == 4]
 32:   return candidates 20:   return candidates
 33:  21: 
 34: def not_in_dict(atoms): 
 35:   """ 
 36:   Writes an error message to `chirality_error` and stderr, if an atom in the topology file is not present in the 
 37:   `defined_valences` dictionary. 
 38:  
 39:   :param atoms: atom graph 
 40:   """ 
 41:   missing = [] 
 42:   for atom in atoms.nodes(): 
 43:     try: 
 44:       _ = defined_valences[atom.element] 
 45:     except KeyError: 
 46:       missing.append(str(atom)) 
 47:   if missing: 
 48:     with open("chirality_warning", "w") as error_file: 
 49:       error_file.write("The following atoms were not present in the defined_valences dictionary\n") 
 50:       error_file.write("of the chirality script, and given a valence of 0 (not bonded):\n\n") 
 51:       sys.stderr.write("The following atoms were not present in the defined_valences dictionary\n") 
 52:       sys.stderr.write("of the chirality script, and given a valence of 0 (not bonded):\n\n") 
 53:       for atom_str in missing: 
 54:         error_file.write(atom_str + "\n") 
 55:       for atom_str in missing: 
 56:         sys.stderr.write(atom_str + "\n") 
 57:  
 58: def multi_bonds(atoms): 22: def multi_bonds(atoms):
 59:   not_in_dict(atoms) 
 60:   multibonded = [atom for atom in atoms.nodes()  23:   multibonded = [atom for atom in atoms.nodes() 
 61:                  if len(nx.edges(atoms, atom)) < valences[atom.element]] 24:                  if len(nx.edges(atoms, atom)) < valences[atom.element]]
 62:   for i, atom in enumerate(multibonded): 25:   for i, atom in enumerate(multibonded):
 63:     paired = False 26:     paired = False
 64:     for other in atoms.neighbors(atom): 27:     for other in atoms.neighbors(atom):
 65:       if isinstance(other, GhostAtom): 28:       if isinstance(other, GhostAtom):
 66:         paired = True 29:         paired = True
 67:         continue 30:         continue
 68:       if len(nx.edges(atoms, other)) < valences[other.element]: 31:       if len(nx.edges(atoms, other)) < valences[other.element]:
 69:         ghost_atom = GhostAtom(**(atom.__dict__)) 32:         ghost_atom = GhostAtom(**(atom.__dict__))
 70:         ghost_atom.name = atom.name + "*" 33:         ghost_atom.name = atom.name + "*"
 71:         ghost_other = GhostAtom(**(other.__dict__)) 34:         ghost_other = GhostAtom(**(other.__dict__))
 72:         ghost_other.name = other.name + "*" 35:         ghost_other.name = other.name + "*"
 73:         atoms.add_edge(other, ghost_atom) 36:         atoms.add_edge(other, ghost_atom)
 74:         atoms.add_edge(atom, ghost_other) 37:         atoms.add_edge(atom, ghost_other)
 75:         paired = True 38:         paired = True
 76:  39: #    if not paired:
  40: #      print "Could not find pair for atom", atom, atom.residue, atom.charge, nx.edges(atoms, atom)
  41:   
 77: def chiral_order(atoms, chiral_atom, depth=6): 42: def chiral_order(atoms, chiral_atom, depth=6):
 78:   # Create a list of ordered atoms to be passed back 43: # Create a list of ordered atoms to be passed back
 79:   ordered = [] 44:   ordered = []
 80:   # Do a quick check whether there are multiple hydrogens 45: # Do a quick check whether there are multiple hydrogens
 81:   neighbors = atoms.neighbors(chiral_atom) 46:   neighbors = atoms.neighbors(chiral_atom)
 82:   hydrogens = [atom for atom in neighbors if atom.element == "H"] 47:   hydrogens = [atom for atom in neighbors if atom.element == "H"]
 83:   if len(hydrogens) < 2: 48:   if len(hydrogens) < 2:
 84:     tree = nx.bfs_tree(atoms, chiral_atom) 49:     tree = nx.bfs_tree(atoms, chiral_atom)
 85:     # Generate the list of shortest paths in the molecule, neglecting the trivial path [chiral_atom] 50:   # Generate the list of shortest paths in the molecule, neglecting the trivial path [chiral_atom]
 86:     paths = sorted(nx.single_source_shortest_path(tree, chiral_atom, depth).values(), reverse = True)[:-1] 51:     paths = sorted(nx.single_source_shortest_path(tree, chiral_atom, depth).values(), reverse = True)[:-1]
 87:     while paths: 52:     while paths:
 88:       # Pop the first element (highest priority path) from the list of paths and remove any duplicates. 53:     # Pop the first element (highest priority path) from the list of paths and remove any duplicates.
 89:       path = paths.pop(0) 54:       path = paths.pop(0)
 90:       paths_no_dups = [unpruned for unpruned in paths if unpruned != path] 55:       paths_no_dups = [unpruned for unpruned in paths if unpruned != path]
 91:       # If there are any duplicates, the paths list will be smaller and we can't resolve a highest priority yet. 56:     # If there are any duplicates, the paths list will be smaller and we can't resolve a highest priority
 92:       if len(paths_no_dups) != len(paths): 57:       if len(paths_no_dups) != len(paths):
 93:         paths = paths_no_dups 58:         paths = paths_no_dups
 94:       # Otherwise, the path is higher priority than all the other paths, so its second atom is the neighbour with 59:       # Otherwise, the path is higher priority than all the other paths, so its second atom is the neighbour with
 95:       # highest priority. 60:       # highest priority.
 96:       else: 61:       else:
 97:         ranked_atom = path[1] 62:         ranked_atom = path[1]
 98:         ordered.append(ranked_atom) 63:         ordered.append(ranked_atom)
 99:         # Drop all the paths containing our ranked atom. 64:       # Drop all the paths containing our ranked atom.
100:         paths = [unpruned for unpruned in paths if unpruned[1] is not ranked_atom] 65:         paths = [unpruned for unpruned in paths if unpruned[1] is not ranked_atom]
101:   return ordered 66:   return ordered
102:  67: 
103: def write_chirality_file(input_filename, output_filename, human_readable_filename): 68: def write_chirality_file(input_filename, output_filename):
104:   molecule = ra.parse_topology_file(input_filename) 69:   molecule = ra.parse_topology_file(input_filename)
105:   atoms = molecule.atoms 70:   atoms = molecule.atoms
106:   chiral_cands = chiral_candidates(atoms) 71:   chiral_cands = chiral_candidates(atoms)
107:   multi_bonds(atoms) 72:   multi_bonds(atoms)
108:   chiral_centres = {} 73:   chiral_centres = {}
109:   for i, chiral_atom in enumerate(chiral_cands): 74:   for i, chiral_atom in enumerate(chiral_cands):
110:     ordered = chiral_order(atoms, chiral_atom) 75:     ordered = chiral_order(atoms, chiral_atom)
111:     if len(ordered) == 4: 76:     if len(ordered) == 4:
112:       chiral_centres[chiral_atom] = ordered 77:       chiral_centres[chiral_atom] = ordered
113:   with open(output_filename, "w") as output_file: 78:   with open(output_filename, "w") as output_file:
114:     for atom in sorted(chiral_centres.keys(), cmp=lambda x, y: cmp(x.index, y.index)): 79:     for atom in sorted(chiral_centres.keys(), cmp=lambda x, y: cmp(x.index, y.index)):
115:       # Write out the list of chiral atoms and their CIP-ranked neighbours. 80:     # Write out the list of chiral atoms and their CIP-ranked neighbours.
116:       output_string = "{0:>8d}{1:>8d}{2:>8d}{3:>8d}{4:>8d}\n".format(atom.index + 1,  81:       output_string = "{0:>8d}{1:>8d}{2:>8d}{3:>8d}{4:>8d}\n".format(atom.index + 1, 
117:                                                                      *[other_atom.index + 1 for other_atom in chiral_centres[atom]]) 82:                                                                      *[other_atom.index + 1 for other_atom in chiral_centres[atom]])
118:       output_file.write(output_string) 83:       output_file.write(output_string)
119:   with open(human_readable_filename, "w") as human_file: 
120:     human_file.write("Atom names given below are in the format: <index (from 0)> <element> <atom name>\n") 
121:     human_file.write("e.g.: 12 C CA is the 13th atom, is a carbon and is named \"CA\"\n\n") 
122:     output_string = "{0:^16s}{1:^16s}{2:^16s}{3:^16s}{4:^16s}\n".format("central atom", "top ranked", "2nd ranked", "3rd ranked", "lowest ranked") 
123:     human_file.write(output_string) 
124:     output_string = "{0:^16s}{1:^16s}{2:^16s}{3:^16s}{4:^16s}\n".format("============", "==========", "==========", "==========", "=============") 
125:     human_file.write(output_string) 
126:     for atom in sorted(chiral_centres.keys(), cmp=lambda x, y: cmp(x.index, y.index)): 
127:       # Write out the list of chiral atoms and their CIP-ranked neighbours - but this time readable for humans. 
128:       output_string = "{0:^16s}{1:^16s}{2:^16s}{3:^16s}{4:^16s}\n".format(str(atom), *[str(other_atom) for other_atom in chiral_centres[atom]]) 
129:       human_file.write(output_string) 
130:  84: 
131: if __name__ == "__main__": 85: if __name__ == "__main__":
132:   write_chirality_file(sys.argv[1], ".chirality_list", ".chirality_readable") 86:   write_chirality_file(sys.argv[1], ".chirality_list")


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0