RDKitで2d_pharmacophore

RDKitにある2Dファーマコフォアを使って
シミラリティマトリックスを作るスクリプトを書いてみた。
 下のデフォルト設定では3000ビット以上のFPになるので
適当なbinに指定し直した方が計算が速くなるかも。
分子によってはGitBitIdxでエラーになるのでそれを回避するため
try/except処理にした。エラーが出る分子については後で原因を追う。
TCの計算はお決まりの書き方でやると、エラーにうまく対処できない。
 見にくいけどdict型でindxとFPの関係を保つようにしてみた。
よくよく考えると、SVMとかするならFPをbitstringで返しておしまいにした方がいいと気がついた。
使い方
$ python p4core.py で入力SDFに関してシミラリティマトリックスを出力される。
ECFP/FCFP/MACCSのTC
諸々と合わせて出力させると違いが分かるかな。

#p4core.py
from rdkit import Chem
from rdkit.Chem import ChemicalFeatures # 2D PharmacoPhore FingerPrints
from rdkit.Chem import Pharm2D
from rdkit.Chem.Pharm2D import Generate
from rdkit.Chem.Pharm2D import SigFactory
from rdkit.Chem.Pharm2D import DefaultSigFactory
from rdkit import RDConfig
from rdkit import DataStructs
import os, sys

dataDir  = os.path.join(RDConfig.RDCodeDir,'Chem/Pharm2D/test_data')
featFact = ChemicalFeatures.BuildFeatureFactory(os.path.join(dataDir,'BaseFeatures.fdef'))
print "########"
sigFact = SigFactory.SigFactory( featFact )
bins = tuple([(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,100)])
sigFact.SetBins(bins)
sigFact.Init()


def calc_p4_fp( mol ):
    fp = Generate.Gen2DFingerprint( mol, sigFact )
    return fp

def calc_tanimoto ( fp1, fp2 ):
    tc = DataStructs.TanimotoSimilarity( fp1,fp2 )
    return tc

if __name__ == "__main__":
    if len(sys.argv) != 3:
        print "You need 2 arguments, input.sdf and output.txt"

    mols = [ mol for mol in Chem.SDMolSupplier(sys.argv[1]) ]

    print str(len(mols)) + "read."

    output = open( sys.argv[2], "w" )
    output.write( "mol1\tmol2\tp4core_tc" )
    output.write( "\n" )

    fps = {}
    error = {}
    for i,mol in enumerate(mols):
        try:
            fp = calc_p4_fp( mol )
            fps[i] = fp
        except:
            error[i] = Chem.MolToSmiles( mol )

    for i in range( len(fps) ):
        for j in range( i ):
            mol1_smi = Chem.MolToSmiles( mols[fps.keys()[i]] )
            mol2_smi = Chem.MolToSmiles( mols[fps.keys()[j]] )
            p4_tanimoto = calc_tanimoto( fps[fps.keys()[i]], fps[fps.keys()[j]] )
            output.write( "%s\t%s\t%s"%( mol1_smi, mol2_smi, p4_tanimoto ) )
            output.write( "\n" )
    output.close()

    error_log = open("errormol.txt","w")
    for i,smi in error.items():
        error_log.write("%s\t%s\n"%(i, smi))
    error_log.close()
    
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s