14_Fig6a_Lanthipeptide_tree.R.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c1782be-f1f6-4dc2-928f-3685e7f799b8",
   "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": "f9b2d90a-661f-4f5b-87f0-ace6484ff413",
   "metadata": {},
   "source": [
    "## File configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a0809be-da58-4c58-9e35-779bbcef9656",
   "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": "a6dc590a-66f8-4ef0-8f3f-7c74f5fc8883",
   "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_6\"\n",
    "FIGURE_TREE = \"Figure_3\"\n",
    "cutoff = \"0.30\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79d7c5e9-1c1c-44c5-8f8d-321ebf18c29a",
   "metadata": {},
   "outputs": [],
   "source": [
    "bigscape_dir = report_dir / \"bigscape/for_cytoscape_antismash_6.1.1/\"\n",
    "bigscape_mapping = [i for i in bigscape_dir.glob(f\"*df_families_{cutoff}.csv\")]\n",
    "bigscape_cluster = [i for i in bigscape_dir.glob(f\"*df_clusters_{cutoff}.csv\")]\n",
    "bigscape_presence = [i for i in bigscape_dir.glob(f\"*df_family_presence_{cutoff}.csv\")]\n",
    "assert len(bigscape_mapping) == 1\n",
    "assert len(bigscape_cluster) == 1\n",
    "assert len(bigscape_presence) == 1\n",
    "df_mapping = pd.read_csv(bigscape_mapping[0], index_col=0)\n",
    "df_cluster = pd.read_csv(bigscape_cluster[0], index_col=0)\n",
    "df_presence = pd.read_csv(bigscape_presence[0], index_col=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab83d4c9-a8b9-44e9-adf2-f91ff3fd5b93",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = list(df_cluster.loc[[i for i in df_cluster.index if \"lanthipeptide\" in df_cluster.loc[i, 'product']], f\"fam_id_{cutoff}\"].unique())\n",
    "fam_mapping = df_mapping.loc[query, \"fam_name\"].to_dict()\n",
    "fam_mapping = {str(k):v for k,v in fam_mapping.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f503e9d-d155-4f82-9bba-e3e57a0a9255",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = [str(i) for i in query]\n",
    "df_presence = df_presence.loc[:, query]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4616316a-a0a9-4817-bbf6-f683a46377ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_presence = df_presence.rename(columns=fam_mapping)\n",
    "df_presence.to_csv(f\"assets/tables/{FIGURE}a_presence_absence_lanthipeptide.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48c1e3b8-879d-4ef7-b0e5-754adda896e2",
   "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": "cc80c1cc-f2ab-4adf-a793-f844cb96f46a",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_data = f\"assets/tables/{FIGURE_TREE}c_ggtree.csv\"\n",
    "tree_file = str((report_dir / \"automlst_wrapper/final_corrected.newick\").resolve())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0b37183-6554-4cd4-8e94-9d2792d99809",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R  -w 800 -h 500 -i tree_data -i tree_file\n",
    "tree <- read.tree(tree_file)\n",
    "#data <- read.csv(\"../automlst_wrapper/df_genomes_tree.csv\")\n",
    "data <- read.csv(tree_data)\n",
    "\n",
    "# midpoint root\n",
    "tree <- phangorn::midpoint(tree)\n",
    "tree <- ladderize(reorder(tree))\n",
    "\n",
    "p <- ggtree(tree)\n",
    "data$tip_label2 <- data$genome_id\n",
    "p <- p %<+% data\n",
    "\n",
    "p2 <- p + geom_tiplab(aes(label=tip_label2, fill=phylogroup),\n",
    "                     size=2.6, hjust=-0.03, family='sans',\n",
    "                    linetype = \"dotted\", linesize = 1) + # size of label border) +\n",
    "        #hexpand(1.5) +\n",
    "        geom_tippoint(size=2.6, alpha=0.8, aes(colour=phylogroup, stroke=1)) + \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",
    "\n",
    "p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "054b0d27-7e57-48d6-a7fe-60902dc262d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_phylogroup = pd.read_csv(f\"assets/tables/{FIGURE_TREE}c_ggtree.csv\", index_col=\"genome_id\")\n",
    "#df_cluster = pd.read_csv(\"../bigscape/for_cytoscape_antismash_6.1.1/2023-01-30 16_15_43_df_clusters_0.30.csv\", index_col=0)\n",
    "for i in df_cluster.index:\n",
    "    genome_id = df_cluster.loc[i, \"genome_id\"]\n",
    "    df_cluster.loc[i, \"phylogroup\"] = df_phylogroup.loc[genome_id, \"phylogroup\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "143cbd2c-8ae0-4add-8541-5ba7fa0bdbba",
   "metadata": {},
   "outputs": [],
   "source": [
    "def map_value_to_color(value):\n",
    "    if value == 1:\n",
    "        return 'black'\n",
    "    else:\n",
    "        return 'white'\n",
    "\n",
    "df = pd.read_csv(f\"assets/tables/{FIGURE}a_presence_absence_lanthipeptide.csv\")\n",
    "df['xindex'] = df.index\n",
    "df = pd.melt(df, [\"genome_id\", \"xindex\"]).rename(columns={\"variable\":\"GCF\"})\n",
    "df['color'] = df['value'].apply(map_value_to_color)\n",
    "\n",
    "df.GCF\n",
    "\n",
    "df_cluster = df_cluster[df_cluster[f'fam_known_compounds_{cutoff}'].isin(df.GCF)]\n",
    "df_cluster['lanthipeptide_class'] = [\"\".join([x for x in i.split(\".\") if 'lanthipeptide' in x]) for i in df_cluster['product']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f6367de-7bc1-411d-ba5d-17f92fdcea05",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_category = {}\n",
    "for cluster in df_cluster[f\"fam_known_compounds_{cutoff}\"].unique():\n",
    "    subset = df_cluster[df_cluster[f\"fam_known_compounds_{cutoff}\"] == cluster]\n",
    "    category_raw = [\"\".join([x for x in i.split(\".\") if \"lanthipeptide\" in x]) for i in subset[\"product\"].to_list()]\n",
    "    category = list(set(category_raw))\n",
    "    if len(category) == 1:\n",
    "        cluster_category[cluster] = \"\".join(category)\n",
    "    else:\n",
    "        if \"\" in category:\n",
    "            category.remove(\"\")\n",
    "        cluster_category[cluster] = \" or \".join(category)\n",
    "cluster_category.values()\n",
    "df[\"class_category\"] = [cluster_category[i] for i in df.GCF]\n",
    "df[\"phylogroup\"] = [df_phylogroup.loc[i, \"phylogroup\"] for i in df.genome_id]\n",
    "\n",
    "for i in df.index:\n",
    "    if df.loc[i, 'value'] == 1:\n",
    "        df.loc[i, 'color'] = f'{df.loc[i, \"phylogroup\"]}_1'\n",
    "    else:\n",
    "        df.loc[i, 'color'] = \"\"\n",
    "\n",
    "        \n",
    "df = df[df.value == 1]\n",
    "xindex = {}\n",
    "for cat in sorted(set(cluster_category.values())):\n",
    "    subset = df[df[\"class_category\"] == cat].GCF.sort_values().drop_duplicates().reset_index(drop=True)\n",
    "    print(cat, len(xindex))\n",
    "    subset.index += (1+len(xindex))\n",
    "    subset = subset.to_dict()\n",
    "    xindex.update({v:k for k,v in subset.items()})\n",
    "xindex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e2a36a3-e22b-4136-bdb0-75a22df442bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# manually annotate as lanthipeptide class iii\n",
    "override = df[df.class_category.isin(['lanthipeptide-class-ii or lanthipeptide-class-iii', 'lanthipeptide-class-iii or lanthipeptide-class-iv'])].index\n",
    "df.loc[override, \"class_category\"] = 'lanthipeptide-class-iii'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1cc0a8b-4bd2-4cd4-b322-45d9595c3067",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.xindex = [xindex[i] for i in df.GCF]\n",
    "\n",
    "presence_absence_table = f\"assets/tables/{FIGURE}a_presence_absence_lanthipeptide_long.csv\"\n",
    "df.to_csv(presence_absence_table, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4a6356a-d723-4630-8b37-f49956391493",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 800 -h 400 -i tree_file -i tree_data -i presence_absence_table\n",
    "df <- read.csv(presence_absence_table,\n",
    "               stringsAsFactors = TRUE)\n",
    "\n",
    "tree <- read.tree(tree_file)\n",
    "#data <- read.csv(\"../automlst_wrapper/df_genomes_tree.csv\")\n",
    "data <- read.csv(tree_data)\n",
    "\n",
    "# midpoint root\n",
    "tree <- phangorn::midpoint(tree)\n",
    "tree <- ladderize(reorder(tree))\n",
    "\n",
    "p <- ggtree(tree)\n",
    "data$tip_label2 <- data$genome_id\n",
    "p <- p %<+% data\n",
    "\n",
    "p2 <- p + geom_tiplab(aes(label=tip_label2, fill=phylogroup),\n",
    "                     size=2.8, hjust=-0.03, family='sans',\n",
    "                    linetype = \"dotted\", linesize = 1) + # size of label border) +\n",
    "        hexpand(0.5) +\n",
    "        geom_tippoint(size=1, alpha=0.8, aes(colour=phylogroup, stroke=1)) + \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",
    "        new_scale_fill() \n",
    "\n",
    "p3 <- p2 #+ theme(legend.position='none')\n",
    "\n",
    "df_text <- df[!duplicated(df$xindex),]\n",
    "p4 <- facet_plot(p3, panel = \"BiG-SCAPE GCFs\", data = df, geom=geom_tile, mapping = aes(x=xindex, y=y, fill = class_category), \n",
    "                 color = \"white\", lwd = 1,\n",
    "                 linetype = 1) \n",
    "#p3 <- facet_plot(p3, panel = \"BiG-SCAPE GCFs\", data = df_text, geom=geom_text, mapping = aes(x=xindex, y = 0, label = GCF), size=2, angle = 90)\n",
    "p4 <- facet_widths(p4, widths = c(1, 1.5)) \n",
    "#p3 + xlim_expand(c(150, 150), \"BiG-SCAPE GCFs\")\n",
    "p4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eb7e502-a547-4470-acc9-30d756dce0ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "from svgutils.compose import *\n",
    "from svgutils.compose import Figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0ffa368-dff1-479e-86d6-e18381aa5b98",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = {i:n for n,i in enumerate(list(df.sort_values(by=\"xindex\").GCF.unique()))}\n",
    "df_region = pd.read_csv(report_dir / \"tables/df_regions_antismash_6.1.1.csv\", index_col=\"bgc_id\")\n",
    "for i in df_cluster.index:\n",
    "    similarity = df_region.loc[i, \"similarity\"]\n",
    "    #print(similarity, similarity > 0.7)\n",
    "    if similarity >= 0.5:\n",
    "        for c in [\"similarity\", \"most_similar_known_cluster_id\", \"most_similar_known_cluster_description\"]:\n",
    "            df_cluster.loc[i, c] = df_region.loc[i, c]\n",
    "df_cluster.loc[:, [\"fam_known_compounds_0.30\", \"most_similar_known_cluster_description\", \"similarity\"]]\n",
    "df_cluster[\"most_similar_known_cluster_description\"].unique()\n",
    "\n",
    "naming_dict = {'venezuelin' : \"Venezuelin\",\n",
    "               'labyrinthopeptin A2 / labyrinthopeptin A1 / labyrinthopeptin A3' : \"Labyrinthopeptin\",\n",
    "               'Ery-9 / Ery-6 / Ery-8 / Ery-7 / Ery-5 / Ery-4 / Ery-3': \"Erythreapeptin\",\n",
    "               'cinnamycin' : \"Kyamicin\", \n",
    "               'A83543A' : \"Spinosyn\", \n",
    "               'anantin C' : \"Anantin\", \n",
    "               'planosporicin' : \"Planosporicin\"}\n",
    "\n",
    "for i in df_cluster.index:\n",
    "    name = df_cluster.loc[i, \"most_similar_known_cluster_description\"]\n",
    "    if name in naming_dict.keys():\n",
    "        df_cluster.loc[i, \"most_similar_known_cluster_description\"] = naming_dict[name]\n",
    "\n",
    "label_mapping = {}\n",
    "        \n",
    "for gcf in df_cluster[\"fam_known_compounds_0.30\"].unique():\n",
    "    subset = df_cluster[df_cluster[\"fam_known_compounds_0.30\"] == gcf]\n",
    "    gcf_number = subset[\"fam_id_0.30\"].unique()\n",
    "    subset = subset.fillna(\"\")\n",
    "    label = list(subset[\"most_similar_known_cluster_description\"].unique())\n",
    "    label = \"\".join(label)\n",
    "    if label == \"ErythreapeptinAnantin\":\n",
    "        label = \"Erythreapeptin|Anantin\"\n",
    "    if label == \"\":\n",
    "        label = f\"GCF_{gcf_number[0]}\"\n",
    "    else:\n",
    "        label = f\"GCF_{gcf_number[0]} ({label})\"\n",
    "    label_mapping[gcf] = label\n",
    "\n",
    "labels = {v:label_mapping[k] for k,v in labels.items()}\n",
    "labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14bc0007-1b6c-46fe-bfac-5be46f45192a",
   "metadata": {},
   "outputs": [],
   "source": [
    "Path(f\"assets/figures/{FIGURE}/\").mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ce6807f-e3d2-453f-a598-7332ce0cf26a",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure_output_svg = f\"assets/figures/{FIGURE}/a.svg\"\n",
    "figure_output_pdf = f\"assets/figures/{FIGURE}/a.pdf\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6d93a86-0ed3-4c9f-a7c6-7105566efece",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 800 -h 400 -i presence_absence_table -i tree_data -i tree_file -i figure_output_svg -i figure_output_pdf \n",
    "df <- read.csv(presence_absence_table,\n",
    "               stringsAsFactors = TRUE)\n",
    "\n",
    "tree <- read.tree(tree_file)\n",
    "#data <- read.csv(\"../automlst_wrapper/df_genomes_tree.csv\")\n",
    "data <- read.csv(tree_data)\n",
    "\n",
    "# midpoint root\n",
    "tree <- phangorn::midpoint(tree)\n",
    "tree <- ladderize(reorder(tree))\n",
    "\n",
    "p <- ggtree(tree)\n",
    "data$tip_label2 <- data$genome_id\n",
    "p <- p %<+% data\n",
    "\n",
    "p2 <- p + geom_tiplab(aes(label=tip_label2, fill=phylogroup),\n",
    "                     size=2.7, hjust=-0.08, family='sans',\n",
    "                    linetype = \"dotted\", linesize = 1) + # size of label border) +\n",
    "        hexpand(1.8) +\n",
    "        geom_tippoint(size=1, alpha=0.8, aes(colour=phylogroup, stroke=1)) + \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",
    "        new_scale_fill() \n",
    "\n",
    "p3 <- p2+ theme(legend.position='none')\n",
    "\n",
    "df_text <- df[!duplicated(df$xindex),]\n",
    "p4 <- facet_plot(p3, panel = \"Lanthipeptide BiG-SCAPE GCFs\", data = df, geom=geom_tile, mapping = aes(x=xindex, y=y, fill = class_category), \n",
    "                 color = \"white\", lwd = 1,\n",
    "                 linetype = 1) \n",
    "#p3 <- facet_plot(p3, panel = \"BiG-SCAPE GCFs\", data = df_text, geom=geom_text, mapping = aes(x=xindex, y = 0, label = GCF), size=2, angle = 90)\n",
    "p4 <- facet_widths(p4, widths = c(1, 3)) \n",
    "#p3 + xlim_expand(c(150, 150), \"BiG-SCAPE GCFs\")\n",
    "ggsave(filename=figure_output_svg, plot=p4, device=svglite, width=2000, height=1400, units=\"px\")\n",
    "ggsave(plot=p4, width=2000, height=1400, units=\"px\", dpi=300, filename=figure_output_pdf, useDingbats=FALSE)\n",
    "p4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54e84199-5f7f-479d-ad87-73e98dbec317",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get coordinates of tiles\n",
    "with open(f'assets/figures/{FIGURE}/a.svg', \"r\") as f:\n",
    "    data = f.readlines()\n",
    "\n",
    "x_points = []\n",
    "for i in data:\n",
    "    if '<rect x=' in i:\n",
    "        if 'stroke-linecap: butt; stroke-linejoin:' in i:\n",
    "            for item in i.split():\n",
    "                if 'x=' in item:\n",
    "                    x_coor = float(item.strip(\"x=\").strip(\"'\"))\n",
    "                    x_points.append(x_coor)\n",
    "x_points = sorted(list(set(x_points)))\n",
    "assert len(x_points) == len(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd557cf8-3e38-43d2-a7d5-a0b42ab79846",
   "metadata": {},
   "outputs": [],
   "source": [
    "import svgutils.transform as sg\n",
    "\n",
    "# load matpotlib-generated figures\n",
    "fig = sg.fromfile(f'assets/figures/{FIGURE}/a.svg')\n",
    "fig2 = sg.fromfile(f'assets/figures/{FIGURE}/a.svg')\n",
    "\n",
    "# get the plot objects\n",
    "plot2 = fig2.getroot()\n",
    "#plot2.moveto(280, 0, scale=0.5)\n",
    "\n",
    "# add text labels\n",
    "for n, label in labels.items():\n",
    "    print(n, x_points[n], label)\n",
    "    x = x_points[n]\n",
    "    y = 330\n",
    "    if 'Erythreapeptin' in label:\n",
    "        weight = \"bold\"\n",
    "    else:\n",
    "        weight = \"normal\"\n",
    "    txt1 = sg.TextElement(x, y, label, size=7.5, weight=weight)\n",
    "    txt1.rotate(45, x, y)\n",
    "    fig.append(txt1)\n",
    "    # save generated SVG files\n",
    "    fig.save(f\"assets/figures/{FIGURE}/a_annotated.svg\")\n",
    "\n",
    "final_figure = Figure(\"650\", \"650\",\n",
    "                      Panel(\n",
    "                          SVG(f\"assets/figures/{FIGURE}/a_annotated.svg\").scale(1.2),\n",
    "                      ))\n",
    "final_figure.save(f\"assets/figures/{FIGURE}/a_annotated.svg\")\n",
    "final_figure"
   ]
  }
 ],
 "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
}