2017-05-31

CDK 2.0

CDK

CDK ver 2.0 が公開されているので、以前作ったプログラムの動作確認も兼ねて、サンプルスクリプトと共に紹介します。
まずは、MACCS keys (166ビット版)から。

○変更点
・IMolecule インターフェースが削除されたので、IAtomContainer に変更
・IteratingSMILESReader 関数の引数に、空の IChemObjectBuilder を追加。

import java.io.*;
import org.openscience.cdk.interfaces.*;
import org.openscience.cdk.io.iterator.IteratingSMILESReader;
import org.openscience.cdk.DefaultChemObjectBuilder;

import java.util.BitSet;
import org.openscience.cdk.fingerprint.MACCSFingerprinter;

class maccs {
    public static void main(String args[]){

        if(args.length!=1){
            System.err.println("maccs <smiles-file>");
            System.exit(1);
        }
        FileInputStream fis = null;
        IteratingSMILESReader isr = null;
        IChemObjectBuilder builder = null;
        try{
            builder = DefaultChemObjectBuilder.getInstance();
            fis = new FileInputStream(new File(args[0]));
            isr = new IteratingSMILESReader(fis, builder);
        } catch (Exception e) {
            e.printStackTrace();
        }

        while( isr.hasNext() ){
            IAtomContainer imol = isr.next();
            String id = (String)imol.getProperty("cdk:Title");

            MACCSFingerprinter ma = new MACCSFingerprinter();
            String fp = "";
            try{
                BitSet fing = ma.getFingerprint(imol);
                for(int i=0; i<ma.getSize(); i++){
                    if( fing.get(i) ){
                        fp = fp + ",1";
                    } else {
                        fp = fp + ",0";
                    }
                }
            } catch (Exception e1) {
                e1.printStackTrace();
            }
            System.out.printf("%s%s\n", id, fp);
        }
    }
}

参考：　旧スクリプト（cdk-1.2）

2012-06-10

PubChem Fingerprint

CDK fingerprint

PubChemデータベースは、独自のフィンガープリントを構築して
Webサイト内での類似度検索に使っています。

そのフィンガープリントの構成は既に公開されているのですが、
881ビットもあるので、ちまちま実装するのは面倒です。

そこでCDKを見てみると…
PubchemFingerprinterクラスが実装されており、パパっと計算できるようです。



import java.io.*;

import org.openscience.cdk.interfaces.*;

import org.openscience.cdk.io.iterator.IteratingSMILESReader;

import org.openscience.cdk.interfaces.IMolecule;
import java.util.BitSet;

import org.openscience.cdk.fingerprint.PubchemFingerprinter;
class pubchemfp {

    public static void main(String args[]){
        if(args.length!=1){

            System.err.println("pubchemfp <smiles-file>");

            System.exit(1);

        }

        FileInputStream fis = null;

        IteratingSMILESReader isr = null;

        try{

            fis = new FileInputStream(new File(args[0]));

            isr = new IteratingSMILESReader(fis);

        } catch (Exception e) {

            e.printStackTrace();

        }
        while( isr.hasNext() ){

            IMolecule imol = (IMolecule)isr.next();

            String id = (String)imol.getProperty("cdk:Title");
            PubchemFingerprinter pcfp = new PubchemFingerprinter();

            String fp = "";

            try{

                BitSet fing = pcfp.getFingerprint(imol);

                for(int i=0; i<pcfp.getSize(); i++){

                    if( fing.get(i) ){

                        fp = fp + ",1";

                    } else {

                        fp = fp + ",0";

                    }

                }

            } catch (Exception e1) {

                e1.printStackTrace();

            }

            System.out.printf("%s%s\n", id, fp);

        }

    }

}

2012-03-29

chemical beauty (2)

(1) のつづき。

Bickertonら*1は、QEDを定義したあと、薬物候補を判別するベンチマークを行いました。
従来から使われているリピンスキーの法則より優れていることを示すためです。

正例：DrugBank登録の薬物
負例：PDB登録の低分子リガンド

