{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Just a simple EDA and clustering embedding ","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n        \nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-05-31T08:15:06.431533Z","iopub.execute_input":"2023-05-31T08:15:06.431955Z","iopub.status.idle":"2023-05-31T08:15:06.493141Z","shell.execute_reply.started":"2023-05-31T08:15:06.431921Z","shell.execute_reply":"2023-05-31T08:15:06.49188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\nfrom Bio import SeqIO\nfrom tqdm import tqdm\n\ndef read_fasta(fastaPath):    \n    fasta_sequences = SeqIO.parse(open(fastaPath), 'fasta')\n    ids = []\n    sequences = []\n    for fasta in fasta_sequences:\n        ids.append(fasta.id)\n        sequences.append(str(fasta.seq))\n    return pd.DataFrame({'Id': ids, 'Sequence': sequences})\n\ndef get_top_go_terms(data, num_terms):\n    term_counts = data['term'].value_counts()\n    freq_counts = term_counts / len(data)\n    freq_top = freq_counts.nlargest(num_terms)\n    return freq_top\n\ntrain_terms = pd.read_csv('/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv', sep='\\t')\ntrain_sequences = read_fasta('/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta')\ntop_terms = get_top_go_terms(train_terms, 10)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:06.495313Z","iopub.execute_input":"2023-05-31T08:15:06.495821Z","iopub.status.idle":"2023-05-31T08:15:14.500603Z","shell.execute_reply.started":"2023-05-31T08:15:06.495792Z","shell.execute_reply":"2023-05-31T08:15:14.499652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_terms","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:14.502047Z","iopub.execute_input":"2023-05-31T08:15:14.502473Z","iopub.status.idle":"2023-05-31T08:15:14.51367Z","shell.execute_reply.started":"2023-05-31T08:15:14.502429Z","shell.execute_reply":"2023-05-31T08:15:14.512561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sequences.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:14.517029Z","iopub.execute_input":"2023-05-31T08:15:14.517565Z","iopub.status.idle":"2023-05-31T08:15:14.546067Z","shell.execute_reply.started":"2023-05-31T08:15:14.517522Z","shell.execute_reply":"2023-05-31T08:15:14.544799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sequences['n_caracs'] = train_sequences['Sequence'].str.len()\ntrain_sequences['n_caracs'].describe()","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:14.54762Z","iopub.execute_input":"2023-05-31T08:15:14.548087Z","iopub.status.idle":"2023-05-31T08:15:14.70198Z","shell.execute_reply.started":"2023-05-31T08:15:14.548047Z","shell.execute_reply":"2023-05-31T08:15:14.700783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"q_low = train_sequences[\"n_caracs\"].quantile(0.01)\nq_hi  = train_sequences[\"n_caracs\"].quantile(0.99)\n\ndf_filtered = train_sequences[(train_sequences[\"n_caracs\"] < q_hi) & (train_sequences[\"n_caracs\"] > q_low)]","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:14.703466Z","iopub.execute_input":"2023-05-31T08:15:14.703899Z","iopub.status.idle":"2023-05-31T08:15:14.726573Z","shell.execute_reply.started":"2023-05-31T08:15:14.703869Z","shell.execute_reply":"2023-05-31T08:15:14.725323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_filtered[['Id', 'n_caracs']].plot.hist(bins= 40)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:14.727973Z","iopub.execute_input":"2023-05-31T08:15:14.728303Z","iopub.status.idle":"2023-05-31T08:15:15.286381Z","shell.execute_reply.started":"2023-05-31T08:15:14.728274Z","shell.execute_reply":"2023-05-31T08:15:15.285277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# results = []\n# for index, row in tqdm(test_data.iterrows(), total=test_data.shape[0], position=0):\n#     for term, freq in top_terms.items():\n#         results.append((row['Id'], term, freq))","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:15.288275Z","iopub.execute_input":"2023-05-31T08:15:15.288615Z","iopub.status.idle":"2023-05-31T08:15:15.293184Z","shell.execute_reply.started":"2023-05-31T08:15:15.288586Z","shell.execute_reply":"2023-05-31T08:15:15.291897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"TEST_SIZE = 0.01","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:16:45.806895Z","iopub.execute_input":"2023-05-31T08:16:45.807558Z","iopub.status.idle":"2023-05-31T08:16:45.812949Z","shell.execute_reply.started":"2023-05-31T08:16:45.807509Z","shell.execute_reply":"2023-05-31T08:16:45.811735Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class search_optimal_cluster():\n    def __init__(self, train_embeds, search_range = range(2, 10)):\n        self.search_range = search_range\n        self.train_embeds = train_embeds\n        self.PCA_n_components = 2\n        self.random_state = 1993\n        \n    def do_PCA(self):\n        pca = PCA(n_components = self.PCA_n_components)\n        pca.fit(self.train_embeds)\n        X = pca.transform(self.train_embeds)\n        print(np.sum(pca.explained_variance_ratio_))\n        \n        return X\n    \n    def calculate_clusters(self, best_k = 2):\n        '''\n        '''\n        self.df_PCA = pd.DataFrame(data=self.do_PCA(), columns=['x1', 'x2'])\n        self.df = pd.DataFrame(self.train_embeds)\n        \n        best_k = best_k\n        best_score = 0\n        sse = []\n        report = {}\n        for k in self.search_range:\n            temp_dict = {}\n            print('***********************************')\n            print(' K: ', k)\n            kmeans = KMeans(init='k-means++',\n                            algorithm='auto', n_clusters=k)\n            labels = kmeans.fit_predict(self.df_PCA)\n            centroids = kmeans.cluster_centers_\n\n            gc.collect()\n\n            # elbow\n            inertia = kmeans.inertia_\n            temp_dict['Sum of squared error'] = inertia\n            chs = calinski_harabasz_score(self.df_PCA, labels)\n            score = silhouette_score(self.df_PCA, labels)\n            if score > best_score:\n                best_k = k\n                best_score = score\n            print(\"For n_clusters = {}, silhouette score is {})\".format(k, score))\n            temp_dict['Calinski Harabasz Score'] = chs\n            temp_dict['Silhouette Score'] = score\n            report[k] = temp_dict\n\n            gc.collect()\n            \n        report_df = pd.DataFrame(report).T\n        display(report_df)\n        self.report = report_df\n\n    def visualize_report(self):\n        '''\n        '''\n        self.report.plot(figsize=(5, 5),\n               xticks=self.search_range,\n               grid=True,\n               title=f'Selecting optimal \"K\"',\n               subplots=True,\n               marker='o',\n               sharex=True)\n        plt.tight_layout()\n        \n    def kelbow_visualizer(self):\n        '''\n        '''\n        from yellowbrick.cluster.elbow import kelbow_visualizer\n\n        kelbow_visualizer(KMeans(random_state= self.random_state),\n                          self.df,\n                          k=self.search_range,\n                          timings=False)\n        \n    def intercluster_distance(self, n = 6):\n        '''\n        '''\n        from yellowbrick.cluster import intercluster_distance\n\n        intercluster_distance(KMeans(n_clusters = n), \n                      self.df_PCA, \n                      embedding='mds', \n                      random_state= self.random_state) # other option for embedding 'tsne'\n        \n    def save_best(self, model_name = '', k = 6):\n        '''\n        '''\n        # \n        n_clusters = k #best_k\n        kmeans = KMeans(n_clusters = n_clusters)\n        kmeans.fit(self.df)\n        labels = kmeans.predict(self.df)\n        centroids = kmeans.cluster_centers_\n        \n        plt.scatter(self.df_PCA.iloc[:, 0], self.df_PCA.iloc[:, 1], c=labels, s=50, cmap='viridis')\n        plt.scatter(centroids[:, 0], centroids[:, 1], c='black', s=200, alpha=0.5);\n        \n        # apply test split\n        IX = np.arange(len(train_embeds))\n        IX_test, IX_train, _, _ = train_test_split( IX, IX, train_size= TEST_SIZE, random_state=42, stratify = labels)\n        \n        np.save('cv_groups_' + model_name + '.npy', labels)\n        np.save('train_split_' + model_name + '.npy', IX_train)\n        np.save('test_split_' + model_name + '.npy', IX_test)","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:27.301406Z","iopub.execute_input":"2023-05-31T08:15:27.301853Z","iopub.status.idle":"2023-05-31T08:15:27.328262Z","shell.execute_reply.started":"2023-05-31T08:15:27.301811Z","shell.execute_reply":"2023-05-31T08:15:27.326804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## T5 ##","metadata":{}},{"cell_type":"code","source":"train_embeds = np.load('/kaggle/input/t5embeds/train_embeds.npy')\ntrain_embeds.shape","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:15.311839Z","iopub.execute_input":"2023-05-31T08:15:15.31225Z","iopub.status.idle":"2023-05-31T08:15:26.244955Z","shell.execute_reply.started":"2023-05-31T08:15:15.312218Z","shell.execute_reply":"2023-05-31T08:15:26.244052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ids= np.load('/kaggle/input/t5embeds/train_ids.npy')\ntrain_ids.shape","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:26.24638Z","iopub.execute_input":"2023-05-31T08:15:26.246979Z","iopub.status.idle":"2023-05-31T08:15:26.302492Z","shell.execute_reply.started":"2023-05-31T08:15:26.246946Z","shell.execute_reply":"2023-05-31T08:15:26.30123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nfrom sklearn.cluster import KMeans\nfrom sklearn.metrics import silhouette_score, calinski_harabasz_score\nfrom sklearn.model_selection import train_test_split\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nimport gc","metadata":{"execution":{"iopub.status.busy":"2023-05-31T08:15:26.304219Z","iopub.execute_input":"2023-05-31T08:15:26.305379Z","iopub.status.idle":"2023-05-31T08:15:27.299697Z","shell.execute_reply.started":"2023-05-31T08:15:26.305334Z","shell.execute_reply":"2023-05-31T08:15:27.29849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nT5_search = search_optimal_cluster(train_embeds, search_range = range(2, 10))\nT5_search.calculate_clusters()","metadata":{"execution":{"iopub.status.busy":"2023-05-26T12:36:18.269289Z","iopub.execute_input":"2023-05-26T12:36:18.269724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"T5_search.visualize_report()","metadata":{"execution":{"iopub.status.busy":"2023-05-26T12:36:07.624324Z","iopub.execute_input":"2023-05-26T12:36:07.624726Z","iopub.status.idle":"2023-05-26T12:36:08.331951Z","shell.execute_reply.started":"2023-05-26T12:36:07.624697Z","shell.execute_reply":"2023-05-26T12:36:08.330737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"T5_search.kelbow_visualizer()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"T5_search.intercluster_distance()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"T5_search.save_best(model_name = 't5')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## EMS_t36  esm2_t36_3B, embeding size = 2560","metadata":{}},{"cell_type":"code","source":"train_embeds = np.load('/kaggle/input/4637427/train_embeds_esm2_t36_3B_UR50D.npy')\ntrain_embeds.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nEMS_t36_search = search_optimal_cluster(train_embeds, search_range = range(2, 10))\nEMS_t36_search.calculate_clusters()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t36_search.visualize_report()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t36_search.kelbow_visualizer()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t36_search.intercluster_distance()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t36_search.save_best(model_name = 'ems_t36')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EMS esm2_t33_650M, embeding size = 1280","metadata":{}},{"cell_type":"code","source":"train_embeds = np.load('/kaggle/input/cafa-5-ems-2-embeddings-numpy/train_embeddings.npy')\ntrain_embeds.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nEMS_t33_search = search_optimal_cluster(train_embeds, search_range = range(2, 10))\nEMS_t33_search.calculate_clusters()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t33_search.visualize_report()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t33_search.kelbow_visualizer()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t33_search.intercluster_distance()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"EMS_t33_search.save_best(model_name = 'ems_t33')","metadata":{},"execution_count":null,"outputs":[]}]}