<jupyter_start><jupyter_text>Using just the affinity values provided in the downloaded sdfs<jupyter_code>important_recs = experimental_df['Receptor'].unique().tolist()

BindingDBAffinities = pd.DataFrame(columns=['Receptor','BindingDB MonomerID','IC50'])
for file in glob('separated_sets/sdf_files/*.sdf'):
    if '3D.sdf' in file:
        continue
    liginfo = []
    receptor = re.findall(r'/(....)_Validation_Affinities.sdf',file)[0]
    if receptor not in important_recs:
        continue
    mono_id,ic="",0
    bdbid,affdata=False,False
    with open(file) as mol_file:
        for line in mol_file:
            if 'BindingDB monomerid' in line:
                bdbid=True
                continue
            if bdbid:
                mono_id = line.strip()
                bdbid=False
                continue
            if 'Enzymologic: Ki nM' in line:
                affdata=True
                continue
            if affdata:
                ic = float(line)
                affdata=False
                continue
            if 'Enzymologic' in line:
                print(line)
            if '$$$$' in line:
                liginfo.append([receptor,mono_id,ic])
                mono_id,ic="",0
    receptor_df = pd.DataFrame(data=liginfo,columns=['Receptor','BindingDB MonomerID','IC50'])
    BindingDBAffinities = BindingDBAffinities.append(receptor_df,ignore_index=True)
        
    
hist_of_diffs = []
for rec in important_recs:
    sdf_info = BindingDBAffinities[BindingDBAffinities['Receptor'] == rec]
    tsv_info = experimental_df[(experimental_df['Receptor'] == rec)]
    tsv_info.head()
    for lig in sdf_info['BindingDB MonomerID'].unique().tolist():
        sdf_val=sdf_info[sdf_info['BindingDB MonomerID']==lig]['IC50'].item()
        tsv_stuff = tsv_info[tsv_info['BindingDB MonomerID']==lig]['IC50 avg']
        if len(tsv_stuff) != 1:
            print(rec)
            print(tsv_stuff.head())
            continue
        tsv_val = tsv_stuff.item()
        diff = np.abs(sdf_val-tsv_val)
        hist_of_diffs.append(diff)

#plt.scatter(range(len(hist_of_diffs)),hist_of_diffs)
arraydif = np.array(hist_of_diffs)
print(len(arraydif))
smul = arraydif[arraydif>5]
zeros = arraydif[arraydif<1e-5]
plt.scatter(range(len(smul)),smul)
print(len(zeros))<jupyter_output><empty_output><jupyter_text>As you can see above, there is about half that are the same between the two data types (tsv and sdf). Though the half that are not identical between the two data sources have very major discrepancies.<jupyter_code>only_sdf_data = experimental_df.copy()
for idx, row in only_sdf_data.iterrows():
    only_sdf_data.at[idx,'IC50 sdf'] = BindingDBAffinities[(BindingDBAffinities['Receptor'] == row['Receptor']) & (BindingDBAffinities['BindingDB MonomerID'] == str(row['BindingDB MonomerID']))]['IC50'].item()
only_bef_lig_check = only_sdf_data[only_sdf_data['IC50 sdf'].notna()]
sdf_only_ic50_vals = only_bef_lig_check.groupby('Receptor').filter(lambda x: len(x) > 1) # make sure receptors have more than 1 ligand
sdf_only_ic50_vals['pIC50'] = sdf_only_ic50_vals['IC50 sdf'].apply(lambda x: -np.log10(x)) #note that this is the negative log10 of the value
print('{} unique receptors with greater than 1 ligand'.format(len(sdf_only_ic50_vals['Receptor'].unique())))
sdf_only_ic50_vals.sort_values(by=['Receptor','file'],inplace=True)
sdf_only_ic50_vals.to_csv('sdf_only_final_ic50_vals.txt', sep=' ', columns=['Receptor','file','BindingDB MonomerID','IC50 sdf','pIC50'],header=False,index=False)<jupyter_output><empty_output><jupyter_text>Make the Model Training Files (Types File)<jupyter_code>#stuff to make the training files
import requests
from bs4 import BeautifulSoup
def searchRCSBStructPage(s):
    ligs = []
    ligand_bigp = s.find(id="smallMoleculespanel")
    ligand_smallt = ligand_bigp.find(id="LigandsMainTable")
    for row in ligand_smallt.find_all(id=re.compile('ligand_row')):
        if 'UNL' in row['id']:
            continue
        data = row.find_all('td')[2]
        ligs.append(data.find_all('br')[1].next_sibling)
    return ligs
def getReferenceLigs(receptors):
    lig_dict = dict()
    for receptor in receptors:
        possible_ligs =[] 
        response = requests.get('https://www.rcsb.org/structure/{}'.format(receptor))
        soupyness = BeautifulSoup(response.content, 'html.parser')
        try:
            lig_dict[receptor] = searchRCSBStructPage(soupyness)
        except:
            try:
                newlig= soupyness.find(id='note_obsoletedBy')
                newreceptor = newlig.find(href=re.compile('/structure/')).text
                response = requests.get('https://www.rcsb.org/structure/{}'.format(newreceptor))
                soupyness = BeautifulSoup(response.content, 'html.parser')
                lig_dict[receptor] = searchRCSBStructPage(soupyness)
            except:
                print(receptor)
    return lig_dict
recs = sdf_only_ic50_vals['Receptor'].unique()
reference_dict = getReferenceLigs(recs)
train_groups = sdf_only_ic50_vals.groupby('Receptor')

sdf_only_ic50_vals['Reference'] = sdf_only_ic50_vals.apply(lambda x: x['Inchi Key'] in reference_dict[x.Receptor],axis=1)
sdf_only_ic50_vals.sort_values(by=['Receptor','Reference'],ascending=[True,False],ignore_index=True,inplace=True)
sdf_only_ic50_vals.to_csv('sdf_only_final_ic50_vals.txt', sep=' ', columns=['Receptor','file','BindingDB MonomerID','IC50 sdf','pIC50','Reference'],header=False,index=False)<jupyter_output><empty_output><jupyter_text>Start here if you want to use the IC50 values from the sdfs downloaded from the BindingDB page and the 475 receptors (maximal amount I have ever gotten)<jupyter_code>sdf_only_ic50_vals = pd.read_csv('sdf_only_final_ic50_vals.txt', sep=' ', header=None)
sdf_only_ic50_vals.columns = ['Receptor','file','BindingDB MonomerID','IC50 sdf','pIC50','Reference']<jupyter_output><empty_output><jupyter_text>maybe using combinations instead of permutations will be better (reduce it to minimum trainable set)<jupyter_code>possible_datas = [sdf_only_ic50_vals]
trainfile_names = ['sdf_only_old_bdb_comb.types']
for i,d in enumerate(possible_datas):
    trainfile = pd.DataFrame(columns=['class','reg','rec','lig1','lig2'])
    grouped = d.groupby('Receptor')
    for rec, group in grouped:
        groupinfo=[]
        for row1,row2 in list(combinations(group.index,2)): ##switched to combinations
            cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))
            dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])
            receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])
            lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])
            lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])
            if dif == np.inf or dif == -np.inf or math.isnan(dif):
                print(lig1,lig2,dif)
                continue
            groupinfo.append([cls,dif,receptor,lig1,lig2])
        group_df = pd.DataFrame(data=groupinfo,columns=['class','reg','rec','lig1','lig2'])
        trainfile = trainfile.append(group_df,ignore_index=True)
    trainfile.to_csv(trainfile_names[i],header=False,index=False,sep=' ',float_format='%.4f')<jupyter_output><empty_output><jupyter_text>Using the same train/test split as the DeltaDelta paperUsing the reference ligand and increasing amounts of additional ligands<jupyter_code>random=45
possible_datas = [sdf_only_ic50_vals]
trfile_names = ['sdfo_train_papersplit_rand{}_allpos'] ##note the allpos in the name because all affinity values are positive (may not make a difference)
ttfile_names = ['sdfo_test_papersplit_rand{}_allpos'] 
trainfile_names = [str_val.format(random) for str_val in trfile_names]
testfile_names = [str_val.format(random) for str_val in ttfile_names]
for i,d in enumerate(possible_datas):
    for addnl in range(1,7):
        trainfile = pd.DataFrame(columns=['class','reg','rec','lig1','lig2'])
        testfile = pd.DataFrame(columns=['class','reg','rec','lig1','lig2'])
        grouped = d.groupby('Receptor')
        for rec, group in grouped:
            traingroupinfo=[]
            testgroupinfo=[]
            if not group['Reference'].sum(): ## have to make sure there is one reference structure per group (arbitrarily assigning the reference to the first row if it isn't already assigned)
                group.at[group.index[0],'Reference'] = True
            if len(group)-1 < addnl:
                train_group = pd.concat([group[(group['Reference'] == True)],group[(group['Reference'] == False)].sample(n=len(group)-1,random_state=random)])
            else:
                train_group = pd.concat([group[(group['Reference'] == True)],group[group['Reference'] == False].sample(n=addnl,random_state=random)])
            ## Training set first
            for row1,row2 in list(combinations(train_group.index,2)): ##switched to combinations
                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])
                if dif == np.inf or dif == -np.inf or math.isnan(dif):
                    continue
