hdiff output

r28072/cistrans.py 2016-07-06 15:37:46.808999174 +0100 r28071/cistrans.py 2016-07-06 15:37:47.161003902 +0100
 15:              "Se": 2, 15:              "Se": 2,
 16:              "P": 5 } 16:              "P": 5 }
 17:  17: 
 18: def peptide_bonds(atoms): 18: def peptide_bonds(atoms):
 19:   peptide_bond_list = [] 19:   peptide_bond_list = []
 20:   amide_carbons = [] 20:   amide_carbons = []
 21: # Identify amide carbons. 21: # Identify amide carbons.
 22:   for carbon in [atom for atom in atoms if atom.element is "C"]: 22:   for carbon in [atom for atom in atoms if atom.element is "C"]:
 23:     if isinstance(carbon, GhostAtom): 23:     if isinstance(carbon, GhostAtom):
 24:       continue 24:       continue
  25:     if carbon.residue.name == 'PRO':
  26:       continue
 25:     neighbors = atoms.neighbors(carbon) 27:     neighbors = atoms.neighbors(carbon)
 26:     nitrogens = [neighbor for neighbor in neighbors if neighbor.element is "N"] 28:     nitrogens = [neighbor for neighbor in neighbors if neighbor.element is "N"]
 27:     oxygens = [neighbor for neighbor in neighbors if neighbor.element is "O" and  29:     oxygens = [neighbor for neighbor in neighbors if neighbor.element is "O" and 
 28:                not isinstance(neighbor, GhostAtom)] 30:                not isinstance(neighbor, GhostAtom)]
 29:     ghost_oxygens = [neighbor for neighbor in neighbors if neighbor.element is "O" and  31:     ghost_oxygens = [neighbor for neighbor in neighbors if neighbor.element is "O" and 
 30:                      isinstance(neighbor, GhostAtom)] 32:                      isinstance(neighbor, GhostAtom)]
 31:     if len(nitrogens) == 1 and len(oxygens) == 1 and len(ghost_oxygens) == 1: 33:     if len(nitrogens) == 1 and len(oxygens) == 1 and len(ghost_oxygens) == 1:
 32:       amide_carbons.append(carbon) 34:       amide_carbons.append(carbon)
 33: # Identify the O - C - N - H atoms. 35: # Identify the Ca - C - N - Ca chain.
 34:   for carbon in amide_carbons: 36:   for carbon in amide_carbons:
 35:     neighbors = atoms.neighbors(carbon) 37:     neighbors = atoms.neighbors(carbon)
 36:     oxygen = [neighbor for neighbor in neighbors if neighbor.element is "O"][0] 38:     c_alpha_1 = [neighbor for neighbor in neighbors if neighbor.element is "C"][0]
 37:     nitrogen = [neighbor for neighbor in neighbors if neighbor.element is "N"][0] 39:     nitrogen = [neighbor for neighbor in neighbors if neighbor.element is "N"][0]
 38:     n_neighbors = atoms.neighbors(nitrogen) 40:     n_neighbors = atoms.neighbors(nitrogen)
 39:     hydrogen = [neighbor for neighbor in n_neighbors if neighbor.element is "H"][0] 41:     c_alpha_2 = [neighbor for neighbor in n_neighbors if list(set(atoms.neighbors(neighbor)) & set(amide_carbons))]
 40:     peptide_bond_list.append((oxygen, carbon, nitrogen, hydrogen)) 42:     if len(c_alpha_2) == 1:
  43:       peptide_bond_list.append((c_alpha_1, carbon, nitrogen, c_alpha_2[0]))
 41:   return peptide_bond_list 44:   return peptide_bond_list
 42:  45: 
 43: def multi_bonds(atoms): 46: def multi_bonds(atoms):
 44:   multibonded = [atom for atom in atoms.nodes()  47:   multibonded = [atom for atom in atoms.nodes() 
 45:                  if len(nx.edges(atoms, atom)) < valences[atom.element]] 48:                  if len(nx.edges(atoms, atom)) < valences[atom.element]]
 46:   for i, atom in enumerate(multibonded): 49:   for i, atom in enumerate(multibonded):
 47:     paired = False 50:     paired = False
 48:     for other in atoms.neighbors(atom): 51:     for other in atoms.neighbors(atom):
 49:       if isinstance(other, GhostAtom): 52:       if isinstance(other, GhostAtom):
 50:         paired = True 53:         paired = True
 61: #      print "Could not find pair for atom", atom, atom.residue, atom.charge, nx.edges(atoms, atom) 64: #      print "Could not find pair for atom", atom, atom.residue, atom.charge, nx.edges(atoms, atom)
 62:    65:   
 63: def write_cis_trans_file(input_filename, output_filename): 66: def write_cis_trans_file(input_filename, output_filename):
 64:   molecule = ra.parse_topology_file(input_filename) 67:   molecule = ra.parse_topology_file(input_filename)
 65:   atoms = molecule.atoms 68:   atoms = molecule.atoms
 66: #  for atom in atoms: 69: #  for atom in atoms:
 67: #    print atom.residue 70: #    print atom.residue
 68:   multi_bonds(atoms) 71:   multi_bonds(atoms)
 69:   with open(output_filename, "w") as output_file: 72:   with open(output_filename, "w") as output_file:
 70:     for bond in sorted(peptide_bonds(atoms), cmp=lambda x, y: cmp(x[0].index, y[0].index)): 73:     for bond in sorted(peptide_bonds(atoms), cmp=lambda x, y: cmp(x[0].index, y[0].index)):
 71:     # Write out the list of atoms in peptide bonds (O - C - N - H). 74:     # Write out the list of atoms in peptide bonds (Ca - C - N - Ca).
 72:       output_string = "{0:>8d}{1:>8d}{2:>8d}{3:>8d}\n".format(*map(lambda x: x.index + 1, bond)) 75:       output_string = "{0:>8d}{1:>8d}{2:>8d}{3:>8d}\n".format(*map(lambda x: x.index + 1, bond))
 73:       output_file.write(output_string) 76:       output_file.write(output_string)
 74:  77: 
 75: if __name__ == "__main__": 78: if __name__ == "__main__":
 76:   write_cis_trans_file(sys.argv[1], ".cis_trans_list") 79:   write_cis_trans_file(sys.argv[1], ".cis_trans_list")


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0