10b1 - MobilomeAnalysis_GenomeSizeProcessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cfd27947-0d7f-4960-9351-7b22db550e32",
   "metadata": {},
   "source": [
    "This part of the pipeline collects the genome sizes of all genomes in the genome set and statistically compares their sizes by rRNA cluster."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "077c164e-5fb5-4b16-a0a5-257226b26a01",
   "metadata": {},
   "source": [
    "### Paths and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e22bb5d-6755-47e6-9469-dfd7e56fcbe9",
   "metadata": {},
   "source": [
    "#### Pipeline input folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd0c9318-45b8-451c-8199-b5152283f244",
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata = \"./genomes_metadata\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ac0eea5-04dc-4cf1-acbd-0af35d634afc",
   "metadata": {},
   "source": [
    "#### Pipeline output folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4edcc105-2f52-43b7-8826-d5fa208d9272",
   "metadata": {},
   "outputs": [],
   "source": [
    "task_root = \"./10-MGEs/genome_sizes\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3a0cf59-aa31-4c78-9f7b-f6cce69034ff",
   "metadata": {},
   "source": [
    "#### Tool pointers and parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87fc1854-5301-4cf8-8463-11090f977531",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "fdcd7105-81a9-494a-9695-b02cdf11131b",
   "metadata": {},
   "source": [
    "#### Libraries and other setups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffbe2bf8-fa41-436b-936a-b295be811f2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import itertools as it\n",
    "import scipy.stats as sts\n",
    "import numpy as np\n",
    "from statannotations.Annotator import Annotator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "273d9ab4-5e97-46ed-8270-d818f5a83c5b",
   "metadata": {},
   "source": [
    "Let's pick some colours from the HUSL colour palette for the plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "985ea8d1-38fd-4ee8-bf18-ecc06aa08bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_palette = sns.husl_palette()\n",
    "custom_palette = [custom_palette[0], custom_palette[2], custom_palette[4], custom_palette[5]]\n",
    "custom_palette"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "478b14bc-4bc6-4023-8440-a0d409289764",
   "metadata": {},
   "source": [
    "### Reading metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3adfecae-e4c5-44e3-b3d3-0e1a2eb6ee0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_table(metadata, usecols = [0,2,4,13], sep = \"\\t\", converters = {4: lambda x: int(x)/1000000})\n",
    "data = data.rename(columns = {'Size (bp)': 'Size (Mb)'})\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0c64529-7899-429e-a498-f05c15098da9",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(data.groupby(['Taxonomic_cluster', 'Assumed_lifestyle']).count())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d53f91a1-e71b-4832-9844-ce7c6491331d",
   "metadata": {},
   "source": [
    "### Stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fc8556c-e1fa-4c9b-b922-987ee17e1936",
   "metadata": {},
   "source": [
    "#### Barplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e71021fd-5c1a-4494-86ed-2f72fe90f836",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (5,2))\n",
    "ax = sns.barplot(ax = ax, data = data, estimator = \"mean\", errorbar = \"se\",\n",
    "                 x = \"Size (Mb)\", y = \"Taxonomic_cluster\", palette = custom_palette,\n",
    "                 width = 0.9, orient = \"h\")\n",
    "plt.xlabel('Genome size (Mb)')\n",
    "plt.ylabel('rRNA cluster')\n",
    "plt.title('Genome sizes')\n",
    "plt.savefig(task_root + \"/\" + \"av_genome_sizes.svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e098fd4d-ff70-48b9-9155-5a49a0e66474",
   "metadata": {},
   "source": [
    "#### Violinplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bc4b631-cb8e-426d-ad07-13c10c9580c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (5,3))\n",
    "ax = sns.violinplot(ax = ax, data = data, x = 'Size (Mb)', y = 'Taxonomic_cluster', palette = custom_palette, orient = 'h', cut = 0)\n",
    "plt.xlabel('Genome size (Mb)')\n",
    "plt.ylabel('rRNA cluster')\n",
    "plt.title('Genome sizes')\n",
    "\n",
    "# Add statistical significance marks\n",
    "pairs = list(it.combinations(data['Taxonomic_cluster'].unique(), 2)) # get all rRNA cluster pairs\n",
    "annotator = Annotator(ax = ax, pairs = pairs, data = data, x = 'Size (Mb)', y = 'Taxonomic_cluster', orient = 'h', cut = 0)\n",
    "annotator.configure(test = 'Mann-Whitney', text_format = 'star', loc = 'inside')\n",
    "annotator.apply_and_annotate()\n",
    "\n",
    "plt.savefig(task_root + \"/\" + 'genome_sizes.svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "626cef6d-73cb-4127-a58c-ff90e46607f9",
   "metadata": {},
   "source": [
    "#### Stats values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d806ce6-4e38-4846-9a50-0c9929f85338",
   "metadata": {},
   "source": [
    "Get all genome sizes and append them to the set of their rRNA cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baecbf96-08d9-4d9f-9060-894fd05cca6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "size = data[['Size (Mb)', 'Taxonomic_cluster']].to_dict(orient = 'list')\n",
    "size_data = list(zip(*size.values()))\n",
    "size_stats = {}\n",
    "for record in size_data:\n",
    "    try:\n",
    "        size_stats[record[1]].append(record[0])\n",
    "    except KeyError:\n",
    "        size_stats[record[1]] = [record[0]]\n",
    "size_stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e12e6210-8be0-445a-b307-a31c79cdc820",
   "metadata": {},
   "source": [
    "General stats and Mann-Whitney U-testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec9ce1c7-ead7-4bf1-81c8-542813da2a70",
   "metadata": {},
   "outputs": [],
   "source": [
    "[(i, [np.mean(j), np.std(j)]) for i,j in size_stats.items()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5feb94ba-ea93-4939-9fe6-a7b71c293ab1",
   "metadata": {},
   "outputs": [],
   "source": [
    "tests = it.combinations(size_stats.keys(), 2) # get all rRNA cluster pairs\n",
    "for comb in tests:\n",
    "    print(str(comb) + ': ' + \n",
    "          str(sts.mannwhitneyu(size_stats[comb[0]], \n",
    "                               size_stats[comb[1]])\n",
    "              [1])\n",
    "         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c1d0866-680f-4dcd-94a1-9eca28e9ce60",
   "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
}