In [1]:
import sys,os
sys.path.insert(0,'/global/project/projectdirs/metatlas/anaconda/lib/python2.7/site-packages' )
from metatlas import metatlas_objects as metob
import csv
curr_ld_lib_path = ''
os.environ['LD_LIBRARY_PATH'] = curr_ld_lib_path + ':/project/projectdirs/openmsi/jupyterhub_libs/boost_1_55_0/lib' + ':/project/projectdirs/openmsi/jupyterhub_libs/lib'
import sys
# sys.path.remove('/anaconda/lib/python2.7/site-packages')
sys.path.append('/global/project/projectdirs/openmsi/jupyterhub_libs/anaconda/lib/python2.7/site-packages')
sys.path.insert(0,'/project/projectdirs/openmsi/projects/meta-iq/pactolus/pactolus' )
from rdkit import Chem
# from rdkit.Chem.rdMolDescriptors import ExactMolWt
from rdkit.Chem import Descriptors
from rdkit.Chem import rdMolDescriptors
from rdkit.Chem import AllChem
from rdkit.Chem.Fingerprints import FingerprintMols
from rdkit.Chem import Draw
from rdkit import DataStructs
import numpy as np
from matplotlib import pylab as plt
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import networkx as nx
In [2]:
allCompounds = metob.retrieve('Compounds')
print len(allCompounds)
In [11]:
for c in allCompounds:
if 'AMP' in c.name:
print c
In [3]:
allCompounds[0]
Out[3]:
In [44]:
#Store the necessary data from the metatlas database for each compound
myNames = [x.name for x in allCompounds]
myMW = [x.MonoIsotopic_molecular_weight for x in allCompounds]
myformula = [x.formula for x in allCompounds]
myID = [x.unique_id for x in allCompounds]
# Convert each InChI to an RDKit molecule
myMols = [Chem.MolFromInchi(str(x.InChI)) for x in allCompounds]
# Calculate Chemical Fingerprints
# fps = [FingerprintMols.FingerprintMol(x) for x in myMols]
# Calculate Chemical Fingerprints Using Fixed BitLength
N = 2048
fps = [FingerprintMols.FingerprintMol(x,minPath=1,maxPath=7,fpSize=N,bitsPerHash=2,useHs=True,tgtDensity=0,minSize=N,branchedPaths=True,
useBondOrder=True,atomInvariants=[],fromAtoms=[],atomBits=[]) for x in myMols]
In [19]:
fp_mat = np.zeros((len(fps),N))
for i,f in enumerate(fps):
for j in range(N):
fp_mat[i,j] = float(DataStructs.BitVectToText(fps[i])[j])
In [20]:
print fp_mat.shape
In [24]:
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import linkage
from scipy.cluster.hierarchy import fcluster
dists = pdist(fp_mat[:,:],'cityblock')
In [26]:
from scipy.spatial.distance import squareform
pd = squareform(dists)
In [25]:
# Y = linkage(dists, method='average')
# cutoff = 0
# metClusterIDs = fcluster(Y,cutoff,criterion='distance')
# print len(np.unique(metClusterIDs))
In [ ]:
from rdkit.DataManip import Metric
import time
start_time = time.time()
d = np.zeros((len(fps),len(fps)))
for i,fps1 in enumerate(fps):
for j,fps2 in enumerate(fps):
if i<j:
pd[i,j] = DataStructs.FingerprintSimilarity(fps[i],fps[j])
### TO
# there is something that doesn't make sense with this.
# probably row column ordering is different than scipy
# d = Metric.GetTanimotoSimMat(fps)
# pd = squareform(d)
e_time = time.time() - start_time
print e_time
In [ ]:
In [18]:
idx1 = 2945
idx2 = 168
test_d = DataStructs.FingerprintSimilarity(fps[idx1],fps[idx2])
print myNames[idx1]," ",myNames[idx2], test_d
print pd[idx1,idx2]
In [46]:
# s = np.sum(pd==0.0,1)
# idx1 = np.argsort(s)[::-1]
# checker = 0+18
# print idx1[checker], myNames[idx1[checker]]
# idx2 = np.argwhere(pd[idx1[checker],:]==0.0)
# for i,idx in enumerate(idx2):
# print i, idx, myNames[idx], pd[checker,idx],pd[idx1[checker],idx[0]]
# # idx2 = np.argwhere(pd[i,:]>=1)
# # for j in idx2:
# # print i,j,myNames[j],pd[i,j]
In [142]:
fid = open('MetAtlas_Network_unity.tab','w')
for i in range(pd.shape[0]):
for j in range(pd.shape[1]):
if i<j:
if pd[i,j]==1:
fid.write('"#%d %s"\t"#%d %s"\t%5.5f\n' % (i,
myNames[i],
j,
myNames[j],
pd[i,j]))
# print i, idx, myNames[idx], pd[checker,idx],pd[idx1[checker],idx[0]]
# pd_row = pd[i,:]
# m = np.max(pd_row)
# if m==1:
# idx = np.argwhere( pd_row == 1.0 )
# for j in idx:
# if j!=i:
# fid.write('"%d%s"\t"%d%s"\t%5.5f\n' % (i,myNames[i],j,myNames[j],pd_row[j]))
fid.close()
In [48]:
fid = open('MetAtlas_Network_Similar.tab','w')
for i in range(pd.shape[0]):
temp = pd[i,0:i]
if len(temp)>1:
m = min(temp)
if m<100000:
for j in range(i-1):
if pd[i,j]==m:
fid.write('"#%d %s %5.4f %s"\t"#%d %s %5.4f %s"\t%5.5f\n' % (i,
myformula[i],
myMW[i],
myNames[i],
j,
myformula[j],
myMW[j],
myNames[j],
pd[i,j]))
# m = max(pd[i,:])
# # if m == 1:
# for j in range(pd.shape[1]):
# if i<j:
# if pd[i,j]==m:
# fid.write('"#%d %s %5.4f %s"\t"#%d %s %5.4f %s"\t%5.5f\n' % (i,
# myformula[i],
# myMW[i],
# myNames[i],
# j,
# myformula[j],
# myMW[j],
# myNames[j],
# pd[i,j]))
# print i, idx, myNames[idx], pd[checker,idx],pd[idx1[checker],idx[0]]
# pd_row = pd[i,:]
# m = np.max(pd_row)
# if m==1:
# idx = np.argwhere( pd_row == 1.0 )
# for j in idx:
# if j!=i:
# fid.write('"%d%s"\t"%d%s"\t%5.5f\n' % (i,myNames[i],j,myNames[j],pd_row[j]))
fid.close()
In [45]:
fid = open('Tanimoto_2048_MetAtlas_Network_Cutoff_Min_Match_100.tab','w')
for i in range(pd.shape[0]):
temp = pd[i,0:i]
if len(temp)>1:
m = min(temp)
if m<100:
for j in range(i-1):
if pd[i,j]==m:
fid.write('"#%d %s %5.4f %s"\t"#%d %s %5.4f %s"\t%5.5f\n' % (i,
myformula[i],
myMW[i],
myNames[i],
j,
myformula[j],
myMW[j],
myNames[j],
pd[i,j]))
# print i, idx, myNames[idx], pd[checker,idx],pd[idx1[checker],idx[0]]
# pd_row = pd[i,:]
# m = np.max(pd_row)
# if m==1:
# idx = np.argwhere( pd_row == 1.0 )
# for j in idx:
# if j!=i:
# fid.write('"%d%s"\t"%d%s"\t%5.5f\n' % (i,myNames[i],j,myNames[j],pd_row[j]))
fid.close()
In [ ]:
In [138]:
print i,j
pd.shape[0]
idx2 = np.where(pd[0,:]==1.0)
print idx2
for i in idx2:
print i
print myNames[i]
print myNames[0]
In [54]:
%%javascript
var nb = IPython.notebook;
var kernel = IPython.notebook.kernel;
var command = "NOTEBOOK_FULL_PATH = '" + nb.base_url + nb.notebook_path + "'";
kernel.execute(command);
In [58]:
filename = os.path.basename(NOTEBOOK_FULL_PATH)
%system cp $filename /project/projectdirs/openmsi/www/
temp = '%s/%s'%('/project/projectdirs/openmsi/www',filename)
%system chmod 775 $temp
print 'http://nbviewer.ipython.org/url/portal.nersc.gov/project/openmsi/%s?flush_cache=true'%filename
In [ ]: