Convert rdkit molecule object to igraph graph object.

Molecules are often handled as graph in chemoinformatics.
There are some libraries for graph analysis in python. Today, I wrote a sample script that convert from molecule to graph.
I used python-igraph and rdkit.
RDkit has method to get adjacency matrix from molecule so, I used the method.
Code is following.

from rdkit import Chem
from rdkit.Chem import rdmolops
import igraph
import numpy as np
# mol2graph function can convert molecule to graph. And add some node and edge attribute.
def mol2graph( mol ):
    admatrix = rdmolops.GetAdjacencyMatrix( mol )
    bondidxs = [ ( b.GetBeginAtomIdx(),b.GetEndAtomIdx() ) for b in mol.GetBonds() ]
    adlist = np.ndarray.tolist( admatrix )
    graph = igraph.Graph()
    g = graph.Adjacency( adlist ).as_undirected()
    for idx in g.vs.indices:
        g.vs[ idx ][ "AtomicNum" ] = mol.GetAtomWithIdx( idx ).GetAtomicNum()
        g.vs[ idx ][ "AtomicSymbole" ] = mol.GetAtomWithIdx( idx ).GetSymbol()
    for bd in bondidxs:
        btype = mol.GetBondBetweenAtoms( bd[0], bd[1] ).GetBondTypeAsDouble()
        g.es[ g.get_eid(bd[0], bd[1]) ][ "BondType" ] = btype
        print( bd, mol.GetBondBetweenAtoms( bd[0], bd[1] ).GetBondTypeAsDouble() )
    return g

Now test it.

# Oseltamivir ( tamiful )
mol2 = Chem.MolFromSmiles( "CCOC(=O)C1=C[C@@H](OC(CC)CC)[C@H](NC(C)=O)[C@@H](N)C1" )
g2=mol2graph(mol2)
(0, 1) 1.0
(1, 2) 1.0
(2, 3) 1.0
(3, 4) 2.0
(3, 5) 1.0
(5, 6) 2.0
(6, 7) 1.0
....

Seems work fine. 😉
Next, get bond infromation.

for e in g2.es:
    print(e)
igraph.Edge(<igraph.Graph object at 0x10ae999a8>, 0, {'BondType': 1.0})
igraph.Edge(<igraph.Graph object at 0x10ae999a8>, 1, {'BondType': 1.0})
igraph.Edge(<igraph.Graph object at 0x10ae999a8>, 2, {'BondType': 1.0})
......

OK next, I tried to plot tamiful as Graph.
Igraph has plot function.
I added some options for plotting.

from igraph import GraphBase
layout = g2.layout_graphopt()
color_dict = {"C": "gray", "N": "blue", "O" :  "white" }
my_plot = igraph.Plot()
my_plot.add(g2,layout=layout,bbox=(400,400),
             margin=90,
             vertex_color = [color_dict[atom] for atom in g2.vs["AtomicSymbole"]],
             vertex_size = [ v.degree()*10 for v in g2.vs ])
my_plot.show()

tamiful

Ummmm @_@ not so good view.
But fun!

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s