ケムインフォマティクスに虚空投げ

Fingerprintについてメモしたいと思います。化学構造のFingerprintでは、ビット列中の各ビットに部分構造(フラグメント)があてはめられています。例えば、1ビット目にあてはめられた部分構造が、化合物中に存在すれば"1"、存在しなければ"0"とセットします。Fingerprintの長さは、通常、数百〜数千ビット程度ではないでしょうか。

代表的なFingerprintの使い方は、類似性検索だと思います。分子Aと分子Bの類似性を評価する場合、Fingerprint A(分子Aより生成)とFingerprint B(分子Bより生成)間の距離を"類似性"の尺度とし評価します。この距離の計算には、Tanimoto係数、コサイン係数など多くの手法が提案されています。

もう1つのFingerprintの作り方として、ハッシングを用いる方法があります。上記方法では、ユーザが各ビットに特定の部分構造を設定しますが、ハッシングでは、与えられた部分構造からハッシュ値を計算し、その値で、何ビット目に1を立てるか決定します。利点としては、はじめに各ビットの部分構造を決定する必要がないこと。また、ハッシュ関数でハッシュ値の上限が設定できるため、ビット列長を自由に設定できることなどが挙げられます。逆に、欠点としては、ハッシュ値の衝突により、特定のビットに2つ以上の部分構造があてはめられることが挙げられます。したがって、類似性評価を行う場合、Fingerprint生成アルゴリズムはチェックすべき項目だと思います。

さて、部分構造検索は、グラフ理論でいうグラフの同形判定であり、この同形を判定する作業は、NP完全問題です。したがって、この部分における高速化は、かなりしんどいタスクとなります。そこで、同形判定をする前に、明らかに部分構造ではない化合物を取り除く工程でFingerprintが利用されています。ここでは、ハッシングによるFingerprintが活躍しています。

例えば、Fingerprint Aをクエリーとして、Fingerprint Bを標的とする場合、Fingerprint AがFingerprint Bの部分構造であるためには、Fingerprint Aの全ての1の立っているビットがBにおいても1である必要があります。この計算はビット演算子を用いて評価できますので、非常に高速に計算できます。そうは言っても、ハッシングによるFingerprintでは衝突が発生しますが、その後に同形判定は必ず行われますので、最終結果には影響を及ぼさないことになります。

次回、Frownsを使って、Fingerprint + SMARTSによる部分構造検索プログラムを作成したいと思います。



人気ブログランキング(クリックして応援してね)

昨日の記事で、CDKを用いてMDL mol形式の入出力を行いました。出力されたmolファイルを分子Viewerに読み込ませても表示されないという方がいるかもしれません。

この原因は、CDKの出力するmol形式のBond Blockにあります。CDKではBond Blockは例のように第7fieldまで出力します。

例：
2 1 1 0 0 0 0

私の記憶では、この第7fieldは“reacting center status”だったと思いますので、CDKの出力は妥当だと思うのですが...?(仕様書確認しないと)

一部の分子Viewerでは、この第7fieldがあると分子が表示されないことがあるようですので、解決策として、以下のように第7fieldを単純に削除すればOKです。

例：
2 1 1 0 0 0

この作業を全て手作業でやるのは大変ですので、スクリプトを書くか、もしくは、babelでmol形式からmol形式に変換することにより解決できます。

$ babel -imol input.mol -omol output.mol

babelも第7field出さないみたい....




人気ブログランキング(クリックして応援してね)

CDKを使って分子ファイルの入出力を行いたいと思います。まずは、ビルドされたCDKのjarファイルを取得し、クラスパスに追加してください。

http://sourceforge.net/project/showfiles.php?group_id=20024

Eclipseを利用されている方は、[Project]->[Properties]->[Java Build Path]を選択->[Libraries]タブをクリック->[Add External JARs]をクリックし、cdk-20050826.jarを追加してください。

今回は、MDL molファイルを入力し、単純に出力することにします。インラインフレームにソースコードを載せていますので、参考にしてください。
molファイルの入力には、MDLReaderクラスを利用します。ここで、ちょっと分かりにくい部分は、readメソッドの使い方だと思います。readメソッドは、引数にChemObjectのサブクラスのオブジェクトを渡す必要があります。今回はMoleculeクラス(ChemObjectのサブクラスの1つです)を使いますので、以下のように記述します。

mol = (Molecule)mr.read(new Molecule());

出力はMDLWriterクラスを使います。

mw.writeMolecule(mol);

その他の分子フォーマットの入出力を行いたい場合、例えば、PDB形式であれば、PDBReader, PDBWriterクラスを使えばできます。利用法もほぼ同じですので、必要に応じてXXXReader, XXXWriterクラスを探してみてください。
CDK API: http://cdk.sourceforge.net/api/






人気ブログランキング(クリックして応援してね)

Frownsを用いて化合物の部分構造検索を行ってみます。
sdfからの分子の入力方法は、6/8の記事と同じ方法で行っています。
今回のプログラムのUsageは以下のとおりです。

　Usage: smartsearcher.py xxx.sdf 'PATTERN' > result.txt

まず、'PATTERN'として入力するSMARTS文字列をcompile関数に渡します。

　pattern = Smarts.compile(sys.argv[2])

compile関数は戻り値として、Matcherクラスのインスタンスを返します。

次にwhile文の中で順次読み込んでいるsdf中の分子とSMARTS文字列の一致判定を行います。この判定には、matchメソッドを利用し、第2引数に1を渡しています。
例えば、分子中にSMART文字列と一致する部分が5箇所あったとしても最大1つだけ見つければよいことを示しています。全ての一致が必要な場合は、-1を渡します。デフォルトが-1ですので、第2引数を渡さなければ全ての一致を検出します。

　match = pattern.match(mol,1)

matchメソッドは一致がない場合Noneオブジェクトを返しますので、以下の式で一致を判定します。

　if match != None:

また、逆に一致する場合、matchメソッドは、SMARTS文字列に一致したatomとbondのPathSetを返します。そして、次のような結果を表示します。

--------------------------------------------------
C1(=C2C(=C(C(=O)C=C2)Br)Oc3c(c(ccc13)O)Br)c4c(C(=O)O)cccc4
match 1
atoms (Atom(18), Atom(17), Atom(16), Atom(15), Atom(20), Atom(19), Atom(22))
bond (Bond(19), Bond(18), Bond(16), Bond(22), Bond(21), Bond(20), Bond(24))
--------------------------------------------------
c1(c(cc(c(c1)O)C(C)(C)C)O)C(C)(C)C
match 1
atoms (Atom(4), Atom(14), Atom(13), Atom(8), Atom(7), Atom(5), Atom(6))
bond (Bond(14), Bond(13), Bond(12), Bond(7), Bond(6), Bond(4), Bond(5))
--------------------------------------------------



Frownsでは、SMARTS検索の一部にgraph matching libraryであるVFLibを利用しています。興味のある方はチェックしてみてください。



人気ブログランキング(クリックして応援してね)

今回は、Frownsを用いたsdfの読み込みについてメモしたいと思います。
FrownsのHPのTutorialに例がありますので、これとは少し異なる方法をとってみます。
ちなみにTutorial中の例は古いバージョンのFrownsを使っている?と思われるため、Ver.0.9aを使う場合は、トップページ上にあるMOL File(example4.py)を利用した方がいいと思います。下記プログラムもVer.0.9aを用いています。

まず、組み込み関数open()に引数としてsdfのファイル名を渡します。そして、その戻り値をMDLモジュール内のsdinクラスの__init__メソッドに渡すことにより、そのインスタンスmdl_infoを生成させます。

次にwhile文を使ってsdfから順次分子を読み込み、それら分子のCanonical Smiles表記を表示します。
ここで、mol_info.next()は3つの戻り値をタプルで返します。
molにはMoleculeクラスのインスタンス、errにはエラーメッセージ、textには、sdfの中身がtextとして収められています。

while文からの脱出は、StopIterationにより行っています。






人気ブログランキング(クリックして応援してね)

ケモメトリックス(化学計量学)を研究現場で実践されている研究者は多いのではないでしょうか。ケモメトリックスは、化学諸問題に対する回帰的学習による問題解決の基盤となる学問・技術であると思います。その背景には、多変量解析、統計解析があるため、ケモメトリックスを学習するには、計量学の理解が必要不可欠となります。

「化学者のための多変量解析」-ケモメトリックス入門-は、統計学の基礎から、回帰分析、主成分分析、PLS回帰分析、そして最後に実践的な適応事例と学部生から研究者まで活用できる書籍だと思います。

