<jupyter_start><jupyter_text>Database extraction<jupyter_code>DF_smiles = pd.read_csv('DF_GI50_IC50_grouped.csv', sep = ',', header = [0, 1])
#DF_smiles['aaa'] = 1
DF_smiles.columns = ['molregNoMin', 'molregNoCount', 'chemblSmiles', 'standardMin', 'standardMax', 'standardMean',
                     'publishedMain', 'publishedMax', 'publishedMean']
DF_smiles.head(10)<jupyter_output><empty_output><jupyter_text>Various methods involving RDkit<jupyter_code>def checkElementPresence(element, smiles):
    mol = Chem.MolFromSmiles(smiles)
    return(element in set([atom.GetSymbol() for atom in mol.GetAtoms()]))
def provideElementsList(smiles):
    mol = Chem.MolFromSmiles(smiles)
    return(list(set([atom.GetSymbol() for atom in mol.GetAtoms()])))
def checkChargePresence(smiles):
    mol = Chem.MolFromSmiles(smiles)
    negativelyCharged = [item.GetFormalCharge() < 0 for item in mol.GetAtoms()]
    positivelyCharged = [item.GetFormalCharge() > 0 for item in mol.GetAtoms()]
    anyCharged = [item1 or item2 for item1, item2 in zip(negativelyCharged, positivelyCharged)]
    return(any(anyCharged))
def provideMoleculeStatistics(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)<jupyter_output><empty_output><jupyter_text>Drawing the molecule<jupyter_code>mol = Chem.MolFromSmiles(DF_smiles.loc[1, 'chemblSmiles'])
AllChem.Compute2DCoords(mol)
AllChem.GenerateDepictionMatching2DStructure(mol, mol)
Draw.MolToImage(mol, fitImage = True, size = (200, 200))<jupyter_output><empty_output><jupyter_text>Number of rings<jupyter_code>mol.GetRingInfo().NumRings()
def checkVariousSmilesProperties(smiles):
    if (checkChargePresence(smiles)):
        response = 100
    elif (checkElementPresence('Br', smiles)):
        response = 10
    elif (checkElementPresence('Cl', smiles)):
        response = 11
    elif (checkElementPresence('I', smiles)):
        response = 12
    elif (checkElementPresence('As', smiles)):
        response = 30
    elif (checkElementPresence('Se', smiles)):
        response = 31
    elif (checkElementPresence('Te', smiles)):
        response = 32
    elif (checkElementPresence('Si', smiles)):
        response = 33
    elif (checkElementPresence('K', smiles)):
        response = 40
    elif (checkElementPresence('Li', smiles)):
        response = 41
    elif (checkElementPresence('Na', smiles)):
        response = 42
    elif (checkElementPresence('Zn', smiles)):
        response = 43
    elif (checkElementPresence('Mg', smiles)):
        response = 44
    elif (checkElementPresence('Au', smiles)):
        response = 50
    elif (checkElementPresence('V', smiles)):
        response = 51
    else:
        response = 0
    return(response)
def checkVariousSmilesProperties(smiles):
    if (checkChargePresence(smiles)):
        response = 100
    else:
        response = 0
    return(response)
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))
len(DF_smiles)
#nTotal = len(DF_smiles)
nTotal = 600000
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 = DF_smiles.loc[iBeg:iEnd, 'chemblSmiles'].apply(canonicalizeSmilesAndProvideDescriptor)
    DF_smiles.loc[iBeg:iEnd, 'canonicalSmiles'] = [item[0] for item in result]
    DF_smiles.loc[iBeg:iEnd, 'negativeCharged'] = [item[1] for item in result]
    DF_smiles.loc[iBeg:iEnd, 'positiveCharged'] = [item[2] for item in result]
    DF_smiles.loc[iBeg:iEnd, 'elementsSet'] = [item[3] for item in result]
    DF_smiles.loc[iBeg:iEnd, 'numberOfRings'] = [item[4] for item in result]<jupyter_output>0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Exception!!! : Clp1(Cl)np(Cl)(Cl)np2(NNP(=O)(NN2)Oc3ccccc3)n1
