#!/usr/bin/env python
# coding: utf-8

# # Thresholds for "random" in fingerprints the RDKit supports
# 
# This is an updated version of a post. The original version of the notebook can be [found in github](https://github.com/greglandrum/rdkit_blog/blob/936b35b77ddca3e843991d1c4b1e3532e754b734/notebooks/Fingerprint%20Thresholds.ipynb).
# 
# A frequent question that comes up when considering fingerprint similarity is: "What threshold should I use to determine what a neighbor is?" The answer is poorly defined. Of course it depends heavily on the details of the fingerprint, but there's also a very subjective component: you want to pick a low enough threshold that you're sure you won't miss anything, but you don't want to pick up too much noise.
# 
# The goal here is to systematically come up with some guidelines that can be used for fingerprints supported within the RDKit. We will do that by looking a similarities between random "drug-like" (MW<600) molecules picked from ChEMBL.
# 
# For the analysis, the 25K similarity values are sorted and the values at particular threshold are examined. 
# 
# There's a fair amount of code and results below, so here's the summary table. To help interpret this: 22500 of the 25000 pairs (90%) have a MACCS keys similarity value less than 0.528. 
# 
# 
# <table>
# <tr><th>Fingerprint</th><th>Metric</th><th>70% level</th><th>80% level</th><th>90% level</th><th>95% level</th><th>99% level</th></tr>
# <tr><td>MACCS</td><td>Tanimoto</td>
#   <td>0.431</td>
#   <td>0.471</td>
#   <td>0.528</td>
#   <td>0.575</td>
#   <td>0.655</td>
# </tr><tr><td>Morgan0 (counts)</td><td>Tanimoto</td>
#   <td>0.429</td>
#   <td>0.471</td>
#   <td>0.525</td>
#   <td>0.568</td>
#   <td>0.651</td>
# </tr><tr><td>Morgan1 (counts)</td><td>Tanimoto</td>
#   <td>0.265</td>
#   <td>0.293</td>
#   <td>0.333</td>
#   <td>0.364</td>
#   <td>0.429</td>
# </tr><tr><td>Morgan2 (counts)</td><td>Tanimoto</td>
#   <td>0.181</td>
#   <td>0.201</td>
#   <td>0.229</td>
#   <td>0.252</td>
#   <td>0.305</td>
# </tr><tr><td>Morgan3 (counts)</td><td>Tanimoto</td>
#   <td>0.141</td>
#   <td>0.156</td>
#   <td>0.178</td>
#   <td>0.196</td>
#   <td>0.238</td>
# </tr><tr><td>Morgan0 (bits)</td><td>Tanimoto</td>
#   <td>0.435</td>
#   <td>0.475</td>
#   <td>0.529</td>
#   <td>0.571</td>
#   <td>0.656</td>
# </tr><tr><td>Morgan1 (bits)</td><td>Tanimoto</td>
#   <td>0.273</td>
#   <td>0.301</td>
#   <td>0.341</td>
#   <td>0.371</td>
#   <td>0.434</td>
# </tr><tr><td>Morgan2 (bits)</td><td>Tanimoto</td>
#   <td>0.197</td>
#   <td>0.217</td>
#   <td>0.246</td>
#   <td>0.269</td>
#   <td>0.322</td>
# </tr><tr><td>Morgan3 (bits)</td><td>Tanimoto</td>
#   <td>0.165</td>
#   <td>0.181</td>
#   <td>0.203</td>
#   <td>0.222</td>
#   <td>0.264</td>
# </tr><tr><td>FeatMorgan0 (counts)</td><td>Tanimoto</td>
#   <td>0.583</td>
#   <td>0.630</td>
#   <td>0.690</td>
#   <td>0.737</td>
#   <td>0.818</td>
# </tr><tr><td>FeatMorgan1 (counts)</td><td>Tanimoto</td>
#   <td>0.390</td>
#   <td>0.425</td>
#   <td>0.474</td>
#   <td>0.511</td>
#   <td>0.581</td>
# </tr><tr><td>FeatMorgan2 (counts)</td><td>Tanimoto</td>
#   <td>0.272</td>
#   <td>0.298</td>
#   <td>0.333</td>
#   <td>0.364</td>
#   <td>0.424</td>
# </tr><tr><td>FeatMorgan3 (counts)</td><td>Tanimoto</td>
#   <td>0.209</td>
#   <td>0.228</td>
#   <td>0.256</td>
#   <td>0.279</td>
#   <td>0.328</td>
# </tr><tr><td>FeatMorgan0 (bits)</td><td>Tanimoto</td>
#   <td>0.583</td>
#   <td>0.630</td>
#   <td>0.690</td>
#   <td>0.737</td>
#   <td>0.818</td>
# </tr><tr><td>FeatMorgan1 (bits)</td><td>Tanimoto</td>
#   <td>0.395</td>
#   <td>0.429</td>
#   <td>0.477</td>
#   <td>0.514</td>
#   <td>0.585</td>
# </tr><tr><td>FeatMorgan2 (bits)</td><td>Tanimoto</td>
#   <td>0.284</td>
#   <td>0.310</td>
#   <td>0.347</td>
#   <td>0.376</td>
#   <td>0.434</td>
# </tr><tr><td>FeatMorgan3 (bits)</td><td>Tanimoto</td>
#   <td>0.228</td>
#   <td>0.248</td>
#   <td>0.276</td>
#   <td>0.299</td>
#   <td>0.349</td>
# </tr><tr><td>RDKit 4 (bits)</td><td>Tanimoto</td>
#   <td>0.209</td>
#   <td>0.239</td>
#   <td>0.285</td>
#   <td>0.325</td>
#   <td>0.426</td>
# </tr><tr><td>RDKit 5 (bits)</td><td>Tanimoto</td>
#   <td>0.197</td>
#   <td>0.219</td>
#   <td>0.253</td>
#   <td>0.287</td>
#   <td>0.368</td>
# </tr><tr><td>RDKit 6 (bits)</td><td>Tanimoto</td>
#   <td>0.230</td>
#   <td>0.250</td>
#   <td>0.280</td>
#   <td>0.308</td>
#   <td>0.369</td>
# </tr><tr><td>RDKit 7 (bits)</td><td>Tanimoto</td>
#   <td>0.313</td>
#   <td>0.346</td>
#   <td>0.389</td>
#   <td>0.429</td>
#   <td>0.507</td>
# </tr><tr><td>linear RDKit 4 (bits)</td><td>Tanimoto</td>
#   <td>0.225</td>
#   <td>0.258</td>
#   <td>0.309</td>
#   <td>0.354</td>
#   <td>0.462</td>
# </tr><tr><td>linear RDKit 5 (bits)</td><td>Tanimoto</td>
#   <td>0.198</td>
#   <td>0.225</td>
#   <td>0.269</td>
#   <td>0.309</td>
#   <td>0.404</td>
# </tr><tr><td>linear RDKit 6 (bits)</td><td>Tanimoto</td>
#   <td>0.187</td>
#   <td>0.210</td>
#   <td>0.246</td>
#   <td>0.282</td>
#   <td>0.365</td>
# </tr><tr><td>linear RDKit 7 (bits)</td><td>Tanimoto</td>
#   <td>0.182</td>
#   <td>0.203</td>
#   <td>0.234</td>
#   <td>0.264</td>
#   <td>0.337</td>
# </tr><tr><td>Atom Pairs (counts)</td><td>Tanimoto</td>
#   <td>0.180</td>
#   <td>0.204</td>
#   <td>0.237</td>
#   <td>0.265</td>
#   <td>0.325</td>
# </tr><tr><td>Torsions (counts)</td><td>Tanimoto</td>
#   <td>0.107</td>
#   <td>0.130</td>
#   <td>0.165</td>
#   <td>0.194</td>
#   <td>0.266</td>
# </tr><tr><td>Atom Pairs (bits)</td><td>Tanimoto</td>
#   <td>0.275</td>
#   <td>0.301</td>
#   <td>0.335</td>
#   <td>0.363</td>
#   <td>0.415</td>
# </tr><tr><td>Torsions (bits)</td><td>Tanimoto</td>
#   <td>0.133</td>
#   <td>0.155</td>
#   <td>0.188</td>
#   <td>0.219</td>
#   <td>0.288</td>
# </tr><tr><td>Avalon 512 (bits)</td><td>Tanimoto</td>
#   <td>0.369</td>
#   <td>0.407</td>
#   <td>0.461</td>
#   <td>0.505</td>
#   <td>0.575</td>
# </tr><tr><td>Avalon 1024 (bits)</td><td>Tanimoto</td>
#   <td>0.269</td>
#   <td>0.297</td>
#   <td>0.340</td>
#   <td>0.375</td>
#   <td>0.449</td>
# </tr><tr><td>Avalon 512 (counts)</td><td>Tanimoto</td>
#   <td>0.300</td>
#   <td>0.333</td>
#   <td>0.379</td>
#   <td>0.418</td>
#   <td>0.491</td>
# </tr><tr><td>Avalon 1024 (counts)</td><td>Tanimoto</td>
#   <td>0.267</td>
#   <td>0.299</td>
#   <td>0.344</td>
#   <td>0.384</td>
#   <td>0.462</td>
# </tr>
# </table>
# 