#                 if dif < 0: ## Lets see if all positive DDG values helps training, so if regression value < 0 then swap order
#                     dif = np.abs(dif)
#                     tmp = row2
#                     row2 = row1
#                     row1 = tmp
                try:
                    cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))
                except:
                    print(train_group)
                    break
                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])
                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])
                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])
                traingroupinfo.append([cls,dif,receptor,lig1,lig2])
            # Test set. Using idea that can only compare to values that are in the training set
            for row1,row2 in list(combinations(group.index,2)):
                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])
                if dif == np.inf or dif == -np.inf or math.isnan(dif):
                    continue
#                 if dif < 0: ## Lets see if all positive DDG values helps training, so if regression value < 0 then swap order
#                     dif = np.abs(dif)
#                     tmp = row2
#                     row2 = row1
#                     row1 = tmp
                if not (bool(row1 in train_group.index) ^ bool(row2 in train_group.index)): ##supposed to be xor
                    continue
                cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))
                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])
                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])
                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])
                testgroupinfo.append([cls,dif,receptor,lig1,lig2])
            group_df = pd.DataFrame(data=traingroupinfo,columns=['class','reg','rec','lig1','lig2'])
            trainfile = trainfile.append(group_df,ignore_index=True)
            tgroup_df = pd.DataFrame(data=testgroupinfo,columns=['class','reg','rec','lig1','lig2'])
            testfile = testfile.append(tgroup_df,ignore_index=True)
        trainfile.to_csv('{}_{}.types'.format(trainfile_names[i],str(addnl)),header=False,index=False,sep=' ',float_format='%.4f')
        testfile.to_csv('{}_{}.types'.format(testfile_names[i],str(addnl)),header=False,index=False,sep=' ',float_format='%.4f')
# Making 5 different train/test splits based on the random seeds
# Additionally making two new columns before receptor for the absolute binding affinity of each ligand
# Use permutations of the ligands rather than combinations (order matters)
randos=[0,1,2,3,4]
trainfile_names = [f'sdfo_train_papersplit_rand{rand}_p' for rand in randos]
validfile_names = [f'sdfo_valid_papersplit_rand{rand}_p' for rand in randos]
testfile_names = [f'sdfo_testv_papersplit_rand{rand}_p' for rand in randos]
for i,random in enumerate(randos):
    for addnl in range(1,7):
        trainfile = pd.DataFrame(columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])
        testfile = pd.DataFrame(columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])
        validfile = pd.DataFrame(columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])
        grouped = sdf_only_ic50_vals.groupby('Receptor')
        for rec, group in grouped:
            traingroupinfo=[]
            testgroupinfo=[]
            validgroupinfo=[]
            if not group['Reference'].sum(): ## have to make sure there is one reference structure per group (arbitrarily assigning the reference to the first row if it isn't already assigned)
                group.at[group.index[0],'Reference'] = True
            if len(group)-1 < addnl:
                train_group = pd.concat([group[(group['Reference'] == True)],group[(group['Reference'] == False)].sample(n=len(group)-1,random_state=random)])
            else:
                train_group = pd.concat([group[(group['Reference'] == True)],group[group['Reference'] == False].sample(n=addnl,random_state=random)])
            not_train_group = set(group.index)-set(train_group.index)
            if len(not_train_group) > 2:
                valid,_ = train_test_split(list(set(group.index)-set(train_group.index)),test_size=0.6,random_state=random)
            else:
                valid = []
            ## Training set first
            for row1,row2 in list(permutations(train_group.index,2)): ##switched to combinations
                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])
                if dif == np.inf or dif == -np.inf or math.isnan(dif):
                    continue
                try:
                    cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))
                except:
                    print(train_group)
                    break
                dg_lig1 = float(group.loc[row1,'pIC50'])
                dg_lig2 = float(group.loc[row2,'pIC50'])
                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])
                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])
                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])
                traingroupinfo.append([cls,dif,dg_lig1,dg_lig2,receptor,lig1,lig2])
            # Test set. Using idea that can only compare to values that are in the training set (i.e. 'Known' values)
            for row1,row2 in list(permutations(group.index,2)):
                if not (bool(row1 in train_group.index) ^ bool(row2 in train_group.index)): ##supposed to be xor
                    continue
                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])
                if dif == np.inf or dif == -np.inf or math.isnan(dif):
                    continue
                dg_lig1 = float(group.loc[row1,'pIC50'])
                dg_lig2 = float(group.loc[row2,'pIC50'])
                cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))
                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])
                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])
                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])
                if row1 in valid or row2 in valid:
                    validgroupinfo.append([cls,dif,dg_lig1,dg_lig2,receptor,lig1,lig2])
                else:
                    testgroupinfo.append([cls,dif,dg_lig1,dg_lig2,receptor,lig1,lig2])
            group_df = pd.DataFrame(data=traingroupinfo,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])
            trainfile = trainfile.append(group_df,ignore_index=True)
            tgroup_df = pd.DataFrame(data=testgroupinfo,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])
            vgroup_df = pd.DataFrame(data=validgroupinfo,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])
            testfile = testfile.append(tgroup_df,ignore_index=True)
            validfile = validfile.append(vgroup_df,ignore_index=True)
        trainfile.to_csv(f'CV_Sets/{trainfile_names[i]}_{addnl}.types',header=False,index=False,sep=' ',float_format='%.4f')
        testfile.to_csv(f'CV_Sets/{testfile_names[i]}_{addnl}.types',header=False,index=False,sep=' ',float_format='%.4f')
        validfile.to_csv(f'CV_Sets/{validfile_names[i]}_{addnl}.types',header=False,index=False,sep=' ',float_format='%.4f')
