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:

Erik Mallinson