# In[1]:


from rdkit import Chem
from rdkit.Chem import rdMolDescriptors
from rdkit.Avalon import pyAvalonTools
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
from rdkit import rdBase
from rdkit import DataStructs
from collections import defaultdict
import pickle,random,gzip
print(rdBase.rdkitVersion)
import time
print(time.asctime())
get_ipython().run_line_magic('pylab', 'inline')


# # Read in the data
# 
# We're using the set of 25K reference pairs generated in an earlier post: http://rdkit.blogspot.ch/2013/10/building-similarity-comparison-set-goal.html
# 
# As a quick reminder: these are pairs of molecules taken from ChEMBL with MW<600 and a count-based MFP0 similarity of at least 0.7 to each other.

# In[2]:


ind = [x.split() for x in gzip.open('../data/chembl16_25K.pairs.txt.gz')]
ms1 = []
ms2 = []
for i,row in enumerate(ind):
    m1 = Chem.MolFromSmiles(row[1])
    ms1.append((row[0],m1))
    m2 = Chem.MolFromSmiles(row[3])
    ms2.append((row[2],m2))
    


# Those pairs are related to each other, but we want random pairs, so shuffle the second list:

# In[3]:


random.seed(23)
random.shuffle(ms2)


# In[4]:


try:
    import ipyparallel as ipp
    rc = ipp.Client()
    dview = rc[:]
    dview.execute('from rdkit import Chem')
    dview.execute('from rdkit import Descriptors')
    dview.execute('from rdkit.Chem import rdMolDescriptors')
    dview.execute('from rdkit.Avalon import pyAvalonTools')
