{"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":"code","source":"import numpy as np \nimport pandas as pd \nfrom IPython.display import clear_output\nimport re\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-05-15T16:06:35.344053Z","iopub.execute_input":"2023-05-15T16:06:35.344403Z","iopub.status.idle":"2023-05-15T16:06:35.383986Z","shell.execute_reply.started":"2023-05-15T16:06:35.344372Z","shell.execute_reply":"2023-05-15T16:06:35.383026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Download diamond","metadata":{}},{"cell_type":"code","source":"!wget http://github.com/bbuchfink/diamond/releases/download/v2.1.6/diamond-linux64.tar.gz\n!tar xzf diamond-linux64.tar.gz\n!rm diamond-linux64.tar.gz","metadata":{"execution":{"iopub.status.busy":"2023-05-15T16:10:23.714454Z","iopub.execute_input":"2023-05-15T16:10:23.714822Z","iopub.status.idle":"2023-05-15T16:10:27.986066Z","shell.execute_reply.started":"2023-05-15T16:10:23.714792Z","shell.execute_reply":"2023-05-15T16:10:27.984705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Run diamond","metadata":{}},{"cell_type":"code","source":"#Run diamond\n!./diamond blastp -d /kaggle/input/deepgoplus-data/data/train_data.dmnd --more-sensitive -q /kaggle/input/cafa-5-protein-function-prediction/Test\\ \\(Targets\\)/testsuperset.fasta \\\n--outfmt 6 qseqid sseqid bitscore > diamond.res","metadata":{"execution":{"iopub.status.busy":"2023-05-15T16:10:27.989248Z","iopub.execute_input":"2023-05-15T16:10:27.990095Z","iopub.status.idle":"2023-05-15T16:27:12.80779Z","shell.execute_reply.started":"2023-05-15T16:10:27.990059Z","shell.execute_reply":"2023-05-15T16:27:12.806491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Run deepgoplus","metadata":{}},{"cell_type":"markdown","source":"We use code from https://github.com/bio-ontology-research-group/deepgoplus. ","metadata":{}},{"cell_type":"markdown","source":"```\n! git clone https://github.com/bio-ontology-research-group/deepgoplus.git\n! gzip diamond.res\n! python deepgoplus/predict.py -if /kaggle/input/cafa-5-protein-function-prediction/Test\\ \\(Targets\\)/testsuperset.fasta \\\n-df /kaggle/working/diamond.res.gz\n```","metadata":{"execution":{"iopub.status.busy":"2023-05-11T13:05:52.810388Z","iopub.execute_input":"2023-05-11T13:05:52.810854Z","iopub.status.idle":"2023-05-11T13:09:58.507197Z","shell.execute_reply.started":"2023-05-11T13:05:52.810811Z","shell.execute_reply":"2023-05-11T13:09:58.505844Z"}}},{"cell_type":"markdown","source":"Instead of cloning the repo we copy all necessary code into the notebook and run it.","metadata":{}},{"cell_type":"code","source":"# from https://github.com/bio-ontology-research-group/deepgoplus/blob/master/utils.py\nfrom collections import deque, Counter\nimport warnings\nimport pandas as pd\nimport numpy as np\nfrom xml.etree import ElementTree as ET\nimport math\n\nBIOLOGICAL_PROCESS = 'GO:0008150'\nMOLECULAR_FUNCTION = 'GO:0003674'\nCELLULAR_COMPONENT = 'GO:0005575'\nFUNC_DICT = {\n    'cc': CELLULAR_COMPONENT,\n    'mf': MOLECULAR_FUNCTION,\n    'bp': BIOLOGICAL_PROCESS}\n\nNAMESPACES = {\n    'cc': 'cellular_component',\n    'mf': 'molecular_function',\n    'bp': 'biological_process'\n}\n\nEXP_CODES = set([\n    'EXP', 'IDA', 'IPI', 'IMP', 'IGI', 'IEP', 'TAS', 'IC',\n    'HTP', 'HDA', 'HMP', 'HGI', 'HEP'])\n\n# CAFA4 Targets\nCAFA_TARGETS = set([\n    '287', '3702', '4577', '6239', '7227', '7955', '9606', '9823', '10090',\n    '10116', '44689', '83333', '99287', '226900', '243273', '284812', '559292'])\n\ndef is_cafa_target(org):\n    return org in CAFA_TARGETS\n\ndef is_exp_code(code):\n    return code in EXP_CODES\n\n\nclass Ontology(object):\n\n    def __init__(self, filename='data/go.obo', with_rels=False):\n        self.ont = self.load(filename, with_rels)\n        self.ic = None\n\n    def has_term(self, term_id):\n        return term_id in self.ont\n\n    def get_term(self, term_id):\n        if self.has_term(term_id):\n            return self.ont[term_id]\n        return None\n\n    def calculate_ic(self, annots):\n        cnt = Counter()\n        for x in annots:\n            cnt.update(x)\n        self.ic = {}\n        for go_id, n in cnt.items():\n            parents = self.get_parents(go_id)\n            if len(parents) == 0:\n                min_n = n\n            else:\n                min_n = min([cnt[x] for x in parents])\n\n            self.ic[go_id] = math.log(min_n / n, 2)\n    \n    def get_ic(self, go_id):\n        if self.ic is None:\n            raise Exception('Not yet calculated')\n        if go_id not in self.ic:\n            return 0.0\n        return self.ic[go_id]\n\n    def load(self, filename, with_rels):\n        ont = dict()\n        obj = None\n        with open(filename, 'r') as f:\n            for line in f:\n                line = line.strip()\n                if not line:\n                    continue\n                if line == '[Term]':\n                    if obj is not None:\n                        ont[obj['id']] = obj\n                    obj = dict()\n                    obj['is_a'] = list()\n                    obj['part_of'] = list()\n                    obj['regulates'] = list()\n                    obj['alt_ids'] = list()\n                    obj['is_obsolete'] = False\n                    continue\n                elif line == '[Typedef]':\n                    if obj is not None:\n                        ont[obj['id']] = obj\n                    obj = None\n                else:\n                    if obj is None:\n                        continue\n                    l = line.split(\": \")\n                    if l[0] == 'id':\n                        obj['id'] = l[1]\n                    elif l[0] == 'alt_id':\n                        obj['alt_ids'].append(l[1])\n                    elif l[0] == 'namespace':\n                        obj['namespace'] = l[1]\n                    elif l[0] == 'is_a':\n                        obj['is_a'].append(l[1].split(' ! ')[0])\n                    elif with_rels and l[0] == 'relationship':\n                        it = l[1].split()\n                        # add all types of relationships\n                        obj['is_a'].append(it[1])\n                    elif l[0] == 'name':\n                        obj['name'] = l[1]\n                    elif l[0] == 'is_obsolete' and l[1] == 'true':\n                        obj['is_obsolete'] = True\n            if obj is not None:\n                ont[obj['id']] = obj\n        for term_id in list(ont.keys()):\n            for t_id in ont[term_id]['alt_ids']:\n                ont[t_id] = ont[term_id]\n            if ont[term_id]['is_obsolete']:\n                del ont[term_id]\n        for term_id, val in ont.items():\n            if 'children' not in val:\n                val['children'] = set()\n            for p_id in val['is_a']:\n                if p_id in ont:\n                    if 'children' not in ont[p_id]:\n                        ont[p_id]['children'] = set()\n                    ont[p_id]['children'].add(term_id)\n        return ont\n\n\n    def get_anchestors(self, term_id):\n        if term_id not in self.ont:\n            return set()\n        term_set = set()\n        q = deque()\n        q.append(term_id)\n        while(len(q) > 0):\n            t_id = q.popleft()\n            if t_id not in term_set:\n                term_set.add(t_id)\n                for parent_id in self.ont[t_id]['is_a']:\n                    if parent_id in self.ont:\n                        q.append(parent_id)\n        return term_set\n\n\n    def get_parents(self, term_id):\n        if term_id not in self.ont:\n            return set()\n        term_set = set()\n        for parent_id in self.ont[term_id]['is_a']:\n            if parent_id in self.ont:\n                term_set.add(parent_id)\n        return term_set\n\n\n    def get_namespace_terms(self, namespace):\n        terms = set()\n        for go_id, obj in self.ont.items():\n            if obj['namespace'] == namespace:\n                terms.add(go_id)\n        return terms\n\n    def get_namespace(self, term_id):\n        return self.ont[term_id]['namespace']\n    \n    def get_term_set(self, term_id):\n        if term_id not in self.ont:\n            return set()\n        term_set = set()\n        q = deque()\n        q.append(term_id)\n        while len(q) > 0:\n            t_id = q.popleft()\n            if t_id not in term_set:\n                term_set.add(t_id)\n                for ch_id in self.ont[t_id]['children']:\n                    q.append(ch_id)\n        return term_set\n\ndef read_fasta(filename):\n    seqs = list()\n    info = list()\n    seq = ''\n    inf = ''\n    with open(filename, 'r') as f:\n        for line in f:\n            line = line.strip()\n            if line.startswith('>'):\n                if seq != '':\n                    seqs.append(seq)\n                    info.append(inf)\n                    seq = ''\n                inf = line[1:]\n            else:\n                seq += line\n        seqs.append(seq)\n        info.append(inf)\n    return info, seqs\n\n\nclass DataGenerator(object):\n\n    def __init__(self, batch_size, is_sparse=False):\n        self.batch_size = batch_size\n        self.is_sparse = is_sparse\n\n    def fit(self, inputs, targets=None):\n        self.start = 0\n        self.inputs = inputs\n        self.targets = targets\n        if isinstance(self.inputs, tuple) or isinstance(self.inputs, list):\n            self.size = self.inputs[0].shape[0]\n        else:\n            self.size = self.inputs.shape[0]\n        self.has_targets = targets is not None\n\n    def __next__(self):\n        return self.next()\n\n    def reset(self):\n        self.start = 0\n\n    def next(self):\n        if self.start < self.size:\n            batch_index = np.arange(\n                self.start, min(self.size, self.start + self.batch_size))\n            if isinstance(self.inputs, tuple) or isinstance(self.inputs, list):\n                res_inputs = []\n                for inp in self.inputs:\n                    if self.is_sparse:\n                        res_inputs.append(\n                            inp[batch_index, :].toarray())\n                    else:\n                        res_inputs.append(inp[batch_index, :])\n            else:\n                if self.is_sparse:\n                    res_inputs = self.inputs[batch_index, :].toarray()\n                else:\n                    res_inputs = self.inputs[batch_index, :]\n            self.start += self.batch_size\n            if self.has_targets:\n                if self.is_sparse:\n                    labels = self.targets[batch_index, :].toarray()\n                else:\n                    labels = self.targets[batch_index, :]\n                return (res_inputs, labels)\n            return res_inputs\n        else:\n            self.reset()\n            return self.next()\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-15T16:27:12.810673Z","iopub.execute_input":"2023-05-15T16:27:12.811637Z","iopub.status.idle":"2023-05-15T16:27:12.854292Z","shell.execute_reply.started":"2023-05-15T16:27:12.811593Z","shell.execute_reply":"2023-05-15T16:27:12.853219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://github.com/bio-ontology-research-group/deepgoplus/blob/master/predict.py\n\nimport numpy as np\nimport pandas as pd\nfrom tensorflow.keras.models import load_model\nfrom subprocess import Popen, PIPE\nimport time\n\nAALETTER = [\n    'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',\n    'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V']\nAANUM = len(AALETTER)\nAAINDEX = dict()\nfor i in range(len(AALETTER)):\n    AAINDEX[AALETTER[i]] = i + 1\nMAXLEN = 2000\n\n\ndef to_onehot(seq, start=0):\n    onehot = np.zeros((MAXLEN, 21), dtype=np.int32)\n    l = min(MAXLEN, len(seq))\n    for i in range(start, start + l):\n        onehot[i, AAINDEX.get(seq[i - start], 0)] = 1\n    onehot[0:start, 0] = 1\n    onehot[start + l:, 0] = 1\n    return onehot\n\n\ndef read_fasta(filename, chunk_size):\n    seqs = list()\n    info = list()\n    seq = ''\n    inf = ''\n    with open(filename, 'rt') as f:\n        for line in f:\n            line = line.strip()\n            if line.startswith('>'):\n                if seq != '':\n                    seqs.append(seq)\n                    info.append(inf)\n                    if len(info) == chunk_size:\n                        yield (info, seqs)\n                        seqs = list()\n                        info = list()\n                    seq = ''\n                inf = line[1:].split()[0]\n            else:\n                seq += line\n        seqs.append(seq)\n        info.append(inf)\n    yield (info, seqs)\n\ndef get_data(sequences):\n    pred_seqs = []\n    ids = []\n    for i, seq in enumerate(sequences):\n        if len(seq) > MAXLEN:\n            st = 0\n            while st < len(seq):\n                pred_seqs.append(seq[st: st + MAXLEN])\n                ids.append(i)\n                st += MAXLEN - 128\n        else:\n            pred_seqs.append(seq)\n            ids.append(i)\n    n = len(pred_seqs)\n    data = np.zeros((n, MAXLEN, 21), dtype=np.float32)\n    \n    for i in range(n):\n        seq = pred_seqs[i]\n        data[i, :, :] = to_onehot(seq)\n    return ids, data","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-15T16:27:12.855865Z","iopub.execute_input":"2023-05-15T16:27:12.856206Z","iopub.status.idle":"2023-05-15T16:27:19.90092Z","shell.execute_reply.started":"2023-05-15T16:27:12.856173Z","shell.execute_reply":"2023-05-15T16:27:19.899981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://github.com/bio-ontology-research-group/deepgoplus/blob/master/predict.py\n\n\nin_file='/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset.fasta'\nout_file='submission.tsv'\ngo_file='/kaggle/input/deepgoplus-data/data/go.obo'\nmodel_file='/kaggle/input/deepgoplus-data/data/model.h5'\nterms_file='/kaggle/input/deepgoplus-data/data/terms.pkl'\nannotations_file='/kaggle/input/deepgoplus-data/data/train_data.pkl'\nchunk_size=1000\ndiamond_file='/kaggle/working/diamond.res'\nthreshold=0.1\nbatch_size=32\nalpha=0.5\n\n# Load GO and read list of all terms\ngo = Ontology(go_file, with_rels=True)\nterms_df = pd.read_pickle(terms_file)\nterms = terms_df['terms'].values.flatten()\n\n# Read known experimental annotations\nannotations = {}\ndf = pd.read_pickle(annotations_file)\nfor row in df.itertuples():\n    annotations[row.proteins] = set(row.prop_annotations)\n\ngo.calculate_ic(annotations.values())\n\ndiamond_preds = {}\nmapping = {}\nwith open(diamond_file, 'rt') as f:\n    for line in f:\n        it = line.strip().split()\n        if it[0] not in mapping:\n            mapping[it[0]] = {}\n        mapping[it[0]][it[1]] = float(it[2])\nfor prot_id, sim_prots in mapping.items():\n    annots = {}\n    allgos = set()\n    total_score = 0.0\n    for p_id, score in sim_prots.items():\n        allgos |= annotations[p_id]\n        total_score += score\n    allgos = list(sorted(allgos))\n    sim = np.zeros(len(allgos), dtype=np.float32)\n    for j, go_id in enumerate(allgos):\n        s = 0.0\n        for p_id, score in sim_prots.items():\n            if go_id in annotations[p_id]:\n                s += score\n        sim[j] = s / total_score\n    for go_id, score in zip(allgos, sim):\n        annots[go_id] = score\n    diamond_preds[prot_id] = annots\n    \n# Load CNN model\nmodel = load_model(model_file)\n# Alphas for the latest model\nalphas = {NAMESPACES['mf']: 0.55, NAMESPACES['bp']: 0.59, NAMESPACES['cc']: 0.46}\n# Alphas for the cafa2 model\n# alphas = {NAMESPACES['mf']: 0.63, NAMESPACES['bp']: 0.68, NAMESPACES['cc']: 0.48}\n    \nstart_time = time.time()\ntotal_seq = 0\nw = open(out_file, 'wt')\nfor prot_ids, sequences in read_fasta(in_file, chunk_size):\n    total_seq += len(prot_ids)\n    deep_preds = {}\n    ids, data = get_data(sequences)\n\n    preds = model.predict(data, batch_size=batch_size)\n    assert preds.shape[1] == len(terms)\n    for i, j in enumerate(ids):\n        prot_id = prot_ids[j]\n        if prot_id not in deep_preds:\n            deep_preds[prot_id] = {}\n        for l in range(len(terms)):\n            if preds[i, l] >= 0.01: # Filter out very low scores\n                if terms[l] not in deep_preds[prot_id]:\n                    deep_preds[prot_id][terms[l]] = preds[i, l]\n                else:\n                    deep_preds[prot_id][terms[l]] = max(\n                        deep_preds[prot_id][terms[l]], preds[i, l])\n    # Combine diamond preds and deepgo\n    for prot_id in prot_ids:\n        annots = {}\n        if prot_id in diamond_preds:\n            for go_id, score in diamond_preds[prot_id].items():\n                annots[go_id] = score * alphas[go.get_namespace(go_id)]\n        for go_id, score in deep_preds[prot_id].items():\n            if go_id in annots:\n                annots[go_id] += (1 - alphas[go.get_namespace(go_id)]) * score\n            else:\n                annots[go_id] = (1 - alphas[go.get_namespace(go_id)]) * score\n        # Propagate scores with ontology structure\n        gos = list(annots.keys())\n        for go_id in gos:\n            for g_id in go.get_anchestors(go_id):\n                if g_id in annots:\n                    annots[g_id] = max(annots[g_id], annots[go_id])\n                else:\n                    annots[g_id] = annots[go_id]\n                \n        sannots = sorted(annots.items(), key=lambda x: x[1], reverse=True)\n        for go_id, score in sannots:\n            if score >= threshold:\n                w.write(prot_id + '\\t' + go_id + '\\t%.3f\\n' % score)\n        w.write('\\n')\nw.close()\ntotal_time = time.time() - start_time\nprint('Total prediction time for %d sequences is %d' % (total_seq, total_time))\n\n","metadata":{"execution":{"iopub.status.busy":"2023-05-15T16:27:19.902869Z","iopub.execute_input":"2023-05-15T16:27:19.903491Z","iopub.status.idle":"2023-05-15T16:40:52.552745Z","shell.execute_reply.started":"2023-05-15T16:27:19.903463Z","shell.execute_reply":"2023-05-15T16:40:52.551498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with open('submission.tsv') as f:\n    for i in range(10):\n        print(f.readline().strip())","metadata":{"execution":{"iopub.status.busy":"2023-05-15T16:40:55.832335Z","iopub.execute_input":"2023-05-15T16:40:55.832766Z","iopub.status.idle":"2023-05-15T16:40:55.839389Z","shell.execute_reply.started":"2023-05-15T16:40:55.832731Z","shell.execute_reply":"2023-05-15T16:40:55.838336Z"},"trusted":true},"execution_count":null,"outputs":[]}]}