回帰分析や主成分分析は、多くの統計関連の書籍で解説されていますが、PLS回帰分析をここまで丁寧かつ分かりやすく記述している書籍は少ないと思います。


PLS回帰分析は、3D-QSARにおいても重要な技術ですので、下記のような専門書を読む前に理解しておくと必ず役に立つと思います。


先日、少し取り上げましたが、PyChemやPyRはケモメトリックスの実践に役に立つライブラリだと思いますので、今後メモしていこうと思っています。




人気ブログランキング(クリックして応援してね)

6月4日の記事で、シクロヘキサンのねじれ角(C1-C2-C3-C4)の計算にOELibで作成したtorsionを利用しました。
そこで、今回は、OELibを用いた分子の入力からねじれ角の計算までをtorsion.cppを通して紹介したいと思います。

torsion.cppのソースコードを下のインラインフレーム内に記載しております。まずは、OELibのヘッダファイルであるopenbabel/mol.hとopenbabel/obconversion.hをincludeし、namespaceを定義します。ちなみに、torsionのUsageは以下のとおりです。
Usage: torsion filename atom1 atom2 atom3 atom4
ここで、filenameは分子ファイル名を示し、atomXは、ねじれ角を構成する原子のIDを示します。
次に引数チェックのあと、分子ファイルの拡張子からファイル形式を判定しています。分子の読み込みは、ファイルから1分子づつfor文の無限ループ内で行っており、mol.Empty()がtrueとなればbreakします。プログラムのメインとなるねじれ角の計算はOBMolクラスのGetTorsionを用いて計算しています。
今回定義したシクロヘキサンのねじれ角は下図のようになります。


torsion.cpp


コンパイル：
$ g++ -I/usr/local/include/openbabel-2.0 -L/usr/local/lib -o torsion　torsion.cpp -llibopenbabel

実行：
$ torsion cyclo.sdf 1 2 3 4

出力：
conf-1 torsion_angle 304.814
conf-2 torsion_angle 291.955
..
..
conf-100 torsion_angle 329.156

無事にねじれ角が計算できると思います。




人気ブログランキング(クリックして応援してね)

先日、Javaのケムインフォマティクス用ライブラリであるCDKを紹介しましたが、もう一つ有名なライブラリとして、JOELibがあります。JOELibは、OpenBabelのライブラリであるOELibのJava版です。とはいえ、OELibを単にJavaに置き換えたというだけでなく、JOELibにしかないクラスもあり、より拡張されたライブラリといえます。

個人的に、CDKは、多機能ではあるが、ややクラスの使い勝手が悪いところ及び一部に実装がまだ十分でないのでは？(新しい機能を積極的に取り入れているからでしょう)というところがあり、一方、JOELibはCDK程の多機能さはないが、クラスの使い勝手や実装も吟味されているという印象をもっています。

JOELibのおもしろそうなところは、Matlab interfaceを持っているところです。Matlabといえば、その汎用シミュレーションパッケージであるSimulinkも便利であり、制御系では必須のツールだと思います。さらに、最近では、システムバイオロジーの分野でも広く利用されていると思いますので、Matlab interfaceに関しては、ちょくちょく調べようと思っています。

ただ、どちらもケムインフォマティクスの基礎となるライブラリは十分に実装されていますし、JOELib中のJOEMolとCDK中のMoleculeクラスは変換可能ですので、2つのライブラリを同時に利用しながら、アプリケーションの開発もできます。

インストー方法も含め今後メモを積み重ねたいと思っています。




人気ブログランキング(クリックして応援してね)

FFEからのTinkerの操作がなかなか便利でしたので、FFEを使った分子の配座探索についてメモすることにしました。
配座解析には様々な方法が提案されていますが、今回は、Tinkerに含まれるDynamicを用いた分子動力学法(MD)により行いたいと思います。

例として、シクロヘキサンの配座解析を行ってみます。
1. 適当なソフトウェアでシクロヘキサンの三次元座標のmolファイルを準備する(cyclo.mol)。そして、babelを使ってtinker入力ファイルを作成する。
$ babel -imol cyclo.mol -otxyz cyclo.xyz

2. FFEを起動し、cyclo.xyzを読み込む。力場はmm2。


3. [Modeling Commands]タブをクリックし、[Optimize]をリストから選択し、[ロケット]ボタンを押す。
これにより、初期配座の最適化を行います。

