import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.cluster import KMeans
from sklearn import metrics
df = pd.read_csv('protein_phi_psi.csv')
df.head()
df.info()
df.isna().sum()
The code shown below reads the phi and psi values from the given dataset and creates a scatter plot based on the two columns. Below that, there is a generated heatmap based on the same values.
plt.figure(figsize=(6,6))
plt.set_cmap("viridis")
plt.scatter(df['phi'], df['psi'], s=0.3)
plt.xticks(np.arange(-180, 181, 30))
plt.yticks(np.arange(-180, 181, 30))
plt.margins(0.001)
plt.xlabel('phi (°)')
plt.ylabel('psi (°)')
plt.title("Phi and Psi Combinations")
plt.plot([-180, 180], [0,0], linestyle = '--', color = 'grey')
plt.plot([0,0], [-180, 180], linestyle = '--', color = 'grey')
plt.show()
plt.figure(figsize=(7,6))
plt.hist2d(df.phi, df.psi, bins=40)
plt.xlabel('phi (°)')
plt.ylabel('psi (°)')
plt.title("Phi and Psi Combinations")
plt.xticks(np.arange(-180, 181, 30))
plt.yticks(np.arange(-180, 181, 30))
cb = plt.colorbar()
cb.set_label('Number of samples')
plt.show()
In order to experiment with different values for K, I created models with k = 2, 3, 4, 5, 6, and 7, and plotted an elbow graph to see which value of k seems to be appropreate.
x = np.array(df['phi'])[..., np.newaxis]
y = np.array(df['psi'])[..., np.newaxis]
X = np.hstack((x, y))
k_list = range(2,7)
distance_list = []
for k in k_list:
kmeans = KMeans(n_clusters=k, random_state=0, n_init='auto').fit(df[["phi","psi"]])
distance_list.append(kmeans.inertia_)
plt.figure(figsize=(8,6))
plt.plot(k_list, distance_list, marker=".")
plt.xticks(np.arange(2,7,1))
plt.title("Elbow graph with K values between 2 to 6", fontsize=15)
plt.xlabel("K values", fontsize=12)
plt.ylabel("Distance", fontsize=12)
plt.show()
kmeans_model_3 = KMeans(n_clusters=3).fit(X)
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6,6))
plt.scatter(x, y, s=0.5, c=kmeans_model_3.labels_)
plt.margins(0.001)
plt.xlabel('phi', fontsize=13)
plt.ylabel('psi', fontsize=13)
plt.xticks(np.arange(-180, 181, 30))
plt.yticks(np.arange(-180, 181, 30))
plt.scatter(kmeans_model_3.cluster_centers_[:,0], kmeans_model_3.cluster_centers_[:,1], c='red', marker='x')
plt.title('Data points and cluster centroids for k=3', fontsize=14)
plt.show()
The code below represents the shifting of the phi and psi angles in order to try and achieve better and more reasonable clustering. Upon experimenting with different angles to shift the phi and psi values by, the following values shown in the code seem to be the combination that provides a higher silhouette score and visually better clustering according to the corresponding plot. With the exception of some points, clusters are more cohesive now and there are seemingly fewer outliers towards the limits of the plot.
df_shifted = df.copy()
df_shifted['phi'] = df_shifted['phi'] + 210
df_shifted['phi'] = df_shifted['phi'] % 360
df_shifted['psi'] = df_shifted['psi'] + 110
df_shifted['psi'] = df_shifted['psi'] % 360
x_shift = np.array(df_shifted['phi'])[..., np.newaxis]
y_shift = np.array(df_shifted['psi'])[..., np.newaxis]
X_shift = np.hstack((x_shift, y_shift))
kmeans_shifted = KMeans(n_clusters=3).fit(X_shift)
centers_shift = kmeans_shifted.cluster_centers_
fig_kmeans_shifted = plt.figure(figsize=(7,7))
plt.scatter(x_shift, y_shift, s=0.3, c=kmeans_shifted.labels_.astype(float))
plt.margins(0.001)
plt.xlabel('phi', fontsize=12)
plt.ylabel('psi', fontsize=12)
plt.xticks(np.arange(0, 361, 30))
plt.yticks(np.arange(0, 361, 30))
plt.scatter(centers_shift[:,0], centers_shift[:,1], c='red', marker='x')
plt.title('Data points and cluster centroids with shifted angles', fontsize=14)
plt.savefig('clusters_shifted.png')
plt.show()
print("Silh. score for k=3 and shifted angles: ", metrics.silhouette_score(X_shift, kmeans_shifted.labels_, metric='euclidean'))
This decision for k=3 is further validated below by calculating the silhouette scores for the 3 k-means models with k=2, 3 and 4 respectively. With the silhouette score ranging from -1 to 1, a score as close to 1 is ideal to represent that the clustering is correct, and that datapoints are dense and close to their corresponding cluster center, without much if any overlap between clusters (this would reduce the score closer to 0). As shown below in the silhouette score prints, the score for k=3 is the highest one out of the 3, and is approaching 0.7, which shows a reasonable, fairly well-delimited clustering.
kmeans_model_2 = KMeans(n_clusters=2).fit(X_shift)
kmeans_model_3 = KMeans(n_clusters=3).fit(X_shift)
kmeans_model_4 = KMeans(n_clusters=4).fit(X_shift)
print("Silh. score for k=2: ", metrics.silhouette_score(X_shift, kmeans_model_2.labels_, metric='euclidean'))
print("Silh. score for k=3: ", metrics.silhouette_score(X_shift, kmeans_model_3.labels_, metric='euclidean'))
print("Silh. score for k=4: ", metrics.silhouette_score(X_shift, kmeans_model_4.labels_, metric='euclidean'))
np.random.seed(0)
size_remove = int(df.shape[0]*0.5)
fig, axes = plt.subplots(2,2, figsize=(12,12))
fig.suptitle("Stability on subsets for the K=3 Clustering", fontsize=17)
for ax in axes.ravel():
indices = np.random.choice(df_shifted.index, size=size_remove)
df_subset = df_shifted.drop(indices)
kmeans = KMeans(n_clusters=3, random_state=0).fit(df_subset[['phi', 'psi']])
ax.scatter(df_subset['phi'], df_subset['psi'], c=kmeans.labels_, s=0.5)
ax.set_xticks(np.arange(0, 361, 30))
ax.set_yticks(np.arange(0, 361, 30))
ax.set_xlabel("phi (°)")
ax.set_ylabel("psi (°)")
ax.margins(0.001)
plt.show()
from sklearn.neighbors import NearestNeighbors
from kneed import KneeLocator
from sklearn.metrics import silhouette_score
from sklearn.cluster import DBSCAN
df_phi_psi = df_shifted[['phi', 'psi']].copy()
df_phi_psi.head()
# Calculate the distance between each point and its 150 nearst neightbors
nearest_neighbors = NearestNeighbors(n_neighbors=150)
neighbors = nearest_neighbors.fit(df_phi_psi)
distances, indeces = neighbors.kneighbors(df_phi_psi)
distances = np.sort(distances[:,99], axis = 0)
# Sort distance values in ascending order
# Plot k-distance elbow plot
points = np.arange(len(distances))
knee = KneeLocator(points, distances, S=1, curve='convex', direction='increasing', interp_method='polynomial')
fig = plt.figure(figsize=(5, 5))
knee.plot_knee()
plt.xlabel("Points")
plt.ylabel("Distance")
plt.title("Elbow Point")
print("number of points is", len(distances))
print("elbow at distance",distances[knee.knee])
eps_list = np.arange(start=5, stop=21, step=1)
minpts_list = np.arange(start=100, stop=201, step=50)
df_silhouette_scores = pd.DataFrame()
for eps in eps_list:
for minpts in minpts_list:
db = DBSCAN(eps=eps, min_samples=minpts)
# limit to the case where DBSCAN clustered the data into 3 clusters (4 clusters including outliers)
if len(np.unique(db.fit_predict(df_phi_psi)))==4:
sil_score=silhouette_score(df_phi_psi, db.fit_predict(df_phi_psi))
else:
continue
trial_parameters="eps:" + str(eps) +" min_sample :" + str(minpts)
df_silhouette_scores=df_silhouette_scores.append(pd.DataFrame(data=[[np.unique(db.fit_predict(df_phi_psi)), sil_score,trial_parameters]], columns=["cluster", "score", "parameters"]))
df_silhouette_scores.sort_values(by='score', ascending=False).head()
# The best hyperparameters with highest Score
df_silhouette_scores.sort_values(by='score', ascending=False).head(1)
db = DBSCAN(eps=17, min_samples=200).fit(df_phi_psi)
core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
core_samples_mask[db.core_sample_indices_] = True
labels = db.labels_
# Number of clusters in labels, ignoring noise if present.
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
n_noise = list(labels).count(-1)
fig_DBSCAN = plt.figure(figsize=(7,7))
plt.scatter(df_phi_psi['phi'], df_phi_psi['psi'], c=labels, s=0.3)
plt.margins(0.001)
plt.xticks(np.arange(0,361,30))
plt.yticks(np.arange(0,361,30))
plt.xlabel("phi", fontsize=13)
plt.ylabel("psi", fontsize=13)
plt.title('Clusters found by DBSCAN', fontsize=15)
plt.show()
print('Estimated number of clusters is %d' % n_clusters)
print('Estimated number of outliers is %d' % n_noise)
fig, axes = plt.subplots(1,2, figsize=(13,6))
axes[0].scatter(x_shift, y_shift, s=0.3, c=kmeans_shifted.labels_.astype(float))
axes[0].margins(0.001)
axes[0].set_xlabel('phi', fontsize=13)
axes[0].set_ylabel('psi', fontsize=13)
axes[0].set_xticks(np.arange(0, 361, 30))
axes[0].set_yticks(np.arange(0, 361, 30))
axes[0].scatter(centers_shift[:,0], centers_shift[:,1], c='red', marker='x')
axes[0].set_title('Clusters found by K-means', fontsize=15)
axes[1].scatter(df_phi_psi['phi'], df_phi_psi['psi'], c=labels, s=0.3)
axes[1].margins(0.001)
axes[1].set_xticks(np.arange(0,361,30))
axes[1].set_yticks(np.arange(0,361,30))
axes[1].set_xlabel("phi", fontsize=13)
axes[1].set_ylabel("psi", fontsize=13)
axes[1].set_title('Clusters found by DBSCAN', fontsize=15)
plt.show()
eps_diff_list = np.arange(-3, 4, 1)
minpts_diff_list = np.arange(-30, 40, 10)
matrix = np.empty((len(eps_diff_list), len(minpts_diff_list)))
for i, eps_diff in enumerate(eps_diff_list):
for j, minpts_diff_diff in enumerate(minpts_diff_list):
db = DBSCAN(eps=17+eps_diff, min_samples=200+minpts_diff_diff).fit(df_phi_psi)
labels = db.labels_
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
matrix[i][j] = n_clusters
plt.figure(figsize=(7,6))
sns.heatmap(matrix, annot=True, cmap='viridis', linewidths=0.01, cbar_kws={"ticks":np.arange(2,5,1)})
plt.xticks(np.arange(0,7,1)+0.5, minpts_diff_list, rotation=0)
plt.yticks(np.arange(0,7,1)+0.5, eps_diff_list, rotation=0)
plt.title("Number of Clusters with different eps and minpts", fontsize=14)
plt.xlabel("Difference in minpts compared to mpts=200", fontsize=13)
plt.ylabel("Difference in minpts comppared to eps=17", fontsize=13)
plt.show()