ここで「QEDの方が高い分離能を示しました。」だけで締めくくらないのが、
この論文のチャームポイントです。

さて、製薬化学者は、化学構造を目で見て、直感的にその良し悪しを判断できます。
これはコンピュータでは対応が困難であり、一種の職人芸と言って良いと思います。
Bickertonらは、これを「Chemical aesthetics」と呼び、
製薬化学者の判断とQEDの出力との関連性をさらに調べました。

すると、

化学者が魅力的と感じた化合物のQED平均値：　0.67
化学者が魅力的でないと感じた化合物のQED平均値：　0.49

という結果を得たそうです。こういう評価を丁寧に行うことも大切ですよね。

*1: Bickerton et al. Quantifying the chemical beauty of drugs. Nat. Chem. (2012) 4, 90-98 PubMed

2012-02-05

芳香環の数

CDK descriptor

芳香環の数は、分子の疎水性に大きく寄与します。
疎水性は、アルブミン結合率・CYP3A4代謝酵素の阻害率・hERGタンパク質阻害率などに強く影響します。
したがって、芳香環を数えることは、「薬らしさ」の指標となります*1。

数えてみましょう。
ただし、ナフタリンはベンゼン環2個と数えます。
すなわち、最小の環構造に分解した後、芳香性を確認します。

これは難解そうに見えますが、
SSSRFinderという便利なクラスで一発解決です。



/* narring.java */
import java.io.*;

import org.openscience.cdk.interfaces.*;

import org.openscience.cdk.io.iterator.IteratingSMILESReader;

import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import org.openscience.cdk.ringsearch.SSSRFinder;

import org.openscience.cdk.CDKConstants;

import org.openscience.cdk.tools.manipulator.RingSetManipulator;
class narring {
    public static void main(String args[]){
        if(args.length!=1){

            System.err.println("narring <smiles-file>");

            System.exit(1);

        }

        FileInputStream fis = null;

        IteratingSMILESReader isr = null;

        try{

            fis = new FileInputStream(new File(args[0]));

            isr = new IteratingSMILESReader(fis);

        } catch (Exception e) {

            e.printStackTrace();

        }

        while( isr.hasNext() ){

            IMolecule imol = (IMolecule)isr.next();

            String id = (String)imol.getProperty("cdk:Title");

            int nar = 0;

            try{

                RingSetManipulator rsm = new RingSetManipulator();

                SSSRFinder sssrf = new SSSRFinder(imol);

                IRingSet rs = sssrf.findSSSR();

                rsm.markAromaticRings(rs);

                IRing r = null;

                for(int i=0;i<rs.getAtomContainerCount();i++){

                    r = (IRing) rs.getAtomContainer(i);

                    boolean isAromatic = r.getFlag(CDKConstants.ISAROMATIC);

                    if (isAromatic) nar++;

                }

            } catch (Exception e1) {

                e1.printStackTrace();

            }

            System.out.printf("%s\t%d\n", id, nar);

        }

    }

}

*1: Ritchie et al. The impact of aromatic ring count on compound developability--are too many aromatic rings a liability in drug design? Drug Discov. Today (2009) 14, 1011-1020 PubMed

2012-01-26

chemical beauty (1)

CDK descriptor

「薬らしさ（drug-likeness）」という指標があります。
これは、物性が承認薬と似ているという基準で候補化合物をざくっと絞り込むときに使います。
有名な指標は、リピンスキーの法則ですが、最近、新たな指標が提案されたので取り上げます。

それは、Bickertonら*1が報告した QED (quantitative estimate of drug-likeness) という計算手法です。
次の8つの記述子の分布をもとに計算されます。

分子量
疎水性（ALOGP）
水素結合ドナーの数
水素結合アクセプターの数
極性表面積
回転可能結合の数
芳香環の数
望ましくない部分構造の数（ALERTS）

さっそく実装に挑戦です。
まず、ALERTSからいきましょう。
smartsパターンをString列で定義して、
SMARTSQueryToolクラスを使って数えます。



import java.io.*;

import org.openscience.cdk.interfaces.*;

