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

隣接行列と結合行列は、計算機上での分子の表現の基礎ともいえる重要な表現方法です。今回はCDKを用いてこれら行列の取得の仕方をメモしたいと思います。

CDKではorg.openscience.cdk.graph.matrixパッケージに含まれるAdjacencyMatrixクラスとConnectionMatrixクラスを用いて取得できます。

　System.out.println("AdjacencyMatrix");
　int am[][] = AdjacencyMatrix.getMatrix(mol);
　for(int i=0;i＜am.length;i++){
　　for(int j=0;j＜am[i].length;j++){
　　　System.out.print(am[i][j]+" ");
　　}
　　System.out.println();
　}

　System.out.println("\nConnectionMatrix");
　double cm[][] = ConnectionMatrix.getMatrix(mol);
　for(int i=0;i＜cm.length;i++){
　　for(int j=0;j＜cm[i].length;j++){
　　　System.out.print(cm[i][j]+" ");
　}
　　System.out.println();
　}

例として、以下の構造を入力してみます。

結果：
AdjacencyMatrix
0 0 0 1 1 1
0 0 0 0 0 1
0 0 0 0 0 1
1 0 0 0 0 0
1 0 0 0 0 0
1 1 1 0 0 0

ConnectionMatrix
0.0 0.0 0.0 1.0 2.0 1.0
0.0 0.0 0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0 0.0 1.0
1.0 0.0 0.0 0.0 0.0 0.0
2.0 0.0 0.0 0.0 0.0 0.0
1.0 1.0 1.0 0.0 0.0 0.0

隣接行列は、原子同士が結合していれば1を、そうでなければ0を当てはめます。一方、結合行列は、原子同士が結合していれば結合次数を、そうでなければ0を当てはめます。結合行列がdouble型の要素をもつ理由は、aromatic bondなどを1.5と表記したいためです。



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

JOELibのAbstractDatabaseクラスを用いたHashcodeの生成についてメモしたいと思います。このクラスにはHashcodeを生成するために以下の3つのメソッドが用意されています。

　getHashcode(JOEMol mol)
　getSMILESHashcode(JOEMol mol)
　getMoleculeHASH(JOEMol mol)

ただし、最後のメソッドは上2つのメソッドの結果を同時に返すメソッドですので実質的には2種類のHashcodeを得ることができます。
それでは、getHashcodeとgetSMILESHashcodeの違いを具体的に見てみます。入力した分子構造は、L体、D体のアラニンの部分構造です。

　JOEMol mol1 = new JOEMol();
　JOEMol mol2 = new JOEMol();
　String smiles1 = "CC(=O)[C@H](C)N";
　String smiles2 = "CC(=O)[C@@H](C)N";
　JOESmilesParser.smiToMol(mol1, smiles1, "mol1");
　JOESmilesParser.smiToMol(mol2, smiles2, "mol2");
　
　int hashcode1 = AbstractDatabase.getHashcode(mol1);
　int hashcode2 = AbstractDatabase.getHashcode(mol2);

　int smileshash1 = AbstractDatabase.getSMILESHashcode(mol1);
　int smileshash2 = AbstractDatabase.getSMILESHashcode(mol2);

　System.out.println("getHashcode()");
　System.out.println(hashcode1);
　System.out.println(hashcode2);

　System.out.println("getSMILESHashcode()");
　System.out.println(smileshash1);
　System.out.println(smileshash2);

結果：
getHashcode()
-715427484
-715427484

getSMILESHashcode()
1803284766
1004112044

getHashcode()ではD体、L体の区別ができないないのに対し、getSMILESHashcode()では区別されています。大きな違いはここにあるのですが、どちらのメソッドも入力する分子の番号付けが変わると異なるHashcodeが生成されてします。したがって、Morgan法などで再番号付けされた分子構造を入力しなければ特定の分子に対してユニークなHashcodeを得ることができませんので注意が必要です。



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

CDKからJmolを利用する方法についてメモしたいと思います。最も簡単な方法は、CDK News, 2.1(2005)で記載されているCdkJmol3DPanelクラスの利用だと思います。ただし、私の利用したCDK及びJmolのバージョン(cdk-20050826,Jmol 10.00)では、CdkJmol3DPanelクラスは動作しませんでした。ネットで調べると、CdkJmol3DPanelクラスとほぼ同等のJmolPanelクラスが公開されており、こちらでは問題なく動作しました。ソースを見るとJmolViewerクラスの使い方が異なっていることが分かりました。また、CDKとJmolの連携に重要なCdkJmolAdapterクラスの使い方もJmolPanelクラスから知ることができます。

