Monthly Archives: 2月 2015

RDKit: Atomに対する基礎的な操作


今回は、Atomに対する基礎的な操作を行ってみます。
まずは、以下のモジュールと分子を読み込みます。

> python
>>> from rdkit import Chem
>>> mol_data = file('alanine.mol','r').read()
>>> mol = Chem.MolFromMolBlock(mol_data)

atom_alanine2
まずは、原子名、原子番号等の情報を取得してみます。

>>> for atom in mol.GetAtoms():
...     print "idx= ",atom.GetIdx()
...     print " symbol= ",atom.GetSymbol()
...     print " atomic_num= ",atom.GetAtomicNum()
...     print " mass= ",atom.GetMass()
...     print " hybridization= ",atom.GetHybridization()
...     print " total_degree= ",atom.GetTotalDegree()
...     print " total_valence= ",atom.GetTotalValence()
...     print " formal_charge= ",atom.GetFormalCharge()
...     print " chiral_tag= ",atom.GetChiralTag()
idx=  0
 symbol=  O
 atomic_num=  8
 mass=  15.999
 hybridization=  SP2
 total_degree=  1
 total_valence=  1
 formal_charge=  -1
 chiral_tag=  CHI_UNSPECIFIED
idx=  1
 symbol=  O
 atomic_num=  8
 mass=  15.999
 hybridization=  SP2
 total_degree=  1
 total_valence=  2
 formal_charge=  0
 chiral_tag=  CHI_UNSPECIFIED
idx=  2
 symbol=  N
 atomic_num=  7
 mass=  14.007
 hybridization=  SP3
 total_degree=  4
 total_valence=  4
 formal_charge=  1
 chiral_tag=  CHI_UNSPECIFIED
idx=  3
 symbol=  C
 atomic_num=  6
 mass=  12.011
 hybridization=  SP3
 total_degree=  4
 total_valence=  4
 formal_charge=  0
 chiral_tag=  CHI_UNSPECIFIED
idx=  4
 symbol=  C
 atomic_num=  6
 mass=  12.011
 hybridization=  SP3
 total_degree=  4
 total_valence=  4
 formal_charge=  0
 chiral_tag=  CHI_UNSPECIFIED
idx=  5
 symbol=  C
 atomic_num=  6
 mass=  12.011
 hybridization=  SP2
 total_degree=  3
 total_valence=  4
 formal_charge=  0
 chiral_tag=  CHI_UNSPECIFIED

次に、注目している原子に結合している原子の情報を取得を行っています。

>>> atom = mol.GetAtomWithIdx(3)
>>> for nbr in atom.GetNeighbors():
...     print nbr.GetIdx(),nbr.GetSymbol()
2 N
4 C
5 C

同じことが、原子からでている結合を使ってもできます。

>>> for bond in atom.GetBonds():
...     oatom = bond.GetOtherAtom(atom)
...     print oatom.GetIdx(),oatom.GetSymbol()
2 N
4 C
5 C

分子と同様に、各原子にも任意のプロパティを設定することができます。

>>> atom.SetProp("PropertyA","123")
>>> print atom.GetProp("PropertyA")
123

利用したソフトウェア:
RDKit_2014_09_2
APIの詳細:
Atom
Bond

RDKit: 直列化した分子の読み込み


先日、cPickleを使って分子の直列化を行いました。
今回は、直列化した分子をファイル出力し、それを入力してみます。
まずは、以下のモジュールを読み込みます。

> python
>>> from rdkit import Chem
>>> import cPickle
>>> import time

通常の分子の読み込み

50,784件の分子の読み込みには、約29.3秒かかっています。

>>> start = time.time()
>>> mols = [mol for mol in Chem.SDMolSupplier('input.sdf')]
>>> end = time.time()
>>> print "Input(normal): %d msec" % int((end-start) * 1000)
>>> print len(mols)
Input(normal): 29313 msec
50784

分子の直列化とファイル出力

直列化には、約10.6秒かかっています。

>>> start = time.time()
>>> outf = open('pickle_mol.txt','w')
>>> for mol in mols:
...     pkl = cPickle.dumps(mol,0)
...     outf.write(pkl+"\n$$$$\n")
>>> outf.close()
>>> end = time.time()
>>> print "Output(serialize): %d msec" % int((end-start) * 1000)
Output(serialize): 10667 msec