except:
    print("could not use ipyparallel")
    dview = None

results_accum = dict()
def compareFPs(ms1,ms2,fpfn,fpName):
    if dview is not None:
        fps = dview.map_sync(lambda x:fpfn(x[1]),ms1)
        fp2s = dview.map_sync(lambda x:fpfn(x[1]),ms2)
    else:
        fps = [fpfn(x[1]) for x in ms1]
        fp2s = [fpfn(x[1]) for x in ms2]
    sims = [DataStructs.TanimotoSimilarity(x,y) for x,y in zip(fps,fp2s)]
    sl = sorted(sims)
    np = len(sl)
    with open('fp_results.txt','a+') as outf:
        outf.write(f'<tr><td>{fpName}</td><td>Tanimoto</td>\n')
        accum = {}
        for bin in (.7,.8,.9,.95,.99):
            simv = sl[int(bin*np)]
            print( bin,simv)
            outf.write(f'  <td>{simv:.3f}</td>\n')
            accum[bin] = simv
        outf.write('</tr>')
        results_accum[fpName] = accum
    hist(sims,bins=20)
    xlabel(fpName)
    


# # MACCS

# In[5]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMACCSKeysFingerprint(x),"MACCS")


# # Morgan FPs

# ## count based

# In[6]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,0),"Morgan0 (counts)")


