Pharmaceutical Developer's Note: 2010-02

2010年2月27日土曜日

PythonでOpenBabelを使う～Pybelを使う

PybelはOpenBabelの様々なクラス・メソッドの中でよく使う物をPython用に定義したお便利環境です。

詳細は
http://openbabel.org/wiki/Using_OpenBabel_from_Python#Pybel
http://openbabel.org/pybel.html

例えばSDファイルを読み込んでMolに変換して構造を表示するのは簡単。

from pybel import *

inSdFilePath = "test.sdf"
mols = readfile("sdf", inSdFilePath)
for mol in mols:
mol.make3D()
mol.draw() #OASAがインストールされていれば構造を表示する

PythonでOpenBabelを使う～SDファイル

とりあえず基本的なやつ。

openbabeltest.py

# coding: shift-jis

from openbabel import *

def ReadSdFile(inSdFilePath="test.sdf"):

#inSdFilePath = "test.sdf"

conv = OBConversion()

answer = conv.SetInAndOutFormats("sdf","can")

mol = OBMol()

end_flag = conv.ReadFile(mol,inSdFilePath)

smiles = []

while end_flag:

smile = conv.WriteString(mol).split('\t')[0] # Smiles文字のみをSplitする

smiles.append(smile)

answer = mol.Clear()

end_flag = conv.Read(mol)

return smiles

# main

smilesStrigList = ReadSdFile()

for smiles in smilesStrigList:

print smiles

コマンドプロンプト

I:\Python>python openbabeltest.py

ClCC(=O)NC1=C(SCC1)C(=O)OC

COC(=O)C1=C(CCS1)NC(=O)c1ccccc1

COC(=O)C1=C(CCS1)NC(=O)c1cccc(C)c1

COC(=O)C1=C(CCS1)NC(=O)c1ccccc1Cl

COC(=O)C1=C(CCS1)NC(=O)c1ccc(C)c(C)c1

PythonでOpenBabelを使う

TIBCO Spotfire3.1が欧米でリリースされ、日本国内でも近々リリースされると思います。

今回のリリースでの目玉はTextView上でのスクリプトの実行が可能になり、ガイドっぽいことができるようになるのとHeatMapのサポート、Bar&LineChartなどいろいろありますが（詳細）、地味ですがIronPythonのScriptのサポートがあります。

これを機にちょっとPythonを勉強しておこうと思って、環境を作ります。詳細はこちら。

Python2.6をインストール　http://www.python.org/download/

最新版は3.1なのですがOpenBabelで構造表示に使うPILとOASAが2.6までしか対応していないので。

Pythonのパスを設定
OpenBabel GUIをインストール
OpenBabelのパスを通す

PATHにI:\Program Files\OpenBabel-2.2.3を追加

環境変数「BABEL_DATADIR」をセット

例えば、set BABEL_DATADIR=I:\Program Files\OpenBabel-2.2.3な感じ

OpenBabel Python bindingsをインストール

PILとOASAをインストール、いずれもPython2.6用を入れます
OpenBabeｌのWebにあるテストスクリプトを動かしてみる

テストスクリプトのうち「mol.make3D()」を実行してみて、

==============================

*** Open Babel Error in OpenBabel::OBForceFieldMMFF94::ParseParamFile

Cannot open parameter file

こんなエラーが出る場合は、環境変数のBABEL_DATADIRが正しくセットされていない時です。コマンドプロンプトでecho %BABEL_DATADIR%とやってみてOpenBabelのGUIのインストールフォルダーが出ることを確認しましょう。

2010年2月21日日曜日

Haskell入門その4

再帰処理
---------------------------------------
-- Listの長さを取得する
---------------------------------------
mylength :: [a] -> Integer -- 型宣言
mylength [] = 0 -- リストの最後の定義
mylength (x:xs) = 1 + mylength xs -- 先頭以外のリストを再帰計算
{-
  x=最初の要素、xs=残りのリスト
  head :: [a] -> a
  head (x:xs) = x
  tail :: [a] -> [a]
  tail (x:xs) = xs

  使い方
  I:\haskell>ghci
  Prelude> :l test012.hs
  *Main> mylength [1..999]
  999
-}