group_rec = sdf_only_ic50_vals.groupby('Receptor')
full_types_list=[]
for rec_name, group in group_rec:
    for idx1, idx2 in list(combinations(group.index,2)):
        assert group.loc[idx1,'Receptor'] == group.loc[idx2,'Receptor']
        regression = float(group.loc[idx1,'pIC50']) - float(group.loc[idx2,'pIC50'])
        if regression == np.inf or regression == -np.inf or math.isnan(regression):
                    continue
#         if regression < 0: ## Lets see if all positive DDG values helps training, so if regression value < 0 then swap order
#             regression = np.abs(regression)
#             tmp = idx2
#             idx2 = idx1
#             idx1 = tmp
        try:
            classification = int(float(group.loc[idx1,'pIC50']) > float(group.loc[idx2,'pIC50']))
        except:
            print(train_group)
            break
        
        receptor = '{0}/{0}_R.pdb'.format(group.loc[idx1,'Receptor'])
        lig1 = '{0}/{1}'.format(group.loc[idx1,'Receptor'],group.loc[idx1,'file'])
        lig2 = '{0}/{1}'.format(group.loc[idx2,'Receptor'],group.loc[idx2,'file'])
        full_types_list.append([classification, regression, receptor,lig1,lig2])
full_types_df = pd.DataFrame(full_types_list,columns=['Classification','∆pIC50','Rec','Lig1','Lig2'])
full_types_df.to_csv('all_bdb_types.types',sep=' ',header=False,index=False)
# Uses permutations
# Outputs cls, reg, absaff1, absaff2, rec, lig1, lig2
group_rec = sdf_only_ic50_vals.groupby('Receptor')
full_types_list=[]
for rec_name, group in group_rec:
    for idx1, idx2 in list(permutations(group.index,2)):
        assert group.loc[idx1,'Receptor'] == group.loc[idx2,'Receptor']
        regression = float(group.loc[idx1,'pIC50']) - float(group.loc[idx2,'pIC50'])
        if regression == np.inf or regression == -np.inf or math.isnan(regression):
                    continue
        try:
            classification = int(float(group.loc[idx1,'pIC50']) > float(group.loc[idx2,'pIC50']))
        except:
            print(train_group)
            break
        dg_lig1 = float(group.loc[idx1,'pIC50'])
        dg_lig2 = float(group.loc[idx2,'pIC50'])
        receptor = '{0}/{0}_R.pdb'.format(group.loc[idx1,'Receptor'])
        lig1 = '{0}/{1}'.format(group.loc[idx1,'Receptor'],group.loc[idx1,'file'])
        lig2 = '{0}/{1}'.format(group.loc[idx2,'Receptor'],group.loc[idx2,'file'])
        full_types_list.append([classification, regression, dg_lig1, dg_lig2, receptor,lig1,lig2])
full_types_df_best = pd.DataFrame(full_types_list,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])
full_types_df_best.to_csv('all_bdb_types_perm_mult.types',sep=' ',header=False,index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15018944177822682, 'fraction_numerical': 0.023844405152816366, 'mean_word_length': 3.2922810060711187, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 35, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '533921', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6766, 'n_tokens_neox': 6458, 'n_words': 1086, 'path': 'OldBindingDB_Investigation.ipynb', 'repo_name': 'drewnutt/DDG', 'size': '19795'}