<jupyter_start><jupyter_text>Database connection<jupyter_code>connection = psycopg2.connect(user = "chembl",
                              password = "<PASSWORD>",
                              host = "127.0.0.1",
                              port = "5432",
                              database = "chembl_25")
cursor = connection.cursor()<jupyter_output><empty_output><jupyter_text>Types of activities<jupyter_code>cursor.execute('select min(AC.standard_type) AS TYPE, count(AC.standard_type) AS COUNT from ACTIVITIES AC GROUP BY standard_type')
activities = pd.DataFrame(cursor.fetchall(), columns = ['type', 'count'])
activities[activities['count'] > 100000]<jupyter_output><empty_output><jupyter_text>Main query<jupyter_code>columns = ['molregno', 'canonical_smiles', 'activity_id', 
           'standard_value', 'standard_units', 'standard_flag', 'standard_type', 'activity_comment',
           'alogp', 'hba', 'hbd', 'psa', 'rtb', 'ro3_pass', 'num_ros_violations', 'molecular_species', 'full_mwt', 'aromatic_rings', 'heavy_atoms', 'qed_weighted']

cursor.execute("select CS.molregno, \
               CS.canonical_smiles, \
               AC.activity_id, \
               AC.standard_value, \
               AC.standard_units, \
               AC.standard_flag, \
               AC.standard_type, \
               AC.activity_comment, \
               CP.ALOGP, \
               CP.HBA, \
               CP.HBD, \
               CP.PSA, \
               CP.RTB, \
               CP.RO3_PASS, \
               CP.NUM_RO5_VIOLATIONS, \
               CP.MOLECULAR_SPECIES, \
               CP.FULL_MWT, \
               CP.AROMATIC_RINGS, \
               CP.HEAVY_ATOMS, \
               CP.QED_WEIGHTED \
               from COMPOUND_STRUCTURES CS \
               inner join ACTIVITIES AC on CS.molregno = AC.molregno \
               inner join COMPOUND_PROPERTIES CP on CS.molregno = CP.MOLREGNO \
               and (AC.standard_type = 'IC50' or AC.standard_type = 'GI50' or AC.standard_type = 'Potency') \
               and (AC.standard_value IS NOT NULL)")
molData = pd.DataFrame(cursor.fetchall(), columns = columns)
connection.close()
molData.shape
molData.head(10).iloc[1,:]<jupyter_output><empty_output><jupyter_text>Data conversion<jupyter_code>floatDescriptors = ['standard_value', 'alogp', 'psa', 'full_mwt', 'qed_weighted']
for moldesc in floatDescriptors:
    molData[moldesc] = molData[moldesc].astype(float)

intDescriptors = ['molregno']#, 'hba', 'hbd', 'rtb', 'aromatic_rings', 'heavy_atoms']
for moldesc in intDescriptors:
    molData[moldesc] = molData[moldesc].astype(int)<jupyter_output><empty_output><jupyter_text>What units do we have?<jupyter_code># mess in units
molData.groupby('standard_units').agg('count')
# take only the entries expressed in nM
molData = molData[molData['standard_units']=='nM'].reset_index()
molData.groupby('standard_units').agg('count')
molData.head(5)<jupyter_output><empty_output><jupyter_text>Data aggregation<jupyter_code>aggFunctions = {
                'molregno': ['min', 'count'], 'canonical_smiles': 'min',
                'standard_value': ['min', 'max'],
                'standard_type': 'min',
                'alogp': ['min', 'max'],
                'hba': ['min', 'max'],
                'hbd': ['min', 'max'],
                'psa': ['min', 'max'],
                'rtb': ['min', 'max'],
                'ro3_pass': 'min',
                'num_ros_violations': 'min',
                'molecular_species': 'min',
                'full_mwt': ['min', 'max'],
                'aromatic_rings': 'min',
                'heavy_atoms': 'min',
                'qed_weighted': ['min', 'max']
                }
grouped = molData.groupby('molregno')
molDataGrouped = grouped.agg(aggFunctions).reset_index()
molDataGrouped.head(5)<jupyter_output><empty_output><jupyter_text>Activity type distribution<jupyter_code># THE GI50/IC50 distribution
molDataGrouped['standard_type'].groupby('min').agg({'min': ['count']})<jupyter_output><empty_output><jupyter_text>Pickling of the data<jupyter_code>import pickle
pcklFile = 'molDataGrouped.pckl'
with open(pcklFile, 'wb') as file:
    pickle.dump(molDataGrouped, file)
