GeneOntology.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8fb377a2",
   "metadata": {},
   "source": [
    "# Download \n",
    "•\tgene2go.gz from ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz  \n",
    "•\tEntrez Gene info from ftp://ftp.ncbi.nih.gov/gene/DATA/gene_info.gz    \n",
    "•\tgo from https://purl.obolibrary.org/obo/go/go-basic.obo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92299cdf",
   "metadata": {},
   "source": [
    "# Get GO_annotations-9606-inferred-allev.tsv from \n",
    "https://github.com/dhimmel/gene-ontology/blob/d57fd938f90c79152449f5cd23d3c438a19ac2f5/code/process.ipynb#L492 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fb2938aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections\n",
    "import os\n",
    "\n",
    "import pandas\n",
    "import networkx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f5d3980f",
   "metadata": {},
   "outputs": [],
   "source": [
    "download_dir ='../data/GeneOntology'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c48de713",
   "metadata": {},
   "outputs": [],
   "source": [
    "remove_subsets = {\n",
    "    'goantislim_grouping', # Grouping classes that can be excluded\n",
    "    'gocheck_do_not_annotate' # Term not to be used for direct annotation\n",
    "    'gocheck_do_not_manually_annotate', # Term not to be used for direct manual annotation\n",
    "}\n",
    "\n",
    "propagate_along = {'is_a', 'part_of'}\n",
    "\n",
    "experimental_codes = {\n",
    "    'EXP', # Inferred from Experiment\n",
    "    'IDA', # Inferred from Direct Assay\n",
    "    'IPI', # Inferred from Physical Interaction\n",
    "    'IMP', # Inferred from Mutant Phenotype\n",
    "    'IGI', # Inferred from Genetic Interaction\n",
    "    'IEP', # Inferred from Expression Pattern\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b6144047",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pip install obonet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3e4d84a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pip install networkx #== 3.2.1 # networkx==2.2 require early version of python but there's no native build for python below 3.8 for Apple Silicon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b37f7495",
   "metadata": {},
   "source": [
    "## Read Gene Ontology graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "88bb27d5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "42255\n",
      "81848\n"
     ]
    }
   ],
   "source": [
    "#https://github.com/dhimmel/gene-ontology/blob/d57fd938f90c79152449f5cd23d3c438a19ac2f5/code/process.ipynb\n",
    "import networkx\n",
    "import obonet # installation: https://github.com/dhimmel/obonet load OBO-formatted ontologies into networkx\n",
    "\n",
    "path = os.path.join(download_dir, 'go-basic.obo') #'https://purl.obolibrary.org/obo/go/go-basic.obo' #data-version: releases/2024-04-24\n",
    "with open(path) as read_file:\n",
    "    graph = obonet.read_obo(read_file)\n",
    "#networkx.info(graph)\n",
    "print(graph.number_of_nodes())\n",
    "print(graph.number_of_edges())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bcadd0ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>go_id</th>\n",
       "      <th>go_name</th>\n",
       "      <th>go_domain</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>GO:0000001</td>\n",
       "      <td>mitochondrion inheritance</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>GO:0000002</td>\n",
       "      <td>mitochondrial genome maintenance</td>\n",
       "      <td>biological_process</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        go_id                           go_name           go_domain\n",
       "0  GO:0000001         mitochondrion inheritance  biological_process\n",
       "1  GO:0000002  mitochondrial genome maintenance  biological_process"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#go_df = utilities.graph_to_dataframe(graph)\n",
    "rows = list()\n",
    "for node, data in graph.nodes(data=True):\n",
    "    rows.append((node, data['name'], data['namespace']))\n",
    "go_df = pandas.DataFrame(rows, columns=['go_id', 'go_name', 'go_domain'])\n",
    "go_df = go_df.sort_values('go_id')\n",
    "    \n",
    "go_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8a4fc80b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "42255\n",
      "73512\n"
     ]
    }
   ],
   "source": [
    "# Remove nodes that should not be annotated\n",
    "remove_nodes = set()\n",
    "for node, data in graph.nodes(data=True):\n",
    "    if remove_subsets & set(data.get('subset', [])):\n",
    "        remove_nodes.add(node)\n",
    "        #graph.remove_node(node)\n",
    "\n",
    "# Remove edges that should not be propagated along\n",
    "remove_edges = []\n",
    "for u, v, key in graph.edges(data=False, keys=True):\n",
    "    if key not in propagate_along:\n",
    "        remove_edges.append((u, v, key))\n",
    "\n",
    "for u, v, key in remove_edges:\n",
    "    graph.remove_edge(u, v, key)\n",
    "\n",
    "assert networkx.is_directed_acyclic_graph(graph)\n",
    "print(graph.number_of_nodes())\n",
    "print(graph.number_of_edges())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "22991fcb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'name': 'isocitrate dehydrogenase (NAD+) activity',\n",
       " 'namespace': 'molecular_function',\n",
       " 'def': '\"Catalysis of the reaction: isocitrate + NAD+ = 2-oxoglutarate + CO2 + NADH.\" [RHEA:23632]',\n",
       " 'synonym': ['\"isocitrate dehydrogenase (NAD) activity\" EXACT [EC:1.1.1.41]',\n",
       "  '\"isocitrate:NAD+ oxidoreductase (decarboxylating)\" EXACT [EC:1.1.1.41]',\n",
       "  '\"NAD dependent isocitrate dehydrogenase activity\" EXACT [EC:1.1.1.41]',\n",
       "  '\"NAD isocitrate dehydrogenase activity\" EXACT [EC:1.1.1.41]',\n",
       "  '\"NAD isocitric dehydrogenase activity\" EXACT [EC:1.1.1.41]',\n",
       "  '\"NAD-linked isocitrate dehydrogenase activity\" EXACT [EC:1.1.1.41]',\n",
       "  '\"NAD-specific isocitrate dehydrogenase activity\" EXACT [EC:1.1.1.41]',\n",
       "  '\"nicotinamide adenine dinucleotide isocitrate dehydrogenase activity\" EXACT [EC:1.1.1.41]'],\n",
       " 'xref': ['EC:1.1.1.41',\n",
       "  'MetaCyc:ISOCITRATE-DEHYDROGENASE-NAD+-RXN',\n",
       "  'Reactome:R-HSA-70967 \"isocitrate + NAD+ => alpha-ketoglutarate + CO2 + NADH + H+ [IDH3]\"',\n",
       "  'RHEA:23632'],\n",
       " 'is_a': ['GO:0004448']}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "graph.nodes['GO:0004449']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "760b19b8",
   "metadata": {},
   "source": [
    "## Read Entrez Gene and annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "fd98ad2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/sg/983t9lp12bsfyjv9ylh6002m0000gt/T/ipykernel_99531/176760046.py:26: DtypeWarning: Columns (15) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  gene_df = pandas.read_table(path, comment='#', names=column_names, na_values=['-'], dtype=dtype, index_col=False)\n",
      "/var/folders/sg/983t9lp12bsfyjv9ylh6002m0000gt/T/ipykernel_99531/176760046.py:26: ParserWarning: Length of header or names does not match length of data. This leads to a loss of data with index_col=False.\n",
      "  gene_df = pandas.read_table(path, comment='#', names=column_names, na_values=['-'], dtype=dtype, index_col=False)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>GeneID</th>\n",
       "      <th>Symbol</th>\n",
       "      <th>type_of_gene</th>\n",
       "      <th>tax_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>11465943</th>\n",
       "      <td>1</td>\n",
       "      <td>A1BG</td>\n",
       "      <td>protein-coding</td>\n",
       "      <td>9606</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11465944</th>\n",
       "      <td>2</td>\n",
       "      <td>A2M</td>\n",
       "      <td>protein-coding</td>\n",
       "      <td>9606</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          GeneID Symbol    type_of_gene  tax_id\n",
       "11465943       1   A1BG  protein-coding    9606\n",
       "11465944       2    A2M  protein-coding    9606"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read Entrez Gene info\n",
    "import pandas \n",
    "import os\n",
    "#download_dir ='../data'\n",
    "path = os.path.join(download_dir, 'gene_info.gz') #'ftp://ftp.ncbi.nih.gov/gene/DATA/gene_info.gz'\n",
    "column_names = [\n",
    "    'tax_id',\n",
    "    'GeneID',\n",
    "    'Symbol',\n",
    "    'LocusTag',\n",
    "    'Synonyms',\n",
    "    'dbXrefs',\n",
    "    'chromosome',\n",
    "    'map_location',\n",
    "    'description',\n",
    "    'type_of_gene',\n",
    "    'Symbol_from_nomenclature_authority',\n",
    "    'Full_name_from_nomenclature_authority',\n",
    "    'Nomenclature_status',\n",
    "    'Other_designations',\n",
    "    'Modification_date',\n",
    "]\n",
    "dtype = {x: str for x in column_names}\n",
    "for column in 'tax_id', 'GeneID':\n",
    "    dtype[column] = int\n",
    "gene_df = pandas.read_table(path, comment='#', names=column_names, na_values=['-'], dtype=dtype, index_col=False)\n",
    "gene_df = gene_df[['GeneID', 'Symbol', 'type_of_gene', 'tax_id']]\n",
    "gene_df = gene_df[gene_df['tax_id']==9606]\n",
    "gene_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "fabe0223",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>tax_id</th>\n",
       "      <th>GeneID</th>\n",
       "      <th>GO_ID</th>\n",
       "      <th>Evidence</th>\n",
       "      <th>Qualifier</th>\n",
       "      <th>GO_term</th>\n",
       "      <th>PubMed</th>\n",
       "      <th>Category</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>22059886</th>\n",
       "      <td>9606</td>\n",
       "      <td>1</td>\n",
       "      <td>GO:0002764</td>\n",
       "      <td>IBA</td>\n",
       "      <td>involved_in</td>\n",
       "      <td>immune response-regulating signaling pathway</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Process</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22059887</th>\n",
       "      <td>9606</td>\n",
       "      <td>1</td>\n",
       "      <td>GO:0003674</td>\n",
       "      <td>ND</td>\n",
       "      <td>enables</td>\n",
       "      <td>molecular_function</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Function</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          tax_id  GeneID       GO_ID Evidence    Qualifier  \\\n",
       "22059886    9606       1  GO:0002764      IBA  involved_in   \n",
       "22059887    9606       1  GO:0003674       ND      enables   \n",
       "\n",
       "                                               GO_term PubMed  Category  \n",
       "22059886  immune response-regulating signaling pathway    NaN   Process  \n",
       "22059887                            molecular_function    NaN  Function  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Download gene2go.gz from ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz\n",
    "#Unzip it\n",
    "#Read annotations #https://github.com/dhimmel/gene-ontology/blob/d57fd938f90c79152449f5cd23d3c438a19ac2f5/code/utilities.py\n",
    "import pandas \n",
    "import os\n",
    "\n",
    "path = os.path.join(download_dir, 'gene2go.gz') #'ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz'\n",
    "column_names = ['tax_id','GeneID','GO_ID','Evidence','Qualifier','GO_term','PubMed','Category']\n",
    "goa_df = pandas.read_table(path, comment='#', names=column_names, na_values=['-'], dtype=None, index_col=False)\n",
    "goa_df = goa_df[goa_df['tax_id']==9606]\n",
    "goa_df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "251bbe94",
   "metadata": {},
   "source": [
    "## Add and propagate annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a5ee7d0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_NOT_qaulifier(qualifier):\n",
    "    \"\"\"\n",
    "    Returns whether a `Qualifier` in gene2go is a NOT qualifier.\n",
    "    http://geneontology.org/page/go-annotation-conventions#not\n",
    "    \"\"\"\n",
    "    if pandas.isnull(qualifier):\n",
    "        return False\n",
    "    if not qualifier:\n",
    "        return False\n",
    "    if qualifier.upper().startswith('NOT'):\n",
    "        return True\n",
    "    return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2c4e2c5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def annotate_graph(graph, goa_df):\n",
    "    \"\"\"Add direct annotations to graph\"\"\"\n",
    "    graph = graph.copy()\n",
    "    \n",
    "    # Add dictionary items for storing annotations\n",
    "    for node, data in graph.nodes.items():\n",
    "        for key in 'direct_annotations', 'direct_not_annotations', 'inferred_annotations':\n",
    "            data[key] = set()\n",
    "\n",
    "    # Populate direct annotations\n",
    "    for i, row in goa_df.iterrows():\n",
    "\n",
    "        go_id = row['GO_ID']\n",
    "        if go_id not in graph:\n",
    "            continue\n",
    "\n",
    "        key = 'direct_not_annotations' if is_NOT_qaulifier(row.Qualifier) else 'direct_annotations'\n",
    "\n",
    "        gene = row['GeneID']\n",
    "        graph.nodes[go_id][key].add(gene) #https://networkx.org/documentation/networkx-2.2/tutorial.html#node-attributes\n",
    "    \n",
    "    return graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8196953e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def propagate_annotations(graph):\n",
    "    \"\"\"Infer annotations through propagations\"\"\"\n",
    "    for node in networkx.topological_sort(graph):\n",
    "        data = graph.nodes[node]\n",
    "        inferred = data['inferred_annotations']\n",
    "        inferred -= data['direct_not_annotations']\n",
    "        inferred |= data['direct_annotations']\n",
    "        for subsuming_node in graph.successors(node):\n",
    "            subsuming_data = graph.nodes[subsuming_node]\n",
    "            subsuming_data['inferred_annotations'] |= inferred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b021d1f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "joiner = lambda x: '|'.join(map(str, x))\n",
    "\n",
    "def extract_annotation_df(graph):\n",
    "    \"\"\"Create an annotation dataframe\"\"\"\n",
    "    rows = list()\n",
    "    for node, data in graph.nodes.items():\n",
    "        if node in remove_nodes:\n",
    "            continue\n",
    "        for kind in 'direct', 'inferred':\n",
    "            for gene in data['{}_annotations'.format(kind)]:\n",
    "                rows.append((node, kind, gene))\n",
    "    \n",
    "    annotation_df = pandas.DataFrame(rows, columns=['go_id', 'kind', 'GeneID'])\n",
    "    annotation_df = annotation_df.merge(gene_df)\n",
    "\n",
    "    rows = list()\n",
    "    for (tax_id, kind), taxon_df in annotation_df.groupby(['tax_id', 'kind']):\n",
    "        for go_id, term_df in taxon_df.groupby('go_id'):\n",
    "            term_df = term_df.sort_values('GeneID')\n",
    "            row = tax_id, go_id, kind, len(term_df), joiner(term_df['GeneID']), joiner(term_df['Symbol'])\n",
    "            rows.append(row)\n",
    "    wide_df = pandas.DataFrame(rows, columns = ['tax_id', 'go_id', 'annotation_type', 'size', 'gene_ids', 'gene_symbols'])\n",
    "    wide_df = go_df.merge(wide_df)\n",
    "    return wide_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e788c6f8",
   "metadata": {},
   "source": [
    "## Extract and save annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "30c53142",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\ngraph_annot = annotate_graph(graph, goa_df)\\n\\npropagate_annotations(graph_annot)\\n\\ngraph_annot.nodes['GO:0004449']['direct_annotations']\\n\\nrows = list()\\nfor node, data in graph_annot.nodes.items():\\n    if node in remove_nodes:\\n        continue\\n    for gene in data['inferred_annotations']:\\n        kind ='inferred_annotations'\\n        rows.append((node, kind, gene))\\n\\nannotation_df = pandas.DataFrame(rows, columns=['go_id', 'kind', 'GeneID'])\\nannotation_df = annotation_df.merge(gene_df)\\nannotation_df\\n\""
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "graph_annot = annotate_graph(graph, goa_df)\n",
    "\n",
    "propagate_annotations(graph_annot)\n",
    "\n",
    "graph_annot.nodes['GO:0004449']['direct_annotations']\n",
    "\n",
    "rows = list()\n",
    "for node, data in graph_annot.nodes.items():\n",
    "    if node in remove_nodes:\n",
    "        continue\n",
    "    for gene in data['inferred_annotations']:\n",
    "        kind ='inferred_annotations'\n",
    "        rows.append((node, kind, gene))\n",
    "\n",
    "annotation_df = pandas.DataFrame(rows, columns=['go_id', 'kind', 'GeneID'])\n",
    "annotation_df = annotation_df.merge(gene_df)\n",
    "annotation_df\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1035c26f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>go_id</th>\n",
       "      <th>go_name</th>\n",
       "      <th>go_domain</th>\n",
       "      <th>tax_id</th>\n",
       "      <th>annotation_type</th>\n",
       "      <th>size</th>\n",
       "      <th>gene_ids</th>\n",
       "      <th>gene_symbols</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>GO:0000002</td>\n",
       "      <td>mitochondrial genome maintenance</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>direct</td>\n",
       "      <td>11</td>\n",
       "      <td>291|1890|4205|4358|4976|9361|10000|55186|80119...</td>\n",
       "      <td>SLC25A4|TYMP|MEF2A|MPV17|OPA1|LONP1|AKT3|SLC25...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>GO:0000002</td>\n",
       "      <td>mitochondrial genome maintenance</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>inferred</td>\n",
       "      <td>26</td>\n",
       "      <td>142|291|1763|1890|2021|4205|4358|4976|5428|624...</td>\n",
       "      <td>PARP1|SLC25A4|DNA2|TYMP|ENDOG|MEF2A|MPV17|OPA1...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>GO:0000009</td>\n",
       "      <td>alpha-1,6-mannosyltransferase activity</td>\n",
       "      <td>molecular_function</td>\n",
       "      <td>9606</td>\n",
       "      <td>direct</td>\n",
       "      <td>2</td>\n",
       "      <td>55650|79087</td>\n",
       "      <td>PIGV|ALG12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>GO:0000009</td>\n",
       "      <td>alpha-1,6-mannosyltransferase activity</td>\n",
       "      <td>molecular_function</td>\n",
       "      <td>9606</td>\n",
       "      <td>inferred</td>\n",
       "      <td>2</td>\n",
       "      <td>55650|79087</td>\n",
       "      <td>PIGV|ALG12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GO:0000012</td>\n",
       "      <td>single strand break repair</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>direct</td>\n",
       "      <td>10</td>\n",
       "      <td>1161|2074|3981|7141|7515|23411|54840|55775|200...</td>\n",
       "      <td>ERCC8|ERCC6|LIG4|TNP1|XRCC1|SIRT1|APTX|TDP1|AP...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40975</th>\n",
       "      <td>GO:2001304</td>\n",
       "      <td>lipoxin B4 metabolic process</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>inferred</td>\n",
       "      <td>1</td>\n",
       "      <td>239</td>\n",
       "      <td>ALOX12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40976</th>\n",
       "      <td>GO:2001306</td>\n",
       "      <td>lipoxin B4 biosynthetic process</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>direct</td>\n",
       "      <td>1</td>\n",
       "      <td>239</td>\n",
       "      <td>ALOX12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40977</th>\n",
       "      <td>GO:2001306</td>\n",
       "      <td>lipoxin B4 biosynthetic process</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>inferred</td>\n",
       "      <td>1</td>\n",
       "      <td>239</td>\n",
       "      <td>ALOX12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40978</th>\n",
       "      <td>GO:2001311</td>\n",
       "      <td>lysobisphosphatidic acid metabolic process</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>direct</td>\n",
       "      <td>2</td>\n",
       "      <td>51205|57406</td>\n",
       "      <td>ACP6|ABHD6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40979</th>\n",
       "      <td>GO:2001311</td>\n",
       "      <td>lysobisphosphatidic acid metabolic process</td>\n",
       "      <td>biological_process</td>\n",
       "      <td>9606</td>\n",
       "      <td>inferred</td>\n",
       "      <td>2</td>\n",
       "      <td>51205|57406</td>\n",
       "      <td>ACP6|ABHD6</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>40980 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            go_id                                     go_name  \\\n",
       "0      GO:0000002            mitochondrial genome maintenance   \n",
       "1      GO:0000002            mitochondrial genome maintenance   \n",
       "2      GO:0000009      alpha-1,6-mannosyltransferase activity   \n",
       "3      GO:0000009      alpha-1,6-mannosyltransferase activity   \n",
       "4      GO:0000012                  single strand break repair   \n",
       "...           ...                                         ...   \n",
       "40975  GO:2001304                lipoxin B4 metabolic process   \n",
       "40976  GO:2001306             lipoxin B4 biosynthetic process   \n",
       "40977  GO:2001306             lipoxin B4 biosynthetic process   \n",
       "40978  GO:2001311  lysobisphosphatidic acid metabolic process   \n",
       "40979  GO:2001311  lysobisphosphatidic acid metabolic process   \n",
       "\n",
       "                go_domain  tax_id annotation_type  size  \\\n",
       "0      biological_process    9606          direct    11   \n",
       "1      biological_process    9606        inferred    26   \n",
       "2      molecular_function    9606          direct     2   \n",
       "3      molecular_function    9606        inferred     2   \n",
       "4      biological_process    9606          direct    10   \n",
       "...                   ...     ...             ...   ...   \n",
       "40975  biological_process    9606        inferred     1   \n",
       "40976  biological_process    9606          direct     1   \n",
       "40977  biological_process    9606        inferred     1   \n",
       "40978  biological_process    9606          direct     2   \n",
       "40979  biological_process    9606        inferred     2   \n",
       "\n",
       "                                                gene_ids  \\\n",
       "0      291|1890|4205|4358|4976|9361|10000|55186|80119...   \n",
       "1      142|291|1763|1890|2021|4205|4358|4976|5428|624...   \n",
       "2                                            55650|79087   \n",
       "3                                            55650|79087   \n",
       "4      1161|2074|3981|7141|7515|23411|54840|55775|200...   \n",
       "...                                                  ...   \n",
       "40975                                                239   \n",
       "40976                                                239   \n",
       "40977                                                239   \n",
       "40978                                        51205|57406   \n",
       "40979                                        51205|57406   \n",
       "\n",
       "                                            gene_symbols  \n",
       "0      SLC25A4|TYMP|MEF2A|MPV17|OPA1|LONP1|AKT3|SLC25...  \n",
       "1      PARP1|SLC25A4|DNA2|TYMP|ENDOG|MEF2A|MPV17|OPA1...  \n",
       "2                                             PIGV|ALG12  \n",
       "3                                             PIGV|ALG12  \n",
       "4      ERCC8|ERCC6|LIG4|TNP1|XRCC1|SIRT1|APTX|TDP1|AP...  \n",
       "...                                                  ...  \n",
       "40975                                             ALOX12  \n",
       "40976                                             ALOX12  \n",
       "40977                                             ALOX12  \n",
       "40978                                         ACP6|ABHD6  \n",
       "40979                                         ACP6|ABHD6  \n",
       "\n",
       "[40980 rows x 8 columns]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Get GO_annotations-9606-inferred-allev.tsv from https://github.com/dhimmel/gene-ontology/blob/d57fd938f90c79152449f5cd23d3c438a19ac2f5/code/process.ipynb#L492 \n",
    "env_type='allev'\n",
    "graph_annot = annotate_graph(graph, goa_df)\n",
    "propagate_annotations(graph_annot)\n",
    "annotation_df = extract_annotation_df(graph_annot)\n",
    "annotation_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1299907e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = annotation_df[(annotation_df['tax_id']==9606)&(annotation_df['annotation_type']=='inferred')]\n",
    "path = os.path.join(download_dir, 'GO_annotations-9606-inferred-allev.tsv') \n",
    "df.to_csv(path, sep='\\t', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9239b1c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "for ev_type in 'allev', 'expev':\n",
    "    goa_subset_df = goa_df\n",
    "    if ev_type == 'expev':\n",
    "        goa_subset_df = goa_subset_df[goa_subset_df.Evidence.isin(experimental_codes)]\n",
    "    graph_annot = annotate_graph(graph, goa_subset_df)\n",
    "    propagate_annotations(graph_annot)\n",
    "    annotation_df = extract_annotation_df(graph_annot)\n",
    "\n",
    "    for (tax_id, annotation_type), df in annotation_df.groupby(['tax_id', 'annotation_type']):\n",
    "        path = utilities.get_annotation_path(annotation_dir, tax_id, annotation_type, ev_type, mkdir=True)\n",
    "        print(path)\n",
    "        df.to_csv(path, sep='\\t', index=False)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27ece0f9",
   "metadata": {},
   "source": [
    "## Compare with previous version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d60f687a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "07830b36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(22283, 8)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# stats\n",
    "import pandas as pd\n",
    "df = pd.read_table('../data/GeneOntology/GO_annotations-9606-inferred-allev.tsv')\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "8dda54b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/sg/983t9lp12bsfyjv9ylh6002m0000gt/T/ipykernel_99531/2484383344.py:4: DtypeWarning: Columns (4,15) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  gene_df = pd.read_table(read_file, na_values=['-'])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "20608"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import gzip\n",
    "path = '../data/EntrezGene/Homo_sapiens.gene_info.gz'\n",
    "read_file = gzip.open(path, 'rb')\n",
    "gene_df = pd.read_table(read_file, na_values=['-'])\n",
    "gene_df = gene_df.rename(columns={'#tax_id':'tax_id'})\n",
    "gene_df = gene_df.query('tax_id == 9606')\n",
    "columns = ['tax_id', 'GeneID', 'Symbol', 'chromosome', 'map_location', 'type_of_gene', 'description']\n",
    "gene_df = gene_df[columns]\n",
    "#gene_df.to_csv('data/genes-human.tsv', sep='\\t', index=False)\n",
    "gene_df = gene_df[gene_df.type_of_gene == 'protein-coding']\n",
    "coding_genes = set(gene_df.GeneID)\n",
    "len(coding_genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ac6ea7b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "go_domain\n",
       "Biological Process    12322\n",
       "Cellular Component     1695\n",
       "Molecular Function     3460\n",
       "Name: go_id, dtype: int64"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "edge\n",
       "Gene–participates–Biological Process    548342\n",
       "Gene–participates–Cellular Component     88885\n",
       "Gene–participates–Molecular Function    104777\n",
       "dtype: int64"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes =[]\n",
    "edges =[]\n",
    "for i, row in df.iterrows():\n",
    "    genes = coding_genes & set(map(int, row.gene_ids.split('|')))\n",
    "    if 2 > len(genes) or len(genes) > 1000:\n",
    "        continue\n",
    "    kind = row['go_domain'].replace('_', ' ').title()\n",
    "    identifier = row['go_id']\n",
    "    name=row['go_name']\n",
    "    nodes.append([identifier, name, kind])\n",
    "    for gene in genes:\n",
    "        source_id = gene\n",
    "        target_id = row['go_id']\n",
    "        edge = 'Gene–participates–' + kind\n",
    "        edges.append([source_id, target_id, edge])\n",
    "        \n",
    "nodes_df = pd.DataFrame(nodes, columns=['go_id', 'go_name', 'go_domain'])    \n",
    "nodes_df.groupby('go_domain')['go_id'].nunique()\n",
    "\n",
    "edges_df = pd.DataFrame(edges, columns=['gene', 'go_id', 'edge'])    \n",
    "edges_df.groupby('edge').size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0355ff0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#compare with hetionet version  438 mesh terms \n",
    "df_hetionet = pd.read_table('https://raw.githubusercontent.com/dhimmel/gene-ontology/87bab297f55db283e65a7a984607316b409415ae/annotations/taxid_9606/GO_annotations-9606-inferred-allev.tsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "b8c4528e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20945"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_df = pd.read_table('https://raw.githubusercontent.com/dhimmel/entrez-gene/a7362748a34211e5df6f2d185bb3246279760546/data/genes-human.tsv')\n",
    "gene_df = gene_df[gene_df.type_of_gene == 'protein-coding']\n",
    "coding_genes = set(gene_df.GeneID)\n",
    "len(coding_genes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "5c363850",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "go_domain\n",
       "Biological Process    11381\n",
       "Cellular Component     1391\n",
       "Molecular Function     2884\n",
       "Name: go_id, dtype: int64"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "edge\n",
       "Gene–participates–Biological Process    559504\n",
       "Gene–participates–Cellular Component     73566\n",
       "Gene–participates–Molecular Function     97222\n",
       "dtype: int64"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes =[]\n",
    "edges =[]\n",
    "for i, row in df_hetionet.iterrows():\n",
    "    genes = coding_genes & set(map(int, row.gene_ids.split('|')))\n",
    "    if 2 > len(genes) or len(genes) > 1000:\n",
    "        continue\n",
    "    kind = row['go_domain'].replace('_', ' ').title()\n",
    "    identifier = row['go_id']\n",
    "    name=row['go_name']\n",
    "    nodes.append([identifier, name, kind])\n",
    "    for gene in genes:\n",
    "        source_id = gene\n",
    "        target_id = row['go_id']\n",
    "        edge = 'Gene–participates–' + kind\n",
    "        edges.append([source_id, target_id, edge])\n",
    "        \n",
    "nodes_df = pd.DataFrame(nodes, columns=['go_id', 'go_name', 'go_domain'])    \n",
    "nodes_df.groupby('go_domain')['go_id'].nunique()\n",
    "\n",
    "edges_df = pd.DataFrame(edges, columns=['gene', 'go_id', 'edge'])    \n",
    "edges_df.groupby('edge').size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62f643bc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}