{"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":"# Interesting finds:\n\n- train and test share ~70,000 protein sequences \n- train and test each have about ~6000 duplicated sequences (different id but identical sequences)\n- in train, ~2300 of the distinct duplicated sequences have different terms (different id, identical sequence, different labels)\n- the most common 6000 terms are used as labels for 90% of the training sequences\n- 99% of the sequences in both train and test are less than 3000 elements long\n- 1/3 of the terms have other relationships besides \"is_a\" (\"part_of\" and \"regulates\")","metadata":{}},{"cell_type":"markdown","source":"# Dependencies","metadata":{}},{"cell_type":"code","source":"from tqdm import tqdm\nimport plotly.express as px\nimport numpy as np\nimport math","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:32:20.257368Z","iopub.execute_input":"2023-06-25T19:32:20.258202Z","iopub.status.idle":"2023-06-25T19:32:21.173024Z","shell.execute_reply.started":"2023-06-25T19:32:20.258167Z","shell.execute_reply":"2023-06-25T19:32:21.171791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Get the data","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nfrom Bio import SeqIO\n\ntrain_sequences = SeqIO.parse('/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta', \"fasta\")\ntrain_sequence_df = pd.DataFrame([{\"id\": s.id, \"seq\": str(s.seq)} for s in train_sequences]).set_index('id')\n\ntrain_terms_df = pd.read_csv('/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv', sep=\"\\t\")\ntrain_terms_df = train_terms_df.rename(columns={\"EntryID\": \"id\"})\ntrain_terms_df_grouped = train_terms_df.groupby(\"id\").aggregate(lambda term: sorted(term.unique().tolist()))[['term']]\n\ntrain_df = train_sequence_df.join(train_terms_df_grouped)\ntrain_df","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:32:21.175562Z","iopub.execute_input":"2023-06-25T19:32:21.176012Z","iopub.status.idle":"2023-06-25T19:32:52.666331Z","shell.execute_reply.started":"2023-06-25T19:32:21.175972Z","shell.execute_reply":"2023-06-25T19:32:52.665087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_sequences = SeqIO.parse('/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset.fasta', \"fasta\")\ntest_df = pd.DataFrame([{\"id\": s.id, \"seq\": str(s.seq)} for s in test_sequences]).set_index('id')\ntest_df","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:32:52.667687Z","iopub.execute_input":"2023-06-25T19:32:52.66811Z","iopub.status.idle":"2023-06-25T19:32:55.475153Z","shell.execute_reply.started":"2023-06-25T19:32:52.66808Z","shell.execute_reply":"2023-06-25T19:32:55.473871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install obonet -q\n\nimport obonet\nterm_graph = obonet.read_obo('/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo')\nterm_df = pd.DataFrame(term_graph.nodes.values())\nterm_df['id'] = term_graph.nodes.keys()\nterm_df = term_df.set_index('id')\nterm_df","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:32:55.476897Z","iopub.execute_input":"2023-06-25T19:32:55.477307Z","iopub.status.idle":"2023-06-25T19:33:31.055037Z","shell.execute_reply.started":"2023-06-25T19:32:55.477269Z","shell.execute_reply":"2023-06-25T19:33:31.05385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Counts","metadata":{}},{"cell_type":"markdown","source":"### number of amino acids in train vs test","metadata":{}},{"cell_type":"code","source":"test_df['seq_len'] = test_df.seq.map(lambda seq: len(seq))\ntrain_df['seq_len'] = train_df.seq.map(lambda seq: len(seq))\n\nprint('total # of amino acids in train:', train_df.seq_len.sum())\nprint('total # of amino acids in test: ', test_df.seq_len.sum())","metadata":{"execution":{"iopub.status.busy":"2023-06-25T19:34:00.276809Z","iopub.execute_input":"2023-06-25T19:34:00.277193Z","iopub.status.idle":"2023-06-25T19:34:00.531246Z","shell.execute_reply.started":"2023-06-25T19:34:00.277164Z","shell.execute_reply":"2023-06-25T19:34:00.52993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### count proteins in train vs test, and overlap","metadata":{}},{"cell_type":"code","source":"print(\"# of proteins in train:       \", len(train_df.index))\nprint(\"# of proteins in test:        \", len(test_df.index))\nprint(\"# of proteins in intersection: \", len(set(train_df.index) & set(test_df.index)))","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:19:06.413991Z","iopub.execute_input":"2023-06-23T22:19:06.415059Z","iopub.status.idle":"2023-06-23T22:19:06.546238Z","shell.execute_reply.started":"2023-06-23T22:19:06.415021Z","shell.execute_reply":"2023-06-23T22:19:06.544965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### count terms in train vs ontology","metadata":{}},{"cell_type":"code","source":"print(\"total # of terms in train:    \", len(train_terms_df))\nprint(\"# of distinct terms in train:   \", len(set(train_terms_df.term)))\nprint(\"# of terms in ontology:         \", len(term_graph))","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:19:06.547858Z","iopub.execute_input":"2023-06-23T22:19:06.54907Z","iopub.status.idle":"2023-06-23T22:19:07.488914Z","shell.execute_reply.started":"2023-06-23T22:19:06.549035Z","shell.execute_reply":"2023-06-23T22:19:07.487706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### count duplicate sequences","metadata":{}},{"cell_type":"code","source":"print(\"# of dupes in train: \", len(train_df.index) - len(set(train_df.seq)))\nprint(\"# of dupes in test:  \", len(test_df.index) - len(set(test_df.seq)))","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:19:07.490479Z","iopub.execute_input":"2023-06-23T22:19:07.491474Z","iopub.status.idle":"2023-06-23T22:19:07.778125Z","shell.execute_reply.started":"2023-06-23T22:19:07.491438Z","shell.execute_reply":"2023-06-23T22:19:07.776875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### analyze duplicates in training\n\ndo the duplicate sequences all have the same terms?","metadata":{}},{"cell_type":"code","source":"train_df['seq_hash'] = train_df.seq.map(lambda seq: hash(seq))\ntrain_df['term_hash'] = train_df.term.map(lambda li: hash(','.join(sorted(li))))\n\ntrain_dupe_seqs = train_df.seq_hash.value_counts()\ntrain_dupe_seqs = list(train_dupe_seqs[train_dupe_seqs > 1].keys())\ntrain_dupe_seqs_df = train_df[train_df.seq_hash.isin(train_dupe_seqs)].sort_values('seq')\ntrain_dupe_seqs_df.to_csv('duplicate_sequences.csv')\n\nprint('total # of duplicated sequences:', len(train_dupe_seqs_df))\nprint('distinct duplicated sequences:  ', len(train_dupe_seqs))\n\ninconsistent_terms_per_seq = train_df.groupby('seq_hash').term_hash.nunique()\ninconsistent_term_counts = inconsistent_terms_per_seq[inconsistent_terms_per_seq > 1]\nseqs_inconsistent_terms_df = train_dupe_seqs_df[train_dupe_seqs_df.seq_hash.isin(inconsistent_term_counts.index)]\n\nprint('total # of duped seqs which have inconsistent terms:    ', len(seqs_inconsistent_terms_df))\nprint('distinct # of duped seqs which have inconsistent terms: ', len(inconsistent_term_counts))\nprint('top 10 duped seqs with the most # of inconsistent terms:', sorted(inconsistent_term_counts, reverse=True)[:10])\n\nseqs_inconsistent_terms_df","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:22:12.66635Z","iopub.execute_input":"2023-06-23T22:22:12.6668Z","iopub.status.idle":"2023-06-23T22:22:14.041369Z","shell.execute_reply.started":"2023-06-23T22:22:12.666768Z","shell.execute_reply":"2023-06-23T22:22:14.040189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Distributions in training data","metadata":{}},{"cell_type":"markdown","source":"### terms per protein","metadata":{}},{"cell_type":"code","source":"train_df['num_terms'] = train_df.term.map(lambda term: len(term))\n\nprint('# terms per protein sequence:')\ntrain_df.num_terms.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:22:16.663407Z","iopub.execute_input":"2023-06-23T22:22:16.663904Z","iopub.status.idle":"2023-06-23T22:22:16.818168Z","shell.execute_reply.started":"2023-06-23T22:22:16.663861Z","shell.execute_reply":"2023-06-23T22:22:16.816956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### proteins per term statistics + percentage plot","metadata":{}},{"cell_type":"code","source":"terms = train_terms_df.groupby('term').id.nunique().to_frame().sort_values('id')\ntotal_seqs = terms.id.sum()\nterms['cumulative_count'] = terms['id'].cumsum()\nterms['cumulative_pct'] = list(100 - terms.cumulative_count.map(lambda c: c / total_seqs * 100).iloc[::-1])\npx.line(terms, x=range(len(terms)), y='cumulative_pct', title=\"cumulative % of proteins using terms\")","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:22:17.106267Z","iopub.execute_input":"2023-06-23T22:22:17.106817Z","iopub.status.idle":"2023-06-23T22:22:22.22968Z","shell.execute_reply.started":"2023-06-23T22:22:17.106764Z","shell.execute_reply":"2023-06-23T22:22:22.228414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### protein sequence length (cumulative percentage plot) in train vs test","metadata":{}},{"cell_type":"code","source":"print('train sequence lengths:')\ntrain_df.seq_len.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:22:22.232083Z","iopub.execute_input":"2023-06-23T22:22:22.232569Z","iopub.status.idle":"2023-06-23T22:22:22.378282Z","shell.execute_reply.started":"2023-06-23T22:22:22.232531Z","shell.execute_reply":"2023-06-23T22:22:22.377143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('test sequence lengths:')\ntest_df.seq_len.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:22:22.379782Z","iopub.execute_input":"2023-06-23T22:22:22.380512Z","iopub.status.idle":"2023-06-23T22:22:22.522731Z","shell.execute_reply.started":"2023-06-23T22:22:22.380477Z","shell.execute_reply":"2023-06-23T22:22:22.521627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_seq_lens = train_df.seq_len.sort_values().to_frame()\ntrain_idxs = np.linspace(0, len(train_seq_lens)-1, 10_000).astype(np.int64)\n\ntest_seq_lens = test_df.seq_len.sort_values().to_frame()\ntest_idxs = np.linspace(0, len(test_seq_lens)-1, 10_000).astype(np.int64)\n\nfig = px.line()\nfig.add_scatter(x=np.arange(10_001) / 100, y=train_seq_lens.seq_len[train_idxs], name=\"train\")\nfig.add_scatter(x=np.arange(10_001) / 100, y=test_seq_lens.seq_len[test_idxs], name=\"test\")\nfig.update_layout(\n    xaxis_title=\"% of sequences shorter than the given length\",\n    yaxis_title=\"sequence length\"\n)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:37:59.375828Z","iopub.execute_input":"2023-06-23T22:37:59.376344Z","iopub.status.idle":"2023-06-23T22:37:59.588241Z","shell.execute_reply.started":"2023-06-23T22:37:59.376287Z","shell.execute_reply":"2023-06-23T22:37:59.58667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Distributions in ontology","metadata":{}},{"cell_type":"markdown","source":"### relational ties (is_a, part_of, …)","metadata":{}},{"cell_type":"code","source":"term_df['relationship_len'] = term_df.relationship.map(lambda r: float('nan') if type(r) is float and math.isnan(r) else len(r))\nterm_df.relationship_len.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:35:33.603097Z","iopub.execute_input":"2023-06-23T22:35:33.603575Z","iopub.status.idle":"2023-06-23T22:35:33.663946Z","shell.execute_reply.started":"2023-06-23T22:35:33.603539Z","shell.execute_reply":"2023-06-23T22:35:33.662504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"term_df['relationship_types'] = term_df.relationship.map(\n    lambda r: \n        float('nan') \n        if type(r) is float and math.isnan(r) else \n        ', '.join(set([e.split(' ')[0] for e in r]))\n)\nterm_df.relationship_types.value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:36:41.042414Z","iopub.execute_input":"2023-06-23T22:36:41.042957Z","iopub.status.idle":"2023-06-23T22:36:41.132163Z","shell.execute_reply.started":"2023-06-23T22:36:41.042919Z","shell.execute_reply":"2023-06-23T22:36:41.130755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### # of \"is_a\" relationships per term","metadata":{}},{"cell_type":"code","source":"term_df['is_a_len'] = term_df.is_a.map(\n    lambda is_a: \n        float('nan') \n        if type(is_a) is float and math.isnan(is_a) else \n        len(is_a)\n)\nterm_df.is_a_len.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:54:46.939002Z","iopub.execute_input":"2023-06-23T22:54:46.939424Z","iopub.status.idle":"2023-06-23T22:54:46.996033Z","shell.execute_reply.started":"2023-06-23T22:54:46.939391Z","shell.execute_reply":"2023-06-23T22:54:46.994933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### distinct by subontology category, vs by # of proteins in train","metadata":{}},{"cell_type":"code","source":"term_df.namespace.value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:56:05.934706Z","iopub.execute_input":"2023-06-23T22:56:05.935223Z","iopub.status.idle":"2023-06-23T22:56:05.955731Z","shell.execute_reply.started":"2023-06-23T22:56:05.935177Z","shell.execute_reply":"2023-06-23T22:56:05.954495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_terms_df.groupby('term').aspect.agg(lambda x: x.iloc[0]).value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T23:02:48.058565Z","iopub.execute_input":"2023-06-23T23:02:48.058953Z","iopub.status.idle":"2023-06-23T23:02:51.049099Z","shell.execute_reply.started":"2023-06-23T23:02:48.058924Z","shell.execute_reply":"2023-06-23T23:02:51.048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_terms_df.aspect.value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:58:38.104067Z","iopub.execute_input":"2023-06-23T22:58:38.104726Z","iopub.status.idle":"2023-06-23T22:58:38.972726Z","shell.execute_reply.started":"2023-06-23T22:58:38.104692Z","shell.execute_reply":"2023-06-23T22:58:38.971629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### todo: ontology chain length","metadata":{}},{"cell_type":"markdown","source":"# Save all the results","metadata":{}},{"cell_type":"code","source":"train_df.to_csv('train_sequences.csv')\nterm_df.to_csv('terms.csv')","metadata":{"execution":{"iopub.status.busy":"2023-06-23T22:19:09.595668Z","iopub.status.idle":"2023-06-23T22:19:09.596206Z","shell.execute_reply.started":"2023-06-23T22:19:09.595905Z","shell.execute_reply":"2023-06-23T22:19:09.595924Z"},"trusted":true},"execution_count":null,"outputs":[]}]}