<jupyter_start><jupyter_text>This isn't normally distrbibuted - therefore, could 'class_weight='balanced' - or, could use tree-based models that perform well on imbalanced datsets as their hierarchcial strucutre alllows them to learn signals from both classses  Lipinski moleculesThis section identifies all molecules that satisfy all criteria of Lipinski's rule of five - for bioavaliable small molecules* 5 or fewer hydrogen bond donors;* 10 or fewer hydrogen bond acceptors;* A molecular weight (MW) of less than 500 Daltons;* An octanol-water partition coefficient (log $P_{o/w}$) of less than 5. Larger log $P_{o/w}$ meansmore lipophilic (i.e., less water soluble).<jupyter_code># Lipinski molecules

def lipinski_mol(df):
    """Determines whether a molecule in dataframe follows Lipinski rules"""
    few_hb_donor = (df["H_donors"] <= 5)
    few_hb_acceptor = (df["H_acceptors"] <= 10)
    low_MW = (df["MW"] < 500)
    low_pc = (df["LogP"] < 5) 
    return (few_hb_acceptor & few_hb_donor & low_MW & low_pc)

def select_lipinski(df, rule = lipinski_mol):
    """Removes compounds from database which don't follow Lipinski rules"""
    df_lip = analyser.analyse_compounds(df, df, 'lipinski_mol', rule)
    # df is passed twice, as it contains the data to operate on, and also is the dataframe to add output to
    df_lip.drop(df_lip[df_lip.lipinski_mol == False].index, inplace=True)
    df_lip.drop(columns = ['lipinski_mol'], inplace=True)  # As all values are now true so column is defunct
    return df_lip


df_lip_compounds = select_lipinski(df_compounds.copy())

print(f"Selected {df_lip_compounds.shape[0]} molecules out of {df_compounds.shape[0]}.")

def select_lipinski_assay(df_c, df_a, rule = lipinski_mol):
    """Removes assay corresponding to compounds from database which don't follow Lipinski rules"""
    df_lip = analyser.analyse_compounds(df_c, df_c, 'lipinski_mol', rule)
    df_a['lipinski_mol'] = df_lip['lipinski_mol']
    df_a.drop(df_a[df_a.lipinski_mol == False].index, inplace=True)
    df_a.drop(columns = ['lipinski_mol'], inplace=True)  # As all values are now true so column is defunct
    return df_a

df_lip_assays = select_lipinski_assay(df_compounds.copy(), df_assay.copy())

print(f"Remaining: {df_lip_compounds.shape[0]} assays out of {df_compounds.shape[0]}.")
# Representation of Lipinski molecules in sample space

def compare_dataframe_pair(df_1, df_2, x_var, y_var, subset_name = 'Lipinski Molecules'):
    try:
        fig, ax = plt.subplots()
        df_1.plot(x=x_var, y=y_var, ax=ax, style='o', color = 'grey', label = 'Molecule Space')
        df_2.plot(x=x_var, y=y_var, ax=ax, style='o', color = 'g', label = subset_name) 
        ax.set_ylabel(y_var); ax.legend()
    except KeyError as k:
        print(f"Column name {k} does not exist in dataframe")
            

compare_dataframe_pair(df_compounds, df_lip_compounds, 'H_acceptors', 'H_donors')
compare_dataframe_pair(df_compounds, df_lip_compounds, 'MW', 'LogP')<jupyter_output><empty_output><jupyter_text>Note that the MW and $P_{o/w}$ are the most significant factors in determining whether a molecule satifes the Lipinski conditions, with the hydrogen bond acceptor and donor limits ruling out very few molecules that have not been ruled out by these factors. Relaxed Lipinski conditionsWe may also consider the subset of molecules that satisfy at least three of the Lipinski conditions, as shown below:<jupyter_code>def relaxed_lipinski(df):
    """Determines whether a molecule in dataframe follows all but one Lipinski rules"""
    few_hb_donor = (df["H_donors"] <= 5) * 1  # (* 1) to convert to int
    few_hb_acceptor = (df["H_acceptors"] <= 10) * 1
    low_MW = (df["MW"] < 500) * 1
    low_pc = (df["LogP"] < 5) * 1
    accepted_cond = few_hb_acceptor + few_hb_donor + low_MW + low_pc
    return (accepted_cond >= 3)