4. 次に、[Dynamic]をリストから選択し、[Number of Steps]に100000を設定、[Temperature]に1000を設定、[Dump]に1を設定(あとはデフォルト値)する。そして、[ロケット]ボタンを押す。
ここでは、1000Kで100psのMDを行い、配座のサンプリングを1psごとに実施しています。

5. [Archive]をリストから選択し、[ロケット]ボタンを押す。
これで、cyclo.arcファイルが生成されます。これをFFEで開くとアニメーション表示で見ることができます。さらに、[Format]タブを選択し、ラジオボタン[X]を選択し、[ロケット]ボタンを押す。これにより、100個の配座が含まれるcyclo.xmolが生成されます。

6. cyclo.xmolをbabelによりsdf形式に変換する。
$ babel -ixyz cyclo.xmol -osdf cyclo.sdf
シクロヘキサンの配座を解析するために、ねじれ角(C1-C2-C3-C4)を計算します。ここでは、OELibを使った自作プログラム"torsion"を使って計算します。
$ torsion cyclo.sdf 1 2 3 4 > result.out
(torsionの詳細は、近々紹介しますね)

結果を下図に示します。横軸にシミュレーション時間、縦軸にねじれ角を示しています。複数のパターンがサンプリングされています。図右に代表的な２つのパターンを示します。


初期配座、シミュレーション時間、アンサンブルの検討など全くやってませんので、結果をまともに解釈はできませんが、FFE(Tinker)の使い方のメモということでご了承ください。
今回は、WindowsでFFEを使ってみましたが、インストールは超簡単ですし、本当に便利だなと思いました。

MDによるシクロヘキサンの配座解析は、下の書籍の中で詳しく紹介されていますので、興味のある方は読んでみてください。
「分子モデリング―基本原理と創薬への応用」江崎俊之訳 地人書館




人気ブログランキング(クリックして応援してね)

タンパク質間相互作用(PPI)は、転写、翻訳、複製、シグナル伝達など様々な細胞内の重要なプロセスを担っています。したがって、PPIに係わる膨大な研究成果をデータベース化することは、これまでに分かっている現象を単純に知るというだけではなく、新しい解釈、系統的な整理、知識化を導く可能性を秘めています。さらに創薬においても、標的タンパク質の選別、副作用の検討など創薬研究者に対する重要な情報源になっていると思います。

現在、多数のPPI Databaseが公開されています(以下のサイトをご参考に)。
http://www.blueprint.org/bind/bind_relateddatabases.html

今回は、PPIで有名なデータベースの１つであるBINDについてメモしたいと思います。私はBINDというと、どうしてもDNSサーバが真っ先に浮かんできてしまいます(笑)。

The Biomolecular Interaction Network Database (BIND)は現在、Commercial版とOpen Access版があります。
http://www.unleashedinformatics.com/index.php?pg=products&refer=bind
ここでは、Open Access版を利用したいと思います。

今回は、p53のPPIについて検索してみることにします。

1. Open AccessのBINDをクリックし、はじめて利用する場合は、ユーザ登録を、そうでない場合はログインする。

2. メイン画面上のSearch Boxにp53と入力し、GOボタンを押す。

3. p53のPPIリストが表示されるので、Mdm2とのPPIを示す"interaction 301661"をクリックする。
"interaction 301661"が表示されていなければ、Search Boxに"301661"と入力し、検索してもOKです。ちなみに、この数値はBIND Idです。


4. BIND id: 301661に関する情報が表示されるので、右側にある"Visualize using..."というリストボックスから"Interaction Network 3.5"を選択する。
すると、ソフトウェアのダウンロード(Java Web Start)がはじまる。

5. ダウンロードが終わり、ソフトウェアが立ち上がると以下の画面が表示される。

p53とMdm2が結合していますね。
BINDのおもしろいところは、"Function","Localization","Binding"に関連する情報を“絵”として表現しているところです。この“絵”の意味は右側に説明がありますので、瞬時にそのタンパク質の概要を知ることができます。

6. 最後に、Mdm2をダブルクリックする。

おお！Mdm2と相互作用するタンパク質がでてきました。
この操作を繰り返すことにより、PPIネットワークをサーフィンできるわけです。

この他にも様々な使い方があるようですので、興味のある方はお試しください。




人気ブログランキング(クリックして応援してね)