41
Exception!!! : CC(Br)C(=O)N=N1=C2C(=Nc3ccccc13)c4cccc5cccc2c45
42
Exception!!! : Clp1(Cl)np(Cl)(Cl)np2(NNP(=S)(NN2)Oc3ccccc3)n1
Exception!!! : CCN1\C(=C\C#C\C(=C/c2sc3ccccc3[n+]2CC)\C)\Sc4ccccc14.[F-][PH2+5]([F-])([F-])([F-])([F-])[F-]
Exception!!! : N1c2ccccc2p3(c4ccccc14)c5ccccc5nc6ccccc36
43
Exception!!! : S=P1(NNp2(NN1)np(np(n2)(N3CC3)N4CC4)(N5CC5)N6CC6)Oc7ccccc7
44
45
46
Exception!!! : N.[Cl-].[Cl-].[Cl-].C1C[O-][Te+4][O-]1
Exception!!! : N.[Cl-].[Cl-].[Cl-].CCCCC1C[O-][Te+4][O-]1
Exception!!! : N.[Cl-].[Cl-].[Cl-].CCCCCCC1C[O-][Te+4][O-]1
47
48
49
50
51
52
53
Exception!!! : CCC1CCCCN1C(=O)N2=NC(=CN2)C(O)(c3ccccc3)c4ccccc4
54
55
56
57
58
Exception!!! : CN1C(=O)NC2=CN3(=C4NC=CC4=C12)CCCCC(C3)N5C(=O)Nc6cnc7[nH]ccc7c56
59<jupyter_text>Get only postprocessed rows<jupyter_code>DF_smiles
DF_smiles = DF_smiles[[not item for item in pd.isna(DF_smiles['canonicalSmiles'])]]
#DF_smiles = DF_smiles[DF_smiles['canonicalSmiles'].notnull()]
DF_smiles.to_csv('DF_GI50_IC50_groupedSanitized.csv', index = False)
DF_smiles = pd.read_csv('DF_GI50_IC50_groupedSanitized.csv')
import ast
DF_smiles.loc[:, 'elementsSet'] = DF_smiles.loc[:, 'elementsSet'].map(lambda x: ast.literal_eval(x))
DF_smiles.head(10)
organicChemistryList = ['B', 'C', 'N', 'O', 'P', 'S', 'F', 'Cl', 'Br', 'I']
organicChemistrySet = set(organicChemistryList)<jupyter_output><empty_output><jupyter_text>Subset check test<jupyter_code>testSet = set(['N', 'C', 'Cl', 'Na'])
testSet
testSet < organicChemistrySet<jupyter_output><empty_output><jupyter_text>Check if smile contains only the elements from the orgranic chemistry set<jupyter_code>DF_smiles.loc[:, 'organicChemistryElementsOnly'] = DF_smiles.loc[:, 'elementsSet'].map(lambda x: set(x) < organicChemistrySet)
DF_smiles.head(10)
plt.hist(DF_smiles['organicChemistryElementsOnly'], bins = 2, width = 0.3)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Length of smiles<jupyter_code>DF_smiles.loc[:, 'smilesLength'] = DF_smiles.loc[:, 'canonicalSmiles'].map(lambda x: len(x))
plt.hist(DF_smiles.loc[:, 'smilesLength'], bins = 50)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Get only the smiles that: are not longer than 100 characters, are not charged and are composed from organic chemistry elements<jupyter_code>limitSmilesLength = 100
DF_chosen = DF_smiles[(DF_smiles['smilesLength'] < limitSmilesLength) & \
                      (DF_smiles['negativeCharged'] == 0) & \
                      (DF_smiles['positiveCharged'] == 0) & \
                      (DF_smiles['numberOfRings'] <= 5) & \
                      (DF_smiles['organicChemistryElementsOnly'])]
DF_chosen<jupyter_output><empty_output><jupyter_text>Encode two-character elements and [nH] pieces<jupyter_code>def encodeElement(series, source, destination):
    return series.map(lambda x: x.replace(source, destination))
DF_chosen.loc[:, 'encodedSmiles'] = DF_chosen.loc[:, 'canonicalSmiles']<jupyter_output>C:\Tools\Anaconda3\lib\site-packages\pandas\core\indexing.py:357: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  self.obj[key] = _infer_fill_value(value)
C:\Tools\Anaconda3\lib\site-packages\pandas\core\indexing.py:537: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  self.obj[item] = s<jupyter_text>'Br' -> 'G'<jupyter_code>DF_chosen.loc[:, 'encodedSmiles'] = encodeElement(DF_chosen.loc[:, 'encodedSmiles'], 'Br', 'G')<jupyter_output>C:\Tools\Anaconda3\lib\site-packages\pandas\core\indexing.py:537: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  self.obj[item] = s<jupyter_text>'Cl' -> 'U'<jupyter_code>DF_chosen.loc[:, 'encodedSmiles'] = encodeElement(DF_chosen.loc[:, 'encodedSmiles'], 'Cl', 'U')<jupyter_output>C:\Tools\Anaconda3\lib\site-packages\pandas\core\indexing.py:537: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  self.obj[item] = s<jupyter_text>[nH] -> 'W'<jupyter_code>DF_chosen.loc[:, 'encodedSmiles'] = encodeElement(DF_chosen.loc[:, 'encodedSmiles'], '[nH]', 'W')<jupyter_output>C:\Tools\Anaconda3\lib\site-packages\pandas\core\indexing.py:537: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  self.obj[item] = s<jupyter_text>Check the charset<jupyter_code>from functools import reduce
def provideSampleCharset(data):
    structures = data.tolist()
    charset = list(reduce(lambda x, y: set(y) | x, structures, set()))
    return(sorted(charset))
print(provideSampleCharset(DF_chosen['encodedSmiles']), len(provideSampleCharset(DF_chosen['encodedSmiles'])))
def findSmilesWithCertainCharacters(data, sequence):
    result = []
    for index, item in enumerate(data):
        if (item.find(sequence) != -1):
            result.append((index, item))
    return(result)<jupyter_output><empty_output><jupyter_text>Remove the '-', the redundant aliphatic bond<jupyter_code>DF_chosen['encodedSmiles'] = DF_chosen['encodedSmiles'].map(lambda x: x.replace('-', ''))
findSmilesWithCertainCharacters(DF_chosen['encodedSmiles'], '-')
print(provideSampleCharset(DF_chosen['encodedSmiles']), len(provideSampleCharset(DF_chosen['encodedSmiles'])))
smilesSample = list(DF_chosen.sample(200)['canonicalSmiles'].values)
mols = []
for smiles in smilesSample:
    mol = Chem.MolFromSmiles(smiles)
    AllChem.Compute2DCoords(mol)
    mols.append(mol)
#Draw.MolsToImage()
Draw.MolsToGridImage(mols, molsPerRow=10)
#Draw.MolsToImage(mols, subImgSize=(200, 200))
maxLenSmiles = np.max(DF_chosen.loc[:, 'smilesLength'])
maxLenSmiles
plt.hist(DF_chosen.loc[:, 'smilesLength'], bins = 20)
plt.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Store on file<jupyter_code>DF_chosen.to_csv('DF_GI50_IC50_groupedSanitizedChosen.csv', index = False)<jupyter_output><empty_output><jupyter_text>Molecular descriptor<jupyter_code>from rdkit import DataStructs
from rdkit import Chem
from rdkit.Chem.Fingerprints import FingerprintMols
from rdkit.Chem.rdMolDescriptors import GetMorganFingerprint
import codecs
encodeToUTF8 = False
def canonicalizeSmiles(molecule):
    rdkitMol = Chem.MolFromSmiles(molecule)
    try:
        newSmiles = Chem.MolToSmiles(rdkitMol)
    except ArgumentError:
        newSmiles = None
    if (encodeToUTF8):
        return(codecs.encode(newSmiles, 'utf-8'))
    else:
        return(newSmiles)
list(DF_smiles.loc[:5, 'canonicalSmiles'])
mols = [Chem.MolFromSmiles(smile) for smile in DF_smiles.loc[:5, 'canonicalSmiles']]
fps = [GetMorganFingerprint(x, 2) for x in mols]
print(DF_smiles.loc[0, 'canonicalSmiles'])
print(Chem.MolToSmiles(mols[0]))<jupyter_output><empty_output><jupyter_text>Check with rdkit<jupyter_code>import codecs
encodeToUTF8 = False
def canonicalizeSmiles(molecule):
    rdkitMol = Chem.MolFromSmiles(molecule)
    newSmiles = Chem.MolToSmiles(rdkitMol)
    if (encodeToUTF8):
        return(codecs.encode(newSmiles, 'utf-8'))
    else:
        return(newSmiles)
smiles_column = 'canonical_smiles_RDKit'
data.loc[:, smiles_column] = data.loc[:, 'canonical_smiles'].apply(canonicalizeSmiles)
data.loc[:10, smiles_column]<jupyter_output><empty_output><jupyter_text>Preparation for Autoencoders Get just structures, left-justified<jupyter_code># left-adjusted, padded with space to equalize the max. length 
structures = data[smiles_column].map(lambda x: list(x.ljust(maxLenSmiles)))
type(structures)
print(structures[0])
iii = 0
for structure in structures:
    isSpace = any([item == ' ' for item in structure])
    if (isSpace):
        print(iii, structure)
    iii += 1
    #any(' ' in structure)
from sklearn.model_selection import train_test_split
trainIdx, testIdx = map(np.array, train_test_split(structures.index, test_size = 0.20))
from functools import reduce<jupyter_output><empty_output><jupyter_text>How reduce works<jupyter_code># The calculation sequence is the following:
# 1. tmp1 = set('a') | set('x')
# 2. tmp2 = set('b') | tmp1
# 3. tmp3 = set('c') | tmp2
tmp = ['a', 'b', 'c']
reduce(lambda x, y: set(y) | x, tmp, set('x'))<jupyter_output><empty_output><jupyter_text>Create the charset<jupyter_code>charset = list(reduce(lambda x, y: set(y) | x, structures, set()))
type(charset)
nCharset = len(charset)
nCharset
print(charset)<jupyter_output><empty_output><jupyter_text>One-hot encoding of the molecules  Auxiliary functions<jupyter_code>def one_hot_array(i, n):
    return map(int, [ix == i for ix in range(n)])

def one_hot_index(vec, charset):
    return map(charset.index, vec)<jupyter_output><empty_output><jupyter_text>The generator that, for given SMILES code, provides the one-hot representation of it* The structures have the same lengths (*maxLenSmiles*)* To get the data one needs to iterate over characters in SMILES<jupyter_code>def oneHotEncodedMolecule(structures, charset):
    for smiles in structures:
        yield(map(lambda x: one_hot_array(x, len(charset)), one_hot_index(smiles, charset)))<jupyter_output><empty_output><jupyter_text>Storing the data on file<jupyter_code>import h5py
outputFile = 'postprocessed.h5'
h5f = h5py.File(outputFile, 'w')
h5f.create_dataset('charset', data = [codecs.encode(item, encoding = 'utf-8') for item in charset])
h5f.close()<jupyter_output><empty_output><jupyter_text>Function that stores the matrices on the HDF5 file<jupyter_code>def storeData(savingStructures, outputFile, dataSetName, chunkSize = 1000):

    h5f = h5py.File(outputFile, 'r+')
    dataset_shape = (len(savingStructures), maxLenSmiles, len(charset))
    newData = h5f.create_dataset(dataSetName, dataset_shape, chunks = tuple([chunkSize]+list(dataset_shape[1:])))

    nChunks = np.ceil(len(savingStructures) / chunkSize)
    #dataChunk = np.zeros(shape = (chunkSize, maxLenSmiles, nCharset))
    
    for jjj in range(int(nChunks)):
        beg = jjj * chunkSize
        if (jjj == nChunks - 1):
            end = len(savingStructures)
        else:
            end = (jjj + 1) * chunkSize
        print(beg, end)
        smilesGenerator = oneHotEncodedMolecule(savingStructures[beg:end], charset)
        iii = 0
        while True:
            try: 
                newData[jjj*chunkSize + iii,] = np.array([list(item) for item in next(smilesGenerator)])
                iii += 1
            except StopIteration:
                break

    h5f.close()
storeData(savingStructures = structures[trainIdx], dataSetName ='data_train', outputFile = outputFile)<jupyter_output><empty_output><jupyter_text>Check if data is written correctly, reading in the charset and molecular matrices<jupyter_code>inputFile = 'postprocessed.h5'
h5f = h5py.File(inputFile, 'r')
readCharset = h5f.get('charset')[:]
data = h5f.get('data_train')[:]
#np.array(data, shape = ())
h5f.close()
data.shape<jupyter_output><empty_output><jupyter_text>Decoding the charset<jupyter_code>import codecs
#codecs.decode(readCharset, 'ascii')
readCharset = [codecs.decode(item, encoding = 'ascii') for item in readCharset]<jupyter_output><empty_output><jupyter_text>Decoding molecules<jupyter_code>def from_one_hot_array(vec):
    oh = np.where(vec == 1)
    if oh[0].shape == (0, ):
        return None
    return int(oh[0][0])

def decode_smiles_from_indexes(vec, charset):
    return "".join(map(lambda x: charset[x], vec)).strip()<jupyter_output><empty_output><jupyter_text>Check if the first read molecule is the same as the first training instance<jupyter_code>smilesCode = [from_one_hot_array(data[0, iii, :]) for iii in range(maxLenSmiles)]
decode_smiles_from_indexes(smilesCode, readCharset)
"".join(structures[trainIdx[0]])<jupyter_output><empty_output><jupyter_text>Same for other test instace<jupyter_code>smilesGenerator = oneHotEncodedMolecule([structures[0]], charset)
dataOne = np.array([list(item) for item in next(smilesGenerator)])
dataOne.shape
smilesCode = [from_one_hot_array(dataOne[iii, :]) for iii in range(maxLenSmiles)]
decode_smiles_from_indexes(smilesCode, charset)
"".join(structures[0])<jupyter_output><empty_output><jupyter_text>RNN-molecule* 1 A. Gupta, <NAME>, <NAME>, <NAME>, <NAME>, and <NAME>, Molecular Informatics 37, 1700111 (2018).<jupyter_code>data
data[data['smiles_length'] > 150]<jupyter_output><empty_output><jupyter_text>Neural Network<jupyter_code>from keras.models import Sequential
from keras.layers import Dense, LSTM, Activation
from keras.optimizers import RMSprop


model = Sequential()
model.add(LSTM(128, input_shape=(maxLenSmiles, len(readCharset)), return_sequences = False))
#model.add(Dropout(0.2))
#model.add(LSTM(128, input_shape=(maxlen, len(chars))))
model.add(Dense(len(readCharset)))
model.add(Activation('softmax'))

optimizer = RMSprop(lr=0.01)
model.compile(loss='categorical_crossentropy', optimizer=optimizer)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13457632414978388, 'fraction_numerical': 0.023175876256444978, 'mean_word_length': 5.452284946236559, 'pattern_counts': {'":': 0, '<': 135, '<?xml version=': 0, '>': 135, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '440214', 'license': 'MIT', 'n_tokens_mistral': 7301, 'n_tokens_neox': 6893, 'n_words': 1340, 'path': 'src/notebooks/smilesAnalysis.ipynb', 'repo_name': 'rafalbachorz/chemicalSmilesSpace', 'size': '19230'}