rdkit-data 201403-1 / usr / share / RDKit / Contrib / mmpa / indexing.py

# Copyright (c) 2012, GlaxoSmithKline Research & Development Ltd.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above
#       copyright notice, this list of conditions and the following
#       disclaimer in the documentation and/or other materials provided
#       with the distribution.
#     * Neither the name of GlaxoSmithKline Research & Development Ltd.
#       nor the names of its contributors may be used to endorse or promote
#       products derived from this software without specific prior written
#       permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# Created by Jameed Hussain, September 2012

import sys
import re
from rdkit import Chem
from optparse import OptionParser

def heavy_atom_count(smi):

    m = Chem.MolFromSmiles(smi)
    return m.GetNumAtoms()

def add_to_index(smi,attachments,cmpd_heavy):

    result = False
    core_size = heavy_atom_count(smi) - attachments

    if(use_ratio):
        core_ratio = float(core_size) / float(cmpd_heavy)
        if(core_ratio <= ratio ):
            result = True
    else:
        if(core_size <= max_size):
            result = True

    return result

def get_symmetry_class(smi):

    symmetry = []

    m = Chem.MolFromSmiles(smi)

    #determine the symmetry class
    #see: http://www.mail-archive.com/rdkit-discuss@lists.sourceforge.net/msg01894.html
    #A thanks to Greg (and Alan)
    Chem.AssignStereochemistry(m,cleanIt=True,force=True,flagPossibleStereoCenters=True)

    #get the symmetry class of the attachements points
    #Note: 1st star is the zero index,
    #2nd star is first index, etc
    for atom in m.GetAtoms():
        if(atom.GetMass() == 0):
            symmetry.append(atom.GetProp('_CIPRank'))

    return symmetry

def cansmirk(lhs,rhs,context):

    #cansmirk algorithm
    #1) cansmi the LHS.
    #2) For the LHS the 1st star will have label 1, 2nd star will have label 2 and so on
    #3) Do a symmetry check of lhs and rhs and use that to decide if the labels on
    #   RHS or/and context need to change.
    #4) For the rhs, if you have a choice (ie. two attachement points are symmetrically
    #   equivalent), always put the label with lower numerical value on the earlier
    #   attachement point on the cansmi-ed smiles

    #print "in: %s,%s" % (lhs,rhs)

    isotope_track={}
    #if the star count of lhs/context/rhs is 1, single cut
    stars = lhs.count("*")

    if(stars > 1):
        #get the symmetry class of stars of lhs and rhs
        lhs_sym = get_symmetry_class(lhs)
        rhs_sym = get_symmetry_class(rhs)

    #deal with double cuts
    if(stars == 2):
        #simple cases
        #unsymmetric lhs and unsymmetric rhs
        if( (lhs_sym[0] != lhs_sym[1]) and (rhs_sym[0] != rhs_sym[1]) ):
            #get 1st and 2nd labels and store the new label for it in isotope_track
            #structure: isotope_track[old_label]=new_label (as strings)
            isotope_track = build_track_dictionary(lhs,stars)

            #switch labels using isotope track
            lhs = switch_labels_on_position(lhs)
            rhs = switch_labels(isotope_track,stars,rhs)
            context = switch_labels(isotope_track,stars,context)
                        
        #symmetric lhs and symmetric rhs
        elif( (lhs_sym[0] == lhs_sym[1]) and (rhs_sym[0] == rhs_sym[1]) ):            
            #the points are all equivalent so change labels on lhs and rhs based on position
            #labels on context don't need to change
            lhs = switch_labels_on_position(lhs)
            rhs = switch_labels_on_position(rhs)
            
        #more difficult cases..
        #symmetric lhs and unsymmetric rhs
        elif( (lhs_sym[0] == lhs_sym[1]) and (rhs_sym[0] != rhs_sym[1]) ):            
            #switch labels lhs based on position
            lhs = switch_labels_on_position(lhs)
            #change labels on rhs based on position but need to record
            #the changes as need to appy them to the context
            isotope_track = build_track_dictionary(rhs,stars)
            rhs = switch_labels_on_position(rhs)
            context = switch_labels(isotope_track,stars,context)
            
        #unsymmetric lhs and symmetric rhs
        elif( (lhs_sym[0] != lhs_sym[1]) and (rhs_sym[0] == rhs_sym[1]) ):            
            #change labels on lhs based on position but need to record
            #the changes as need to appy them to the context
            isotope_track = build_track_dictionary(lhs,stars)
            lhs = switch_labels_on_position(lhs)
            context = switch_labels(isotope_track,stars,context)
            #as rhs is symmetric, positions are equivalent so change labels on position
            rhs = switch_labels_on_position(rhs)
            
    #deal with triple cut
    #unwieldy code but most readable I can make it
    elif(stars == 3):
        #simple cases
        #completely symmetric lhs and completely symmetric rhs
        if( ( (lhs_sym[0] == lhs_sym[1]) and (lhs_sym[1] == lhs_sym[2]) and (lhs_sym[0] == lhs_sym[2]) ) and
        ( (rhs_sym[0] == rhs_sym[1]) and (rhs_sym[1] == rhs_sym[2]) and (rhs_sym[0] == rhs_sym[2]) ) ):
            #the points are all equivalent so change labels on lhs and rhs based on position
            #labels on context don't need to change
            lhs = switch_labels_on_position(lhs)
            rhs = switch_labels_on_position(rhs)            
            
        #completely symmetric lhs and completely unsymmetric rhs
        elif( ( (lhs_sym[0] == lhs_sym[1]) and (lhs_sym[1] == lhs_sym[2]) and (lhs_sym[0] == lhs_sym[2]) ) and
        ( (rhs_sym[0] != rhs_sym[1]) and (rhs_sym[1] != rhs_sym[2]) and (rhs_sym[0] != rhs_sym[2]) ) ):
            #alter lhs in usual way
            lhs = switch_labels_on_position(lhs)
            #change labels on rhs based on position but need to record
            #the changes as need to appy them to the context
            isotope_track = build_track_dictionary(rhs,stars)
            rhs = switch_labels_on_position(rhs)
            context = switch_labels(isotope_track,stars,context)
            
        #completely unsymmetric lhs and completely unsymmetric rhs
        elif( ( (lhs_sym[0] != lhs_sym[1]) and (lhs_sym[1] != lhs_sym[2]) and (lhs_sym[0] != lhs_sym[2]) ) and 
        ( (rhs_sym[0] != rhs_sym[1]) and (rhs_sym[1] != rhs_sym[2]) and (rhs_sym[0] != rhs_sym[2]) ) ):
            #build the isotope track
            isotope_track = build_track_dictionary(lhs,stars)
            #alter lhs in usual way
            lhs = switch_labels_on_position(lhs)
            #change rhs and context based on isotope_track
            rhs = switch_labels(isotope_track,stars,rhs)
            context = switch_labels(isotope_track,stars,context)
                        
        #completely unsymmetric lhs and completely symmetric rhs
        elif( ( (lhs_sym[0] != lhs_sym[1]) and (lhs_sym[1] != lhs_sym[2]) and (lhs_sym[0] != lhs_sym[2]) ) and
        ( (rhs_sym[0] == rhs_sym[1]) and (rhs_sym[1] == rhs_sym[2]) and (rhs_sym[0] == rhs_sym[2]) ) ):
            #build isotope trach on lhs
            isotope_track = build_track_dictionary(lhs,stars)
            #alter lhs in usual way
            lhs = switch_labels_on_position(lhs)
            #change labels on context
            context = switch_labels(isotope_track,stars,context)
            #all positions on rhs equivalent so add labels on position
            rhs = switch_labels_on_position(rhs)            

        #more difficult cases, partial symmetry
        #completely unsymmetric on lhs and partial symmetry on rhs
        elif( (lhs_sym[0] != lhs_sym[1]) and (lhs_sym[1] != lhs_sym[2]) and (lhs_sym[0] != lhs_sym[2]) ):            
            #build the isotope track
            isotope_track = build_track_dictionary(lhs,stars)
            #alter lhs in usual way
            lhs = switch_labels_on_position(lhs)
            #change rhs and context based on isotope_track
            rhs = switch_labels(isotope_track,stars,rhs)
            context = switch_labels(isotope_track,stars,context)

            #tweak positions on rhs based on symmetry
            #rhs 1,2 equivalent
            if(rhs_sym[0] == rhs_sym[1]):
                #tweak rhs position 1 and 2 as they are symmetric
                rhs = switch_specific_labels_on_symmetry(rhs,rhs_sym,1,2)

            #rhs 2,3 equivalent
            elif(rhs_sym[1] == rhs_sym[2]):
                #tweak rhs position 1 and 2 as they are symmetric
                rhs = switch_specific_labels_on_symmetry(rhs,rhs_sym,2,3)
                
            #rhs 1,3 equivalent - try for larger set in future
            elif(rhs_sym[0] == rhs_sym[2]):
                #tweak rhs position 1 and 2 as they are symmetric
                rhs = switch_specific_labels_on_symmetry(rhs,rhs_sym,1,3)
                
        #now we are left with things with partial symmetry on lhs and not completely symmetric or unsymmetric on rhs
        else:
            #lhs 1,2,3 equivalent and any sort of partial symmetry on rhs
            if( (lhs_sym[0] == lhs_sym[1]) and (lhs_sym[1] == lhs_sym[2]) and (lhs_sym[0] == lhs_sym[2]) ):
                
                #alter lhs in usual way
                lhs = switch_labels_on_position(lhs)
                #change labels on rhs based on position but need to record
                #the changes as need to appy them to the context
                isotope_track = build_track_dictionary(rhs,stars)
                rhs = switch_labels_on_position(rhs)
                context = switch_labels(isotope_track,stars,context)
                

            #now deal partial symmetry on lhs or rhs.
            #Cases where:
            #lhs 1,2 equivalent
            #lhs 2,3 equivalent
            #lhs 1,3 equivalent
            else:
                #build isotope track on lhs
                isotope_track = build_track_dictionary(lhs,stars)
                #alter lhs in usual way
                lhs = switch_labels_on_position(lhs)
                #change rhs and context based on isotope_track
                rhs = switch_labels(isotope_track,stars,rhs)
                context = switch_labels(isotope_track,stars,context)

                #tweak positions on rhs based on symmetry

                #lhs 1,2 equivalent
                if(lhs_sym[0] == lhs_sym[1]):
                    #tweak rhs position 1 and 2 as they are symmetric on lhs
                    rhs = switch_specific_labels_on_symmetry(rhs,rhs_sym,1,2)                    

                #lhs 2,3 equivalent
                elif(lhs_sym[1] == lhs_sym[2]):
                    #tweak rhs position 1 and 2 as they are symmetric on lhs
                    rhs = switch_specific_labels_on_symmetry(rhs,rhs_sym,2,3)
                    
                #lhs 1,3 equivalent - try for larger set in future
                elif(lhs_sym[0] == lhs_sym[2]):
                    #tweak rhs position 1 and 2 as they are symmetric on lhs
                    rhs = switch_specific_labels_on_symmetry(rhs,rhs_sym,1,3)

    smirk = "%s>>%s" % (lhs,rhs)

    return smirk,context

def switch_specific_labels_on_symmetry(smi,symmetry_class,a,b):

    #check if a and b positions are symmetrically equivalent
    #if equivalent, swap labels if the lower numerical label is not on the
    #1st symmetrically equivalent attachment points in the smi

    if(symmetry_class[a-1] == symmetry_class[b-1]):
        #what are the labels on a and b

        matchObj = re.search( r'\[\*\:([123])\].*\[\*\:([123])\].*\[\*\:([123])\]', smi )
        if matchObj:
            #if the higher label comes first, fix
            if(int(matchObj.group(a)) > int(matchObj.group(b))):
            #if(int(matchObj.group(1)) > int(matchObj.group(2))):
                smi = re.sub(r'\[\*\:'+matchObj.group(a)+'\]', '[*:XX' + matchObj.group(b) + 'XX]' , smi)
                smi = re.sub(r'\[\*\:'+matchObj.group(b)+'\]', '[*:XX' + matchObj.group(a) + 'XX]' , smi)
                smi = re.sub('XX', '' , smi)

    return smi

def switch_labels_on_position(smi):
    
    #move the labels in order of position
    smi = re.sub(r'\[\*\:[123]\]', '[*:XX1XX]' , smi, 1)
    smi = re.sub(r'\[\*\:[123]\]', '[*:XX2XX]' , smi, 1)
    smi = re.sub(r'\[\*\:[123]\]', '[*:XX3XX]' , smi, 1)
    smi = re.sub('XX', '' , smi)
    
    return smi

def switch_labels(track,stars,smi):

    #switch labels based on the input dictionary track
    if(stars > 1):
        #for k in track:
        #        print "old: %s, new: %s" % (k,track[k])

        if(track['1'] != '1'):
            smi = re.sub(r'\[\*\:1\]', '[*:XX' + track['1'] + 'XX]' , smi)

        if(track['2'] != '2'):
            smi = re.sub(r'\[\*\:2\]', '[*:XX' + track['2'] + 'XX]' , smi)

        if(stars == 3):
            if(track['3'] != '3'):
                smi = re.sub(r'\[\*\:3\]', '[*:XX' + track['3'] + 'XX]' , smi)

        #now remove the XX
        smi = re.sub('XX', '' , smi)

    return smi

def build_track_dictionary(smi,stars):

    isotope_track = {}

    #find 1st label, record it in isotope_track as key, with value being the
    #new label based on its position (1st star is 1, 2nd star 2 etc.)
    if(stars ==2):
        matchObj = re.search( r'\[\*\:([123])\].*\[\*\:([123])\]', smi )
        if matchObj:
            isotope_track[matchObj.group(1)] = '1'
            isotope_track[matchObj.group(2)] = '2'

    elif(stars ==3):
        matchObj = re.search( r'\[\*\:([123])\].*\[\*\:([123])\].*\[\*\:([123])\]', smi )
        if matchObj:
            isotope_track[matchObj.group(1)] = '1'
            isotope_track[matchObj.group(2)] = '2'
            isotope_track[matchObj.group(3)] = '3'

    return isotope_track

def index_hydrogen_change():
    #Algorithm details
    #have an index of common fragment(key) => fragments conected to it (values)
    #Need to add *-H to the values where appropriate - and its
    #appropriate when the key is what you would get if you chopped a H off a cmpd.
    #Therefore simply need to check if key with the * replaced with a H is
    #the same as any full smiles in the set
    #
    #Specific details:
    #1) Loop through keys of index
    #2) If key is the result of a single cut (so contains only 1 *) replace the * with H, and cansmi
    #3) If full smiles matches key in hash above, add *-H to that fragment index.

    for key in index:

        attachments = key.count('*')
        #print attachments

        if(attachments==1):

            smi = key

            #simple method
            smi = re.sub(r'\[\*\:1\]', '[H]' , smi)
            
            #now cansmi it
            temp = Chem.MolFromSmiles(smi)

            if(temp == None):
                sys.stderr.write('Error with key: %s, Added H: %s\n' %(key,smi) )
            else:
                c_smi = Chem.MolToSmiles( temp, isomericSmiles=True )

                if(c_smi in smi_to_id):
                    core = "[*:1][H]"
                    id = smi_to_id[c_smi]

                    value = "%s;t%s" % (id,core)
                    #add to index
                    index[key].append(value)
                

