<jupyter_start><jupyter_text>Imports Colab Should take ~5-10 minIf are not running this in Colab, and you would like to, click this:and run the following cell once per kernel runtime reset, if resetting the kernel, no need to run.<jupyter_code>import os
import subprocess
import sys

GIT_NAME='melting_points_ml'
GIT_URL='https://beangoben:testPASS12@github.com/beangoben/{}.git'.format(GIT_NAME)
IN_COLAB = 'google.colab' in sys.modules

def run_cmd(cmd):
    print('Output of "{}":'.format(cmd))
    print(subprocess.run(cmd,stdout=subprocess.PIPE, shell=True).stdout.decode('utf-8'))
        
if IN_COLAB:
    SRC_DIR='.'
    run_cmd('sudo apt-get install git-lfs')
    run_cmd('git lfs install')
    run_cmd('rm -rf sample_data')
    run_cmd('rm -rf {}'.format(GIT_NAME))
    run_cmd('git clone --verbose --progress {}'.format(GIT_URL))
    run_cmd('mv {}/* . '.format(GIT_NAME))
    run_cmd('rm -rf {}'.format(GIT_NAME))
    sys.path.append(os.path.join(SRC_DIR,'code'))
else:
    SRC_DIR='..'
    
print('Using colab? {}, using root directory "{}"'.format(IN_COLAB,SRC_DIR))<jupyter_output>Using colab? False, using root directory ".."<jupyter_text>install extra packages in colab<jupyter_code>import sys
import os
if IN_COLAB:
    SRC_DIR='.'
    sys.path.append(os.path.join(SRC_DIR,'code'))
    import colab_helpers
    colab_helpers.install_rdkit()
else:
    SRC_DIR='..'<jupyter_output><empty_output><jupyter_text>if in colab, we need to add the custom anaconda env (with rdkit) to the path<jupyter_code>if IN_COLAB:
    sys.path.append('/usr/local/lib/python3.7/site-packages/')<jupyter_output><empty_output><jupyter_text>Regular imports<jupyter_code>import json
import os
import random
from collections import OrderedDict
# scientific python
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy as sp
import seaborn as sns
# Chemistry
import rdkit
from rdkit import Chem
from rdkit.Chem.Draw import MolsToGridImage, MolToImage
from rdkit import DataStructs
from rdkit.Chem.rdMolDescriptors import GetMorganFingerprint
# ML
import sklearn
from sklearn import preprocessing<jupyter_output><empty_output><jupyter_text>Repo specific imports<jupyter_code>sys.path.append(os.path.join(SRC_DIR, 'code'))
import utils
utils.plot_settings()<jupyter_output>/home/beangoben/anaconda3/lib/python3.7/site-packages/tqdm/autonotebook/__init__.py:14: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  " (e.g. in jupyter console)", TqdmExperimentalWarning)<jupyter_text>Data wranglingFirst we define our folder for the original data<jupyter_code>data_path = os.path.join(SRC_DIR, 'data/MP_47K/')
results_path = os.path.join(SRC_DIR,'results/eda')<jupyter_output><empty_output><jupyter_text>One JSON/molecule<jupyter_code>json_files = utils.get_all_jsons(data_path)
filename = os.path.join(data_path, random.choice(json_files))
with open(filename) as afile:
    sample = json.load(afile)