ちょっとわかってきたような気になってきたことが示唆される

2010年2月20日土曜日

Haskell入門その3

module Main

where

import IO

import Random

main = do

hSetBuffering stdin LineBuffering

num <- randomRIO (1::Int, 100)

putStrLn "I'm thinking of a number between 1 and 100"

doGuessing num

doGuessing num = do

putStrLn "Enter your guess:"

guess <- getLine

let guessNum = read guess

if guessNum < num

then do putStrLn "Too low!"

doGuessing num

else if read guess > num

then do putStrLn "Too high!"

doGuessing num

else do putStrLn "You Win!"

これをコンパイルしてみた

I:\haskell>ghc --make test011.hs -o test011.exe

なんじゃ、このバカでかいファイルは！？

Haskell入門その2

回帰

{----------------------------
  N!
-----------------------------}
my_factorial 1 = 1
my_factorial n = n * my_factorial (n-1)
{----------------------------
  べき乗
  a^b
-----------------------------}
my_exponent a 1 = a
my_exponent a b = a * my_exponent a (b-1)
{----------------------------
  フィナボッチ数列
  n番目のフィナボッチ数を求める
-----------------------------}
fina 1 = 1
fina 2 = 1
fina n = fina (n-2) - fina (n-1)
{----------------------------
  a : b
-----------------------------}
mult a 0 = 0
mult a b = a + mult a (b-1)

GHCIからhsファイルを読み込む

I:\haskell>ghci
Prelude> :l test008.hs
*Main>

Haskell入門その1

test001.hs

module BMI where
type BMI = Double
stdBMI :: BMI
stdBMI = 22.0
type Height = Double
type Weight = Double
bmi :: (Height, Weight) -> BMI
bmi (h,w) = w / (h ^ 2)

コマンドプロンプト

I:\haskell>ghci test002.hs
GHCi, version 6.12.1: http://www.haskell.org/ghc/ :? for help
Loading package ghc-prim ... linking ... done.
Loading package integer-gmp ... linking ... done.
Loading package base ... linking ... done.
Loading package ffi-1.0 ... linking ... done.
[1 of 1] Compiling BMI ( test002.hs, interpreted )
Ok, modules loaded: BMI.
*BMI> bmi(1.71,79.5)
27.18785267261722

う～む、まだ楽しくない

2010年2月13日土曜日

IUPAC名からSmilesを起こす

ということで苦労しまくったIUPAC名からSmilesを起こすプロジェクトですが、最終的には以下のような流れになりました。

１）Opsin0.5.3でIUPAC名をCMLに変換
２）CDKでCMLをIMoleculeに変換
３）CDKでIMoleculeからMolファイル形式に変換
４）OpenBabelでMolファイルからSmilesに変換

フォームアプリケーション的はコードを書きます。Smilesにするだけなら2次元構造は必要なないのですが、忘記録的に入れときます。この流れがdelegateでも動いたらOkですな。

using System;
using System.Windows.Forms;
using System.IO;
using System.Diagnostics;

using org.openscience.cdk.io;
using org.openscience.cdk;
using org.openscience.cdk.interfaces;
using org.openscience.cdk.layout;
using java.io;

using uk.ac.cam.ch.wwmm.opsin;
using OpenBabel;

namespace OpsinTest
{
   public partial class Form1 : Form
   {
   NameToStructure nts = null;

   public Form1()
   {
   InitializeComponent();
   nts = new NameToStructure();
   }

   private void button3_Click(object sender, EventArgs e)
   {
   string IUPAC = textBox1.Text;

   try
   {
   textBox2.Text = "";
   Application.DoEvents();

   // Create CML from IUPAC name (Opsin)
   string cml = nts.parseToCML(IUPAC).toXML();

   // Convert CML to IMolecule (CDK)
   StringBufferInputStream str_stream = new StringBufferInputStream(cml);

   CMLReader cmlr = new CMLReader();
   cmlr.setReader(str_stream);

   ChemFile chemFile = new ChemFile();
   ChemFile chem = (ChemFile)cmlr.read(chemFile);

   IMolecule mol = chem.getChemSequence(0).getChemModel(0).getMoleculeSet().getMolecule(0);

   // Convert from IMolecule to SD (CDK)
   java.io.StringWriter sww = new java.io.StringWriter();
   MDLWriter mw = new MDLWriter(sww);
   StructureDiagramGenerator sdg = new StructureDiagramGenerator();

   sdg.setMolecule(mol);
   sdg.generateCoordinates();
   mol = sdg.getMolecule();
   mw.write(mol);
   string sd = sww.toString();

   // Convert from SD to smiles (OpenBabel)
   OBMol obMol = new OBMol();
   OBConversion obConv = new OBConversion();
   if (!obConv.SetInAndOutFormats("mol", "can")) { return; }
   if (!obConv.ReadString(obMol, sd)) { return; }

   textBox2.Text = obConv.WriteString(obMol);

   obMol.Dispose();
   obConv.Dispose();
   obMol = null;
   obConv = null;
   }
   catch (Exception ex)
   {
   string mes = ex.Message;
   textBox2.Text = mes;

   }
   }

   }
}

OBConversionでメモリーリーク？

OpenBabelの構造コンバーターの「OBConversion」はめちゃくちゃたくさんの構造フォーマットに対応した便利なConverterなんですが、どうも変な動きをすることがあります。

たとば、下記のようにMolの配列にSmiles文字列の配列をConvertしながら入れる場合、下記のコードはVisualStudioでは特に問題ありません。

（コードA）
for(i=0; i < 3; i++)
{
   OBConversion conv = new OBConversion();
   conv.SetInAndOutFormats("smi", "can");
   conv.ReadString(obMol[i], obSmi[i]);
}

しかし下記のコードのように、CML文字列を変換する場合、実行中は何も問題ないのですが、終了するとメモリーエラーを起こします。

（コードB）

for(i=0; i < 3; i++)

{

OBConversion conv = new OBConversion();

conv.SetInAndOutFormats("cml", "can");

conv.ReadString(obMol[i], obCml[i]);

}

これは正しくは、こうするべきです。

（コードC）

for(i=0; i < 3; i++)

{

OBConversion conv = new OBConversion();

conv.SetInAndOutFormats("cml", "can");

conv.ReadString(obMol[i], obCml[i]);

conv.Dispose();

conv = null;

}

つまり、使用後のOBConversionインスタンスはちゃんと破棄する必要があります。このコードでは終了しても何もエラーがなく正常に動きます。

ただ不思議なのはコードAではメモリーエラーを起こさず、コードBではエラーを起こす点です。

想像するにCMLのConvertで何か内部的なメモリーの破棄忘れがあって、破棄しない場合にメモリーリークを起こしているではと思ってます。

さらにはコードCでもメモリーエラーを起こす場合があります。これはコードCを含む関数をdelegateで使う場合で、OBConversionのインスタンスが正常に破棄されず、ReadStringの箇所でメモリーエラーを起こします。

ということで、実用問題としてOpenBabelのCML変換は使えない、という結論に達しました。

これがわかるのに2日もかかりました。

ということでCML変換はCDKでやることにしました。その内容はまた別に書きます。

OpsinをC#で使う

以前に.NETでCDKを使うで書いたように、JAVAの実行ファイル（.jarファイル）を.NETのclass library（dllファイル）にConvertする方法を使えば、いろいろなJavaの便利環境を.NETで使うことができます。