ファイル入力と分子の非直列化

直列化された分子のファイルを読み込んで、非直列化した場合の時間は、約7.4秒となります。
直列化の時間は必要ですが、約4倍速くなっています。
今回は、テストのために、分子の区切りとして”$$$$”を用いて、ファイル出力しましたが、本来は、1分子ごとにデータベースに格納して利用した方がいいと思います。

>>> start = time.time()
>>> mols = []
>>> pkl = ""
>>> inf = open('pickle_mol.txt')
>>> while 1:
...    line = inf.readline()
...    if not line:
...        break
...    if line[:4] == "$$$$":
...        mol = cPickle.loads(pkl)
...        mols.append(mol)
...        pkl = ""
...    else:
...        pkl = pkl+line
>>> inf.close()
>>> end = time.time()
>>> print  "Input(deselialize): %d msec" % int((end-start) * 1000)
>>> print len(mols)
Input(deselialize): 7440 msec
50784

利用したソフトウェア:
RDKit_2014_09_2

RDKit: Moleculeに対する基礎的な操作


今回は、Moleculeに対する基礎的な操作を行ってみます。
まずは、以下のモジュールと分子を読み込みます。

> python
>>> from rdkit import Chem
>>> mol_data = file('acetaminofen.mol','r').read()
>>> mol = Chem.MolFromMolBlock(mol_data)

原子数、結合の数、そして重原子数の取得

>>> nAtom = mol.GetNumAtoms()
>>> nBond = mol.GetNumBonds()
>>> nHAtom = mol.GetNumHeavyAtoms()
>>> print "nAtom=",nAtom
>>> print "nBond=",nBond
nAtom= 11
nBond= 11
nHAtom= 11

Conformer情報の取得

配座情報を取得してみます。

>>> confs = mol.GetConformers()
>>> print "nConfs=",len(confs)
>>> conf = confs[0]
>>> for idx in range(conf.GetNumAtoms()):
...     pos = conf.GetAtomPosition(idx)
...     print "  %2d" % idx,
...     print " %8.3f%8.3f%8.3f" % (pos.x,pos.y,pos.z)
nConfs= 1
   0     2.866  -2.595   0.000
   1     4.598   1.405   0.000
   2     2.866   1.405   0.000
   3     2.866   0.405   0.000
   4     2.000  -0.095   0.000
   5     3.732  -0.095   0.000
   6     2.000  -1.095   0.000
   7     3.732  -1.095   0.000
   8     2.866  -1.595   0.000
   9     3.732   1.905   0.000
  10     3.732   2.905   0.000

分子の直列化

cPickleを使えばPythonのオブジェクトを直列化することができます。
ここでは、分子情報の入ったmolを直列化しています。
dumpsの第二引数は、直列化のためのプロトコルを示します。
バージョン0より2の方が、サイズが小さくなることがわかります。

>>> import cPickle
>>> pkl1 = cPickle.dumps(mol,0)
>>> print len(pkl1)
1100
>>> pkl2 = cPickle.dumps(mol,2)
>>> print len(pkl2)
318

直列化されたmolを復元してみます。

>>> molpkl1 = cPickle.loads(pkl1)
>>> print Chem.MolToSmiles(molpkl1)
CC(=O)Nc1ccc(O)cc1
>>> molpkl2 = cPickle.loads(pkl2)
>>> print Chem.MolToSmiles(molpkl2)
CC(=O)Nc1ccc(O)cc1

次は、molをbinaryに変換して、復元してみます。

>>> bina = mol.ToBinary()
>>> print len(bina)
318
>>> molbina = Chem.Mol(bina)
>>> print Chem.MolToSmiles(molbina)
CC(=O)Nc1ccc(O)cc1

binaryが一番サイズが小さいことが分かります。


利用したソフトウェア:
RDKit_2014_09_2

APIの詳細:
Mol
Conformer

RDKit: Fraggleで類似構造検索


今回は、Fraggleで類似構造検索を行ってみます。
Fraggleは、RDKit_2014_09_02\Contrib\fraggleフォルダにあるスクリプトを使って実行できます。

