<jupyter_start><jupyter_text>Gather data from various source into the analysis repository<jupyter_code># imports
import os
import shutil
import numpy as np
import pandas as pd
import yaml
import git
import tempfile
import sys
sys.path.append(os.path.join(os.getcwd(), '..'))

import pint
unit_registry = pint.UnitRegistry()

from PLBenchmarks import targets, ligands, edges
from IPython.core.display import HTML

from tqdm.notebook import tqdm

import benchmarkpl
path = benchmarkpl.__path__[0]<jupyter_output>Warning: Unable to load toolkit 'OpenEye Toolkit'. The Open Force Field Toolkit does not require the OpenEye Toolkits, and can use RDKit/AmberTools instead. However, if you have a valid license for the OpenEye Toolkits, consider installing them for faster performance and additional file format support: https://docs.eyesopen.com/toolkits/python/quickstart-python/linuxosx.html OpenEye offers free Toolkit licenses for academics: https://www.eyesopen.com/academic-licensing<jupyter_text>Set path of data directory<jupyter_code>targets.set_data_dir(os.path.join(path))
# directory name where results for each target are stored
results_dir = '10_results'<jupyter_output><empty_output><jupyter_text>Number of targets, ligands and edges in the data set<jupyter_code>nligs, nedgs = 0, 0
print(f'{"Target":10s} {"Num Ligs":>10s} {"Num Edges":>10s}')
print(33 * '-')
for target in tqdm(targets.target_dict):
    print(f'{target:10s} {len(ligands.LigandSet(target)):10d} {len(edges.EdgeSet(target)):10d}')
    nligs += len(ligands.LigandSet(target))
    nedgs += len(edges.EdgeSet(target))
print(33 * '-')
print(f'{"total":10s} {nligs:10d} {nedgs:10d}')<jupyter_output>Target       Num Ligs  Num Edges
---------------------------------<jupyter_text>Experimental values stored in repositoryRetrieve the experimental values stored in the dataset<jupyter_code># function to retrieve exp. data from PLBenchmarks 
def getExpResults(target):
    edg = edges.EdgeSet(target)
    df = edg.get_dataframe(columns=[0,1, 'exp. DeltaG [kcal/mol]', 'exp. Error [kcal/mol]'])
    df.index = pd.Series([f'{target}_edge_' + str(lig1) + '_' + str(lig2) for lig1, lig2 in zip(df[0].values, df[1].values)])
   
    # remove unit of exp. values
    df['exp_DDG'] = df['exp. DeltaG [kcal/mol]'].apply(lambda x: x.magnitude)
    df['exp_dDDG'] = df['exp. Error [kcal/mol]'].apply(lambda x: x.magnitude)

    # filter only to relevant columns
    df = df.filter(items = [0,1,'exp_DDG', 'exp_dDDG'])
    df[0]= df[0].astype(str)
    df[1]= df[1].astype(str)
    return df
getExpResults('tyk2').head()
author = "hahn"
software = "experiment"
for target in targets.target_dict:
    df = getExpResults(target)
    if df is None:
        continue
    os.makedirs(os.path.join(path, targets.get_target_dir(target)), exist_ok=True)
    os.makedirs(os.path.join(path, targets.get_target_dir(target), results_dir), exist_ok=True)
    with open(os.path.join(path, targets.get_target_dir(target), results_dir, f'{target}_{software}_{author}.yaml'), 'w') as file:
        df = df.filter([0,1,'exp_DDG', 'exp_dDDG'])
        df.rename(columns={0: 'ligandA', 1: 'ligandB', 'exp_DDG': 'DDG', 'exp_dDDG': 'dDDG'}, inplace=True)
        df['unit']='kilocalories / mole'
        yaml.dump(df.T.to_dict(), file)<jupyter_output><empty_output><jupyter_text>pmx calculations with openFF parameters (Hahn et al.)<jupyter_code># function to retrieve data from PLBenchmarks calculations, soted in 00_data/input 
def getRawResults(target, forcefield='openff-1.0.0.offxml'):
    file_path = os.path.join(path, 
                                   '..', 
                                   '00_data', 
                                   'input',  
                                   f'{target}_{forcefield}.dat'
                                  )
    if not os.path.exists(file_path):
        print(f'File {file_path} does not exist.')
        return
    # read in result file
    res = pd.read_csv(file_path,
                      header=None, 
                      comment='#', 
                      skipinitialspace=True, 
                      names=['edge', 'calc DDG', 'calc dDDG' , 'add dDDG'], 
                      sep=' ')
    res.index = res['edge']
    # read in exp. data
    edg = edges.EdgeSet(target)
    df = edg.get_dataframe(columns=[0,1, 'exp. DeltaG [kcal/mol]', 'exp. Error [kcal/mol]'])
    df.index = pd.Series(['edge_' + str(lig1) + '_' + str(lig2) for lig1, lig2 in zip(df[0].values, df[1].values)])
   
    # copy in calculated values
    df['calc_DDG'] = res.loc[:,'calc DDG']
    df['calc_dDDG'] = res.loc[:,'calc dDDG']
    df['calc_dDDG(additional)'] = res.loc[:, 'add dDDG']
    
    # remove unit of calculated values
    df['exp_DDG'] = df['exp. DeltaG [kcal/mol]'].apply(lambda x: x.magnitude)
    df['exp_dDDG'] = df['exp. Error [kcal/mol]'].apply(lambda x: x.magnitude)

    # filter only to relevant columns
    df = df.filter(items = [0,1,'exp_DDG', 'exp_dDDG', 'calc_DDG', 'calc_dDDG', 'calc_dDDG(additional)'])
    df[0]= df[0].astype(str)
    df[1]= df[1].astype(str)
    df.index = pd.Series([f'{target}_edge_' + str(lig1) + '_' + str(lig2) for lig1, lig2 in zip(df[0].values, df[1].values)])
    return df
getRawResults('tyk2').head()
author = "hahn"
software = "pmx"
forcefield = "openff-1.0.0.offxml"
for target in targets.target_dict:
    df = getRawResults(target, forcefield)
    if df is None:
        continue
    os.makedirs(os.path.join(path, targets.get_target_dir(target), results_dir), exist_ok=True)
    with open(os.path.join(path, targets.get_target_dir(target), results_dir, f'{target}_{software}_{forcefield}_{author}.yaml'), 'w') as file:
        df = df.filter([0,1,'calc_DDG', 'calc_dDDG'])
        df.rename(columns={0: 'ligandA', 1: 'ligandB', 'calc_DDG': 'DDG', 'calc_dDDG': 'dDDG'}, inplace=True)
        df['unit']='kilocalories / mole'
        yaml.dump(df.T.to_dict(), file)<jupyter_output><empty_output><jupyter_text>Gather data from Gapsys et al.retrieve from https://github.com/deGrootLab/pmx<jupyter_code>temp_directory = tempfile.mkdtemp()

git.Repo.clone_from('https://github.com/deGrootLab/pmx', temp_directory, branch='master', depth=1)
def getGapsysResults(target):  
    file_name = f'{temp_directory}/protLig_benchmark/ddg_data/{target}.dat'
    if not os.path.exists(file_name):
        print(f'File {file_name} does not exist.')
        return None
    data = pd.read_csv(f'{temp_directory}/protLig_benchmark/ddg_data/{target}.dat', sep='\s+', header=None, comment='#',
                   names=['edge', 'exp', 'gaff', 'dgaff', 'cgenff', 'dcgenff', 'cons', 'dcons', 'fep5', 'dfep5', 'fep1', 'dfep1'])
    data['dexp']=pd.Series([0.0]*data.shape[0])
    df = edges.EdgeSet(target).get_dataframe()
    if target == 'jnk1':
        df.index = pd.Series([f'{str(a).split("-")[0]}_{str(b).split("-")[0]}' for a, b in zip(df[0].values, df[1].values)])
    else:
        df.index = pd.Series([f'{a}_{b}' for a, b in zip(df[0].values, df[1].values)])
    
    newdata = data.copy()
    newdata.index=newdata['edge']
    newdata.columns = pd.MultiIndex.from_arrays([np.array(newdata.columns), ['', 'exp'] + ['pmx'] * 6 + ['fep'] * 4 + ['exp'], [''] + ['kj/mol'] * 12], names=['forcefield', 'method', 'unit'])
    
    newdata.loc[:,('ligandA', '', '')] = df[0].apply(lambda x: str(int(float(x))) if str(x).endswith('.0') else x).astype(str)
    newdata.loc[:,('ligandB', '', '')] = df[1].apply(lambda x: str(int(float(x))) if str(x).endswith('.0') else x).astype(str)
    newdata = newdata.drop(columns=('edge', '', ''))
    newdata.sort_index(axis=1, level=1, inplace=True, sort_remaining=False)
    newdata.index=[f'{target}_edge_{x[("ligandA", "", "")]}_{x[("ligandB", "", "")]}' for i, x in newdata.iterrows()]
    
    return newdata
getGapsysResults('jnk1').head()
author = "gapsys"
for target in targets.target_dict:
    original_df = getGapsysResults(target)
    if original_df is None:
        continue
    os.makedirs(os.path.join(path, targets.get_target_dir(target), results_dir), exist_ok=True)
    for forcefield in ['fep5', 'fep1', 'gaff', 'cgenff']:
        df = original_df.loc[:, np.in1d(original_df.columns.get_level_values(0), ['ligandA', 'ligandB', f'{forcefield}', f'd{forcefield}'])].copy()
        software = list(df.columns.get_level_values(1))[list(df.columns.get_level_values(0)).index(f'{forcefield}')]
        unit = list(df.columns.get_level_values(2))[list(df.columns.get_level_values(0)).index(f'{forcefield}')]
        if unit == 'kj/mol':
            unit = 'kilojoules / mole'
        df.columns = df.columns.get_level_values(0)
        df.rename(columns={0: 'ligandA', 1: 'ligandB', f'{forcefield}': 'DDG', f'd{forcefield}': 'dDDG'}, inplace=True)
        df['unit'] = unit
        if forcefield.startswith('fep'):
            forcefield = f'opls3e_{forcefield[-1]}'
        with open(os.path.join(path, targets.get_target_dir(target), results_dir, f'{target}_{software}_{forcefield}_{author}.yaml'), 'w') as file:
            yaml.dump(df.T.to_dict(), file)
    
    forcefield = 'exp'
    df = original_df.loc[:, np.in1d(original_df.columns.get_level_values(0), ['ligandA', 'ligandB', f'{forcefield}', f'd{forcefield}'])].copy()
    software = list(df.columns.get_level_values(1))[list(df.columns.get_level_values(0)).index(f'{forcefield}')]
    unit = list(df.columns.get_level_values(2))[list(df.columns.get_level_values(0)).index(f'{forcefield}')]
    if unit == 'kj/mol':
        unit = 'kilojoules / mole'
    df.columns = df.columns.get_level_values(0)
    df.rename(columns={0: 'ligandA', 1: 'ligandB', f'{forcefield}': 'DDG', f'd{forcefield}': 'dDDG'}, inplace=True)
    df['unit'] = unit
    with open(os.path.join(path, targets.get_target_dir(target), results_dir, f'{target}_experiment_{author}.yaml'), 'w') as file:
        yaml.dump(df.T.to_dict(), file)<jupyter_output><empty_output><jupyter_text>Compare experimental values stored in repository with Gapsys exp. data<jupyter_code>for target in targets.target_dict:
    author = "hahn"
    software = "experiment"
    file_name = os.path.join(path, targets.get_target_dir(target), results_dir, 
                                   f'{target}_{software}_{author}.yaml'
                       )
    if os.path.exists(file_name):
        with open(file_name, 'r') as file:
            data1 = yaml.safe_load(file)
    else:
        print(f"File {file_name} for target {target} not available")
        continue
        
    author = "gapsys"
    software = "experiment"
    file_name = os.path.join(path, targets.get_target_dir(target), results_dir, 
                                   f'{target}_{software}_{author}.yaml'
                       )
    if os.path.exists(file_name):
        with open(file_name, 'r') as file:
            data2 = yaml.safe_load(file)
    else:
        print(f"File {file_name} for target {target} not available")
        continue
        
    for e, edata in data1.items():
        if e in data2:
            v1 = unit_registry.Quantity(edata['DDG'], edata['unit'])
            v2 = unit_registry.Quantity(data2[e]['DDG'], data2[e]['unit'])
            if not np.isclose(v1.to('kilocalories / mole').magnitude, 
                              v2.to('kilocalories / mole').magnitude, 
                              atol=.05,
                              equal_nan=False):
                print(target, e, v1.to('kilocalories / mole'), v2.to('kilocalories / mole'))
        else:
            print(target, e)<jupyter_output>bace bace_edge_CAT-24_CAT-17e 1.74 kilocalorie / mole 1.32887189292543 kilocalorie / mole
bace bace_edge_CAT-24_CAT-17i 2.29 kilocalorie / mole 1.8809751434034416 kilocalorie / mole<jupyter_text>Get results from Wang et. al., JACS 2015Input file taken from https://pubs.acs.org/doi/suppl/10.1021/ja512751q/suppl_file/ja512751q_si_003.xlsx (retrieval 2020-09-21) and converted to csv file<jupyter_code>def getWangResults(target):    
    file_name = f'{path}/../00_data/input/ja512751q_si_003.csv'
    if not os.path.exists(file_name):
        print(f'File {file_name} does not exist.')        
        return None
    data = pd.read_csv(file_name, sep=',')
    data['system'] = data['system'].apply(lambda x: str(x).lower())
    indeces = data.loc[data['system']!='nan', :].index
    data_per_target = [data.loc[indeces[i]:indeces[i+1]-1,:] for i in range(indeces.shape[0]-1)]
    data_per_target = {d['system'].iloc[0]: d for d in data_per_target}
    for t, d in data_per_target.items():
        d.index = [f'{target}_edge_{x["Ligand1"]}_{x["Ligand2"]}' for i, x in d.iterrows()]
        d = d[['Ligand1', 'Ligand2', 'bennett_ddG', 'bennett_error']]
        d = d.rename(columns={'Ligand1': 'LigandA',
                             'Ligand2': 'LigandB',
                             'bennett_ddG': 'DDG',
                             'bennett_error': 'dDDG'})
        d['unit'] = 'kilocalories / mole'
        data_per_target[t] = d
    if target in data_per_target:
        return data_per_target[target]
    else:
        return None
getWangResults('jnk1').head()
author = 'wang'
software = 'fep+'
forcefield = 'opls2.1'
for target in targets.target_dict:
    df = getWangResults(target)
    if df is None:
        continue
    os.makedirs(os.path.join(path, targets.get_target_dir(target), results_dir), exist_ok=True)
    with open(os.path.join(path, targets.get_target_dir(target), results_dir, f'{target}_{software}_{forcefield}_{author}.yaml'), 'w') as file:
        yaml.dump(df.T.to_dict(), file)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16500184979652238, 'fraction_numerical': 0.018275989641139474, 'mean_word_length': 3.3983078424991864, 'pattern_counts': {'":': 4, '<': 30, '<?xml version=': 0, '>': 32, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '137098', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4979, 'n_tokens_neox': 4795, 'n_words': 924, 'path': '00_data/00b_gather_data.ipynb', 'repo_name': 'dfhahn/protein-ligand-benchmark-analysis', 'size': '13515'}