23_FigS8_prokka_feat.R.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1227e32d-12e8-45f9-b6c7-5ccf80b0eb75",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import yaml\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "466530c3-a6d9-434a-b40c-dd81f1ea9d19",
   "metadata": {},
   "source": [
    "## File configurations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0126f790-ba27-43ef-b009-e318dd837d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"config.yaml\", \"r\") as f:\n",
    "    notebook_configuration = yaml.safe_load(f)\n",
    "notebook_configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4dd12fa-3e81-4548-ac7e-4f04afb9125c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bgcflow_dir = Path(notebook_configuration[\"bgcflow_dir\"])\n",
    "project_name = \"mq_saccharopolyspora\"\n",
    "report_dir = bgcflow_dir / f\"data/processed/{project_name}\"\n",
    "FIGURE = \"Figure_S8\"\n",
    "FIGURE_TREE = \"Figure_3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ff1eb93-9238-4e24-90ce-5df1d31551c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "R_feature_data = f\"assets/tables/{FIGURE}_features_stacked.csv\"\n",
    "R_bases_data = f\"assets/tables/{FIGURE}_bases_stacked.csv\"\n",
    "R_contigs_data = f\"assets/tables/{FIGURE}_contigs_stacked.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "660d2b25-bbe1-48d1-a287-4ff1c9576b0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = []\n",
    "for item in (report_dir / \"genbank\").glob(\"*.txt\"):\n",
    "    df = pd.read_csv(item, sep=\":\", skiprows=1, header=None, index_col=0).rename(columns={1:item.stem}).T\n",
    "    server_path = \"<a href='{{ project().file_server() }}/genbank/\"\n",
    "    df.loc[item.stem, \"gbk file\"] = server_path + f\"{item.stem}.gbk' target='_blank''>{item.stem}.gbk</a>\"\n",
    "    df.loc[item.stem, \"CDS table\"] = server_path + f\"{item.stem}.tsv' target='_blank''>{item.stem}.tsv</a>\"\n",
    "    summary.append(df)\n",
    "df = pd.concat(summary).sort_index()\n",
    "df_table = df.loc[:, ['contigs', 'bases', 'CDS', 'rRNA', 'tRNA', 'tmRNA', 'repeat_region']].fillna(0)\n",
    "\n",
    "result = {}\n",
    "for gbk_file in (report_dir / \"genbank\").glob(\"*.tsv\"):\n",
    "    gbk_file = Path(gbk_file)\n",
    "    gbk_table = pd.read_csv(gbk_file, sep=\"\\t\")\n",
    "    genome_id = gbk_file.stem\n",
    "    gbk_table['genome_id'] = genome_id\n",
    "    unknown_cds = gbk_table['product'].value_counts()['hypothetical protein']\n",
    "    known_cds = gbk_table.shape[0] - unknown_cds\n",
    "    result[genome_id] = {\"CDS (Unknown)\" : unknown_cds, \n",
    "                         \"CDS (Known)\" : known_cds}\n",
    "df_table = pd.concat([df_table, pd.DataFrame.from_dict(result).T], axis=1)\n",
    "\n",
    "df_table.to_csv(f\"assets/tables/{FIGURE}_prokka.csv\", index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45afa348-c66d-4bbb-b28c-dc13bc74f64f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_features = df_table.loc[:, ['rRNA', 'tRNA', 'tmRNA', 'repeat_region', \"CDS (Unknown)\", \"CDS (Known)\"]].reset_index(drop=False).melt(id_vars='index').rename(columns={\"index\":\"genome_id\", 'variable':\"category\", \"value\":\"count\"})\n",
    "df_features.to_csv(R_feature_data, index=False)\n",
    "\n",
    "df_bases = df_table.loc[:, ['bases']].reset_index(drop=False).melt(id_vars='index').rename(columns={\"index\":\"genome_id\", 'variable':\"category\", \"value\":\"count\"})\n",
    "df_bases.to_csv(R_bases_data, index=False)\n",
    "\n",
    "df_contigs = df_table.loc[:, ['contigs']].reset_index(drop=False).melt(id_vars='index').rename(columns={\"index\":\"genome_id\", 'variable':\"category\", \"value\":\"count\"})\n",
    "df_contigs.to_csv(R_contigs_data, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3973e1e2-fe32-4956-bf9d-c2a661aa0e56",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "library(\"treeio\")\n",
    "library(\"ggtree\")\n",
    "library(\"tidyverse\")\n",
    "library(\"ggstar\")\n",
    "library(\"ggnewscale\")\n",
    "library(\"ggtreeExtra\")\n",
    "library(\"phangorn\")\n",
    "library(\"svglite\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edc46555-9553-4fe2-9fc3-9c1c772afc93",
   "metadata": {},
   "outputs": [],
   "source": [
    "R_tree = str(report_dir / \"automlst_wrapper/final_corrected.newick\")\n",
    "R_data = f\"assets/tables/{FIGURE_TREE}c_ggtree.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1d183ce-c791-42a6-8d9f-b1ff2fc38788",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R  -w 1200 -h 1200 -i R_tree -i R_data\n",
    "tree <- read.tree(R_tree)\n",
    "#data <- read.csv(\"../automlst_wrapper/df_genomes_tree.csv\")\n",
    "data <- read.csv(R_data)\n",
    "\n",
    "# midpoint root\n",
    "tree <- phangorn::midpoint(tree)\n",
    "tree <- ladderize(reorder(tree))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc78d391-4510-4d01-89a4-7673a24a9e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R  -w 600 -h 300\n",
    "\n",
    "p <- ggtree(tree)\n",
    "p <- p %<+% data\n",
    "#p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7f49830-398f-4a64-9b84-fc860a3028d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "p2 = p + geom_tiplab(aes(label=tip_label, fill=phylogroup),\n",
    "                     size=2.6, hjust=-0.03, family='sans',\n",
    "                    linetype = \"dotted\", linesize = 1) + # size of label border) +\n",
    "        hexpand(5) +\n",
    "        geom_tippoint(size=2.6, alpha=0.8, aes(colour=phylogroup, stroke=1, shape=sequence_quality)) + \n",
    "        scale_color_manual(labels=c(\"P1\", \"P2\", \"P3\", \"P4\", \"P5\", \"P6\", \"P7\", \"P8\"),\n",
    "                           values=c('#264653','#e9c46a','#808080', '#808080','#f4a261','#808080', '#e76f51', '#2a9d8f'))\n",
    "p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be3c499f-f49d-462a-8220-4bba400ca2c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i R_feature_data\n",
    "# add feature\n",
    "feature_data <- read.csv(R_feature_data)\n",
    "p3 = p2 + new_scale_fill() + geom_facet(panel='Feature count (log scale)',\n",
    "                     geom=geom_bar,\n",
    "                     data=feature_data,\n",
    "                     mapping=aes(x=count, fill=category),\n",
    "                     position=\"stack\",\n",
    "                     orientation=\"y\", \n",
    "                     stat=\"identity\") +\n",
    "        scale_x_log10() + theme(text=element_text(size=10, family=\"sans\")) + theme_bw() #+\n",
    "        #scale_fill_manual(labels=c(\"bgcs_on_contig_edge\", \"complete_bgcs\"),\n",
    "        #                  values=c(\"#cc2936\", \"#08415c\")) \n",
    "p3 = facet_widths(p3, widths = c(3, 2))\n",
    "p3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e3d5f74-3b11-4b00-8aba-c16cdf40ccdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "R_fig3c_svg = f\"assets/figures/{FIGURE}/{FIGURE}.svg\"\n",
    "R_fig3c_pdf = f\"assets/figures/{FIGURE}/{FIGURE}.pdf\"\n",
    "Path(R_fig3c_svg).parent.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29857d29-57af-4067-9e98-f9a34827ad83",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i R_fig3c_svg -i R_fig3c_pdf\n",
    "ggsave(file=R_fig3c_svg , plot=p3, device=svglite, width=4000, height=1800, units=\"px\")\n",
    "ggsave(plot=p3, width=4000, height=1600, units=\"px\", dpi=300, filename=R_fig3c_pdf, useDingbats=FALSE)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "857b42e1-d791-4cd1-88b1-e281b66fa471",
   "metadata": {},
   "source": [
    "%%R -i R_bases_data\n",
    "# add feature\n",
    "bases_data <- read.csv(R_bases_data)\n",
    "p4 = p3 + new_scale_fill() + geom_facet(panel='Bases',\n",
    "                     geom=geom_bar,\n",
    "                     data=bases_data,\n",
    "                     mapping=aes(x=count, fill=category),\n",
    "                     position=\"stack\",\n",
    "                     orientation=\"y\", \n",
    "                     stat=\"identity\") #+\n",
    "        #scale_fill_manual(labels=c(\"bgcs_on_contig_edge\", \"complete_bgcs\"),\n",
    "        #                  values=c(\"#cc2936\", \"#08415c\")) \n",
    "p4"
   ]
  },
  {
   "cell_type": "raw",
   "id": "d386143b-42e1-4737-901b-20c57d131dbe",
   "metadata": {},
   "source": [
    "%%R -i R_contigs_data\n",
    "# add feature\n",
    "contigs_data <- read.csv(R_contigs_data)\n",
    "p5 = p4 + new_scale_fill() + geom_facet(panel='Contigs',\n",
    "                     geom=geom_bar,\n",
    "                     data=contigs_data,\n",
    "                     mapping=aes(x=count, fill=category),\n",
    "                     position=\"stack\",\n",
    "                     orientation=\"y\", \n",
    "                     stat=\"identity\") #+\n",
    "        #scale_fill_manual(labels=c(\"bgcs_on_contig_edge\", \"complete_bgcs\"),\n",
    "        #                  values=c(\"#cc2936\", \"#08415c\")) \n",
    "p5"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}