Molecular Graph-ics with Pybel

Graphs are great. There are books and books of algorithms written by generations of computer scientists that take graphs and do interesting things. And better still, there are open source programming libraries available that implement many of these algorithms. So, given a Pybel Molecule, how can it be converted for use by these libraries?

After some googling around, I found three graph libraries accessible from Python (on Windows): networkx, igraph and the Boost Graph Library (BGL). Whatever library is used, the solution is pretty much the same; iterate over all of the atoms and bonds of the molecule, and add nodes and edges to the graph.

Here's a function that takes a Pybel Molecule and returns a networkx graph (remember to install networkx first...):

def mol_to_networkxgraph(mol):
edges = []
bondorders = []
for bond in ob.OBMolBondIter(mol.OBMol):
bondorders.append(bond.GetBO())
edges.append( (bond.GetBeginAtomIdx() - 1, bond.GetEndAtomIdx() - 1) )
g = networkx.Graph()
g.add_edges_from(edges)
return g

What about making an igraph graph?

def mol_to_igraph(mol):
edges = []
bondorders = []
for bond in ob.OBMolBondIter(mol.OBMol):
bondorders.append(bond.GetBO())
edges.append( (bond.GetBeginAtomIdx() - 1, bond.GetEndAtomIdx() - 1) )

atomtypes = [atom.type for atom in mol]

g = igraph.Graph(edges=edges,
vertex_attrs={'atomtype':atomtypes},
edge_attrs={'bondorder': bondorders})
return g

And finally, making a BGL graph:

def mol_to_boostgraph(mol):
edges = []
bondorders = []
for bond in ob.OBMolBondIter(mol.OBMol):
bondorders.append(bond.GetBO())
edges.append( (bond.GetBeginAtomIdx() - 1, bond.GetEndAtomIdx() - 1) )

g = boost.graph.Graph(edges)
bondordermap = g.edge_property_map("integer")
for edge, bondorder in zip(g.edges, bondorders):
bondordermap[edge] = bondorder
g.edge_properties["bondorder"] = bondordermap

atomtypemap = g.vertex_property_map("string")
atomtypes = [atom.type for atom in mol]
for vertex, atomtype in zip(g.vertices, atomtypes):
atomtypemap[vertex] = atomtype
g.vertex_properties["atomtype"] = atomtypemap

return g

Now that you have a graph you can use any of the algorithms provided by the library.

Image credit: Creative Type websites as graphs by Erik Mallinson (CC BY-NC-SA 2.0)