import org.openscience.cdk.io.iterator.IteratingSMILESReader;

import org.openscience.cdk.smiles.smarts.SMARTSQueryTool;
public class alerts {
    private static final String[] smarts = {

        "*1[O,S,N]*1", "[S,C](=[O,S])[F,Br,Cl,I]",

        "[CX4][Cl,Br,I]", "[C,c]S(=O)(=O)O[C,c]",

        "[$([CH]),$(CC)]#CC(=O)[C,c]", "[$([CH]),$(CC)]#CC(=O)O[C,c]",

        "n[OH]", "[$([CH]),$(CC)]#CS(=O)(=O)[C,c]",

        "C=C(C=O)C=O", "n1c([F,Cl,Br,I])cccc1",

        "[CH1](=O)", "[O,o][O,o]", "[C;!R]=[N;!R]", "[N!R]=[N!R]",

        "[#6](=O)[#6](=O)", "[S,s][S,s]", "[N,n][NH2]", "C(=O)N[NH2]", "[C,c]=S",

        "[$([CH2]),$([CH][CX4]),$(C([CX4])[CX4])]=[$([CH2]),$([CH][CX4]),$(C([CX4])[CX4])]",

        "C1(=[O,N])C=CC(=[O,N])C=C1", "C1(=[O,N])C(=[O,N])C=CC=C1",

        "a21aa3a(aa1aaaa2)aaaa3", "a31a(a2a(aa1)aaaa2)aaaa3",

        "a1aa2a3a(a1)A=AA=A3=AA=A2", "c1cc([NH2])ccc1",

        "[Hg,Fe,As,Sb,Zn,Se,se,Te,B,Si,Na,Ca,Ge,Ag,Mg,K,Ba,Sr,Be,Ti,Mo,Mn,Ru,Pd,Ni,Cu,Au,Cd,Al,Ga,Sn,Rh,Tl,Bi,Nb,Li,Pb,Hf,Ho]", "I",

        "OS(=O)(=O)[O-]", "[N+](=O)[O-]",

        "C(=O)N[OH]", "C1NC(=O)NC(=O)1", "[SH]", "[S-]",

        "c1ccc([Cl,Br,I,F])c([Cl,Br,I,F])c1[Cl,Br,I,F]", "c1cc([Cl,Br,I,F])cc([Cl,Br,I,F])c1[Cl,Br,I,F]",

        "[CR1]1[CR1][CR1][CR1][CR1][CR1][CR1]1", "[CR1]1[CR1][CR1]cc[CR1][CR1]1",

        "[CR2]1[CR2][CR2][CR2][CR2][CR2][CR2][CR2]1",

        "[CR2]1[CR2][CR2]cc[CR2][CR2][CR2]1",

        "[CH2R2]1N[CH2R2][CH2R2][CH2R2][CH2R2][CH2R2]1",

        "[CH2R2]1N[CH2R2][CH2R2][CH2R2][CH2R2][CH2R2][CH2R2]1",

        "C#C", "[OR2,NR2]@[CR2]@[CR2]@[OR2,NR2]@[CR2]@[CR2]@[OR2,NR2]",

        "[$([N+R]),$([n+R]),$([N+]=C)][O-]", "[C,c]=N[OH]",

        "[C,c]=NOC=O", "[C,c](=O)[CX4,CR0X3,O][C,c](=O)",

        "c1ccc2c(c1)ccc(=O)o2", "[O+,o+,S+,s+]",

        "N=C=O", "[NX3,NX4][F,Cl,Br,I]",

        "c1ccccc1OC(=O)[#6]", "[CR0]=[CR0][CR0]=[CR0]",

        "[C+,c+,C-,c-]", "N=[N+]=[N-]",

        "C12C(NC(N1)=O)CSC2", "c1c([OH])c([OH,NH2,NH])ccc1",

        "P", "[N,O,S]C#N", "C=C=O", "[Si][F,Cl,Br,I]",

        "[SX2]O", "[SiR0,CR0](c1ccccc1)(c2ccccc2)(c3ccccc3)",

        "O1CCCCC1OC2CCC3CCCCC3C2", "N=[CR0][N,n,O,S]",

        "[cR2]1[cR2][cR2]([Nv3X3,Nv4X4])[cR2][cR2][cR2]1[cR2]2[cR2][cR2][cR2]([Nv3X3,Nv4X4])[cR2][cR2]2",

        "C=[C!r]C#N",

        "[cR2]1[cR2]c([N+0X3R0,nX3R0])c([N+0X3R0,nX3R0])[cR2][cR2]1",

        "[cR2]1[cR2]c([N+0X3R0,nX3R0])[cR2]c([N+0X3R0,nX3R0])[cR2]1",

        "[cR2]1[cR2]c([N+0X3R0,nX3R0])[cR2][cR2]c1([N+0X3R0,nX3R0])",

        "[OH]c1ccc([OH,NH2,NH])cc1", "c1ccccc1OC(=O)O",

        "[SX2H0][N]", "c12ccccc1(SC(S)=N2)", "c12ccccc1(SC(=S)N2)", "c1nnnn1C=O",

        "s1c(S)nnc1NC=O", "S1C=CSC1=S", "C(=O)Onnn", "OS(=O)(=O)C(F)(F)F",

        "N#CC[OH]", "N#CC(=O)", "S(=O)(=O)C#N", "N[CH2]C#N",

        "C1(=O)NCC1", "S(=O)(=O)[O-,OH]", "NC[F,Cl,Br,I]", "C=[C!r]O",

        "[NX2+0]=[O+0]", "[OR0,NR0][OR0,NR0]",

        "(C(=O)O[C,H1]).(C(=O)O[C,H1]).(C(=O)O[C,H1])", "[CX2R0][NX3R0]",

        "c1ccccc1[C;!R]=[C;!R]c2ccccc2",

        "[NX3R0,NX4R0,OR0,SX2R0][CX4][NX3R0,NX4R0,OR0,SX2R0]",

        "[s,S,c,C,n,N,o,O]~[n+,N+](~[s,S,c,C,n,N,o,O])(~[s,S,c,C,n,N,o,O])~[s,S,c,C,n,N,o,O]",

        "[s,S,c,C,n,N,o,O]~[nX3+,NX3+](~[s,S,c,C,n,N])~[s,S,c,C,n,N]",

        "[*]=[N+]=[*]", "[SX3](=O)[O-,OH]", "N#N",

        "F.F.F.F", "[R0;D2][R0;D2][R0;D2][R0;D2]",

        "[cR,CR]~C(=O)NC(=O)~[cR,CR]", "C=!@CC=[O,S]",

        "[#6,#8,#16][C,c](=O)O[C,c]", "c[C;R0](=[O,S])[C,c]",

        "c[SX2][C;!R]", "C=C=C", "c1nc([F,Cl,Br,I,S])ncc1", "c1ncnc([F,Cl,Br,I,S])c1",

        "c1nc(c2c(n1)nc(n2)[F,Cl,Br,I])", "[C,c]S(=O)(=O)c1ccc(cc1)F"};

