03 - fastANI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2ae56af3-7ed0-4d0b-955d-fb536e0a6ffd",
   "metadata": {},
   "source": [
    "This part of the pipeline carries out a fastANI analysis and immediately processes the results into a clustermap. It also collects the assemblies that had an ANI similarity higher the recommended threshold of 95% to delineate species."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93b405a3-6acb-4065-884d-b6fa9d6492d2",
   "metadata": {},
   "source": [
    "### Checking dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1689eb8-0591-48b3-9b1e-5c461d7500db",
   "metadata": {},
   "outputs": [],
   "source": [
    "!fastANI --version"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e491e6e9-c341-4f9a-aece-6f3b2bfdce3d",
   "metadata": {},
   "source": [
    "### Paths and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed9022d7-c198-4d64-97f8-ff79e9525586",
   "metadata": {},
   "source": [
    "#### Pipeline input folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97785fec-da2e-4282-bc4e-17ea9d4d9d43",
   "metadata": {},
   "outputs": [],
   "source": [
    "genomes = \"./02-QC/data/genomes\"\n",
    "metadata = \"./genomes_metadata\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1151b578-60bc-4c29-9f01-6254d50fc429",
   "metadata": {},
   "source": [
    "#### Pipeline output folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6a641aa-d67e-4772-baa8-1d8492c76721",
   "metadata": {},
   "outputs": [],
   "source": [
    "task_root = \"./03-fastANI\"\n",
    "\n",
    "!mkdir -p $task_root"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "378b7f8b-ae90-44de-aba1-a07b4c72844f",
   "metadata": {},
   "source": [
    "#### Tool pointers and parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4281e12-76ed-48b0-b3de-0543fe904f5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "runner_script = \"./utils/run_fastANI.sh\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e4e27c8-7fe9-4081-8793-3acc82b87ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_color_palette = {'b': '#36489a',\n",
    "                    'g': '#228833',\n",
    "                    'r': '#ee6677',\n",
    "                    'k': '#ccbb44'}\n",
    "cp = lambda x: my_color_palette[x]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec12d706-00f9-41d5-9066-ba44aef10885",
   "metadata": {},
   "source": [
    "### Run the fastANI analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b447cf0a-8e73-4e8c-8bf7-3050ae3cf964",
   "metadata": {},
   "outputs": [],
   "source": [
    "!bash $runner_script $genomes $task_root"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4b97d54-57aa-4e8a-a9f5-61c8c895c30d",
   "metadata": {},
   "source": [
    "### Process the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6ca30cc-b7ce-4d2b-8a9d-b353f4324f68",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from scipy.cluster import hierarchy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75663c44-41ee-4ffb-a47c-91d4a53bcc8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fastani_file = task_root+\"/out\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a6fb426-e416-44c3-bae8-9c4259de3bf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fastani = pd.read_table(fastani_file, usecols = [0,1,2], header = None, names = ['Genome 1', 'Genome 2', 'ANI'])\n",
    "fastani['Genome 1'] = fastani['Genome 1'].apply(lambda x: x.split('/')[-1][:-4])\n",
    "fastani['Genome 2'] = fastani['Genome 2'].apply(lambda x: x.split('/')[-1][:-4])\n",
    "fastani = fastani.pivot(index = 'Genome 1', columns = 'Genome 2', values = 'ANI').fillna(70)\n",
    "fastani = (fastani + fastani.T)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b564cc7-24ae-4966-a443-fdd4588ea649",
   "metadata": {},
   "outputs": [],
   "source": [
    "fastani"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1d6aa83-e804-462f-84a0-857dd80e7e27",
   "metadata": {},
   "outputs": [],
   "source": [
    "fastani.to_csv(task_root+'/processed_out.tsv', sep = '\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19b196df-307c-47a1-920f-761c5b6f9468",
   "metadata": {},
   "source": [
    "### Define rRNA group colouring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc216ae9-0adb-422c-ba10-281ad9fb0cfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "groups = pd.read_table(metadata, usecols = [1,2]).set_index('Genome_accession').squeeze()\n",
    "group_colours = groups.map(dict(zip(['1','4','14a','14b'],[cp(c) for c in 'rbgk'])))\n",
    "group_colours.name = \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed678bce-01fe-44f0-987c-d1e24c1b3ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "group_colours"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f31c833d-2e38-40cd-a68f-5619bc70dd0f",
   "metadata": {},
   "source": [
    "### Making ANI clustermap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a8822d1-d586-403c-8c8b-ab6fdbc91b82",
   "metadata": {},
   "outputs": [],
   "source": [
    "col_linkage = hierarchy.linkage(fastani, method = \"average\", metric = \"euclidean\", optimal_ordering = True)\n",
    "row_linkage = hierarchy.linkage(fastani.T, method = \"average\", metric = \"euclidean\", optimal_ordering = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe468c96-c225-457c-9eba-c82d8183ac67",
   "metadata": {},
   "outputs": [],
   "source": [
    "#sns.set(font_scale=1)\n",
    "fig = sns.clustermap(fastani, col_linkage = col_linkage, row_linkage = row_linkage, col_cluster = True, row_cluster = True,\n",
    "                     row_colors = group_colours, cmap=sns.color_palette('magma_r', as_cmap=True), dendrogram_ratio = 0.15,\n",
    "                     figsize = (8,8), cbar_kws = {'orientation': 'horizontal'}, cbar_pos = (0.01, 0.89, 0.16, 0.05),\n",
    "                     xticklabels = False, yticklabels = False)\n",
    "fig.ax_cbar.set_title('ANI (%)')\n",
    "ax = fig.ax_heatmap\n",
    "ax.set(xlabel = None, ylabel = None)\n",
    "ax.text(65,65, \"ANI group 1\", fontsize = 14, ha = \"center\", va = \"center\")\n",
    "ax.text(435,435, \"ANI group 2\", fontsize = 14, ha = \"center\", va = \"center\")\n",
    "ax.text(340,340, \"ANI group 3\", fontsize = 14, ha = \"center\", va = \"center\")\n",
    "ax.text(215,215, \"ANI group 4\", fontsize = 14, ha = \"center\", va = \"center\")\n",
    "plt.savefig(task_root+\"/ANI_heatmap.svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f62e7747-6560-450e-8f4c-ba7c7bb3528c",
   "metadata": {},
   "source": [
    "### Super-threshold ANIs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6ab6711-f972-441e-9f56-7646ac2c8b16",
   "metadata": {},
   "outputs": [],
   "source": [
    "fastani_molten = fastani.where(np.triu(np.ones(fastani.shape), k = 1).astype(bool)).stack().reset_index().rename(columns = {0: 'ANI'})\n",
    "fastani_molten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bc27748-1a7f-423b-a1ad-17b72d2f2333",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Super-threshold ANIs omitting identical genomes\n",
    "super_anis = fastani_molten[(fastani_molten['ANI'] >= 95) & (fastani_molten['ANI'] < 100)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "437f656a-e11d-4303-811a-07fdc709edbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "super_anis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b94ab6e-9e5f-42fb-acde-0806d37e3821",
   "metadata": {},
   "outputs": [],
   "source": [
    "super_anis.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24e8e0c5-7719-4276-a309-a0be5072e66b",
   "metadata": {},
   "outputs": [],
   "source": [
    "super_anis.to_csv(task_root+'/super_threshold', sep = '\\t', index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0505ad77-7403-4612-ba97-b83a2acef0ef",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}