with open(pcklFile, 'rb') as file:
    molDataGrouped = pickle.load(file)
molDataGrouped.columns
molDataGrouped.loc[:, ['canonical_smiles', 'min']].values
molDataGrouped.head(10)<jupyter_output><empty_output><jupyter_text>Smiles analysis<jupyter_code>from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import Chem
from rdkit.Chem import Draw
def provideMoleculeStatistics(smiles):
    #print(smiles)
    mol = Chem.MolFromSmiles(smiles)
    
    newSmiles = Chem.MolToSmiles(mol, canonical = True, isomericSmiles = False)
    negativeCharged = sum([ item.GetFormalCharge() if item.GetFormalCharge() < 0 else 0 for item in mol.GetAtoms() ])
    positiveCharged = sum([ item.GetFormalCharge() if item.GetFormalCharge() > 0 else 0 for item in mol.GetAtoms() ])
    #anyCharged = any([item1 or item2 for item1, item2 in zip(negativelyCharged, positivelyCharged)])
    
    elementsList = list(set([atom.GetSymbol() for atom in mol.GetAtoms()]))
    numberOfRings = mol.GetRingInfo().NumRings()
    
    return(newSmiles, negativeCharged, positiveCharged, elementsList, numberOfRings)

import codecs
encodeToUTF8 = False
def canonicalizeSmilesAndProvideDescriptor(smiles):
    #rdkitMol = Chem.MolFromSmiles(molecule)
    #Chem.SanitizeMol(rdkitMol)
    try:
        #newSmiles = Chem.MolToSmiles(rdkitMol, canonical = True, isomericSmiles = False)
        newSmiles, negativeCharged, positiveCharged, elementsList, numberOfRings = provideMoleculeStatistics(smiles)
        #smilesDescription = checkVariousSmilesProperties(newSmiles)
        #elementsSet = provideElementsList(newSmiles)        
    except:
        newSmiles, negativeCharged, positiveCharged, elementsList, numberOfRings = (None, None, None, None, None)
        # There was a trouble in catching the ArgumentError exception (originatefd most likely in Boost.Python 
        # therefore any exceptio s caught here)
        print('Exception!!! :', smiles)
        
    if (encodeToUTF8):
        return((codecs.encode(newSmiles, 'utf-8'), negativeCharged, positiveCharged, elementsList, numberOfRings))
    else:
        return((newSmiles, negativeCharged, positiveCharged, elementsList, numberOfRings))
print(pd.__version__)
def printRow(row):
    print('in')
    print(row)
    mol = Chem.MolFromSmiles(row)

molDataGrouped.loc[:2, ('canonical_smiles', 'min')].apply(printRow)
 #.apply(canonicalizeSmilesAndProvideDescriptor)


len(molDataGrouped)
sourceColumn = ('canonical_smiles', 'min')
nTotal = len(molDataGrouped)
nStart = 0
nSize = 10000
nBatch = np.ceil((nTotal - nStart)/nSize).astype(int)
for iii in range(nBatch):
    iBeg = nStart + iii * nSize
    if (iii == nBatch - 1):
        iEnd = nTotal
    else:
        iEnd = nStart + (iii + 1) * nSize
    print(iii)
    result = molDataGrouped.loc[iBeg:iEnd, sourceColumn].apply(canonicalizeSmilesAndProvideDescriptor)
    molDataGrouped.loc[iBeg:iEnd, 'canonicalSmiles'] = [item[0] for item in result]
    molDataGrouped.loc[iBeg:iEnd, 'negativeCharged'] = [item[1] for item in result]
    molDataGrouped.loc[iBeg:iEnd, 'positiveCharged'] = [item[2] for item in result]
    molDataGrouped.loc[iBeg:iEnd, 'elementsSet'] = [item[3] for item in result]
    molDataGrouped.loc[iBeg:iEnd, 'numberOfRings'] = [item[4] for item in result]
