2011-12-10

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

PyMOL

ひょんなことから、
「科学・技術フェスタ in 京都2011」の出展資料を準備することになりました。

研究室のボスによると、
「くすりがタンパク質に結合する様子を動画で示したい」
とのこと。いま、PyMOLという立体構造表示ソフトで動画を作成しているところです。

簡単に足あとを残しておきます。

PDBで結晶構造を探す( http://www.rcsb.org/pdb/explore.do?structureId=2ITY )
- ページ右上の「Download Files」から、タンパク質のPDBファイルをダウンロード
- ページ中ほどの「Ligand Chemical Component」から、リガンドのSDファイルをダウンロード
リガンドのxyz座標を少しずつ移動させながら、ひとつのSDファイルに出力
PyMOLでPDBファイルと上記SDファイルを開き、再生ボタンをクリック
- 化学構造を順番に表示する → 動いているように見える

という流れです。
落ち着いたときに、プログラムや動画も併せて更新しようと思います。

2011-11-04

部分構造検索 (2)

OpenBabelに「obgrep」という部分構造検索コマンドがあるのですが、
CDKのUniversalIsomorphismTesterクラスとどちらが優れているか検証してみようという企画です。

【比較条件】

QueryはSMILES形式の環構造。
検索対象DBは、DrugBankの6620化合物（SDF形式）
CDK ver1.4.2
OpenBabel ver2.3.0

CASE1: βラクタム [N1C(=O)CC1]
まずは小さめの環から。
無事に同じ結果（76 件）が出力されました。以降もそうあって欲しいのですけれども。
計算速度は、CDK: 12秒, OpenBabel: 15秒でした。

CASE2: 1,4-ジヒドロピリジン [C1C=CNC=C1]
これも単環なので大丈夫でしょう、と思いきや
・CDK: 35 件
・obgrep: 14 件
という異なる出力。
分子を見比べてみると、どうやらCDKでは、ニューキノロン系抗菌薬もヒットするようです。
確かにジヒドロピリジンを部分構造に含んでいます。
二重結合が芳香性に寄与しているという点が実装の分かれ目ですね。

CASE3: クマリン [c1ccc2c(c1)ccc(=O)o2]
共鳴するケトンがいやらしい感じですが、
・CDK: 0 件
・obgrep: 28 件
ここでは、CDKがダメですね。原因は不明です。

このシリーズは、もう少し続けてみます。

2011-10-09

CPSA記述子

（編集中）

2011-10-01

ファンデルワールス体積

CDK descriptor

わざわざ3次元座標を生成しなくても、
各原子の原子タイプから回帰分析により近似できる記述子があります。
たとえば、疎水性（XLogP）、表面積（TPSA）、そして今回、
分子の体積を求めます。

Zhaoら*1によると、
以下の式でファンデルワールス体積Vが近似できるそうです。
　V = Σ(原子の寄与) - 5.92 x (結合の数) - 14.7 x (芳香環の数) - 3.8 x (非芳香環の数)

CDKにも、VABCDescriptorクラスとして実装されています。



/* vabc.java */
import java.io.*;

import org.openscience.cdk.interfaces.*;

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

import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import java.util.HashMap;

import java.util.Map;

import org.openscience.cdk.CDKConstants;

import org.openscience.cdk.aromaticity.CDKHueckelAromaticityDetector;

import org.openscience.cdk.config.AtomTypeFactory;

import org.openscience.cdk.exception.CDKException;

import org.openscience.cdk.ringsearch.SSSRFinder;


class vabc {
    private static Map<String,Double> bondiiVolumes = new HashMap<String, Double>() {{

        put("H", 7.2382293504);

        put("C", 20.5795259250667);

        put("N", 15.5985308577667);

        put("O", 14.7102267005611);

        put("Cl", 22.4492971208333);

        put("Br", 26.5218483279667);

        put("F", 13.3057882007064);

        put("I", 32.5150310206656);

        put("S",24.4290240576);

        put("P",24.4290240576);

        put("As", 26.5218483279667);

        put("B", 40.48);

        put("Se", 28.7309115245333);

        put("Si", 38.7923854248);

    }};

    private static AtomTypeFactory atomTypeList = null;
    public static void main(String args[]){
        if(args.length!=1){

            System.err.println("vabc <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");
            double sum = 0.0;

            try{

                if (atomTypeList == null) {

                    atomTypeList = AtomTypeFactory.getInstance(

                             "org/openscience/cdk/dict/data/cdk-atom-types.owl",

                             imol.getBuilder() // take whatever we got first

                    );

                }

                int totalHCount = 0;

                for (IAtom atom : imol.atoms()) {

                    Double bondiiVolume = bondiiVolumes.get(atom.getSymbol());

                    if (bondiiVolume == null)

                        throw new CDKException("Unsupported element.");
                    sum += bondiiVolume;
                    // add volumes of implicit hydrogens?

                    IAtomType type =  atomTypeList.getAtomType(atom.getAtomTypeName());

                    if (type == null)

                        throw new CDKException("Unknown atom type: " + atom.getSymbol());

                    if (type.getFormalNeighbourCount() == null)

                        throw new CDKException("Formal neighbor count not given for : "

                                               + type.getAtomTypeName());

                    int hCount = type.getFormalNeighbourCount() -

                        imol.getConnectedAtomsCount(atom);

                    sum += (hCount * bondiiVolumes.get("H"));

                    totalHCount += hCount;

                }

                sum -= 5.92 * (imol.getBondCount() + totalHCount);
                boolean[] originalFlags = imol.getFlags();

                CDKHueckelAromaticityDetector.detectAromaticity(imol);

                SSSRFinder ringFinder = new SSSRFinder(imol);

                IRingSet ringSet = ringFinder.findSSSR();

                if (ringSet.getAtomContainerCount() > 0) {

                    int aromRingCount = 0;

                    int nonAromRingCount = 0;

                    for (IAtomContainer ring : ringSet.atomContainers()) {

                        if (ringIsAromatic(ring)) {

                            aromRingCount++;

                        } else {

                            nonAromRingCount++;

                        }

                    }

                    imol.setFlags(originalFlags);

                    sum -= 14.7 * aromRingCount;

                    sum -= 3.8 * nonAromRingCount;

                }

            } catch (Exception e1) {

                e1.printStackTrace();

            }

            System.out.printf("%s\t%.3f\n", id, sum);

        }

    }


    private static boolean ringIsAromatic(IAtomContainer ring) {

        for (IAtom atom : ring.atoms()) {

            if (!atom.getFlag(CDKConstants.ISAROMATIC)) return false;

        }

        return true;

    }
}

*1: Zhao et al. Fast Calculation of van der Waals Volume as a Sum of Atomic and Bond Contributions and Its Application to Drug Compounds. J. Org. Chem. (2003) 68, 7368-7373 PubMed

2011-09-23

E-state指標

CDK descriptor

先週までのcanonicalな話は、
ECFPというフィンガープリントを実装するための準備だったのですが、
長引きそうなので、少し休題。

さて、気をとり直して次の話題に。
電子の分布を数値で表す「E-state指標(electrotopological index)」*1です。

まず、有機化学の基礎から。
共有結合は、σ結合とπ結合から構成されています。
π結合はσ結合より弱いため、
π結合上の電子（π電子）は、孤立電子対のように、
求電子剤と反応しやすい性質を持ちます。
では、そのような反応性電子の存在状況を原子ごとに数値化してみます。

ある原子をi、その価電子の数をδ_i、σ電子の数をσ_i、主量子数をn_i としたときに、
　I_i = ( (2/L_i)^2 x δ_i + 1 ) / σ_i
を「固有状態(intrinsic state)」と呼びます。
σ電子の数に比べてπ電子や孤立電子対の数が多い原子ほど、
すなわち電子豊富な原子ほど、I_iが大きくなります。
たとえば、メチル基(-CH3)の炭素では、I_i = 2 なのに対して、
カルボニル酸素(=O)では、I_i = 7 となります。

そして、原子間の電子による影響ΔI_iを考慮して、
E-state指標 S_i を以下のように定義します。
　S_i = I_i + ΔI_i = I_i + Σ{ (I_i - I_j) / (d_ij + 1)^2 }
ここで、d_ij は、原子i,j間の距離（経路長）です。
Σは、原子jについて足しあわせます。

分子単位で数値化するときには、
すべての原子のE-state指標 S_i を足しあわせます。
これは固有状態和（intrinsic state sum）と呼ばれているようです。

CDKでは、固有状態の計算までは実装されておらず、
電子の状況で分類された原子タイプの頻度を数えるKierHallSmartsDescriptorクラスが
あるのみです。その名の通り、SMARTSパターンで検索しています。
E-stateの実装は難しくないのですが、
単に頻度を出力する方が実用的なのかもしれませんね。

*1: Hall and Kier, Electrotopological State Indices for Atom Types: A Novel Combination of Electronic, Topological, and Valence State Information. J. Chem. Inf. Comput. Sci. (1995) 35, 1039-1045

2011-09-17

Canonical label その3

CDK

各原子に割り当てられた素数の積と
並び替えらた原子の中間順位を出力してみます。



/* inv_label_EC.java */
import java.io.*;

import org.openscience.cdk.interfaces.*;

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

import org.openscience.cdk.DefaultChemObjectBuilder;
import org.openscience.cdk.smiles.SmilesGenerator;

import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import org.openscience.cdk.smiles.InvPair;

import org.openscience.cdk.tools.periodictable.PeriodicTable;

import org.openscience.cdk.CDKConstants;

import java.util.*;


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

            System.err.println("inv_label <sd-file>");

            System.exit(1);

        }

        FileInputStream fis = null;

        IteratingMDLReader isr = null;

        try{

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

            isr = new IteratingMDLReader(fis, DefaultChemObjectBuilder.getInstance() );

        } catch (Exception e) {

            e.printStackTrace();

        }
        while( isr.hasNext() ){

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

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

            try {

                AtomContainerManipulator acm = new AtomContainerManipulator();

                IAtomContainer atoms = acm.removeHydrogens(imol);

                vect = createInvarLabel(atoms);
                System.out.printf("%s", id);

                for (int i = 0 ; i < vect.size() ; i++){

                    System.out.printf(" %s", atoms.getAtom(i).getSymbol() );

                }

                System.out.print("\n");

                for (int i = 0 ; i < vect.size() ; i++){

                    System.out.printf(" %s",

                                      atoms.getAtom(i).getProperty("InvariancePair"));

                }

                System.out.print("\n");
                step3(vect, atoms);
            }catch (Exception e1) {

                e1.printStackTrace();

            }

        }

    }
    private static void primeProduct(ArrayList v, IAtomContainer atomContainer) {

        Iterator it = v.iterator();

        Iterator n;

        InvPair inv;

        IAtom a;

        long summ;

        while (it.hasNext()) {

            inv = (InvPair) it.next();

            List neighbour = atomContainer.getConnectedAtomsList(inv.getAtom());

            n = neighbour.iterator();

            summ = 1;

            while (n.hasNext()) {

                a = (IAtom) n.next();

                int next = ( (InvPair)a.getProperty(InvPair.INVARIANCE_PAIR)).getPrime();

                summ = summ * next;

            }

            inv.setLast(inv.getCurr());

            inv.setCurr(summ);

        }

    }
    private static void sortArrayList(ArrayList v) {

        Collections.sort(v, new Comparator() {

                public int compare(Object o1, Object o2) {

                    return (int) ( ( (InvPair) o1).getCurr() - ( (InvPair)o2).getCurr());

                }

            });

        Collections.sort(v, new Comparator() {

                public int compare(Object o1, Object o2) {

                    return (int) ( ( (InvPair) o1).getLast() - ( (InvPair)o2).getLast());

                }

            });

    }
    private static void rankArrayList(ArrayList v) {

        int num = 1;

        int[] temp = new int[v.size()];

        InvPair last = (InvPair) v.get(0);

        Iterator it = v.iterator();

        InvPair curr;

        for (int x = 0; it.hasNext(); x++) {

            curr = (InvPair) it.next();

            if (!last.equals(curr)) {

                num++;

            }

            temp[x] = num;

            last = curr;

        }

        it = v.iterator();

        for (int x = 0; it.hasNext(); x++) {

            curr = (InvPair) it.next();

            curr.setCurr(temp[x]);

            curr.setPrime();

        }

    }
    private static boolean isInvPart(ArrayList v) {

        if ( ( (InvPair) v.get(v.size()-1)).getCurr() == v.size())

            return true;

        Iterator it = v.iterator();

        InvPair curr;

        while (it.hasNext()) {

            curr = (InvPair) it.next();

            if (curr.getCurr() != curr.getLast())

                return false;

        }

        return true;

    }
    private static void breakTies(ArrayList v) {

        Iterator it = v.iterator();

        InvPair curr;

        InvPair last = null;

        int tie = 0;

        boolean found = false;

        for (int x = 0; it.hasNext(); x++) {

            curr = (InvPair) it.next();

            curr.setCurr(curr.getCurr() * 2);

            curr.setPrime();

            if (x != 0 && !found && curr.getCurr() == last.getCurr()) {

                tie = x - 1;

                found = true;

            }

            last = curr;

        }

        curr = (InvPair) v.get(tie);

        curr.setCurr(curr.getCurr() - 1);

        curr.setPrime();

    }
    private static ArrayList createInvarLabel(IAtomContainer ac) {

        StringBuffer inv;

        ArrayList vect = new ArrayList();

        for (int i = 0 ; i < ac.getAtomCount(); i++){

            IAtom a = ac.getAtom(i);

            inv = new StringBuffer();

            Integer ddd = a.getImplicitHydrogenCount();
            // 1. Num connections (incl H)

            inv.append( ac.getConnectedAtomsList(a).size() +

                        (a.getImplicitHydrogenCount() == CDKConstants.UNSET ?

                         0 : a.getImplicitHydrogenCount()) );

            // 2. Num of non H bonds

            //System.out.printf(" %d\n", ac.getConnectedAtomsList(a).size() );

            inv.append( String.format("%02d",

                                      ac.getConnectedAtomsList(a).size()) );

            // 3. Atomic num

            inv.append( String.format("%02d",

                                      PeriodicTable.getAtomicNumber(a.getSymbol()) ) );

            // 4. Sign of charge

            Integer fcharge = a.getFormalCharge();

            if (fcharge == CDKConstants.UNSET) fcharge = 0;

            if (fcharge < 0){

                inv.append(1);

            } else {

                inv.append(0);

            }

            // 5. Absolute charge

            inv.append( (int)Math.abs( (a.getFormalCharge() ==

                                CDKConstants.UNSET ? 0.0 : a.getFormalCharge()) ) );

            // 6. Hydrogen count

            inv.append( (a.getImplicitHydrogenCount() == CDKConstants.UNSET ?

                         0 : a.getImplicitHydrogenCount()));
            vect.add(new InvPair(Long.parseLong(inv.toString()), a));

        }

        return vect;

    }

}

けみげの