df_rlip_compounds = select_lipinski(df_compounds.copy(), rule=relaxed_lipinski)

print(f"Selected {df_rlip_compounds.shape[0]} molecules out of {df_compounds.shape[0]}.")<jupyter_output>Selected 1950 molecules out of 2037.<jupyter_text>We see that this new criterion includes a further 422 molecules, compared to the stricter condition where all Lipinski conditions must be passed. We expect from our previous work that most of these new molecules will have broken the MW or $P_{o/w}$ conditions, but may verify this graphically.<jupyter_code>def compare_dataframe_trio(df_1, df_2, df_3, x_var, y_var):
    try:
        fig, ax = plt.subplots()
        df_1.plot(x=x_var, y=y_var, ax=ax, style='o', color = 'grey', label = 'Molecule Space')
        df_2.plot(x=x_var, y=y_var, ax=ax, style='o', color = 'blue', label = 'Relaxed Lipinski Molecules') 
        df_3.plot(x=x_var, y=y_var, ax=ax, style='o', color = 'g', label = 'Lipinski Molecules') 
        ax.set_ylabel(y_var); ax.legend()
    except KeyError as k:
        print(f"Column name {k} does not exist in dataframe")
            

compare_dataframe_trio(df_compounds, df_rlip_compounds, df_lip_compounds, 'H_acceptors', 'H_donors')
compare_dataframe_trio(df_compounds, df_rlip_compounds, df_lip_compounds, 'MW', 'LogP')<jupyter_output><empty_output><jupyter_text>As expected, we see two new tranches of molecules (in blue) on the MW/$P_{o/w}$ plot, that have been accepted under the new relaxed conditions. In general, outliers that lie further from the main cluster have not been accepted - these heavier molecules are more likely to fail the hydrogben bond donor/acceptor conditions, as they contain more functional groups. We may also consider the distribution of these compounds across subspaces of other assay data:<jupyter_code>compare_dataframe_trio(df_compounds, df_rlip_compounds, df_lip_compounds, 'NHOH_count', 'NO_count')
compare_dataframe_trio(df_compounds, df_rlip_compounds, df_lip_compounds, 'Heavy_atoms', 'Heavy_atom_mw')
compare_dataframe_trio(df_compounds, df_rlip_compounds, df_lip_compounds, 'Rotatable_bonds', 'Valence_electonrs')
compare_dataframe_trio(df_compounds, df_rlip_compounds, df_lip_compounds, 'Rings', 'TPSA')

#tpsa is good, rings are not - compare tpsa to valence e?
#nhoh is h bond donor, no is hydrogen bond acceptor
print(df_lip_compounds.columns)<jupyter_output>Index(['SMILES', 'H_acceptors', 'H_donors', 'MW', 'LogP', 'Heavy_atoms',
       'Heavy_atom_mw', 'NHOH_count', 'NO_count', 'Rotatable_bonds',
       'Valence_electonrs', 'Rings', 'TPSA', 'fingerprint'],
      dtype='object')<jupyter_text>Clustering We have computed ECFP (Morgan) fingerprints with radius 2 and 2048 bits for each of the Moonshot compounds. We will then cluster the compounds using the Tanimoto (a.k.a. Jaccard) index. The remove zero size clusters?or remove zero variance bits and conduct pca analysis<jupyter_code>#Define clustering setup
def ClusterFps(fps,cutoff=0.2):
    from rdkit import DataStructs
    from rdkit.ML.Cluster import Butina

    # first generate the distance matrix:
    dists = []
    nfps = len(fps)
    for i in range(1,nfps):
        sims = DataStructs.BulkTanimotoSimilarity(fps[i],fps[:i])
        dists.extend([1-x for x in sims])

    # now cluster the data:
    cs = Butina.ClusterData(dists,nfps,cutoff,isDistData=True)
    return cs
clusters=ClusterFps(df_compounds['fingerprint'],cutoff=0.322)  # for 95% confidence that clusters are different

def cluster_size_dist(clusters):
    """ Records the size of each cluster.
    """
    cluster_sizes = []
    for cluster in clusters:
        cluster_sizes.append(len(cluster))
    return cluster_sizes

cluster_size_dist = cluster_size_dist(clusters)
print(f"Generated {len(clusters)} with mean occupancy of {sum(cluster_size_dist)/len(clusters):.2f} molecules")<jupyter_output>Generated 911 with mean occupancy of 2.24 molecules<jupyter_text>Our clustering algorithm has generated a large number of clusters, with very few molecules per cluster. What does this mean though - are our molecules evenly spread between clusters?<jupyter_code>plt.figure(figsize=(14,6))
plt.subplot(121)
hist, bins, _ = plt.hist(cluster_size_dist, bins=12)
plt.xlabel('Cluster Size'); plt.ylabel('Frequency')
plt.title("Linear Scale")

plt.subplot(122)
logbins = np.logspace(np.log10(bins[0]),np.log10(bins[-1]),len(bins))
plt.hist(cluster_size_dist, bins=logbins)
plt.xscale('log'); plt.yscale('log')
plt.xlabel('Cluster Size'); plt.ylabel('Frequency')
plt.title("Log Scale")

plt.show()<jupyter_output><empty_output><jupyter_text>From the linear histogram, it is clear that the vast majority of the clusters are small; in fact the log scale shows the frequency of clusters decreases almost exponentially with size. We believe this over-representation of small clusters is a result of an extremely diverse dataset, with many unique molecules unrelated to others in the dataset (and a small number of well clustered molecules, ie derived from the same scaffold).<jupyter_code># find features associated with the most common clusters, and display graphically.
# remove molecules in smallest clusters and recluster?
# or remove zero variance bits in fingerprint and repeat.
def find_max_cluster(clusters):
    """ Returns the largest cluster
    """
    max_cluster = ()
    for cluster in clusters:
        if len(cluster) > len(max_cluster):
            max_cluster = cluster
    return max_cluster

max_cluster = find_max_cluster(clusters)
print(f"The largest cluster has {len(max_cluster)} compounds.")<jupyter_output>The largest cluster has 143 compounds.<jupyter_text>This cluster can be visualised within the sample space, for example in terms of MD and $P_{o/w}$:<jupyter_code>compare_dataframe_pair(df_compounds, df_compounds.iloc[list(max_cluster),], 'MW', 'LogP', subset_name='Largest Cluster')
compare_dataframe_pair(df_compounds, df_compounds.iloc[list(max_cluster),], 'NHOH_count', 'NO_count', subset_name='Largest Cluster')
compare_dataframe_pair(df_compounds, df_compounds.iloc[list(max_cluster),], 'TPSA', 'Valence_electonrs', subset_name='Largest Cluster')

'NHOH_count', 'NO_count'<jupyter_output><empty_output><jupyter_text>We may also visualise molecules from the cluster, using RDKit's built-in drawer:<jupyter_code>n_molecule = 3
cluster_mols=[]
for i in range(n_molecule):
    m = Chem.MolFromSmiles(df_compounds['SMILES'][max_cluster[i]])
    cluster_mols.append(m)

Draw.MolsToGridImage(cluster_mols)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10286433133346236, 'fraction_numerical': 0.009386491194116508, 'mean_word_length': 5.395420792079208, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 35, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '259875', 'license': 'MIT', 'n_tokens_mistral': 3622, 'n_tokens_neox': 3342, 'n_words': 1002, 'path': 'populate_df.ipynb', 'repo_name': 'KCGallagher/moonshot-predictor', 'size': '10334'}