10b3 - MobilomeAnalysis_PseudogenesProcessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e1c54f2f-bb47-4daf-99a9-897cb8d612b4",
   "metadata": {},
   "source": [
    "This part of the pipeline processes the raw Pseudofinder output and statistically compares the normalised pseudogene counts by rRNA cluster."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb305754-9ec9-475e-b61c-7e201aa4f5b4",
   "metadata": {},
   "source": [
    "### Paths and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63302c23-06b8-41c2-aa15-10a2ac191800",
   "metadata": {},
   "source": [
    "#### Pipeline input folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42447fe5-7531-43ef-91c1-47dc1d3ee561",
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata = \"./genomes_metadata\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d98356c2-9941-414c-beaa-4536446e3c93",
   "metadata": {},
   "source": [
    "#### Pipeline output folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77f247fb-d398-419f-af52-1fb617e5dafa",
   "metadata": {},
   "outputs": [],
   "source": [
    "task_root = \"./10-MGEs/pseudogenes\"\n",
    "output_folder = task_root+\"/output\"\n",
    "results_folder = task_root+\"/processed_output\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e549094e-9f5c-481d-a75c-a82d74fe9157",
   "metadata": {},
   "source": [
    "#### Tool pointers and parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c5059be-13a7-4243-8da9-7562680b4f1d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "228b435e-e759-4dce-9967-217a9adcc8a7",
   "metadata": {},
   "source": [
    "#### Libraries and other setups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70d0125b-e5d6-4372-9bf0-9826f70c4d98",
   "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": "code",
   "execution_count": null,
   "id": "d5a4488f-5f45-4a9a-b2c8-9f88bb36af93",
   "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 = {'1': '#ee6677',\n",
    "                  '4': '#4477aa',\n",
    "                  '14a': '#228833',\n",
    "                  '14b': '#ccbb44'}\n",
    "custom_palette"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebff6dcd-5cc8-4311-a334-99c724b71706",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs(results_folder, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d86c84f-5098-4770-9680-e3ef76629cc6",
   "metadata": {},
   "source": [
    "## Reading input files"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92bc04e2-7ba3-422a-99ec-54dd9cf766b2",
   "metadata": {},
   "source": [
    "### Parsing pseudofinder results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef794b67-2b4d-413d-8498-3b61b6a1f837",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.makedirs(results_folder, exist_ok=True)\n",
    "result_dirs = os.listdir(output_folder)\n",
    "hits = []\n",
    "# Pseudogene hits are listed in the \"*_pseudos.fasta\" file, so counting the records ('>') in those files\n",
    "for dir in result_dirs:\n",
    "    dir_conts = os.listdir(output_folder + '/' + dir)\n",
    "    pseudo_fastas = [f for f in dir_conts if '_pseudos.fasta' in f]\n",
    "    for pfa in pseudo_fastas:\n",
    "        with open(output_folder + '/' + dir + '/' + pfa, \"r\") as handle:\n",
    "            cont = handle.read()\n",
    "            counts = cont.count('>')\n",
    "        record = {'assembly_ID': dir, 'counts': counts}\n",
    "        hits.append(record)\n",
    "hits = pd.DataFrame(hits)\n",
    "hits"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df8a048-ddf0-4f74-ad5a-3db88b16b2d2",
   "metadata": {},
   "source": [
    "### Loading the cluster annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73aa898f-c33a-4a03-bbd6-09d1b6797f04",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cluster_annotations_0 = pd.read_table(metadata, sep = '\\t', usecols = [1,2])\n",
    "cluster_annotations_0.columns = ['assemblyID', 'cluster']\n",
    "cluster_annotations = cluster_annotations_0.to_dict(orient = 'list')\n",
    "cluster_annotations = dict(zip(*cluster_annotations.values()))\n",
    "cluster_annotations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e135c41b-7b45-4ccf-9988-36243b7addae",
   "metadata": {},
   "source": [
    "### Reading genome sizes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54c695e3-db1a-4076-a6dd-ef1b531aa7e2",
   "metadata": {},
   "source": [
    "Necessary for normalising the pseudogene counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36eece84-9e98-4fcb-b252-cf1fae8fd89d",
   "metadata": {},
   "outputs": [],
   "source": [
    "genome_sizes_0 = pd.read_table(metadata, sep = '\\t', usecols = [1,4])\n",
    "genome_sizes_0.columns = ['assemblyID', 'size']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64d9ee3d-0102-40f1-acd1-2a9b85474ee8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "genome_sizes = genome_sizes_0.to_dict(orient = 'list')\n",
    "genome_sizes = dict(zip(*genome_sizes.values()))\n",
    "genome_sizes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc986ef1-5d90-41af-8247-3803dbffbdc2",
   "metadata": {},
   "source": [
    "### Adding metadata columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "545d79fc-8694-4f5c-9bb4-c4daf02595ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "hits['cluster'] = hits['assembly_ID'].apply(lambda x: cluster_annotations[x])\n",
    "hits['size'] = hits['assembly_ID'].apply(lambda x: genome_sizes[x])\n",
    "hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acce77ea-d680-4ea0-9ac7-b22f3dac2649",
   "metadata": {},
   "outputs": [],
   "source": [
    "hits.to_csv(results_folder + \"/counts\", sep = '\\t', index = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68bb373d-d8c1-4559-bcaf-dd72234d0fa7",
   "metadata": {},
   "source": [
    "### General count plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a404119e-bd9e-4ad6-85c8-20ee1b74ef56",
   "metadata": {},
   "source": [
    "Normalise by genome size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cd4aa1d-22e6-40d1-ad6d-4c2c47a967d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_counts = hits.rename(columns = {'counts': 'No. pseudogenes'})\n",
    "cluster_counts['Norm. no. pseudogenes'] = cluster_counts['No. pseudogenes']/cluster_counts['size']*1000000\n",
    "cluster_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b1269b5-c262-421b-8790-7be5063eeca9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_counts.to_csv(results_folder + \"/\" + 'counts_per_assembly_cluster', sep = '\\t', index = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60df7a1e-5e7f-45d1-b9fb-52f03a39504d",
   "metadata": {},
   "source": [
    "#### Barplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd8859e7-edd0-4030-9b56-948957efd10d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (5,2))\n",
    "ax = sns.barplot(ax = ax, data = cluster_counts, estimator = \"mean\", errorbar = \"se\",\n",
    "                 x = \"Norm. no. pseudogenes\", y = \"cluster\", palette = custom_palette,\n",
    "                 width = 0.9, orient = \"h\", order = ['1','14a','4','14b'])\n",
    "plt.xlabel('Avg. norm. no. pseudogenes')\n",
    "plt.ylabel('rRNA cluster')\n",
    "plt.title('Pseudogenes')\n",
    "plt.savefig(results_folder + \"/\" + \"av_counts_pseudogenes_cluster_bar.svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2d46245-2b57-46f0-88d4-e421dd018c77",
   "metadata": {},
   "source": [
    "#### Violinplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de635de0-7b80-42b4-a88a-9a69447c8c3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (5,3))\n",
    "ax = sns.violinplot(ax = ax, data = cluster_counts, x = 'Norm. no. pseudogenes', y = 'cluster', \n",
    "                    palette = custom_palette, orient = 'h', cut = 0, order = ['1','14a','4','14b'])\n",
    "plt.xlabel('Norm. no. pseudogenes')\n",
    "plt.ylabel('rRNA cluster')\n",
    "plt.title('Pseudogenes')\n",
    "\n",
    "# Add statistical significance markers\n",
    "pairs = list(it.combinations(cluster_counts['cluster'].unique(), 2))\n",
    "annotator = Annotator(ax = ax, pairs = pairs, data = cluster_counts, x = 'Norm. no. pseudogenes', y = 'cluster', orient = 'h', cut = 0,\n",
    "                      order = ['1','14a','4','14b'])\n",
    "annotator.configure(test = 'Mann-Whitney', text_format = 'star', loc = 'inside')\n",
    "annotator.apply_and_annotate()\n",
    "\n",
    "plt.savefig(results_folder + \"/\" + 'counts_pseudogenes_cluster_violin.svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faf731b5-c8e6-4045-a4dd-24fc5bfa0a11",
   "metadata": {},
   "source": [
    "#### Exact stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3126789d-b7f3-409f-9f78-ba67e1ff2c99",
   "metadata": {},
   "source": [
    "Getting all the counts grouped by rRNA cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bb15c2a-12b9-40a3-80f6-cf1aa0cd4c38",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cluster_counts_stats = cluster_counts[['Norm. no. pseudogenes', 'cluster']].to_dict(orient = 'list')\n",
    "cluster_counts_stats = list(zip(*cluster_counts_stats.values()))\n",
    "counts_stats = {}\n",
    "for record in cluster_counts_stats:\n",
    "    try:\n",
    "        counts_stats[record[1]].append(record[0])\n",
    "    except KeyError:\n",
    "        counts_stats[record[1]] = [record[0]]\n",
    "counts_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cbe8944-6f20-4d50-94cc-b262ef7c0492",
   "metadata": {},
   "outputs": [],
   "source": [
    "[(i, [np.mean(j), np.std(j)]) for i,j in counts_stats.items()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6665dad-88ec-42db-9899-ac8b6e1517ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "tests = it.combinations(counts_stats.keys(), 2)\n",
    "for comb in tests:\n",
    "    print(str(comb) + ': ' + \n",
    "          str(sts.mannwhitneyu(counts_stats[comb[0]], \n",
    "                           counts_stats[comb[1]])\n",
    "              [1])\n",
    "         )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f7a301-0fbe-4cad-9adb-4421506bf081",
   "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
}