Monthly Archives: 3月 2015

RDKit: pyAvalonTools


今回は、RDKitから利用できるpyAvalonToolsを使ってみます。
まずは、以下のモジュールと分子を読み込みます。

> python
>>> from rdkit import Chem
>>> from rdkit.Avalon import pyAvalonTools
>>> from rdkit.Chem import Draw
>>> mol = Chem.MolFromSmiles('Cc1c(cc(cc1)NC(=O)c1ccc(cc1)CN1CCN(CC1)C)Nc1nccc(n1)c1cnccc1')

Fingerprintの生成

GetAvalonFPを使って、ここでは128-bitのFingerprintを生成しています。
また、GetAvalonFPAsWordsを使うと、32-bitの整数型でFingerprintを得ることもできます。

>>> fp = pyAvalonTools.GetAvalonFP(mol,128)
>>> print fp.ToBitString()
10111101011010111110100111000101111101110111100011110110110011110011110110101011111011111001111110101101111110011010111111111100
>>> fpw = pyAvalonTools.GetAvalonFPAsWords(mol,128)
>>> print fpw
[4294967229L, 4294967279L, 4294967228L, 4294967221L]

Canonical smilesの生成

rdkitとpyAvalonToolsでcanonical smilesを出力してみます。
同じ分子のsmilesですが、文字列としては、異なったsmilesを出力しています。
smilesの規範化は、利用するソフトウェアによって、異なることが分かります。

>>> cansmi_ava = pyAvalonTools.GetCanonSmiles(mol)
>>> cansmi_rdk = Chem.MolToSmiles(mol)
>>> print cansmi_ava
>>> print cansmi_rdk
Cc1ccc(cc1Nc2nccc(n2)c3cccnc3)NC(=O)c5ccc(CN4CCN(C)CC4)cc5
Cc1ccc(NC(=O)c2ccc(CN3CCN(C)CC3)cc2)cc1Nc1nccc(-c2cccnc2)n1

2D構造の生成

rdkitと同様に、pyAvalonToolsでも2D構造の生成ができます。

>>> pyAvalonTools.Generate2DCoords(mol)
>>> Draw.MolToFile(mol,'ava.png',size=(200,150))

ava


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

APIの詳細:
pyAvalonTools

scikit-learnで機械学習


今回は、scikit-learnを使って機械学習を行ってみます。

scikit-learnのインストール

ここでは、Windows8にscikit-learnをインストールします。

> pip install scikit-learn

Scipyをインストールしていない場合は、以下のサイトからインストーラを取得し、インストールします。
http://sourceforge.net/projects/scipy/files/scipy/0.12.0/scipy-0.12.0-win32-superpack-python2.7.exe/download

学習用データとテスト用データの読み込み

学習用とテスト用のデータを読み込みます。
Fingerprintは、ConvertToNumpyArrayを使ってnumpyの配列に変換しています。

>>> from rdkit import Chem
>>> from rdkit import DataStructs
>>> from rdkit.Chem import AllChem
>>> from sklearn.ensemble import RandomForestClassifier
>>> import numpy as np
>>> 
>>> ms_train = [m for m in Chem.SDMolSupplier('train.sdf')]
>>> fps_train = [AllChem.GetMorganFingerprintAsBitVect(m,2) for m in ms_train]
>>> nfp_train = []
>>> for fp in fps_train:
...     nfp = np.zeros((1,))
...     DataStructs.ConvertToNumpyArray(fp,nfp)
...     nfp_train.append(nfp)
>>> 
>>> ms_pred = [m for m in Chem.SDMolSupplier('pred.sdf')]
>>> fps_pred = [AllChem.GetMorganFingerprintAsBitVect(m,2)  for m in ms_pred]
>>> nfp_pred = []
>>> for fp in fps_pred:
...     nfp = np.zeros((1,))
...     DataStructs.ConvertToNumpyArray(fp,nfp)
...     nfp_pred.append(nfp)