Fraggleの特徴ですが、一般的なFingerprintでは、クエリー分子とターゲット分子の骨格において、骨格中の窒素の位置がわずかに異なるだけで、その類似度に大きな影響を及ぼします。

一方、Fraggleは、分子中のフラグメントの一致を重要視することで、骨格中の窒素の位置のわずかな違い等が類似度に大きく反映しないように工夫されています。

以下の資料に詳細が説明がされています。
https://github.com/rdkit/UGM_2013/blob/master/Presentations/Hussain.Fraggle.pdf

Step1. Query Fragmentation

クエリー分子のフラグメント化を行います。
クエリー分子は、SMILESで入力しますが、ここでは、SMILESの後にスペースを1つ入れ、IDをつけています。

query.smi:

c1(nc2c(n1CC)ccnc2)c1ccccc1 query

fraggle_query

python fraggle.py < query.smi > query_frag.csv

4パターンのフラグメントが生成されました。

query_frag.csv:

c1(nc2c(n1CC)ccnc2)c1ccccc1,query,[*]c1nc(-c2ccccc2)n(CC)c1[*]
c1(nc2c(n1CC)ccnc2)c1ccccc1,query,[*]n1c2ccncc2nc1-c1ccccc1
c1(nc2c(n1CC)ccnc2)c1ccccc1,query,[*]c1ccccc1.[*]c1ccncc1[*]
c1(nc2c(n1CC)ccnc2)c1ccccc1,query,[*]c1ccncc1[*]

fraggle_frags2

Step2. Tversky Search

Step1で生成したフラグメントとターゲット分子間の類似度をTverskey係数で算出し、閾値を設定することにより選別を行います。ここでは、閾値を0.6としています。
target.smiの形式は、クエリー分子と同じです。

target.smi:

n1(c(nc2c1ncnc2N)c1cc(F)ccc1)C target01
n1cnc2c(c1N(C)C)nc(n2Cc1ccccc1)C target02
n1cnc2c(c1N1CCCC1)nc([nH]2)c1cnccc1 target03

target

> python rdkit_tversky.py -f query_frag.csv -c 0.6 < target.smi > frag_tversky.out

クエリー分子のfragment1と2が、target01と類似度0.6以上を示しています。
また、クエリー分子のfragment2が、target03と類似度0.6以上を示しています。
target02と類似度0.6以上をもつフラグメントは無かったようです。

frag_tversky.out:

[*]c1nc(-c2ccccc2)n(CC)c1[*],c1(nc2c(n1CC)ccnc2)c1ccccc1,query,n1(c(nc2c1ncnc2N)c1cc(F)ccc1)C,target01,0.651376146789
[*]n1c2ccncc2nc1-c1ccccc1,c1(nc2c(n1CC)ccnc2)c1ccccc1,query,n1(c(nc2c1ncnc2N)c1cc(F)ccc1)C,target01,0.641025641026
[*]n1c2ccncc2nc1-c1ccccc1,c1(nc2c(n1CC)ccnc2)c1ccccc1,query,n1cnc2c(c1N1CCCC1)nc([nH]2)c1cnccc1,target03,0.676923076923

Step3. Post-Processing

Step2.で選別されたクエリー分子のフラグメントを基に、ターゲット分子とターゲット分子間の類似度を計算します。

> python atomcontrib.py -c 0.6 < frag_tversky.out > result.csv

result.csv:

SMILES,ID,QuerySMI,QueryID,Fraggle_Similarity,RDK5_Similarity
n1cnc2c(c1N1CCCC1)nc([nH]2)c1cnccc1,target03,c1(nc2c(n1CC)ccnc2)c1ccccc1,query,0.603773584906,0.386920980926
n1(c(nc2c1ncnc2N)c1cc(F)ccc1)C,target01,c1(nc2c(n1CC)ccnc2)c1ccccc1,query,0.759398496241,0.534591194969

Fraggleで算出された類似度とRDK5の類似度を示しています。
Fingerprintの種類により、ベースラインの類似度の値が異なりますので比較は難しいのですが、以下の化合物間においては、Fraggleの方が、RDK5よりも高い類似度を示しています。

