05_Fig4b_query_bigfam.python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f32bf8b3-dd6a-47c0-8ec9-518679e96b2c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import sqlite3\n",
    "\n",
    "import math\n",
    "import json, yaml\n",
    "import altair as alt\n",
    "import ast\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import networkx as nx\n",
    "import plotly.graph_objects as go\n",
    "import altair as alt\n",
    "\n",
    "def shannon_div(data):\n",
    "    total = sum(data.values())\n",
    "    container = []\n",
    "    for k, v in data.items():\n",
    "        p = v / total\n",
    "        val = p*(math.log(p))\n",
    "        container.append(val)\n",
    "    shannon = abs(sum(container))\n",
    "    return shannon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d801d33d-f29e-46e6-9d73-7012a271a128",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_bigfam_class(gcf_id, conn, threshold=\"<= 900\", rank=0):\n",
    "    \"\"\"\n",
    "    Given a gcf_id, returns chemical class and chemical subclass of the model\n",
    "    \"\"\"\n",
    "    sql_query = f\"\"\"\n",
    "    SELECT \n",
    "        \"gcf_membership\".\"gcf_id\" AS \"gcf_id\", \n",
    "        \"gcf_membership\".\"bgc_id\" AS \"bgc_id\",      \n",
    "        \"gcf_membership\".\"membership_value\" AS \"membership_value\", \n",
    "        \"gcf_membership\".\"rank\" AS \"rank\", \n",
    "        \"Bgc Class - Bgc\".\"bgc_id\" AS \"Bgc Class - Bgc__bgc_id\", \n",
    "        \"Bgc Class - Bgc\".\"chem_subclass_id\" AS \"Bgc Class - Bgc__chem_subclass_id\", \n",
    "        \"Chem Subclass\".\"id\" AS \"Chem Subclass__id\", \n",
    "        \"Chem Subclass\".\"class_id\" AS \"Chem Subclass__class_id\", \n",
    "        \"Chem Subclass\".\"name\" AS \"Chem Subclass__name\", \n",
    "        \"Chem Class - Class\".\"id\" AS \"Chem Class - Class__id\", \n",
    "        \"Chem Class - Class\".\"name\" AS \"Chem Class - Class__name\"\n",
    "    FROM \"gcf_membership\"\n",
    "    LEFT JOIN \"bgc_class\" \"Bgc Class - Bgc\" ON \"gcf_membership\".\"bgc_id\" = \"Bgc Class - Bgc\".\"bgc_id\" LEFT JOIN \"chem_subclass\" \"Chem Subclass\" ON \"Bgc Class - Bgc\".\"chem_subclass_id\" = \"Chem Subclass\".\"id\" LEFT JOIN \"chem_class\" \"Chem Class - Class\" ON \"Chem Subclass\".\"class_id\" = \"Chem Class - Class\".\"id\"\n",
    "    WHERE (\"gcf_membership\".\"gcf_id\" = {gcf_id}\n",
    "       AND \"gcf_membership\".\"membership_value\" {threshold} AND \"gcf_membership\".\"rank\" = {rank})\n",
    "    \"\"\"\n",
    "    df = pd.read_sql_query(sql_query, conn)\n",
    "    non_singleton_map = df.bgc_id.value_counts().to_dict()\n",
    "    bgc = [k for k, v in non_singleton_map.items() if v > 1]\n",
    "    df.loc[df[~df.bgc_id.isin(bgc)].index, \"score\"] = 1\n",
    "    for i in df[df.bgc_id.isin(bgc)].index:\n",
    "        score = 1 / non_singleton_map[df.loc[i, \"bgc_id\"]]\n",
    "        df.loc[i, \"score\"] = score\n",
    "        \n",
    "    chemical_class = df.groupby(\"Chem Class - Class__name\").sum().score / len(df.bgc_id.unique())\n",
    "    chemical_class = chemical_class.to_dict()\n",
    "    chemical_subclass = df.groupby(\"Chem Subclass__name\").sum().score / len(df.bgc_id.unique())\n",
    "    chemical_subclass = chemical_subclass.to_dict()\n",
    "    top_chemical_class = max(chemical_class, key=chemical_class.get)\n",
    "    top_chemical_subclass = max(chemical_subclass, key=chemical_subclass.get)\n",
    "    bgc_member = len(df.bgc_id.unique())\n",
    "    \n",
    "    result = {\"bgc_member\" : bgc_member,\n",
    "              \"chemical_class_hits\" : len(df),\n",
    "              \"top_chemical_class\" : top_chemical_class,\n",
    "              \"top_chemical_class_proportion\" : chemical_class[top_chemical_class],\n",
    "              \"top_chemical_subclass\" : top_chemical_subclass,\n",
    "              \"top_chemical_subclass_proportion\" : chemical_subclass[top_chemical_subclass],\n",
    "              \"chemical_class\" : chemical_class,\n",
    "              \"chemical_subclass\" : chemical_subclass,\n",
    "              }\n",
    "    return result\n",
    "\n",
    "def get_bigfam_taxa(gcf_id, conn, threshold=\"<= 900\", rank=0, level=5):\n",
    "    \"\"\"\n",
    "    Given a gcf_id, returns chemical class and chemical subclass of the model\n",
    "    \"\"\"\n",
    "    sql_all_hits = f\"\"\"\n",
    "    SELECT \n",
    "        \"gcf_membership\".\"gcf_id\" AS \"gcf_id\", \n",
    "        \"gcf_membership\".\"bgc_id\" AS \"bgc_id\", \n",
    "        \"gcf_membership\".\"membership_value\" AS \"membership_value\", \n",
    "        \"gcf_membership\".\"rank\" AS \"rank\"\n",
    "    FROM \"gcf_membership\"\n",
    "    WHERE (\"gcf_membership\".\"gcf_id\" = {gcf_id}\n",
    "        AND \"gcf_membership\".\"membership_value\" {threshold} AND \"gcf_membership\".\"rank\" = {rank})\n",
    "    \"\"\"\n",
    "    \n",
    "    sql_query = f\"\"\"\n",
    "    SELECT \n",
    "        \"gcf_membership\".\"gcf_id\" AS \"gcf_id\", \n",
    "        \"gcf_membership\".\"bgc_id\" AS \"bgc_id\", \n",
    "        \"gcf_membership\".\"membership_value\" AS \"membership_value\", \n",
    "        \"gcf_membership\".\"rank\" AS \"rank\", \n",
    "        \"Bgc Taxonomy - Bgc\".\"bgc_id\" AS \"Bgc Taxonomy - Bgc__bgc_id\", \n",
    "        \"Bgc Taxonomy - Bgc\".\"taxon_id\" AS \"Bgc Taxonomy - Bgc__taxon_id\", \n",
    "        \"Taxon\".\"id\" AS \"Taxon__id\", \"Taxon\".\"level\" AS \"Taxon__level\", \n",
    "        \"Taxon\".\"name\" AS \"Taxon__name\", \"Taxon Class - Level\".\"id\" AS \"Taxon Class - Level__id\", \n",
    "        \"Taxon Class - Level\".\"level\" AS \"Taxon Class - Level__level\", \n",
    "        \"Taxon Class - Level\".\"name\" AS \"Taxon Class - Level__name\"\n",
    "    FROM \"gcf_membership\"\n",
    "    LEFT JOIN \"bgc_taxonomy\" \"Bgc Taxonomy - Bgc\" ON \"gcf_membership\".\"bgc_id\" = \"Bgc Taxonomy - Bgc\".\"bgc_id\" LEFT JOIN \"taxon\" \"Taxon\" ON \"Bgc Taxonomy - Bgc\".\"taxon_id\" = \"Taxon\".\"id\" LEFT JOIN \"taxon_class\" \"Taxon Class - Level\" ON \"Taxon\".\"level\" = \"Taxon Class - Level\".\"level\"\n",
    "    WHERE (\"gcf_membership\".\"gcf_id\" = {gcf_id}\n",
    "        AND \"gcf_membership\".\"membership_value\" {threshold}\n",
    "        AND \"gcf_membership\".\"rank\" = {rank}\n",
    "        AND \"Taxon\".\"level\" = {level})\n",
    "    \"\"\"\n",
    "    df_all_hits = pd.read_sql_query(sql_all_hits, conn)\n",
    "    df = pd.read_sql_query(sql_query, conn)\n",
    "    taxonomic_level = df['Taxon Class - Level__name'].unique()[0]\n",
    "    df = df_all_hits.merge(df, on=[\"bgc_id\", \"gcf_id\", \"membership_value\", \"rank\"], how=\"outer\").fillna(\"Unassigned\")\n",
    "    taxa_distribution = df.Taxon__name.value_counts().to_dict()\n",
    "    top_genus = max(taxa_distribution, key=taxa_distribution.get)\n",
    "    taxonomic_hits = len(df.bgc_id.unique())\n",
    "    result = {\"taxonomic_hits\" : taxonomic_hits,\n",
    "              \"taxonomic_level\" : taxonomic_level,\n",
    "              \"H-index\" : shannon_div(taxa_distribution),\n",
    "              \"richness\" : len(taxa_distribution.keys()),\n",
    "              \"top_taxa\" : top_genus,\n",
    "              \"top_taxa_proportion\" : taxa_distribution[top_genus] / taxonomic_hits,\n",
    "              \"taxa_distribution\" : taxa_distribution,\n",
    "             }\n",
    "    return result\n",
    "\n",
    "def get_bigfam_summary(gcf_id, conn, threshold=\"<= 900\", rank=0, level=5):\n",
    "    result = get_bigfam_class(gcf_id, conn)\n",
    "    result.update(get_bigfam_taxa(gcf_id, conn))\n",
    "    return {gcf_id : result}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4132d92-557e-4053-8382-a59d8739f66f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def generate_bigfam_network(df, rank=0):\n",
    "    \"\"\"\n",
    "    Build a networkx graph from bigscape df network\n",
    "    \"\"\"\n",
    "    df[df['rank'] == rank]\n",
    "    G = nx.from_pandas_edgelist(df, source='bgc_id', target='gcf_id', edge_attr=['membership_value',\n",
    "           'rank'], edge_key='bigfam_edge_id')\n",
    "    return G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c95c3d6b-d48c-4874-9c39-02a5aa2f27cb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def annotate_bigfam_models(G, df, columns = ['bgc_member', 'chemical_class_hits', 'top_chemical_class',\n",
    "                                             'top_chemical_class_proportion', 'top_chemical_subclass',\n",
    "                                             'top_chemical_subclass_proportion', 'taxonomic_hits', \n",
    "                                             'taxonomic_level', 'H-index', 'richness', 'top_taxa', \n",
    "                                             'top_taxa_proportion']):\n",
    "    for bgc in G.nodes:\n",
    "        if bgc in df.index:\n",
    "            G.nodes[bgc]['node_type'] = \"BiG-FAM GCFs\"\n",
    "            G.nodes[bgc]['text'] = f\"GCF {bgc}<br>size: {df.loc[bgc, 'bgc_member']}<br>top_chemical_class: {df.loc[bgc, 'top_chemical_class']} ({df.loc[bgc, 'top_chemical_class_proportion']:.0%})\\\n",
    "<br>top_chemical_subclass: {df.loc[bgc, 'top_chemical_subclass']} ({df.loc[bgc, 'top_chemical_subclass_proportion']:.0%})\\\n",
    "<br>top_taxa: {df.loc[bgc, 'top_taxa']} ({df.loc[bgc, 'top_taxa_proportion']:.0%})\"\n",
    "            for c in columns:\n",
    "                G.nodes[bgc][c] = df.loc[bgc, c]\n",
    "    return G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf9fc3e7-8e46-4ac6-bd1b-f9c4383f02e2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def annotate_bigfam_antismash(G, antismash_region_path, columns = [\"genome_id\", \"product\", \"contig_edge\", \"region_length\", \"most_similar_known_cluster_id\", \"most_similar_known_cluster_description\", \"most_similar_known_cluster_type\", \"similarity\"]):\n",
    "    df_antismash = pd.read_csv(antismash_region_path, index_col=0)\n",
    "    df_antismash = df_antismash.loc[:, columns]\n",
    "    for bgc in G.nodes:\n",
    "        if bgc in df_antismash.index:\n",
    "            G.nodes[bgc]['node_type'] = \"BGC\"\n",
    "            G.nodes[bgc]['text'] = f\"{bgc}<br>{df_antismash.loc[bgc, 'product']}<br>{df_antismash.loc[bgc, 'most_similar_known_cluster_description']}\"\n",
    "            for c in df_antismash.columns:\n",
    "                G.nodes[bgc][c] = df_antismash.loc[bgc, c]\n",
    "    return G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fefefba-8c6a-491f-91f8-c11b85dba89c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def create_edge_trace(G, name):\n",
    "    edge_trace = go.Scatter(\n",
    "        x=[],\n",
    "        y=[],\n",
    "        name=name,\n",
    "        line=dict(width=0.5,color='#888'),\n",
    "        hoverinfo='none',\n",
    "        mode='lines')\n",
    "\n",
    "    for edge in G.edges():\n",
    "        x0, y0 = G.nodes[edge[0]]['pos']\n",
    "        x1, y1 = G.nodes[edge[1]]['pos']\n",
    "        edge_trace['x'] += tuple([x0, x1, None])\n",
    "        edge_trace['y'] += tuple([y0, y1, None])\n",
    "    return edge_trace\n",
    "\n",
    "def create_node_trace(G, node_trace_category, color, showtextlabel=False, nodesize=10, nodeopacity=0.8, nodesymbol=\"circle\", linewidth=1, linecolor=\"black\", textposition=\"top center\"):\n",
    "    if showtextlabel:\n",
    "        markermode = \"markers+text\"\n",
    "    else:\n",
    "        markermode = \"markers\"\n",
    "    node_trace = go.Scatter(\n",
    "            x=[],\n",
    "            y=[],\n",
    "            text=[],\n",
    "            textposition=textposition,\n",
    "            mode=markermode,\n",
    "            hoverinfo='text',\n",
    "            name=node_trace_category,\n",
    "            marker=dict(\n",
    "                symbol=nodesymbol,\n",
    "                opacity=nodeopacity,\n",
    "                showscale=False,\n",
    "                color=color,\n",
    "                size=nodesize,\n",
    "                line=dict(width=linewidth, color=linecolor)))\n",
    "\n",
    "    for node in G.nodes():\n",
    "        if G.nodes[node][\"node_trace\"] == node_trace_category:\n",
    "            x, y = G.nodes[node]['pos']\n",
    "            node_trace['x'] += tuple([x])\n",
    "            node_trace['y'] += tuple([y])\n",
    "            node_trace['text'] +=tuple([G.nodes[node]['text']])\n",
    "    return node_trace"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c2b4173-f7e2-4436-9b4f-c27a9953eaf7",
   "metadata": {},
   "source": [
    "## File configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e31764e5-0fd3-4190-82b1-474df301ccfb",
   "metadata": {
    "tags": []
   },
   "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": "be30544d-7c98-4e18-9eda-2bbc6fc6f959",
   "metadata": {
    "tags": []
   },
   "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_4\"\n",
    "MASH_FIGURE = \"Figure_3\"\n",
    "\n",
    "query_dir = report_dir / \"bigslice/query_as_6.1.1/\"\n",
    "df_bigfam_model = pd.read_csv(query_dir / \"gcf_summary.csv\")\n",
    "df_bigfam_hits = pd.read_csv(report_dir / \"bigslice/query_as_6.1.1/query_network.csv\")\n",
    "antismash_region_path = report_dir / \"tables/df_regions_antismash_6.1.1.csv\"\n",
    "mash_path = f\"assets/tables/{MASH_FIGURE}b_mash_hcluster.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9265868-76ee-4301-b85a-a3d0e3e39964",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "output = Path(f\"assets/tables/{FIGURE}_bigfam_models.csv\")\n",
    "\n",
    "if output.is_file():\n",
    "    df = pd.read_csv(output, index_col=0)\n",
    "    pass\n",
    "else:\n",
    "    bigfam_db = bgcflow_dir / \"resources/bigslice/full_run_result/result/data.db\"\n",
    "    conn = sqlite3.connect(bigfam_db)\n",
    "\n",
    "    result = {}\n",
    "    for i in df_bigfam_model.gcf_id:\n",
    "        value = get_bigfam_summary(i, conn)\n",
    "        result.update(value)\n",
    "\n",
    "    df = pd.DataFrame.from_dict(result).T\n",
    "    df.to_csv(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67ea9fd7-65a1-41e9-8a1a-39660113a0df",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bigfam_network_path = report_dir / \"bigslice/query_as_6.1.1/query_network.csv\"\n",
    "\n",
    "df_network = pd.read_csv(Path(bigfam_network_path))\n",
    "df_cut = df[df.top_taxa_proportion <= 0.3]\n",
    "#df_network = pd.read_csv(Path(bigfam_network_path))\n",
    "#df_network = df_network[df_network.gcf_id.isin(df_cut.index)]\n",
    "\n",
    "G = generate_bigfam_network(df_network)\n",
    "G = annotate_bigfam_antismash(G, antismash_region_path)\n",
    "G = annotate_bigfam_models(G, df)\n",
    "G_raw = G.copy()\n",
    "\n",
    "for n in df_cut.index:\n",
    "    print(f\"removing {n}\")\n",
    "    G.remove_node(n)\n",
    "\n",
    "# position nodes\n",
    "pos = nx.nx_agraph.graphviz_layout(G)\n",
    "for n, p in pos.items():\n",
    "    G.nodes[n]['pos'] = p\n",
    "    \n",
    "# Add phylogorup information\n",
    "df_mash = pd.read_csv(mash_path, index_col=0)\n",
    "for node in G.nodes:\n",
    "    if 'genome_id' in G.nodes[node].keys():\n",
    "        G.nodes[node]['phylogroup'] = df_mash.loc[G.nodes[node]['genome_id'], \"phylogroup\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a88b66a2-8154-4aa7-b9a2-549c864e8c79",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Calculate the number of GCFs with representation across different genera\n",
    "\n",
    "\n",
    "df_genus_dist = pd.DataFrame(index=df.index)\n",
    "for gcf_id in df.index:\n",
    "    taxa = str(df.loc[gcf_id, \"taxa_distribution\"])\n",
    "    taxa_dict = ast.literal_eval(taxa)\n",
    "    for genus in taxa_dict.keys():\n",
    "        df_genus_dist.loc[gcf_id, genus] = taxa_dict[genus]\n",
    "df_genus_dist.fillna(0, inplace=True)\n",
    "\n",
    "# Genera represented by most number of GCFs\n",
    "df_genus_dist_binary = df_genus_dist > 0 \n",
    "df_genus_dist_binary = df_genus_dist_binary * 1\n",
    "df_genus_dist_binary.sum().sort_values(ascending=False)[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cef6d45-0ba7-46ef-88bf-2652b1b73e74",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#node_color = phylogroup\n",
    "phylomap = {}\n",
    "color = [\"#264653\", \"#287271\", \"#2a9d8f\", \"#8ab17d\", \"#e9c46a\", \"#f4a261\", \"#ee8959\", \"#e76f51\"]\n",
    "for num, phylogroup in enumerate(df_mash.phylogroup.unique()):\n",
    "    phylomap[phylogroup] = color[num]\n",
    "phylomap[\"BiG-FAM GCFs\"] = \"blue\""
   ]
  },
  {
   "cell_type": "raw",
   "id": "8c29e19b-7a28-41e0-a058-e779470d7171",
   "metadata": {
    "tags": []
   },
   "source": [
    "node_trace_category = {\"BGC\" : {\"nodecolor\" : \"blue\",\n",
    "                                \"nodeshape\" : \"circle\",\n",
    "                                \"linecolor\" : \"black\"},\n",
    "                       \"BiG-FAM GCFs\" : {\"nodecolor\" : \"red\",\n",
    "                                         \"nodeshape\" : \"square\",\n",
    "                                         \"linecolor\" : \"black\"},\n",
    "                      }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec8b061f-f4f7-494f-a00a-d71f98d1cf55",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "color_category = 'phylogroup'\n",
    "usecolormap = phylomap\n",
    "node_trace_category = {}\n",
    "for node in G.nodes:\n",
    "    if G.nodes[node]['node_type'] == \"BiG-FAM GCFs\":\n",
    "        nodeshape = \"diamond\"\n",
    "        color = usecolormap[\"BiG-FAM GCFs\"]\n",
    "        cat = \"BiG-FAM GCFs\"\n",
    "    else:            \n",
    "        cat = G.nodes[node][color_category]\n",
    "        color = usecolormap[cat]\n",
    "        if G.nodes[node]['similarity'] > 0.8:\n",
    "            cat = cat + \" \" + \"Known (antiSMASH) > 80% similarity\"\n",
    "            nodeshape = \"circle\"\n",
    "        elif G.nodes[node]['similarity'] > 0.4:\n",
    "            cat = cat + \" \" + \"Known (antiSMASH) > 40% similarity\"\n",
    "            nodeshape = \"triangle-up\"\n",
    "        elif G.nodes[node]['similarity'] < 0.4:\n",
    "            cat = cat + \" \" + \"Unknown (antiSMASH) < 40% similarity\"\n",
    "            nodeshape = \"triangle-down\"\n",
    "        else:\n",
    "            cat = cat + \" \" + \"Unknown\"\n",
    "            nodeshape = \"star\"\n",
    "        \n",
    "        linecolor = \"black\"\n",
    "    \n",
    "    G.nodes[node]['node_trace'] = cat\n",
    "    node_trace_category[cat] = {\"nodecolor\" : color,\n",
    "                                \"nodeshape\" : nodeshape,\n",
    "                                \"linecolor\" : linecolor}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa3e9479-2b91-480e-b639-e53e106fffac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "edge_trace = create_edge_trace(G, \"bigslice_similarity\")\n",
    "traces = [edge_trace]\n",
    "for cat in node_trace_category.keys():\n",
    "    nodeopacity = 0.8\n",
    "    color = node_trace_category[cat][\"nodecolor\"]\n",
    "    nodeshape = node_trace_category[cat][\"nodeshape\"]\n",
    "    linecolor = node_trace_category[cat][\"linecolor\"]\n",
    "    node_trace = create_node_trace(G, cat, color, nodesymbol=nodeshape, nodeopacity=nodeopacity, nodesize=10, linewidth=1, linecolor=linecolor)\n",
    "    traces.append(node_trace)\n",
    "    \n",
    "fig2 = go.Figure(data=traces,\n",
    "                layout=go.Layout(\n",
    "                    paper_bgcolor='rgba(0,0,0,0)',\n",
    "                    plot_bgcolor='rgba(0,0,0,0)',\n",
    "                    showlegend=True,\n",
    "                    hovermode='closest',\n",
    "                    margin=dict(b=20,l=5,r=5,t=40),\n",
    "                    xaxis=dict(showgrid=False, zeroline=False, showticklabels=False),\n",
    "                    yaxis=dict(showgrid=False, zeroline=False, showticklabels=False),\n",
    "                    width=1200, height=800))\n",
    "fig2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eecd074a-876d-43b5-b251-ff08269b5d9b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bgc_hits = len([n for n in G_raw.nodes if G_raw.nodes[n]['node_type'] == 'BGC'])\n",
    "GCF_hits = len([n for n in G_raw.nodes if G_raw.nodes[n]['node_type'] != 'BGC'])\n",
    "top_taxa_saccharopolyspora = df[df.top_taxa_proportion == 1].top_taxa.value_counts().to_dict()['Saccharopolyspora']\n",
    "shannon_non_zero = df[df['H-index'] != 0]\n",
    "shannon_non_zero[shannon_non_zero['H-index'] > 2].shape\n",
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d4ebbb1-d5bb-466e-8b3d-5f81c6632c15",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#GCF_201888\n",
    "\n",
    "shannon_min_non_zero_tax = ast.literal_eval(str(df.loc[shannon_non_zero['H-index'].astype(float).idxmin(), \"taxa_distribution\"]))\n",
    "shannon_min_non_zero_tax = {k : v/sum(shannon_min_non_zero_tax.values()) for k,v in shannon_min_non_zero_tax.items()}\n",
    "shannon_min_non_zero_tax_clean = [f\"{k} ({v:.1%})\" for k,v in sorted(shannon_min_non_zero_tax.items(), key=lambda x:x[1], reverse=True)]\n",
    "df.loc[220008, 'top_taxa']\n",
    "df.loc[220008, 'top_taxa_proportion'] * df.loc[220008, 'bgc_member']\n",
    "def get_top_taxa(df, gcf):\n",
    "    result = f\"{df.loc[gcf, 'top_taxa']} ({df.loc[gcf, 'top_taxa_proportion'] * df.loc[gcf, 'bgc_member']:.0f} genomes)\"\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43ff49d9-f981-46d7-84ea-35874c21cad0",
   "metadata": {},
   "outputs": [],
   "source": [
    "shannon_non_zero['H-index'].astype(float).idxmin()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cdf15bd-d939-417e-884c-113a2aefd303",
   "metadata": {},
   "source": [
    "The comparison of BGCs against BiG-FAM-DB can be used to investigate if detected BGCs are widely spread across other genomes from the public datasets. The BiG-SLICE-based query resulted in 389 BGCs having similarities against 130 GCF models from the BiGFAM database. From the 130 detected GCFs, there were 71 BiG-FAM GCFs that are specifically distributed in the genus Saccharopolyspora. Whereas several other GCFs included BGCs from different genera such as Streptomyces (47 genomes), Amycolatopsis (31 genomes), Kitasatospora (29 genomes), Nocardia (25 genomes), Pseudomonas_E (22 genomes), and many others. We further calculated the Shannon diversity index (H) for each of the GCFs representing the distribution of BGCs across different genera (Figure 4b). Of the 59 GCFs with positive Shannon index, 33 were highly distributed across many genera with a Shannon index of greater than 2. The GCF GCF_201888 with the lowest positive Shannon index (~0.3) contained 12,444 BGCs distributed across 43 genera with the majority belonging to Staphylococcus (~94.2%) followed by Acinetobacter (~4.2%). The known BGC in this BiG-FAM GCF coded for the biosynthesis of staphylobactin (also known as staphyloferrin B), which is a siderophore with a role in the virulence of Staphylococcus aureus (84, 85). The detailed comparative analysis of the predicted BGCs against the MIBIG entry showed that the Saccharopolypora genomes indeed possess a BGC that is very similar to the staphylobactin BGC from Staphylococcus aureus (Figure S11). Some of the other GCFs like GCF_215691  (Shannon index: 0.35), GCF_220008 (Shannon index: 0.41), and GCF_203048 (Shannon index: 0.52) were found predominantly in genera such as Psudeomonas_E (1241 genomes), Mycobacterium (227 genomes), and Streptomyces (114 genomes), respectively (Figure 4B). These examples particularly highlight the BGCs that are potentially transferred across different phylogenetic groups through horizontal gene transfer events."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a9a3298-4e88-49c9-a35a-c04b7244672c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "text1 = f\"The comparison of BGCs against BiG-FAM-DB can be used to investigate if detected BGCs are widely spread across other genomes from the public datasets. The BiG-SLICE-based query resulted in {bgc_hits} BGCs having similarities against {GCF_hits} GCF models from the BiGFAM database.\"\n",
    "text2 = f\"From the {GCF_hits} detected GCFs, there were {top_taxa_saccharopolyspora} BiG-FAM GCFs that are specifically distributed in the genus Saccharopolyspora. Whereas several other GCFs included BGCs from different genera such as Streptomyces (47 genomes), Amycolatopsis (31 genomes), Kitasatospora (29 genomes), Nocardia (25 genomes), Pseudomonas_E (22 genomes), and many others.\"\n",
    "text3 = f\"We further calculated the Shannon diversity index (H) for each of the GCFs representing the distribution of BGCs across different genera (Figure 4b). Of the {shannon_non_zero.shape[0]} GCFs with positive Shannon index, {shannon_non_zero[shannon_non_zero['H-index'] > 2].shape[0]} were highly distributed across many genera with a Shannon index of greater than 2. The GCF GCF_{shannon_non_zero['H-index'].idxmin()} with the lowest positive Shannon index (~{shannon_non_zero['H-index'].min():.1f}) contained {df.loc[shannon_non_zero['H-index'].idxmin(), 'bgc_member']} BGCs distributed across {len(shannon_min_non_zero_tax.keys())} genera with the majority belonging to {shannon_min_non_zero_tax_clean[0]}) followed by {shannon_min_non_zero_tax_clean[1]}.\"\n",
    "text4 = f\"The known BGC in this BiG-FAM GCF coded for the biosynthesis of staphylobactin (also known as staphyloferrin B), which is a siderophore with a role in the virulence of Staphylococcus aureus (84, 85). The detailed comparative analysis of the predicted BGCs against the MIBIG entry showed that the Saccharopolypora genomes indeed possess a BGC that is very similar to the staphylobactin BGC from Staphylococcus aureus (Figure S11).\"\n",
    "text5 = f\"Some of the other GCFs like GCF_215691 (Shannon index: {df.loc[215691, 'H-index']:.2f}), GCF_220008 (Shannon index: {df.loc[220008, 'H-index']:.2f}), and GCF_203048 (Shannon index: {df.loc[203048, 'H-index']:.2f}) were found predominantly in genera such as {get_top_taxa(df, 215691)}, {get_top_taxa(df, 220008)}, and {get_top_taxa(df, 203048)}, respectively (Figure 4B). These examples particularly highlight the BGCs that are potentially transferred across different phylogenetic groups through horizontal gene transfer events.\"\n",
    "final_text = \" \".join([text1, text2, text3, text4, text5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b8a1ecf-8d04-4a4f-bbeb-f23ab095c28e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with open(f\"assets/figures/{FIGURE}/{FIGURE}b_text.txt\", \"w\") as f:\n",
    "    f.write(final_text)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "441848a3-b079-45fc-a31b-b418d77c2a04",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mapping = df_bigfam_hits.gcf_id.value_counts()\n",
    "for gcf in df.index:\n",
    "    df.loc[gcf, \"dataset_hits\"] = mapping[gcf]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb7827b3-dfc2-4bf6-bf6a-ec6847fe37dc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "domain = list(df[df[\"top_taxa_proportion\"] > 0.5][\"top_taxa\"].value_counts().to_dict().keys())\n",
    "domain.append(\"Other\")\n",
    "print(domain)\n",
    "r = []\n",
    "range_ = [\"#264653\", \"#287271\", \"#2a9d8f\", \"#8ab17d\", \"#e9c46a\", \"#f4a261\", \"#ee8959\", \"#e76f51\", \"#edede9\"]\n",
    "for num, d in enumerate(domain):\n",
    "    if num < len(range_):\n",
    "        r.append(range_[num])\n",
    "    else:\n",
    "        r.append(\"white\")\n",
    "#range_ = [\"#264653\", \"white\", \"white\", \"white\", \"white\", \"white\", \"white\", \"white\", \"white\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cc8dc6a-b181-42c6-ab41-7053d7ab2d26",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "source = df.copy()\n",
    "source = source.reset_index().rename(columns={\"index\":\"BiG-FAM_id\"})\n",
    "for i in source.index:\n",
    "    if source.loc[i, \"top_taxa_proportion\"] <= 0.5:\n",
    "        source.loc[i, \"top_taxa\"] = \"Other\"\n",
    "        \n",
    "chart_one = alt.Chart(source).mark_point().encode(\n",
    "    alt.Y('H-index:Q',\n",
    "          scale=alt.Scale(domain=(-0.5, 5)),\n",
    "          axis=alt.Axis(title=\"Shannon Index (H')\")\n",
    "         ),\n",
    "    alt.X('bgc_member:Q',\n",
    "          scale=alt.Scale(type=\"log\"),\n",
    "          axis=alt.Axis(title=\"Member Size\")\n",
    "         ),\n",
    "    alt.Size('dataset_hits',\n",
    "             scale=alt.Scale(type='pow', domain=(1, 30)), \n",
    "             title=\"Number of hits in dataset\"\n",
    "            ),\n",
    "    alt.Color(\"top_taxa:N\", scale=alt.Scale(domain=domain, range=r), title=\"Top Genus (>=50%)\"),\n",
    "    tooltip=['BiG-FAM_id', 'bgc_member', 'chemical_class_hits', 'top_chemical_class', alt.Tooltip('top_chemical_class_proportion', format='.1%'), \n",
    "             'top_chemical_subclass', alt.Tooltip('top_chemical_subclass_proportion', format='.1%'),\n",
    "             'taxonomic_level', 'richness', 'top_taxa', alt.Tooltip('top_taxa_proportion', format='.1%'), \n",
    "             alt.Tooltip('H-index:Q', format='.2f')],\n",
    ").mark_point(\n",
    "    filled=True,\n",
    "    stroke='black',\n",
    "    strokeWidth=0.5,\n",
    "    opacity=0.8,\n",
    "    size=1000\n",
    ").configure_header(\n",
    "    title=None,\n",
    "    labels=False\n",
    ").configure_axis(\n",
    "    labelFontSize=10,\n",
    "    titleFontSize=12\n",
    ").configure_legend(\n",
    "    labelFontSize=10,\n",
    "    titleFontSize=12,\n",
    ").configure_view(\n",
    "    continuousHeight=250,\n",
    "    continuousWidth=250,\n",
    ").configure_legend(\n",
    "  orient='right'\n",
    ")\n",
    "outfile = Path(f\"assets/figures/{FIGURE}/{FIGURE}b.svg\")\n",
    "outfile_html = Path(f\"assets/figures/{FIGURE}/{FIGURE}b.html\")\n",
    "outfile.parent.mkdir(parents=True, exist_ok=True)\n",
    "chart_one.save(outfile)\n",
    "chart_one.interactive().save(outfile_html)\n",
    "chart_one.interactive()"
   ]
  }
 ],
 "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}