banana_conjecture/The_Banana_Conjecture.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "- Description of Notebook: This analysis compares the number of orthologs between banana and human using OMA, OrthoInspector, and best-bidirectional hit.\n",
    "- Last updated: 8 Dec 2020\n",
    "- By: Natasha Glover"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-08T16:13:22.679716Z",
     "start_time": "2020-12-08T16:13:20.694663Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#import libraries\n",
    "from omadb import Client\n",
    "import requests\n",
    "import json\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "#seaborn options\n",
    "sns.set(rc={'axes.facecolor':'white', 'figure.facecolor':'white'}, font_scale=1.5)\n",
    "sns.set_style('whitegrid')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Human-Banana orthologs using OrthoInspector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "*OrthoInspector is a software suite for inference of orthologous relationships between protein coding-genes and an online resource to access and query precomputed orthology databases.* https://lbgi.fr/orthoinspectorv3/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-08T16:13:22.685409Z",
     "start_time": "2020-12-08T16:13:22.682007Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "human_ncbi_taxid = str(9606)\n",
    "banana_ncbi_taxid = str(214687)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-08T16:13:27.699950Z",
     "start_time": "2020-12-08T16:13:22.688455Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 7668 banana genes with a human ortholog.\n",
      "There are 36440 total number of genes in banana.\n",
      "Using OrthoInspector, 21% of the banana genome's proteins have an ortholog in human.\n"
     ]
    }
   ],
   "source": [
    "#get all the banana genes with an ortholog in human\n",
    "response = requests.get(\"http://lbgi.fr/orthoinspectorv3/api/Eukaryota/species/\"+ banana_ncbi_taxid +\\\n",
    "                        \"/orthologs/\"+ human_ncbi_taxid)\n",
    "\n",
    "#count the number of banana genes with a human ortholog\n",
    "banana_genes_with_orthologs = response.json()\n",
    "\n",
    "print(\"There are {} banana genes with a human ortholog.\".format(len(banana_genes_with_orthologs.keys())))\n",
    "\n",
    "#count the number of banana genes\n",
    "banana_genome = requests.get(\"http://lbgi.fr/orthoinspectorv3/api/Eukaryota/species/214687/proteins\").json()\n",
    "\n",
    "print(\"There are {} total number of genes in banana.\".format(len(banana_genome)))\n",
    "\n",
    "#percentage of banana genes with a human ortholog\n",
    "percentage = len(banana_genes_with_orthologs.keys())/len(banana_genome)*100\n",
    "print(\"Using OrthoInspector, {:.0f}% of the banana genome's proteins have an ortholog in human.\".format(percentage))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Same thing with human."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-08T16:13:32.340301Z",
     "start_time": "2020-12-08T16:13:27.703489Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 5058 human genes with a banana ortholog.\n",
      "There are 21006 total number of genes in human.\n",
      "Using OrthoInspector, 24% of the human genome's proteins have an ortholog in banana.\n"
     ]
    }
   ],
   "source": [
    "#get all the human genes with an ortholog in banana\n",
    "response = requests.get(\"http://lbgi.fr/orthoinspectorv3/api/Eukaryota/species/\"+ human_ncbi_taxid +\\\n",
    "                        \"/orthologs/\"+ banana_ncbi_taxid)\n",
    "\n",
    "#count the number of human  genes with a banana ortholog\n",
    "human_genes_with_orthologs = response.json()\n",
    "\n",
    "print(\"There are {} human genes with a banana ortholog.\".format(len(human_genes_with_orthologs.keys())))\n",
    "\n",
    "#count the number of human genes\n",
    "human_genome = requests.get(\"http://lbgi.fr/orthoinspectorv3/api/Eukaryota/species/\"+ human_ncbi_taxid +\"/proteins\").json()\n",
    "\n",
    "print(\"There are {} total number of genes in human.\".format(len(human_genome)))\n",
    "\n",
    "#percentage of human genes with a banana ortholog\n",
    "percentage = len(human_genes_with_orthologs.keys())/len(human_genome)*100\n",
    "print(\"Using OrthoInspector, {:.0f}% of the human genome's proteins have an ortholog in banana.\".format(percentage))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Human-Banana orthologs using OMA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Here, I use PyOMADB, the user-friendly Python wrapper around the OMA REST API. You can install it with \"pip install omadb\" from the command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-08T16:13:32.352111Z",
     "start_time": "2020-12-08T16:13:32.342888Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "c = Client()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "First, get the number of genes in the human and banana genomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-08T16:13:53.639384Z",
     "start_time": "2020-12-08T16:13:32.355962Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 20159 total number of genes in human.\n"
     ]
    }
   ],
   "source": [
    "#count number of genes in the human genome\n",
    "human_genome = c.genomes.proteins(\"HUMAN\").as_dataframe()\n",
    "human_genome = human_genome[human_genome['is_main_isoform']==True]\n",
    "\n",
    "print(\"There are {} total number of genes in human.\".format(len(human_genome)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-12-08T16:14:19.296882Z",
     "start_time": "2020-12-08T16:13:53.642944Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 36439 total number of genes in banana.\n"
     ]
    }
   ],
   "source": [
    "#count number of genes in the banana genome\n",
    "banana_genome = c.genomes.proteins(\"MUSAM\").as_dataframe()\n",
    "banana_genome = banana_genome[banana_genome['is_main_isoform']==True]\n",
    "\n",
    "print(\"There are {} total number of genes in banana.\".format(len(banana_genome)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Next, get the orthologs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.722Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#get the pairwise orthologs between human and banana\n",
    "human_banana_pairs = list(c.pairwise('HUMAN', 'MUSAM'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Finally, compute the percentage of genes with an ortholog in each species:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.725Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#count the number and percent of human genes with an ortholog in banana\n",
    "\n",
    "#take the set to get a unique list, as a given gene might have multiple orthologs\n",
    "human_genes_with_orthologs = set([x['entry_1']['omaid'] for x in human_banana_pairs if x['entry_1']['is_main_isoform']==True])\n",
    "\n",
    "human_percentage = len(human_genes_with_orthologs)/len(human_genome)* 100\n",
    "\n",
    "print(\"Using OMA, {:.0f}% of the human genome's proteins have an ortholog in banana.\".format(human_percentage), \"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.729Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#count the number and percent of banana genes with an ortholog in human\n",
    "\n",
    "banana_genes_with_orthologs = set([x['entry_2']['omaid'] for x in human_banana_pairs if x['entry_2']['is_main_isoform']==True])\n",
    "\n",
    "banana_percentage = len(banana_genes_with_orthologs)/len(banana_genome)* 100\n",
    "\n",
    "print(\"Using OMA, {:.0f}% of the banana genome's proteins have an ortholog in human.\".format(banana_percentage))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Human-banana orthologs using Best Bidirectional Hit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Update:\n",
    "- Downloaded human genome from UniProt: ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/reference_proteomes/Eukaryota/UP000005640_9606.fasta.gz\n",
    "- This reference proteome has 20609 genes, one representative protein per locus.\n",
    "\n",
    "- Downloaded the banana genome from UniProt: ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/reference_proteomes/Eukaryota/UP000012960_214687.fasta.gz\n",
    "- 36474 genes (1 protein sequence per gene)\n",
    "\n",
    "Commands for running blast:\n",
    "- source /dcsrsoft/spack/bin/setup_dcsrsoft\n",
    "- module load gcc\n",
    "- module load blast-plus/2.9.0\n",
    "- makeblastdb -in human_genome.faa -dbtype prot\n",
    "- makeblastdb -in banana_genome.faa -dbtype prot\n",
    "- mkdir banana_input\n",
    "- mkdir human_input\n",
    "- pyfasta split -n 100 human_genome.faa\n",
    "- mv human_genome.[0-9]*.faa human_input/\n",
    "- pyfasta split -n 100 banana_genome.faa\n",
    "- mv banana_genome.[0-9]*.faa banana_input/\n",
    "\n",
    "The next part uses the SLURM script.\n",
    "- sbatch blast_banana_vs_human.sh\n",
    "- sbatch blast_human_vs_banana.sh\n",
    "\n",
    "- cat banana_vs_human_output/* > banana_vs_human.blastp\n",
    "- cat human_vs_banana_output/* > human_vs_banana.blastp\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.733Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "header = [\"qseqid\",\"sseqid\",\"pident\",\"length\",\"mismatch\",\"gapopen\",\"qstart\",\"qend\",\"qlen\", \\\n",
    "          \"sstart\",\"send\",\"slen\",\"evalue\",\"bitscore\"]\n",
    "human_df = pd.read_csv(\"human_vs_banana.blastp\", sep=\"\\t\", header=None)\n",
    "human_df.columns = header\n",
    "human_df = human_df.drop_duplicates(subset=\"qseqid\", keep=\"first\")\n",
    "human_df.rename({\"qseqid\":\"human_gene\", \"sseqid\":\"banana_gene\"}, axis=1, inplace=True)\n",
    "print(len(human_df))\n",
    "      \n",
    "human_df[:3]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.738Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "banana_df = pd.read_csv(\"banana_vs_human.blastp\", sep=\"\\t\", header=None)\n",
    "banana_df.columns = header\n",
    "banana_df = banana_df.drop_duplicates(subset=\"qseqid\", keep=\"first\")\n",
    "banana_df.rename({\"qseqid\":\"banana_gene\", \"sseqid\":\"human_gene\"}, axis=1, inplace=True)\n",
    "print(len(banana_df))\n",
    "banana_df[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.742Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bbh_df = pd.merge(left=human_df, right=banana_df, on=['human_gene', 'banana_gene'], how=\"inner\")\n",
    "bbh_df[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.745Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "nb_bbhs = len(bbh_df)\n",
    "nb_human_genes_uniprot = 20609\n",
    "nb_banana_genes_uniprot = 36474\n",
    "print(\"There are {} genes in the uniprot reference human proteome.\"\\\n",
    "      .format(nb_human_genes_uniprot))\n",
    "print(\"There are {} genes in the uniprot reference banana proteome.\"\\\n",
    "      .format(nb_banana_genes_uniprot))\n",
    "print(\"There are {} best bidirectional hits between human and banana.\"\\\n",
    "      .format(nb_bbhs))\n",
    "print(\"This is {:.0f}% of the human genome, and {:.0f}% of the banana genome.\"\\\n",
    "      .format(nb_bbhs/nb_human_genes_uniprot*100, nb_bbhs/nb_banana_genes_uniprot*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.748Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "summary_df = pd.DataFrame(data={\"method\": [\"OMA\",\"OrthoInspector\",\"BBH\",\"Internet memes\"], \"human\": [17,24,17,50], \\\n",
    "                                \"banana\":[14,21,10,50]})\n",
    "df = pd.melt(summary_df, id_vars=\"method\", value_vars=[\"human\", \"banana\"])\n",
    "df.rename({\"variable\":\"species\", \"value\": \"% of genome shared with the other\"}, axis=1, inplace=True)\n",
    "\n",
    "sns.barplot(x=\"method\", y=\"% of genome shared with the other\", data=df, hue=\"species\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.751Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "myplot = sns.barplot(x=\"method\", y=\"human\", data=summary_df)\n",
    "for p in myplot.patches:\n",
    "    myplot.annotate(format(p.get_height(), '.0f'), \n",
    "                   (p.get_x() + p.get_width() / 2., p.get_height()), \n",
    "                   ha = 'center', va = 'center', \n",
    "                   xytext = (0, 9), \n",
    "                   size = 16,\n",
    "                   textcoords = 'offset points')\n",
    "plt.title(\"Percentage of Human Genes Orthologous to Banana\", size=26)\n",
    "plt.ylabel(\"%\")\n",
    "plt.savefig('./Percentage_human_banana_orthologs_barplot.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "# Bonus! GO enrichment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "So what are the functions of the human-banana orthologs? We can find out by doing a Gene Ontology (GO) enrichment analysis. \n",
    "\n",
    "To do this, I will get the canonical ID of each human gene with an ortholog and use this as the \"study set\" for a GOEA using the [panther tool](http://pantherdb.org/webservices/go/overrep.jsp).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.755Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "#write file for upload to go enrichment tool\n",
    "canonical_ids = human_genome[human_genome['omaid'].isin(human_genes_with_orthologs)]['canonicalid']\n",
    "\n",
    "with open(\"human_orthologs_for_go_enrichment.txt\", \"w\") as outfile:\n",
    "    for humanid in canonical_ids:\n",
    "        outfile.write(humanid + \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Use the list of genes in the previous step to upload to the online go enrichment webserver. The following code works on the resulting table file, which I named \"GOanalysis.txt.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.758Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "go_df = pd.read_csv(\"./GOanalysis.txt\", sep=\"\\t\", skiprows=11)\n",
    "\n",
    "#do some manipulation to get the go term and the p-value of enriched terms\n",
    "go_df['go_term'] = go_df['GO biological process complete'].apply(lambda x: x.split(\"(\")[-1])\n",
    "go_df['go_term'] = go_df['go_term'].apply(lambda x: x.split(\")\")[0])\n",
    "go_df = go_df[go_df['go_term'].str.contains(\"GO\")]\n",
    "\n",
    "#print for revigo\n",
    "#print(go_df[['go_term','human_orthologs_for_go_enrichment.txt (FDR)']].to_string(index=False, justify=\"left\"))\n",
    "\n",
    "#write to file\n",
    "go_df[['go_term', 'human_orthologs_for_go_enrichment.txt (FDR)']].to_csv(\"./enrichedGOterms_bonferroni.tsv\", sep=\"\\t\", \\\n",
    "                                                                         index=False, header=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "source": [
    "Summarize and visualize the resulting list of enriched GO terms with GO-Figure! https://www.biorxiv.org/content/10.1101/2020.12.02.408534v1.full.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "start_time": "2020-12-08T16:13:20.762Z"
    },
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "python ~/Applications/GO-Figure/gofigure.py -i ~/DessimozRepos/blogpost-code/enrichedGOterms.tsv \\\n",
    "-o gofigure/ --legend full --title \"Enriched GO terms for human-banana orthologs\" --cluster_labels description-numbered"
   ]
  }
 ],
 "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.0"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "102px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}