molDataGrouped.head(4)
molDataGrouped[molDataGrouped['elementsSet'] == None]
with open('molDataGroupedDesc.pckl', 'wb') as file:
    pickle.dump(molDataGrouped, file)
with open('molDataGroupedDesc.pckl', 'rb') as file:
    molDataGrouped = pickle.load(file)
molDataGrouped = molDataGrouped[~molDataGrouped['elementsSet'].isnull()]
organicChemistryList = ['B', 'C', 'N', 'O', 'P', 'S', 'F', 'Cl', 'Br', 'I']
organicChemistrySet = set(organicChemistryList)
testSet = set(['N', 'C', 'Cl'])
testSet < organicChemistrySet
molDataGrouped['organicChemistryElementsOnly'] = molDataGrouped['elementsSet'].apply(lambda x: set(x) < organicChemistrySet)
molDataGrouped.groupby('organicChemistryElementsOnly').count()
molDataGrouped['canonicalSmilesLength'] = molDataGrouped['canonicalSmiles'].apply(lambda x: len(x))
plt.hist(molDataGrouped['canonicalSmilesLength'], bins = 50)
limitSmilesLength = 100
molDataGroupedChosen = molDataGrouped[(molDataGrouped['canonicalSmilesLength'] < limitSmilesLength) & \
                      (molDataGrouped['negativeCharged'] == 0) & \
                      (molDataGrouped['positiveCharged'] == 0) & \
                      (molDataGrouped['numberOfRings'] <= 5) & \
                      (molDataGrouped['organicChemistryElementsOnly'])].reset_index()
replacementDict = {'Br': 'G', 'Cl': 'U', '[nH]': 'W'}
molDataGroupedChosen['encodedSmiles'] = molDataGroupedChosen['canonicalSmiles'].replace(replacementDict, regex=True)
molDataGroupedChosen.head(5)
molDataGroupedChosen.shape
# molDataGroupedFinal_100.pckl: smiles length < 100
with open('molDataGroupedFinal_100.pckl', 'wb') as file:
    pickle.dump(molDataGroupedChosen, file)
with open('molDataGroupedFinal.pckl', 'rb') as file:
    molDataGroupedChosen = pickle.load(file)

molDataGroupedChosen.head()
nSmilesCodes = 50000
mask = random.randint(0, molDataGroupedChosen.shape[0], size=nSmilesCodes)
staticFeatures = pd.DataFrame()
toBeAveraged = ['standard_value', 'alogp', 'hba', 'hbd', 'psa', 'rtb', 'full_mwt', 'qed_weighted']
for quantity in toBeAveraged:
    staticFeatures.loc[:, quantity] = (molDataGroupedChosen.loc[mask, (quantity, 'min')] + molDataGroupedChosen.loc[mask, (quantity, 'max')])/2
    staticFeatures.loc[:, quantity].astype(float)
toBeTaken = ['aromatic_rings', 'heavy_atoms']
for quantity in toBeTaken:
    staticFeatures.loc[:, quantity] = molDataGroupedChosen.loc[mask, (quantity, 'min')]
    staticFeatures.loc[:, quantity].astype(float)
staticFeatures.loc[:, 'number_of_rings'] = molDataGroupedChosen.loc[mask, 'numberOfRings'].astype(float)
staticFeatures.loc[:, 'number_of_rings'] = staticFeatures.loc[:, 'number_of_rings'].astype(float)
print(staticFeatures.head(2))
#staticFeatures = staticFeatures.values
staticFeatures['full_mwt'] = staticFeatures.full_mwt.astype(float)
staticFeatures['qed_weighted'] = staticFeatures.qed_weighted.astype(float)
staticFeatures['aromatic_rings'] = staticFeatures.qed_weighted.astype(float)

thres = 100000
staticFeatures[staticFeatures['standard_value'] < thres].shape[0] / staticFeatures['standard_value'].shape[0]
staticFeatures = staticFeatures[staticFeatures['standard_value'] < thres]
staticFeatures.shape
allDescriptors = ['standard_value', 'alogp', 'hba', 'hbd', 'psa', 'rtb', 'full_mwt', 'qed_weighted', 'aromatic_rings', 'heavy_atoms', 'number_of_rings']
#allDescriptors = ['standard_value']
#quantity = 'alogp'
plotIdx = 1
nRows = np.ceil(len(allDescriptors) / 2)
fig = plt.figure(figsize=(16, 16)) 
for quantity in allDescriptors:
    print(quantity)
    plt.subplot(nRows, 2, plotIdx)
    plt.hist(staticFeatures[~staticFeatures[quantity].isnull()][quantity], bins = 10)
    plt.title(quantity)
    plotIdx += 1<jupyter_output><empty_output><jupyter_text>Dynamic features<jupyter_code>smilesCodes = molDataGroupedChosen.loc[staticFeatures.index, 'encodedSmiles']