今回は私の心の師匠のkさんのhttp://blog.kzfmix.com/entry/1214827204で使っていたOpsinを.NETで使えるか試してみました。

Opsinは化合物のIUPAC名をCMLに変換するモジュールで、CMLからさらにSDやSmilesなどに変換することができます。これを使えば、論文や特許に記載されている化合物の文字列から構造を起こすことができるという優れものです。

使ったOpsinはOscar3に入っているOpsin-0.5.3.jar。実は最初に0.1.0を試したのですが、どうもIUPACの認識率が低く、使い物にならんなーと思って諦めかけたのですが、このバージョンになってだいぶ良くなったので使うことにしました。

いつもどおりにIKVMでjarをdllに変換します。このときkeyfileがあれば指定することもできます。keyfile付きで変換することで「厳格な名前」付きdllにすることができるというわけです。

出来上がったdllをVisualStudioのC#のプロジェクトにいれます。そして

using OpenBabel;
using uk.ac.cam.ch.wwmm.opsin;

NameToStructure nts = new NameToStructure();
string cml = nts.parseToCML("4-iodobenzoic acid").toXML();

OBMol obMol = new OBMol();
OBConversion conv = new OBConversion();
conv.SetInAndOutFormats("cml", "can");
conv.ReadString(obMol, cml);
string smiles = conv.WriteString(obMol);
string[] smile = smiles.Split(new string[] { "\t" }, StringSplitOptions.None);
string result = smile[0];

こんな流れのコードを使えばIUPAC名からSmiles文字列を作ることができます。

ただ、特許や論文のIUPAC名はけっこういい加減で、" ' -　などの使い方、あるいは半角スペースの入り方でOpsinがエラーを吐くことが多いです。Opsinに渡す前にIUPAC名のCleanが必要なようです。

またここで書いたOBConversionクラスですが、どうもメモリーリークがありそうなので、CMLの変換には使わない方が良いでしょう。このあたりについては次に・・・・

2010年2月4日木曜日

Spotfire DecisionSiteで.NET

ずっとTIBCO Spotfire+C#で仕事していると、たまにSpotfire DecisionSiteのスクリプトをいじらなきゃいけなくなったときに超めんどくさくなります。

なんといってもデバッグできない。昔はしこしこJavascriptでAlertとかログファイルに書き出したりしてなんとか構築したのに、もう戻れなくなっている。

ということで、DecisionSiteでVisualStudio C#を使う方法を検討しました。もちろん直接Javascriptのデバッグもできるとは思いますが、それじゃ面白くないので、C#でクラスライブラリーを作ってそれを参照させることにしました。

VisualStudio C#でクラスライブラリーを作ります。このとき「COM参照」ができるようにPropertiesを設定しておきます。

ClassはPublicで作って普通にコンパイルするとdllとtblファイルが出来ます。この2つのファイルをDecisionSiteのガイドやツールに登録してやればJavascriptのnew ActiveX(<namespace.class名>)で使えるようになります。あとは登録したpublic methodを呼べば使えるというわけです。

思っていたより結構気軽なので、よく使うレコードセット処理や文字列処理を入れて使う予定です。

2010年2月2日火曜日

TIBCO SpotfireとOracle Procedure

TIBCO SpotfireではOracleのテーブル、ビューに加えてProcedureに対応しています。Procedureというより、REF CURSORで返されるデータセットを受け取れるといったほうが正確ですね。

具体的には、Packageとして作成したFunctionでREF CURSORとして返すような動的なSQLを作ることで、柔軟なOracleクエリーの結果を受け取ることができます。

これが意外と便利で、作り方によっては同じクエリーでもユーザごとに異なる結果をとったり、テーブルをまたいだ結果をとったり、ユーザ操作をハックして選択したデータに応じた結果を取得できます。

もちろん通常の操作で、データの更新や追加、消去などをProcedureに入れることで、SpotfireからOracleのデータテーブルの操作もできます。

一度トライしてみてください。

登録: 投稿 (Atom)