select diverse molecules from dataset.

To select subset from compound library, I think about molecular diversity.
So, I often use clustering .
There are many clustering method in chemo-informatics area.
But, today I used another approach. Use GA.
My strategy is that, evaluate diversity of molecules is sum of dissimilarity.
Let’s start.
For convenience, I used rdkit first_prop_200.sdf as dataset.

import scipy.misc as scm
import numpy as np
import random
from rdkit import Chem
from rdkit.Chem import DataStructs
from rdkit.Chem import AllChem
from rdkit.Chem import DataStructs
from deap import creator, base
from deap import algorithms
from deap import tools

mols = [mol for mol in Chem.SDMolSupplier('first_200.props.sdf')]
n_mols = len( mols )

# define diversity as sum of dissimilarity
def eval_diversity(individual):
    combination = scm.comb( len(individual),2 )
    dis_sim_mat=[]
    fps =[ AllChem.GetMorganFingerprintAsBitVect(mols[i], 2) for i in individual ]
    for idx in range(1,len(fps)):
        # return distance returnDistance=1
        dis_sim = DataStructs.BulkCosineSimilarity(fps[idx], fps[:idx], returnDistance=1)
        dis_sim_mat.extend(dis_sim)
    return sum(dis_sim_mat)/combination,

Points.
1. I want to maximise dissimilarity
2. each chromosome has ten molecules.
3. population is 20.

creator.create( 'Fitness', base.Fitness, weights=(1.0,) )
creator.create( 'Individual', list, fitness = creator.Fitness )
toolbox = base.Toolbox()
#generate random int from range(200). It means select molecules in random.
toolbox.register('attr_item', random.randint, 0,n_mols)
# select 10 molecules from library and maximize diversity.
toolbox.register('individual', tools.initRepeat, creator.Individual, toolbox.attr_item, 10)
toolbox.register('population', tools.initRepeat, list, toolbox.individual)

toolbox.register('evaluate', eval_diversity)
toolbox.register('mate', tools.cxTwoPoint)
toolbox.register('mutate', tools.mutShuffleIndexes, indpb = 0.05)
toolbox.register('select', tools.selTournament, tournsize=3)

Then write main function.

def main():
    random.seed(64)
    #set cross over and mutation probability and generation.
    CXPB, MUTPB, NGEN= 0.5, 0.2, 40
    # use all molecures.
    pop = toolbox.population(n=20)
    print('start evaluation !!!')
    fitnesses = list( map(toolbox.evaluate, pop))
    for ind, fit in zip(pop, fitnesses):
        ind.fitness.values = fit
    print('evaled %i inds' %len(pop))
    for g in range(NGEN):
        print('--GENERATION %i--' % g)
        # select the next generation individuals
        offspring = toolbox.select( pop, len(pop) )
        # clone the selected individuals
        offspring = list( map(toolbox.clone, offspring) )
        # Apply crossover and mutaion on the offspring
        for  child1, child2 in zip( offspring[::2], offspring[1::2]):
            if random.random() < CXPB:
                toolbox.mate( child1, child2 )
                del child1.fitness.values
                del child2.fitness.values
        for mutant in offspring:
            if random.random() < MUTPB:
                toolbox.mutate(mutant)
                del mutant.fitness.values
        # Evaluate the individuals with an invalid fitness
        invalid_ind = [ ind for ind in offspring if not ind.fitness.valid ]
        fitnesses = map(toolbox.evaluate, invalid_ind)
        for ind, fit in zip(invalid_ind, fitnesses):
            ind.fitness.values = fit
        print(' evaled %i individuals '%len(invalid_ind))
        pop[:] = offspring
        fits = [ind.fitness.values[0] for ind in pop]
        length = len(pop)
        mean = sum(fits) / length
        sum2 = sum(x*x for x in fits)
        std = abs(sum2/length-mean**2)**0.5
        print('min %s' % min(fits))
        print ('max %s' % max(fits))
        print ('avg %s' % mean)
        print('std %s' % std)
    print('--end of evolution--')
    best_ind = tools.selBest(pop, 1)[0]
    print('bset individual is %s, %s' % (best_ind, best_ind.fitness.values))

Ready.
Then run script.

main()
<out put from here>
start evaluation !!!
evaled 20 inds
--GENERATION 0--
 evaled 18 individuals 
min 0.765494194159
max 0.850630899091
avg 0.818827150217
std 0.0248685727781
........
 evaled 10 individuals 
min 0.912074571976
max 0.912074571976
avg 0.912074571976
std 2.58095682795e-08
--end of evolution--
bset individual is [41, 130, 1, 113, 5, 198, 94, 130, 93, 151], (0.91207457197601005,)
<end>

Now I got answer.
Check molecules.

res = [mols[i] for i in [41, 130, 1, 113, 5, 198, 94, 130, 93, 151]]
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
img = Draw.MolsToGridImage( res,molsPerRow=5,subImgSize=(200,200))
img

Ah …
Screen Shot 2015-02-08 at 10.48.01 PM

Now I got best dissimilarity score was 0.91207457197601005.
Compare with random fashion.

sample = random.sample(range(200),10)
sample
>[37, 160, 157, 30, 163, 38, 112, 118, 135, 85]
eval_diversity(sample)
>(0.82878571084365449,)

😉
Select molecules using GA is useful because, I could make evaluate function in myself.
Today I used only maximise dissimilarity, but I can also use maximise drug likeness or activity for target, or any more…
So I’ll study more deeply about deap.

広告

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中