sci_train
sci_test

予測モデルの構築

ここでは、学習用データを用いて、Random Forestにより予測モデルの構築を行っています。
教師データとして、アミノ酸を”0″、核酸を”1″として、与えています。

rf = RandomForestClassifier(n_estimators=100, random_state=1)
act_train = [0,0,0,0,1,1,1,1]
rf.fit(nfp_train,act_train)

予測モデルの評価

テスト用データに含まれている化合物がアミノ酸”0″であるか、核酸”1″であるか予測してみます。

>>> for nfp in nfp_pred:
...     pred = rf.predict(nfp)
...     prob = rf.predict_proba(nfp)
...     print pred,prob
[0] [[ 0.77  0.23]]
[1] [[ 0.22  0.78]]

1番目の分子は、アミノ酸”0″、2番目の分子は、核酸”1″と予測されました。

Similarity Mapによる可視化

1番目の分子を予測モデルとSimilarity Mapを用いて可視化してみます。

>>> from rdkit.Chem.Draw import SimilarityMaps
>>> def getScore(fp, function):
...     score=function(fp)
...     return score[0][0]
>>> fig, maxweight = SimilarityMaps.GetSimilarityMapForModel(ms_pred[0], SimilarityMaps.GetMorganFingerprint, lambda fp: getScore(fp, rf.predict_proba))
>>> fig.savefig("m0_amino.png",bbox_inches='tight')

主鎖の部分が、アミノ酸”0″の予測に強く寄与していることが分かります。
m0_amino2


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

参考にしたサイト
RDKit Online Documentation

APIの詳細:
ConvertToNumpyArray
RandomForestClassifier

RDKit: フラグメントからFingerprintの生成


今回は、フラグメントからFingerprintの生成を行ってみます。
まずは、以下のモジュールと分子を2つ読み込みます。

> python
>>> from rdkit import Chem
>>> from rdkit.Chem import FragmentCatalog
>>> import os
>>> ms = [m for m in Chem.SDMolSupplier('input.sdf')]

Pharma2D_input

官能基情報の読み込み

官能基の定義は、FunctionalGroups.txtに記載されています。
例えば、methyl amideは以下のようにSMARTSで定義されています。

-NC(=O)CH3 *-[N;D2]-[C;D3](=O)-[C;D1;H3]

>>> fName = os.path.join(Chem.RDConfig.RDDataDir,'FunctionalGroups.txt')
>>> fparams = FragmentCatalog.FragCatParams(1,4,fName)
>>> print "Number of functional groups: ",fparams.GetNumFuncGroups()
Number of functional groups:  39

フラグメントの生成

入力した2つの分子からフラグメントを生成します。

>>> fcat = FragmentCatalog.FragCatalog(fparams)
>>> fcgen = FragmentCatalog.FragCatGenerator()
>>> for m in ms:
...     fcgen.AddFragsFromMol(m,fcat)
>>> n = fcat.GetNumEntries()
>>> print "Number of fragments: ",n
Number of fragments:  203
>>> for i in range(n):
...     print fcat.GetEntryDescription(i)
CN
C<=O>N
cN
.

フラグメントからFingerprintの生成

フラグメントからFingerprintの生成は、FragFPGeneratorを用いて行います。

>>> fpgen = FragmentCatalog.FragFPGenerator()
>>> fps = [fpgen.GetFPForMol(m,fcat) for m in ms]

共通のビットを算出し、それが表すフラグメントを確認してみます。

>>> fp_common = fps[0]&fps[1]
>>> onbit = list(fp_common.GetOnBits())
>>> for i in onbit:
...     print "Common fragment:",fcat.GetEntryDescription(i)
Common fragment: CN
Common fragment: C<=O>N
.

最後にTanimoto係数を算出してみます。

>>> a = float(fps[0].GetNumOnBits())
>>> b = float(fps[1].GetNumOnBits())
>>> c = float((fps[0]&fps[1]).GetNumOnBits())
>>> tc = c/(a+b-c)
>>> print tc
0.285714285714