if __name__=='__main__':

    #note max heavy atom count does not
    #include the attachement points (*)
    max_size = 10
    ratio = 0.3
    use_ratio = False

    index={}
    smi_to_id={}
    id_to_smi={}

    id_to_heavy={}

    #set up the command line options
    #parser = OptionParser()
    parser = OptionParser(description="Program to generate MMPs")
    parser.add_option('-s', '--symmetric', default=False, action='store_true', dest='sym',
                      help='Output symmetrically equivalent MMPs, i.e output both cmpd1,cmpd2, SMIRKS:A>>B and cmpd2,cmpd1, SMIRKS:B>>A')
    parser.add_option('-m','--maxsize',action='store', dest='maxsize', type='int',
                      help='Maximum size of change (in heavy atoms) allowed in matched molecular pairs identified. DEFAULT=10. \
                      Note: This option overrides the ratio option if both are specified.')
    parser.add_option('-r','--ratio',action='store', dest='ratio', type='float',
                      help='Maximum ratio of change allowed in matched molecular pairs identified. The ratio is: size of change / \
                      size of cmpd (in terms of heavy atoms). DEFAULT=0.3. Note: If this option is used with the maxsize option, the maxsize option will be used.')

    #parse the command line options
    (options, args) = parser.parse_args()

    #print options
    if(options.maxsize != None):
        max_size = options.maxsize
    elif(options.ratio != None):
        ratio = options.ratio
        if(ratio >= 1):
            print "Ratio specified: %s. Ratio needs to be less than 1."
            sys.exit(1)
        use_ratio = True

    #read the STDIN
    for line in sys.stdin:

        line = line.rstrip()
        smi,id,core,context = line.split(',')

        #fill in dictionaries
        smi_to_id[smi]=id
        id_to_smi[id]=smi

        #if using the ratio option, check if heavy atom
        #of mol already calculated. If not, calculate and store
        cmpd_heavy = None
        if(use_ratio):
            if( (id in id_to_heavy) == False):
                id_to_heavy[id] = heavy_atom_count(smi)

            cmpd_heavy = id_to_heavy[id]

        #deal with cmpds that have not been fragmented
        if(len(core) == 0) and (len(context) == 0):
            continue

        #deal with single cuts
        if(len(core) == 0):
            side_chains = context.split('.')

            #minus 1 for the attachement pt
            if( add_to_index(side_chains[1],1,cmpd_heavy)==True ):
                context = side_chains[0]
                core = side_chains[1]

                value = "%s;t%s" % (id,core)

                #add the array if no key exists
                #add the context with id to index
                index.setdefault(context, []).append(value)

            #minus 1 for the attachement pt
            if( add_to_index(side_chains[0],1,cmpd_heavy)==True ):
                context = side_chains[1]
                core = side_chains[0]

                value = "%s;t%s" % (id,core)

                #add the array if no key exists
                #add the context with id to index
                index.setdefault(context, []).append(value)

        #double or triple cut
        else:

            attachments = core.count('*')

            if( add_to_index(core,attachments,cmpd_heavy)==True ):
                value = "%s;t%s" % (id,core)

                #add the array if no key exists
                #add the context with id to index
                index.setdefault(context, []).append(value)

    #index the H change
    index_hydrogen_change()

    #Now index is ready

    #loop through the index
    for key in index:

        total = len(index[key])

        #check if have more than one value
        if(total == 1):
            continue

        for xa in xrange(total):

            for xb in xrange(xa, total):

                if(xa != xb):
                    #now generate the pairs

                    id_a,core_a = index[key][xa].split(";t")
                    id_b,core_b = index[key][xb].split(";t")

                    #make sure pairs are not same molecule
                    if(id_a != id_b):

                        #make sure LHS and RHS of SMIRKS are not the same
                        if(core_a != core_b):

                            smirks,context = cansmirk(core_a,core_b,key)
                            print "%s,%s,%s,%s,%s,%s" % ( id_to_smi[id_a], id_to_smi[id_b], id_a, id_b, smirks, context )

                            #deal with symmetry switch
                            if(options.sym == True):
                                smirks,context = cansmirk(core_b,core_a,key)
                                print "%s,%s,%s,%s,%s,%s" % ( id_to_smi[id_b], id_to_smi[id_a], id_b, id_a, smirks, context )