2006年06月30日

mol2ps:rotate機能

mol2psの便利なオプションとして"--rotate"があります。これはx,y,z軸に対して任意の角度で回転させることのできる機能です。
それではpiperidineをY軸に対して0,30,60,90°と回転させてみます。

 $ mol2ps --rotate=0,0,0 piperidine.mol > angle0.ps
 $ mol2ps --rotate=0,30,0 piperidine.mol > angle30.ps
 $ mol2ps --rotate=0,60,0 piperidine.mol > angle60.ps
 $ mol2ps --rotate=0,90,0 piperidine.mol > angle90.ps

上記コマンドの実行時には、"--linewidth=1.5 --fontsize=16"も付加していますが、表示幅の都合上削除しています。

rotateP.png

このような機能は、紙面上で立体化学の解説をしたいときなど役にたつと思います。


banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 06:29| Comment(0) | TrackBack(0) | 分子グラフィクス | このブログの読者になる | 更新情報をチェックする

2006年06月29日

mol2ps:分子ファイルをpsファイルへ

mol2psは、sdf形式等の分子ファイルをpsファイルに変換するソフトウェアです。データベースに画像ファイルとして構造式を格納したり、印刷物に構造式を載せたりと様々な用途で利用できると思います。

インストールはとても簡単で、Windowsの場合ですと、mol2psのHPから"mol2ps-latest-win32.zip"をdownloadし、展開するだけです。
展開するとmol2ps.exeができます。

mol2psには以下の様な便利なオプションがありますので、好みの画像になるように調節できます。

Usage: mol2ps [options]
where is the file containing the molecular structure
(supported formats: MDL *.mol or *.sdf, Alchemy *.mol, Sybyl *.mol2)
if is "-" (without quotes), the program reads from standard input