今回は、2つの分子からフラグメントを生成しましたが、通常は、複数の分子からフラグメントを生成し、Fingerprintを生成することが多いと思います。


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

APIの詳細:
FragCatParams
FragCatalog
FragCatGenerator
AddFragsFromMol
FragFPGenerator
GetEntryDescription

RDKit: 2D Pharmacophore Fingerprints


今回は、2D Pharmacophore Fingerprintの生成と類似度の算出を行ってみます。
まずは、以下のモジュールと分子を読み込みます。
input.sdfには2つの分子が格納されています。

> python
>>> from rdkit import Chem
>>> from rdkit.Chem.Pharm2D import Generate
>>> from rdkit.Chem.Pharm2D.SigFactory import SigFactory
>>> from rdkit import DataStructs
>>> from rdkit.Chem import ChemicalFeatures
>>> from rdkit import RDConfig
>>> import os
>>> ms = [m for m in Chem.SDMolSupplier('input.sdf')]

Pharma2D_input

薬理作用団の定義の読み込み

ファーマコフォアを形成する薬理作用団の定義されたファイルを読み込みます。

>>> fdefName = os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef')
>>> featFactory = ChemicalFeatures.BuildFeatureFactory(fdefName)

Fingerprint作成のための準備

SigFactoryでFingerprint生成のためのパラメータ設定を行います。
ここでは、薬理作用団の間の距離は、(0,2),(2,5),(5,8)の3つの範囲で分割しています。

>>> sigFactory = SigFactory(featFactory)
>>> sigFactory.SetBins([(0,2),(2,5),(5,8)])
>>> sigFactory.Init()

Fingerprintの生成と類似度の評価

Gen2DFingerprintによりFingerprintを生成し、Tanimoto係数を算出しています。

>>> fps = [Generate.Gen2DFingerprint(m,sigFactory) for m in ms]
>>> sim = DataStructs.TanimotoSimilarity(fps[0],fps[1])
>>> print sim
0.373831775701

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

APIの詳細:
BuildFeatureFactory
SigFactory
Gen2DFingerprint

RDKit: Similarity mapの作成


分子間の類似度は、Fingerprint間のTanimoto係数等により定量化することができますが、具体的に分子のどの領域が類似しているのか知りたいことがあります。Similarity mapを使うと、類似している領域を可視化して、確認することができます。
まずは、以下のモジュールと分子を読み込みます。

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

Similarity mapを作成してみます。

>>> ms = [m for m in Chem.SDMolSupplier('input.sdf')]
>>> 
>>> weights = SimilarityMaps.GetAtomicWeightsForFingerprint(ms[0], ms[1], SimilarityMaps.GetMorganFingerprint)
>>> fig = SimilarityMaps.GetSimilarityMapFromWeights(ms[1], weights, size=(150, 150))
>>> fig.savefig("ms1.png",bbox_inches='tight')
>>> 
>>> weights = SimilarityMaps.GetAtomicWeightsForFingerprint(ms[1], ms[0], SimilarityMaps.GetMorganFingerprint)
>>> fig = SimilarityMaps.GetSimilarityMapFromWeights(ms[0], weights, size=(150, 150))
>>> fig.savefig("ms0.png",bbox_inches='tight')

simMap

weightsの代わりに原子の電荷を用いれば、部分電荷の度合いを可視化することもできます。

>>> AllChem.ComputeGasteigerCharges(ms[0])
>>> charges = [float(atom.GetProp('_GasteigerCharge')) for atom in ms[0].GetAtoms()]
>>> fig = SimilarityMaps.GetSimilarityMapFromWeights(ms[0],charges, size=(150, 150),scale=10)
>>> fig.savefig("ms0_charge.png",bbox_inches='tight')

ms0_charge


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

