<jupyter_start><jupyter_text>1. Name and properties of compounds in small clusters of manifold plots2. Use clustering algorithm to find clusters3. Plot molecules based on how large they are. Use features such as molecular weight, volume, number of atoms. What defines the size of molecules?4. EDA - find class specific features5. General toxicity vs. specific toxicity<jupyter_code>from google.colab import drive
drive.mount('/gdrive')
%cd /gdrive
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
os.listdir('My Drive/Colab Notebooks/Others/')
from sklearn import manifold
from sklearn.utils import shuffle
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from collections import OrderedDict
from functools import partial
import rdkit.Chem as Chem
from rdkit.Chem import AllChem

def smiles_to_fps(data, fp_length, fp_radius):
    return stringlist2intarray(np.array([smile_to_fp(s, fp_length, fp_radius) for s in data]))

def smile_to_fp(s, fp_length, fp_radius):
    m = Chem.MolFromSmiles(s)
    return (AllChem.GetMorganFingerprintAsBitVect(
        m, fp_radius, nBits=fp_length)).ToBitString()

def stringlist2intarray(A):
    '''This function will convert from a list of strings "10010101" into in integer numpy array.'''
    return np.array([list(s) for s in A], dtype=int)
# load the feature dataset
features_df = pd.read_csv("My Drive/Colab Notebooks/Others/total_df.csv", index_col=[0])
features_df.head()
index = features_df.index
features_df = shuffle(features_df, random_state=43)
features_df.index = index
features_df.head()
feature_compound_name = features_df['compound'].values
class_label_name = ['Ototoxic' if y == 1.0 else 'Non-toxic' for y in features_df['class_label'].values]
Y_feature = features_df['class_label'].values
X_feature = features_df.drop(columns=['class_label', 'compound']).values
X_feature = StandardScaler().fit_transform(X_feature)
Y_feature.shape, X_feature.shape
import plotly.express as px
def run_manifold(X, Y, name, n_components=2):
    # setup manifold methods
    methods = OrderedDict()
#     methods['MDS'] = manifold.MDS(n_components, max_iter=100, n_init=1)
    methods['t-SNE'] = manifold.TSNE(n_components=n_components, init='pca', random_state=0)
#     methods['PCA'] = PCA(n_components=n_components)
    
    for i, (label, method) in enumerate(methods.items()):
        Y_ = method.fit_transform(X)
        if n_components == 2:
            fig = px.scatter(Y_, x=0, y=1, color=class_label_name, hover_name=feature_compound_name, 
                            title=f'Manifold method {label}', labels={'0': 'PC 1', '1': 'PC 2'})
        elif n_components == 3:
            fig = px.scatter_3d(Y_, x=0, y=1, z=2, color=class_label_name, hover_name=feature_compound_name, 
                        title=f'Manifold method {label}', labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'})
        fig.show()
run_manifold(X_feature, Y_feature, "PyBioMed Features", n_components=2)<jupyter_output>/usr/local/lib/python3.7/dist-packages/sklearn/manifold/_t_sne.py:793: FutureWarning: The default learning rate in TSNE will change from 200.0 to 'auto' in 1.2.
  FutureWarning,
/usr/local/lib/python3.7/dist-packages/sklearn/manifold/_t_sne.py:986: FutureWarning: The PCA initialization in TSNE will change to have the standard deviation of PC1 equal to 1e-4 in 1.2. This will ensure better convergence.
  FutureWarning,<jupyter_text>Clustering using the t-SNE components. Note that t-SNE components changes upon successive runs.<jupyter_code>sne_components = manifold.TSNE(n_components=2, init='pca', random_state=0).fit_transform(X_feature)
from sklearn.cluster import KMeans
from sklearn.metrics.pairwise import pairwise_distances_argmin
# run clustering
n_clusters=20
kmeans = KMeans(init='k-means++', n_clusters=n_clusters, n_init=10, random_state=0)
kmeans.fit(sne_components)
# we want to have same colors for the same clusters. 
# let's pair the cluster centers per closest one. 
kmeans_cluster_center = kmeans.cluster_centers_
kmeans_labels = pairwise_distances_argmin(sne_components, kmeans_cluster_center)
kmeans_labels_str = ["Cluster "+str(p) for p in kmeans_labels]
fig = px.scatter(sne_components, x=0, y=1, 
                 color=kmeans_labels_str, 
                 hover_name=[pp + " " + qq for pp, qq in zip(feature_compound_name, class_label_name)],
                            title=f'k-Means clustering of t-SNE components', labels={'0': 'PC 1', '1': 'PC 2'})
fig.show()
# plot the clusters
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 5))
colors = ["#4EACC5", "#FF9C34", "#d2c596", "#a254a2", "#289ec1", "#18e4ac", "#072f8b", "#b46c9a", "#9f0c05", 
         "#5974d1", "#0dccb9", "#f264d1", "#b0b071", "#da4634", "#704e31", "#ee8b6a", "#eae883", "#116549", 
         "#f8a82a", "#6c776f"]


for k, col in zip(range(n_clusters), colors):
    my_members = kmeans_labels == k
    cluster_center = kmeans_cluster_center[k]
    ax.plot(sne_components[my_members, 0], sne_components[my_members, 1], "w", markerfacecolor=col, marker=".", markersize=10)
    ax.plot(
        cluster_center[0],
        cluster_center[1],
        "o",
        markerfacecolor=col,
        markeredgecolor="k",
        markersize=6,
    )
ax.set_title("KMeans Clustering of t-SNE components")
ax.set_xticks(())
ax.set_yticks(())<jupyter_output><empty_output><jupyter_text>Now, take the clusters and study what features the molecules in the clusters have in common.<jupyter_code><jupyter_output><empty_output><jupyter_text>What defines the size of the molecules?Ototoxicity is caused by blocking of the MET-channel. Is there any correlation between the size of the molecules and the class labels. We will use the following features as proxy for molecular volume (or size) to do some analysis. 1. Molecular weight (Weight)2. Number of all atoms (nta)3. Number of O atoms (noxy)4. Number of C atoms (ncarb)5. Number of single bonds (nsb)We also have counts of various atoms such as1. Hydrogen2. Halogens3. Heavy atoms4. F atoms5. Cl atoms6. Br atoms7. I atoms8. C atoms9. P atoms10. S atoms11. O atoms12. N atomsOther measures for molecular size are bond length, position of atoms, volume which can be computed from Molar mass and density.<jupyter_code>molecular_size_features = ['Weight', 'nta', 'noxy', 'ncarb', 'nsb']
for indicator in molecular_size_features:
    fig = px.scatter(
        features_df,
        y=indicator,
        x=np.arange(0, features_df.shape[0]),
        color = class_label_name,
        hover_name=features_df['compound'],
        title=f'Scatter plot of {indicator} feature as proxy for molecular size',
        width=700, 
        height=500
    )
    fig.show()<jupyter_output><empty_output>