    // "[15N,13C,18O,2H,34S]"};
    public static void main(String args[]){
        if(args.length!=1){

            System.err.println("alerts <smiles-file>");

            System.exit(1);

        }

        FileInputStream fis = null;

        IteratingSMILESReader isr = null;

        try{

            fis = new FileInputStream(new File(args[0]));

            isr = new IteratingSMILESReader(fis);

        } catch (Exception e) {

            e.printStackTrace();

        }

        while( isr.hasNext() ){

            IMolecule imol = (IMolecule)isr.next();

            String id = (String)imol.getProperty("cdk:Title");
            int nalert = 0;
            try {

                SMARTSQueryTool sqt = new SMARTSQueryTool("C");

                for (int i = 0; i < smarts.length; i++) {

                    sqt.setSmarts(smarts[i]);

                    if ( sqt.matches(imol) ) {

                        nalert = nalert + 1;

                    }

                }

            } catch (Exception e1) {

                e1.printStackTrace();

            }

            System.out.printf("%s\t%d\n", id, nalert);

        }

    }

}

最後の放射性元素のパターンは、実装が面倒だったのでコメントアウトしました。
放射性化合物・・・ふつうは入っていませんよね。

*1: Bickerton et al. Quantifying the chemical beauty of drugs. Nat. Chem. (2012) 4, 90-98 PubMed