JmolPanelクラスさえ使えればあとは簡単です。

　JmolPanel jp = new JmolPanel();
　jp.setMol(mol);
　JFrame frame = new JFrame();
　Container cont = frame.getContentPane();
　cont.add(jp);
　frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
　frame.pack();
　frame.setBounds(10,10,300,300);
　frame.setVisible(true);

こんな感じでJmolが利用できます。




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

PythonからOpenBabelのクラスを利用して、簡易版obgrepを作成したいと思います。

simple_obgrep.pyはsdf形式のファイルを入力し、部分構造としてベンゼン(c1ccccc1)を含んでいる分子を検出するシンプルなプログラムです。

実行：
　$ python simple_obgrep.py > output.sdf

simple_obgrep.py:




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

OpenBabelでは、3種類のfingerprintを利用してsimilarity searchを行うことができます。fingerprintの比較はTanimoto coefficientで行われています。

それでは、実際にsimilarity searchをしてみます。

　$ babel 'query.smi' 'mydata.sdf' -ofpt -xfFP2

　結果：
　> Tanimoto from first mol = 0.166667
　> Tanimoto from first mol = 0.285714
　> Tanimoto from first mol = 0.27027
　・
　・
　> Tanimoto from first mol = 0.0802469
　> Tanimoto from first mol = 0.149425
　> Tanimoto from first mol = 0.122449
　39 molecules converted

ここでは、query.smiとmydata.sdfに含まれる全ての分子との間のTanimoto coefficientが計算されています。-xfFP2はfingerprintにFP2を利用することを示しています。fingerprintの種類は以下のコマンドで知ることができます。

　> babel -F
　FP2 -- Indexes linear fragments up to 7 atoms.
　FP3 -- SMARTS patterns specified in the file patterns.txt
　FP4 -- SMARTS patterns specified in the file SMARTS_InteLigand.txt

patterns.txtとSMARTS_InteLigand.txtは、通常、/usr/local/share/openbabel/2.0.2などにインストールされていると思います。このファイルに独自にpatternを追加すればオリジナルのfingerprintで探索ができます。



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

Lazy Structure-Activity Relationships(lazar)という毒性予測ツールが公開されています。JOELibが利用されているとのことです。Web版も公開されているため試しに使ってみるにも便利だと思います。Rat,Mouse,Hamsterに対するCarcinogenicityなどの予測ができるようです。論文のpreprintもPDFで公開されていますので、内部的な処理も理解でき、勉強になると思います。


話は変わりますが、今日、“旅の指さし会話帳DS ドイツ”を買いました。少し遊んでみましたがなかなか面白いソフトです。ただ、私には実際に旅先でDSを相手に向けて質問することはできそうにありません(笑)。





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

OpenBabelではSWIGを用いて、そのほとんどクラスをPythonから利用できるようにしています(Perlからも利用できるみたいです)。

OpenBabelを普通にインストールしただけではPythonから利用できないので、ちょっとしたインストール作業が必要です。

1. openbabel-2.0.1.tar.gzを展開したディレクトリに移動します。

　＞cd $OpenBabelSrc/openbabel-2.0.2/scripts/python

2. ビルド、インストール作業を行います。

　＞python setup.py build
　＞python setup.py install

これで、PythonからOpenBabelのクラスが利用できるようになります。

3. 簡単なサンプルプログラムを動かしてみます。

　sample.py：
　import openbabel

　mol = openbabel.OBMol()
　conv = openbabel.OBConversion()
　conv.SetInAndOutFormats("smi","mdl")
　conv.ReadString(mol,"c1ccccc1O")
　print mol.NumAtoms()
　mol.AddHydrogens()
　print mol.NumAtoms()
　conv.WriteFile(mol,'phenol.mol')

実行:
　$python sample.py

出力：
　7　<-水素なしの時の原子数
　13　<-水素ありの時の原子数

出力ファイル(phenol.mol):

OpenBabel