print(json.dumps(sample, indent=2))<jupyter_output>{
  "smiles": "SC1=NN=C(N1C1=CC=CC=C1)N1CCCCC1",
  "Tm": 290.15,
  "num_atms": 34,
  "gaussian": {
    "optdone": 1,
    "sigma_profile": "yes",
    "theory": "DFT",
    "functional": "B3LYP",
    "basiset": "6-31G**",
    "properties": {
      "R2extent": 7048.3766,
      "electric_dipole_moment_norm": 5.08583331323,
      "electric_quadrupole_moment_norm": 17.0597151377,
      "total_energy": -1122.19654864,
      "HOMO_LUMO_gap": 0.17583,
      "cosmo_polarized_energy": -14.55,
      "cosmo_cavitation_energy": 37.26,
      "cosmo_dispersion_energy": -26.38,
      "cosmo_repulsion_energy": 2.02,
      "cosmo_total_nonelectrostatic_energy": 12.9,
      "SMD_solvation_energy": 2.87
    },
    "coords": [
      {
        "x": -3.4092609553,
        "y": 3.4251511016,
        "z": -2.6974217351,
        "element": "S"
      },
      {
        "x": -2.5153130047,
        "y": 2.8910717272,
        "z": -1.2658043313,
        "element": "C"
      },
      {
        "x": -2.160112709,
     [...]<jupyter_text>SMILES to rdkit mol<jupyter_code>mol = Chem.MolFromSmiles(sample['smiles'])
MolToImage(mol)<jupyter_output><empty_output><jupyter_text>3D geometry<jupyter_code>sample['gaussian']['coords']<jupyter_output><empty_output><jupyter_text>Some electronic properties<jupyter_code>sample['gaussian']['properties']<jupyter_output><empty_output><jupyter_text>Valence orbitals<jupyter_code>sample['gaussian']['valence_orbital']<jupyter_output><empty_output><jupyter_text>Cosmo profile<jupyter_code>cosmo_profile = sample['gaussian']["cosmo_sigma_profile"]
plt.plot([float(i) for i in cosmo_profile.keys()], cosmo_profile.values())
plt.xlabel('Sigma')
plt.show()<jupyter_output><empty_output><jupyter_text>Wrangling a whole databaseUse the utility function `utils.read_json_data` along with `utils.get_all_jsons` with the **suffixes** argument (list of strings) to extract a subset of the data (i.g. *"OCHEM", "Bradley", "Enamine", "Bergstrom"*). If used with `None`, will read all data.<jupyter_code>def get_info(sample):
    """Function that process one json sample and returns a list of useful data."""
    index = sample['filename'].replace('.json', '').split('_')[-2]
    db = sample['filename'].replace('.json', '').split('_')[-1]
    return [sample['smiles'], sample['Tm'], db, index]

json_files = utils.get_all_jsons(data_path,suffixes=None)
data = utils.read_json_data(json_files,data_path, get_info)
df = pd.DataFrame(data, columns=['smiles', 'mp', 'db', 'db_index'])
print(df.shape)
df<jupyter_output><empty_output><jupyter_text>Add rdkit mol column<jupyter_code>df['mol'] =df['smiles'].apply(Chem.MolFromSmiles)<jupyter_output><empty_output><jupyter_text>Exploratory Data Analysis (EDA) Sets and colorsHere we define the database sets and associate colors to each label<jupyter_code>mp_dict = OrderedDict([("All", ["OCHEM", "Bradley", "Enamine", "Bergstrom"]),
                       ("Bradley + Bergstrom", ["Bradley", "Bergstrom"]),
                       ("OCHEM", ["OCHEM"]), ("Enamine", ["Enamine"])])

colors = sns.color_palette('Dark2', len(mp_dict))
db2color = {db: colors[indx] for indx, db in enumerate(mp_dict.keys())}

sns.palplot(colors)
plt.xticks(list(range(len(mp_dict))), list(mp_dict.keys()), rotation=90)
plt.show()<jupyter_output><empty_output><jupyter_text>Summary statistics<jupyter_code>stats_list = []
for label, db_list in mp_dict.items():
    sub_df = df[df['db'].isin(db_list)]
    stats = OrderedDict([('Label', label)])
    result_stats = sp.stats.describe(sub_df['mp'])._asdict()
    result_stats['std'] = np.sqrt(result_stats['variance'])
    stats.update(result_stats)
    stats_list.append(stats)
stats_df = pd.DataFrame(stats_list)
utils.save_latextable(stats_df,results_path,'stats')
stats_df<jupyter_output><empty_output><jupyter_text>Distribution plots<jupyter_code>alpha = 0.75
bins = 50

for index, (label, db_list) in enumerate(mp_dict.items()):
    y = df[df['db'].isin(db_list)]['mp'].values
    col = db2color[label]
    plt_label = utils.legend_stats_label(label, y,20)
    sns.distplot(y,
                 bins=bins,
                 kde=False,
                 hist_kws={
                     "histtype": 'step',
                     "color": col,
                     "alpha": alpha,
                     'linewidth': 3
                 },
                 kde_kws={
                     "color": col,
                     "alpha": alpha,
                     'shade': False,
                     'linewidth': 3
                 },
                 label=plt_label)

for xc in [50 + 273.15, 250 + 273.15]:
    plt.axvline(x=xc, c='k', linestyle='--', linewidth=3)

plt.xlabel('Melting Temperature (K)')
plt.ylabel('Samples')
utils.header_legend('',20, prop={'family': 'monospace', 'size': 'xx-small'}, title='Databases', loc='upper right',bbox_to_anchor=(1.3,1.0))
utils.save_figure(results_path,'hist')
plt.show()
plot_df = pd.DataFrame(
    dict([(label, pd.Series(df[df['db'].isin(db_list)]['mp'].values))
          for label, db_list in mp_dict.items()]))
sns.violinplot(data=plot_df,
               scale="area",
               inner='box',
               palette=colors,
               orient="v",
               linewidth=2.5)
plt.ylabel('Melting Temperature (K)')
utils.save_figure(results_path,'violin')
plt.show()<jupyter_output><empty_output><jupyter_text>Representative sample across MP<jupyter_code>bins = 20
encode = preprocessing.KBinsDiscretizer(n_bins=bins, encode='ordinal')
df['bin'] = encode.fit_transform(df['mp'].values.reshape(-1, 1))
mol_list, mp_list = [], []
smi_list = []
for index, grp in df.groupby(by='bin'):
    row = grp.sample(n=1).iloc[0]
    smi_list.append(row['smiles'])
    mol_list.append(row['mol'])
    mp_list.append('MP={:.2f}'.format(row['mp']))
svg= MolsToGridImage(mol_list, legends=mp_list, molsPerRow=5,useSVG=True)
utils.save_molgrid(svg,results_path,'sample')
utils.display_svg(svg)<jupyter_output><empty_output><jupyter_text>Finding similar molecules<jupyter_code>def fp_distance(fp1,fp2):
    return 1-DataStructs.DiceSimilarity(fp1,fp2)

radius=4
tol=0.4
get_fp = lambda mol: GetMorganFingerprint(mol,radius)

similar_df = df.copy()
similar_df['ref']= False
similar_df.loc[df.sample(n=1).index,'ref']=True
ref_row = similar_df.query('ref').iloc[0]
ref_mol,ref_mp = ref_row['mol'],ref_row['mp']
ref_fp = get_fp(ref_mol)
similar_df['dist'] = similar_df['mol'].apply(lambda mol: fp_distance(ref_fp, get_fp(mol)))
similar_df['similar'] = np.logical_and(similar_df['dist'] < tol, np.logical_not(similar_df['ref']))
print('Found {} similar molecules to the reference'.format(sum(similar_df['similar'])))
mol_list = [ref_mol]+similar_df.query('similar')['mol'].tolist()
legends=['ref: {:.1f}'.format(ref_mp)]+df['mp'].apply(lambda x: '{:.1f}'.format(x)).tolist()
svg= MolsToGridImage(mol_list, legends=legends, molsPerRow=4,useSVG=True)
utils.display_svg(svg)<jupyter_output><empty_output><jupyter_text>Structural Analysis via MP heatmap Reduce fingerprints via importance<jupyter_code>def counted_fp(mol):
    
from sklearn.ensemble import RandomForestRegressor

def rdkit_fp(mol,info={}):
    return Chem.RDKFingerprint(mol, minPath=3,maxPath=7,fpSize=2048, bitInfo=info)

def rf_bit_importances(mol_list,y,fp_func):
    model = RandomForestRegressor()
    x = np.vstack([fp_func(mol) for mol in mol_list])
    model.fit(x,y)
    return x, model.feature_importances_
y = df['mp'].values
x, bit_importances = rf_bit_importances(df['mol'].tolist(),y,rdkit_fp)
thresholds = np.logspace(0,-4,100)
threholded_bits = [sum(bit_importances > thres) for thres in thresholds]
elbow_index = utils.get_curve_elbow(threholded_bits)
mask = bit_importances > thresholds[elbow_index]
x = x[:,mask]
print('Using importance threshold of {}'.format(thresholds[elbow_index]))
print('FP is reduced to {} bits'.format(threholded_bits[elbow_index]))
n_bins=100
encode = preprocessing.KBinsDiscretizer(n_bins=n_bins, encode='ordinal')
bins = encode.fit_transform(y.reshape(-1, 1)).ravel()
x_avg = np.zeros((n_bins,sum(mask)))
for i in range(n_bins):
    x_avg[i,:]=np.mean(x[bins==i],axis=0)
sns.heatmap(x_avg.T,center=0)
def in_range(a,lower_lim,upper_lim):
    return lower_lim <= a < upper_lim

n_bins=100
percentile=False
n_bits = X_all.shape[1]

sub_df = df[df[calib_label].apply(in_range,args=(-350,-100))]
if percentile:
    bins = np.percentile(sub_df[calib_label],np.linspace(0.0,100.0,n_bins+1))
else:
    bins = np.linspace(-350,-100.0,n_bins+1)



plt.figure(figsize=(10,2))
sns.distplot(sub_df[calib_label].values,bins=bins)
plt.ylabel('Density')
plt.yticks([])
sns.despine()
plt.savefig('Redox_dist.svg',dpi=300)
plt.show()
sns.heatmap(X_avg.T,center=0)
plt.ylabel('RxnFP Bits')
plt.yticks([0,n_bits-1],[0,n_bits-1],rotation='horizontal')

if percentile:
    plt.xlabel('Percentile (%)')
    plt.xticks(np.arange(0,101,25),np.arange(0,101,25),rotation='vertical')
else:
    plt.xlabel('$E^{{m}^\prime}$ (pH=7) PM7 GP-calibrated ($mV$)')
    value_ticks = [ int(bins[math.ceil(i*n_bins)]) for i in np.linspace(0,1,5)]
    plt.xticks(np.arange(0,101,25),value_ticks,rotation='vertical')

plt.savefig('Rxn_heatmap.svg',dpi=300)
plt.savefig('Rxn_heatmap.png',dpi=300)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15874826827479424, 'fraction_numerical': 0.026648194931138458, 'mean_word_length': 5.003913894324853, 'pattern_counts': {'":': 37, '<': 81, '<?xml version=': 0, '>': 80, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '222130', 'license': 'MIT', 'n_tokens_mistral': 4921, 'n_tokens_neox': 4580, 'n_words': 756, 'path': 'notebooks/Exploratory_Data_Analysis.ipynb', 'repo_name': 'qchem64/melting_points_ml', 'size': '12271'}