2012-01-22

化合物が流れる動画 (3)

OpenBabel

前回の続き。
こんどは、構造の一部のみを"回転"させます。
OpenbabelのRotateメソッドを使ってみます。
回転可能な結合を適当に見つけてきて、
二面角をなす4原子の番号を関数に入力します。
（この番号は、SDFにおける登場順で付けられます）



#include <stdio.h>

#include <iostream>

#include <fstream>

#include <string>

#include <openbabel/mol.h>

#include <openbabel/obconversion.h>
using namespace std;

using namespace OpenBabel;
#define PI 3.14159265358979323846
void mol_rotate(int a1, int b1, int c1, int d1, int rot_ang, OBMol* comp,

                OBConversion* conv, ofstream* comp_out);
int main(int argc,char *argv[]){

  ifstream comp_in( argv[1] );

  ofstream comp_out( argv[2] );
  OBConversion conv;

  conv.SetInFormat("SDF");

  conv.SetOutFormat("SDF");
  OBMol comp;

  conv.Read(&comp, &comp_in);
  mol_rotate(7, 8, 9, 10, 120, &comp, &conv, &comp_out);

  mol_rotate(22, 23, 24, 25, 180, &comp, &conv, &comp_out);

  mol_rotate(7, 8, 9, 10, 240, &comp, &conv, &comp_out);
}
void mol_rotate(int a1, int b1, int c1, int d1, int rot_ang, OBMol* comp,

                 OBConversion* conv, ofstream* comp_out){

  OBAtom *a, *b, *c, *d;

  a=(*comp).GetAtom(a1);

  b=(*comp).GetAtom(b1);

  c=(*comp).GetAtom(c1);

  d=(*comp).GetAtom(d1);
  double ang = (*comp).GetTorsion(a, b, c, d) * PI / 180.0;
  for (int i=0; i<=rot_ang; i=i+5) {

    double nang= ang + PI / 180 * i;

    char name[10];

    sprintf(name, "%d", i);

    (*comp).SetTorsion(a, b, c, d, nang);

    (*comp).SetTitle( name );

    (*conv).Write(comp, comp_out);

  }

}

2012-01-08

化合物が流れる動画 (2)

OpenBabel

あけましておめでとうございます。
本年もよろしくお願いいたします。

今年こそは、新しいことに挑戦しようと思います。

さて、化合物の動画について、まとめてみます。
まずは、単純な平行移動から。
Openbabel（ver2.3.0）を使います。Translateメソッドで移動できます。
下の例では、ベクトル(1,1,1)の方向に移動しつつ、20コマ分の座標を保存します。



#include <stdio.h>

#include <iostream>

#include <fstream>

#include <string>

#include <openbabel/mol.h>

#include <openbabel/obconversion.h>

#include <openbabel/math/vector3.h>
using namespace std;

using namespace OpenBabel;
int main(int argc,char *argv[]){

  ifstream comp_in( argv[1] );

  ofstream comp_out( argv[2] );
  OBConversion conv;

  conv.SetInFormat("SDF");

  conv.SetOutFormat("SDF");
  OBMol comp;

  conv.Read(&comp, &comp_in);
  for(int i=1; i<=20; i++){

    comp.Translate(VX + VY + VZ);

    conv.Write(&comp, &comp_out);

  }

}

VX ... (x,y,z) = (1,0,0)
VY ... (x,y,z) = (0,1,0)
VZ ... (x,y,z) = (0,0,1)

次に、このSDFをPyMOL(ver.1.1r1)で読み込み、Movieとして保存します。
連番のPNG形式ファイルが出力されます。

最後に、PNG→GIFアニメ変換を行います。
ImageMagickのconvertコマンドが使いやすいです。

convert -delay 10 *.png moving.gif