13 13 0 0 0 0 0 0 0 0999 V2000
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 O 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
1 6 2 0 0 0
1 2 1 0 0 0
1 8 1 0 0 0
2 3 2 0 0 0
2 9 1 0 0 0
3 4 1 0 0 0
3 10 1 0 0 0
4 5 2 0 0 0
4 11 1 0 0 0
5 6 1 0 0 0
5 12 1 0 0 0
6 7 1 0 0 0
7 13 1 0 0 0
M END

OpenBabelでは、まだ2次元座標の生成機能がないこともあり、xyzが全て0となります。
ScriptからOpenBabelを利用できるとコンパイル作業もなく便利ですね。



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

PymolのpluginにGlobal Protein Surface Survey(GPSS)というpluginが公開されていますが、かなり便利だと思います。タンパク質のpocketやcavityの提示だけではなく、ligand binding surfaceやDNA binding surfaceなども提示してくれます。実際の計算はpluginで行っているのではなく、CASTpやCSAのデータを利用しているようです。

インストールは簡単で、GPSSのサイトからGPSSpyMOL.pyを取得して、$PYMOLPATH/module/pmg_tk/startup/にコピーするだけです。

具体的には、以下の様なグラフィクスが描画され、pocketの形状がよく分かります。





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

CDKでは、Fingerprinterクラスを用いてfingerprintの生成を行います。先日メモしたJOELibでのfingerprint生成法とは異なり各ビットの特性を事前に定義する必要はありません。

CDKのfingerprint生成の概要ですが、まず、ある原子をrootとし、そこから任意の深さ(searchDepth)まで、DFS(Depth-First Search)を行い、経由した経路をSMILESなどの線形表記で文字列に置き換えます。次に、任意の最大値(size)を設定したハッシュ関数により文字列をハッシュ値に変換します。最後に、このハッシュ値に位置するビットを1にセットすることによりfingerprintを生成します。sizeとsearchDepthを適切に設定しないとハッシュ値の衝突や異なる分子で同じfingerprintが多発するという現象が生じます。

Fingerprinterクラスのデフォルトでは、size=1024, searchDepth=6となっています。JOELibで生成した葉酸のfingerprintと比較するためにsize=54, searchDepth=3としてfingerprintを生成したいと思います。

　mol = (Molecule)mr.read(new Molecule());
　BitSet finger = Fingerprinter.getFingerprint(mol,54,3);
　for(int i=0;i＜54;i++){
　　if(finger.get(i)){
　　　System.out.print("1");
　}else{
　　　System.out.print("0");
　}


　　　　　　　　　　　　　　↓
Fingerprint:
111001011110111110111111111110111111101111101101111111



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

JOELibを使って分子のfingerprintを作成したいと思います。fingerprintは以前pythonのcheminformaticsライブラリであるfrownsの記事でも取り上げました。JOELibでは、SSKey3DSクラスを用いてfingerprintを生成できます。生成されるfingerprintは54bitです。各bitが何を示すのかはココをご覧ください。

まずは、SSKey3DSクラスのgetDescInfo()メソッドを用いて、Descriptorの情報を調べます。

　SSKey3DS ssk = new SSKey3DS();
　DescriptorInfo di = ssk.getDescInfo();
　System.out.println(di.getName());
　System.out.println(di.getResult());

出力：
　Pharmacophore_fingerprint_1
　joelib.desc.result.BitResult

SSKey3DSはBitResultを返すことが分かりました。

それでは、実際にfingerprintの取得を行ってみます。fingerprintの計算は、calculate(JOEMol mol)メソッドで行います。戻り値であるDescResultはBitResultにキャストします。そして、maxBitSizeでfingerprintの長さ及びgetBinaryValue()メソッドで生成されたfingerprintを取得します。

　SSKey3DS ssk = new SSKey3DS();
　sreader.readNext(mol);
　DescResult result = ssk.calculate(mol);
　int fp_size = ((BitResult)result).maxBitSize;
　JOEBitVec br = ((BitResult)result).getBinaryValue();
　for(int i=0;i＜fp_size;i++){
　　if(br.get(i)){
　　　System.out.print("1");
　　}else{
　　　System.out.print("0");
　　}
　}

最後に、例として葉酸のfingerprintを生成したいと思います。


　　　　　　　　　　　　　　　↓
Fingerprint:
111001110000101100100010000100101010001011111111111111

以前メモしたScreeningAssistantのdiversity計算には、このSSKey3DSクラスで作成したfingerprintが利用されているようです。



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