APIの詳細:
GetAtomicWeightsForFingerprint
GetSimilarityMapFromWeights

RDKit: 化学反応


今回は、ReactionFromSmartsを使って、2つのアミノ酸をペプチド(アミド)結合で繋げてみたいと思います。
まずは、以下のモジュールを読み込みます。

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

アミド結合を記述し、確認のため、ReactionToImageを使って描画してみます。

>>> rxn = AllChem.ReactionFromSmarts('[O:1]=[C:2][OH:3].[NH2:4]>>[O:1]=[C:2][N:4].[O:3]')
>>> img = Draw.ReactionToImage(rxn,subImgSize=(100,100))
>>> img.save('reaction.png')

reaction2

AlaとValをreactantとして入力します。

>>> ala = Chem.MolFromSmiles('NC(C)C(=O)O')
>>> val = Chem.MolFromSmiles('NC(C(C)C)C(=O)O')
>>> ps = rxn.RunReactants((ala,val))
>>> for p in ps:
...     for pp in p:
...         smi = Chem.MolToSmiles(pp)
...         print smi,
...     print
CC(C)C(NC(=O)C(C)N)C(=O)O O

出力として、Ala-Valと水分子が得られます。


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

APIの詳細:
ReactionFromSmarts
ReactionToImage
RunReactants

RDKit: RECAPでフラグメントの生成


今回は、RECAPを用いて、分子のフラグメント化を行ってみます。
まずは、以下のモジュールと分子を読み込みます。

> python
>>> from rdkit import Chem
>>> from rdkit.Chem import Recap
>>> from rdkit.Chem import Draw
>>> mol = Chem.MolFromSmiles('O=c1nc([nH]c2NCC(Nc12)CN(c1ccc(cc1)C(=O)NC(CCC(=O)O)C(=O)O)C=O)N')

RECAP_input

RecapDecomposeによりフラグメント化します。
結果は、木構造で得られます。

>>> node = Recap.RecapDecompose(mol)
>>> print node

葉ノード(フラグメント)を取得し、画像で出力してみます。

>>> leaves = [leaf.mol for leaf in node.GetLeaves().values()]
>>> img = Draw.MolsToGridImage(leaves)
>>> img.save('leaves.png')

RECAP_leaves

全てのノード(フラグメント)を取得し、画像で出力してみます。

>>> all_nodes = [node.mol for node in node.GetAllChildren().values()]
>>> img = Draw.MolsToGridImage(all_nodes)
>>> img.save('all_nodes.png')

RECAP_all_nodes


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

APIの詳細:
RecapDecompose
GetLeaves
GetAllChildren

RDKit: BRICSでバーチャルライブラリの構築


今回は、BRICSを使ってバーチャルライブラリの構築を行ってみます。
まずは、以下のモジュールを読み込みます。

>>> from rdkit import Chem
>>> from rdkit.Chem import AllChem
>>> from rdkit.Chem import BRICS
>>> from rdkit.Chem import rdMolDescriptors
>>> import random
>>> random.seed(1)

分子からフラグメントを生成

ここでは、ChEMBLからBACE阻害剤を取得して、入力ファイルとしています。