smilesCodes
maxlen = -1
for code in smilesCodes:
    if len(code) > maxlen:
        maxlen = len(code)
maxlen
minlen = 1e6
for code in smilesCodes:
    if len(code) < minlen:
        minlen = len(code)
minlen
def pad_smile(string, max_len, padding='right'):
    if len(string) <= max_len:
        if padding == 'right':
            return string + " " * (max_len - len(string))
        elif padding == 'left':
            return " " * (max_len - len(string)) + string
        elif padding == 'none':
            return string

smilesCodes = smilesCodes.apply(lambda x: pad_smile(x, max_len=maxlen, padding='right'))
chars = sorted(list(set(smilesCodes.str.cat(sep=''))))
print('total chars:', len(chars))
char2indices = dict((c, i) for i, c in enumerate(chars))
indices2char = dict((i, c) for i, c in enumerate(chars))
dynamicFeatures = np.zeros((len(smilesCodes), maxlen, len(chars)), dtype=np.float)
dynamicFeatures.shape
for codeidx, code in enumerate(smilesCodes):
    for charidx, char in enumerate(code):
        dynamicFeatures[codeidx, charidx, char2indices[char]] = 1
sums = []
for idx in range(dynamicFeatures.shape[0]):
    sums.append(np.sum(dynamicFeatures[idx, :, :]))
plt.hist(sums)
staticChosen = ['alogp', 'full_mwt']
scaler = StandardScaler()
scaler.fit(staticFeatures[staticChosen])
staticFeaturesStandard = scaler.transform(staticFeatures[staticChosen])<jupyter_output><empty_output><jupyter_text>Autoencoder architecture<jupyter_code>from keras.layers import LSTM, TimeDistributed, concatenate, Input, Dense, RepeatVector, Lambda
from keras.models import Model
from keras.activations import relu, sigmoid, tanh
from keras.callbacks import EarlyStopping, ModelCheckpoint
from keras.optimizers import Adam, RMSprop
import keras.backend as K
from keras.utils import plot_model
from keras import losses
import numpy.random as rnd
def sampling(args):
    z_mean, z_log_var = args
    batch = K.shape(z_mean)[0]
    dim = K.int_shape(z_mean)[1]
    # by default, random_normal has mean=0 and std=1.0
    epsilon = K.random_normal(shape=(batch, dim))
    return z_mean + K.exp(0.5 * z_log_var) * epsilon