# In[7]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,1),"Morgan1 (counts)")


# In[8]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,2),"Morgan2 (counts)")


# In[9]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,3),"Morgan3 (counts)")


# ## bit-vector based

# In[10]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,0,1024),"Morgan0 (bits)")


# In[11]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,1,1024),"Morgan1 (bits)")


# In[12]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,2,1024),"Morgan2 (bits)")


# In[13]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,3,1024),"Morgan3 (bits)")


# # FeatMorgan

# ## count based

# In[14]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,0,useFeatures=True),"FeatMorgan0 (counts)")


# In[15]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,1,useFeatures=True),"FeatMorgan1 (counts)")


# In[16]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,2,useFeatures=True),"FeatMorgan2 (counts)")


# In[17]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetMorganFingerprint(x,3,useFeatures=True),"FeatMorgan3 (counts)")


# ## bit vectors

# In[18]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,0,1024,useFeatures=True),"FeatMorgan0 (bits)")


# In[19]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,1,1024,useFeatures=True),"FeatMorgan1 (bits)")


# In[20]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,2,1024,useFeatures=True),"FeatMorgan2 (bits)")


# In[21]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedMorganFingerprint(x,3,1024,useFeatures=True),"FeatMorgan3 (bits)")


# # RDKit

# ## Branched (default)

# In[22]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=4),"RDKit 4 (bits)")


# In[23]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=5),"RDKit 5 (bits)")


# In[24]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=6),"RDKit 6 (bits)")


# In[25]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=7),"RDKit 7 (bits)")


# ## linear

# In[26]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=4,branchedPaths=False),"linear RDKit 4 (bits)")


# In[27]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=5,branchedPaths=False),"linear RDKit 5 (bits)")


# In[28]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=6,branchedPaths=False),"linear RDKit 6 (bits)")


# In[29]:


compareFPs(ms1,ms2,lambda x:Chem.RDKFingerprint(x,maxPath=7,branchedPaths=False),"linear RDKit 7 (bits)")


# # Atom pairs and torsions

# ## count-based

# In[30]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetAtomPairFingerprint(x),"Atom Pairs (counts)")


# In[31]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetTopologicalTorsionFingerprint(x),"Topological Torsions (counts)")


# ## bit vectors

# In[32]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedAtomPairFingerprintAsBitVect(x),"Atom Pairs (bits)")


# In[33]:


compareFPs(ms1,ms2,lambda x:rdMolDescriptors.GetHashedTopologicalTorsionFingerprintAsBitVect(x),"Topological Torsions (bits)")


# # Avalon

# In[34]:


compareFPs(ms1,ms2,lambda x:pyAvalonTools.GetAvalonFP(x,512),"Avalon 512 (bits)")


# In[35]:


compareFPs(ms1,ms2,lambda x:pyAvalonTools.GetAvalonFP(x,1024),"Avalon 1024 (bits)")


# # Avalon Counts

# In[36]:


compareFPs(ms1,ms2,lambda x:pyAvalonTools.GetAvalonCountFP(x,512),"Avalon 512 (counts)")


# In[37]:


compareFPs(ms1,ms2,lambda x:pyAvalonTools.GetAvalonCountFP(x,1024),"Avalon 1024 (counts)")


# In[38]:


pickle.dump(results_accum,open('./results/fp_thresholds_random_accum.pkl','wb+'))