valid options are:
--font=, default: Helvetica
--fontsize=, default: 14
--fontsizesmall=, default: 9 (for subscripts)
--linewidth=, default: 1.0 (linewidth in points; use 1 decimal)
--rotate=, default: auto (n,n,n specifies the
angles to rotate the molecule around the X, Y, and Z axis (in degrees)
--autoscale=, default: on (scales the molecule to fit the natural
C-C bond length)
--striphydrogen=, default: on (strips all explicit H atoms)
--hydrogenonhetero=, default: on (adds H to all hetero atoms)
--hydrogenonmethyl=, default: on (adds H to all methyl C atoms)
--hydrogenonstereo=, default: on (shows H if bond is "up" or "down")
--showmolname=, default: off (prints name above the structure)

それでは、実際に実行してみます。

$ mol2ps --fontsize=16 --linewidth=1.5 compound.mol > compound.ps

pngへの変換は、ImageMagickのconvertコマンドなどを用いて変換できます。
$ convert compound.ps compound.png

compound.png

かなりきれいな画像だと思います。また、--hydrogenonheteroオプション等で水素原子の付加の制御ができるところも魅力ですね。


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


posted by わばのり at 00:21| Comment(2) | TrackBack(0) | 分子グラフィクス | このブログの読者になる | 更新情報をチェックする

2006年06月28日

JOELibでRDFコードの計算

RDFコードは3次元descriptorの1つで、RDF(動径分布関数)に基づいて考案されたものです。JOELibではRDFコードを計算するためにRadialDistributionFunctionクラスが用意されています。
基本的な使い方は次のようになります。

RadialDistributionFunction rd = new RadialDistributionFunction();
DescResult result = rd.calculate(mol);
double v[] = ((APropDoubleArrResult)result).value;

最終的には、double[]として結果を受け取るのですが、その際、
(APropDoubleArrResult)をどのようにして知るかというと、RadialDistributionFunctionクラスのgetDescInfo()メソッドとDescriptorInfoクラスのgetResult()メソッドを次のように使うとできます。

DescriptorInfo di = rd.getDescInfo();
System.out.println(di.getResult());
出力:joelib.desc.result.APropDoubleArrResult

RDFコードで使われている原子の物性はデフォルトではGasteiger_Marsiliになっています。
RDFコードは、分子の回転や並進に依存しない3次元descriptorですので、分子の重ね合わせをしなくとも分子間のRDFを比較することができ、便利です。
個人的にお気に入りのdescriptorですので、またメモしたいと思ってます。


banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 08:32| Comment(0) | TrackBack(0) | JOELib | このブログの読者になる | 更新情報をチェックする

2006年06月27日

Chemical MIME

MIME(Multipurpose Internet Mail Extension)とは、従来電子メールで半角英文字といくつかの記号しか利用できなかった規格を画像データなど文字以外のデータを扱えるように拡張した規格です。現在では、Webサーバーとブラウザ間においてもMIMEを用いてデータタイプの判別が行われています。
MIMEタイプは「タイプ名/サブタイプ名」の文字列で表されています。

分子データファイルにもMIMEタイプが存在し、Chemical MIMEと呼ばれているようです。いくつかの例を示してみます。

chemical/x-pdb
chemical/x-mdl-molfile
chemical/x-xyz
chemical/x-gaussian

ここで、x-は正式なサブタイプ名が与えられていないことを示しています。

MozillaのPluginDocを見てみますと、以下のような対応になっていました。

MIMEタイプ プラグイン   
chemical/x-cdx CambridgeSoft ChemDraw
chemical/x-chemdraw CambridgeSoft ChemDraw
chemical/x-csml MDL Chime
chemical/x-gaussian-cube MDL Chime
chemical/x-gaussian-input MDL Chime
chemical/x-jcamp-dx MDL Chime
CambridgeSoft ChemDraw
chemical/x-mdl-molfile MDL Chime
chemical/x-mdl-rxnfile MDL Chime
chemical/x-mdl-tgf MDL Chime
chemical/x-mopac-input MDL Chime
chemical/x-pdb MDL Chime
chemical/x-spectra CambridgeSoft ChemDraw
chemical/x-xyz

分子データファイルを自分のサーバ上に置く場合などにChemical MIMEの設定に注意しようと思っています。


banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 08:46| Comment(0) | TrackBack(0) | その他ツール | このブログの読者になる | 更新情報をチェックする

2006年06月26日

Bioclipse

BioclipseEclipseをベースとしたケム/バイオインフォマティクスのビジュアルプラットフォームです。Eclipseと同様に、新たな機能をpluginとして組み込みことができ、CDK-plugin、Jmol-pluginなどがすでに含まれています。なかなか便利です。

bioclipse.PNG

それにしてもEclipseはすばらしいですね。私はEclipseがないとJavaでプログラムする気すら起こりません。

近いうちに、Bioclipseのpluginを作成したいと思っています。


banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 08:40| Comment(0) | TrackBack(0) | 分子グラフィクス | このブログの読者になる | 更新情報をチェックする

2006年06月25日

JOELibの分子Viewer

JOELibでは、Viewerクラスを用いることで簡単に分子Viewerを利用することができます。Viewerクラスのインスタンスを生成する際、コンストラクタに分子ファイル名を渡すことにより、分子を表示させることができます。ちなみにUIManagerはアプリケーションの見た目を変えるクラスです。

 try
 {
  UIManager.setLookAndFeel
   (UIManager.getSystemLookAndFeelClassName());
 }
 catch (Exception e)
 {
  e.printStackTrace();
 }
 Viewer view = new Viewer(args[0]);

次のような概観となります。
viewer1.png

menuからViewを選ぶと、表示形式及び各種物性情報などを原子の色として視覚的に表現することができます。私の環境では、若干動作が不安定なところが見られますが、何かと使える場面があると思います。


viewer2.png


banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 18:21| Comment(0) | TrackBack(0) | JOELib | このブログの読者になる | 更新情報をチェックする

2006年06月24日

Smalltalkベースの分子グラフィクスライブラリ

5,6年前でしょうか、「じゅん」というSmalltalkのグラフィクスライブラリをHPで見ました。かっこいいグラフィクスライブラリだな、自分でも使いこなしたいなと思い、まずはSmalltalkを覚えなきゃと「Smalltalk イディオム」を買いました。



本を買ったのはいいのですが、そのうちやろうという弱い決意は全く遂行されず、新品に近い状態で本棚に眠っています。

ところが、昨年末に国立情報学研究所から「ケモじゅん」という「じゅん」をベースとした化学系グラフィックスライブラリが公開されていました。
使用手引書やマニュアルを読んでみましたが、丁寧に書かれており、すばらしいと思いました。しかも国産ですしね。

これを機会に、Smalltalkをちゃんと覚えようと新たな強い?決意が芽生えました。あと、Jun for Javaというのもあるみたいです。こちらも楽しそう。



banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 16:04| Comment(0) | TrackBack(0) | 分子グラフィクス | このブログの読者になる | 更新情報をチェックする

2006年06月23日

Pymolの小技

タンパク質の活性部位に2つの異なる化合物を同時に表示させたいことがよくあると思います。例えば、PDBファイルに

 タンパク質
 TER
 化合物1
 TER
 化合物2

と記述し、PymolでOpenしてみると下図のように2つの化合物が融合してしまい何じゃこれっとなってしまいます。他の分子Viewerの中にはTERで区切ってさえいればきちんと認識してくれるものもあります。

sp2.png

単純な解決方法として、

 タンパク質
 TER
 分子A

というファイルと、

 分子B

というファイルを2つ作り、1つめのファイルをOpenしたあと、続けてもう1つのファイルをOpenすれば2つの分子がきちんと分離されます。
sp1.png


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

posted by わばのり at 08:18| Comment(0) | TrackBack(0) | 分子グラフィクス | このブログの読者になる | 更新情報をチェックする

2006年06月22日

JOELibでDescriptorの計算

JOELibのDescriptorsは、"Native values","Atom properties","Fingerprints","Transformations"の4つのカテゴリーに分けられています。"Native values"は、Drug-likenessフィルタやQSARで利用可能な基本的なDescriptorsの集まりですので、今回は、この"Native values"の取得方法についてメモしたいと思います。

"Native values"に属するDescriptor名を知るには、DescriptorHelperクラスのインスタンスを生成し、getNativeDescs()メソッドを実行します。

 Vector nativeDes = DescriptorHelper.instance().getNativeDescs();
 String desName;
 for (int i = 0; i < nativeDes.size(); i++)
 {
  desName = (String)nativeDes.get(i);
  System.out.println(desName);
 }

ここで得られるdesNameがDescriptorの値を取得する際に必要となります。

実際に、Descriptorの値を取得するには、DescriptorFactoryクラスのgetDescriptor(String name)メソッドでDescriptorを取得し、calculate(JOEMol mol)メソッドによりDescriptorの値の計算を行います。

 Descriptor descHBA = DescriptorFactory.getDescriptor("Number_of_HBA_2");
 DescResult results = descHBA.calculate(mol);

Descriptorの値をint型にするには次のように記述します。

 int hba =((IntResult)results).getInt();

logPのようにdouble型のDescriptorの値をもつものは、次のように記述します。
 
 double logP =((DoubleResult)results).getDouble();

CDKと同様に、JOELibでもDescriptorの値を取得すれば、簡単にrule of fiveフィルタが作成でき便利です。
ちなみにJOELibでのlogPの計算は、以下の論文の手法が使われています。

 J. Chem. Inf. Comput. Sci., 1999, 39, 868-873.


banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 08:29| Comment(0) | TrackBack(0) | JOELib | このブログの読者になる | 更新情報をチェックする

2006年06月21日

JOELibでPOVRay

JOELibを用いて分子をPOVRay形式で出力してみます。POVRayは、レイトレーシング法により高品質な3次元画像を生成するフリーな3Dグラフィックスツールです。JOELibはPOVRayのシーン記述ファイルである「POVファイル」を出力することができるため、センスがあれば(私にはありませんが)、非常にかっこいい分子画像を作成することができます。

まずは、POVRayクラスのインスタンスを生成し、initWriter(OutputStream os)メッソドを実行します。

 joelib.io.types.POVRay ray = new joelib.io.types.POVRay();
 ray.initWriter(output);

次に分子の表示形式のセットを行います。
 
 ray.setOutputType(joelib.io.types.POVRay.SPHERE);

ここでは、SPHERE表示を選択しています。その他の表示形式は、POVRayクラスのField Summaryをご覧ください。
POVファイルの出力はwrite(JOEMol mol)メソッドで行います。
 
 ray.write(mol);

出力ファイルをPOVRayにより描画してみます。

SPHERE形式
SPHERE.png

Ball and Stick形式
BandS.png
 
POVRayで描画させると真っ先にやりたくなる操作は、縮小・拡大/回転だと思います。POVファイルの一番下にある以下のような部分にrotate<>及びscale<>を書き加えることにより実現できます。

object { Molecule1
   rotate<90,0,0>
   scale<0.8,0.8,0.8>
}


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



posted by わばのり at 07:51| Comment(0) | TrackBack(0) | JOELib | このブログの読者になる | 更新情報をチェックする

2006年06月20日

PDBの二次データベース

現在、Protein Data Bank(PDB)には、約37000の立体構造が登録されているようです。年ごとの登録数をみると指数関数的に増加していますね。PDBの成熟に伴い、多くの二次データベースが派生していると思いますが、私自身は低分子とタンパク質の複合体構造に興味があるので、これに関連しそうなデータベースをメモしておきます。

MSDChem
 PDB中のLigandを部分構造,Fragments,Fingerprintなどで検索ができる。

The PDBbind Database
 Protein-ligand複合体の実測binding affinityデータが収められている。

Computed Ligand Binding Energy(CLiBE) Database
 コンピュータで計算したprotein-ligand複合体の相互作用エネルギーが収められている。

HIV Protease Database:
 HIV-1,HIV-2,SIVの立体構造が整理されて収められている。

Protein Ligand Database (PLD):
 Protein-ligand複合体のbinding affinityデータが実測/計算値共に収められている。



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









posted by わばのり at 08:21| Comment(0) | TrackBack(0) | その他ツール | このブログの読者になる | 更新情報をチェックする

2006年06月19日

JOElibで水素原子の付加

JOElibでは、JOEMolクラスのaddHydrogens()メソッドにより、H原子を付加することができます。H原子を付加する際、基本的には、以下のコードにより、H原子の個数を決定しています(実際のJOELib中のコードです)。

 hcount = atom.getImplicitValence() - atom.getValence();

JOELibの便利な機能として、pHを考慮したH原子の付加ができます。具体的には、joelib/data/plain/phmodel.txtに従っています。例えば、carboxylic acidの場合、以下のように記述されています。

TRANSFORM O=CO[#1:1] >> O=CO
TRANSFORM O=C[OQ1-0:1] >> O=C[O-:1]

この記述によれば、COOHはCOO-に変換されるため、負に荷電しており、H原子は付加されません。したがって、phmodel.txtをカスタマイズすることにより、H原子の付加を自由にコントロールできます。記述方法は、Reaction SMARTSだと思います。

もう1つ便利な機能として、polar hydrogensのみの付加があります。例えば、水素結合のパターンを目視で認識する際に便利ですし、Amberにもこれに対応した力場パラメータが存在するので、分子モデリングでも利用することができます。

それでは、JOEMolクラスのaddHydrogens(boolean polaronly, boolean correctForPH)メソッドを用いて、上述した内容をアラニンを用いて確認したいと思います。第一引数が、polar hydrogensのみのH原子付加を行うか、第二引数は、phmodel.txtを用いたH原子付加のコントロールを行うかを示しております。

hbond.png
上図中のtrue,falseは第一引数、第二引数の順番になっております。ちなみに、addHydrogens()メソッドでは、false,trueとなっております。


banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 08:10| Comment(0) | TrackBack(0) | JOELib | このブログの読者になる | 更新情報をチェックする

2006年06月18日

頭の体操

昨日本屋でなんとなく本を見ていると多湖輝さんの「頭の体操」シリーズが平積みで置いてあった。小中学生の頃、このシリーズが大好きで、難解な問題をクリアーできたときの爽快感がたまらなかった。それまで知っていた“なぞなぞ”や“パズル”とは、明らかに一線をきす問題に、これを考えつく多湖さんは、当時、私の尊敬する大人のNo.1であった。もう1つこのシリーズの印象に残っているのは、本の表紙で、具体的には思い出せなかったのでネットで探してみると、ココにあった。このなにか特異な絵柄も惹きつけられる要素だったのだと思う。第10集あたりまでは読んだ記憶があるが、だんだん興味が他に移っていき読まなくなった。

あれから20年程たったが、多湖さんにもう一度挑戦してみようと思う。DS脳トレで鍛えているおじさんは、20年たっても脳年齢は21才なはず(笑)なので8割以上はクリアーできるだろう(といいんなあ)と信じている。





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


posted by わばのり at 07:28| Comment(0) | TrackBack(0) | 日記 | このブログの読者になる | 更新情報をチェックする

2006年06月17日

JOELibでバベる!

JOELibを使って簡易babelを作成してみます。元々babel自体がOELibで作られていますので、JOELibでbabelっぽいプログラムを作るのは容易だろうと想像できます。

JOELibではSimpleReaderクラスを用いてファイルの入力を行います。
SimpleReaderクラスのインスタンスの生成時に、ファイル形式の情報が必要となるため、IOTypeHolderクラスを用いて、その情報を取得します。

 IOType inType = IOTypeHolder.instance().filenameToType(inputf);

ここでは、インスタンスの数を1つに保つために、GoFのSingletonパターンが使われていると思います。

 sreader = new SimpleReader(inputs,inType);

また、ファイルの出力は、SimpleWriterクラスを用いて行います。

 IOType outType = IOTypeHolder.instance().filenameToType(outputf);
 swriter = new SimpleWriter(outputs,outType);

そして、JOEMolクラスのインスタンスの生成を次のように行います。
 
 JOEMol mol = new JOEMol(inType,outType);

分子の入力は、readNextメソッドで行い、出力はwriteNextメソッドで行っています。

 while(true){
  if(!sreader.readNext(mol))
    break;
  swriter.writeNext(mol);
 }

実行
# java jbabel inputfile outputfile

拡張子で、ファイル形式の判定を行っていますので、とてもシンプルになりました。ファイル変換だけをみればbabelと同等の機能をもっていると思います。

以下に今回のソースコードを示しております。変換可能なファイル形式は、ココに記載されています。

簡単にbabelもどきができましたね(当たり前ですか..)。個人的にJOELibのクラス構成をとても気に入っていますので、ちょくちょくメモする予定です。

Jbabel.java



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


posted by わばのり at 07:16| Comment(0) | TrackBack(0) | JOELib | このブログの読者になる | 更新情報をチェックする

2006年06月16日

CDKで"rule of five"フィルタの作成 2

RuleOfFiveDescriptorクラスには、以下6つのメソッドがあります。

 calculate(AtomContainer mol)
 getParameterNames()
 getParameters()
 getParameterType(java.lang.String name)
 getSpecification()
 setParameters(java.lang.Object[] params)

APIの説明を読むと、getParameterNames()及びTypeあたりでクラス内で使われているdescriptorを知り、setParameters(java.lang.Object[] params)でそのdescriptorの値を任意に指定できるのかなと思いました。デフォルトでは、Lipinski's Ruleだろうけど、Mw<350などに設定したい場合があるので結構親切じゃんっと思い込んだわけです。ただ、calculate(AtomContainer mol)の説明に以下のような不可解な文書が記載されています。

 the method take a boolean checkAromaticity: if the boolean is  
 true, it means that aromaticity has to be checked

checkAromaticityがtrueならaromaticityのチェックをするんだなということは分かりますが、どうやって設定するの?と思ってしまいます。aromaticityのチェックはしておかないと特にXLogPの計算結果に影響しそうなので...もしかしてデフォルトでFALSEではないよな...など心配になってきたので、ソースコードを見ることにしました。

まずは、setParameters(java.lang.Object[] params)を見ると次のようなコードが書かれていました。

 if (params.length > 1) {
   throw new CDKException("RuleOfFiveDescriptor only expects one
   parameter");
  }
if (!(params[0] instanceof Boolean)) {
throw new CDKException("The first parameter must be of type
Boolean");
}
checkAromaticity = ((Boolean) params[0]).booleanValue();

えっ!このメソッドからcheckAromaticityがセットできることは分かりましたが、いろんなパラメータのセットができるという私の推測は見事にはずれてしまいました。ということで、以下のコードでcheckAromaticityをTRUEにできます。

 // checkAromaticityをTRUEにする
 Boolean param[] = new Boolean[1];
 param[0] = Boolean.TRUE;
 try {
  rule5.setParameters(param);
 } catch (CDKException e1) {
  System.err.println(e1.toString());
 }

本当は、コード作成者は、複数のパラメータをセットできるようにしたかったのではないかな.....
ちなみにcheckAromaticityはデフォルトでは、falseですので、注意してください。あと、当然ですが、getParameterNames()でも返ってくるのは、"checkAromaticity"だけです。

ということで、RuleOfFiveDescriptorクラスは、checkAromaticityをTRUEかFALSEに設定して、calculate()すればよい(それしかない)と思われます。

前回も記載しましたが、"rule of five"フィルタのコアになるコードは以下のようになります。
 
 //RuleOf5の計算
 dv = rule5.calculate((AtomContainer)mol);
 IntegerResult ir = (IntegerResult)dv.getValue();

 //RuleOf5を満たすものをファイル出力
 if(ir.intValue() == 0){
  mw.writeMolecule(mol);
  count++;
 }

これで、"rule of five"フィルタの完成です。
注意点としては、CDKにおける水素結合のドナー、アクセプターの定義を確認しないとせっかくの有望な化合物がこぼれてしまう可能性があること。また、CDKのlogPの計算方法(精度)も把握しておかないと同様の心配がありますので、このあたりもそのうち取り上げてみます。

CDKex2.java



banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 08:03| Comment(0) | TrackBack(0) | CDK | このブログの読者になる | 更新情報をチェックする

2006年06月15日

CDKで"rule of five"フィルタの作成 1

CDKには、"rule of five"フィルタを簡単に構築できるRuleOfFiveDescriptorクラスがあります。"rule of five"とは、Lipinskiにより提案された化合物のDrug-likeness(薬らしさ)を評価する方法の1つです。

RuleOfFiveDescriptorクラスでは、以下のように"rule of five"の各項目を満たさなければlipinskifailuresに1を加えています。すなわち、"rule of five"を全て満たした化合物は、lipinskifailuresが0となるわけです。

 if(xlogPvalue > 5.0) { lipinskifailures += 1; }
 if(acceptors > 10) { lipinskifailures += 1; }
 if(donors > 5) { lipinskifailures += 1; }
 if(mwvalue > 500.0) { lipinskifailures += 1; }
 if(rotatablebonds > 10.0) { lipinskifailures += 1; }

今回は、sdf形式のファイルを入力し、"rule of five"を満たす化合物のファイルを出力するプログラムを作成したいと思います。

 Usage: CDKex2 data.sdf data_passed.sdf

まずは、概略ですが、data.sdfに含まれる化合物を順次取り込むために、前回用いたMDLReaderクラスではなく、IteratingMDLReaderクラスを用います。

 IteratingMDLReader itr = new IteratingMDLReader(fr);
 while(itr.hasNext())
   mol =(Molecule)itr.next();

sdf形式のファイル中にH原子が付加されていない場合は、HydrogenAdderクラスを用いて付加できます。

 HydrogenAdder adder = new HydrogenAdder();
 adder.addExplicitHydrogensToSatisfyValency(mol);

"rule of five"の評価はRuleOfFiveDescriptorクラスのcalculate()メソッドにより実施し、getValue()メソッドを用いて前述したlipinskifailuresの値を取得します。

 RuleOfFiveDescriptor rule5 = new RuleOfFiveDescriptor();
 dv = rule5.calculate((AtomContainer)mol);
 IntegerResult ir = (IntegerResult)dv.getValue();
 if(ir.intValue() == 0)
   mw.writeMolecule(mol);

RuleOfFiveDescriptorクラスはちょっとクセ(問題)のあるクラスですので、これも含め、次回、実際のコードを用いてプログラムの解説をしたいと思います。



banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 08:49| Comment(0) | TrackBack(0) | CDK | このブログの読者になる | 更新情報をチェックする

2006年06月14日

分子動力学

普段、Amber/CHARMM/NAMDなどを研究に利用していると、一度は、分子動力学(MD)を自力でプログラミングしたい。もしくは、MDの中身をより詳しく理解したい。このような話をよく聞きます。
私は、学生のときMDを自力で開発していましたが、結局、私の実力ではAmberには勝てないなっと心底悟らされ、今ではMDを与えられた課題を解くツールとして、いかに使いこなすかで満足しています。パラメータを検討した(いじりまわした)末に、実測値と相関が取れたりするとマジっか?と思いながらもうれしいものです。

MD関連の書籍は、多くでていると思いますが、MDの理論体系を見通せる書籍としては、以下の2冊がお勧めです。




実際に、MDをプログラミングするとっかかりとして、非常に参考になったのは、以下の書籍でした。FORTRANで書かれたソースプログラムが掲載されていますので、Javaとかで書き直しただけで、何となくMDを自作した気分になれます(笑)。ただし、実際にプログラミングするといろいろな問題点が把握できると思いますので、MDを利用する上でも役に立つと思います。



banner_02.gif
人気ブログランキング(クリックして応援してね)
posted by わばのり at 07:23| Comment(0) | TrackBack(0) | 書籍紹介 | このブログの読者になる | 更新情報をチェックする

2006年06月13日

Pythonでケムインフォ:Fingerprint 2

FrownsでのFingerprintの使い方をメモしたいと思います。
 
 import frowns.Fingerprint

によりモジュールを読み込み、generateFingerprint関数によりFingerprintを生成します。

 def generateFingerprint(molecule, numInts=32, pathLength=7):

デフォルトでは、ビット長は32、部分構造のパス長が7となっています。
ここでは、生成したFingerprintを部分構造検索の前工程に利用するプログラムを作成致します。

プログラムの概要ですが、

1.
クエリー構造のFingerprintを作成する。ただし、ここでは、Smiles.smilinで変換可能なSMARTSパターンを入力すること。

2. #Fingerprint generation
標的となる化合物郡のFingerprintを全て生成し、dbリストに収める。

3. #Smarts Search
第3引数であるflagが1ならSMARTS検索の前工程としてFingerprint同士の一致を評価し、0なら評価しない。
その後、SMARTS検索。

実行例1:Fingerprintを利用しない。
$ python smartsearch2.py data.sdf "c1ccccc1OCCCCC" 0

出力:
C(=O)(O)c1c(OC(=O)CCCCC)cccc1
O(CCCCCCCCCCCC)c1ccc(C(=O)O)cc1
Time 4.43799996376
2 Hit compounds

実行例2:Fingerprintを利用する。
$ python smartsearch2.py data.sdf "c1ccccc1OCCCCC" 1

出力:
C(=O)(O)c1c(OC(=O)CCCCC)cccc1
O(CCCCCCCCCCCC)c1ccc(C(=O)O)cc1
Time 1.84400010109
2 Hit compounds

この例では、Fingerprintの利用により約2.4倍、検索速度のUPが確認できました。実際の検索システムでは、Fingerprintは化合物のデータベースへの登録時に、ビット文字列などで格納されることが多いと思います。

smartsearch2.py



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



 
posted by わばのり at 08:06| Comment(0) | TrackBack(0) | Frowns | このブログの読者になる | 更新情報をチェックする

2006年06月12日

Pythonでケムインフォ:Fingerprint

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による部分構造検索プログラムを作成したいと思います。


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





posted by わばのり at 07:11| Comment(0) | TrackBack(0) | Frowns | このブログの読者になる | 更新情報をチェックする

2006年06月11日

CDK:出力形式について

昨日の記事で、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出さないみたい....



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

posted by わばのり at 07:29| Comment(0) | TrackBack(0) | CDK | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。


×

この広告は1年以上新しい記事の投稿がないブログに表示されております。