def prepare_model(static, dynamic, k, window, charsetLen, lr, lossFunction, showArch):
    input_dynamic = Input(shape=(window, charsetLen), name="input_dynamic")
    input_static = Input(shape=(static,), name="input_static")
    latent = Dense(k[0], activation=relu)(input_static)
    dense_h = Dense(k[0])(latent)
    dense_c = Dense(k[0])(latent)
    lstm_layer, state_h, state_c = LSTM(k[0], return_sequences=True, return_state=True)(input_dynamic,
                                                                                        initial_state=[dense_h,
                                                                                                       dense_c])

    for x in k[1:-1]:
        concat_h = concatenate([dense_h, state_h])
        dense_h = Dense(x)(concat_h)
        concat_c = concatenate([dense_c, state_c])
        dense_c = Dense(x)(concat_c)
        lstm_layer, state_h, state_c = LSTM(x, return_sequences=True, return_state=True)(lstm_layer,
                                                                                         initial_state=[dense_h,
                                                                                                        dense_c])
    x = k[-1]
    concat_h = concatenate([dense_h, state_h])
    dense_h = Dense(x)(concat_h)
    concat_c = concatenate([dense_c, state_c])
    dense_c = Dense(x)(concat_c)
    lstm_layer, state_h, state_c = LSTM(x, return_state=True)(lstm_layer, initial_state=[dense_h, dense_c])
    concat = concatenate([lstm_layer, latent])

    # autoencoder
    z_mean = Dense(x, name='z_mean')(concat)
    z_log_var = Dense(x, name='z_log_var')(concat)
    
    
    z = Lambda(sampling, output_shape=(x,), name='z')([z_mean, z_log_var])
    
    
    
    state_h = Dense(k[-2], activation=relu)(z)
    dense_h = Dense(k[-2], activation=relu)(z)
    state_c = Dense(k[-2], activation=relu)(z)
    dense_c = Dense(k[-2], activation=relu)(z)
    lstm_layer = RepeatVector(window)(z)

    for x in np.flip(k[:-1]):
        concat_h = concatenate([dense_h, state_h])
        dense_h = Dense(x)(concat_h)
        concat_c = concatenate([dense_c, state_c])
        dense_c = Dense(x)(concat_c)
        lstm_layer, state_h, state_c = LSTM(x, return_sequences=True, return_state=True)(lstm_layer,
                                                                                         initial_state=[dense_h,
                                                                                                        dense_c])

    #result_series = TimeDistributed(Dense(charsetLen))(lstm_layer)
    result_series = LSTM(charsetLen, return_sequences=True, activation='softmax')(lstm_layer)
    concat = concatenate([state_h, state_c])
    #result_sigmoid = Dense(static-3, activation=sigmoid)(concat)
    result_relu = Dense(static, activation=sigmoid)(concat)

    #model = Model(inputs=[input_dynamic, input_static], outputs=[result_series, result_sigmoid, result_relu])
    model = Model(inputs=[input_dynamic, input_static], outputs=[result_series, result_relu])
    optimizer = RMSprop(lr=lr)
    model.compile(optimizer=optimizer, loss=lossFunction, metrics=['binary_crossentropy', 'mean_absolute_error'])
    if (showArch):
        print(model.summary())
    return model
def fit(staticFeatures, dynamicFeatures, model, step=1):
    #dynamic_data = np.empty((0, window, 1), np.float)
    #helper = []
    #for d in dynamic:
    #    new_data = rolling_window(d, window, step)
    #    helper.append(len(new_data))
    #    dynamic_data = np.append(dynamic_data, new_data, axis=0)
    #print(len(helper))
    #static_data = np.repeat(static, helper, axis=0)
    order = rnd.permutation(len(staticFeatures))

    early_stopping = EarlyStopping(monitor='val_loss', patience=5)
    bst_model_path = 'autoencoder.h5'
    checkpoint = ModelCheckpoint(bst_model_path, save_best_only=True, save_weights_only=True, monitor='val_loss')

    size = int(staticFeatures.shape[0] * 0.9)
    training_dynamic, training_static = dynamicFeatures[order[:size]], staticFeatures[order[:size]]
    testing_dynamic, testing_static = dynamicFeatures[order[size:]], staticFeatures[order[size:]]
    print(training_dynamic.shape, training_static.shape)
    print(testing_dynamic.shape, testing_static.shape)
    model.fit([training_dynamic, training_static], 
              [training_dynamic, training_static],
                   epochs=10,
                   batch_size=64,
                   callbacks=[early_stopping, checkpoint],
                   validation_data=([testing_dynamic, testing_static], 
                                    [testing_dynamic, testing_static]))
lr = 0.001
model = prepare_model(staticFeaturesStandard.shape[1], 1, [64,64], dynamicFeatures.shape[1], dynamicFeatures.shape[2], lr, ['binary_crossentropy', 'mean_absolute_error'], True)
fit(staticFeaturesStandard, dynamicFeatures, model)
dynamicFeatures[0,:,:].reshape(-1, dynamicFeatures.shape[1], dynamicFeatures.shape[2])
staticFeaturesStandard[0,:].reshape(-1, staticFeaturesStandard.shape[1]).shape
prediction = model.predict([dynamicFeatures[0,:,:].reshape(-1, dynamicFeatures.shape[1], dynamicFeatures.shape[2]), staticFeaturesStandard[0,:].reshape(-1, staticFeaturesStandard.shape[1]) ])
prediction[0]
import SmilesEnumerator<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup