\n\nThe Critical Assessment of Functional Annotation (CAFA) is an ongoing global, community-wide effort designed to provide a large-scale assessment of computational methods dedicated to predicting protein function. \nProteins are responsible for many activities in our tissues, organs, and bodies and they also play a central role in the structure and function of cells. Proteins are large molecules composed of 20 types of building-blocks known as amino acids. The human body makes tens of thousands of different proteins, and each protein is composed of dozens or hundreds of amino acids that are linked sequentially. This amino-acid sequence determines the 3D structure and conformational dynamics of the protein, and that, in turn, determines its biological function. The accurate assignment of biological function to the protein is the key to understanding life at the molecular level. More knowledge of the functions assigned to proteins could lead to curing diseases and improving human and animal health and wellness in areas as varied as medicine and agriculture. \n
\n \n
\n3D structure of the protein myoglobin showing turquoise α-helices. This protein was the first to have its structure solved by X-ray crystallography \n\n This kernel has been influenced by great public kernels posted in this contest. Thank you all fellow kagglers for sharing all this great content, very helpful especially for people like me that do not have a strong background in the science of biology. If you find something interesting please do not forget to UPVOTE! ","metadata":{}},{"cell_type":"code","source":"!pip install --force-reinstall --no-deps numpy==1.22.3 \n","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:35:24.614459Z","iopub.execute_input":"2023-08-18T17:35:24.615622Z","iopub.status.idle":"2023-08-18T17:35:31.293306Z","shell.execute_reply.started":"2023-08-18T17:35:24.615573Z","shell.execute_reply":"2023-08-18T17:35:31.291296Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install --force-reinstall --no-deps goatools==1.1.6","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:35:31.29573Z","iopub.execute_input":"2023-08-18T17:35:31.296111Z","iopub.status.idle":"2023-08-18T17:36:01.136558Z","shell.execute_reply.started":"2023-08-18T17:35:31.296083Z","shell.execute_reply":"2023-08-18T17:36:01.135738Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 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\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n for filename in filenames:\n print(os.path.join(dirname, filename))","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:01.137698Z","iopub.execute_input":"2023-08-18T17:36:01.138043Z","iopub.status.idle":"2023-08-18T17:36:01.158066Z","shell.execute_reply.started":"2023-08-18T17:36:01.138011Z","shell.execute_reply":"2023-08-18T17:36:01.157342Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import goatools\nfrom goatools.semantic import TermCounts, get_info_content\nfrom goatools.associations import dnld_assc\nfrom goatools import obo_parser\nfrom goatools.obo_parser import GODag\nfrom goatools.semantic import resnik_sim\nfrom IPython.display import Image, display","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:01.160623Z","iopub.execute_input":"2023-08-18T17:36:01.161066Z","iopub.status.idle":"2023-08-18T17:36:01.185198Z","shell.execute_reply.started":"2023-08-18T17:36:01.161028Z","shell.execute_reply":"2023-08-18T17:36:01.183863Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.cm as cm\nfrom sklearn.cluster import KMeans\nfrom sklearn.metrics import silhouette_samples, silhouette_score,roc_auc_score\nfrom sklearn import preprocessing,metrics\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.decomposition import PCA\nfrom sklearn.manifold import trustworthiness as trust \nfrom yellowbrick.cluster.elbow import KElbowVisualizer\nfrom termcolor import colored\nimport seaborn as sns\nimport gc","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:01.186861Z","iopub.execute_input":"2023-08-18T17:36:01.187193Z","iopub.status.idle":"2023-08-18T17:36:02.772631Z","shell.execute_reply.started":"2023-08-18T17:36:01.187163Z","shell.execute_reply":"2023-08-18T17:36:02.771279Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pip install pygraphviz ","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:02.774993Z","iopub.execute_input":"2023-08-18T17:36:02.775731Z","iopub.status.idle":"2023-08-18T17:36:16.90286Z","shell.execute_reply.started":"2023-08-18T17:36:02.775686Z","shell.execute_reply":"2023-08-18T17:36:16.901164Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ** Training data **\n\nFor this contest, we are given protein sequences on which we have to predict Gene Ontology (GO) terms in each of the three subontologies: Molecular Function (MF), Biological Process (BP), and Cellular Component (CC). Gene Ontology is used to annotate proteins and describes molecules functions in the style of a directed acyclic graph (DAG). A very helpful guide for Gene Ontology can be found here [2]. \nThe nodes in this graph are functional descriptors (terms or classes) connected by relational ties between them (is_a, part_of, etc.). For example, terms 'protein binding activity' and 'binding activity' are related by an is_a relationship; however, the edge in the graph is often reversed to point from binding towards protein binding. The protein's function is therefore represented by a subset of one or more of the subontologies. We ll use GOATools,a library for GO term manipulation, GOEA testing, and custom ontology visualization, to load our GO DAG and plot a random GO term up to root.","metadata":{}},{"cell_type":"code","source":"godag = GODag(\"/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\")\nfin_gaf = os.path.join(os.getcwd(), \"tair.gaf\")\nassociations = dnld_assc(fin_gaf, godag)\ntermcounts = TermCounts(godag, associations) \ngo_id2 = 'GO:0097192'\nrec = godag[go_id2]\nlineage_png = '/kaggle/working/GO_lineage.png'\ngodag.draw_lineage([rec] )\ndisplay(Image(lineage_png, width=600, height=300))","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:16.904805Z","iopub.execute_input":"2023-08-18T17:36:16.905182Z","iopub.status.idle":"2023-08-18T17:36:29.421981Z","shell.execute_reply.started":"2023-08-18T17:36:16.905145Z","shell.execute_reply":"2023-08-18T17:36:29.421082Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"GO terms are generally defined to be taxon neutral. As a result of research experiments, scientists might discover new biological processes, locations or functions that a protein performs. When this happens proteins are annotated with the corresponding GO term and these annotations are then incorporated to the many protein databases such as UniProtKB and UniProt-GOA. Our training data file \"Train_terms.tsv\" contains protein sequences that have at least one experimentally determined GO term in at least one subontology with relevant annotations in three columns: \n\n- **EntryID** is the protein's UniProt accession ID\n- **term** is the GO term ID\n- **aspect** indicates in which ontology the term appears\n-------------------------------------------------------------","metadata":{}},{"cell_type":"code","source":"Train_terms = pd.read_csv('../input/cafa-5-protein-function-prediction/Train/train_terms.tsv',sep=\"\\t\")\nprint(\"\\033[1m Training data available :\",Train_terms.shape[0])\nNumProteins=Train_terms[\"EntryID\"].unique()\npd.DataFrame(Train_terms.head(5))","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:29.423415Z","iopub.execute_input":"2023-08-18T17:36:29.42373Z","iopub.status.idle":"2023-08-18T17:36:32.596967Z","shell.execute_reply.started":"2023-08-18T17:36:29.423705Z","shell.execute_reply":"2023-08-18T17:36:32.595527Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(6, 3))\nax = sns.countplot(data=Train_terms, x='aspect')\nax.set_yticks(ax.get_yticks().astype(int)) \nax.set_yticklabels(ax.get_yticks())\nplt.title('Go Domains')\nplt.ylabel('No. of Occurances')","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:32.598243Z","iopub.execute_input":"2023-08-18T17:36:32.599218Z","iopub.status.idle":"2023-08-18T17:36:35.094551Z","shell.execute_reply.started":"2023-08-18T17:36:32.599167Z","shell.execute_reply":"2023-08-18T17:36:35.093043Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"These are the three ways that a gene functions :\n- **BPO** is the larger process, or ‘biological programs’ accomplished by multiple molecular activities\n- **CCO** location relative to the cell where activity takes place\n- **MFO** molecular-level process or activity the gene carries out","metadata":{}},{"cell_type":"code","source":"print('\\033[1m \\033[95mUnique proteins in training data:',len(NumProteins))\nGoterms=Train_terms[\"term\"].unique()\nprint('\\033[92m Unique Go Terms in training data:',len(Goterms))","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:35.099935Z","iopub.execute_input":"2023-08-18T17:36:35.100336Z","iopub.status.idle":"2023-08-18T17:36:35.519347Z","shell.execute_reply.started":"2023-08-18T17:36:35.100303Z","shell.execute_reply":"2023-08-18T17:36:35.517789Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ** Go-terms selection **\nWe have to select only some of the 31466 Go terms to work with , with max limit to 1500 as contest suggests. Selecting Go terms based on their frequency is a method widely used in previous CAFA contests and will be used in this solution but first lets see how are GO terms assigned to target proteins.","metadata":{"execution":{"iopub.status.busy":"2023-07-26T07:23:09.976554Z","iopub.execute_input":"2023-07-26T07:23:09.977111Z","iopub.status.idle":"2023-07-26T07:23:09.985676Z","shell.execute_reply.started":"2023-07-26T07:23:09.977069Z","shell.execute_reply":"2023-07-26T07:23:09.983967Z"}}},{"cell_type":"code","source":"Gos_per_prot = Train_terms.groupby('EntryID')['term'].apply(list).reset_index()\nGos_per_prot.columns = ['EntryID', 'term']\nGos_per_prot['Go_Terms_Count'] = Gos_per_prot['term'].apply(len)\nprint('\\033[1m Go terms related to each protein')\nGos_per_prot","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:35.52113Z","iopub.execute_input":"2023-08-18T17:36:35.521547Z","iopub.status.idle":"2023-08-18T17:36:39.872566Z","shell.execute_reply.started":"2023-08-18T17:36:35.521513Z","shell.execute_reply":"2023-08-18T17:36:39.871659Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(4,3))\nn, bins, patches = plt.hist(x=Gos_per_prot.Go_Terms_Count, bins='auto',alpha=0.7, rwidth=0.85)\nplt.xlim(right=150, left=0)\nplt.grid(axis='y', alpha=0.75)\nplt.xlabel('Number of Go terms')\nplt.ylabel('Proteins')\nplt.title('Go Terms per Protein')","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:39.874027Z","iopub.execute_input":"2023-08-18T17:36:39.874567Z","iopub.status.idle":"2023-08-18T17:36:40.766626Z","shell.execute_reply.started":"2023-08-18T17:36:39.874536Z","shell.execute_reply":"2023-08-18T17:36:40.765184Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Go_freqs = (Train_terms['term'].value_counts())\ntargetGos=list(Go_freqs.index[:1500]) ","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:36:40.768255Z","iopub.execute_input":"2023-08-18T17:36:40.768588Z","iopub.status.idle":"2023-08-18T17:36:41.270846Z","shell.execute_reply.started":"2023-08-18T17:36:40.768558Z","shell.execute_reply":"2023-08-18T17:36:41.269297Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Embeddings
\n
\n\nProtein language models enable the encoding of protein sequences in \nnumerical vectors so they can be used in machine learning. The overall idea of a using a protein language model is to \"translate\" protein sequences as sentences and amino acids as words in order to extract useful information and representation. These sequence representations are called embeddings and they are used to predict global or local protein properties. There are many models available like Esm/Esm1b, ProtT5, and SeqVec with different characteristics and performance [3], this kernel uses ProtT5 model embeddings thanks to Sergei Fironov's work, available to download from \nhttps://www.kaggle.com/datasets/sergeifironov/t5embeds. ","metadata":{}},{"cell_type":"code","source":"X = np.load('/kaggle/input/t5embeds/train_embeds.npy')\ntrain_ids = np.load('/kaggle/input/t5embeds/train_ids.npy')\ncolumn_names_A = ['Protein_Emb_' + str(i) for i in range(X.shape[1])]\nEmb_df =pd.DataFrame(X, columns=column_names_A)\nprint('\\033[1m Statistics of protein embeddings from ProtT5 model ')\nEmb_df.describe().T","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:54:59.503243Z","iopub.execute_input":"2023-08-18T17:54:59.503694Z","iopub.status.idle":"2023-08-18T17:55:09.642505Z","shell.execute_reply.started":"2023-08-18T17:54:59.50366Z","shell.execute_reply":"2023-08-18T17:55:09.641016Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Biopython
\n
\n\nWe are given two files , train_sequences.fasta which contains the protein sequences for the training dataset and testsuperset.fasta that contains protein sequences on which predictions will be estimated. The header for each sequence contains the protein's UniProt accession ID and the Taxon ID of the species this protein belongs to. We shall use Biopython , a set of freely available tools for biological computation. Bio.SeqIO provides a simple uniform interface to input and output assorted sequence file formats as SeqRecord objects. In bioinformatics and biochemistry, the FASTA format is a text-based format for representing either nucleotide sequences or amino acid (protein) sequences, in which nucleotides or amino acids are represented using single-letter codes. We can use information from .fasta in many ways , for example to create embeddings, study sequence lengths or as we ll see later to extract protein sequences for similarity detection in formed protein clusters.","metadata":{"execution":{"iopub.status.busy":"2023-07-31T14:31:24.476094Z","iopub.execute_input":"2023-07-31T14:31:24.476458Z","iopub.status.idle":"2023-07-31T14:31:24.485156Z","shell.execute_reply.started":"2023-07-31T14:31:24.476409Z","shell.execute_reply":"2023-07-31T14:31:24.482999Z"}}},{"cell_type":"code","source":"from Bio import SeqIO\n\npathbio='/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'\nsequences = SeqIO.parse(pathbio, \"fasta\")\nseq_list = []\nfor seq in sequences:\n seq_list.append(seq)\n \nprint('\\033[1m \\033[91m There are ', {len(seq_list)}, 'protein Sequences in trainining .fasta file')\nprint('\\033[1m \\033[94m' ,f'\\n Sample Sequence: \\n{seq_list[20]}')\nfrom Bio.SeqIO.FastaIO import SimpleFastaParser\nwith open(pathbio) as fasta_file: \n identifiers = []\n lengths = []\n seqs =[]\n for title, sequence in SimpleFastaParser(fasta_file):\n identifiers.append(title.split(None, 1)[0]) # First word is ID\n lengths.append(len(sequence))\n seqs.append(sequence)\nIdsSeq=pd.DataFrame({'EntryID1':identifiers, 'Sequences' : seqs, 'length':lengths})\nprint('\\033[1m \\033[92m ')\nprint(' Protein Ids and sequence lengths')\nprint('\\033[1m \\033[92m', IdsSeq[:5])","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:37:02.777277Z","iopub.execute_input":"2023-08-18T17:37:02.777532Z","iopub.status.idle":"2023-08-18T17:37:07.543965Z","shell.execute_reply.started":"2023-08-18T17:37:02.777508Z","shell.execute_reply":"2023-08-18T17:37:07.542781Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plt.figure(figsize=(3,3))\nn, bins, patches = plt.hist(x=IdsSeq.length, bins='auto',alpha=0.7, rwidth=0.85)\nplt.xlim(right=2500, left=0)\nplt.grid(axis='y', alpha=0.75)\nplt.xlabel('Length')\nplt.ylabel('Frequency')\nplt.title('Protein Sequence Lengths')\n","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:37:07.54518Z","iopub.execute_input":"2023-08-18T17:37:07.545514Z","iopub.status.idle":"2023-08-18T17:37:10.631763Z","shell.execute_reply.started":"2023-08-18T17:37:07.545484Z","shell.execute_reply":"2023-08-18T17:37:10.630301Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Clustering Embeddings
\n
\nEmbeddings are representations of protein sequences (composition of amino acids and order in the sequence) so clustering them is expected to group sequences in such a way that each cluster has proteins that share similar features and functionality. This sequence similarity can indicate ”homologous” proteins or genes that share common ancestry. The main idea of this work is that GO terms that appear in a cluster a specific cluster might be assigned on uknown sequences that belong to the same cluster, transfering successfully relevant predicted function annotations. The results of course depend on the quality of clustering. \n\n## ** Dimensionality Reduction **\n\nOur embeddings have high dimensionality so trying to apply clustering analysis we 'll face the problem of the “Curse of Dimensionality”. At so high dimensions measures like as Euclidean and Manhattan distances , used for clustering become less effective. Moreover dimensionality reduction before clustering is a technique expected to improve final results since it removes redundant features, noise and irrelevant data. PCA and UMAP are two well known techniques that will be used before our clustering methods.\n## ** K-means **\n K-Means clustering is a very polular unsupervised learning algorithm approach that aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean (cluster centers or cluster centroid), serving as a prototype of the cluster. We ll use the Elbow method to determine the optimal number of clusters, always keeping in mind that K-means faces problems with clusters of unequal sizes and densities and also outliers. The silhouette coefficient will be used to analyze the quality of our clusters. So let's start with PCA to reduce data dimensionality keeping 95% of the variance.\n ","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\nscaler = StandardScaler()\nscaler.fit(Emb_df)\nEmb_dfstd = scaler.transform(Emb_df)\npca = PCA(n_components= 0.95 ,svd_solver = \"full\",random_state = 20) \ncomponents_PCA=pca.fit_transform(Emb_dfstd)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:40:11.717377Z","iopub.execute_input":"2023-08-18T17:40:11.71772Z","iopub.status.idle":"2023-08-18T17:40:34.838261Z","shell.execute_reply.started":"2023-08-18T17:40:11.717696Z","shell.execute_reply":"2023-08-18T17:40:34.836577Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":" #### ** Elbow Method **\n The KElbowVisualizer implements the “Elbow” method to select the optimal number of clusters by fitting the kmeans model with a range of values for K. The \"Elbow\" method is a a heuristic method that tests various values for K and calculates WCSS the sum of the squared distance between each point and the centroid in a cluster.","metadata":{}},{"cell_type":"code","source":"k_values = range(10,300)\nclusterer = KMeans(init='k-means++', n_init='auto' ,\n max_iter =500, random_state=21)\nvisualizer = KElbowVisualizer(clusterer, k=k_values)\nvisualizer.fit(components_PCA[:40000])\nvisualizer.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:40:47.312909Z","iopub.execute_input":"2023-08-18T17:40:47.313297Z","iopub.status.idle":"2023-08-18T17:40:48.062014Z","shell.execute_reply.started":"2023-08-18T17:40:47.313266Z","shell.execute_reply":"2023-08-18T17:40:48.060756Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#### ** Silhouette score **\n\nSilhouette Coefficient is a metric used to calculate the quality of a clustering technique. It is a measure of how similar an object is to its own cluster (cohesion) compared to other clusters (separation), giving an average value for all the samples. The silhouette ranges from −1 to +1, where a high value indicates that the object is well matched to its own cluster and poorly matched to neighboring clusters. If most objects have a high value, then the clustering configuration is appropriate. If many points have a low or negative value, then the clustering configuration may have too many or too few clusters.","metadata":{}},{"cell_type":"code","source":"n_clusters = 54\nplt.figure(figsize=(5, 4))\nax = plt.gca()\nax.set_xlim([-0.1, 1])\nax.set_ylim([0, len(components_PCA) + (n_clusters + 1) * 10])\nclusterer = KMeans(n_clusters=n_clusters,init='k-means++', n_init='auto' ,\n max_iter =300, random_state=21)\ncluster_labels = clusterer.fit_predict(components_PCA)\nsilhouette_avg = silhouette_score(components_PCA, cluster_labels)\nprint(\"For\",n_clusters,\"clusters the average silhouette_score is :\",silhouette_avg, )\nprint(' ')\nprint(' ')\nsample_silhouette_values = silhouette_samples(components_PCA, cluster_labels)\ny_lower = 10\nfor i in range(n_clusters):\n ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]\n ith_cluster_silhouette_values.sort()\n size_cluster_i = ith_cluster_silhouette_values.shape[0]\n y_upper = y_lower + size_cluster_i\n color = cm.nipy_spectral(float(i) / n_clusters)\n ax.fill_betweenx(np.arange(y_lower, y_upper), 0,ith_cluster_silhouette_values, facecolor=color, edgecolor=color,alpha=0.7,)\n ax.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))\n y_lower = y_upper + 10 # 10 for the 0 samples\nax.set_xlabel(\"The silhouette coefficient values\")\nax.set_ylabel(\"Cluster number\")\nax.axvline(x=silhouette_avg, color=\"red\", linestyle=\"--\")\nax.set_yticks([]) # Clear the yaxis labels / ticks\nax.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])\nplt.show()","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"As we see clusters have many outliers and overall quality is low, experimenting with K values did not improve results. ","metadata":{}},{"cell_type":"code","source":"#no_clusters =visualizer.elbow_value_\nno_clusters =54\nclusterer = KMeans(n_clusters= no_clusters, max_iter =500, n_init='auto', random_state=21)\nlabels= clusterer.fit_predict(components_PCA)\nEmb_df[\"ClustersK\"] = labels\nEmb_df['seq']=IdsSeq.Sequences\nEmb_df['EntryID1']=IdsSeq.EntryID1\nunique_values, counts1 = np.unique(labels, return_counts=True)\nplt.figure(figsize=(6,4))\nplt.bar(unique_values, counts1)\nplt.xlabel('K-Means Clusters')\nplt.ylabel('Proteins in Cluster')\nplt.ylim(0, 5500)\nplt.title('Protein distriblution in clusters for K-Means')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:40:54.095047Z","iopub.execute_input":"2023-08-18T17:40:54.095407Z","iopub.status.idle":"2023-08-18T17:40:55.501605Z","shell.execute_reply.started":"2023-08-18T17:40:54.095378Z","shell.execute_reply":"2023-08-18T17:40:55.499833Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## ** HDBSCAN **\nDBSCAN – is a density-based clustering non-parametric algorithm: given a set of points in some space, it groups together points that are closely packed together (points with many nearby neighbors), marking as outliers points in low-density regions.We shall use HDBSCAN,a modification of DBSCAN that allows clusters WITH varying density . The metric used to evaluate the quality of results is Density Based Clustering Validation (DBCV) that considers both density and shape properties of cluster. But first,let's start with UMAP, a dimensionality reduction technique that preserves data's global structure. UMAP scales efficiently in terms of both dataset size and dimensionality,improving clustering algorithms [4]. It can be controlled using n_components for the number of dimensions after the projection and min_dist the minimum distance apart that points are allowed to be in the low dimensional representation. A range of values between (10-50) for n_components was tested using Trustworthiness as a metric. Based on data complexity and volume n_components was set to 15. \n#### ** Trustworthiness criterion **\n\nThis metric (Jarkko Venna and Samuel Kaski. 2001) is used to indicate to what extent the local structure is retained after applying UMAP. Values >0.90 are considered sufficient. Due to memory limitations it was calculated over 30000 samples of ProtT5 embeddings. Since we do not want to focus on local structure k=15 is selected.","metadata":{}},{"cell_type":"code","source":"import umap \nfrom umap import validation \nXum=Emb_df.drop(['ClustersK','seq','EntryID1'], axis=1)\nclust = [5,7,10,15,20,25,30,35,40,45,50]\nfor K in clust: \n um = umap.UMAP(n_components=K) \n components_umap = um.fit_transform(Xum[:30000 ]) \n trust_score = trust(components_umap[:30000 ],Xum[:30000]) \n print('K=',K,'=',trust_score)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:41:00.017793Z","iopub.execute_input":"2023-08-18T17:41:00.018145Z","iopub.status.idle":"2023-08-18T17:42:00.397445Z","shell.execute_reply.started":"2023-08-18T17:41:00.018116Z","shell.execute_reply":"2023-08-18T17:42:00.396513Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"n_components=7\num = umap.UMAP(n_components=7 ,random_state=42) \ncomponents_umap = um.fit_transform(Xum) ","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:42:02.995086Z","iopub.execute_input":"2023-08-18T17:42:02.995981Z","iopub.status.idle":"2023-08-18T17:47:12.737666Z","shell.execute_reply.started":"2023-08-18T17:42:02.995939Z","shell.execute_reply":"2023-08-18T17:47:12.736341Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pip install hdbscan","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:47:19.961774Z","iopub.execute_input":"2023-08-18T17:47:19.96215Z","iopub.status.idle":"2023-08-18T17:48:24.552202Z","shell.execute_reply.started":"2023-08-18T17:47:19.96212Z","shell.execute_reply":"2023-08-18T17:48:24.550802Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import hdbscan\nfrom hdbscan import HDBSCAN\nfrom hdbscan.flat import (HDBSCAN_flat,approximate_predict_flat,\n membership_vector_flat, all_points_membership_vectors_flat)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:48:28.370599Z","iopub.execute_input":"2023-08-18T17:48:28.371047Z","iopub.status.idle":"2023-08-18T17:48:28.394552Z","shell.execute_reply.started":"2023-08-18T17:48:28.371008Z","shell.execute_reply":"2023-08-18T17:48:28.393244Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn.model_selection import RandomizedSearchCV\nimport logging\nimport hdbscan\nfrom sklearn.metrics import make_scorer\nlogging.captureWarnings(True)\nclusterer = HDBSCAN_flat(components_umap,prediction_data=True,gen_min_span_tree=True) \n\ngrid = {\"min_samples\" : [10 ,20,30 ,50], \"min_cluster_size\":[20 ,30,50,100 ], \n \"cluster_selection_method\" : [\"eom\",\"leaf\"], \"metric\" : [\"euclidean\",\"manhattan\"] }\n\nvalidity_scorer = make_scorer(hdbscan.validity.validity_index,greater_is_better=True)\nrandom_search = RandomizedSearchCV(clusterer,param_distributions=grid,n_iter=3,scoring=validity_scorer,random_state=22)\nrandom_search.fit(components_umap)\nprint(f\"Best Parameters {random_search.best_params_}\")\nprint(f\"DBCV score :{random_search.best_estimator_.relative_validity_}\")","metadata":{"execution":{"iopub.status.busy":"2023-08-18T15:03:48.886827Z","iopub.execute_input":"2023-08-18T15:03:48.887601Z","iopub.status.idle":"2023-08-18T15:03:48.96181Z","shell.execute_reply.started":"2023-08-18T15:03:48.88756Z","shell.execute_reply":"2023-08-18T15:03:48.960043Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"clusterer = HDBSCAN_flat(components_umap,min_cluster_size=30,min_samples= 30 , metric= \"manhattan\", cluster_selection_method= \"leaf\" ,prediction_data=True)\ncluster_labels = clusterer.labels_\nunique_values, counts = np.unique(cluster_labels, return_counts=True)\n\nprint('\\033[1m \\033[94m Number of HDBSCAN clusters:',len(unique_values))\nprint( ' ------------------------------------------')\nprint(' Unclustered proteins:',counts[0])\n","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:48:55.559734Z","iopub.execute_input":"2023-08-18T17:48:55.560144Z","iopub.status.idle":"2023-08-18T17:50:15.209938Z","shell.execute_reply.started":"2023-08-18T17:48:55.560112Z","shell.execute_reply":"2023-08-18T17:50:15.208313Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nplt.figure(figsize=(6,4))\nplt.bar(unique_values[1:], counts[1:])\nplt.xlabel('HDBSCAN Clusters')\nplt.ylabel('Proteins in Cluster')\nplt.ylim(0, 800)\nplt.title('Protein distriblution in clusters')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:51:14.225916Z","iopub.execute_input":"2023-08-18T17:51:14.226345Z","iopub.status.idle":"2023-08-18T17:51:15.009992Z","shell.execute_reply.started":"2023-08-18T17:51:14.226304Z","shell.execute_reply":"2023-08-18T17:51:15.008692Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Emb_df[\"ClustersH\"] = cluster_labels\nEmb_df['seq']=IdsSeq.Sequences\nEmb_df['EntryID1']=IdsSeq.EntryID1","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:51:20.967361Z","iopub.execute_input":"2023-08-18T17:51:20.967812Z","iopub.status.idle":"2023-08-18T17:51:20.998937Z","shell.execute_reply.started":"2023-08-18T17:51:20.967776Z","shell.execute_reply":"2023-08-18T17:51:20.997696Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Similarity in clusters
\n
\n \n### ** Pairwise Local Sequence Alignment Scores **\nIn bioinformatics, protein sequence alignment is a way of arranging the sequences of DNA, RNA, or protein to identify regions of similarity that may indicate similar function or structure shared between the sequences. Local alignments are more useful for dissimilar sequences that are suspected to contain regions of similarity or similar sequence motifs within their larger sequence context. Pairwise sequence alignment will be examined here for 10 random protein sequences in 10 different clusters in order to have a quick look at the level of similarity in the cluster (computed vs 30 other proteins in cluster). Note that there are also many online similarity searching programs like eg BLAST that could be used in this project for sequence analysis. A similarity score of 30% (% identical over their entire lengths) is an accepted limit indicating that two sequences might be homologous. Results below show that there is high similarity between clustered embeddings (degree of similarity varying from cluster to cluster). We have though to achieve high level of similarity since \"......for very similar proteins (about 50 % identical residues) the chance of completely incorrect annotation is low < 6%\"\" , so if we manage to cluster proteins with high similarity we could transfer Go terms and increase accuracy. [5].","metadata":{}},{"cell_type":"code","source":"from Bio import pairwise2\nimport random\n\ndef cl_similarity(Clusters):\n simdfHDB = pd.DataFrame()\n i=0\n for cl in Clusters:\n similarities=[]\n for prot in range(10):\n sp=[]\n p=random.randint(1, 30)\n sequence1 =cl.seq.iloc[p]\n for second_prot in range(30):\n sequence2 = (cl.seq.iloc[second_prot])\n if sequence1 !=sequence2:\n alignments = pairwise2.align.localms(sequence1, sequence2, 2, -1, -0.5, -0.1)\n best_alignment = alignments[0]\n aligned_sequence1 = best_alignment[0]\n aligned_sequence2 = best_alignment[1]\n total_length = max(len(sequence1), len(sequence2))\n num_matches = sum(a == b for a, b in zip(aligned_sequence1, aligned_sequence2))\n similarity_percentage = (num_matches / total_length) * 100\n sp.append(similarity_percentage)\n similarities.append(sum(sp)/len(sp))\n simdfHDB[f'Cluster_{i}'] = similarities\n i+=1\n return(simdfHDB)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:51:31.215625Z","iopub.execute_input":"2023-08-18T17:51:31.216027Z","iopub.status.idle":"2023-08-18T17:51:31.230979Z","shell.execute_reply.started":"2023-08-18T17:51:31.215991Z","shell.execute_reply":"2023-08-18T17:51:31.229675Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"clustsKM = []\nclustsH = []\nfor value in range(0, 10):\n clustr = Emb_df[Emb_df['ClustersK'] == value].copy()\n clustsKM.append(clustr)\n clustr2 = Emb_df[Emb_df['ClustersH'] == value].copy()\n clustsH.append(clustr2)\n \nprint ('\\033[1m \\033 In-cluster similarity (%) for 10 random proteins (vs 30 other) in 10 random clusters ')\nprint(' ')\nprint('\\033[1m \\033[92m ------------------------------ K-means ---------------------------- ')\ncl_similarity(clustsKM)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:51:34.018533Z","iopub.execute_input":"2023-08-18T17:51:34.018979Z","iopub.status.idle":"2023-08-18T17:52:03.471899Z","shell.execute_reply.started":"2023-08-18T17:51:34.018948Z","shell.execute_reply":"2023-08-18T17:52:03.470538Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('\\033[1m \\033[94m ------------------------------ HDBSCAN ----------------------------')\ncl_similarity(clustsH)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:52:08.59105Z","iopub.execute_input":"2023-08-18T17:52:08.591413Z","iopub.status.idle":"2023-08-18T17:52:17.18463Z","shell.execute_reply.started":"2023-08-18T17:52:08.591383Z","shell.execute_reply":"2023-08-18T17:52:17.183985Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Results show how HDBSCAN outperforms K-means clustering so that will be the choice, despite marking so many proteins as noise.","metadata":{}},{"cell_type":"markdown","source":"### ** Cosine similarity between all pairs of embeddings in clusters 0-5 **\n\nCosine similarity is a measure of similarity between two vectors, two proportional vectors have a cosine similarity of 1 (perfectly similar), two orthogonal vectors have a similarity of 0 (no similarity), and two opposite vectors have a similarity of -1. Lets create cosine similarity matrices for all embeddings in clusters (testing only for clusters 0-5 due to computational limitations) to check how similar embeddings in our formed clusters are. ","metadata":{}},{"cell_type":"code","source":"from sklearn.metrics.pairwise import cosine_similarity\n\ndef plot_sim(ax, sim, title):\n lower_triangular_mask = np.triu(np.ones(sim.shape), k=1)\n sim = np.ma.array(sim, mask=lower_triangular_mask)\n im = ax.imshow(sim, cmap='viridis')\n ax.set_title(title, fontsize=14, weight='bold')\n return im\nsim=[]\nfor i in range(6):\n a=clustsH[i].drop(['ClustersH','seq','EntryID1','ClustersK'], axis=1)\n #Cos_sim.append(a)\n sim.append(cosine_similarity(a))","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:52:50.290453Z","iopub.execute_input":"2023-08-18T17:52:50.290904Z","iopub.status.idle":"2023-08-18T17:52:50.321448Z","shell.execute_reply.started":"2023-08-18T17:52:50.290869Z","shell.execute_reply":"2023-08-18T17:52:50.32058Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"fig, axs = plt.subplots(2, 3, figsize=(15, 10))\nim0=plot_sim(axs[0, 0],sim[0],'Cluster 0')\nim1=plot_sim(axs[0, 1],sim[1],'Cluster 1')\nim2=plot_sim(axs[0, 2],sim[2],'Cluster 2')\nim3=plot_sim(axs[1, 0],sim[3],'Cluster 3')\nim4=plot_sim(axs[1, 1],sim[4],'Cluster 4')\nim5=plot_sim(axs[1, 2],sim[5],'Cluster 5')\nplt.tight_layout(rect=[0, 0.03, 1, 0.95])\nfig.suptitle('Similarity Matrices for embeddings in clusters 0-5', fontsize=18,fontweight='bold')\ncax = plt.axes([0.98, 0.1, 0.02, 0.8]) \ncbar = plt.colorbar(im5, cax=cax, orientation='vertical')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:52:54.229908Z","iopub.execute_input":"2023-08-18T17:52:54.230337Z","iopub.status.idle":"2023-08-18T17:52:55.812286Z","shell.execute_reply.started":"2023-08-18T17:52:54.230303Z","shell.execute_reply":"2023-08-18T17:52:55.810868Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"del sim,im0,im1,im2,im3,im4,im5,Xum\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:52:59.998105Z","iopub.execute_input":"2023-08-18T17:52:59.99858Z","iopub.status.idle":"2023-08-18T17:53:01.512335Z","shell.execute_reply.started":"2023-08-18T17:52:59.998539Z","shell.execute_reply":"2023-08-18T17:53:01.511332Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
GO terms assignment
\n
\n\nHaving clustered our proteins, we could use GO terms that appear in each cluster and try to select the most relevant and similar ones to assign on unseen proteins that belong to this cluster and should show similar function. There are various methods like Resnik or Jiang one could use based on GO distances to the closest common ancestor term to detect relationships between them based on GO hierarchy. They are refered to be problematic though in the DAG structure of GO and got low scoring in some testing I did for this project. I also tried to perform GO Enrichment Analysis to calculate p-values for GO terms that appear in cluster using various significance level (alpha) values but with no success (maybe due to the small size of the clusters or because tehre was no significant association of proteins with particular GO terms). So finally I decided to examine frequency of GO terms in each cluster and in full dataset as well as selection of GOs overrepresented in the cluster compared to the whole dataset (as three thresholds). I did not have time to examine thresholds thoroughly so here is a solution selecting GOs based of frequences and cluster sizes. It is a very basic solution since number of GOs to assign should be selected per protein (maybe examining sequence lengths too) and thresholds used are for now based on best scores, anyway statistical analysis is certainly not enough to guarantee biological relevance. There is an interesting relevant publication for all the above here [6]. Now as for Noise Cluster that contains almost half of the samples, an MLP model will be trained to make predictions for unseen data that belong to noise. ","metadata":{}},{"cell_type":"code","source":"Goterms_in_clusters=[]\ncluster_sizes = counts[1:]\naverage_cluster_size = sum(cluster_sizes) / len(cluster_sizes)\nscaling = 0.2\nfull_protein_set_freq = pd.DataFrame({'GO_full': Go_freqs.index, 'Frequency': Go_freqs.values})\n\nfor i in range(len(unique_values)-1):\n clustHDB=Emb_df[Emb_df[\"ClustersH\"]==i] \n prots=list(clustHDB.EntryID1) \n prot_gos_df = Gos_per_prot[Gos_per_prot['EntryID'].isin(prots)] \n T = [item for sublist in prot_gos_df['term'] for item in sublist] \n ValidGos = [element for element in T if element in targetGos] \n frequency_counts = pd.value_counts(ValidGos).reset_index()\n frequency_counts.columns = ['Go_cluster', 'Frequency_cl']\n merged_df = pd.merge(frequency_counts, full_protein_set_freq, left_on='Go_cluster', right_on='GO_full', how='left')\n merged_df['Freq_full'] = merged_df['Frequency']\n result_df = merged_df.drop(['GO_full', 'Frequency'], axis=1)\n cluster_threshold = scaling * cluster_sizes[i]\n best_mean_Gos = [term for term, freq_cluster, freq_global in zip(result_df['Go_cluster'], result_df['Frequency_cl'], result_df['Freq_full'])\n if freq_cluster >= cluster_threshold]\n Goterms_in_clusters.append(best_mean_Gos)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:53:03.881924Z","iopub.execute_input":"2023-08-18T17:53:03.882334Z","iopub.status.idle":"2023-08-18T17:53:17.976938Z","shell.execute_reply.started":"2023-08-18T17:53:03.8823Z","shell.execute_reply":"2023-08-18T17:53:17.975447Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Multitarget dataframe
\n
\n\nWe shall use existing dataframe with proteins and corresponding Go terms and transform it to a dataframe with binary values. There will be 1500 columns corresponding to selected target Go Terms. Remember that each row corresponds to a protein so values of 1 denote that this Go Term is assigned to the protein. ","metadata":{}},{"cell_type":"code","source":"q=list(train_ids)\nrow_K=Gos_per_prot.EntryID\nsorted_row_K = row_K.sort_values(key=lambda x: x.map({v: i for i, v in enumerate(q)}))\nGos_per_protOK = Gos_per_prot.reindex(sorted_row_K.index).reset_index()\nGos_per_protOK.drop(['index'], axis=1)\n\nY = pd.DataFrame(np.zeros((Gos_per_protOK.shape[0], len(targetGos))), columns=targetGos)\nfor i, row in Gos_per_protOK.iterrows():\n for item in row['term']:\n if item in targetGos:\n Y.at[i, item] = 1 ","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:53:23.512088Z","iopub.execute_input":"2023-08-18T17:53:23.512479Z","iopub.status.idle":"2023-08-18T17:53:36.374226Z","shell.execute_reply.started":"2023-08-18T17:53:23.512447Z","shell.execute_reply":"2023-08-18T17:53:36.372733Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Prediction models
\n
\nAt this point we shall train a simple KNN classifier on our labels from HDBSCAN in order to predict new ones for unseen data later. We could use built-in approximate_predict method after HDBSCAN but it seems that it is not so stable and doesn't give good results. Moreover, we need a MLP to deal with noise cluster (-1) that has almost half of the proteins , to make predictions on unseen data.\n\n### ** KNN classifier**","metadata":{}},{"cell_type":"code","source":"from sklearn.neighbors import KNeighborsClassifier\n\nX_train = Emb_df[column_names_A]\ny_train = Emb_df['ClustersH']\n\nknn_classifier = KNeighborsClassifier(n_neighbors=3)\nknn_classifier.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:53:40.717674Z","iopub.execute_input":"2023-08-18T17:53:40.718109Z","iopub.status.idle":"2023-08-18T17:53:41.065583Z","shell.execute_reply.started":"2023-08-18T17:53:40.718076Z","shell.execute_reply":"2023-08-18T17:53:41.064228Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### ** Multilayer Perceptron Model**","metadata":{}},{"cell_type":"code","source":"unclustered = Emb_df[Emb_df['ClustersH'] == -1]\nY_nonclustered = Y.loc[unclustered.index]\nunclustered.reset_index(drop=True)\nY_nonclustered.reset_index(drop=True)\nXuncl= unclustered.drop([\"ClustersH\",\"seq\",\"EntryID1\",\"ClustersK\"], axis=1) ","metadata":{},"outputs":[],"execution_count":null},{"cell_type":"code","source":"del IdsSeq,X ,sorted_row_K,Gos_per_protOK,Gos_per_prot ,Train_terms,unclustered\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:53:52.545728Z","iopub.execute_input":"2023-08-18T17:53:52.546141Z","iopub.status.idle":"2023-08-18T17:53:53.69861Z","shell.execute_reply.started":"2023-08-18T17:53:52.546106Z","shell.execute_reply":"2023-08-18T17:53:53.697673Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import keras.utils \nfrom keras.utils import np_utils\nfrom keras.optimizers import Adam \nfrom keras.layers import Activation, Dense,Flatten ,Dropout, BatchNormalization, Input,LeakyReLU\nfrom keras.callbacks import EarlyStopping ,ModelCheckpoint,ReduceLROnPlateau\nfrom keras.models import Sequential \nfrom keras import models, utils,backend,layers\nfrom sklearn.model_selection import train_test_split\n\ndef build_model(X):\n keras.backend.clear_session()\n model = models.Sequential()\n n_cols = X.shape[1]\n model.add(layers.Dense(256,activation='relu', input_shape=(n_cols,)))\n model.add(layers.Dropout(0.2))\n model.add(layers.Dense(1500, activation='sigmoid'))\n return model \n\ndef plotmodel(hist):\n train_loss = hist.history['loss']\n val_loss = hist.history['val_loss']\n train_acc = hist.history['accuracy']\n val_acc = hist.history['val_accuracy']\n plt.figure(figsize=(12, 4))\n plt.subplot(1, 2, 1)\n plt.plot(range(1, len(train_loss) + 1), train_loss, label='Training Loss')\n plt.plot(range(1, len(val_loss) + 1), val_loss, label='Validation Loss')\n plt.xlabel('Epochs')\n plt.ylabel('Loss')\n plt.title('Training and Validation Loss')\n plt.legend()\n plt.subplot(1, 2, 2)\n plt.plot(range(1, len(train_acc) + 1), train_acc, label='Training Accuracy')\n plt.plot(range(1, len(val_acc) + 1), val_acc, label='Validation Accuracy')\n plt.xlabel('Epochs')\n plt.ylabel('Accuracy')\n plt.title('Training and Validation Accuracy')\n plt.legend()\n # Show the plot\n plt.tight_layout()\n plt.show()\n\ndef train_model(X,Y):\n X_train, X_test, Y_train, Y_test = train_test_split( X, Y, test_size=0.3, random_state=42)\n reduce_lr = keras.callbacks.ReduceLROnPlateau(monitor='val_logloss', factor=0.3, patience=5, mode='min', min_lr=1E-5)\n early_stopping = keras.callbacks.EarlyStopping(monitor='val_logloss', min_delta=1E-5, patience=15, mode='min',restore_best_weights=True)\n modelkeras=build_model(X)\n modelkeras.compile(optimizer=keras.optimizers.Adam(learning_rate=1e-3, decay=1e-3 / 200), loss=keras.losses.BinaryCrossentropy(), metrics=['accuracy',keras.metrics.AUC(multi_label=True)])\n hist = modelkeras.fit(X_train,Y_train, batch_size=128, epochs=45,verbose=0,validation_data = (X_test,Y_test),callbacks=[reduce_lr, early_stopping])\n loss, accuracy ,AUC = modelkeras.evaluate(X_test,Y_test)\n Y_pred = modelkeras.predict(X_test)\n roc_auc_scores = []\n for i in range(Y.shape[1]):\n unique_classes = np.unique(Y_test.iloc[:, i])\n if len(unique_classes) > 1: # Check if more than one class is present for the label\n roc_auc = roc_auc_score(Y_test.iloc[:, i], Y_pred[:, i])\n roc_auc_scores.append(roc_auc)\n # Step 5: Calculate the average ROC AUC score for the multi-label prediction task\n avg_roc_auc = np.mean(roc_auc_scores)\n \n print(\"ROC AUC Scores for Each Label:\",np.mean(roc_auc_scores))\n plotmodel(hist)\n return modelkeras \n\nmodelkeras_clustered=train_model(Xuncl,Y_nonclustered)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:56:44.689368Z","iopub.execute_input":"2023-08-18T17:56:44.689862Z","iopub.status.idle":"2023-08-18T17:57:07.311271Z","shell.execute_reply.started":"2023-08-18T17:56:44.689822Z","shell.execute_reply":"2023-08-18T17:57:07.310168Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Predict and Submit
\n
","metadata":{}},{"cell_type":"code","source":"X_emb2 = np.load('/kaggle/input/t5embeds/test_embeds.npy')\ntest_ids = np.load('/kaggle/input/t5embeds/test_ids.npy')\nEmb_df2 =pd.DataFrame(X_emb2, columns=column_names_A)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:57:14.256922Z","iopub.execute_input":"2023-08-18T17:57:14.257366Z","iopub.status.idle":"2023-08-18T17:57:26.65973Z","shell.execute_reply.started":"2023-08-18T17:57:14.257328Z","shell.execute_reply":"2023-08-18T17:57:26.657881Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"predicted_clusters_D = knn_classifier.predict(Emb_df2)\nEmb_df2['Clusters'] = predicted_clusters_D","metadata":{"execution":{"iopub.status.busy":"2023-08-18T17:57:35.502728Z","iopub.execute_input":"2023-08-18T17:57:35.503208Z","iopub.status.idle":"2023-08-18T18:06:13.588681Z","shell.execute_reply.started":"2023-08-18T17:57:35.50317Z","shell.execute_reply":"2023-08-18T18:06:13.587316Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Emb_df2['GOlists'] = [Goterms_in_clusters[index] if 0 <= index < len(Goterms_in_clusters) else [] for index in Emb_df2['Clusters']]\nmultiple_prots = [e for e in test_ids for _ in range(len(targetGos))]\nsolution_df = pd.DataFrame({'Proteins': multiple_prots})","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:08:07.431651Z","iopub.execute_input":"2023-08-18T18:08:07.432224Z","iopub.status.idle":"2023-08-18T18:08:07.647057Z","shell.execute_reply.started":"2023-08-18T18:08:07.432187Z","shell.execute_reply":"2023-08-18T18:08:07.645758Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"del multiple_prots,Emb_df,X_emb2, predicted_clusters_D, Y ,Emb_dfstd, X_train\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:08:11.697335Z","iopub.execute_input":"2023-08-18T18:08:11.697675Z","iopub.status.idle":"2023-08-18T18:08:11.746056Z","shell.execute_reply.started":"2023-08-18T18:08:11.69765Z","shell.execute_reply":"2023-08-18T18:08:11.744242Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"Y2 = pd.DataFrame(np.zeros((Emb_df2.shape[0], len(targetGos))), columns=targetGos)\nfor i, row in Emb_df2.iterrows():\n for item in row['GOlists']:\n if item in targetGos:\n Y2.at[i, item] = 1 ","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:12:44.631363Z","iopub.execute_input":"2023-08-18T18:12:44.631771Z","iopub.status.idle":"2023-08-18T18:13:06.176287Z","shell.execute_reply.started":"2023-08-18T18:12:44.631721Z","shell.execute_reply":"2023-08-18T18:13:06.174077Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"solution_df['GO Terms'] = targetGos* (Y2.shape[0])","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:13:14.83764Z","iopub.execute_input":"2023-08-18T18:13:14.83808Z","iopub.status.idle":"2023-08-18T18:13:14.859398Z","shell.execute_reply.started":"2023-08-18T18:13:14.838042Z","shell.execute_reply":"2023-08-18T18:13:14.857302Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"unclustered2 = Emb_df2[Emb_df2['Clusters'] == -1]\nunclustered2= unclustered2.drop(['Clusters','GOlists'], axis=1)\na=unclustered2.index\nY2.loc[a]= modelkeras_clustered.predict(unclustered2.values)","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:15:13.622176Z","iopub.execute_input":"2023-08-18T18:15:13.622614Z","iopub.status.idle":"2023-08-18T18:15:19.22697Z","shell.execute_reply.started":"2023-08-18T18:15:13.622583Z","shell.execute_reply":"2023-08-18T18:15:19.225319Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"del unclustered2 ,Emb_df2\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:15:24.459823Z","iopub.execute_input":"2023-08-18T18:15:24.460234Z","iopub.status.idle":"2023-08-18T18:15:24.504547Z","shell.execute_reply.started":"2023-08-18T18:15:24.460201Z","shell.execute_reply":"2023-08-18T18:15:24.503101Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"solution_df['Predictions'] =Y2.values.ravel()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:15:27.723055Z","iopub.execute_input":"2023-08-18T18:15:27.724003Z","iopub.status.idle":"2023-08-18T18:15:27.730165Z","shell.execute_reply.started":"2023-08-18T18:15:27.723954Z","shell.execute_reply":"2023-08-18T18:15:27.728841Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"del Go_freqs,NumProteins,Y2 , components_umap,y_train,clustHDB ,test_ids,train_ids,seqs,seq_list\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-08-18T18:19:37.192624Z","iopub.execute_input":"2023-08-18T18:19:37.193186Z","iopub.status.idle":"2023-08-18T18:19:38.478227Z","shell.execute_reply.started":"2023-08-18T18:19:37.193149Z","shell.execute_reply":"2023-08-18T18:19:38.47624Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"solution_df.to_csv(\"submission.tsv\",header=False, index=False, sep=\"\\t\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"
\n
Methodology and possible improvements
\n
\n\n#### ** Various techniques and parameters were tested during this work, some of them improved my LB significally.**\n\n**1) Target Go Terms selection**: The Naive method of selecting 1500 Gos based on frequency of appearance was used. Selecting only those GOs that appeared on all three aspects did not show better LB scores. I also tried to select GO terms based on significance (.obo graph degree centrality or IA scores) but it did not work out. There are other methods for the choice of target GO terms that could lead to more important and relevant GO terms selection,influencing final results. \n**2) Embeddings selection**: ProtT5 is presented here but using embeddings from other models such as SeqVec, or combining Esm with T5 could improve LB. \n**3) Clustering method**: When using PCA/K-means for clustering ,pairwise local sequence alignment average scores were always low for small and large number of clusters (<35%). UMAP/HDBSCAN clustering was always more successful showing higher in-cluster similarity and so it why was prefered. Fine tuning UMAP/HDBSCAN could lead to better clustering reducing proteins characterized as noise leading to higher scores \n**4) Go terms assignment** : In-cluster GO terms selection was based on frequency of appearance both in cluster and in full set. This is a very basic appoach but working ...., anyway more sophisticated techniques could be used for Go terms assignment. \n**5) Classifiers** : Due to contest resources limitations ,I only tried MLP and KNN classifiers. There is a plethora of well known algorithms that could handle this multilabel task differently and could be used. Also imbalanced target label distributions and dependencies between labels should be examined. Dealing with this aspect should also improve model performance. ","metadata":{}},{"cell_type":"markdown","source":"
\n
References
\n
\nThere are really many inspiring notebooks and discussions published for this contest covering various aspects, that helped me a lot in to build this kernel either using code or ideas (and some help from ChatGPT as well!!) Just to mentiom a few :\n\nhttps://www.kaggle.com/code/alexandervc/baseline-multilabel-to-multitarget-binary \nhttps://https://www.kaggle.com/code/shtrausslearning/biopython-bioinformatics-basics#5-|-PAIRWISE-SEQUENCE-ALIGNMENT \nhttps://www.kaggle.com/code/gusthema/cafa-5-protein-function-with-tensorflow \nhttps://www.kaggle.com/code/leonidkulyk/eda-cafa5-pfp-interactive-dags-plotly \nhttps://www.youtube.com/watch?v=h4zTjPzPud4 Andrey Shevtsov \"Introduction to CAFA5 Protein function prediction Kaggle competition\" \n\n### **Relevant Publications** \n\n[1] Littmann, M., Heinzinger, M., Dallago, C. et al. Embeddings from deep learning transfer GO annotations beyond homology. Sci Rep 11, 1160 (2021). https://doi.org/10.1038/s41598-020-80786-0 \n[2] Vesztrocy AW, Dessimoz C. A Gene Ontology Tutorial in Python. Methods Mol Biol. 2017;1446:221-229. doi: 10.1007/978-1-4939-3743-1_16. PMID: 27812946. \n[3] Tran C, Khadkikar S, Porollo A. Survey of Protein Sequence Embedding Models. Int J Mol Sci. 2023 Feb 14;24(4):3775. doi: 10.3390/ijms24043775. PMID: 36835188; PMCID: PMC9963412. \n[4] Allaoui M, Kherfi ML, Cheriet A. Considerably Improving Clustering Algorithms Using UMAP Dimensionality Reduction Technique: A Comparative Study. Image and Signal Processing. 2020 Jun 5;12119:317–25. doi: 10.1007/978-3-030-51935-3_34. PMCID: PMC7340901. \n[5] Sangar, V., Blankenberg, D.J., Altman, N. et al. Quantitative sequence-function relationships in proteins based on gene ontology. BMC Bioinformatics 8, 294 (2007). https://doi.org/10.1186/1471-2105-8-294 \n[6] Warith Eddine Djeddi, Sadok Ben Yahia, Engelbert Mephu Nguifo. A novel approach for predicting\nprotein functions by transferring annotation via alignment networks. 2019. hal-02070419\n\nThank you very much for your time reading this kernel. \n#### ** If you found something that you liked or gave you an idea........ don't forget to Upvote! **\n\n#### ** Good luck to everyone that dedicated time to this competition! **","metadata":{}}]}
![](\"https://upload.wikimedia.org/wikipedia/commons/6/60/Myoglobin.png\")
\n