fraggs_result3


利用したソフトウェア:
RDKit_2014_09_2

RDKit: Ring情報の取得


今回は、分子の中のring情報を取得してみます。
まずは、以下のモジュールと分子を読み込みます。

> python
>>> from rdkit import Chem
>>> mol = Chem.MolFromSmiles('C1(N=CC=C2)=C2C=CC=C1')

quinoline2

RingInfoクラスの使い方

RingInfoクラスを使って、各原子が何個のringに属しているのか調べてみます。

>>> rinfo = mol.GetRingInfo()
>>> for idx in range(mol.GetNumAtoms()):
...    print idx,rinfo.NumAtomRings(idx)
0 2
1 1
2 1
3 1
4 1
5 2
6 1
7 1
8 1
9 1

0番と5番の原子は、2つのringに属し、それ以外は1つのringに属していることが分かります。
入力として使ったキノリンは、10個のAtomから構成されています。
従って、NumAtomRingsの引数は、0~9までが対象となるのですが、10以上の値を入れてもエラーとはならず、0が戻ってきますので、注意が必要です。

分子の中に任意のサイズのringが存在するか否かも確認できます。
ringサイズは、範囲で指定できます。

>>> print rinfo.IsAtomInRingOfSize(3,6)
True
>>> print rinfo.IsAtomInRingOfSize(7,12)
False

ここでは、Smallest rings(最小環)を判定の対象としています。
従って、キノリンには、ringサイズが3から6のringはありますが、7から12は、無いことが分かります。
構造式を見て分かるように、この結果は、キノリン骨格の中に6員環が存在することに起因しています。

原子ごとにring情報の取得

RingInfoクラスを使わなくても、ring情報の取得はできます。

>>> for atom in mol.GetAtoms():
...     print atom.GetIdx(),atom.IsInRing(),atom.IsInRingSize(5),atom.IsInRingSize(6)
0 True False True
1 True False True
2 True False True
3 True False True
4 True False True
5 True False True
6 True False True
7 True False True
8 True False True
9 True False True

Smallest Set of Smallest Rings(SSSR)の取得

SSSRを取得してみます。

>>> sssr = Chem.GetSymmSSSR(mol)
>>> for ring in sssr:
...     print list(ring)
[1, 2, 3, 4, 5, 0]
[6, 7, 8, 9, 0, 5]

SSSRとして、6員環が2つ見つかりました。


利用したソフトウェア:
RDKit_2014_09_2

APIの詳細:
GetRingInfo
GetSymmSSSR

RDKit: 水素原子の付加と削除


今回は水素原子の付加と削除を行ってみます。
AddHsとRemoveHsを使います。

hatom.py:

from rdkit import Chem

def countH(mol):
    """ H原子のカウント
    """
    Hs = 0
    for atom in mol.GetAtoms():
        if atom.GetAtomicNum() == 1:
            Hs += 1

    return Hs

# SMILESをMolに変換
mol = Chem.MolFromSmiles('Oc1ccc(NC(=O)C)cc1')

# H原子のカウントと表示
print Chem.MolToSmiles(mol)
Hs=countH(mol)
print "Num of Hs=",Hs

# H原子の付加
print "\nAdd Hs"
mol = Chem.AddHs(mol)
print Chem.MolToSmiles(mol)

# H原子のカウントと表示
Hs=countH(mol)
print "Num of Hs=",Hs

# H原子の削除
print "\nRemove Hs"
mol=Chem.RemoveHs(mol)
print Chem.MolToSmiles(mol)

# H原子のカウントと表示
Hs=countH(mol)
print "Num of Hs=",Hs
> python hatom.py
CC(=O)Nc1ccc(O)cc1
Num of Hs= 0

Add Hs
[H]Oc1c([H])c([H])c(N([H])C(=O)C([H])([H])[H])c([H])c1[H]
Num of Hs= 9

Remove Hs
CC(=O)Nc1ccc(O)cc1
Num of Hs= 0

利用したソフトウェア:
RDKit_2014_09_2

APIの詳細:
AddHs
RemoveHs

RDKit: プロパティの操作


