Canonical label その2
前回のつづき。
原子のラベル付けを何回か繰り返してcanonical labelを生成する
というプロセスについて書きかけたので、続けてみます。
ペンタン分子（C-C-C-C-C）の場合、まず原子の状態を数値化して、
10106003--20206002--20206002--20206002--10106003
となります。前回のとおりです。ここで、値の小さいものから順位をつけると、
1--2--2--2--1
となります。しかし、同じ2番のついている炭素原子でも、
中心の炭素から見れば両脇にエチル基(-CC)が位置しているのに対し、
その隣の炭素から見ればメチル基(-C)とプロピル基(-CCC)が隣接します。
この違いを順位に反映させるために、「Extended Connectivity」と呼ばれる
考え方を以下で導入します。
基本的な方針は、隣接している原子の順位を数値化することです。
最もシンプルなものは、隣接原子の順位を単純に足し合わせるものです。
ペンタンの場合、
(2)--(1+2)--(2+2)--(2+1)--(2)
これに再び順位を付けると、
1--2--3--2--1
となり、中心の炭素を差別化することができました。
ただし、足し算の場合、和が衝突するという欠点があります。
たとえば、隣接原子の順位が(5, 5, 2)である原子と
(3, 3, 6)である原子を比較した場合、
5+5+2 = 12
3+3+6 = 12
となり、隣接要素が違うにもかかわらず同順位になってしまいます。
二乗和としても、
25+25+4 = 54
9+9+36 = 54
となり、やはり衝突することがあります。
そこで「素数」の登場です。
素数を小さい順に並べた数列 (2, 3, 5, 7, 11, 13, 17, ...) を準備し、
先ほどの順位列をこの数列に置き換えます。
そして、和ではなく、積をとってみましょう。すると、
11x11x3 = 363
5x5x13 = 325
となり、区別できます。素数の特性を活かしたエレガントな方法ですね。
このプロセスを繰り返せば、２つ隣、３つ隣の原子の順位も伝播させることができます。
順位の入れ替わりがなくなった段階で、最終的なcanonical labelが決定します。
詳細は、Weiningerらの論文を参考に*1。次回は実装です。
*1: Weininger et al. SMILES. 2. Algorithm for Generation of Unique SMILES Notation. J. Chem. Inf. Comput. Sci. (1989) 29, 97-101
Canonical label その1
前回のつづき。
Canonical SMILES文字列を生成するためには、
まず、文字列の開始点となる原子を決定しなければなりません。
そのために、それぞれの原子がどのような環境にいるかを数値化して、
その順番に原子を並べ替える、というのがひとつの方法です。
Weiningerら*1は、以下の6種類の基準で数値化し、それらの数値を連結して
8桁の整数を生成して昇順に並べました。
こうすることで、たとえば、エーテル分子（C-C-O-C-C）の場合、
10106003--20206002--20208000--20206002--10106003
という数値列に置き換えられるので、
最小値の"10106003"に対応する原子から始めれば良いということになります。
ただし、最小値が2ヶ所以上にある場合は、
Morgan法のように、隣接原子の順位を反映させつつ並び替えを繰り返します。
（この話は長くなりそうなので、今回は一旦ここで止めておきます）
さて、CDKでは、CanonicalLabelerクラスとして、一連のラベル付け機能が実装されています。
Weiningerらの表現が曖昧なので、実装によって結果も変わってきます。
CDKとOpenBabelの出力結果の違いも、この辺に原因がありそうです。
以下は、CDKを土台としたサンプルコードですが、少し変更を加えています。
/* inv_label.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 org.openscience.cdk.tools.manipulator.AtomContainerManipulator;import java.util.*;
class inv_label {
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");StringBuffer inv;
ArrayList vect = new ArrayList();
IAtomContainer ac;try {
AtomContainerManipulator acm = new AtomContainerManipulator();
ac = acm.removeHydrogens(imol);
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
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));
}
}catch (Exception e1) {
e1.printStackTrace();
}
for (int i = 0 ; i < vect.size() ; i++){
System.out.printf("%s\t%d\t%s\n", id, i, vect.get(i));
}
}
}
}
*1: Weininger et al. SMILES. 2. Algorithm for Generation of Unique SMILES Notation. J. Chem. Inf. Comput. Sci. (1989) 29, 97-101
Canonical SMILES
ふたつの化学構造が同一構造なのか、手っ取り早く知りたいとき、
"Canonical (正準化された)" SMILES形式を使うと便利です。
化学構造を"通常の"SMILES形式で表現しようとすると、一意に定まらないことがあります。
たとえば、エーテル（diethyl ether）の場合、
CCOCC
C(C)OCC
C(OCC)C
O(CC)CC
のように、記述の開始点と順序により、４通りの表現が可能です。
そこで、これらの中からひとつ代表SMILESに統一しよう、そうすれば
文字列の一致が化学構造の一致に対応付けられる、というのが"Canonical"の考え方です。
さて、ここからは実装の話で、私も最近知ったことなのですが、
あるソフトで"Canonical"であったとしても、別のソフトでそうとは限りません。
ソフトによって、それぞれ"Canonical"の定義が存在するようなのです。
（これに気づいたのは、OpenBabelのver.2.2.0 と ver.2.3.0 の間の相違でした）
たとえば、OpenBabelとCDKの間でも、出力結果が異なります。
CDKの場合、以下のスクリプトで、SDF形式からCanonical SMILES形式に変換できます。
/* sdf2can.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;class sdf2can {
public static void main(String args[]){
if(args.length!=1){
System.err.println("sdf2can <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");
SmilesGenerator sg = new SmilesGenerator();
AtomContainerManipulator acm = new AtomContainerManipulator();
try {
sg.setUseAromaticityFlag(true);
String can = sg.createSMILES(acm.removeHydrogens(imol));
System.out.printf("%s\t%s\n", can, id);
} catch (Exception e1) {
e1.printStackTrace();
}
}
}
}
次回では、CDKにおける"Canonical"の定義について、追いかけてみようと思います。
BEDROC
CRANのenrichvsパッケージにメールアドレスを載せていた影響なのか、
先日質問メールがとどきました。どうやら、Truchonらの論文*1を読んで
ExcelでBEDROCを計算しようとして躓いているご様子。
従来より、スクリーニングの性能評価には、
エンリッチメント曲線やROC曲線が用いられてきました。
しかし、もっと曲線の立ち上がり（すなわち、スコア上位化合物の活性／非活性）に
焦点を当てた尺度が必要、とのことで提案されたのが、BEDROCです。
通常、ROC曲線下面積を求める時には、そのままTP(真陽性) / FP(偽陽性)を足しあわせていきますが、
BEDROCの場合、FPが大きくなるにつれ、TPの重みが小さくなっていきます。
enrichvsにもBEDROC関数を実装していたので、質問メールが届いたのでした。
いちおう返信しましたが、大丈夫ですかな。
CRANには、pROCも実装されていますが、これも同じような考え方ですね。
横軸が対数スケールになるので、その曲線下面積は、FPが小さい範囲での
予測性能に左右されるように設計されています。
ファーマコフォア検索
ファーマコフォア (pharmacophore)とは、分子構造の抽象的な特徴を３次元空間上に配置したものです。
たとえば、２つの窒素原子と１つの酸素原子が以下のように配置した
ファーマコフォアを定義してみます。
CDKのPharmacophoreMatcherクラスを使って、この条件にあう分子を検索できるようです。
以下、サンプルコードです。
（編集中）
/* ph4_search.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.isomorphism.matchers.QueryAtomContainer;
import org.openscience.cdk.pharmacophore.*;
class ph4_search {
public static void main(String args[]){
if(args.length!=1){
System.err.println("ph4 <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();
}PharmacophoreQuery query = new PharmacophoreQuery();
PharmacophoreQueryAtom o = new PharmacophoreQueryAtom("D", "[OX1]");
PharmacophoreQueryAtom n1 = new PharmacophoreQueryAtom("A", "[N]");
PharmacophoreQueryAtom n2 = new PharmacophoreQueryAtom("A", "[N]");
query.addAtom(o);
query.addAtom(n1);
query.addAtom(n2);
PharmacophoreQueryBond b1 = new PharmacophoreQueryBond(o, n1, 4.0, 4.5);
PharmacophoreQueryBond b2 = new PharmacophoreQueryBond(o, n2, 4.0, 5.0);
PharmacophoreQueryBond b3 = new PharmacophoreQueryBond(n1, n2, 5.4, 5.8);
query.addBond(b1);
query.addBond(b2);
query.addBond(b3);while( isr.hasNext() ){
IMolecule imol = (IMolecule)isr.next();
String id = (String)imol.getProperty("cdk:Title");
PharmacophoreMatcher matcher = new PharmacophoreMatcher();
matcher.setPharmacophoreQuery(query);
try{
if( matcher.matches(imol, true) ){
System.out.printf("%s\n", id);
}
} catch (Exception e1) {
e1.printStackTrace();
}
}
}
}
これを応用すれば、
何千種類ものファーマコフォア・パターンを予め定義しておいて、
各分子がそれぞれのパターンに適合するかを{0, 1}のベクトルで表現する
「ファーマコフォア・フィンガープリント」も実装可能です。
（注：PharmacophoreQueryAtomのSMARTSパターンに制約があるようなので、
市販ソフトほどの使い勝手はありません）
Mendeleyで文献情報共有
最近は、論文誌の種類が増えて、
先行文献を網羅的に把握することが難しくなってきています。
また、過去に読んだ論文でも、
あの文章はどの論文のくだりだったか、思い出すだけでもひと苦労です。
このブログで紹介した文献も10種類を超えましたので、これを機に
Mendeleyというフリーの文献管理ソフトを導入してみます。
Mendeleyの特徴は、iPaperのように、PDFをドラッグ＆ドロップで論文を登録できる点です。
また、横断的にPDF全文検索できる点もありがたいですね。
なお、容量に制限（500MB）がありますが、
Webサイトを介して、他の研究者と論文情報を共有することもできます。
Webサイト側の論文リストは、１クリックでローカルPC側の論文リストと同期します。
作ってみました。
・けみげの専用、文献リスト
http://www.mendeley.com/groups/1295393/%E3%81%91%E3%81%BF%E3%81%92%E3%81%AE/papers/
（こんな使い方して良いのかな・・・少なくとも英語で記述すべきだったかも）
分子骨格と作用選択性
さきの記事で、「分子骨格(molecular framework)」を定義しました。
今回はその応用編です。
Yangらは、標的タンパク質選択性（最近流行のキーワードでもあります）との関連性に注目しました*1。
- f_MF = (「分子骨格」内の重原子数) / (全体の重原子数)
と定義するときに、f_MF = 0.65を超えると比較的多くのタンパク質に作用する、と言うのです。
確かに、骨格以外のニョキニョキした構造が
ちょうど鍵のギザギザしたイメージと重なり、
特定のタンパク質にだけピッタリ嵌りそうな気もします。
当たり前のようで誰も計算してみなかった統計。
J. Med. Chem.誌に載るためには、
そのようなコロンブスの卵的な要素が重要になってきそうです。
さて、f_MFを算出するプログラムは以下のとおりです。
環構造を１つしか持たないときは、さらに条件分岐して計算します。
/* fmf.java */import java.io.*;
import org.openscience.cdk.interfaces.*;
import org.openscience.cdk.io.iterator.IteratingSMILESReader;
import org.openscience.cdk.interfaces.IMolecule;
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import org.openscience.cdk.fragment.MurckoFragmenter;class fmf {
public static void main(String args[]){
if(args.length!=1){
System.err.println("fmf <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");MurckoFragmenter fragmenter = new MurckoFragmenter(true, 3);
double result = 0.0;
try {
fragmenter.generateFragments(imol);
IAtomContainer[] framework = fragmenter.getFrameworksAsContainers();
IAtomContainer[] ringSystems = fragmenter.getRingSystemsAsContainers();
if (framework.length == 1) {
result = framework[0].getAtomCount() /
(double) imol.getAtomCount();
} else if (framework.length == 0 && ringSystems.length == 1) {
result = ringSystems[0].getAtomCount() /
(double) imol.getAtomCount();
} else result = 0.0;
} catch (Exception e1) {
e1.printStackTrace();
}
System.out.printf("%s\t%.3f\n", id, result);
}
}
}
CDKでは、FMFDescriptorクラスとしてまとめられています。
