<jupyter_start><jupyter_text>This notebook is an analysis of the errors in relative free energy calculations from the Schrödinger JACS dataset:> <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., ... & <NAME>. (2015). Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field. Journal of the American Chemical Society, 137(7), 2695-2703.http://doi.org/10.1021/ja512751q<jupyter_code>%matplotlib inline
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pylab as plt
# Read the Excel sheet
df = pd.read_excel('ja512751q_si_003.xlsx', sheet_name='ddG')
# Delete rows with summary statistics
rows_to_drop = list()
for i, row in df.iterrows():
    if str(df.loc[i,'Ligand1']) == 'nan':
        rows_to_drop.append(i)
print('dropping rows: {}'.format(rows_to_drop))
df = df.drop(index=rows_to_drop);
# Populate 'system' field for each entry
system = df.loc[0,'system']
for i, row in df.iterrows():
    if str(df.loc[i,'system']) == 'nan':
        df.loc[i, "system"] = system
    else:
        system = df.loc[i, "system"]
df
def bootstrap_statistic(errors, ci=0.95, statistic='rmse', nbootstrap = 1000):
    """Compute mean and confidence intervals of specified statistic.
    
    Parameters
    ----------
    errors : ndarray of shape (N,)
        The errors to bootstrap
    ci : float, optional, default=0.95
        Interval for CI
    statistic : str, optional, default='rmse'
        Statistic to compute; one of ['rmse', 'mue']
    nbootstrap : int, optional. default=10000
        Number of bootstrap samples
        
    Returns
    -------
    mean : float
        The mean statistic for the whole dataset
    stderr : float
        The standard error
    low, high : float
        Low and high ends of CI    
    """
    
    def compute_statistic(errors, statistic):
        if statistic == 'rmse':
            return np.sqrt(np.mean(errors**2))
        elif statistic == 'mue':
            return np.mean(np.abs(errors))
        else:
            raise Exception("unknown statistic '{}'".format(statistic))
            
    sample_size = len(errors)
    rmse_n = np.zeros([nbootstrap], np.float64)
    for replicate in range(nbootstrap):
        bootstrapped_sample = np.random.choice(np.arange(sample_size), size=[sample_size])
        rmse_n[replicate] = compute_statistic(errors[bootstrapped_sample], statistic)
        
    rmse_n = np.sort(rmse_n)
        
    rmse = compute_statistic(errors, statistic)
    drmse = np.std(rmse_n)
    
    low_frac = (1.0-ci)/2.0
    high_frac = 1.0 - low_frac
    rmse_low = rmse_n[int(np.floor(nbootstrap*low_frac))]
    rmse_high = rmse_n[int(np.ceil(nbootstrap*high_frac))]

    return rmse, drmse, rmse_low, rmse_high
# Collect data by system
def plot_data(rows):
    exp_ddg = rows['exp_ddG'].values
    ccc_ddg = rows['ccc_ddG'].values
    ccc_error = rows['ccc_error'].values
    error = rows['Error_ccc'].values
    xmin = min(exp_ddg.min(), ccc_ddg.min()) - 0.5
    xmax = min(exp_ddg.max(), ccc_ddg.max()) + 0.5
    plt.xlabel('experiment $\Delta \Delta G$ (kcal/mol)');
    plt.ylabel('cycle closure $\Delta \Delta G$ (kcal/mol)');
    plt.plot([xmin, xmax], [xmin, xmax], 'b-');
    plt.plot(exp_ddg, ccc_ddg, 'ko');
    plt.plot([exp_ddg, exp_ddg], [ccc_ddg-2*ccc_error, ccc_ddg+2*ccc_error], 'k-', alpha=0.3);

    rmse, drmse, rmse_low, rmse_high = bootstrap_statistic(error, statistic='rmse')
    mue, dmue, mue_low, mue_high = bootstrap_statistic(error, statistic='mue')
    title = """{} (N = {})
    RMSE {:.1f} [95%: {:.1f}, {:.1f}] kcal/mol
    MUE {:.1f} [95%: {:.1f}, {:.1f}] kcal/mol
    """.format(system, len(exp_ddg), rmse, rmse_low, rmse_high, mue, mue_low, mue_high)
    plt.title(title);
    
    plt.axis('square');
    plt.xlim(xmin, xmax);
    plt.ylim(xmin, xmax);

systems = df['system'].unique()
nsystems = len(systems)
nx = int(np.ceil(np.sqrt(nsystems)))
ny = int(np.ceil(np.sqrt(nsystems)))
fig = plt.figure(figsize=[12,12])
for plot_index, system in enumerate(systems):
    plt.subplot(ny, nx, plot_index+1)
    rows = df.query("system == '{}'".format(system))
    plot_data(rows)
    
plt.subplot(ny, nx, nsystems+1)
system = 'all'
plot_data(df)
    
fig.tight_layout()
fig.savefig('jacs-analysis.pdf');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14982817869415807, 'fraction_numerical': 0.01993127147766323, 'mean_word_length': 3.3572854291417165, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '462610', 'license': 'MIT', 'n_tokens_mistral': 1648, 'n_tokens_neox': 1544, 'n_words': 372, 'path': 'notebook.ipynb', 'repo_name': 'jchodera/jacs-dataset-analysis', 'size': '4378'}