今回は、sdファイル内のプロパティの情報を操作してみます。
化合物名や、物性値等、様々な情報を分子に持たせることができるので、よく使う機能だと思います。

まずは、以下のモジュールを読み込みます。

> python
>>> from rdkit import Chem

input.sdfを読み込んで、sdファイルのタイトル行を表示してみます。
GetPropを用いて、プロパティ名は’_Name’で取得できます。

>>> mols = [mol for mol in Chem.SDMolSupplier('input.sdf') if mol is not None]
>>> for mol in mols:
>>>     print mol.GetProp('_Name')
comp1
comp2
comp3

どのようなプロパティがあるのかを知りたい場合は、GetPropNamesを用います。

>>> mol = mols[0]
>>> for key in mol.GetPropNames():
>>>     print key,mols[0].GetProp(key)
EXACT_MASS 151.063329
INCHI InChI=1S/C8H9NO2/c1-6(10)9-7-2-4-8(11)5-3-7/h2-5,11H,1H3,(H,9,10)
INCHIKEY RZVAJINKPMORJF-UHFFFAOYSA-N
IUPAC__NAME N-(4-hydroxyphenyl)acetamide
XLOGP 0.5

ClearPropにより、プロパティは削除することもできます。
また、SetPropにより、新しいプロパティを設定することもできます。

>>> mol.ClearProp('INCHI')
>>> print mol.SetProp('SMILES','CC(=O)Nc1ccc(O)cc1')
CC(=O)Nc1ccc(O)cc1


利用したソフトウェア:
RDKit_2014_09_2

APIの詳細:
GetProp
GetPropNames
ClearProp
SetProp

RDKit: 分子の出力


今回は、分子の出力を行ってみます。
まずは以下のモジュールをインポートします。

> python
>>> from rdkit import Chem
>>> from rdkit.Chem import AllChem

1つの分子のファイル出力

ここでは、Compute2DCoordsを使って2Dの座標データを取得しています。
MolToMolBlockでmol形式のテキストを生成し、ファイル出力しています。

>>> mol = Chem.MolFromSmiles('c1ccccc1OC')
>>> AllChem.Compute2DCoords(mol)
>>> mol_data = Chem.MolToMolBlock(mol)
>>> file('output.mol','w').write(mol_data)

複数の分子のファイル出力

sd形式で複数の分子を出力するためには、SDWriterを使います。

>>> mol1 = Chem.MolFromSmiles('CC(=O)Nc1ccc(O)cc1')
>>> mol2 = Chem.MolFromSmiles('CC(=O)Nc1ccccc1')
>>> mols = [mol1,mol2]
>>> outf = Chem.SDWriter('output.sdf')
>>> for mol in mols:
...    AllChem.Compute2DCoords(mol)
...    outf.write(mol)

利用したソフトウェア:
 RDKit_2014_09_2
APIの詳細:
 SDWriter
 Compute2DCoords

RDKit: 分子の入力


今回は、RDKitで分子の入力を行ってみます。
まずは、以下のモジュールをインポートします。

> python
>>> from rdkit import Chem
>>> import gzip

1つの分子の入力

ファイル(input.mol)を読み込んで、原子数を表示しています。

>>> mol = Chem.MolFromMolFile('input.mol')
>>> print mol.GetNumAtoms()
11

ファイル(input.mol)を通常のテキストファイルと同じ方法で読み込み、Mol型に変換することもできます。

>>> mol_data = file('input.mol','r').read()
>>> mol = Chem.MolFromMolBlock(mol_data)
>>> print mol.GetNumAtoms()
11

複数の分子の入力

sdファイルを読み込む方法です。ここでは、3つの分子がsdファイルに収まっています。
また、不適切な分子の場合、Noneが返されるため、Noneではない分子だけ原子数を表示しています。
ちなみに、ランダムアクセスもできます。

>>> sdsuppl = Chem.SDMolSupplier('input.sdf')
>>> for mol in sdsuppl:
...    if mol is not None:
...       print mol.GetNumAtoms()
11
10
8
>>> sdsuppl[2].GetNumAtoms()
8

リスト内包表記を使えば、シンプルになります。

>>>sdsuppl = Chem.SDMolSupplier('input.sdf')
>>>mols = [mol for mol in sdsuppl if mol is not None]

MolFromMolBlockを使っても複数分子の読み込みはできます。
ただし、プロパティの情報は、読み込まれません。

>>> inf = open('input.sdf','r')
>>> mol_data =""
>>> while 1:
...    line = inf.readline()
...    if not line:
...       break
...    if line[:4] == "$$$$":
...       mol = Chem.MolFromMolBlock(mol_data)
...       print Chem.MolToSmiles(mol)
...       mol_data = ""
...    else:
...       mol_data = mol_data+line
CC(=O)Nc1ccc(O)cc1
CC(=O)Nc1ccccc1
Oc1ccc(O)cc1

大きなsdファイルを読み込む場合、特に分子に対し逐次的に何らかの計算を行う場合は、ForwardSDMolSupplierが適していると思います。

>>> for mol in Chem.ForwardSDMolSupplier('input.sdf'):
...    if mol is not None:
...       print mol.GetNumAtoms()
11
10
8

gzで圧縮されたsdfを読み込みたい場合にも使えます。

>>> inf = gzip.open('input.sdf.gz')
>>> fsdsuppl = Chem.ForwardSDMolSupplier(inf)
>>> mols = [mol for mol in fsdsuppl if mol is not None]

複数のSMILESが格納されたsmiファイルを読み込むこともできます。
ここでは、タイトル行なし、タブ区切りの以下のファイルを読み込んでみます。

input.smi:

Oc1ccc(NC(=O)C)cc1 compA
O=C(Nc1ccccc1)C compB
Oc1ccc(O)cc1 compC

>>> smisuppl = Chem.SmilesMolSupplier('input.smi',delimiter='\t',titleLine=False)
>>> mols = [mol for mol in smisuppl if mol is not None]


利用したソフトウェア:
RDKit_2014_09_2

APIの詳細:
MolFromMolFile
MolFromMolBlock
MolFromSmiles
SDMolSupplier
ForwardSDMolSupplier
SmilesMolSupplier

RDKitのインストール (Windows)


RDKitは、代表的なChemoinformaticsのツールキットです。
RDKitを利用すれば、創薬に関連する多くのアプリケーションを迅速に開発することができます。

今回は、Windows8にRDKitをインストールしたいと思います。

RDKitに必要なソフトウェアのインストール

・Python2.7のインストール

入手先:https://www.python.org/
ファイル名:python-2.7.9.msi

入手したファイルをダブルクリックし、インストーラの指示に従ってください。
インストールが終わったら、環境変数PATHに以下を追加してください。

C:\Python27;C:\Python27\Scripts

Scriptsフォルダには、後で利用するpipとeasy_installがインストールされています。

・Microsoft Visual C++ Compiler for Python 2.7 のインストール

入手先:http://www.microsoft.com/en-us/download/details.aspx?id=44266
ファイル名:VCForPython27.msi

入手したファイルをダブルクリックし、インストーラの指示に従ってください。
本ソフトウェアは、NumPyのインストールの際に、必要となります。

・NumPyのインストール
pipを使ってインストールします。

> pip install numpy

・Pillow(PIL)のインストール
easy_installを使ってインストールします。

> easy_install pillow

RDKitのインストール

入手先:http://sourceforge.net/projects/rdkit/files/
ファイル名:RDKit_2014_09_2.win32.py27.zip

ファイルを展開して、RDKit_2014_09_2フォルダを任意の場所に置きます。
ここでは、Cドライブ直下におき、C:\RDKit_2014_09_2としました。

環境変数の設定を行います。

RDBASEに以下を設定

C:\RDKit_2014_09_2

PYTHONPATHに以下を設定

%RDBASE%

PATHに以下を追加

;%RDBASE%\lib

動作確認

それでは、簡単なプログラムで動作を確認してみます。
ここでは、分子をSMILESで入力し、comp.pngとしてファイル出力をしています。

> python
>>> from rdkit import Chem
>>> from rdkit.Chem import Draw
>>> mol = Chem.MolFromSmiles(‘c1ccccc1O’)
>>> Draw.MolToImageFile(mol,’comp.png’,size=(150,150))

comp.png:
comp