6. Molecular Similarity#
Objectives#
Generete molecular fingerprints for a given molecule.
Evaluate structural similarity between molecules using different molecular fingerpints and similarity metrics.
Many useful documents/papers describe various aspects of molecular similarity, including molecular fingerprints and similarity measures. Please read these if you need more details.
Getting Started with the RDKit in Python
(https://www.rdkit.org/docs/GettingStartedInPython.html#fingerprinting-and-molecular-similarity)Fingerprint Generation, GraphSim Toolkit 2.4.2
(https://docs.eyesopen.com/toolkits/python/graphsimtk/fingerprint.html)Chemical Fingerprints
(https://docs.chemaxon.com/display/docs/Chemical+Fingerprints)Extended-Connectivity Fingerprints
(https://doi.org/10.1021/ci100050t)
1. Fingerprint Generation#
from rdkit import Chem
mol = Chem.MolFromSmiles('CC(C)C1=C(C(=C(N1CC[C@H](C[C@H](CC(=O)O)O)O)C2=CC=C(C=C2)F)C3=CC=CC=C3)C(=O)NC4=CC=CC=C4')
1-(1) MACCS keys#
The MACCS key is a binary fingerprint (a string of 0’s and 1’s). Each bit position represents the presence (=1) or absence (=0) of a pre-defined structural feature. The feature definitions for the MACCS keys are available at:
rdkit/rdkit
from rdkit.Chem import MACCSkeys
fp = MACCSkeys.GenMACCSKeys(mol)
print(type(fp))
for i in range(len(fp)):
print(fp[i], end='')
fp.ToBitString() # Alternative, easier way to convert it to a bitstring.
<class 'rdkit.DataStructs.cDataStructs.ExplicitBitVect'>
00000000000000000000000000000000000000000010000000000100000000100100000000110000100101010111100011001000100110110000011001110100110111111101101011111111111111111111110
'00000000000000000000000000000000000000000010000000000100000000100100000000110000100101010111100011001000100110110000011001110100110111111101101011111111111111111111110'
len(fp)
167
Note that the MACCS key is 166-bit-long, but RDKit generates a 167-bit-long fingerprint. It is because the index of a list/vector in many programming languages (including python) begins at 0. To use the original numbering of the MACCS keys (1-166) (rather than 0-165), the MACCS keys were implemented to be 167-bit-long, with Bit 0 being always zero. Because Bit 0 is set to OFF for all compounds, it does not affect the evaluation of molecular similarity.
These are some methods that allow you to get some additional information on the MACCS Keys.
print(fp.GetNumBits())
print(fp.GetNumOffBits())
print(fp.GetNumOnBits())
print(fp.ToBinary())
167
105
62
b'\xe0\xff\xff\xff\xa7\x00\x00\x00>\x00\x00\x00T\x14\x10\x04\x10\x00\x08\x04\x02\x02\x02\x00\x00\x00\x06\x00\x04\x06\x04\x00\x02\x00\n\x00\x04\x00\x00\x02\x04\x00\x02\x00\x00\x00\x00\x00\x00\x02\x00\x02\x02\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x02'
Exercise 1a: Generate the MACCS keys for the molecules represented by the following SMILES, and get the positions of the bits set to ON in each of the three fingerprints. What fragments do these bit positions correspond to? (the bit definitions are available at
smiles = [ 'C1=CC=CC=C1', # Benzene (Kekule)
'c1ccccc1', # Benzene ("Aromatized" carbons)
'C1CCCCC1'] # Cyclohexene
# Write your code in this cell.
Write the fragment definition of the bits ON (one is already provided for you as an example).
1-(2) Circular Fingerprints#
Circular fingerprints are hashed fingerprints. They are generated by exhaustively enumerating “circular” fragments (containing all atoms within a given radius from each heavy atom of the molecule) and then hashing these fragments into a fixed-length bitstring. (Here, the “radius” from an atom is measured by the number of bonds that separates two atoms).
Examples of circular fingerprints are the extended-connectivity fingerprint (ECFPs) and their variant called FCFPs (Functional-Class Fingerprints), originally described in a paper by Rogers and Hahn (https://doi.org/10.1021/ci100050t). The RDKit implementation of these fingerprints are called “Morgan Fingerprints” (https://www.rdkit.org/docs/GettingStartedInPython.html#morgan-fingerprints-circular-fingerprints).
from rdkit.Chem import AllChem
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=1024).ToBitString()
print(fp)
0100000000000000000000000000010001100000000000000000010000000000100000000000000010000000001000000000000000000100000001000100000010000000100000000000000010000000000000000100000100000000000000010000000000000000000000000000000000010000000000000000000000000000000000000010000000000000000100000100000000000000000000000000000000000000000000000000000000000001000010010000000000000000000000000010010000000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100000000000100000000000000000000000001000000000000000000000100010000000000000000000000000010010000000000000000001000100000000010000000001001000000000100100000000000000000000000000000000000000000000000000001000000000100001000000000000000000000000000000000000000001000000000000000000000000000000000000000000010010000000100000000000000000000000000000010000000000000000001000000000000000000100000000000000000000000000000000000000001000000000
[10:36:26] DEPRECATION WARNING: please use MorganGenerator
When comparing the RDK’s Morgan fingerprints with the ECFP/FCFP fingerprints, it is important to remember that the name of ECFP/FCFP fingerprints are suffixed with the diameter of the atom environments considered, while the Morgan Fingerprints take a radius parameter (e.g., the second argument “2” of GetMorganFingerprintAsBitVect() in the above code cell). The Morgan fingerprint generated above (with a radius of 2) is comparable to the ECFP4 fingerprint (with a diameter of 4).
Exercise 1b: For the moleculess below, generate the 512-bit-long Morgan Fingeprint comparable to the FCFP6 fingerprint.
Search for the compounds by name and get their SMILES strings.
Generate the molecular fingerprints from the SMILES strings.
Print the generated fingerprints.
To generate FCFP (not ECFP), read the following document: https://www.rdkit.org/docs/GettingStartedInPython.html#morgan-fingerprints-circular-fingerprints
synonyms = [ 'diphenhydramine', 'cetirizine', 'fexofenadine', 'loratadine' ]
# Write your code in this cell
1-(3) Path-Based Fingeprints#
Path-based fingerprints are also hashed fingerprints. They are generated by enumerating linear fragments of a given length and hashing them into a fixed-length bitstring. An example is the RDKit’s topological fingeprint. As described in the RDK documentation (https://www.rdkit.org/docs/GettingStartedInPython.html#topological-fingerprints), while this fingerprint can be generated using FingerprintMols.FingerprintMol(), it is recommended to use rdmolops.RDKFingerprint() to generate the fingerprint using non-default parameter values.
from rdkit.Chem import rdmolops
fp = rdmolops.RDKFingerprint(mol, fpSize=2048, minPath=1, maxPath=7).ToBitString()
print(fp)
11011101111101010111110010111110101111011110100111111111110111001111111101101111011111111111101111111011111101111101010010110111011111011111111111011011011111101110001101111111000110010101111111110110111011010111110101111101111001111111011111101111110011101110111010010110001110001111011001111110001111111001111011111110000100001111101111110111111111000111111111111110111001010111011111111101111011111010000110110101111111110111111011100111011110110111010111111011111111101111110111011011101111111111011111111111111011110100100011111111110110101110110001101111101110111101011110101111111101010010011110010111101110110011000011110111011111111110111111110111101111001000011101001111111110011101101011111001011101001111011011011111011111111111111111011101110000111011111011110111110101001111001110111111111111111101101111110110101111111001111001101101100111111100111101101010111111000101111101101111111110110111111111011111011011111111011111101011101010010011101101101110111110011111000011101011110011111011100100111011011010111111111010010101111111100111101110111010111110001101011111110101110111111101111111111111010110100011101101110011011010011110100111111111101111101110101111010011110001110101101101011101110001111111111011101001110110001101111110111011100011010101111110111111111011110110111111101110001011011110000111111001111100011111011011111101110010001001010011111110110011011011111000011101111101110111111111101111111100111110100010111011000110111101100110111111111100111101111101111101111111100011111101111111011100101110110110010011110111110110111011101110101001110110111101011111011011100101111111110100110100111100110111010011001111010110010111011111111010101100111111110111001000111110111011111101110001110101111111111010001101101111111101110111111101011011111011110111111011101110101011011011101101111101110101111101011011110101111111111111010011110011111111111101011110111111110001011111111111111111000111101111110110111111101111111111111001011011110111111111111011111110101111101110101111011010010101010110010011111011011110111011
1-(4) PubChem Fingerprint#
The PubChem Fingerprint is a 881-bit-long binary fingerprint (ftp://ftp.ncbi.nlm.nih.gov/pubchem/specifications/pubchem_fingerprints.pdf). Similar to the MACCS keys, it uses a pre-defined fragment dictionary. The PubChem fingerprint for each compound in PubChem can be downloaded from PubChem. However, because they are base64-encoded, they should be decoded into binary bitstrings or bitvectors.
Details about how to decode base64-encoded PubChem fingerprints is described on the last page of the PubChem Fingerprint specification (ftp://ftp.ncbi.nlm.nih.gov/pubchem/specifications/pubchem_fingerprints.pdf). Below is a user-defined function that decodes a PubChem fingerprint into a bit string.
from base64 import b64decode
def PCFP_BitString(pcfp_base64) :
pcfp_bitstring = "".join( ["{:08b}".format(x) for x in b64decode( pcfp_base64 )] )[32:913]
return pcfp_bitstring
pcfps = 'AAADcYBgAAAAAAAAAAAAAAAAAAAAAAAAAAAwAAAAAAAAAAABAAAAGAAAAAAACACAEAAwAIAAAACAACBCAAACAAAgAAAIiAAAAIgIICKAERCAIAAggAAIiAcAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA=='
print( len(PCFP_BitString(pcfps)) )
print(PCFP_BitString(pcfps))
881
10000000011000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001100000000000000000000000000000000000000000000000000000000000000000000000000010000000000000000000000000001100000000000000000000000000000000000000000000000100000000000100000000001000000000000001100000000000010000000000000000000000000000000100000000000000000100000010000100000000000000000000000100000000000000000001000000000000000000000000010001000100000000000000000000000000010001000000010000010000000100010100000000001000100010000100000000010000000000000001000001000000000000000000010001000100000000111000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
The generated bitstring can be converted to a bitvector that can be used for molecular similarity computation in RDKit (to be discussed in the next section).
from rdkit import DataStructs
bitvect = DataStructs.CreateFromBitString(PCFP_BitString(pcfps))
type(bitvect)
rdkit.DataStructs.cDataStructs.ExplicitBitVect
2. Computation of similarity scores#
import requests
import time
cids = [ 54454, # Simvastatin (Zocor)
54687, # Pravastatin (Pravachol)
60823, # Atorvastatin (Lipitor)
446155, # Fluvastatin (Lescol)
446157, # Rosuvastatin (Crestor)
5282452, # Pitavastatin (Livalo)
97938126 ] # Lovastatin (Altoprev)
Let’s get the SMILES strings from PubChem, generate Mol objects from them, and draw their chemical structures.
prolog = "https://pubchem.ncbi.nlm.nih.gov/rest/pug"
str_cid = ",".join([ str(x) for x in cids])
url = prolog + "/compound/cid/" + str_cid + "/property/isomericsmiles/txt"
res = requests.get(url)
smiles = res.text.split()
from rdkit import Chem
from rdkit.Chem import Draw
mols = [ Chem.MolFromSmiles(x) for x in smiles ]
Chem.Draw.MolsToGridImage(mols, molsPerRow=4, subImgSize=(200,200), legends=[str(x) for x in cids] )
[10:36:27] SMILES Parse Error: syntax error while parsing: Status:
[10:36:27] SMILES Parse Error: check for mistakes around position 2:
[10:36:27] Status:
[10:36:27] ~^
[10:36:27] SMILES Parse Error: Failed parsing SMILES 'Status:' for input: 'Status:'
[10:36:27] SMILES Parse Error: syntax error while parsing: 400
[10:36:27] SMILES Parse Error: check for mistakes around position 1:
[10:36:27] 400
[10:36:27] ^
[10:36:27] SMILES Parse Error: Failed parsing SMILES '400' for input: '400'
[10:36:27] SMILES Parse Error: syntax error while parsing: Code:
[10:36:27] SMILES Parse Error: check for mistakes around position 3:
[10:36:27] Code:
[10:36:27] ~~^
[10:36:27] SMILES Parse Error: Failed parsing SMILES 'Code:' for input: 'Code:'
[10:36:27] SMILES Parse Error: syntax error while parsing: PUGREST.BadRequest
[10:36:27] SMILES Parse Error: check for mistakes around position 2:
[10:36:27] PUGREST.BadRequest
[10:36:27] ~^
[10:36:27] SMILES Parse Error: Failed parsing SMILES 'PUGREST.BadRequest' for input: 'PUGREST.BadRequest'
[10:36:27] SMILES Parse Error: syntax error while parsing: Message:
[10:36:27] SMILES Parse Error: check for mistakes around position 1:
[10:36:27] Message:
[10:36:27] ^
[10:36:27] SMILES Parse Error: Failed parsing SMILES 'Message:' for input: 'Message:'
[10:36:27] SMILES Parse Error: syntax error while parsing: Invalid
[10:36:27] SMILES Parse Error: check for mistakes around position 3:
[10:36:27] Invalid
[10:36:27] ~~^
[10:36:27] SMILES Parse Error: Failed parsing SMILES 'Invalid' for input: 'Invalid'
[10:36:27] SMILES Parse Error: syntax error while parsing: property
[10:36:27] SMILES Parse Error: check for mistakes around position 2:
[10:36:27] property
[10:36:27] ~^
[10:36:27] SMILES Parse Error: Failed parsing SMILES 'property' for input: 'property'
Now generate MACCS keys for each compound.
from rdkit import DataStructs
from rdkit.Chem import MACCSkeys
fps = [ MACCSkeys.GenMACCSKeys(x) for x in mols ]
---------------------------------------------------------------------------
ArgumentError Traceback (most recent call last)
Cell In[21], line 4
1 from rdkit import DataStructs
2 from rdkit.Chem import MACCSkeys
----> 4 fps = [ MACCSkeys.GenMACCSKeys(x) for x in mols ]
Cell In[21], line 4, in <listcomp>(.0)
1 from rdkit import DataStructs
2 from rdkit.Chem import MACCSkeys
----> 4 fps = [ MACCSkeys.GenMACCSKeys(x) for x in mols ]
ArgumentError: Python argument types in
rdkit.Chem.rdMolDescriptors.GetMACCSKeysFingerprint(NoneType)
did not match C++ signature:
GetMACCSKeysFingerprint(RDKit::ROMol mol)
Now let’s compute the pair-wise similarity scores among them. To make higher scores easier to find, they are indicated with the “*” character(s).
for i in range(0, len(fps)) :
for j in range(i+1, len(fps)) :
score = DataStructs.FingerprintSimilarity(fps[i], fps[j])
print(cids[i], "vs.", cids[j], ":", round(score,3), end='')
if ( score >= 0.85 ):
print(" ****")
elif ( score >= 0.75 ):
print(" ***")
elif ( score >= 0.65 ):
print(" **")
elif ( score >= 0.55 ):
print(" *")
else:
print(" ")
54454 vs. 54687 : 0.812 ***
54454 vs. 60823 : 0.354
54454 vs. 446155 : 0.379
54454 vs. 446157 : 0.307
54454 vs. 5282452 : 0.4
54454 vs. 97938126 : 0.867 ****
54687 vs. 60823 : 0.387
54687 vs. 446155 : 0.397
54687 vs. 446157 : 0.287
54687 vs. 5282452 : 0.421
54687 vs. 97938126 : 0.8 ***
60823 vs. 446155 : 0.662 **
60823 vs. 446157 : 0.535
60823 vs. 5282452 : 0.507
60823 vs. 97938126 : 0.365
446155 vs. 446157 : 0.526
446155 vs. 5282452 : 0.735 **
446155 vs. 97938126 : 0.393
446157 vs. 5282452 : 0.473
446157 vs. 97938126 : 0.298
5282452 vs. 97938126 : 0.444
By default, the similarity score is generated using the Tanimoto equation. RDKit also supports other similarity metrics, including Dice, Cosine, Sokal, Russel, Kulczynski, McConnaughey, and Tversky. The definition of these metrics is available at the LibreTexts page (https://bit.ly/2kx9NCd).
print("Tanimoto :", round(DataStructs.TanimotoSimilarity(fps[0], fps[1]), 4))
print("Dice :", round(DataStructs.DiceSimilarity(fps[0], fps[1]), 4))
print("Cosine :", round(DataStructs.CosineSimilarity(fps[0], fps[1]), 4))
print("Sokal :", round(DataStructs.SokalSimilarity(fps[0], fps[1]), 4))
print("McConnaughey:", round(DataStructs.McConnaugheySimilarity(fps[0], fps[1]), 4))
Tanimoto : 0.8125
Dice : 0.8966
Cosine : 0.8971
Sokal : 0.6842
McConnaughey: 0.7952
The Tversky score is an asymmetric similarity measure, and its computation requires the weightings of the two molecules being compared.
for i in range(0,11) :
alpha = round( i * 0.1, 1 )
beta = round( 1 - alpha, 1 )
print("(alpha, beta) = (", alpha, ",", beta, ") : ", end='')
print(round(DataStructs.TverskySimilarity( fps[0], fps[1], alpha, beta ), 4))
(alpha, beta) = ( 0.0 , 1.0 ) : 0.9286
(alpha, beta) = ( 0.1 , 0.9 ) : 0.922
(alpha, beta) = ( 0.2 , 0.8 ) : 0.9155
(alpha, beta) = ( 0.3 , 0.7 ) : 0.9091
(alpha, beta) = ( 0.4 , 0.6 ) : 0.9028
(alpha, beta) = ( 0.5 , 0.5 ) : 0.8966
(alpha, beta) = ( 0.6 , 0.4 ) : 0.8904
(alpha, beta) = ( 0.7 , 0.3 ) : 0.8844
(alpha, beta) = ( 0.8 , 0.2 ) : 0.8784
(alpha, beta) = ( 0.9 , 0.1 ) : 0.8725
(alpha, beta) = ( 1.0 , 0.0 ) : 0.8667
Exercise 2a: Compute the Tanimoto similarity scores between the seven compounds used in this section, using the PubChem fingerprints
Download the PubChem Fingerprint for the seven CIDs.
Convert the downloaded fingerprints into bit vectors.
Compute the pair-wise Tanimoto scores using the bit vectors.
# Write your code in this cell
3. Interpretation of similarity scores#
Using molecular fingeprints. we can compute the similarity scores between molecules. However, how should these scores be interpreted? For example, the Tanimoto score between CID 60823 and CID 446155 is computed to be 0.662, but does it mean that the two compounds are similar? How similar is similar? The following analysis would help answer these questions.
Step 1. Randomly select 1,000 compounds from PubChem and download their SMILES strings.
import random
random.seed(0)
cid_max = 138962044 # The maximum CID in PubChem as of September 2019
cids = []
for x in range(1000):
cids.append(random.randint(1, cid_max + 1))
chunk_size = 100
if len(cids) % chunk_size == 0 :
num_chunks = int( len(cids) / chunk_size )
else :
num_chunks = int( len(cids) / chunk_size ) + 1
smiles = []
for i in range(num_chunks):
if (i == 0):
print("Processing chunk ", end='')
print(i, end=' ')
idx1 = chunk_size * i
idx2 = chunk_size * (i + 1)
str_cids = ",".join([ str(x) for x in cids[idx1:idx2]])
url = prolog + "/compound/cid/" + str_cids + "/property/isomericsmiles/txt"
res = requests.get(url)
if ( res.status_code == 200) :
smiles.extend( res.text.split() )
else :
print("Chunk", i, "Failed to get SMILES.")
time.sleep(0.2)
print("Done!")
print("# Number of SMILES : ", len(smiles))
Processing chunk 0 1 2 3 4 5 6 7 8 9 Done!
# Number of SMILES : 1000
Step 2. Generate the MACCSKeys for each compound.
from rdkit import Chem
mols = [ Chem.MolFromSmiles(x) for x in smiles if x != None ]
fps = [ MACCSkeys.GenMACCSKeys(x) for x in mols if x != None ]
print("# Number of compounds:", len(mols))
print("# Number of fingerprints:", len(fps))
# Number of compounds: 1000
# Number of fingerprints: 1000
# Run this cell if the number of compounds != the number of fingerprints.
#if ( len(cids) != len(fps) ):
# print("SMILES at index", mols.index(None), ":", smiles[ mols.index(None) ])
Step 3. Compute the Tanimoto scores between compounds.
print("# The number of compound pairs:", (len(fps) * (len(fps) - 1))/2 )
# The number of compound pairs: 499500.0
scores = []
for i in range(0, len(fps)) :
if (i == 0) :
print("Processing compound ", end='')
if (i % 100 == 0) :
print(i, end=' ')
for j in range(i+1, len(fps)) :
scores.append(DataStructs.FingerprintSimilarity(fps[i], fps[j]))
print("Done!")
print("# Number of scores : ", len(scores))
Processing compound 0 100 200 300 400 500 600 700 800 900 Done!
# Number of scores : 499500
Step 4. Generate a histogram that shows the distribution of the pair-wise scores.
import matplotlib.pyplot as plt
%matplotlib inline
mybins = [ x * 0.01 for x in range(101)]
fig = plt.figure(figsize=(8,4), dpi=300)
plt.subplot(1, 2, 1)
plt.title("Distribution")
plt.hist(scores, bins=mybins)
plt.subplot(1, 2, 2)
plt.title("Cumulative Distribution")
plt.hist(scores, bins=mybins, density=True, cumulative=1)
plt.plot([0,1],[0.95,0.95]);
for i in range(21) :
thresh = i / 20
num_similar_pairs = len([x for x in scores if x >= thresh])
prob = num_similar_pairs / len(scores) * 100
print("%.3f %8d (%8.4f %%)" % (thresh, num_similar_pairs, round(prob,4)))
0.000 499500 (100.0000 %)
0.050 497010 ( 99.5015 %)
0.100 488714 ( 97.8406 %)
0.150 469193 ( 93.9325 %)
0.200 435672 ( 87.2216 %)
0.250 385198 ( 77.1167 %)
0.300 318637 ( 63.7912 %)
0.350 245621 ( 49.1734 %)
0.400 175158 ( 35.0667 %)
0.450 111668 ( 22.3560 %)
0.500 66599 ( 13.3331 %)
0.550 32209 ( 6.4482 %)
0.600 13897 ( 2.7822 %)
0.650 4668 ( 0.9345 %)
0.700 1351 ( 0.2705 %)
0.750 355 ( 0.0711 %)
0.800 76 ( 0.0152 %)
0.850 24 ( 0.0048 %)
0.900 6 ( 0.0012 %)
0.950 1 ( 0.0002 %)
1.000 0 ( 0.0000 %)
print("Average:", sum(scores)/len(scores))
Average: 0.3488195436976387
From the distribution of the similarity scores among 1,000 compounds, we observe the following:
If you randomly select two compounds from PubChem, the similarity score between them (computed using the Tanimoto equation and MACCS keys) is ~0.35 on average.
About %5 of randomly selected compound pairs have a similarity score greater than 0.55.
About %1 of randomly selected compound pairs have a similarity score greater than 0.65.
If two compounds have a Tanimoto score of 0.35, it is close to the avaerage Tanimoto score between randomly selected compounds and there is a 50% chance that you will get a score of 0.35 or greater just by selecting two compounds from PubChem. Therefore, it is reasonable to consider the two compounds are not similar.
The Tanimoto index may have a value ranging from 0 (for no similarity) to 1 (for identical molecules) and the midpoint of this value range is 0.5. Because of this, a Tanimoto score of 0.55 may not sound great enough to consider two compounds to be similar. However, according to the score distribution curve generated here, only ~5% of randomly selected compound pairs will have a score greater than this.
In the previous section, we computed the similarity scores between some cholesterol-lowering drugs, and CID 60823 and CID 446155 had a Tanimoto score of 0.662. Based on the score distribution curve generated in the second section, we can say that the probablilty of two randomly selected compounds from PubChem having a Tanimoto score greater than 0.662 is less than 1%.
The following code cell demonstrates how to find an appropriate similarity score threshold above which a given percentage of the compound pairs will be considered to be similar to each other.
scores.sort() # Sort the scores in an increasing order.
# to find a threshold for top 3% compound pairs (i.e., 97% percentile)
print("# total compound pairs: ", len(scores))
print("# 95% of compound pairs: ", len(scores) * 0.97)
print("# score at 95% percentile:", scores[ round(len(scores) * 0.97) ] )
# total compound pairs: 499500
# 95% of compound pairs: 484515.0
# score at 95% percentile: 0.5945945945945946
Exercise 3a: In this exercise, we want to generate the distribution of the similarity scores among 1,000 compounds randomly selected from PubChem, using different molecular fingeprints and similarity metrics.
For molecular fingerprints, use the following:
PubChem Fingerprint
MACCS keys
Morgan Fingerprint (ECFP4 analogue, 1024-bit-long)
For similarity metrics, use the following:
Tanimoto similarity
Dice similarity
Cosine similarity
As a result, a total of 9 distribution curves need to be generated.
Here are additional instructions to follow:
When generating the histograms, bin the scores from 0 to 1 with an increment of 0.01.
For each distribution curve, determine the similarity score threshold so that 1% of the compound pairs have a similarity score greater than or equal to this threshold.
Use RDKit to generate the MACCS keys and Morgan fingerprint and download the PubChem fingerprints from PubChem.
For reproducibility, use random.seed(2019) before you generate random CIDs.
Step 1: Generate 1,000 random CIDs, download the isomeric SMILES for them, and create the RDKit mol objects from the downloaded SMILES strings.
# Write your code in this cell
Step 2: Generate the fingerprints, compute the similarity scores, determine similarity thresholds, and make histograms.
# Write your code in this cell