AddDB_updater.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AddDB Updater\n",
    "\n",
    "This is a Jupyter notebook to build a python script that can parse and create a fullly annotated tsv file to replace gene_xref for ANNOVAR. This script could be useful to update these database and be added to the Achabilarity container.\n",
    "\n",
    "# Import library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parser function\n",
    "\n",
    "These functions will parse them for merging.\n",
    "\n",
    "## HGNC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hgnc(file):\n",
    "    hgnc = pd.read_csv(file, sep='\\t')\n",
    "    hgnc.rename(columns={'Approved symbol': '#Gene_name'}, inplace=True)\n",
    "    hgnc = hgnc.sort_values(by=['#Gene_name'])\n",
    "    hgnc = hgnc.reset_index(drop=True)\n",
    "    return(hgnc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  #Gene_name                   Approved name\n",
      "0       A1BG          alpha-1-B glycoprotein\n",
      "1   A1BG-AS1            A1BG antisense RNA 1\n",
      "2       A1CF  APOBEC1 complementation factor\n",
      "3        A2M           alpha-2-macroglobulin\n",
      "4    A2M-AS1             A2M antisense RNA 1\n",
      "\n",
      "Info\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 41782 entries, 0 to 41781\n",
      "Data columns (total 2 columns):\n",
      "#Gene_name       41782 non-null object\n",
      "Approved name    41782 non-null object\n",
      "dtypes: object(2)\n",
      "memory usage: 652.9+ KB\n",
      "None\n",
      "\n",
      "Duplicates if exist\n",
      "Empty DataFrame\n",
      "Columns: [#Gene_name, Approved name]\n",
      "Index: []\n"
     ]
    }
   ],
   "source": [
    "hgnc_list = hgnc('~/Kevin/AddDB_updater/data/hgnc.tsv')\n",
    "print(hgnc_list.head())\n",
    "print(\"\\nInfo\")\n",
    "print(hgnc_list.info())\n",
    "print(\"\\nDuplicates if exist\")\n",
    "print(hgnc_list[hgnc_list.duplicated() == True])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GnomAD constraint score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each gene, if there is multiples entries because of multiples transcripts, we will only keep the highest predicted impact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gnomad_score(file):\n",
    "    gnomad = pd.read_csv(file, sep='\\t',header=0)\n",
    "    gnomad.rename(columns={'gene': '#Gene_name'}, inplace=True)\n",
    "    gnomad_select = gnomad[['#Gene_name', 'oe_lof_upper_rank',\n",
    "       'oe_lof_upper_bin','oe_lof','oe_lof_lower','oe_lof_upper','oe_mis','oe_mis_lower', 'oe_mis_upper','oe_syn','oe_syn_lower', 'oe_syn_upper','constraint_flag']]\n",
    "    gnomad_select = gnomad_select.sort_values(by=['#Gene_name'])\n",
    "    gnomad_select = gnomad_select.groupby('#Gene_name', as_index=False).min()\n",
    "    gnomad_select = gnomad_select.reset_index(drop=True)\n",
    "    return(gnomad_select)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  #Gene_name  oe_lof_upper_rank  oe_lof_upper_bin   oe_lof  oe_lof_lower  \\\n",
      "0       A1BG            13015.0               6.0  0.78457         0.524   \n",
      "1       A1CF             9254.0               4.0  0.60537         0.425   \n",
      "2        A2M             5366.0               2.0  0.40526         0.305   \n",
      "3      A2ML1            10116.0               5.0  0.77171         0.629   \n",
      "4    A3GALT2            17288.0               9.0  1.17270         0.771   \n",
      "\n",
      "   oe_lof_upper   oe_mis  oe_mis_lower  oe_mis_upper   oe_syn  oe_syn_lower  \\\n",
      "0         1.208  1.01410         0.922         1.116  1.02990         0.897   \n",
      "1         0.880  0.84521         0.765         0.934  1.03550         0.895   \n",
      "2         0.544  0.81065         0.758         0.866  0.87995         0.796   \n",
      "3         0.952  0.96729         0.911         1.027  1.00480         0.916   \n",
      "4         1.764  0.99002         0.873         1.124  0.98196         0.816   \n",
      "\n",
      "   oe_syn_upper  \n",
      "0         1.184  \n",
      "1         1.201  \n",
      "2         0.973  \n",
      "3         1.103  \n",
      "4         1.187  \n",
      "     #Gene_name  oe_lof_upper_rank  oe_lof_upper_bin   oe_lof  oe_lof_lower  \\\n",
      "1181       AQP1             9257.0               4.0  0.34085         0.154   \n",
      "\n",
      "      oe_lof_upper   oe_mis  oe_mis_lower  oe_mis_upper  oe_syn  oe_syn_lower  \\\n",
      "1181         0.881  0.88463         0.776          1.01  1.0934         0.945   \n",
      "\n",
      "      oe_syn_upper  \n",
      "1181         1.262  \n",
      "\n",
      "Info\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 19658 entries, 0 to 19657\n",
      "Data columns (total 12 columns):\n",
      "#Gene_name           19658 non-null object\n",
      "oe_lof_upper_rank    19155 non-null float64\n",
      "oe_lof_upper_bin     19155 non-null float64\n",
      "oe_lof               19155 non-null float64\n",
      "oe_lof_lower         19155 non-null float64\n",
      "oe_lof_upper         19155 non-null float64\n",
      "oe_mis               19657 non-null float64\n",
      "oe_mis_lower         19658 non-null float64\n",
      "oe_mis_upper         19658 non-null float64\n",
      "oe_syn               19657 non-null float64\n",
      "oe_syn_lower         19658 non-null float64\n",
      "oe_syn_upper         19658 non-null float64\n",
      "dtypes: float64(11), object(1)\n",
      "memory usage: 1.8+ MB\n",
      "None\n",
      "\n",
      "Duplicates if exist\n",
      "Empty DataFrame\n",
      "Columns: [#Gene_name, oe_lof_upper_rank, oe_lof_upper_bin, oe_lof, oe_lof_lower, oe_lof_upper, oe_mis, oe_mis_lower, oe_mis_upper, oe_syn, oe_syn_lower, oe_syn_upper]\n",
      "Index: []\n"
     ]
    }
   ],
   "source": [
    "gnomad_score_list =  gnomad_score('~/Kevin/AddDB_updater/data/gnomad.v2.1.1.lof_metrics.by_gene.txt')\n",
    "print(gnomad_score_list.head())\n",
    "print(gnomad_score_list[gnomad_score_list['#Gene_name'] == 'AQP1'])\n",
    "print(\"\\nInfo\")\n",
    "print(gnomad_score_list.info())\n",
    "print(\"\\nDuplicates if exist\")\n",
    "print(gnomad_score_list[gnomad_score_list.duplicated() == True])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OMIM genemap2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As gene could have multiples entries, we keep the one with the highest Phenotypes informations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def omim(file):\n",
    "    omim = pd.read_csv(file, sep='\\t',skiprows=3)\n",
    "    omim_select = omim[['Approved Symbol','Phenotypes']]\n",
    "    omim_select.rename(columns={'Approved Symbol': '#Gene_name'}, inplace=True)\n",
    "    omim_select = omim_select.dropna(subset=['#Gene_name']) \n",
    "    omim_select = omim_select.sort_values(by=['#Gene_name'])\n",
    "    omim_select = omim_select.fillna('')\n",
    "    omim_select = omim_select.groupby('#Gene_name', as_index=False)['Phenotypes'].max()\n",
    "    omim_select = omim_select.drop_duplicates()\n",
    "    omim_select= omim_select.reset_index(drop=True)\n",
    "    return(omim_select)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  #Gene_name                                         Phenotypes\n",
      "0       A1BG                                                   \n",
      "1       A1CF                                                   \n",
      "2        A2M  Alpha-2-macroglobulin deficiency, 614036 (1), ...\n",
      "3      A2ML1  {Otitis media, susceptibility to}, 166760 (3),...\n",
      "4     A4GALT  [Blood group, P1Pk system, P(2) phenotype], 11...\n",
      "5      A4GNT                                                   \n",
      "6       AAAS  Achalasia-addisonianism-alacrimia syndrome, 23...\n",
      "7       AACS                                                   \n",
      "8      AADAC                                                   \n",
      "9      AADAT                                                   \n",
      "\n",
      "Info\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 15878 entries, 0 to 15877\n",
      "Data columns (total 2 columns):\n",
      "#Gene_name    15878 non-null object\n",
      "Phenotypes    15878 non-null object\n",
      "dtypes: object(2)\n",
      "memory usage: 248.2+ KB\n",
      "None\n",
      "\n",
      "Duplicates if exist\n",
      "Empty DataFrame\n",
      "Columns: [#Gene_name, Phenotypes]\n",
      "Index: []\n",
      "     #Gene_name                                         Phenotypes\n",
      "2914     CSF2RA  Surfactant metabolism dysfunction, pulmonary, ...\n"
     ]
    }
   ],
   "source": [
    "omim_list = omim('~/Kevin/AddDB_updater/data/genemap2.txt').reset_index(drop=True)\n",
    "print(omim_list.head(n=10))\n",
    "print(\"\\nInfo\")\n",
    "print(omim_list.info())\n",
    "print(\"\\nDuplicates if exist\")\n",
    "print(omim_list[omim_list.duplicated() == True])\n",
    "print(omim_list[omim_list['#Gene_name'] == 'CSF2RA'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## UniProt database"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As Uniprot give multiples genes names, we need to split them into single rows.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def explode(df, lst_cols, fill_value='', preserve_index=False):\n",
    "    # make sure `lst_cols` is list-alike\n",
    "    if (lst_cols is not None\n",
    "        and len(lst_cols) > 0\n",
    "        and not isinstance(lst_cols, (list, tuple, np.ndarray, pd.Series))):\n",
    "        lst_cols = [lst_cols]\n",
    "    # all columns except `lst_cols`\n",
    "    idx_cols = df.columns.difference(lst_cols)\n",
    "    # calculate lengths of lists\n",
    "    lens = df[lst_cols[0]].str.len()\n",
    "    # preserve original index values    \n",
    "    idx = np.repeat(df.index.values, lens)\n",
    "    # create \"exploded\" DF\n",
    "    res = (pd.DataFrame({\n",
    "                col:np.repeat(df[col].values, lens)\n",
    "                for col in idx_cols},\n",
    "                index=idx)\n",
    "             .assign(**{col:np.concatenate(df.loc[lens>0, col].values)\n",
    "                            for col in lst_cols}))\n",
    "    # append those rows that have empty lists\n",
    "    if (lens == 0).any():\n",
    "        # at least one list in cells is empty\n",
    "        res = (res.append(df.loc[lens==0, idx_cols], sort=False)\n",
    "                  .fillna(fill_value))\n",
    "    # revert the original index order\n",
    "    res = res.sort_index()\n",
    "    # reset index if requested\n",
    "    if not preserve_index:        \n",
    "        res = res.reset_index(drop=True)\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Process file. As some gene have multiples entries in Uniprot, we only keep the first iteration (usually the longest transcript)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def uniprot(file):\n",
    "    uniprot = pd.read_csv(file, sep='\\t')\n",
    "    uniprot_select = uniprot.iloc[:,3:]\n",
    "    uniprot_select.rename(columns={'Gene names  (primary )': 'Gene_name'}, inplace=True)\n",
    "    uniprot_select = uniprot_select.dropna(subset=['Gene_name']) \n",
    "    uniprot_select = explode(uniprot_select.assign(Gene_name=uniprot_select['Gene_name'].str.split(';')),lst_cols=['Gene_name'])\n",
    "    uniprot_select.rename(columns={'Gene_name': '#Gene_name'}, inplace=True)\n",
    "    uniprot_select = uniprot_select.groupby('#Gene_name', as_index=False).min()\n",
    "    uniprot_select = uniprot_select.sort_values(by=['#Gene_name'])\n",
    "    uniprot_select = uniprot_select.drop_duplicates()\n",
    "    uniprot_select = uniprot_select.reset_index(drop=True)\n",
    "    return(uniprot_select)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mobidic-umai/anaconda3/lib/python3.7/site-packages/pandas/core/frame.py:4025: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  return super(DataFrame, self).rename(**kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  #Gene_name                                      Function [CC]  \\\n",
      "0      AMY1B                                               None   \n",
      "1      AMY1C                                               None   \n",
      "2     ARL17B  FUNCTION: GTP-binding protein that functions a...   \n",
      "3     BOLA2B  FUNCTION: Acts as a cytosolic iron-sulfur (Fe-...   \n",
      "4      BPY2B                                               None   \n",
      "5      BPY2C                                               None   \n",
      "6      C4B_2  FUNCTION: Non-enzymatic component of the C3 an...   \n",
      "7     CCL3L3  FUNCTION: Chemotactic for lymphocytes and mono...   \n",
      "8     CCL4L2  FUNCTION: Chemokine that induces chemotaxis of...   \n",
      "9      CDY1B  FUNCTION: Has histone acetyltransferase activi...   \n",
      "\n",
      "                              Involvement in disease  \\\n",
      "0                                               None   \n",
      "1                                               None   \n",
      "2                                               None   \n",
      "3                                               None   \n",
      "4                                               None   \n",
      "5                                               None   \n",
      "6  DISEASE: Systemic lupus erythematosus (SLE) [M...   \n",
      "7                                               None   \n",
      "8                                               None   \n",
      "9                                               None   \n",
      "\n",
      "                                  Tissue specificity  \n",
      "0                                               None  \n",
      "1                                               None  \n",
      "2                                               None  \n",
      "3                                               None  \n",
      "4  TISSUE SPECIFICITY: Expressed exclusively in t...  \n",
      "5  TISSUE SPECIFICITY: Expressed exclusively in t...  \n",
      "6  TISSUE SPECIFICITY: Complement component C4 is...  \n",
      "7                                               None  \n",
      "8  TISSUE SPECIFICITY: Detected in B-cells. {ECO:...  \n",
      "9  TISSUE SPECIFICITY: Testis-specific. Detected ...  \n",
      "\n",
      "Info\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 20228 entries, 0 to 20227\n",
      "Data columns (total 4 columns):\n",
      "#Gene_name                20228 non-null object\n",
      "Function [CC]             16288 non-null object\n",
      "Involvement in disease    4246 non-null object\n",
      "Tissue specificity        9806 non-null object\n",
      "dtypes: object(4)\n",
      "memory usage: 632.2+ KB\n",
      "None\n",
      "\n",
      "Duplicates if exist\n",
      "Empty DataFrame\n",
      "Columns: [#Gene_name, Function [CC], Involvement in disease, Tissue specificity]\n",
      "Index: []\n"
     ]
    }
   ],
   "source": [
    "uniprot_list = uniprot('~/Kevin/AddDB_updater/data/uniprot.tsv')\n",
    "print(uniprot_list.head(n=10))\n",
    "print(\"\\nInfo\")\n",
    "print(uniprot_list.info())\n",
    "print(\"\\nDuplicates if exist\")\n",
    "print(uniprot_list[uniprot_list.duplicated() == True])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     #Gene_name                                      Function [CC]  \\\n",
      "1606       BBC3  FUNCTION: Essential mediator of p53/TP53-depen...   \n",
      "\n",
      "     Involvement in disease Tissue specificity  \n",
      "1606                    NaN                NaN  \n"
     ]
    }
   ],
   "source": [
    "print(uniprot_list[uniprot_list['#Gene_name'] == 'BBC3'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merge into one file\n",
    "\n",
    "This function will merge all databases into HGNC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_db(hgnc_file,omim_file,gnomad_score_file,uniprot_file):\n",
    "    hgnc_list = hgnc(hgnc_file)\n",
    "    omim_list = omim(omim_file)\n",
    "    gnomad_score_list = gnomad_score(gnomad_score_file)\n",
    "    uniprot_list = uniprot(uniprot_file)\n",
    "    gene_fullxref = hgnc_list.merge(omim_list,on='#Gene_name').merge(gnomad_score_list,on='#Gene_name').merge(uniprot_list,on='#Gene_name')\n",
    "    gene_fullxref = gene_fullxref.sort_values(by=['#Gene_name'])\n",
    "    gene_fullxref = gene_fullxref.drop_duplicates()\n",
    "    gene_fullxref = gene_fullxref.fillna('')\n",
    "    return gene_fullxref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mobidic-umai/anaconda3/lib/python3.7/site-packages/pandas/core/frame.py:4025: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  return super(DataFrame, self).rename(**kwargs)\n",
      "/home/mobidic-umai/anaconda3/lib/python3.7/site-packages/pandas/core/frame.py:4025: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  return super(DataFrame, self).rename(**kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  #Gene_name                                    Approved name  \\\n",
      "0       A1BG                           alpha-1-B glycoprotein   \n",
      "1       A1CF                   APOBEC1 complementation factor   \n",
      "2        A2M                            alpha-2-macroglobulin   \n",
      "3      A2ML1                     alpha-2-macroglobulin like 1   \n",
      "4     A4GALT  alpha 1,4-galactosyltransferase (P blood group)   \n",
      "\n",
      "                                          Phenotypes oe_lof_upper_rank  \\\n",
      "0                                                                13015   \n",
      "1                                                                 9254   \n",
      "2  Alpha-2-macroglobulin deficiency, 614036 (1), ...              5366   \n",
      "3  {Otitis media, susceptibility to}, 166760 (3),...             10116   \n",
      "4  [Blood group, P1Pk system, P(2) phenotype], 11...             16517   \n",
      "\n",
      "  oe_lof_upper_bin   oe_lof oe_lof_lower oe_lof_upper   oe_mis  oe_mis_lower  \\\n",
      "0                6  0.78457        0.524        1.208   1.0141         0.922   \n",
      "1                4  0.60537        0.425         0.88  0.84521         0.765   \n",
      "2                2  0.40526        0.305        0.544  0.81065         0.758   \n",
      "3                5  0.77171        0.629        0.952  0.96729         0.911   \n",
      "4                8  0.94609        0.553        1.654   1.0345         0.930   \n",
      "\n",
      "   oe_mis_upper   oe_syn  oe_syn_lower  oe_syn_upper  \\\n",
      "0         1.116   1.0299         0.897         1.184   \n",
      "1         0.934   1.0355         0.895         1.201   \n",
      "2         0.866  0.87995         0.796         0.973   \n",
      "3         1.027   1.0048         0.916         1.103   \n",
      "4         1.151   1.1197         0.967         1.300   \n",
      "\n",
      "                                       Function [CC]  \\\n",
      "0                                                      \n",
      "1  FUNCTION: Essential component of the apolipopr...   \n",
      "2  FUNCTION: Is able to inhibit all four classes ...   \n",
      "3  FUNCTION: Is able to inhibit all four classes ...   \n",
      "4  FUNCTION: Necessary for the biosynthesis of th...   \n",
      "\n",
      "                              Involvement in disease  \\\n",
      "0                                                      \n",
      "1                                                      \n",
      "2                                                      \n",
      "3  DISEASE: Otitis media (OM) [MIM:166760]: An in...   \n",
      "4                                                      \n",
      "\n",
      "                                  Tissue specificity  \n",
      "0                        TISSUE SPECIFICITY: Plasma.  \n",
      "1  TISSUE SPECIFICITY: Widely expressed with high...  \n",
      "2  TISSUE SPECIFICITY: Secreted in plasma. {ECO:0...  \n",
      "3  TISSUE SPECIFICITY: In the epidermis, expresse...  \n",
      "4  TISSUE SPECIFICITY: Ubiquitous. Highly express...  \n",
      "\n",
      "Info\n",
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 14065 entries, 0 to 14064\n",
      "Data columns (total 17 columns):\n",
      "#Gene_name                14065 non-null object\n",
      "Approved name             14065 non-null object\n",
      "Phenotypes                14065 non-null object\n",
      "oe_lof_upper_rank         14065 non-null object\n",
      "oe_lof_upper_bin          14065 non-null object\n",
      "oe_lof                    14065 non-null object\n",
      "oe_lof_lower              14065 non-null object\n",
      "oe_lof_upper              14065 non-null object\n",
      "oe_mis                    14065 non-null object\n",
      "oe_mis_lower              14065 non-null float64\n",
      "oe_mis_upper              14065 non-null float64\n",
      "oe_syn                    14065 non-null object\n",
      "oe_syn_lower              14065 non-null float64\n",
      "oe_syn_upper              14065 non-null float64\n",
      "Function [CC]             14065 non-null object\n",
      "Involvement in disease    14065 non-null object\n",
      "Tissue specificity        14065 non-null object\n",
      "dtypes: float64(4), object(13)\n",
      "memory usage: 1.9+ MB\n",
      "None\n",
      "\n",
      "Duplicates if exist\n",
      "Empty DataFrame\n",
      "Columns: [#Gene_name, Approved name, Phenotypes, oe_lof_upper_rank, oe_lof_upper_bin, oe_lof, oe_lof_lower, oe_lof_upper, oe_mis, oe_mis_lower, oe_mis_upper, oe_syn, oe_syn_lower, oe_syn_upper, Function [CC], Involvement in disease, Tissue specificity]\n",
      "Index: []\n"
     ]
    }
   ],
   "source": [
    "gene_fullxref_list = merge_db(hgnc_file='~/Kevin/AddDB_updater/data/hgnc.tsv',omim_file='~/Kevin/AddDB_updater/data/genemap2.txt', gnomad_score_file='~/Kevin/AddDB_updater/data/gnomad.v2.1.1.lof_metrics.by_gene.txt', uniprot_file='~/Kevin/AddDB_updater/data/uniprot.tsv')\n",
    "gene_fullxref_list.to_csv('~/Kevin/AddDB_updater/data/test2.txt',sep='\\t',index=False)\n",
    "print(gene_fullxref_list.head())\n",
    "print(\"\\nInfo\")\n",
    "print(gene_fullxref_list.info())\n",
    "print(\"\\nDuplicates if exist\")\n",
    "print(gene_fullxref_list[gene_fullxref_list.duplicated() == True])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove characters that are not permitted in regular expression into \"_\" (used in ANNOVAR for -gx). Note that for OMIM file, there is some brackets that is missing, thus give problems for regex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_fullxref_list = pd.read_csv('~/Kevin/AddDB_updater/data/test2.txt',sep='\\t')\n",
    "#COMMAND = \" sed 's/+/plus/g' | awk 'BEGIN{FS=OFS=\"\\t\"} {for (i=6;i<=7;i++) gsub(/-/,\"_\",$i)}1' |  awk 'BEGIN{FS=OFS=\"\\t\"} {gsub(/-/,\"_\",$(NF-1))}1'\" \n",
    "#test = subprocess.call(['sed', 's/+/plus/g', test])\n",
    "gene_fullxref_list = gene_fullxref_list.replace('\\+','plus',regex=True)\n",
    "gene_fullxref_list.loc[:,['Approved name','Phenotypes','Function [CC]','Involvement in disease','Tissue specificity']] = gene_fullxref_list.loc[:,['Approved name','Phenotypes','Function [CC]','Involvement in disease','Tissue specificity']].replace('-','_',regex=True)\n",
    "gene_fullxref_list.loc[:,['Phenotypes']] = gene_fullxref_list.loc[:,['Phenotypes']].replace('\\(','_',regex=True)\n",
    "gene_fullxref_list.loc[:,['Phenotypes']] = gene_fullxref_list.loc[:,['Phenotypes']].replace('\\)','_',regex=True)\n",
    "gene_fullxref_list.columns = gene_fullxref_list.columns.str.replace(' ','_')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  #Gene_name                                    Approved_name  \\\n",
      "0       A1BG                           alpha_1_B glycoprotein   \n",
      "1       A1CF                   APOBEC1 complementation factor   \n",
      "2        A2M                            alpha_2_macroglobulin   \n",
      "3      A2ML1                     alpha_2_macroglobulin like 1   \n",
      "4     A4GALT  alpha 1,4_galactosyltransferase (P blood group)   \n",
      "\n",
      "                                          Phenotypes  oe_lof_upper_rank  \\\n",
      "0                                                NaN            13015.0   \n",
      "1                                                NaN             9254.0   \n",
      "2  Alpha_2_macroglobulin deficiency, 614036 _1_, ...             5366.0   \n",
      "3  {Otitis media, susceptibility to}, 166760 _3_,...            10116.0   \n",
      "4  [Blood group, P1Pk system, P_2_ phenotype], 11...            16517.0   \n",
      "\n",
      "   oe_lof_upper_bin   oe_lof  oe_lof_lower  oe_lof_upper   oe_mis  \\\n",
      "0               6.0  0.78457         0.524         1.208  1.01410   \n",
      "1               4.0  0.60537         0.425         0.880  0.84521   \n",
      "2               2.0  0.40526         0.305         0.544  0.81065   \n",
      "3               5.0  0.77171         0.629         0.952  0.96729   \n",
      "4               8.0  0.94609         0.553         1.654  1.03450   \n",
      "\n",
      "   oe_mis_lower  oe_mis_upper   oe_syn  oe_syn_lower  oe_syn_upper  \\\n",
      "0         0.922         1.116  1.02990         0.897         1.184   \n",
      "1         0.765         0.934  1.03550         0.895         1.201   \n",
      "2         0.758         0.866  0.87995         0.796         0.973   \n",
      "3         0.911         1.027  1.00480         0.916         1.103   \n",
      "4         0.930         1.151  1.11970         0.967         1.300   \n",
      "\n",
      "                                       Function_[CC]  \\\n",
      "0                                                NaN   \n",
      "1  FUNCTION: Essential component of the apolipopr...   \n",
      "2  FUNCTION: Is able to inhibit all four classes ...   \n",
      "3  FUNCTION: Is able to inhibit all four classes ...   \n",
      "4  FUNCTION: Necessary for the biosynthesis of th...   \n",
      "\n",
      "                              Involvement_in_disease  \\\n",
      "0                                                NaN   \n",
      "1                                                NaN   \n",
      "2                                                NaN   \n",
      "3  DISEASE: Otitis media (OM) [MIM:166760]: An in...   \n",
      "4                                                NaN   \n",
      "\n",
      "                                  Tissue_specificity  \n",
      "0                        TISSUE SPECIFICITY: Plasma.  \n",
      "1  TISSUE SPECIFICITY: Widely expressed with high...  \n",
      "2  TISSUE SPECIFICITY: Secreted in plasma. {ECO:0...  \n",
      "3  TISSUE SPECIFICITY: In the epidermis, expresse...  \n",
      "4  TISSUE SPECIFICITY: Ubiquitous. Highly express...  \n"
     ]
    }
   ],
   "source": [
    "print(gene_fullxref_list.head())\n",
    "gene_fullxref_list.to_csv('~/Kevin/AddDB_updater/data/aaa.txt',sep='\\t',index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}