>>> ms = [m for m in Chem.SmilesMolSupplier('BACE.smi',delimiter='\t',titleLine=False) if m is not None]
>>> print len(ms)
461
>>> frags = set()
>>> for m in ms:
...     frag = BRICS.BRICSDecompose(m)
...     frags.update(frag)
>>> print frags
set(['[15*]C1CCC2(Cc3ccc([16*])cc3C23N=C(N)N(C)C3=O)CC1', '[16*]c1ccccc1F',...........
>>> print len(frags)
322

フラグメントから分子の生成

分子量が300から500の化合物を10個、ファイル出力しています。

>>> ms_out = BRICS.BRICSBuild([Chem.MolFromSmiles(smi) for smi in frags])
>>> fw = Chem.SDWriter('new_comps.sdf')
>>> nMol = 0
>>> for m in ms_out:
...     m.UpdatePropertyCache(True)
...     mw = rdMolDescriptors.CalcExactMolWt(m)
...     if mw >= 300 and mw < 500:
...         AllChem.Compute2DCoords(m)
...         fw.write(m)
...         nMol += 1
...         if nMol >= 10:
...             break
>>> fw.close()

少ない化合物数ですが、バーチャルライブラリを生成できました。
new_comps


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

APIの詳細:
BRICSDecompose
BRICSBuild

RDKit: Maximum Common Substructure


今回は、Maximum Common Substructure(MCS)の検出を行ってみます。
まずは、以下のモジュールを読み込みます。

> python
>>> from rdkit import Chem
>>> from rdkit.Chem import rdFMCS
>>> from rdkit.Chem import Draw

ここでは、3つの分子の間のMCSの検出を行っています。

>>> ms = [m for m in Chem.SDMolSupplier('input.sdf') if m is not None]
>>> mcs = rdFMCS.FindMCS(ms)
>>> #mcs = rdFMCS.FindMCS(ms,ringMatchesRingOnly=True)
>>> #mcs = rdFMCS.FindMCS(ms,completeRingsOnly=True)
>>> mcs_smarts = mcs.smartsString
>>> mcs_mol = Chem.MolFromSmarts(mcs_smarts)
>>> for i,m in enumerate(ms):
...     match_atoms = m.GetSubstructMatch(mcs_mol)
...     print match_atoms
...     Draw.MolToFile(m,'comp_%d.png' % i,highlightAtoms=match_atoms)
(1, 9, 11, 6, 4, 8, 3, 5, 7, 13, 17, 16, 12, 10, 15)
(17, 14, 13, 12, 3, 6, 2, 4, 5, 9, 11, 10, 8, 7, 15)
(0, 8, 13, 11, 5, 4, 2, 1, 3, 10, 16, 15, 9, 6, 14)

MCSをハイライトすると以下のようになります。
mcs1

FindMCSにringMatchesRingOnly=Trueとして値を渡してみます。
mcs2

FindMCSにcompleteRingsOnly=Trueとして値を渡してみます。
mcs3
オプションによって、MCSが変化することがよく分かると思います。


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

APIの詳細:
FindMCS

MolVSによるによる互変異性体の生成


今回は、MolVS(Molecule Validation and Standardization)を使ってtautomerの生成をしたいと思います。

MolVSのインストール

事前にRDKitのインストールが必要です。
pipを使ってインストールします。

> pip install molvs

MolVSの実行

ここでは、cytosineのtautomerを生成してみます。
cytosine
enumerate_tautomers_smilesで複数のtautomerを生成することができます。
また、canonicalize_tautomer_smilesにより、規範化されたtautomerも生成できます。

>>> from rdkit import Chem
>>> from molvs import canonicalize_tautomer_smiles
>>> from molvs import enumerate_tautomers_smiles
>>> from rdkit.Chem import Draw
>>> 
>>> mol_data = file('cytosine.mol','r').read()
>>> m = Chem.MolFromMolBlock(mol_data)
>>> smi = Chem.MolToSmiles(m)
>>> tautomers = enumerate_tautomers_smiles(smi)
>>> can_tautomer = canonicalize_tautomer_smiles(smi)
>>> print tautomers
set(['N=c1ccnc(O)[nH]1', 'N=c1cc[nH]c(=O)[nH]1', 'Nc1cc[nH]c(=O)n1', 'N=c1cc[nH]c(O)n1', 'Nc1ccnc(=O)[nH]1', 'Nc1ccnc(O)n1'])
>>> print can_tautomer
N=c1cc[nH]c(=O)[nH]1

結果

Tautomers
tautomers
Canonicalize Tautomer
can_tautomers
利用したソフトウェア:

RDKit_2014_09_2
MolVS-0.0.3

APIの詳細:
MolVS_API