notebooks/Data_Preparation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "***\n",
    "***\n",
    "\n",
    "<img width='700' src=\"https://user-images.githubusercontent.com/8030363/108961534-b9a66980-7634-11eb-96e2-cc46589dcb8c.png\" style=\"vertical-align:middle\">\n",
    "\n",
    "## Pre-Knowledge Graph Build Data Preparation\n",
    "***\n",
    "\n",
    "**Author:** [TJCallahan](https://mail.google.com/mail/u/0/?view=cm&fs=1&tf=1&to=callahantiff@gmail.com)  \n",
    "**GitHub Repository:** [PheKnowLator](https://github.com/callahantiff/PheKnowLator/wiki)  \n",
    "**Release:** **[v2.0.0](https://github.com/callahantiff/PheKnowLator/wiki/v2.0.0)**\n",
    "  \n",
    "<br>  \n",
    "  \n",
    "**Purpose:** This notebook serves as a script to download and process data in order to generate mapping and filtering data needed to build edges for the PheKnowLator knowledge graph. For more information on the data sources utilize within this script, please see the [Data Sources](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources) Wiki page.\n",
    "\n",
    "<br>\n",
    "\n",
    "**Assumptions:**   \n",
    "- Raw data downloads ➞ `./resources/processed_data/unprocessed_data`    \n",
    "- Processed data write location ➞ `./resources/processed_data`  \n",
    "\n",
    "<br>\n",
    "\n",
    "**Dependencies:**   \n",
    "- **Scripts**: This notebook utilizes several helper functions, which are stored in the [`data_utils.py`](https://github.com/callahantiff/PheKnowLator/blob/master/pkt_kg/utils/data_utils.py) and [`kg_utils.py`](https://github.com/callahantiff/PheKnowLator/blob/master/pkt_kg/utils/kg_utils.py) scripts.  \n",
    "- **Data**: Hyperlinks to all downloaded and generated data sources are provided through [Zenodo](https://zenodo.org/doi/10.5281/zenodo.7030039). <u>This notebook will download everything that is needed for you</u>.  \n",
    "_____\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "***\n",
    "\n",
    "### [Create Identifier Maps ](#create-identifier-maps)  \n",
    "- [HUMAN TRANSCRIPT, GENE, AND PROTEIN IDENTIFIER MAPPING](#human-transcript,-gene,-and-protein-identifier-mapping)\n",
    "  - [Entrez Gene-Ensembl Transcript](#entrezgene-ensembltranscript)  \n",
    "  - [Entrez Gene-Protein Ontology](#entrezgene-proteinontology)  \n",
    "  - [Ensembl Gene-Entrez Gene](#ensemblgene-entrezgene)\n",
    "  - [Gene Symbol-Ensembl Transcript](#genesymbol-ensembltranscript)  \n",
    "  - [STRING-Protein Ontology](#string-proteinontology)  \n",
    "  - [Uniprot Accession-Protein Ontology](#uniprotaccession-proteinontology)\n",
    "  \n",
    "\n",
    "- [OTHER IDENTIFIER MAPPING](#other-identifier-mapping) \n",
    "  - [ChEBI Identifiers](#mesh-chebi) \n",
    "  - [Human Disease and Phenotype Identifiers](#disease-identifiers)\n",
    "  - [Human Protein Atlas Tissue and Cell Types](#hpa-uberon)  \n",
    "  - [Reactome Pathways - Pathway Ontology](#reactome-pw)  \n",
    "  - [Genomic Identifiers - Sequence Ontology](#genomic-soo)  \n",
    "\n",
    "\n",
    "### [Create Edge Datasets](#create-edge-datasets)\n",
    "- [ONTOLOGIES](#ontologies)  \n",
    "  - [Protein Ontology](#protein-ontology)  \n",
    "  - [Relations Ontology](#relations-ontology)  \n",
    "\n",
    "\n",
    "- [LINKED DATA](#linked-data)  \n",
    "  - [Clinvar Variant-Diseases and Phenotypes](#clinvar-variant)\n",
    "  - [Uniprot Protein-Cofactor and Protein-Catalyst](#uniprot-protein-cofactorcatalyst)  \n",
    "\n",
    "\n",
    "### [Create Instance Data and/or Subclass Metadata](#create-instance-metadata)  \n",
    "- [Genes/RNA](#gene-and-rna-metadata)\n",
    "- [Pathways](#pathway-metadata)\n",
    "- [Variants](#variant-metadata) \n",
    "- [Relations](#relations-metadata) \n",
    "\n",
    "____"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set-Up Environment\n",
    "_____"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # uncomment and run to install any required modules from notebooks/requirements.txt\n",
    "# import sys\n",
    "# !{sys.executable} -m pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if running a local version of pkt_kg, uncomment the code below\n",
    "# import sys\n",
    "# sys.path.append('../')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import needed libraries\n",
    "import datetime\n",
    "import glob\n",
    "import itertools\n",
    "import networkx\n",
    "import numpy\n",
    "import os\n",
    "import openpyxl\n",
    "import pandas\n",
    "import pickle\n",
    "import re\n",
    "import requests\n",
    "import sys\n",
    "\n",
    "from collections import Counter\n",
    "from functools import reduce\n",
    "from rdflib import Graph, Namespace, URIRef, BNode, Literal\n",
    "from rdflib.namespace import OWL, RDF, RDFS\n",
    "from reactome2py import content\n",
    "from tqdm import tqdm\n",
    "from typing import Dict\n",
    "\n",
    "from pkt_kg.utils import *  # import pkt_kg utility script containing helper functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define Global Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# directory to use for processing data\n",
    "unprocessed_data_location = '../resources/processed_data/unprocessed_data/'\n",
    "processed_data_location = '../resources/processed_data/'\n",
    "\n",
    "# directory to write relations data to\n",
    "relations_data_location = '../resources/relations_data/'\n",
    "\n",
    "# directory to write node metadata to\n",
    "node_data_location = '../resources/node_data/'\n",
    "\n",
    "# directory to write kg construction approach dictionary to\n",
    "construction_approach_location = '../resources/construction_approach/'\n",
    "\n",
    "# directory to write ontology data to\n",
    "ontology_data_location = '../resources/ontologies/'\n",
    "\n",
    "# owltools location\n",
    "owltools_location = '../pkt_kg/libs/owltools'\n",
    "\n",
    "# obo spacespace\n",
    "obo = Namespace('http://purl.obolibrary.org/obo/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "***\n",
    "### CREATE MAPPING DATASETS  <a class=\"anchor\" id=\"create-identifier-maps\"></a>\n",
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Human Transcript, Gene, and Protein Identifier Mapping  <a class=\"anchor\" id=\"human-transcript,-gene,-and-protein-identifier-mapping\"></a>\n",
    "***\n",
    "\n",
    "**Data Source Wiki Pages:**   \n",
    "- [Ensembl](https://uswest.ensembl.org/)  \n",
    "- [Uniprot Knowledgebase](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/#uniprot-knowledgebase)  \n",
    "- [HGNC](ftp://ftp.ebi.ac.uk/pub/databases/genenames/new/tsv/hgnc_complete_set.txt) \n",
    "- [NCBI Gene](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/#ncbi-gene) \n",
    "- [Protein Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/#protein-ontology)\n",
    "\n",
    "<br>\n",
    "\n",
    "**Purpose:** To map create `protein-coding gene`-`protein` edges and mappings between the identifiers types listed below. The edges types produced from each of these mappings will be further described within each of the subsequent identifier mapping sections:  \n",
    "- [Entrez Gene-Ensembl Transcript](#entrezgene-ensembltranscript)  \n",
    "- [Entrez Gene-Protein Ontology](#entrezgene-proteinontology)  \n",
    "- [Ensembl Gene-Entrez Gene](#ensemblgene-entrezgene)\n",
    "- [Gene Symbol-Ensembl Transcript](#genesymbol-ensembltranscript)  \n",
    "- [STRING-Protein Ontology](#string-proteinontology)  \n",
    "- [Uniprot Accession-Protein Ontology](#uniprotaccession-proteinontology)\n",
    "\n",
    "<br>\n",
    "\n",
    "**Gene and Transcript Types:** The transcript and gene/locus types were reviewed by a PhD Molecular biologist to confirm whether or not they should be classified as `protein-coding` or not, which is useful for creating `genomic`-`rna`, `genomic`-`protein`, and `rna`-`protein` edges in the knowledge graph. For more information on this classification, please see the table below. Definitions of concepts in the table have been taken from [HGNC](https://www.genenames.org/help/symbol-report/), [Ensembl](https://uswest.ensembl.org/info/genome/genebuild/biotypes.html), [NCBI](https://www.ncbi.nlm.nih.gov/IEB/ToolBox/CPP_DOC/lxr/source/src/objects/entrezgene/entrezgene.asn), and Wikipedia.\n",
    "\n",
    "<table>\n",
    "<th align=\"center\">Gene and Transcript Type</th>  \n",
    "<th align=\"center\">Definition</th>\n",
    "<th align=\"center\">Type</th>\n",
    "<th align=\"center\">Genomic material <i>transcribed_to</i> RNA</th>\n",
    "<th align=\"center\">RNA <i>translated_to</i> Protein</th>\n",
    "<th align=\"center\">Genomic material <i>has_gene_product</i> Protein</th>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">biological-region</td> \n",
    "  <td rowspan=\"2\">Biological_region (SO:0001411); Special note: This is a parental feature spanning all other feature annotation on each RefSeq Functional Element record. It is a 'misc_feature' in GenBank flat files but a 'Region' feature in ASN.1 and GFF3 formats</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td>transcript</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_C_gene</td> \n",
    "  <td rowspan=\"2\">Constant chain immunoglobulin gene that undergoes somatic recombination before transcription</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_C_pseudogene</td> \n",
    "  <td rowspan=\"2\">Inactivated immunoglobulin gene. Immunoglobulin gene segments that are inactivated due to frameshift mutations and/or stop codons in the open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\t \t \n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_D_gene</td> \n",
    "  <td rowspan=\"2\">Diversity chain immunoglobulin gene that undergoes somatic recombination before transcription</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_J_gene</td> \n",
    "  <td rowspan=\"2\">IG J gene: Joining chain immunoglobulin gene that undergoes somatic recombination before transcription</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_J_pseudogene</td> \n",
    "  <td rowspan=\"2\">Inactivated immunoglobulin gene. Immunoglobulin gene segments that are inactivated due to frameshift mutations and/or stop codons in the open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_pseudogene</td> \n",
    "  <td rowspan=\"2\">Inactivated immunoglobulin gene. Immunoglobulin gene segments that are inactivated due to frameshift mutations and/or stop codons in the open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_V_gene</td> \n",
    "  <td rowspan=\"2\">Variable chain immunoglobulin gene that undergoes somatic recombination before transcription</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">IG_V_pseudogene</td> \n",
    "  <td rowspan=\"2\">Inactivated immunoglobulin gene. Immunoglobulin gene segments that are inactivated due to frameshift mutations and/or stop codons in the open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">lncRNA</td> \n",
    "  <td rowspan=\"2\">RNA, long non-coding - non-protein coding genes that encode long non-coding RNAs (lncRNAs) (SO:0001877); these are at least 200 nt in length. Subtypes include intergenic (SO:0001463), intronic (SO:0001903) and antisense (SO:0001904)</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">miRNA</td> \n",
    "  <td rowspan=\"2\">RNA, micro - non-protein coding genes that encode microRNAs (miRNAs) (SO:0001265)</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">misc_RNA</td> \n",
    "  <td rowspan=\"2\">Non-protein coding genes that encode miscellaneous types of small ncRNAs, such as vault (SO:0000404) and Y (SO:0000405) RNA genes</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">Mt_rRNA</td> \n",
    "  <td rowspan=\"2\">Mitochondrial rRNA</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">Mt_tRNA</td> \n",
    "  <td rowspan=\"2\">Mitochondrial tRNA</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">ncRNA</td> \n",
    "  <td rowspan=\"2\">Noncoding RNA</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">non_stop_decay</td> \n",
    "  <td rowspan=\"2\">Transcripts that have polyA features (including signal) without a prior stop codon in the CDS, i.e. a non-genomic polyA tail attached directly to the CDS without 3' UTR. These transcripts are subject to degradation</td>\n",
    "  <td>gene</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">nonsense_mediated_decay</td> \n",
    "  <td rowspan=\"2\">If the coding sequence (following the appropriate reference) of a transcript finishes >50bp from a downstream splice site then it is tagged as NMD. If the variant does not cover the full reference coding sequence then it is annotated as NMD if NMD is unavoidable i.e. no matter what the exon structure of the missing portion is the transcript will be subject to NMD</td>\n",
    "  <td>gene</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">other</td> \n",
    "  <td rowspan=\"2\">other</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">phenotype</td> \n",
    "  <td rowspan=\"2\"> Mapped phenotypes where the causative gene has not been identified (SO:0001500) </td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">polymorphic_pseudogene</td> \n",
    "  <td rowspan=\"2\">Pseudogene owing to a SNP/DIP but in other individuals/haplotypes/strains the gene is translated</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">processed_pseudogene</td> \n",
    "  <td rowspan=\"2\">Pseudogene that lack introns and is thought to arise from reverse transcription of mRNA followed by reinsertion of DNA into the genome</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">processed_transcript</td> \n",
    "  <td rowspan=\"2\">Gene/transcript that doesn't contain an open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">protein_coding</td> \n",
    "  <td rowspan=\"2\">Contains an open reading frame (ORF)</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>yes</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">pseudogene</td> \n",
    "  <td rowspan=\"2\">Have homology to proteins but generally suffer from a disrupted coding sequence and an active homologous gene can be found at another locus. Sometimes these entries have an intact coding sequence or an open but truncated open reading frame, in which case there is other evidence used (for example genomic polyA stretches at the 3' end) to classify them as a pseudogene</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">retained_intron</td> \n",
    "  <td rowspan=\"2\">Has an alternatively spliced transcript believed to contain intronic sequence relative to other, coding, variants</td>\n",
    "  <td>gene</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">ribozyme</td> \n",
    "  <td rowspan=\"2\">Ribozymes are RNA molecules that have the ability to catalyze specific biochemical reactions, including RNA splicing in gene expression, similar to the action of protein enzymes</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">rRNA</td> \n",
    "  <td rowspan=\"2\">RNA, ribosomal - non-protein coding genes that encode ribosomal RNAs (rRNAs) (SO:0001637)</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">rRNA_pseudogene</td> \n",
    "  <td rowspan=\"2\">A gene that has homology to known protein-coding genes but contain a frameshift and/or stop codon(s) which disrupts the open reading frame. Thought to have arisen through duplication followed by loss of function</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">scaRNA</td> \n",
    "  <td rowspan=\"2\">Small Cajal body-specific RNAs are a class of small nucleolar RNAs that specifically localize to the Cajal body, a nuclear organelle involved in the biogenesis of small nuclear ribonucleoproteins/td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">scRNA</td> \n",
    "  <td rowspan=\"2\">RNA, small cytoplasmic - non-protein coding genes that encode small cytoplasmic RNAs (scRNAs) (SO:0001266)</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">snoRNA</td> \n",
    "  <td rowspan=\"2\">RNA, small nucleolar - non-protein coding genes that encode small nucleolar RNAs (snoRNAs) containing C/D or H/ACA box domains (SO:0001267)</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">snRNA</td> \n",
    "  <td rowspan=\"2\">RNA, small nuclear - non-protein coding genes that encode small nuclear RNAs (snRNAs) (SO:0001268)</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">sRNA</td> \n",
    "  <td rowspan=\"2\">Bacterial small RNAs (sRNA) are small RNAs produced by bacteria; they are 50- to 500-nucleotide non-coding RNA molecules, highly structured and containing several stem-loops</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">TEC</td> \n",
    "  <td rowspan=\"2\">TEC (To be Experimentally Confirmed). This is used for non-spliced EST clusters that have polyA features. This category has been specifically created for the ENCODE project to highlight regions that could indicate the presence of protein coding genes that require experimental validation, either by 5' RACE or RT-PCR to extend the transcripts, or by confirming expression of the putatively-encoded peptide with specific antibodies</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>yes</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">TR_C_gene</td> \n",
    "  <td rowspan=\"2\">Constant chain T cell receptor gene that undergoes somatic recombination before transcription/td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">TR_D_gene</td> \n",
    "  <td rowspan=\"2\">Diversity chain T cell receptor gene that undergoes somatic recombination before transcription/td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">TR_J_gene</td> \n",
    "  <td rowspan=\"2\">Joining chain T cell receptor gene that undergoes somatic recombination before transcription</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">TR_J_pseudogene</td> \n",
    "  <td rowspan=\"2\">T cell receptor pseudogene - T cell receptor gene segments that are inactivated due to frameshift mutations and/or stop codons in the open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">TR_V_gene</td> \n",
    "  <td rowspan=\"2\">Variable chain T cell receptor gene that undergoes somatic recombination before transcription</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">TR_V_pseudogene</td> \n",
    "  <td rowspan=\"2\">T cell receptor pseudogene - T cell receptor gene segments that are inactivated due to frameshift mutations and/or stop codons in the open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">transcribed_processed_pseudogene</td> \n",
    "  <td rowspan=\"2\">Pseudogene where protein homology or genomic structure indicates a pseudogene, but the presence of locus-specific transcripts indicates expression. These can be classified into 'Processed', 'Unprocessed' and 'Unitary'</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">transcribed_unitary_pseudogene</td> \n",
    "  <td rowspan=\"2\">Pseudogene where protein homology or genomic structure indicates a pseudogene, but the presence of locus-specific transcripts indicates expression. These can be classified into 'Processed', 'Unprocessed' and 'Unitary'</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">transcribed_unprocessed_pseudogene</td> \n",
    "  <td rowspan=\"2\">Pseudogene where protein homology or genomic structure indicates a pseudogene, but the presence of locus-specific transcripts indicates expression. These can be classified into 'Processed', 'Unprocessed' and 'Unitary'</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "<tr>\n",
    "<td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">translated_processed_pseudogene</td> \n",
    "  <td rowspan=\"2\">Pseudogenes that have mass spec data suggesting that they are also translated. These can be classified into 'Processed', 'Unprocessed'</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">translated_unprocessed_pseudogene</td> \n",
    "  <td rowspan=\"2\">Inactivated immunoglobulin gene. Immunoglobulin gene segments that are inactivated due to frameshift mutations and/or stop codons in the open reading frame</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">tRNA</td> \n",
    "  <td rowspan=\"2\">Transfer RNA</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">unitary_pseudogene</td> \n",
    "  <td rowspan=\"2\">A species specific unprocessed pseudogene without a parent gene, as it has an active orthologue in another species</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">unknown</td> \n",
    "  <td rowspan=\"2\">Entries where the locus type is currently unknown</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "  <td> --- </td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\">unprocessed_pseudogene</td> \n",
    "  <td rowspan=\"2\">Pseudogene that can contain introns since produced by gene duplication</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "<tr>\n",
    "  <td rowspan=\"2\" align=\"center\">vaultRNA</td> \n",
    "  <td rowspan=\"2\" align=\"center\">Short non coding RNA genes that form part of the vault ribonucleoprotein complex</td>\n",
    "  <td>gene</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td>\n",
    "</tr>\n",
    "  <td>transcript</td> \n",
    "  <td>yes</td> \n",
    "  <td>no</td> \n",
    "  <td>no</td> \n",
    "</tr>\n",
    "</table> \n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "**Output:** This script downloads and saves the following data:  \n",
    "- Human Ensembl Gene Set ➞ `Homo_sapiens.GRCh38.<<release>>.gtf`\n",
    "- Human Ensembl-UniProt Identifiers ➞ `Homo_sapiens.GRCh38.<<release>>.uniprot.tsv` \n",
    "- Human Ensembl-Entrez Identifiers ➞ `Homo_sapiens.GRCh38.<<release>>.entrez.tsv` \n",
    "- Human Gene Identifiers ➞ [`Homo_sapiens.gene_info`](ftp://ftp.ncbi.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz), [`hgnc_complete_set.txt`](ftp://ftp.ebi.ac.uk/pub/databases/genenames/new/tsv/hgnc_complete_set.txt)  \n",
    "- Human Protein Identifiers ➞ [`promapping.txt`](https://proconsortium.org/download/current/promapping.txt)  \n",
    "- UniProt Identifiers ➞ [`uniprot_identifier_mapping.tab`](https://www.uniprot.org/uniprot/?query=&fil=organism%3A%22Homo%20sapiens%20(Human)%20%5B9606%5D%22&columns=id%2Cdatabase(GeneID)%2Cdatabase(Ensembl)%2Cdatabase(HGNC)%2Cgenes(PREFERRED)%2Cgenes(ALTERNATIVE))\n",
    "\n",
    "*All Merged Data Sets:*  \n",
    "- `Merged_Human_Ensembl_Entrez_HGNC_Uniprot_Identifiers.txt` \n",
    "- `Merged_gene_rna_protein_identifiers.pkl`  \n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Genomic Typing Dictionary**  \n",
    "Read in the  `genomic_typing_dict.pkl` dictionary, which is needed in order to preprocess the genomic identifier datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://zenodo.org/records/10056198/files/genomic_typing_dict.pkl.zip?download=1'\n",
    "if not os.path.exists(unprocessed_data_location + 'genomic_typing_dict.pkl'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "genomic_type_mapper = pickle.load(open(unprocessed_data_location + 'genomic_typing_dict.pkl', 'rb'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**HGNC Data** \n",
    "\n",
    "_Human Gene Set Data_ - `hgnc_complete_set.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'http://ftp.ebi.ac.uk/pub/databases/genenames/new/tsv/hgnc_complete_set.txt'\n",
    "if not os.path.exists(unprocessed_data_location + 'hgnc_complete_set.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "hgnc = pandas.read_csv(unprocessed_data_location + 'hgnc_complete_set.txt', header=0, delimiter='\\t', low_memory=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preprocess Data_  \n",
    "Data file needs to be lightly cleaned before it can be merged with other data. This light cleaning includes renaming columns, replacing `NaN` with `None`, updating data types (i.e. making all columns type `str`), and unnesting `|` delimited data. The final step is to update the gene_type variable such that each of the variable values is re-grouped to be protein-coding, other or ncRNA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hgnc = hgnc.loc[hgnc['status'].apply(lambda x: x == 'Approved')]\n",
    "hgnc = hgnc[['hgnc_id', 'entrez_id', 'ensembl_gene_id', 'uniprot_ids', 'symbol', 'locus_type', 'alias_symbol', 'name', 'location', 'alias_name']]\n",
    "hgnc.rename(columns={'uniprot_ids': 'uniprot_id', 'location': 'map_location', 'locus_type': 'hgnc_gene_type'}, inplace=True)\n",
    "hgnc['hgnc_id'].replace('.*\\:', '', inplace=True, regex=True)  # strip 'HGNC' off of the identifiers\n",
    "hgnc.fillna('None', inplace=True)  # replace NaN with 'None'\n",
    "hgnc['entrez_id'] = hgnc['entrez_id'].apply(lambda x: str(int(x)) if x != 'None' else 'None')  # make col str\n",
    "\n",
    "# combine certain columns into single column\n",
    "hgnc['name'] = hgnc['name'] + '|' + hgnc['alias_name']\n",
    "hgnc['synonyms'] = hgnc['alias_symbol'] + '|' + hgnc['alias_name'] + '|' + hgnc['name']\n",
    "hgnc['symbol'] = hgnc['symbol'] + '|' + hgnc['alias_symbol']\n",
    "\n",
    "# explode nested data and reformat values in preparation for combining it with other gene identifiers\n",
    "explode_df_hgnc = explodes_data(hgnc.copy(), ['ensembl_gene_id', 'uniprot_id', 'symbol', 'name', 'synonyms'], '|')\n",
    "\n",
    "# reformat hgnc gene type\n",
    "for val in genomic_type_mapper['hgnc_gene_type'].keys():\n",
    "    explode_df_hgnc['hgnc_gene_type'].replace(val, genomic_type_mapper['hgnc_gene_type'][val], inplace=True)\n",
    "\n",
    "# reformat master hgnc gene type\n",
    "explode_df_hgnc['master_gene_type'] = explode_df_hgnc['hgnc_gene_type']\n",
    "master_dict = genomic_type_mapper['hgnc_master_gene_type']\n",
    "for val in master_dict.keys():\n",
    "    explode_df_hgnc['master_gene_type'].replace(val, master_dict[val], inplace=True)\n",
    "\n",
    "# post-process reformatted data\n",
    "explode_df_hgnc.drop(['alias_symbol', 'alias_name'], axis=1, inplace=True)  # remove original gene type column\n",
    "explode_df_hgnc.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "explode_df_hgnc.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Ensembl Data**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Human Gene Set Data_ - `Homo_sapiens.GRCh38.102.gtf.gz`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'ftp://ftp.ensembl.org/pub/release-102/gtf/homo_sapiens/Homo_sapiens.GRCh38.102.gtf.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'Homo_sapiens.GRCh38.102.gtf'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "ensembl_geneset = pandas.read_csv(unprocessed_data_location + 'Homo_sapiens.GRCh38.102.gtf',\n",
    "                                  header = None, delimiter='\\t', skiprows=5, usecols=[8], low_memory=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preprocess Data_  \n",
    "Data file needs to be reformatted in order for it to be able to be merged with the other gene, RNA, and protein identifier data. To do this, we iterate over each row of the data and extract the fields shown below in `column_names`, making each of these extracted fields their own column. The final step is to update the gene_type variable such that each of the variable values is re-grouped to be `protein-coding`, `other` or `ncRNA`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ensembl_data = list(ensembl_geneset[8]); ensembl_df_data = []\n",
    "for i in tqdm(range(0, len(ensembl_data))):\n",
    "    if 'gene_id' in ensembl_data[i] and 'transcript_id' in ensembl_data[i]:\n",
    "        row_dict = {x.split(' \"')[0].lstrip(): x.split(' \"')[1].strip('\"') for x in ensembl_data[i].split(';')[0:-1]}\n",
    "        ensembl_df_data += [(row_dict['gene_id'], row_dict['transcript_id'], row_dict['gene_name'],\n",
    "                           row_dict['gene_biotype'], row_dict['transcript_name'], row_dict['transcript_biotype'])]\n",
    "# convert to data frame\n",
    "ensembl_geneset = pandas.DataFrame(ensembl_df_data,\n",
    "                                   columns=['ensembl_gene_id', 'transcript_stable_id', 'symbol',\n",
    "                                            'ensembl_gene_type', 'transcript_name', 'ensembl_transcript_type'])\n",
    "\n",
    "# reformat ensembl gene type\n",
    "gene_dict = genomic_type_mapper['ensembl_gene_type']\n",
    "for val in gene_dict.keys(): ensembl_geneset['ensembl_gene_type'].replace(val, gene_dict[val], inplace=True)\n",
    "# reformat master gene type\n",
    "ensembl_geneset['master_gene_type'] = ensembl_geneset['ensembl_gene_type']\n",
    "gene_dict = genomic_type_mapper['ensembl_master_gene_type']\n",
    "for val in gene_dict.keys(): ensembl_geneset['master_gene_type'].replace(val, gene_dict[val], inplace=True)\n",
    "# reformat master transcript type\n",
    "ensembl_geneset['ensembl_transcript_type'].replace('vault_RNA', 'vaultRNA', inplace=True, regex=False)\n",
    "ensembl_geneset['master_transcript_type'] = ensembl_geneset['ensembl_transcript_type']\n",
    "trans_dict = genomic_type_mapper['ensembl_master_transcript_type']\n",
    "for val in trans_dict.keys(): ensembl_geneset['master_transcript_type'].replace(val, trans_dict[val], inplace=True)\n",
    "\n",
    "# post-process reformatted data\n",
    "ensembl_geneset.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "ensembl_geneset.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Ensembl Annotation Data**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Ensembl-UniProt_ - `Homo_sapiens.GRCh38.102.uniprot.tsv`  \n",
    "Once the main ensembl gene set has been read in, the next step is to read in the `ensembl-uniprot` mapping file. These files are vital for successfully merging the ensembl identifiers with the uniprot data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url_uniprot = 'ftp://ftp.ensembl.org/pub/release-102/tsv/homo_sapiens/Homo_sapiens.GRCh38.102.uniprot.tsv.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'Homo_sapiens.GRCh38.102.uniprot.tsv'):\n",
    "    data_downloader(url_uniprot, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "ensembl_uniprot = pandas.read_csv(unprocessed_data_location + 'Homo_sapiens.GRCh38.102.uniprot.tsv', header=0, delimiter='\\t', low_memory=False)\n",
    "\n",
    "# preprocess data\n",
    "ensembl_uniprot.rename(columns={'xref': 'uniprot_id', 'gene_stable_id': 'ensembl_gene_id'}, inplace=True)\n",
    "ensembl_uniprot.replace('-', 'None', inplace=True)\n",
    "ensembl_uniprot.fillna('None', inplace=True)\n",
    "ensembl_uniprot = ensembl_uniprot.loc[ensembl_uniprot['xref_identity'].apply(lambda x: x != 'None')]\n",
    "ensembl_uniprot = ensembl_uniprot.loc[ensembl_uniprot['uniprot_id'].apply(lambda x: '-' not in x)]  # remove isoforms\n",
    "ensembl_uniprot = ensembl_uniprot.loc[ensembl_uniprot['info_type'].apply(lambda x: x == 'DIRECT')]\n",
    "# ensembl_uniprot['master_gene_type'] = ['protein-coding'] * len(ensembl_uniprot)\n",
    "# ensembl_uniprot['master_transcript_type'] = ['protein-coding'] * len(ensembl_uniprot)\n",
    "ensembl_uniprot.drop(['db_name', 'info_type', 'source_identity', 'xref_identity', 'linkage_type'], axis=1, inplace=True)\n",
    "ensembl_uniprot.drop_duplicates(subset=None, keep='first', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Ensembl-Entrez_ - `Homo_sapiens.GRCh38.102.entrez.tsv`  \n",
    "Once the main ensembl gene set has been read in, the next step is to read in the `ensembl-entrez` mapping file. These files are vital for successfully merging the ensembl identifiers with the entrez data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url_entrez = 'ftp://ftp.ensembl.org/pub/release-102/tsv/homo_sapiens/Homo_sapiens.GRCh38.102.entrez.tsv.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'Homo_sapiens.GRCh38.102.entrez.tsv'):\n",
    "    data_downloader(url_entrez, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "ensembl_entrez = pandas.read_csv(unprocessed_data_location + 'Homo_sapiens.GRCh38.102.entrez.tsv', header=0, delimiter='\\t', low_memory=False)\n",
    "\n",
    "# preprocess data\n",
    "ensembl_entrez.rename(columns={'xref': 'entrez_id', 'gene_stable_id': 'ensembl_gene_id'}, inplace=True)\n",
    "ensembl_entrez = ensembl_entrez.loc[ensembl_entrez['db_name'].apply(lambda x: x == 'EntrezGene')]\n",
    "ensembl_entrez = ensembl_entrez.loc[ensembl_entrez['info_type'].apply(lambda x: x == 'DEPENDENT')]\n",
    "ensembl_entrez.replace('-', 'None', inplace=True)\n",
    "ensembl_entrez.fillna('None', inplace=True)\n",
    "ensembl_entrez.drop(['db_name', 'info_type', 'source_identity', 'xref_identity', 'linkage_type'], axis=1, inplace=True)\n",
    "ensembl_entrez.drop_duplicates(subset=None, keep='first', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Merge Annotation Data_ - `ensembl_uniprot` + `ensembl_entrez`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merge_cols = list(set(ensembl_entrez).intersection(set(ensembl_uniprot)))\n",
    "ensembl_annot = pandas.merge(ensembl_uniprot, ensembl_entrez, on=merge_cols, how='outer')\n",
    "ensembl_annot.fillna('None', inplace=True)\n",
    "\n",
    "# preview data\n",
    "ensembl_annot.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Merge Ensembl Annotation and Gene Set Data_ - `ensembl_geneset` + `ensembl_annot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merge_cols = list(set(ensembl_annot).intersection(set(ensembl_geneset)))\n",
    "ensembl = pandas.merge(ensembl_geneset, ensembl_annot, on=merge_cols, how='outer')\n",
    "ensembl.fillna('None', inplace=True)\n",
    "ensembl.replace('NA','None', inplace=True, regex=False)\n",
    "\n",
    "# preview data\n",
    "ensembl.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Save Cleaned Ensembl Data_  \n",
    "Save the cleaned Ensembl data so that it can be used when generating node metadata for transcript identifiers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ensembl.to_csv(processed_data_location + 'ensembl_identifier_data_cleaned.txt', header=True, sep='\\t', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**UniProt Data**   \n",
    "_Human Gene Set Data_ - `uniprot_identifier_mapping.tab`\n",
    "\n",
    "This data was obtained by querying the [UniProt Knowledgebase](https://www.uniprot.org/uniprot/) using the *organism:\"Homo sapiens (Human) [9606]\"* keyword and including the following columns:\n",
    "- Entry (Standard)    \n",
    "- GeneID (*Genome Annotation*)  \n",
    "- Ensembl (*Genome Annotation*)  \n",
    "- HGNC (*Organism-specific*)  \n",
    "- Gene names (primary) (*Names & Taxonomy*)    \n",
    "- Gene synonym (primary) (*Names & Taxonomy*)    \n",
    "\n",
    "The URL to access the results of this query is obtained by clicking on the share symbol and copying the free-text from the box. To obtain the data in a tab-delimited format the following string is appended to the end of the URL: \"&format=tab\".\n",
    "\n",
    "**NOTE.** Be sure to obtain a new URL from the [UniProt Knowledgebase](https://www.uniprot.org/uniprot/) when rebuilding to ensure you are getting the most up-to-date data. This query was last generated on `12/02/2020`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://www.uniprot.org/uniprot/?query=&fil=organism%3A%22Homo%20sapiens%20(Human)%20%5B9606%5D%22&columns=id%2Creviewed%2Cdatabase(GeneID)%2Cdatabase(Ensembl)%2Cdatabase(HGNC)%2Cgenes(ALTERNATIVE)%2Cgenes(PREFERRED)&format=tab'\n",
    "if not os.path.exists(unprocessed_data_location + 'uniprot_identifier_mapping.tab'):\n",
    "    data_downloader(url, unprocessed_data_location, 'uniprot_identifier_mapping.tab')\n",
    "\n",
    "# load data\n",
    "uniprot = pandas.read_csv(unprocessed_data_location + 'uniprot_identifier_mapping.tab', header=0, delimiter='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preprocess Data_  \n",
    "Data file needs to be lightly cleaned before it can be merged with other data. This light cleaning includes renaming columns, replacing `NaN` with `None`, and unnesting `\"|\"` delimited data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uniprot.fillna('None', inplace=True)  # replace NaN with 'None'\n",
    "uniprot.rename(columns={'Entry': 'uniprot_id',\n",
    "                        'Cross-reference (GeneID)': 'entrez_id',\n",
    "                        'Ensembl transcript': 'transcript_stable_id',\n",
    "                        'Cross-reference (HGNC)': 'hgnc_id',\n",
    "                        'Gene names  (synonym )': 'synonyms',\n",
    "                        'Gene names  (primary )' :'symbol'}, inplace=True)\n",
    "\n",
    "# update space-delimited synonyms to a pipe (i.e. '|')\n",
    "uniprot['synonyms'] = uniprot['synonyms'].apply(lambda x: '|'.join(x.split()) if x.isupper() else x)\n",
    "\n",
    "# only keep reviewed entries\n",
    "uniprot = uniprot.loc[uniprot['Status'].apply(lambda x: x != 'unreviewed')]\n",
    "\n",
    "# explode nested data\n",
    "explode_df_uniprot = explodes_data(uniprot.copy(), ['transcript_stable_id', 'entrez_id', 'hgnc_id'], ';')\n",
    "explode_df_uniprot = explodes_data(explode_df_uniprot.copy(), ['symbol', 'synonyms'], '|')\n",
    "\n",
    "# strip out uniprot names\n",
    "explode_df_uniprot['transcript_stable_id'].replace('\\s.*','', inplace=True, regex=True)\n",
    "\n",
    "# remove duplicates\n",
    "explode_df_uniprot.drop(['Status'], axis=1, inplace=True)\n",
    "explode_df_uniprot.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "explode_df_uniprot.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**NCBI Data**   \n",
    "_Human Gene Set Data_ - `Homo_sapiens.gene_info`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'ftp://ftp.ncbi.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'Homo_sapiens.gene_info'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "ncbi_gene = pandas.read_csv(unprocessed_data_location + 'Homo_sapiens.gene_info', header=0, delimiter='\\t', low_memory=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preprocess Data_  \n",
    "Data file needs to be lightly cleaned before it can be merged with other data. This light cleaning includes renaming columns, replacing `NaN` with `None`, updating data types (i.e. making all columns type `str`), and unnesting `|` delimited data. Then, the `gene_type` variable is cleaned such that each of the variable's values are re-grouped to be `protein-coding`, `other` or `ncRNA`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocess data\n",
    "ncbi_gene = ncbi_gene.loc[ncbi_gene['#tax_id'].apply(lambda x: x == 9606)]  # remove non-human rows\n",
    "ncbi_gene.replace('-', 'None', inplace=True)\n",
    "ncbi_gene.rename(columns={'GeneID': 'entrez_id', 'Symbol': 'symbol', 'Synonyms': 'synonyms'}, inplace=True)\n",
    "ncbi_gene['synonyms'] = ncbi_gene['synonyms'] + '|' + ncbi_gene['description'] + '|' + ncbi_gene['Full_name_from_nomenclature_authority'] + '|' + ncbi_gene['Other_designations']\n",
    "ncbi_gene['symbol'] = ncbi_gene['Symbol_from_nomenclature_authority'] + '|' + ncbi_gene['symbol']\n",
    "ncbi_gene['name'] = ncbi_gene['Full_name_from_nomenclature_authority'] + '|' + ncbi_gene['description']\n",
    "\n",
    "# explode nested data\n",
    "explode_df_ncbi_gene = explodes_data(ncbi_gene.copy(), ['symbol', 'synonyms', 'name', 'dbXrefs'], '|')\n",
    "\n",
    "# clean up results\n",
    "explode_df_ncbi_gene['entrez_id'] = explode_df_ncbi_gene['entrez_id'].astype(str)\n",
    "explode_df_ncbi_gene = explode_df_ncbi_gene.loc[explode_df_ncbi_gene['dbXrefs'].apply(lambda x: x.split(':')[0] in ['Ensembl', 'HGNC', 'IMGT/GENE-DB'])]\n",
    "explode_df_ncbi_gene['hgnc_id'] = explode_df_ncbi_gene['dbXrefs'].loc[explode_df_ncbi_gene['dbXrefs'].apply(lambda x: x.startswith('HGNC'))]\n",
    "explode_df_ncbi_gene['ensembl_gene_id'] = explode_df_ncbi_gene['dbXrefs'].loc[explode_df_ncbi_gene['dbXrefs'].apply(lambda x: x.startswith('Ensembl'))]\n",
    "explode_df_ncbi_gene.fillna('None', inplace=True)\n",
    "\n",
    "# reformat entrez gene type\n",
    "explode_df_ncbi_gene['entrez_gene_type'] = explode_df_ncbi_gene['type_of_gene']\n",
    "gene_dict = genomic_type_mapper['entrez_gene_type']\n",
    "for val in gene_dict.keys(): explode_df_ncbi_gene['entrez_gene_type'].replace(val, gene_dict[val], inplace=True)\n",
    "# reformat master gene type\n",
    "explode_df_ncbi_gene['master_gene_type'] = explode_df_ncbi_gene['entrez_gene_type']\n",
    "gene_dict = genomic_type_mapper['master_gene_type']\n",
    "for val in gene_dict.keys(): explode_df_ncbi_gene['master_gene_type'].replace(val, gene_dict[val], inplace=True)\n",
    "\n",
    "# post-process reformatted data\n",
    "explode_df_ncbi_gene.drop(['type_of_gene', 'dbXrefs', 'description', 'Nomenclature_status', 'Modification_date',\n",
    "                           'LocusTag', '#tax_id', 'Full_name_from_nomenclature_authority', 'Feature_type',\n",
    "                           'Symbol_from_nomenclature_authority'], axis=1, inplace=True)\n",
    "explode_df_ncbi_gene['hgnc_id'] = explode_df_ncbi_gene['hgnc_id'].replace('HGNC:', '', regex=True)\n",
    "explode_df_ncbi_gene['ensembl_gene_id'] = explode_df_ncbi_gene['ensembl_gene_id'].replace('Ensembl:', '', regex=True)\n",
    "explode_df_ncbi_gene.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "explode_df_ncbi_gene.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Protein Ontology Identifier Mapping Data**   \n",
    "_Protein Ontology Identifier Data_ - `promapping.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://proconsortium.org/download/current/promapping.txt'\n",
    "if not os.path.exists(unprocessed_data_location + 'promapping.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "pro_map = pandas.read_csv(unprocessed_data_location + 'promapping.txt', header=None, names=['pro_id', 'entry', 'pro_mapping'], delimiter='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preprocess Data_  \n",
    "Basic filtering to to include `Protein Ontology` mappings to `Uniprot` identifiers and cleaning to update formatting of accession values (i.e. removing `UniProtKB:`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pro_map = pro_map.loc[pro_map['entry'].apply(lambda x: x.startswith('Uni') and '_VAR' not in x and ', ' not in x)]  # keep 'UniProtKB' rows\n",
    "pro_map = pro_map.loc[pro_map['pro_mapping'].apply(lambda x: x.startswith('exact'))] # keep exact mappings\n",
    "pro_map['pro_id'].replace('PR:','PR_', inplace=True, regex=True)  # replace PR: with PR_\n",
    "pro_map['entry'].replace('(^\\w*\\:)','', inplace=True, regex=True)  # remove id prefixes\n",
    "pro_map = pro_map.loc[pro_map['pro_id'].apply(lambda x: '-' not in x)] # remove isoforms\n",
    "pro_map.rename(columns={'entry': 'uniprot_id'}, inplace=True)  # rename columns before merging\n",
    "pro_map.drop(['pro_mapping'], axis=1, inplace=True)  # remove uneeded columns\n",
    "pro_map.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "pro_map.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### Merging Processed Genomic Identifier Data Sources  \n",
    "Merging all of the genomic identifier data sources is needed in order to create a map that can be used to integrate the different genomic data sources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "*Data Sources:* `hgnc` + `ensembl`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "merge_cols = list(set(explode_df_hgnc.columns).intersection(set(ensembl.columns)))\n",
    "ensembl_hgnc_merged_data = pandas.merge(ensembl, explode_df_hgnc, on=merge_cols, how='outer')\n",
    "\n",
    "# clean up merged data\n",
    "ensembl_hgnc_merged_data.fillna('None', inplace=True)\n",
    "ensembl_hgnc_merged_data.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "ensembl_hgnc_merged_data.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "*Data Sources:* `ensembl_hgnc_merged_data` + `explode_df_uniprot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merge_cols = list(set(ensembl_hgnc_merged_data.columns).intersection(set(explode_df_uniprot.columns)))\n",
    "ensembl_hgnc_uniprot_merged_data = pandas.merge(ensembl_hgnc_merged_data, explode_df_uniprot, on=merge_cols, how='outer')\n",
    "\n",
    "# clean up merged data\n",
    "ensembl_hgnc_uniprot_merged_data.fillna('None', inplace=True)\n",
    "ensembl_hgnc_uniprot_merged_data.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "ensembl_hgnc_uniprot_merged_data.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "*Data Sources:* `ensembl_hgnc_uniprot_merged_data` + `Homo_sapiens.gene_info`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "merge_cols = merge_cols = list(set(ensembl_hgnc_uniprot_merged_data).intersection(set(explode_df_ncbi_gene.columns)))\n",
    "ensembl_hgnc_uniprot_ncbi_merged_data = pandas.merge(ensembl_hgnc_uniprot_merged_data, explode_df_ncbi_gene, on=merge_cols, how='outer')\n",
    "\n",
    "# clean up merged data\n",
    "ensembl_hgnc_uniprot_ncbi_merged_data.fillna('None', inplace=True)\n",
    "ensembl_hgnc_uniprot_ncbi_merged_data.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "ensembl_hgnc_uniprot_ncbi_merged_data.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "*Data Sources:* `ensembl_hgnc_uniprot_ncbi_merged_data` + `promapping.txt`  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "merged_data = pandas.merge(ensembl_hgnc_uniprot_ncbi_merged_data, pro_map, on='uniprot_id', how='outer')\n",
    "\n",
    "# clean up merged data\n",
    "merged_data.fillna('None', inplace=True)\n",
    "merged_data.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "\n",
    "# preview data\n",
    "merged_data.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Fix Symbol Formatting_  \n",
    "some genes are formatted similarly to dates (e.g. `DEC1`), which can be erroneously re-formatted during input as a date value (i.e. `1-DEC`). In order for the data to be successfully merged with other data sources, all date-formatted genes need to be resolved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clean_dates = []\n",
    "for x in tqdm(list(merged_data['symbol'])):\n",
    "    if '-' in x and len(x.split('-')[0]) < 3 and len(x.split('-')[1]) == 3:\n",
    "        clean_dates.append(x.split('-')[1].upper() + x.split('-')[0])\n",
    "    else: clean_dates.append(x)\n",
    "\n",
    "# add cleaned date var back to data set\n",
    "merged_data['symbol'] = clean_dates\n",
    "merged_data.fillna('None', inplace=True)\n",
    "\n",
    "# make sure that all gene and transcript type colunmns have none recoded to unknown or not protein-coding\n",
    "merged_data['hgnc_gene_type'].replace('None', 'unknown', inplace=True, regex=False)\n",
    "merged_data['ensembl_gene_type'].replace('None', 'unknown', inplace=True, regex=False)\n",
    "merged_data['entrez_gene_type'].replace('None', 'unknown', inplace=True, regex=False)\n",
    "merged_data['master_gene_type'].replace('None', 'unknown', inplace=True, regex=False)\n",
    "merged_data['master_transcript_type'].replace('None', 'not protein-coding', inplace=True, regex=False)\n",
    "merged_data['ensembl_transcript_type'].replace('None', 'unknown', inplace=True, regex=False)\n",
    "\n",
    "# remove duplicates\n",
    "merged_data_clean = merged_data.drop_duplicates(subset=None, keep='first')\n",
    "\n",
    "# write data\n",
    "merged_data_clean.to_csv(processed_data_location + 'Merged_Human_Ensembl_Entrez_HGNC_Uniprot_Identifiers.txt', header=True, sep='\\t', index=False)\n",
    "    \n",
    "# preview data\n",
    "merged_data_clean.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Create a Master Mapping Dictionary**  \n",
    "Although the above steps result in a `pandas.Dataframe` of the merged identifiers, there is still work needed in order to be able to obtain a complete mapping between the identifiers. For example, if you were to search for Entrez gene identifier `entrez_259234` you would find the following mappings: `entrez_259234-ENSG00000233316`, `entrez_259234-DSCR10`. If you only had `ENSG00000233316`, with the current data you would be unable to obtain the gene symbol without first mapping to the Entrez gene identifier. \n",
    "\n",
    "To solve this problem, we build a master dictionary where the keys are `ensembl_gene_id`, `transcript_stable_id`, `protein_stable_id`, `uniprot_id`, `entrez_id`, `hgnc_id`, `pro_id`, and `symbol` identifiers and values are the list of genomic identifiers that match to each identifier. It's important to note that there are several labeling identifiers (i.e. `name`, `chromosome`, `map_location`, `Other_designations`, `synonyms`, `transcript_name`, `*_gene_types`, and `trasnscript_type_update`), which will only be mapped when clustered against one of the primary identifier types (i.e. the keys described above).\n",
    "\n",
    "_Note_. The next chunk does a lot of heavy lifting and takes approximately ~40 minutes to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reformat data to convert all nones, empty values, and unknowns to NaN\n",
    "for col in merged_data_clean.columns:\n",
    "    merged_data_clean[col] = merged_data_clean[col].apply(lambda x: '|'.join([i for i in x.split('|') if i != 'None']))\n",
    "merged_data_clean.replace(to_replace=['None', '', 'unknown'], value=numpy.nan, inplace=True)\n",
    "identifiers = [x for x in merged_data_clean.columns if x.endswith('_id')] + ['symbol']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert data to dictionary\n",
    "master_dict = {}\n",
    "for idx in tqdm(identifiers):\n",
    "    grouped_data = merged_data_clean.groupby(idx)\n",
    "    grp_ids = set([x for x in list(grouped_data.groups.keys()) if x != numpy.nan])\n",
    "    for grp in grp_ids:\n",
    "        df = grouped_data.get_group(grp).dropna(axis=1, how='all')\n",
    "        df_cols, key = df.columns, idx + '_' + grp\n",
    "        val_df = [[col + '_' + x for x in set(df[col]) if isinstance(x, str)] for col in df_cols if col != idx]\n",
    "        if len(val_df) > 0:\n",
    "            if key in master_dict.keys(): master_dict[key] += [i for j in val_df for i in j if len(i) > 0]\n",
    "            else: master_dict[key] = [i for j in val_df for i in j if len(i) > 0]  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Finalizing Master Mapping Dictionary_  \n",
    "Then, we need to identify a master gene and transcript type for each entity because the last ran code chunk can result in several genes and transcripts with differing types (i.e. `protein-coding` or `not protein-coding`). The next step collects all information for each gene and transcript and performs a voting procedure to select a single primary gene and transcript type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reformatted_mapped_identifiers = dict()\n",
    "for key, values in tqdm(master_dict.items()):\n",
    "    identifier_info = set(values); gene_prefix = 'master_gene_type_'; trans_prefix = 'master_transcript_type_'\n",
    "    if key.split('_')[0] in ['protein', 'uniprot', 'pro']: pass\n",
    "    elif 'transcript' in key:\n",
    "        trans_match = [x.replace(trans_prefix, '') for x in values if trans_prefix in x]\n",
    "        if len(trans_match) > 0:\n",
    "            t_type_list = ['protein-coding' if ('protein-coding' in trans_match or 'protein_coding' in trans_match) else 'not protein-coding']\n",
    "            identifier_info |= {'transcript_type_update_' + max(set(t_type_list), key=t_type_list.count)}\n",
    "    else:\n",
    "        gene_match = [x.replace(gene_prefix, '') for x in values if x.startswith(gene_prefix) and 'type' in x]\n",
    "        if len(gene_match) > 0:\n",
    "            g_type_list = ['protein-coding' if ('protein-coding' in gene_match or 'protein_coding' in gene_match) else 'not protein-coding']\n",
    "            identifier_info |= {'gene_type_update_' + max(set(g_type_list), key=g_type_list.count)}\n",
    "    reformatted_mapped_identifiers[key] = identifier_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save a copy of the dictionary\n",
    "# output > 4GB requires special approach: https://stackoverflow.com/questions/42653386/does-pickle-randomly-fail-with-oserror-on-large-files\n",
    "filepath = processed_data_location + 'Merged_gene_rna_protein_identifiers.pkl'\n",
    "\n",
    "# defensive way to write pickle.write, allowing for very large files on all platforms\n",
    "max_bytes, bytes_out = 2**31 - 1, pickle.dumps(reformatted_mapped_identifiers)\n",
    "n_bytes = sys.getsizeof(bytes_out)\n",
    "\n",
    "with open(filepath, 'wb') as f_out:\n",
    "    for idx in range(0, n_bytes, max_bytes):\n",
    "        f_out.write(bytes_out[idx:idx+max_bytes])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# # load data\n",
    "# filepath = processed_data_location + 'Merged_gene_rna_protein_identifiers.pkl'\n",
    "\n",
    "# # defensive way to write pickle.load, allowing for very large files on all platforms\n",
    "# max_bytes = 2**31 - 1\n",
    "# input_size = os.path.getsize(filepath)\n",
    "# bytes_in = bytearray(0)\n",
    "\n",
    "# with open(filepath, 'rb') as f_in:\n",
    "#     for _ in range(0, input_size, max_bytes):\n",
    "#         bytes_in += f_in.read(max_bytes)\n",
    "\n",
    "# # load ickled data\n",
    "# reformatted_mapped_identifiers = pickle.loads(bytes_in)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Ensembl Gene-Entrez Gene <a class=\"anchor\" id=\"ensemblgene-entrezgene\"></a>\n",
    "\n",
    "\n",
    "**Purpose:** To map Ensembl gene identifiers to Entrez gene identifiers when creating `gene`-`gene` edges\n",
    "\n",
    "**Output:** `ENSEMBL_GENE_ENTREZ_GENE_MAP.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomic_id_mapper(reformatted_mapped_identifiers,\n",
    "                  processed_data_location + 'ENSEMBL_GENE_ENTREZ_GENE_MAP.txt',\n",
    "                  'ensembl_gene_id', 'entrez_id', 'ensembl_gene_type', 'entrez_gene_type',\n",
    "                  'gene_type_update', 'gene_type_update')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print the number of rows, and preview it\n",
    "egeg_data = pandas.read_csv(processed_data_location + 'ENSEMBL_GENE_ENTREZ_GENE_MAP.txt',\n",
    "                            header=None, delimiter='\\t', low_memory=False,\n",
    "                            names=['Ensembl_Gene_IDs', 'Entrez_Gene_IDs',\n",
    "                                   'Ensembl_Gene_Type', 'Entrez_Gene_Type',\n",
    "                                   'Master_Gene_Type1', 'Master_Gene_Type2'])\n",
    "\n",
    "print('There are {edge_count} ensembl gene-entrez gene edges'.format(edge_count=len(egeg_data)))\n",
    "egeg_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Ensembl Transcript-Protein Ontology <a class=\"anchor\" id=\"ensembltranscript-proteinontology\"></a>\n",
    "\n",
    "**Purpose:** To map Ensembl transcript identifiers to Protein Ontology identifiers when creating `rna`-`protein` edges\n",
    "\n",
    "**Output:** `ENSEMBL_TRANSCRIPT_PROTEIN_ONTOLOGY_MAP.txt`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomic_id_mapper(reformatted_mapped_identifiers,\n",
    "                  processed_data_location + 'ENSEMBL_TRANSCRIPT_PROTEIN_ONTOLOGY_MAP.txt',\n",
    "                  'transcript_stable_id', 'pro_id', 'ensembl_transcript_type', None,\n",
    "                  'transcript_type_update', None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print the number of rows, and preview it\n",
    "etpr_data = pandas.read_csv(processed_data_location + 'ENSEMBL_TRANSCRIPT_PROTEIN_ONTOLOGY_MAP.txt',\n",
    "                            header=None, delimiter='\\t', low_memory=False, usecols=[0, 1, 2, 4],\n",
    "                            names=['Ensembl_Transcript_IDs', 'Protein_Ontology_IDs',\n",
    "                                   'Ensembl_Transcript_Type', 'Master_Transcript_Type'])\n",
    "\n",
    "print('There are {edge_count} ensembl transcript-protein ontology edges'.format(edge_count=len(etpr_data)))\n",
    "etpr_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Entrez Gene-Ensembl Transcript <a class=\"anchor\" id=\"entrezgene-ensembltranscript\"></a>\n",
    "\n",
    "**Purpose:** To map entrez gene identifiers to Ensembl transcript identifiers when creating `gene`-`rna` edges\n",
    "\n",
    "**Output:** `ENTREZ_GENE_ENSEMBL_TRANSCRIPT_MAP.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomic_id_mapper(reformatted_mapped_identifiers,\n",
    "                  processed_data_location + 'ENTREZ_GENE_ENSEMBL_TRANSCRIPT_MAP.txt',\n",
    "                  'entrez_id', 'transcript_stable_id', 'entrez_gene_type', 'ensembl_transcript_type',\n",
    "                  'gene_type_update', 'transcript_type_update')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print the number of rows, and preview it\n",
    "eet_data = pandas.read_csv(processed_data_location + 'ENTREZ_GENE_ENSEMBL_TRANSCRIPT_MAP.txt',\n",
    "                           header=None, delimiter='\\t', low_memory=False,\n",
    "                           names=['Entrez_Gene_IDs', 'Ensembl_Transcript_IDs',\n",
    "                                  'Entrez_Gene_Type', 'Ensembl_Transcript_Type',\n",
    "                                  'Master_Gene_Type', 'Master_Transcript_Type'])\n",
    "\n",
    "print('There are {edge_count} entrez gene identifiers-ensembl transcript edges'.format(edge_count=len(eet_data)))\n",
    "eet_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Entrez Gene-Protein Ontology <a class=\"anchor\" id=\"entrezgene-proteinontology\"></a>\n",
    "\n",
    "**Purpose:** To map Protein Ontology identifiers to Ensembl transcript identifiers when creating the following edges:   \n",
    "- chemical-protein  \n",
    "- gene-protein\n",
    "\n",
    "**Output:** `ENTREZ_GENE_PRO_ONTOLOGY_MAP.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomic_id_mapper(reformatted_mapped_identifiers,\n",
    "                  processed_data_location + 'ENTREZ_GENE_PRO_ONTOLOGY_MAP.txt',\n",
    "                  'entrez_id', 'pro_id', 'entrez_gene_type', None,\n",
    "                  'gene_type_update', None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print the number of rows, and preview it\n",
    "egpr_data = pandas.read_csv(processed_data_location + 'ENTREZ_GENE_PRO_ONTOLOGY_MAP.txt',\n",
    "                            header=None, delimiter='\\t', low_memory=False, usecols = [0, 1, 2, 4],\n",
    "                            names=['Gene_IDs', 'Protein_Ontology_IDs',\n",
    "                                   'Entrez_Gene_Type', 'Master_Gene_Type'])\n",
    "\n",
    "print('There are {edge_count} entrez gene-protein ontology edges'.format(edge_count=len(egpr_data)))\n",
    "egpr_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Gene Symbol-Ensembl Transcript <a class=\"anchor\" id=\"genesymbol-ensembltranscript\"></a>\n",
    "\n",
    "**Purpose:** To map gene symbols to Ensembl transcript identifiers when creating the following edges: \n",
    "- chemical-rna  \n",
    "- rna-anatomy  \n",
    "- rna-cell  \n",
    "\n",
    "**Output:** `GENE_SYMBOL_ENSEMBL_TRANSCRIPT_MAP.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomic_id_mapper(reformatted_mapped_identifiers,\n",
    "                  processed_data_location + 'GENE_SYMBOL_ENSEMBL_TRANSCRIPT_MAP.txt',\n",
    "                  'symbol', 'transcript_stable_id', 'master_gene_type', 'ensembl_transcript_type',\n",
    "                  'gene_type_update', 'transcript_type_update')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print the number of rows, and preview it\n",
    "set_data = pandas.read_csv(processed_data_location + 'GENE_SYMBOL_ENSEMBL_TRANSCRIPT_MAP.txt',\n",
    "                           header=None, delimiter='\\t', low_memory=False,\n",
    "                           names=['Gene_Symbols', 'Ensembl_Transcript_IDs',\n",
    "                                  'Gene_Type', 'Ensembl_Transcript_Type',\n",
    "                                  'Master_Gene_Type', 'Master_Transcript_Type'])\n",
    "\n",
    "print('There are {edge_count} gene symbol-ensembl transcript edges'.format(edge_count=len(set_data.drop_duplicates())))\n",
    "set_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "### STRING-Protein Ontology <a class=\"anchor\" id=\"string-proteinontology\"></a>\n",
    "\n",
    "**Purpose:** To map STRING identifiers to Protein Ontology identifiers when creating `protein`-`protein` edges \n",
    "\n",
    "**Output:** `STRING_PRO_ONTOLOGY_MAP.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomic_id_mapper(reformatted_mapped_identifiers,\n",
    "                  processed_data_location + 'STRING_PRO_ONTOLOGY_MAP.txt',\n",
    "                  'protein_stable_id', 'pro_id', None, None, None, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print the number of rows, and preview it\n",
    "stpr_data = pandas.read_csv(processed_data_location + 'STRING_PRO_ONTOLOGY_MAP.txt',\n",
    "                            header=None, delimiter='\\t', low_memory=False, usecols=[0, 1],\n",
    "                            names=['STRING_IDs', 'Protein_Ontology_IDs'])\n",
    "\n",
    "print('There are {edge_count} string-protein ontology edges'.format(edge_count=len(stpr_data.drop_duplicates())))\n",
    "stpr_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "### Uniprot Accession-Protein Ontology <a class=\"anchor\" id=\"uniprotaccession-proteinontology\"></a>\n",
    "\n",
    "**Purpose:** To map Uniprot accession identifiers to Protein Ontology identifiers when creating the following edges:  \n",
    "- protein-gobp  \n",
    "- protein-gomf  \n",
    "- protein-gocc  \n",
    "- protein-cofactor  \n",
    "- protein-catalyst \n",
    "- protein-pathway\n",
    "\n",
    "**Output:** `UNIPROT_ACCESSION_PRO_ONTOLOGY_MAP.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genomic_id_mapper(reformatted_mapped_identifiers,\n",
    "                  processed_data_location + 'UNIPROT_ACCESSION_PRO_ONTOLOGY_MAP.txt',\n",
    "                  'uniprot_id', 'pro_id', None, None, None, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print the number of rows, and preview it\n",
    "uapr_data = pandas.read_csv(processed_data_location + 'UNIPROT_ACCESSION_PRO_ONTOLOGY_MAP.txt',\n",
    "                            header=None, delimiter='\\t', low_memory=False, usecols=[0, 1],\n",
    "                            names=['Uniprot_Accession_IDs', 'Protein_Ontology_IDs'])\n",
    "\n",
    "print('There are {edge_count} uniprot accession-protein ontology edges'.format(edge_count=len(uapr_data.drop_duplicates())))\n",
    "uapr_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "***\n",
    "***\n",
    "### Other Identifier Mapping <a class=\"anchor\" id=\"other-identifier-mapping\"></a>\n",
    "***\n",
    "* [ChEBI Identifiers](#mesh-chebi)  \n",
    "* [Human Protein Atlas Tissue and Cell Types](#hpa-uberon) \n",
    "* [Human Disease and Phenotype Identifiers](#disease-identifiers) \n",
    "* [Reactome Pathways and the Pathway Ontology](#reactome-pw)  \n",
    "* [Genomic Identifiers and the Sequence Ontology](#genomic-so)  \n",
    "\n",
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ChEBI-MeSH Identifiers <a class=\"anchor\" id=\"mesh-chebi\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [mapping-mesh-to-chebi](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#mapping-mesh-identifiers-to-chebi-identifiers)  \n",
    "\n",
    "**Purpose:** Map MeSH identifiers to ChEBI identifiers when creating the following edges:  \n",
    "- chemical-gene  \n",
    "- chemical-disease\n",
    "\n",
    "**Dependencies:** Recapitulates the [`LOOM`](https://www.bioontology.org/wiki/BioPortal_Mappings) algorithm implemented by BioPortal when creating mappings between resources. The procedure is relatively straightforward and consists of the following:\n",
    "- For all MeSH `SCR Chemicals`, obtain the following information:  \n",
    "  - <u>Identifiers</u>: MeSH identifiers     \n",
    "  - <u>Labels</u>: string labels using the `RDFS:label` object property  \n",
    "  - <u>Synonyms</u>: track down all synonyms using the `vocab:concept` and `vocab:preferredConcept` object properties   \n",
    "- For all ChEBI classes, obtain the following information:  \n",
    "  - <u>Labels</u>: string labels using the `RDFS:label` object property  \n",
    "  - <u>Synonyms</u>: track down all synonyms using all `synonym` object properties \n",
    "  \n",
    "*Alternatively:* You can use the [`ncbo_rest_api.py`](https://gist.github.com/callahantiff/a28fb3160782f42f104e9ec41553af0d) script to pull mappings from the BioPortal API, but note that it takes >2 days for it to finish.\n",
    "\n",
    "**Output:** `CHEBI_MESH_MAP.txt`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***  \n",
    "**MeSH**  \n",
    "Downloads the `nt`-formatted version of the current MeSH vocabulary. Preprocesing is then performed in order to reformat the data so that it can be converted into a Pandas DataFrame in preparation of merging it with `ChEBI` in order to identify overlapping concepts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'ftp://nlmpubs.nlm.nih.gov/online/mesh/rdf/2021/mesh2021.nt'\n",
    "if not os.path.exists(unprocessed_data_location + 'mesh2021.nt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "    \n",
    "# load data\n",
    "mesh = [x.split('> ') for x in tqdm(open(unprocessed_data_location + 'mesh2021.nt').readlines())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocess data\n",
    "mesh_dict, results = {}, []\n",
    "for row in tqdm(mesh):\n",
    "    dbx, lab, msh_type = None, None, None\n",
    "    s, p, o = row[0].split('/')[-1], row[1].split('#')[-1], row[2]  \n",
    "    if s[0] in ['C', 'D'] and ('.' not in s and 'Q' not in s) and len(s) >= 5:\n",
    "        s = 'MESH_' + s\n",
    "        if p == 'preferredConcept' or p == 'concept': dbx = 'MESH_' + o.split('/')[-1]\n",
    "        if 'label' in p.lower(): lab = o.split('\"')[1]\n",
    "        if 'type' in p.lower(): msh_type = o.split('#')[1]\n",
    "        if s in mesh_dict.keys():\n",
    "            if dbx is not None: mesh_dict[s]['dbxref'].add(dbx)\n",
    "            if lab is not None: mesh_dict[s]['label'].add(lab)\n",
    "            if msh_type is not None: mesh_dict[s]['type'].add(msh_type)\n",
    "        else:\n",
    "            mesh_dict[s] = {'dbxref': set() if dbx is None else {dbx},\n",
    "                            'label': set() if lab is None else {lab},\n",
    "                            'type': set() if msh_type is None else {msh_type},\n",
    "                            'synonym': set()}\n",
    "\n",
    "# fine tune dictionary - obtain labels for each entry's synonym identifiers\n",
    "for key in tqdm(mesh_dict.keys()):\n",
    "    for i in mesh_dict[key]['dbxref']:\n",
    "        if len(mesh_dict[key]['dbxref']) > 0 and i in mesh_dict.keys():\n",
    "            mesh_dict[key]['synonym'] |= mesh_dict[i]['label']\n",
    "\n",
    "# expand data and convert to pandas DataFrame\n",
    "for key, value in tqdm(mesh_dict.items()):\n",
    "    results += [[key, list(value['label'])[0], 'NAME']]\n",
    "    if len(value['synonym']) > 0:\n",
    "        for i in value['synonym']:\n",
    "            results += [[key, i, 'SYNONYM']]\n",
    "mesh_filtered = pandas.DataFrame({'CODE': [x[0] for x in results],\n",
    "                                  'TYPE': [x[2] for x in results],\n",
    "                                  'STRING': [x[1] for x in results]})\n",
    "\n",
    "# lowercase all strings and remove white space and punctuation\n",
    "mesh_filtered['STRING'] = mesh_filtered['STRING'].str.lower()\n",
    "mesh_filtered['STRING'] = mesh_filtered['STRING'].str.replace('[^\\w]','')\n",
    "\n",
    "# preview data\n",
    "mesh_filtered.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***  \n",
    "**ChEBI**  \n",
    "Downloads the flat-file containing labels and synonyms for all classes in the `ChEBI` ontology. Preprocessing is then performed in order to reformat the data so that it can be converted into a Pandas DataFrame in preparation of merging it with `MeSH` in order to identify overlapping concepts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'ftp://ftp.ebi.ac.uk/pub/databases/chebi/Flat_file_tab_delimited/names.tsv.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'names.tsv'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "    \n",
    "# load data\n",
    "chebi = pandas.read_csv(unprocessed_data_location + 'names.tsv', header=0, delimiter='\\t')\n",
    "\n",
    "# preprocess data\n",
    "chebi_filtered = chebi[['COMPOUND_ID', 'TYPE', 'NAME']]\n",
    "chebi_filtered.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "chebi_filtered.columns = ['CODE', 'TYPE', 'STRING']\n",
    "\n",
    "# append CHEBI to the number in each code\n",
    "chebi_filtered['CODE'] = chebi_filtered['CODE'].apply(lambda x: \"{}{}\".format('CHEBI_', x))\n",
    "\n",
    "# lowercase all strings and remove white space and punctuation\n",
    "chebi_filtered['STRING'] = chebi_filtered['STRING'].str.lower()\n",
    "chebi_filtered['STRING'] = chebi_filtered['STRING'].str.replace('[^\\w]','')\n",
    "\n",
    "# preview data\n",
    "chebi_filtered.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***  \n",
    "**Merge Identifier Data**  \n",
    "Performs an inner merge of the `MeSH` and `ChEBI` Pandas DataFrames in order to find concepts that exist in both DataFrames. Results are then written out to a text file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# merge data\n",
    "chem_merge = pandas.merge(chebi_filtered[['STRING', 'CODE']], mesh_filtered[['STRING', 'CODE']], on='STRING', how='inner')\n",
    "\n",
    "# filter results\n",
    "mesh_edges = set()\n",
    "for idx, row in chem_merge.drop_duplicates().iterrows():\n",
    "    mesh, chebi = row['CODE_y'], row['CODE_x']\n",
    "    syns = [x for x in mesh_dict[mesh]['dbxref'] if 'C' in x or 'D' in x]\n",
    "    mesh_edges.add(tuple([mesh, chebi]))\n",
    "    if len(syns) > 0:\n",
    "        for x in syns:\n",
    "            mesh_edges.add(tuple([x, chebi]))\n",
    "\n",
    "# write resulting mappings\n",
    "with open(processed_data_location + 'MESH_CHEBI_MAP.txt', 'w') as out:\n",
    "    for pair in mesh_edges:\n",
    "        out.write(pair[0] + '\\t' + pair[1] + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "data = pandas.read_csv(processed_data_location + 'MESH_CHEBI_MAP.txt', header=None, names=['MESH_ID', 'CHEBI_ID'], delimiter='\\t')\n",
    "\n",
    "# preview mapping results\n",
    "print('There are {} MeSH-ChEBI Edges'.format(len(data)))\n",
    "data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "### Disease and Phenotype Identifiers <a class=\"anchor\" id=\"disease-identifiers\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [DisGeNET](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#disgenet)  \n",
    "\n",
    "**Purpose:** This script downloads the Human Phenotype Ontology (HPO), the MonDO Disease Ontology (MONDO), and [disease_mappings.tsv](https://www.disgenet.org/static/disgenet_ap1/files/downloads/disease_mappings.tsv.gz) in order to map UMLS identifiers to HPO and MONDO identifiers when creating the following edges:  \n",
    "- chemical-disease  \n",
    "- disease-phenotype  \n",
    "- chemical-phenotype  \n",
    "- gene-phenotype  \n",
    "- variant-phenotype  \n",
    "\n",
    "**Output:**   \n",
    "- Human Disease Ontology Mappings ➞ `DISEASE_MONDO_MAP.txt`\n",
    "- Human Phenotype Ontology Mappings ➞ `PHENOTYPE_HPO_MAP.txt`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**MONDO Identifiers**  \n",
    "`MONDO` contains DbXRef mappings to other disease terminology identifiers. To make this useful, we will store the DbXRefs as a dictionary with `MONDO` identifiers as the values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download ontology\n",
    "if not os.path.exists(unprocessed_data_location + 'mondo_with_imports.owl'):\n",
    "    command = '{} {} --merge-import-closure -o {}'\n",
    "    os.system(command.format(owltools_location, 'http://purl.obolibrary.org/obo/mondo.owl',\n",
    "                             unprocessed_data_location + 'mondo_with_imports.owl'))\n",
    "    \n",
    "# read data into RDFLib graph object\n",
    "mondo_graph = Graph().parse(unprocessed_data_location + 'mondo_with_imports.owl')\n",
    "print('There are {} axioms in the ontology (date: {})'.format(len(mondo_graph), datetime.datetime.now().strftime('%m/%d/%Y')))\n",
    "\n",
    "# get dbxrefs for all MONDO classes\n",
    "dbxref_res = gets_ontology_class_dbxrefs(mondo_graph)[0]\n",
    "mondo_dict = {str(k).lower().split('/')[-1]: {str(i).split('/')[-1].replace('_', ':') for i in v} for k, v in dbxref_res.items() if 'MONDO' in str(v)}\n",
    "\n",
    "# pickle dictionary\n",
    "pickle.dump(mondo_dict, open(processed_data_location + 'Mondo_Identifier_Map.pkl', 'wb'), protocol=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**HPO Identifiers**  \n",
    "`HPO` contains DbXRef mappings to other disease terminology identifiers. To make this useful, we will store the DbXRefs as a dictionary with `HPO` identifiers as the values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download ontology\n",
    "if not os.path.exists(unprocessed_data_location + 'hp_with_imports.owl'):\n",
    "    command = '{} {} --merge-import-closure -o {}'\n",
    "    os.system(command.format(owltools_location, 'http://purl.obolibrary.org/obo/hp.owl',\n",
    "                             unprocessed_data_location + 'hp_with_imports.owl'))\n",
    "\n",
    "# read data into RDFLib graph object\n",
    "hp_graph = Graph().parse(unprocessed_data_location + 'hp_with_imports.owl')\n",
    "print('There are {} axioms in the ontology (date: {})'.format(len(hp_graph), datetime.datetime.now().strftime('%m/%d/%Y')))\n",
    "\n",
    "# get dbxrefs for all HPO classes\n",
    "dbxref_res = gets_ontology_class_dbxrefs(hp_graph)[0]\n",
    "hp_dict = {str(k).lower().split('/')[-1]: {str(i).split('/')[-1].replace('_', ':') for i in v} for k, v in dbxref_res.items() if 'HP' in str(v)}\n",
    "\n",
    "# pickle dictionary\n",
    "pickle.dump(hp_dict, open(processed_data_location + 'HPO_Identifier_Map.pkl', 'wb'), protocol=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**DisGeNET Disease Mappings**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://www.disgenet.org/static/disgenet_ap1/files/downloads/disease_mappings.tsv.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'disease_mappings.tsv'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "    \n",
    "# load data\n",
    "disease_data = pandas.read_csv(unprocessed_data_location + 'disease_mappings.tsv', header=0, delimiter='\\t')\n",
    "\n",
    "# reformat data\n",
    "disease_data['vocabulary'] = disease_data['vocabulary'].str.lower()\n",
    "disease_data['diseaseId'] = disease_data['diseaseId'].str.lower()\n",
    "disease_data['vocabulary'] = ['doid' if x == 'do' else 'ordoid' if x == 'ordo' else x for x in disease_data['vocabulary']]\n",
    "\n",
    "# preview data\n",
    "disease_data.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Build Disease Identifier Dictionary_  \n",
    "In order to improve efficiency when mapping different disease terminology identifiers to the [MonDO Disease Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#mondo-disease-ontology) and [Human Phenotype Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#human-phenotype-ontology), we create a dictionary of disease identifiers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all CUIs found with HPO and MONDO\n",
    "disease_data_keep = disease_data.query('vocabulary == \"hpo\" | vocabulary == \"mondo\"')\n",
    "\n",
    "# create mondo and hpo dictionary\n",
    "hp_mondo_dict = {}\n",
    "for idx, row in tqdm(disease_data_keep.iterrows(), total=disease_data_keep.shape[0]):\n",
    "    if row['vocabulary'] == 'mondo': key, value = 'umls:' + row['diseaseId'], 'MONDO:' + row['code']\n",
    "    else: key, value = 'umls:' + row['diseaseId'], row['code']\n",
    "    if key in hp_mondo_dict.keys(): hp_mondo_dict[key] |= {value}\n",
    "    else: hp_mondo_dict[key] = {value}\n",
    "# add ontology mappings from MONDO and HPO\n",
    "for key in tqdm(hp_mondo_dict.keys()):\n",
    "    if key in mondo_dict.keys():\n",
    "        hp_mondo_dict[key] = set(list(hp_mondo_dict[key]) + list(mondo_dict[key]))\n",
    "    if key in hp_dict.keys():\n",
    "        hp_mondo_dict[key] = set(list(hp_mondo_dict[key]) + list(hp_dict[key]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get all rows for HPO/MONDO CUIs to obtain mappings to other disease identifiers\n",
    "disease_data_other = disease_data[disease_data.diseaseId.isin(disease_data_keep['diseaseId'])]\n",
    "\n",
    "# get all other codes that map to MONDO or HPO by hopping through MONDO/HPO relevant CUIs\n",
    "disease_dict = {}\n",
    "for idx, row in tqdm(disease_data_other.iterrows(), total=disease_data_other.shape[0]):\n",
    "    if row['vocabulary'] == 'mondo' or row['vocabulary'] == 'hpo':\n",
    "        key, value = 'umls:' + row['diseaseId'].lower(), row['code']\n",
    "        if key in disease_dict.keys(): disease_dict[key] |= {value}\n",
    "        else: disease_dict[key] = {value}\n",
    "    else:\n",
    "        if 'mondo' not in row['code'] or 'hp' not in row['code']:\n",
    "            if ':' not in row['code']: key, value = row['vocabulary'] + ':' + row['code'], hp_mondo_dict['umls:' + row['diseaseId']]\n",
    "            else: key, value = row['code'], hp_mondo_dict['umls:' + row['diseaseId']]\n",
    "            if key in disease_dict.keys(): disease_dict[key] |= value\n",
    "            else: disease_dict[key] = value\n",
    "\n",
    "# add ontology dictionaries\n",
    "disease_dict = {**disease_dict, **mondo_dict, **hp_dict}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Write Mapping Data_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(processed_data_location + 'DISEASE_MONDO_MAP.txt', 'w') as outfile1, open(processed_data_location + 'PHENOTYPE_HPO_MAP.txt', 'w') as outfile2:\n",
    "    for k, v in tqdm(disease_dict.items()):\n",
    "        if any(x for x in v if x.startswith('MONDO')):\n",
    "            for idx in [x.replace(':', '_') for x in v if 'MONDO' in x]:\n",
    "                outfile1.write(k.upper().split(':')[-1] + '\\t' + idx + '\\n')\n",
    "        if any(x for x in v if x.startswith('HP')):\n",
    "            for idx in [x.replace(':', '_') for x in v if 'HP' in x]:\n",
    "                outfile2.write(k.upper().split(':')[-1]  + '\\t' + idx + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preview Processed MONDO Disease Ontology Mappings_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print row count, and preview it\n",
    "dis_data = pandas.read_csv(processed_data_location + 'DISEASE_MONDO_MAP.txt', header=None, names=['Disease_IDs', 'MONDO_IDs'], delimiter='\\t')\n",
    "\n",
    "print('There are {} disease-MONDO edges'.format(len(dis_data)))\n",
    "dis_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preview Processed Human Phenotype Mappings_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print row count, and preview it\n",
    "hp_data = pandas.read_csv(processed_data_location + 'PHENOTYPE_HPO_MAP.txt', header=None, names=['Disease_IDs', 'HP_IDs'], delimiter='\\t')\n",
    "\n",
    "print('There are {} phenotype-HPO edges'.format(len(hp_data)))\n",
    "hp_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "### Human Protein Atlas/GTEx Tissue/Cells - UBERON + Cell Ontology + Cell Line Ontology <a class=\"anchor\" id=\"hpa-uberon\"></a>\n",
    "\n",
    "**Data Source Wiki Page:**  \n",
    "- [human-protein-atlas](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/#human-protein-atlas) \n",
    "- [genotype-tissue-expression-project](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#the-genotype-tissue-expression-gtex-project)  \n",
    "\n",
    "<br>\n",
    "\n",
    "**Purpose:** Downloads a query for cell, tissue, and blood types with overexpressed protein-coding genes in the human proteome ([`proteinatlas_search.tsv`](https://www.proteinatlas.org/api/search_download.php?search=&columns=g,eg,up,pe,rnatsm,rnaclsm,rnacasm,rnabrsm,rnabcsm,rnablsm,scl,t_RNA_adipose_tissue,t_RNA_adrenal_gland,t_RNA_amygdala,t_RNA_appendix,t_RNA_basal_ganglia,t_RNA_bone_marrow,t_RNA_breast,t_RNA_cerebellum,t_RNA_cerebral_cortex,t_RNA_cervix,_uterine,t_RNA_colon,t_RNA_corpus_callosum,t_RNA_ductus_deferens,t_RNA_duodenum,t_RNA_endometrium_1,t_RNA_epididymis,t_RNA_esophagus,t_RNA_fallopian_tube,t_RNA_gallbladder,t_RNA_heart_muscle,t_RNA_hippocampal_formation,t_RNA_hypothalamus,t_RNA_kidney,t_RNA_liver,t_RNA_lung,t_RNA_lymph_node,t_RNA_midbrain,t_RNA_olfactory_region,t_RNA_ovary,t_RNA_pancreas,t_RNA_parathyroid_gland,t_RNA_pituitary_gland,t_RNA_placenta,t_RNA_pons_and_medulla,t_RNA_prostate,t_RNA_rectum,t_RNA_retina,t_RNA_salivary_gland,t_RNA_seminal_vesicle,t_RNA_skeletal_muscle,t_RNA_skin_1,t_RNA_small_intestine,t_RNA_smooth_muscle,t_RNA_spinal_cord,t_RNA_spleen,t_RNA_stomach_1,t_RNA_testis,t_RNA_thalamus,t_RNA_thymus,t_RNA_thyroid_gland,t_RNA_tongue,t_RNA_tonsil,t_RNA_urinary_bladder,t_RNA_vagina,t_RNA_B-cells,t_RNA_dendritic_cells,t_RNA_granulocytes,t_RNA_monocytes,t_RNA_NK-cells,t_RNA_T-cells,t_RNA_total_PBMC,cell_RNA_A-431,cell_RNA_A549,cell_RNA_AF22,cell_RNA_AN3-CA,cell_RNA_ASC_diff,cell_RNA_ASC_TERT1,cell_RNA_BEWO,cell_RNA_BJ,cell_RNA_BJ_hTERT+,cell_RNA_BJ_hTERT+_SV40_Large_T+,cell_RNA_BJ_hTERT+_SV40_Large_T+_RasG12V,cell_RNA_CACO-2,cell_RNA_CAPAN-2,cell_RNA_Daudi,cell_RNA_EFO-21,cell_RNA_fHDF/TERT166,cell_RNA_HaCaT,cell_RNA_HAP1,cell_RNA_HBEC3-KT,cell_RNA_HBF_TERT88,cell_RNA_HDLM-2,cell_RNA_HEK_293,cell_RNA_HEL,cell_RNA_HeLa,cell_RNA_Hep_G2,cell_RNA_HHSteC,cell_RNA_HL-60,cell_RNA_HMC-1,cell_RNA_HSkMC,cell_RNA_hTCEpi,cell_RNA_hTEC/SVTERT24-B,cell_RNA_hTERT-HME1,cell_RNA_HUVEC_TERT2,cell_RNA_K-562,cell_RNA_Karpas-707,cell_RNA_LHCN-M2,cell_RNA_MCF7,cell_RNA_MOLT-4,cell_RNA_NB-4,cell_RNA_NTERA-2,cell_RNA_PC-3,cell_RNA_REH,cell_RNA_RH-30,cell_RNA_RPMI-8226,cell_RNA_RPTEC_TERT1,cell_RNA_RT4,cell_RNA_SCLC-21H,cell_RNA_SH-SY5Y,cell_RNA_SiHa,cell_RNA_SK-BR-3,cell_RNA_SK-MEL-30,cell_RNA_T-47d,cell_RNA_THP-1,cell_RNA_TIME,cell_RNA_U-138_MG,cell_RNA_U-2_OS,cell_RNA_U-2197,cell_RNA_U-251_MG,cell_RNA_U-266/70,cell_RNA_U-266/84,cell_RNA_U-698,cell_RNA_U-87_MG,cell_RNA_U-937,cell_RNA_WM-115,blood_RNA_basophil,blood_RNA_classical_monocyte,blood_RNA_eosinophil,blood_RNA_gdT-cell,blood_RNA_intermediate_monocyte,blood_RNA_MAIT_T-cell,blood_RNA_memory_B-cell,blood_RNA_memory_CD4_T-cell,blood_RNA_memory_CD8_T-cell,blood_RNA_myeloid_DC,blood_RNA_naive_B-cell,blood_RNA_naive_CD4_T-cell,blood_RNA_naive_CD8_T-cell,blood_RNA_neutrophil,blood_RNA_NK-cell,blood_RNA_non-classical_monocyte,blood_RNA_plasmacytoid_DC,blood_RNA_T-reg,blood_RNA_total_PBMC,brain_RNA_amygdala,brain_RNA_basal_ganglia,brain_RNA_cerebellum,brain_RNA_cerebral_cortex,brain_RNA_hippocampal_formation,brain_RNA_hypothalamus,brain_RNA_midbrain,brain_RNA_olfactory_region,brain_RNA_pons_and_medulla,brain_RNA_thalamus&format=tsv)) via [API](https://www.proteinatlas.org/about/help/dataaccess) and median gene-level TPM by tissue for all genes that are not protein-coding ([`GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct`](https://storage.googleapis.com/gtex_analysis_v8/rna_seq_data/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct.gz)) in order to create mappings between cell and tissue type strings to the Uber-Anatomy, Cell Ontology, and Cell Line Ontology concepts (see [human-protein-atlas](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#human-protein-atlas) for details on the mapping process). The mappings are then used to create the following edge types:  \n",
    "- rna-cell line  \n",
    "- rna-tissue type   \n",
    "- protein-cell line  \n",
    "- protein-tissue type  \n",
    "\n",
    "\n",
    "**Output:**  \n",
    "- All HPA tissue and cell type strings ➞ `HPA_tissues.txt`  \n",
    "- Mapping HPA strings to ontology concepts (documentation) ➞ `zooma_tissue_cell_mapping_04JAN2020.xlsx` \n",
    "- Final HPA-ontology mappings ➞ `HPA_GTEx_TISSUE_CELL_MAP.txt`\n",
    "- HPA Edges ➞ `HPA_GTEX_RNA_GENE_PROTEIN_EDGES.txt`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Human Protein Atlas**  \n",
    "To expedite the mapping process, all HPA tissues, cells, cell lines, and fluid types are extracted from the HPA data columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://www.proteinatlas.org/api/search_download.php?search=&columns=g,eg,up,pe,rnatsm,rnaclsm,rnacasm,rnabrsm,rnabcsm,rnablsm,scl,t_RNA_adipose_tissue,t_RNA_adrenal_gland,t_RNA_amygdala,t_RNA_appendix,t_RNA_basal_ganglia,t_RNA_bone_marrow,t_RNA_breast,t_RNA_cerebellum,t_RNA_cerebral_cortex,t_RNA_cervix,_uterine,t_RNA_colon,t_RNA_corpus_callosum,t_RNA_ductus_deferens,t_RNA_duodenum,t_RNA_endometrium_1,t_RNA_epididymis,t_RNA_esophagus,t_RNA_fallopian_tube,t_RNA_gallbladder,t_RNA_heart_muscle,t_RNA_hippocampal_formation,t_RNA_hypothalamus,t_RNA_kidney,t_RNA_liver,t_RNA_lung,t_RNA_lymph_node,t_RNA_midbrain,t_RNA_olfactory_region,t_RNA_ovary,t_RNA_pancreas,t_RNA_parathyroid_gland,t_RNA_pituitary_gland,t_RNA_placenta,t_RNA_pons_and_medulla,t_RNA_prostate,t_RNA_rectum,t_RNA_retina,t_RNA_salivary_gland,t_RNA_seminal_vesicle,t_RNA_skeletal_muscle,t_RNA_skin_1,t_RNA_small_intestine,t_RNA_smooth_muscle,t_RNA_spinal_cord,t_RNA_spleen,t_RNA_stomach_1,t_RNA_testis,t_RNA_thalamus,t_RNA_thymus,t_RNA_thyroid_gland,t_RNA_tongue,t_RNA_tonsil,t_RNA_urinary_bladder,t_RNA_vagina,t_RNA_B-cells,t_RNA_dendritic_cells,t_RNA_granulocytes,t_RNA_monocytes,t_RNA_NK-cells,t_RNA_T-cells,t_RNA_total_PBMC,cell_RNA_A-431,cell_RNA_A549,cell_RNA_AF22,cell_RNA_AN3-CA,cell_RNA_ASC_diff,cell_RNA_ASC_TERT1,cell_RNA_BEWO,cell_RNA_BJ,cell_RNA_BJ_hTERT+,cell_RNA_BJ_hTERT+_SV40_Large_T+,cell_RNA_BJ_hTERT+_SV40_Large_T+_RasG12V,cell_RNA_CACO-2,cell_RNA_CAPAN-2,cell_RNA_Daudi,cell_RNA_EFO-21,cell_RNA_fHDF/TERT166,cell_RNA_HaCaT,cell_RNA_HAP1,cell_RNA_HBEC3-KT,cell_RNA_HBF_TERT88,cell_RNA_HDLM-2,cell_RNA_HEK_293,cell_RNA_HEL,cell_RNA_HeLa,cell_RNA_Hep_G2,cell_RNA_HHSteC,cell_RNA_HL-60,cell_RNA_HMC-1,cell_RNA_HSkMC,cell_RNA_hTCEpi,cell_RNA_hTEC/SVTERT24-B,cell_RNA_hTERT-HME1,cell_RNA_HUVEC_TERT2,cell_RNA_K-562,cell_RNA_Karpas-707,cell_RNA_LHCN-M2,cell_RNA_MCF7,cell_RNA_MOLT-4,cell_RNA_NB-4,cell_RNA_NTERA-2,cell_RNA_PC-3,cell_RNA_REH,cell_RNA_RH-30,cell_RNA_RPMI-8226,cell_RNA_RPTEC_TERT1,cell_RNA_RT4,cell_RNA_SCLC-21H,cell_RNA_SH-SY5Y,cell_RNA_SiHa,cell_RNA_SK-BR-3,cell_RNA_SK-MEL-30,cell_RNA_T-47d,cell_RNA_THP-1,cell_RNA_TIME,cell_RNA_U-138_MG,cell_RNA_U-2_OS,cell_RNA_U-2197,cell_RNA_U-251_MG,cell_RNA_U-266/70,cell_RNA_U-266/84,cell_RNA_U-698,cell_RNA_U-87_MG,cell_RNA_U-937,cell_RNA_WM-115,blood_RNA_basophil,blood_RNA_classical_monocyte,blood_RNA_eosinophil,blood_RNA_gdT-cell,blood_RNA_intermediate_monocyte,blood_RNA_MAIT_T-cell,blood_RNA_memory_B-cell,blood_RNA_memory_CD4_T-cell,blood_RNA_memory_CD8_T-cell,blood_RNA_myeloid_DC,blood_RNA_naive_B-cell,blood_RNA_naive_CD4_T-cell,blood_RNA_naive_CD8_T-cell,blood_RNA_neutrophil,blood_RNA_NK-cell,blood_RNA_non-classical_monocyte,blood_RNA_plasmacytoid_DC,blood_RNA_T-reg,blood_RNA_total_PBMC,brain_RNA_amygdala,brain_RNA_basal_ganglia,brain_RNA_cerebellum,brain_RNA_cerebral_cortex,brain_RNA_hippocampal_formation,brain_RNA_hypothalamus,brain_RNA_midbrain,brain_RNA_olfactory_region,brain_RNA_pons_and_medulla,brain_RNA_thalamus&format=tsv'\n",
    "if not os.path.exists(unprocessed_data_location + 'proteinatlas_search.tsv'):\n",
    "    data_downloader(url, unprocessed_data_location, 'proteinatlas_search.tsv.gz')\n",
    "\n",
    "# load data\n",
    "hpa = pandas.read_csv(unprocessed_data_location + 'proteinatlas_search.tsv', header=0, delimiter='\\t')\n",
    "hpa.fillna('None', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# retrieve terms to map and write results\n",
    "with open(unprocessed_data_location + 'HPA_tissues.txt', 'w') as outfile:\n",
    "    for x in tqdm(list(hpa.columns)):\n",
    "        if x.endswith('[NX]'):\n",
    "            outfile.write(x.split('RNA - ')[-1].split(' [NX]')[:-1][0] + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Genotype-Tissue Expression Project**  \n",
    "Import the tissues, cells, cell lines, and fluids that we externally mapped from HPA and GTEx data to [UBERON](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#uber-anatomy-ontology), the [Cell Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#cell-ontology), and the [Cell Line Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#cell-line-ontology)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "url='https://storage.googleapis.com/gtex_analysis_v8/rna_seq_data/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "gtex = pandas.read_csv(unprocessed_data_location + 'GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct', header=0, skiprows=2, delimiter='\\t')\n",
    "gtex.fillna('None', inplace=True)  # replace NaN with 'None'\n",
    "gtex['Name'].replace('(\\..*)','', inplace=True, regex=True)  # remove identifier type, which appears after '.'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url='https://zenodo.org/records/10056198/files/zooma_tissue_cell_mapping_04JAN2020.xlsx.zip?download=1'\n",
    "if not os.path.exists(unprocessed_data_location + 'zooma_tissue_cell_mapping_04JAN2020.xlsx'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "    \n",
    "# load ontology mapping data\n",
    "mapping_data = pandas.read_excel(open(unprocessed_data_location + 'zooma_tissue_cell_mapping_04JAN2020.xlsx', 'rb'),\n",
    "                                 sheet_name='Concept_Mapping - 04JAN2020', header=0, engine='openpyxl')\n",
    "mapping_data.fillna('None', inplace=True)  # convert NaN to None\n",
    "\n",
    "# preview data\n",
    "mapping_data.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Write HPA and GTEx Mapping Data_  \n",
    "The HPA and GTEx mapping data is written locally so that it can be used by the `PheKnowLator` algorithm when creating the knowledge graph edge lists. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(processed_data_location + 'HPA_GTEx_TISSUE_CELL_MAP.txt', 'w') as out:\n",
    "    for idx, row in tqdm(mapping_data.iterrows(), total=mapping_data.shape[0]):\n",
    "        if row['UBERON'] != 'None': out.write(str(row['TERM']).strip() + '\\t' + str(row['UBERON']).strip() + '\\n')\n",
    "        if row['CL'] != 'None': out.write(str(row['TERM']).strip() + '\\t' + str(row['CL']).strip() + '\\n')\n",
    "        if row['CLO'] != 'None': out.write(str(row['TERM']).strip() + '\\t' + str(row['CLO']).strip() + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load mapping data\n",
    "mapping_data = pandas.read_csv(processed_data_location + 'HPA_GTEx_TISSUE_CELL_MAP.txt', header=None, names=['TISSUE_CELL_TERM', 'ONTOLOGY_IDs'], delimiter='\\t')\n",
    "\n",
    "# preview data\n",
    "mapping_data.head(n=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Create Edge Data Set**\n",
    "\n",
    "_Human Protein Atlas_  \n",
    "The `HPA` data is looped over and reformatted such that all tissue, cell, cell lines, and fluid types are stored as a nested list. The anatomy type is specified as an item in the list according to its type in order to make mapping more efficient while building the knowledge graph edge list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hpa_results = []\n",
    "for idx, row in tqdm(hpa.iterrows(), total=hpa.shape[0]):\n",
    "    ens, gene, uniprot, evid = str(row['Ensembl']), str(row['Gene']), str(row['Uniprot']), str(row['Evidence'])\n",
    "    if row['RNA tissue specific NX'] != 'None':\n",
    "        for x in row['RNA tissue specific NX'].split(';'):\n",
    "            hpa_results += [[ens, gene, uniprot, evid, 'anatomy', str(x.split(':')[0])]]\n",
    "    if row['RNA cell line specific NX'] != 'None':\n",
    "        for x in row['RNA cell line specific NX'].split(';'):\n",
    "            hpa_results += [[ens, gene, uniprot, evid, 'cell line', str(x.split(':')[0])]]\n",
    "    if row['RNA brain regional specific NX'] != 'None':\n",
    "        for x in row['RNA brain regional specific NX'].split(';'):\n",
    "            hpa_results += [[ens, gene, uniprot, evid, 'anatomy', str(x.split(':')[0])]]\n",
    "    if row['RNA blood cell specific NX'] != 'None':\n",
    "        for x in row['RNA blood cell specific NX'].split(';'):\n",
    "            hpa_results += [[ens, gene, uniprot, evid, 'anatomy', str(x.split(':')[0])]]\n",
    "    if row['RNA blood lineage specific NX'] != 'None':\n",
    "        for x in row['RNA blood lineage specific NX'].split(';'):\n",
    "            hpa_results += [[ens, gene, uniprot, evid, 'anatomy', str(x.split(':')[0])]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Genotype-Tissue Expression Project_  \n",
    "The `GTEx` edge data is created by first filtering out all _protein-coding_ genes that appear in the `HPA` cell transcriptome data set. Once filter so that we are only left noncoding genes, we perform an additional filtering step to only add genes and their corresponding tissue, cell, or fluid, if the median expression is `>= 1.0`. The `GTEx` is formatted such that all anatomical entities occur as their own column and all unique genes occur as a row, thus the expression filtering step is performed while also reformatting the file. The genes and tissues/cells/fluids that meet criteria are stored as a nested list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove rows that contain protein coding genes already in the hpa data\n",
    "hpa_genes = list(hpa['Ensembl'].drop_duplicates(keep='first', inplace=False))\n",
    "gtex = gtex.loc[gtex['Name'].apply(lambda x: x not in hpa_genes)]\n",
    "\n",
    "# loop over data and re-organize - only keep results with tpm >= 1 and if gene symbol is not a protein-coding gene\n",
    "gtex_results = []\n",
    "for idx, row in tqdm(gtex.iterrows(), total=gtex.shape[0]):\n",
    "    for col in list(gtex.columns)[2:]:\n",
    "        typ = 'cell line' if 'Cells' in col else 'anatomy'\n",
    "        if row[col] >= 1.0:\n",
    "            evidence = 'Evidence at transcript level'\n",
    "            gtex_results += [[str(row['Name']), str(row['Description']), 'None', evidence, typ, str(col)]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Writes Edge Data*  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(processed_data_location + 'HPA_GTEX_RNA_GENE_PROTEIN_EDGES.txt', 'w') as out:\n",
    "    for x in tqdm(hpa_results + gtex_results):\n",
    "        out.write(x[0] + '\\t' + x[1] + '\\t' + x[2] + '\\t' + x[3] + '\\t' + x[4] + '\\t' + x[5] + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, return edge count, and preview it\n",
    "hpa_edges = pandas.read_csv(processed_data_location + 'HPA_GTEX_RNA_GENE_PROTEIN_EDGES.txt',\n",
    "                           header=None, low_memory=False, sep='\\t',\n",
    "                           names=['Ensembl_IDs', 'Gene_Symbols', 'Uniprot_IDs', 'Evidence', 'Anatomy_Type', 'Anatomy'])\n",
    "\n",
    "print('There are {edge_count} edges'.format(edge_count=len(hpa_edges)))\n",
    "hpa_edges.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "### Mapping Reactome Pathways to the Pathway Ontology <a class=\"anchor\" id=\"reactome-pw\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [Pathway Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/#pathway-ontology)  \n",
    "\n",
    "**Purpose:** This script downloads the [canonical pathways](http://compath.scai.fraunhofer.de/export_mappings) and [kegg-reactome pathway mappings](https://github.com/ComPath/resources/blob/master/mappings/kegg_reactome.csv) files from the [ComPath Ecosystem](https://github.com/ComPath) in order to create the following identifier mappings:  \n",
    "- `Reactome Pathway Identifiers`  ➞ `KEGG Pathway Identifiers` ➞ `Pathway Ontology Identifiers` \n",
    "\n",
    "**Output:**  \n",
    "- `REACTOME_PW_GO_MAPPINGS.txt`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Pathway Ontology**   \n",
    "Use [OWL Tools](https://github.com/owlcollab/owltools/wiki) to download the [Pathway Ontology](http://www.obofoundry.org/ontology/pw.html). Once downloaded, we read the ontology in as a `RDFLib` graph object so that we can query it to obtain all `DbXRefs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download ontology\n",
    "if not os.path.exists(unprocessed_data_location + 'pw_with_imports.owl'):\n",
    "    command = '{} {} --merge-import-closure -o {}'\n",
    "    os.system(command.format(owltools_location, 'http://purl.obolibrary.org/obo/pw.owl',\n",
    "                             unprocessed_data_location + 'pw_with_imports.owl'))\n",
    "\n",
    "# load the knowledge graph\n",
    "pw_graph = Graph().parse(unprocessed_data_location + 'pw_with_imports.owl')\n",
    "print('There are {} axioms in the ontology (date: {})'.format(len(pw_graph), datetime.datetime.now().strftime('%m/%d/%Y')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Reformat Mapping Results_  \n",
    "Create a dictionary of mapping results where pathway ontology identifiers are values and the keys are `DbXRef` identifiers.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get dbxref results\n",
    "dbxref_res = gets_ontology_class_dbxrefs(pw_graph)[0]\n",
    "dbxref_dict = {str(k).lower().split('/')[-1]: {str(i).split('/')[-1].replace('_', ':') for i in v} for k, v in dbxref_res.items() if 'PW_' in str(v)}\n",
    "\n",
    "# get synonym results\n",
    "syn_res = gets_ontology_class_synonyms(pw_graph)[0]\n",
    "synonym_dict = {str(k).lower().split('/')[-1]: {str(i).split('/')[-1].replace('_', ':') for i in v} for k, v in syn_res.items() if 'PW_' in str(v)}\n",
    "\n",
    "# combine results into single dictionary\n",
    "id_mappings = {**dbxref_dict, **synonym_dict}\n",
    "\n",
    "print('There are {} results (date: {})'.format(len(id_mappings), datetime.datetime.now().strftime('%m/%d/%Y')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Reactome Pathways**  \n",
    "Download a file of all [Reactome Pathways](https://reactome.org/download/current/ReactomePathways.txt), [Reactome's GO Annotations]('https://reactome.org/download/current/gene_association.reactome.gz'), and [Reactome's mappings to CHEBI](https://reactome.org/download/current/ChEBI2Reactome_All_Levels.txt). This file will be filtered to only include human pathways."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Reactome Pathway Stable Identifiers_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://reactome.org/download/current/ReactomePathways.txt'\n",
    "if not os.path.exists(unprocessed_data_location + 'ReactomePathways.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "reactome_pathways = pandas.read_csv(unprocessed_data_location + 'ReactomePathways.txt', header=None, delimiter='\\t', low_memory=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove all non-human pathways and save as list\n",
    "reactome_pathways = reactome_pathways.loc[reactome_pathways[2].apply(lambda x: x == 'Homo sapiens')] \n",
    "reactome_map = {x:set(['PW_0000001']) for x in set(list(reactome_pathways[0]))}     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Reactome's Mappings to GO Annotations_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://reactome.org/download/current/gene_association.reactome.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'gene_association.reactome'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "reactome_pathways2 = pandas.read_csv(unprocessed_data_location + 'gene_association.reactome', header=None, delimiter='\\t', skiprows=3, low_memory=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove all non-human pathways and save as list\n",
    "reactome_pathways2 = reactome_pathways2.loc[reactome_pathways2[12].apply(lambda x: x == 'taxon:9606')] \n",
    "reactome_map.update({x.split(':')[-1]:set(['PW_0000001']) for x in set(list(reactome_pathways2[5]))})     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Reactome's Mappings to ChEBI_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://reactome.org/download/current/ChEBI2Reactome_All_Levels.txt'\n",
    "if not os.path.exists(unprocessed_data_location + 'ChEBI2Reactome_All_Levels.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "reactome_pathways3 = pandas.read_csv(unprocessed_data_location + 'ChEBI2Reactome_All_Levels.txt', header=None, delimiter='\\t', low_memory=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove all non-human pathways and save as list\n",
    "reactome_pathways3 = reactome_pathways3.loc[reactome_pathways3[5].apply(lambda x: x == 'Homo sapiens')] \n",
    "reactome_map.update({x:set(['PW_0000001']) for x in set(list(reactome_pathways3[1]))})     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**ComPath Reactome Pathway Mappings**  \n",
    "Use [ComPath Mappings](https://github.com/ComPath/resources/tree/master/mappings) to obtain the following mappings:  `Reactome Pathways`  ➞ `KEGG Pathways` ➞ `Pathway Ontology` "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Canonical Pathways_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url1 = 'http://compath.scai.fraunhofer.de/export_mappings'\n",
    "if not os.path.exists(unprocessed_data_location + 'compath_canonical_pathway_mappings.txt'):\n",
    "    data_downloader(url1, unprocessed_data_location, 'compath_canonical_pathway_mappings.txt')\n",
    "\n",
    "# load data\n",
    "compath_cannonical = pandas.read_csv(unprocessed_data_location + 'compath_canonical_pathway_mappings.txt', header=None, delimiter='\\t', low_memory=False)\n",
    "compath_cannonical.fillna('None', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx, row in tqdm(compath_cannonical.iterrows(), total=compath_cannonical.shape[0]):\n",
    "    if row[6] == 'kegg' and 'kegg:' + row[5].strip('path:hsa') in id_mappings.keys() and row[2] == 'reactome':\n",
    "        for x in id_mappings['kegg:' + row[5].strip('path:hsa')]:\n",
    "            if row[1] in reactome_map.keys(): reactome_map[row[1]] |= set([x.split('/')[-1]])\n",
    "            else: reactome_map[row[1]] = set([x.split('/')[-1]])\n",
    "    if (row[2] == 'kegg' and 'kegg:' + row[1].strip('path:hsa') in id_mappings.keys()) and row[6] == 'reactome':\n",
    "        for x in id_mappings['kegg:' + row[1].strip('path:hsa')]:\n",
    "            if row[5] in reactome_map.keys(): reactome_map[row[5]] |= set([x.split('/')[-1]])\n",
    "            else: reactome_map[row[5]] = set([x.split('/')[-1]])         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_KEGG - Reactome Mappings_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url2 = 'https://raw.githubusercontent.com/ComPath/resources/master/mappings/kegg_reactome.csv'\n",
    "if not os.path.exists(unprocessed_data_location + 'kegg_reactome.csv'):\n",
    "    data_downloader(url2, unprocessed_data_location, 'kegg_reactome.csv')\n",
    "\n",
    "# load data\n",
    "kegg_reactome_map = pandas.read_csv(unprocessed_data_location + 'kegg_reactome.csv', header=0, delimiter=',', low_memory=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx, row in tqdm(kegg_reactome_map.iterrows(), total=kegg_reactome_map.shape[0]):\n",
    "    if row['Source Resource'] == 'reactome' and 'kegg:' + row['Target ID'].strip('path:hsa') in id_mappings.keys():\n",
    "        for x in id_mappings['kegg:' + row['Target ID'].strip('path:hsa')]:\n",
    "            if row['Source ID'] in reactome_map.keys(): reactome_map[row['Source ID']] |= set([x.split('/')[-1]])\n",
    "            else: reactome_map[row['Source ID']] = set([x.split('/')[-1]])\n",
    "    if row['Target Resource'] == 'reactome' and 'kegg:' + row['Source Resource'].strip('path:hsa') in id_mappings.keys():\n",
    "        for x in id_mappings['kegg:' + row['Source ID'].strip('path:hsa')]:\n",
    "            if row['Target ID'] in reactome_map.keys(): reactome_map[row['Target ID']] |= set([x.split('/')[-1]])\n",
    "            else: reactome_map[row['Target ID']] = set([x.split('/')[-1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Reactome Pathway GO Annotation Mappings**  \n",
    "Use Reactome's [API](https://reactome.org/dev/content-service) to obtain the following mappings: `Reactome Pathway Identifiers`  ➞ `Gene Ontology Identifiers`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "for request_ids in tqdm(list(chunks(list(reactome_map.keys()), 20))):\n",
    "    result, key = content.query_ids(ids=','.join(request_ids)), 'goBiologicalProcess'\n",
    "    if result is not None and (isinstance(result, List) or result['code'] != 404):\n",
    "        for res in result:\n",
    "            if key in res.keys():\n",
    "                if res['stId'] in reactome_map.keys(): reactome_map[res['stId']] |= {'GO_' + res[key]['accession']}\n",
    "                else: reactome_map[res['stId']] = {'GO_' + res[key]['accession']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Write Data*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reformat identifiers -- replacing ontology concepts with ':' to '_'\n",
    "temp_dict = dict()\n",
    "for key, value in tqdm(reactome_map.items()):\n",
    "    temp_dict[key] = set(x.replace(':', '_') for x in value)\n",
    "\n",
    "# overwrite original reactome dict with cleaned mappings\n",
    "reactome_map = temp_dict\n",
    "\n",
    "# output data\n",
    "with open(processed_data_location + 'REACTOME_PW_GO_MAPPINGS.txt', 'w') as out:\n",
    "    for key in tqdm(reactome_map.keys()):\n",
    "        for x in reactome_map[key]:\n",
    "            if x.startswith('PW') or x.startswith('GO'): out.write(key + '\\t' + x + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print row count, and preview it\n",
    "pw_data = pandas.read_csv(processed_data_location + 'REACTOME_PW_GO_MAPPINGS.txt', header=None, names=['Pathway_IDs', 'Mapping_IDs'], delimiter='\\t')\n",
    "\n",
    "print('There are {edge_count} pathway ontology mappings'.format(edge_count=len(pw_data)))\n",
    "pw_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "***\n",
    "\n",
    "### Mapping Genomic Identifiers to the Sequence Ontology <a class=\"anchor\" id=\"genomic-soo\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [Sequence Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/_edit#sequence-ontology)  \n",
    "\n",
    "**Purpose:** This script downloads the `genomic_sequence_ontology_mappings.xlsx` file in order to create the following identifier mappings:  \n",
    "- `Gene BioTypes`  ➞ `Sequence Ontology Identifiers`  \n",
    "- `RNA BioTypes`  ➞ `Sequence Ontology Identifiers`  \n",
    "- `variant Types`  ➞ `Sequence Ontology Identifiers`\n",
    "\n",
    "**Output:**  \n",
    "- `SO_GENE_TRANSCRIPT_VARIANT_TYPE_MAPPING.txt`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url='https://zenodo.org/records/10056198/files/genomic_sequence_ontology_mappings.xlsx.zip?download=1'\n",
    "if not os.path.exists(unprocessed_data_location + 'genomic_sequence_ontology_mappings.xlsx'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "mapping_data = pandas.read_excel(open(unprocessed_data_location + 'genomic_sequence_ontology_mappings.xlsx', 'rb'),\n",
    "                                 sheet_name='GenomicType_SO_Map_09Mar2020', header=0, engine='openpyxl')\n",
    "\n",
    "# convert data to dictionary\n",
    "genomic_type_so_map = {}\n",
    "for idx, row in tqdm(mapping_data.iterrows(), total=mapping_data.shape[0]):\n",
    "    genomic_type_so_map[row['source_*_type'] + '_' + row['Genomic']] = row['SO ID']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Genes**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in genomic mapping data\n",
    "genomic_mapped_ids = pickle.load(open(processed_data_location + 'Merged_gene_rna_protein_identifiers.pkl', 'rb'))\n",
    "\n",
    "sequence_map = {}\n",
    "for identifier in tqdm(genomic_mapped_ids.keys()):    \n",
    "    if identifier.startswith('entrez_id_') and identifier.replace('entrez_id_', '') != 'None':\n",
    "        id_clean = identifier.replace('entrez_id_', '')\n",
    "        \n",
    "        # get identifier types\n",
    "        ensembl = [x.replace('ensembl_gene_type_', '') for x in genomic_mapped_ids[identifier] if x.startswith('ensembl_gene_type') and x != 'ensembl_gene_type_unknown']\n",
    "        hgnc = [x.replace('hgnc_gene_type_', '')  for x in genomic_mapped_ids[identifier] if x.startswith('hgnc_gene_type') and x != 'hgnc_gene_type_unknown']\n",
    "        entrez = [x.replace('entrez_gene_type_', '')  for x in genomic_mapped_ids[identifier] if x.startswith('entrez_gene_type') and x != 'entrez_gene_type_unknown']\n",
    "        \n",
    "        # determine gene type\n",
    "        if len(ensembl) > 0: gene_type = genomic_type_so_map[ensembl[0].replace('ensembl_gene_type_', '') + '_Gene']\n",
    "        elif len(hgnc) > 0: gene_type = genomic_type_so_map[hgnc[0].replace('hgnc_gene_type_', '') + '_Gene']\n",
    "        elif len(entrez) > 0: gene_type = genomic_type_so_map[entrez[0].replace('entrez_gene_type_', '') + '_Gene']\n",
    "        else: gene_type = 'SO_0000704'  \n",
    "        \n",
    "        # update sequence map\n",
    "        if id_clean in sequence_map.keys(): sequence_map[id_clean] += [gene_type]\n",
    "        else: sequence_map[id_clean] = [gene_type]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Transcripts**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in processed Ensembl Transcript data \n",
    "transcript_data = pandas.read_csv(processed_data_location + 'ensembl_identifier_data_cleaned.txt', header=0, delimiter='\\t', low_memory=False)\n",
    "\n",
    "# convert to dictionary\n",
    "transcripts = {}\n",
    "for idx, row in tqdm(transcript_data.iterrows(), total=transcript_data.shape[0]):\n",
    "    if row['transcript_stable_id'] != 'None':\n",
    "        if row['transcript_stable_id'].replace('transcript_stable_id_', '') in transcripts.keys():\n",
    "            transcripts[row['transcript_stable_id'].replace('transcript_stable_id_', '')] += [row['ensembl_transcript_type']]\n",
    "        else: transcripts[row['transcript_stable_id'].replace('transcript_stable_id_', '')] = [row['ensembl_transcript_type']]\n",
    "            \n",
    "# update so map dictionary\n",
    "for identifier in tqdm(transcripts.keys()):\n",
    "    if transcripts[identifier][0] == 'protein_coding': trans_type = genomic_type_so_map['protein-coding_Transcript']\n",
    "    elif transcripts[identifier][0] == 'misc_RNA': trans_type = genomic_type_so_map['miscRNA_Transcript']\n",
    "    else: trans_type = genomic_type_so_map[list(set(transcripts[identifier]))[0] + '_Transcript']\n",
    "    sequence_map[identifier] = [trans_type, 'SO_0000673']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Variants**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# read in variant summary data \n",
    "url = 'ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'variant_summary.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "    \n",
    "# load data    \n",
    "variant_data = pandas.read_csv(unprocessed_data_location + 'variant_summary.txt', header=0, delimiter='\\t', low_memory=False)\n",
    "\n",
    "# convert to dictionary\n",
    "variants = {}\n",
    "for idx, row in tqdm(variant_data.iterrows(), total=variant_data.shape[0]):\n",
    "    if row['Assembly'] == 'GRCh38' and row['RS# (dbSNP)'] != -1:\n",
    "        if 'rs' + str(row['RS# (dbSNP)']) in variants.keys(): variants['rs' + str(row['RS# (dbSNP)'])] |= set([row['Type']])\n",
    "        else: variants['rs' + str(row['RS# (dbSNP)'])] = set([row['Type']])\n",
    "\n",
    "# update so map dictionary\n",
    "for identifier in tqdm(variants.keys()):\n",
    "    for typ in variants[identifier]:\n",
    "        var_type = genomic_type_so_map[typ.lower() + '_Variant']\n",
    "        if identifier in sequence_map.keys(): sequence_map[identifier] += [var_type]\n",
    "        else: sequence_map[identifier] = [var_type]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*** \n",
    "**Write Data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reformat data and write it out\n",
    "with open(processed_data_location + 'SO_GENE_TRANSCRIPT_VARIANT_TYPE_MAPPING.txt', 'w') as outfile:\n",
    "    for key in tqdm(sequence_map.keys()):\n",
    "        for map_type in sequence_map[key]:\n",
    "            outfile.write(key + '\\t' + map_type + '\\n')\n",
    "\n",
    "# load data, print row count, and preview it\n",
    "so_data = pandas.read_csv(processed_data_location + 'SO_GENE_TRANSCRIPT_VARIANT_TYPE_MAPPING.txt', header=None, delimiter='\\t', names=['Identifier', 'Sequence_Ontology_ID'])\n",
    "\n",
    "print('There are {edge_count} sequence ontology mappings'.format(edge_count=len(so_data)))\n",
    "so_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Combine Pathway and Sequence Ontology Mapping Data in Dictionary**  \n",
    "Combine the pathway and sequence mapping data into a dictionary and output it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine genomic and pathway maps\n",
    "subclass_mapping = {}  \n",
    "sequence_map.update(reactome_map)\n",
    "\n",
    "# iterate over pathway lists and combine them\n",
    "for key in tqdm(sequence_map.keys()):\n",
    "    subclass_mapping[key] = sequence_map[key]\n",
    "\n",
    "# save a copy of the dictionary\n",
    "pickle.dump(subclass_mapping, open(construction_approach_location + 'subclass_construction_map.pkl', 'wb'), protocol=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "***\n",
    "***\n",
    "### CREATE EDGE DATASETS  <a class=\"anchor\" id=\"create-edge-datasets\"></a>\n",
    "***\n",
    "***\n",
    "\n",
    "### Ontologies  <a class=\"anchor\" id=\"ontologies\"></a>\n",
    "***\n",
    "- [Protein Ontology](#protein-ontology)  \n",
    "- [Relations Ontology](#relations-ontology)  \n",
    "\n",
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### Protein Ontology <a class=\"anchor\" id=\"protein-ontology\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [protein-ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#human-phenotype-ontology)  \n",
    "\n",
    "**Purpose:** This script uses [OWLTools](https://github.com/owlcollab/owltools) to download the [pr.owl](http://purl.obolibrary.org/obo/pr.owl) (with imports) file from [ProConsortium.org](https://proconsortium.org/) in order to create a version of the ontology that contains only human proteins. This is achieved by performing forward and reverse breadth first search over all proteins which are `owl:subClassOf` [Homo sapiens protein](https://proconsortium.org/app/entry/PR%3A000029067/).\n",
    "\n",
    "<br>\n",
    "\n",
    "**Output:**  \n",
    "- Human Protein Ontology ➞ `human_pro.owl`\n",
    "- Classified Human Protein Ontology (Hermit) ➞ `human_pro_closed.owl`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download ontology\n",
    "if not os.path.exists(unprocessed_data_location + 'pr_with_imports.owl'):\n",
    "    command = '{} {} --merge-import-closure -o {}'\n",
    "    os.system(command.format(owltools_location, 'http://purl.obolibrary.org/obo/pr.owl',\n",
    "                             unprocessed_data_location + 'pr_with_imports.owl'))\n",
    "    \n",
    "# read in ontology as graph (the ontology is large so this takes ~60 minutes)\n",
    "print('Loading Protein Ontology')\n",
    "pr_graph = Graph().parse(unprocessed_data_location + 'pr_with_imports.owl')\n",
    "print('There are {} axioms in the ontology (date: {})'.format(len(pr_graph), datetime.datetime.now().strftime('%m/%d/%Y')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Convert Ontology to Directed MulitGraph_  \n",
    "In order to create a version of the ontology which includes all relevant human edges, we need to first convert the KG to a [directed multigraph](https://networkx.github.io/documentation/stable/reference/classes/multidigraph.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "networkx_mdg: networkx.MultiDiGraph = networkx.MultiDiGraph()\n",
    "    \n",
    "for s, p, o in tqdm(pr_graph):\n",
    "    networkx_mdg.add_edge(s, o, **{'key': p})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Identify Human Proteins_   \n",
    "A list of human proteins is obtained by querying the ontology to return all ontology classes `only_in_taxon some Homo sapiens`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Approach 1 - Query Loaded Graph to Obtain Human Protein Classes*  \n",
    "Does not require using external resources or SPARQL Endpoints. This is the preferred approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "human_classes_restriction = list(pr_graph.triples((None, OWL.someValuesFrom, obo.NCBITaxon_9606)))\n",
    "human_classes = [list(pr_graph.subjects(RDFS.subClassOf, x[0])) for x in human_classes_restriction]\n",
    "human_pro_classes = list(str(i) for j in human_classes for i in j if 'PR_' in str(i))\n",
    "\n",
    "print('There are {} edges in the ontology (date:{})'.format(len(human_pro_classes), datetime.datetime.now().strftime('%m/%d/%Y')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Approach 2 - Query PRO Consortium SPARQL Endpoint to Obtain Human Protein Classes*  \n",
    "This approach should only be used when the PRO endpoint is not limiting the number of results that are returned. As of `October 2021`, this was happening so please use *Approach 1* which is guaranteed to return the correct results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# # download data\n",
    "# url = 'https://sparql.proconsortium.org/virtuoso/sparql?query=PREFIX+obo%3A+%3Chttp%3A%2F%2Fpurl.obolibrary.org%2Fobo%2F%3E%0D%0A%0D%0ASELECT+%3FPRO_term%0D%0AFROM+%3Chttp%3A%2F%2Fpurl.obolibrary.org%2Fobo%2Fpr%3E%0D%0AWHERE+%7B%0D%0A+++++++%3FPRO_term+rdf%3Atype+owl%3AClass+.%0D%0A+++++++%3FPRO_term+rdfs%3AsubClassOf+%3Frestriction+.%0D%0A+++++++%3Frestriction+owl%3AonProperty+obo%3ARO_0002160+.%0D%0A+++++++%3Frestriction+owl%3AsomeValuesFrom+obo%3ANCBITaxon_9606+.%0D%0A%0D%0A+++++++%23+use+this+to+filter-out+things+like+hgnc+ids%0D%0A+++++++FILTER+%28regex%28%3FPRO_term%2C%22http%3A%2F%2Fpurl.obolibrary.org%2Fobo%2F*%22%29%29+.%0D%0A%7D&format=text%2Fhtml&debug='\n",
    "# if not os.path.exists(unprocessed_data_location + 'human_pro_classes.html'):\n",
    "#     data_downloader(url, unprocessed_data_location, 'human_pro_classes.html')\n",
    "\n",
    "# # load data\n",
    "# df_list = pandas.read_html(unprocessed_data_location + 'human_pro_classes.html')\n",
    "\n",
    "# # extract data from html table - pro classes only_in_taxon some Homo sapiens\n",
    "# human_pro_classes = list(df_list[-1]['PRO_term'])\n",
    "# print('There are {} edges in the ontology (date:{})'.format(len(human_pro_classes), datetime.datetime.now().strftime('%m/%d/%Y')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Construct Human PRO_   \n",
    "Now that we have all of the paths from the original graph that are relevant to humans, we can construct a human-only version of the PRotein Ontology. After building the human subset, we verify the number of connected components and get 1. However, after reformatting the graph using [OWLTools](https://github.com/owlcollab/owltools) you will see that there are 3 connected components: component 1 (n=`1051673`); component 2 (n=`12`); and component 3 (n=`2`). The contents of components 2 and 3 are shown below:\n",
    "\n",
    "```python\n",
    "[{'http://purl.obolibrary.org/obo/IAO_0000115',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#hasAlternativeId',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#hasBroadSynonym',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#hasDbXref',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#hasExactSynonym',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#hasNarrowSynonym',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#hasOBONamespace',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#hasRelatedSynonym',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#id',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#is_transitive',\n",
    "  'http://www.geneontology.org/formats/oboInOwl#shorthand',\n",
    "  'http://www.w3.org/2002/07/owl#AnnotationProperty'},\n",
    " \n",
    " {'N41f0be4cf00c48929605b1e69a09f326',\n",
    "  'http://www.w3.org/2002/07/owl#Ontology'}]\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new graph using bfs paths\n",
    "human_pro_graph = Graph()\n",
    "human_networkx_mdg = networkx.MultiDiGraph()\n",
    "\n",
    "for node in tqdm(human_pro_classes):\n",
    "    forward = list(networkx.edge_bfs(networkx_mdg, URIRef(node), orientation='original'))\n",
    "    reverse = list(networkx.edge_bfs(networkx_mdg, URIRef(node), orientation='reverse'))\n",
    "    \n",
    "    # add edges from forward and reverse bfs paths\n",
    "    for path in set(forward + reverse):\n",
    "        human_pro_graph.add((path[0], path[2], path[1]))\n",
    "        human_networkx_mdg.add_edge(path[0], path[1], **{'key': path[2]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get connected component information\n",
    "print('Finding Connected Components')\n",
    "components = list(networkx.connected_components(human_networkx_mdg.to_undirected()))\n",
    "component_dict = sorted(components, key=len, reverse=True)\n",
    "\n",
    "# if more than 1 connected component, only keep the biggest\n",
    "if len(component_dict) > 1:\n",
    "    print('Cleaning Graph: Removing Small Disconnected Components')\n",
    "    for node in tqdm([x for y in component_dict[1:] for x in list(y)]):\n",
    "        human_pro_graph.remove((node, None, None))\n",
    "\n",
    "# save data\n",
    "print('Saving Human Subset of the Protein Ontology')\n",
    "human_pro_graph.serialize(destination=unprocessed_data_location + 'human_pro.owl', format='xml')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Classify Ontology_  \n",
    "To ensure that we have correctly built the new ontology, we run the hermit reasoner over it to ensure that there are no incomplete triples or inconsistent classes. In order to do this, we will call the reasoner using [OWLTools](https://github.com/owlcollab/owltools), which this script assumes has already been downloaded to the `./resources/lib` directory. The following arguments are then called to run the reasoner (from the command line):  \n",
    "\n",
    "___\n",
    "\n",
    "```bash\n",
    "../pkt_kg/libs/owltools ./resources/processed_data/unprocessed_data/human_pro.owl --reasoner elk --run-reasoner --assert-implied -o ./resources/processed_data/human_pro_closed.owl\n",
    "```\n",
    "___\n",
    "\n",
    "\n",
    "_**Note.** This step takes around 5 minutes to run. When run from the command line the reasoner determined that the ontology was consistent and 200 new axioms were inferred (12/01/2020)._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run reasoner\n",
    "command = '{} {} --reasoner {} --run-reasoner --assert-implied -o {}'\n",
    "os.system(command.format(owltools_location, unprocessed_data_location + 'human_pro.owl', 'elk',\n",
    "                         ontology_data_location + 'pr_with_imports.owl'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Examine Cleaned Human PRO_  \n",
    "Once we have cleaned the ontology we can get counts of components, nodes, and edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gets_ontology_statistics(ontology_data_location + 'pr_with_imports.owl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "***\n",
    "\n",
    "### Relations Ontology <a class=\"anchor\" id=\"relations-ontology\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [RO](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#relation-ontology)  \n",
    "\n",
    "**Purpose:** This script downloads the [ro.owl](http://purl.obolibrary.org/obo/ro.owl) file from [obofoundry.org](http://www.obofoundry.org/) in order to obtain all `ObjectProperties` and their inverse relations.  \n",
    "\n",
    "**Output:** \n",
    "- Relations and Inverse Relations ➞ `INVERSE_RELATIONS.txt`\n",
    "- Relations and Labels ➞ `RELATIONS_LABELS.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download ontology\n",
    "if not os.path.exists(unprocessed_data_location + 'ro_with_imports.owl'):\n",
    "    command = '{} {} --merge-import-closure -o {}'\n",
    "    os.system(command.format(owltools_location, 'http://purl.obolibrary.org/obo/ro.owl',\n",
    "                             unprocessed_data_location + 'ro_with_imports.owl'))\n",
    "# load graph\n",
    "ro_graph = Graph().parse(unprocessed_data_location + 'ro_with_imports.owl')\n",
    "print('There are {} edges in the ontology (date:{})'.format(len(ro_graph), datetime.datetime.now().strftime('%m/%d/%Y')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Identify Relations and Inverse Relations**  \n",
    "Identify all relations and their inverse relations using the `owl:inverseOf` property. To make it easier to look up the inverse relations, each pair is listed twice, for example:  \n",
    "- [location of](http://www.ontobee.org/ontology/RO?iri=http://purl.obolibrary.org/obo/RO_0001015) `owl:inverseOf` [located in](http://www.ontobee.org/ontology/RO?iri=http://purl.obolibrary.org/obo/RO_0001025)  \n",
    "- [located in](http://www.ontobee.org/ontology/RO?iri=http://purl.obolibrary.org/obo/RO_0001025) `owl:inverseOf` [location of](http://www.ontobee.org/ontology/RO?iri=http://purl.obolibrary.org/obo/RO_0001015)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(relations_data_location + 'INVERSE_RELATIONS.txt', 'w') as outfile:\n",
    "    outfile.write('Relation' + '\\t' + 'Inverse_Relation' + '\\n')\n",
    "    for s, p, o in tqdm(ro_graph):\n",
    "        if 'owl#inverseOf' in str(p):\n",
    "            if 'RO' in str(s) and 'RO' in str(o):\n",
    "                outfile.write(str(s.split('/')[-1]) + '\\t' + str(o.split('/')[-1]) + '\\n')\n",
    "                outfile.write(str(o.split('/')[-1]) + '\\t' + str(s.split('/')[-1]) + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preview Processed Data_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print row count, and preview it\n",
    "ro_data = pandas.read_csv(relations_data_location + 'INVERSE_RELATIONS.txt', header=0, delimiter='\\t')\n",
    "\n",
    "print('There are {edge_count} RO Relations and Inverse Relations'.format(edge_count=len(ro_data)))\n",
    "ro_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Get Relations Labels**  \n",
    "Identify all relations and their labels for use when building the knowledge graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {str(x[2]).lower(): str(x[0]) for x in ro_graph if '/RO_' in str(x[0]) and 'label' in str(x[1]).lower()}\n",
    "\n",
    "# write data to file\n",
    "with open(relations_data_location + 'RELATIONS_LABELS.txt', 'w') as outfile:\n",
    "    outfile.write('Label' + '\\t' + 'Relation' + '\\n')\n",
    "    for k, v in results.items():\n",
    "        outfile.write(str(v).split('/')[-1] + '\\t' + str(k) + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preview Processed Data_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print row count, and preview it\n",
    "ro_data_label = pandas.read_csv(relations_data_location + 'RELATIONS_LABELS.txt', header=0, delimiter='\\t')\n",
    "\n",
    "print('There are {edge_count} RO Relations and Labels'.format(edge_count=len(ro_data_label)))\n",
    "ro_data_label.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "***\n",
    "***\n",
    "### Linked Data <a class=\"anchor\" id=\"linked-data\"></a>\n",
    "***\n",
    "* [Clinvar Variant-Diseases and Phenotypes](#clinvar-variant) \n",
    "* [Uniprot Protein-Cofactor and Protein-Catalyst](#uniprot-protein-cofactorcatalyst)  \n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "***\n",
    "### Clinvar Variant-Diseases and Phenotypes <a class=\"anchor\" id=\"clinvar-variant\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [Clinvar](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#clinvar)  \n",
    "\n",
    "**Purpose:** This script downloads the [variant_summary.txt](ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz) file from [ClinVar](https://www.ncbi.nlm.nih.gov/clinvar/) in order to create the following edges:  \n",
    "- gene-variant  \n",
    "- variant-disease  \n",
    "- variant-phenotype  \n",
    "\n",
    "**Output:** `CLINVAR_VARIANT_GENE_DISEASE_PHENOTYPE_EDGES.txt`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'variant_summary.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "clinvar_data = pandas.read_csv(unprocessed_data_location + 'variant_summary.txt', header=0, delimiter='\\t', low_memory=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Preprocess Data_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace NaN with 'None'\n",
    "clinvar_data.fillna('None', inplace=True)\n",
    "\n",
    "# explode nested data\n",
    "explode_df_clinvar = explodes_data(clinvar_data.copy(), ['PhenotypeIDS'], ';')\n",
    "explode_df_clinvar = explodes_data(explode_df_clinvar.copy(), ['PhenotypeIDS'], ',')\n",
    "\n",
    "# edit column formatting\n",
    "explode_df_clinvar['PhenotypeIDS'].replace('Orphanet:ORPHA','ORPHA:', inplace=True, regex=True)\n",
    "explode_df_clinvar['PhenotypeIDS'].replace('Human Phenotype Ontology:HP:','HP_', inplace=True, regex=True)\n",
    "\n",
    "# write data\n",
    "explode_df_clinvar.to_csv(processed_data_location + 'CLINVAR_VARIANT_GENE_DISEASE_PHENOTYPE_EDGES.txt', header=True, sep='\\t', encoding='utf-8', index=False)\n",
    "\n",
    "# print row count and preview data\n",
    "print('There are {edge_count} variant edges'.format(edge_count=len(explode_df_clinvar)))\n",
    "explode_df_clinvar.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<br>\n",
    "\n",
    "***\n",
    "\n",
    "### Uniprot  Protein-Cofactor and Protein-Catalyst <a class=\"anchor\" id=\"uniprot-protein-cofactorcatalyst\"></a>\n",
    "\n",
    "**Data Source Wiki Page:** [Uniprot](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/#uniprot-knowledgebase)  \n",
    "\n",
    "**Purpose:** This script downloads the [uniprot-cofactor-catalyst.tab](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources/#uniprot-knowledgebase) file from the [Uniprot Knowledge Base](https://www.uniprot.org) in order to create the following edges:  \n",
    "- protein-cofactor  \n",
    "- protein-catalyst  \n",
    "\n",
    "**Data:** This data was obtained by querying the [UniProt Knowledgebase](https://www.uniprot.org/uniprot/) using the *reviewed:yes AND organism:\"Homo sapiens (Human) [9606]\"\"* keyword and including the following columns:\n",
    "- Entry (Standard) \n",
    "- Status (Standard) \n",
    "- PRO (*Miscellaneous*)  \n",
    "- ChEBI (Cofactor) (*Chemical entities*)   \n",
    "- ChEBI (Catalytic activity) (*Chemical entities*)  \n",
    "\n",
    "The URL to access the results of this query is obtained by clicking on the share symbol and copying the free-text from the box. To obtain the data in a tab-delimited format the following string is appended to the end of the URL: \"&format=tab\".\n",
    "\n",
    "**NOTE.** Be sure to obtain a new URL from the [UniProt Knowledgebase](https://www.uniprot.org/uniprot/) when rebuilding to ensure you are getting the most up-to-date data. This query was last generated on `12/02/2020`.\n",
    "\n",
    "<br>\n",
    "\n",
    "**Output:**  \n",
    "- protein-cofactor ➞ `UNIPROT_PROTEIN_COFACTOR.txt`\n",
    "- protein-catalyst ➞ `UNIPROT_PROTEIN_CATALYST.txt`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'https://www.uniprot.org/uniprot/?query=&fil=organism%3A%22Homo%20sapiens%20(Human)%20%5B9606%5D%22&columns=id%2Creviewed%2Centry%20name%2Cdatabase(PRO)%2Cchebi(Cofactor)%2Cchebi(Catalytic%20activity)&format=tab'\n",
    "if not os.path.exists(unprocessed_data_location + 'uniprot-cofactor-catalyst.tab'):\n",
    "    data_downloader(url, unprocessed_data_location, 'uniprot-cofactor-catalyst.tab')\n",
    "\n",
    "# upload datta\n",
    "data = open(unprocessed_data_location + 'uniprot-cofactor-catalyst.tab').readlines()\n",
    "\n",
    "# reformat data and write it out\n",
    "with open(processed_data_location + 'UNIPROT_PROTEIN_COFACTOR.txt', 'w') as outfile1, open(processed_data_location + 'UNIPROT_PROTEIN_CATALYST.txt', 'w') as outfile2:\n",
    "    for line in tqdm(data):\n",
    "        # get cofactors\n",
    "        if 'CHEBI' in line.split('\\t')[4]: \n",
    "            for i in line.split('\\t')[4].split(';'):\n",
    "                chebi = i.split('[')[-1].replace(']', '').replace(':', '_')\n",
    "                outfile1.write('PR_' + line.split('\\t')[3].strip(';') + '\\t' + chebi + '\\n')\n",
    "        # get catalysts\n",
    "        if 'CHEBI' in line.split('\\t')[5]:       \n",
    "            for i in line.split('\\t')[5].split(';'):\n",
    "                chebi = i.split('[')[-1].replace(']', '').replace(':', '_')\n",
    "                outfile2.write('PR_' + line.split('\\t')[3].strip(';') + '\\t' + chebi + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Cofactor Data**  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print row count, and preview it\n",
    "pcp1_data = pandas.read_csv(processed_data_location + 'UNIPROT_PROTEIN_COFACTOR.txt', header=None, names=['Protein_Ontology_IDs', 'CHEBI_IDs'], delimiter='\\t')\n",
    "\n",
    "print('There are {edge_count} protein-cofactor edges'.format(edge_count=len(pcp1_data)))\n",
    "pcp1_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "\n",
    "**Catalyst Data**  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data, print row count, and preview it\n",
    "pcp2_data = pandas.read_csv(processed_data_location + 'UNIPROT_PROTEIN_CATALYST.txt', header=None, names=['Protein_Ontology_IDs', 'CHEBI_IDs'], delimiter='\\t')\n",
    "\n",
    "print('There are {edge_count} protein-catalyst edges'.format(edge_count=len(pcp2_data)))\n",
    "pcp2_data.head(n=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "***\n",
    "***\n",
    "### INSTANCE AND/OR SUBCLASS (NON-ONTOLOGY CLASS) METADATA <a class=\"anchor\" id=\"create-instance-metadata\"></a>\n",
    "***\n",
    "\n",
    "**Data Source Wiki Page:** [Dependencies](https://github.com/callahantiff/PheKnowLator/wiki/Dependencies/#node-metadata) \n",
    "\n",
    "**Purpose:** The goal of this section is to obtain metadata for each non-ontology instance and/or subclass data source and all relations used in the knowledge graph. For **[`Release V2.0.0`](https://github.com/callahantiff/PheKnowLator/wiki/v2.0.0)**, the following are non-ontology instance and/or subclass data and require the compiling of metadata:\n",
    "- [Genes](#gene-metadata)\n",
    "- [RNA](#rna-metadata)\n",
    "- [Variants](#variant-metadata)  \n",
    "- [Pathways](#pathway-metadata)\n",
    "- [Relations](#relations-metadata)\n",
    "\n",
    "<br>\n",
    "\n",
    "**Metadata:** The <u>metadata</u> we will gather includes:  \n",
    "\n",
    "| **Metadata Type** | **Definition** | **Example Node**  | **Example Node Metadata** | \n",
    "| :---: | :---: | :---: | :---: | \n",
    "| Label | The primary label or name for the node | `R-HSA-1006173` | \"CFH:Host cell surface\" |       \n",
    "| Description | A definition or other useful details about the node | `rs794727058` | This `germline` `single nucleotide variant` located on chromosome `5 (GRCh38: NC_000005.10, start/stop positions (126555930/126555930))` with `pathogenic` clinical significance and a last review date of `2/23/2015` (review status: `criteria provided, single submitter`). |        \n",
    "| Synonym | Alternative terms used for a node | `81399` | \"OR1-1, OR7-21\" |           \n",
    "\n",
    "The metadata information will be used to create the following edges in the knowledge graph:  \n",
    "- **Label** ➞ node `rdfs:label`  \n",
    "- **Description** ➞ node `obo:IAO_0000115` description \n",
    "- **Synonyms** ➞ node `oboInOwl:hasExactSynonym` synonym \n",
    "\n",
    "<br>\n",
    "\n",
    "*<b>NOTE.</b> All node metadata are written to the `node_data` directory as a `pickled` dictionary called `node_metadata_dict.pkl`. The algorithm will look for this dictionary in the `node_data` directory and if it is not there, then no node metadata will be created.*\n",
    "\n",
    "<br>\n",
    "\n",
    "### Prepare Metadata Dictionaries\n",
    "***\n",
    "\n",
    "**Purpose:** To create the resources needed in order to create metadata dictionaries, which are in turn used to obtain metadata for instance and/or subclass data nodes. This process has the following steps:\n",
    "\n",
    "**1. [Generate Metadata Dictionaries](#generate-metadata-dictionaries):** In order to efficiently obtain metadata for all non-ontology instance and/or subclass data nodes and all relations, we first read in the data for each type (i.e. genes, rna, pathways, variants, and relations) and convert them into a dictionary. Then, each metadata dictionary is merged together and saved to a `master_metadata_dictionary`, keyed by identifier.\n",
    "  - <u>Input Datasets</u>:  \n",
    "    - Genes ➞ `Homo_sapiens.gene_info`   \n",
    "    - RNA ➞ `ensembl_identifier_data_cleaned.txt` \n",
    "    - Pathways ➞ [`reactome2py API`](https://github.com/reactome/reactome2py) ; `ReactomePathways.txt`; `gene_association.reactome.gz`; `ChEBI2Reactome_All_Levels.txt`; `kegg_reactome.csv`   \n",
    "    - Variants ➞ `variant_summary.txt`  \n",
    "    - Relations ➞ `ro_with_imports.owl`  \n",
    "    \n",
    "Example Metadata Dictionary Output:\n",
    "\n",
    "```python\n",
    "{\n",
    "    'nodes': {\n",
    "        'http://www.ncbi.nlm.nih.gov/gene/1': {\n",
    "            'Label': 'A1BG',\n",
    "            'Description': \"A1BG has locus group protein-coding' and is located on chromosome 19 (19q13.43).\",\n",
    "            'Synonym': 'HYST2477alpha-1B-glycoprotein|HEL-S-163pA|ABG|A1B|GAB'} ... },\n",
    "    'relations': {\n",
    "        'http://purl.obolibrary.org/obo/RO_0002533': {\n",
    "            'Label': 'sequence atomic unit',\n",
    "            'Description': 'Any individual unit of a collection of like units arranged in a linear order',\n",
    "            'Synonym': 'None'} ... }\n",
    "}\n",
    "```\n",
    "\n",
    "<br>\n",
    "\n",
    "**2. [Write Metadata Files](#write-metadata-files):** The `master_metadata_dictionary` dictionary from _Step 1_ is `pickled` and saved to the `resources/node_data/` directory.\n",
    "\n",
    "<br>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Metadata Dictionaries  <a class=\"anchor\" id=\"generate-metadata-dictionaries\"></a>\n",
    "In this step, the goal is to create a metadata dictionary for each node type that does not rely on API data. In this case, only the **Gene**, **RNA**, and **Variant** nodes require data that is not from an API.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### Genes Metadata Dictionary <a class=\"anchor\" id=\"gene-metadata\"></a>\n",
    "\n",
    "The nested dictionary of gene metadata is created by looping over the merged data described in the prior column. The `keys` of the dictionary are `Entrez gene identifiers` and the `values` are dictionaries for each metadata type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# entrez gene data\n",
    "entrez_gene_data = pandas.read_csv(unprocessed_data_location + 'Homo_sapiens.gene_info', header=0, delimiter='\\t', low_memory=False)\n",
    "\n",
    "# remove all rows that are not human\n",
    "entrez_gene_data = entrez_gene_data.loc[entrez_gene_data['#tax_id'].apply(lambda x: x == 9606)]\n",
    "\n",
    "# replace NaN and '-' with 'None'\n",
    "entrez_gene_data.fillna('None', inplace=True)\n",
    "entrez_gene_data.replace('-','None', inplace=True, regex=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create metadata\n",
    "genes, lab, desc, syn = [], [], [], []\n",
    "for idx, row in tqdm(entrez_gene_data.iterrows(), total=entrez_gene_data.shape[0]):\n",
    "    gene_id, sym, defn, gene_type = row['GeneID'], row['Symbol'], row['description'], row['type_of_gene']\n",
    "    chrom, map_loc, s1, s2 = row['chromosome'], row['map_location'], row['Synonyms'], row['Other_designations']\n",
    "    if gene_id != 'None':\n",
    "        genes.append('http://www.ncbi.nlm.nih.gov/gene/' + str(gene_id))\n",
    "        if sym != 'None' or sym != '': lab.append(sym)\n",
    "        else: lab.append('Entrez_ID:' + gene_id)\n",
    "        if 'None' not in [defn, gene_type, chrom, map_loc]:\n",
    "            desc_str = \"{} has locus group '{}' and is located on chromosome {} ({}).\"\n",
    "            desc.append(desc_str.format(sym, gene_type, chrom, map_loc))\n",
    "        else: desc.append(\"{} locus group '{}'.\".format(sym, gene_type))\n",
    "        if s1 != 'None' and s2 != 'None': syn.append('|'.join(set([x for x in (s1 + s2).split('|') if x != 'None' or x != ''])))\n",
    "        elif s1 != 'None': syn.append('|'.join(set([x for x in s1.split('|') if x != 'None' or x != ''])))\n",
    "        elif s2 != 'None': syn.append('|'.join(set([x for x in s2.split('|') if x != 'None' or x != ''])))\n",
    "        else: syn.append('None')\n",
    "\n",
    "# combine into new data frame\n",
    "metadata = pandas.DataFrame(list(zip(genes, lab, desc, syn)), columns=['ID', 'Label', 'Description', 'Synonym'])\n",
    "metadata = metadata.astype(str)\n",
    "metadata.drop_duplicates(subset='ID', keep='first', inplace=True)\n",
    "\n",
    "# convert df to dictionary\n",
    "metadata.set_index('ID', inplace=True)\n",
    "gene_metadata_dict = metadata.to_dict('index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### RNA Metadata Dictionary  <a class=\"anchor\" id=\"rna-metadata\"></a>\n",
    "\n",
    "The nested dictionary of rna metadata is created by looping over the cleaned human [Ensembl](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#ensembl) gene, RNA, and protein identifier data set (`ensembl_identifier_data_cleaned.txt`). The `keys` of the dictionary are `Ensembl transcript identifiers` and the `values` are dictionaries for each metadata type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "rna_gene_data = pandas.read_csv(processed_data_location + 'ensembl_identifier_data_cleaned.txt', header=0, delimiter='\\t', low_memory=False)\n",
    "\n",
    "# remove rows without identifiers\n",
    "rna_gene_data = rna_gene_data.loc[rna_gene_data['transcript_stable_id'].apply(lambda x: x != 'None')]\n",
    "\n",
    "# remove unneeded columns\n",
    "rna_gene_data.drop(['ensembl_gene_id', 'symbol', 'protein_stable_id', 'uniprot_id', 'master_transcript_type',\n",
    "                    'entrez_id', 'ensembl_gene_type', 'master_gene_type', 'symbol'], axis=1, inplace=True)\n",
    "\n",
    "# remove duplicates\n",
    "rna_gene_data.drop_duplicates(subset=['transcript_stable_id', 'transcript_name', 'ensembl_transcript_type'], keep='first', inplace=True)\n",
    "\n",
    "# replace NaN with 'None'\n",
    "rna_gene_data.fillna('None', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create metadata\n",
    "rna, lab, desc, syn = [], [], [], []\n",
    "for idx, row in tqdm(rna_gene_data.iterrows(), total=rna_gene_data.shape[0]):\n",
    "    rna_id, ent_type, nme = row['transcript_stable_id'], row['ensembl_transcript_type'], row['transcript_name']\n",
    "    rna.append('https://uswest.ensembl.org/Homo_sapiens/Transcript/Summary?t=' + rna_id)\n",
    "    if nme != 'None':\n",
    "        lab.append(nme)\n",
    "    else:\n",
    "        lab.append('Ensembl_Transcript_ID:' + rna_id)\n",
    "        nme = 'Ensembl_Transcript_ID:' + rna_id\n",
    "    if ent_type != 'None': desc.append(\"Transcript {} is classified as type '{}'.\".format(nme, ent_type))\n",
    "    else: desc.append('None')\n",
    "    syn.append('None')\n",
    "\n",
    "# combine into new data frame\n",
    "metadata = pandas.DataFrame(list(zip(rna, lab, desc, syn)), columns=['ID', 'Label', 'Description', 'Synonym'])\n",
    "metadata = metadata.astype(str)\n",
    "metadata.drop_duplicates(subset='ID', keep='first', inplace=True)\n",
    "\n",
    "# convert df to dictionary\n",
    "metadata.set_index('ID', inplace=True)\n",
    "rna_metadata_dict = metadata.to_dict('index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### Variant Metadata Dictionary <a class=\"anchor\" id=\"variant-metadata\"></a>  \n",
    "\n",
    "The nested dictionary of rna metadata is created by looping over the human [ClinVar Variant](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#clinvar) identifier data set (`variant_summary.txt`). The `keys` of the dictionary are `dbSNP identifiers` and the `values` are dictionaries for each metadata type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download data\n",
    "url = 'ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'variant_summary.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "\n",
    "# load data\n",
    "var_data = pandas.read_csv(unprocessed_data_location + 'variant_summary.txt', header=0, delimiter='\\t', low_memory=False)\n",
    "\n",
    "# remove rows without identifiers\n",
    "var_data = var_data.loc[var_data['Assembly'].apply(lambda x: x == 'GRCh38')]\n",
    "var_data = var_data.loc[var_data['RS# (dbSNP)'].apply(lambda x: x != -1)]\n",
    "\n",
    "# de-dup data\n",
    "var_metadata = var_data[['#AlleleID', 'Type', 'Name', 'ClinicalSignificance', 'RS# (dbSNP)', 'Origin',\n",
    "                         'ChromosomeAccession', 'Chromosome', 'Start', 'Stop', 'ReferenceAllele',\n",
    "                         'Assembly', 'AlternateAllele','Cytogenetic', 'ReviewStatus', 'LastEvaluated']] \n",
    "\n",
    "# replace NaN with 'None'\n",
    "var_metadata.replace('na', 'None', inplace=True)\n",
    "var_metadata.fillna('None', inplace=True)\n",
    "\n",
    "# remove duplicate dbSNP ids by choosing the most recent reviewed variant\n",
    "var_metadata.sort_values('LastEvaluated', ascending=False, inplace=True)\n",
    "var_metadata.drop_duplicates(subset='RS# (dbSNP)', keep='first', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create metadata\n",
    "variant, label, desc, syn = [], [], [], []\n",
    "for idx, row in tqdm(var_metadata.iterrows(), total=var_metadata.shape[0]):\n",
    "    var_id, lab = row['RS# (dbSNP)'], row['Name']\n",
    "    if var_id != 'None':\n",
    "        variant.append('https://www.ncbi.nlm.nih.gov/snp/rs' + str(var_id))\n",
    "        if lab != 'None': label.append(lab)\n",
    "        else: label.append('dbSNP_ID:rs' + str(var_id))\n",
    "        sent = \"This variant is a {} {} located on chromosome {} ({}, start:{}/stop:{} positions, \" +\\\n",
    "               \"cytogenetic location:{}) and has clinical significance '{}'. \" +\\\n",
    "               \"This entry is for the {} and was last reviewed on {} with review status '{}'.\"\n",
    "        desc.append(sent.format(row['Origin'].replace(';', '/'), row['Type'].replace(';', '/'), row['Chromosome'], row['ChromosomeAccession'],\n",
    "                                row['Start'], row['Stop'], row['Cytogenetic'], row['ClinicalSignificance'],\n",
    "                                row['Assembly'], row['LastEvaluated'], row['ReviewStatus']).replace('None', 'UNKNOWN'))\n",
    "        syn.append('None')\n",
    "    \n",
    "# combine into new data frame\n",
    "var_metadata_final = pandas.DataFrame(list(zip(variant, label, desc, syn)), columns =['ID', 'Label', 'Description', 'Synonym'])\n",
    "var_metadata_final.drop_duplicates(subset=None, keep='first', inplace=True)\n",
    "var_metadata_final = var_metadata_final.astype(str)\n",
    "\n",
    "# convert df to dictionary\n",
    "var_metadata_final.set_index('ID', inplace=True)\n",
    "var_metadata_dict = var_metadata_final.to_dict('index') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### Pathway Metadata Dictionary <a class=\"anchor\" id=\"pathway-metadata\"></a>  \n",
    "\n",
    "The nested dictionary of pathway metadata is created by looping over the human [Reactome Pathway Database](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#reactome-pathway-database) identifier data set (`ReactomePathways.txt`); Reactome-Gene Association data (`gene_association.reactome.gz`), and Reactome-ChEBI data (`ChEBI2Reactome_All_Levels.txt`). The `keys` of the dictionary are `Reactome identifiers` and the `values` are dictionaries for each metadata type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download reactome pathways data\n",
    "url = 'https://reactome.org/download/current/ReactomePathways.txt'\n",
    "if not os.path.exists(unprocessed_data_location + 'ReactomePathways.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "# load data\n",
    "reactome_pathways = pandas.read_csv(unprocessed_data_location + 'ReactomePathways.txt', header=None, delimiter='\\t', low_memory=False)\n",
    "reactome_pathways = reactome_pathways.loc[reactome_pathways[2].apply(lambda x: x == 'Homo sapiens')] \n",
    "\n",
    "# reactome gene association data\n",
    "url = 'https://reactome.org/download/current/gene_association.reactome.gz'\n",
    "if not os.path.exists(unprocessed_data_location + 'gene_association.reactome'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "# load data\n",
    "reactome_pathways2 = pandas.read_csv(unprocessed_data_location + 'gene_association.reactome', header=None, delimiter='\\t', skiprows=3, low_memory=False)\n",
    "reactome_pathways2 = reactome_pathways2.loc[reactome_pathways2[12].apply(lambda x: x == 'taxon:9606')]\n",
    "reactome_pathways2[5].replace('REACTOME:','', inplace=True, regex=True) \n",
    "\n",
    "# reactome CHEBI data\n",
    "url = 'https://reactome.org/download/current/ChEBI2Reactome_All_Levels.txt'\n",
    "if not os.path.exists(unprocessed_data_location + 'ChEBI2Reactome_All_Levels.txt'):\n",
    "    data_downloader(url, unprocessed_data_location)\n",
    "# load data\n",
    "reactome_pathways3 = pandas.read_csv(unprocessed_data_location + 'ChEBI2Reactome_All_Levels.txt', header=None, delimiter='\\t', low_memory=False)\n",
    "# remove all non-human pathways and save as list\n",
    "reactome_pathways3 = reactome_pathways3.loc[reactome_pathways3[5].apply(lambda x: x == 'Homo sapiens')] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get metadata\n",
    "nodes = list(set(reactome_pathways[0]) | set(reactome_pathways2[5]) | set(reactome_pathways3[1]))\n",
    "pathway_metadata_final = metadata_api_mapper(nodes)\n",
    "\n",
    "# update dictionary\n",
    "pathway_metadata_final['ID'] = pathway_metadata_final['ID'].map('https://reactome.org/content/detail/{}'.format)\n",
    "pathway_metadata_final.set_index('ID', inplace=True)\n",
    "\n",
    "# convert df to dictionary\n",
    "pathway_metadata_dict = pathway_metadata_final.to_dict('index') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "#### Relations Metadata Dictionary <a class=\"anchor\" id=\"relations-metadata\"></a>  \n",
    "\n",
    "The nested dictionary of relation metadata is created by looping over the human [Relations Ontology](https://github.com/callahantiff/PheKnowLator/wiki/v2-Data-Sources#relations-ontology) identifier data set (`ro_with_imports.owl`). The `keys` of the dictionary are `Relations Ontology identifiers` and the `values` are dictionaries for each metadata type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download ontology\n",
    "if not os.path.exists(unprocessed_data_location + 'ro_with_imports.owl'):\n",
    "    command = '{} {} --merge-import-closure -o {}'\n",
    "    os.system(command.format(owltools_location, 'http://purl.obolibrary.org/obo/ro.owl',\n",
    "                             unprocessed_data_location + 'ro_with_imports.owl'))\n",
    "# load graph\n",
    "ro_graph = Graph().parse(unprocessed_data_location + 'ro_with_imports.owl')\n",
    "print('There are {} edges in the ontology (date:{})'.format(len(ro_graph), datetime.datetime.now().strftime('%m/%d/%Y')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get metadata\n",
    "relation_metadata_dict, obo = {}, Namespace('http://purl.obolibrary.org/obo/')\n",
    "\n",
    "# get ontology information\n",
    "cls = [x for x in gets_ontology_classes(ro_graph) if '/RO_' in str(x)] +\\\n",
    "      [x for x in gets_object_properties(ro_graph) if '/RO_' in str(x)]\n",
    "master_synonyms = [x for x in ro_graph if 'synonym' in str(x[1]).lower() and isinstance(x[0], URIRef)]\n",
    "\n",
    "for x in tqdm(cls):\n",
    "    # labels\n",
    "    cls_label = [x for x in ro_graph.objects(x, RDFS.label) if '@' not in n3(x) or '@en' in n3(x)]\n",
    "    labels = str(cls_label[0]) if len(cls_label) > 0 else 'None'\n",
    "    # synonyms\n",
    "    cls_syn = [str(i[2]) for i in master_synonyms if x == i[0]]\n",
    "    synonym = str(cls_syn[0]) if len(cls_syn) > 0 else 'None'\n",
    "    # description\n",
    "    cls_desc = [x for x in ro_graph.objects(x, obo.IAO_0000115) if '@' not in n3(x) or '@en' in n3(x)]\n",
    "    desc = '|'.join([str(cls_desc[0])]) if len(cls_desc) > 0 else 'None'\n",
    "    \n",
    "    relation_metadata_dict[str(x)] = {\n",
    "        'Label': labels, 'Description': desc, 'Synonym': synonym\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Create Master Metadata Dictionary** \n",
    "\n",
    "To make it easier to navigate the mapping of each instance node in an edge, a master dictionary is created and keyed by node type. This is most useful when both nodes in an edge are instances, but of different data types (e.g. `gene-rna`).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine all metadata dictionaries\n",
    "master_metadata_dictionary = {'nodes': {**gene_metadata_dict,\n",
    "                                        **rna_metadata_dict,\n",
    "                                        **var_metadata_dict,\n",
    "                                        **pathway_metadata_dict},\n",
    "                              'relations': relation_metadata_dict}\n",
    "\n",
    "# verify metadata strings are properly formatted\n",
    "temp_copy = master_metadata_dictionary.copy(); master_metadata_dictionary = dict()\n",
    "for key, value in tqdm(temp_copy.items()):\n",
    "    master_metadata_dictionary[key] = {}\n",
    "    for ent_key, ent_value in value.items():\n",
    "        updated_inner_dict = {k: re.sub('\\s\\s+', ' ', v.replace('\\n', ' '))\n",
    "                              if v is not None else v for k, v in ent_value.items()}\n",
    "        master_metadata_dictionary[key][ent_key] = updated_inner_dict\n",
    "del temp_copy\n",
    "\n",
    "# save dictionary locally\n",
    "pickle.dump(master_metadata_dictionary, open(node_data_location + 'node_metadata_dict.pkl', 'wb'), protocol=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<br>\n",
    "\n",
    "***\n",
    "***\n",
    "\n",
    "```\n",
    "@misc{callahan_tj_2019_3401437,\n",
    "  author       = {Callahan, TJ},\n",
    "  title        = {PheKnowLator},\n",
    "  month        = mar,\n",
    "  year         = 2019,\n",
    "  doi          = {10.5281/zenodo.3401437},\n",
    "  url          = {https://doi.org/10.5281/zenodo.3401437}\n",
    "}\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}