01_Fig3a_seqfu.python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aa0f0d95-5e6b-436f-a971-9486d0841329",
   "metadata": {},
   "source": [
    "# SeqFu2\n",
    "Summary of [Seqfu](https://github.com/telatin/seqfu2) results from project: `[{{ project().name }}]` \n",
    "\n",
    "## Description\n",
    "[Seqfu2](https://github.com/telatin/seqfu2) provides a sequence statistic overview of the genomes in the dataset.\n",
    "\n",
    "## Genome Statistics Overview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d21448bb-afe9-4fdb-b964-2143a0ccd1e4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pathlib import Path\n",
    "import altair as alt\n",
    "import yaml\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "974303df-fefc-4d21-82dc-420074563bc8",
   "metadata": {
    "tags": []
   },
   "source": [
    "%%bash\n",
    "#pip install -U altair_viewer\n",
    "#pip install altair_saver\n",
    "pip install vl-convert-python\n",
    "pip install selenium==4.2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a46880a2-ae00-465d-a4c0-5b2fc4c41391",
   "metadata": {},
   "source": [
    "## File Configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bab09367-1019-4c03-9401-9693de605473",
   "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": "bb9c9dce-20fe-4e6c-99fd-ce6731167069",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Set up paths and input file\n",
    "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_3\"\n",
    "\n",
    "# Seqfu result\n",
    "#seqfu_table = report_dir / \"tables/df_seqfu_stats.csv\"\n",
    "seqfu_table = bgcflow_dir / \"data/processed/qc_saccharopolyspora/tables/df_seqfu_stats.csv\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53fb8b15-379c-4be1-827e-243602f89d95",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Taxonomic placement result - merges GTDB API and GTDB-tk result\n",
    "gtdb_table = report_dir / \"tables/df_gtdb_meta.csv\"\n",
    "#gtdb_table = \"../../qc_saccharopolyspora/tables/df_gtdb_gtdbtk_meta.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bd02500-c7bb-4e77-a4f7-01bde6380c55",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# NCBI report for later use (automatically build when using samples from NCBI)\n",
    "df_ncbi = pd.read_csv(report_dir / \"tables/df_ncbi_meta.csv\").set_index(\"genome_id\")\n",
    "\n",
    "# Integrate both result into single table\n",
    "df_seqfu = pd.read_csv(seqfu_table)\n",
    "df_seqfu = df_seqfu.rename(columns={'File' : 'genome_id'}).set_index('genome_id')\n",
    "df_gtdb = pd.read_csv(gtdb_table).set_index('genome_id')\n",
    "df = pd.concat([df_seqfu, df_gtdb], join=\"inner\", axis=1).reset_index()\n",
    "df = df.set_index(\"genome_id\", drop=False)\n",
    "\n",
    "#df_all = pd.read_csv(\"../../qc_saccharopolyspora/tables/df_seqfu_stats.csv\").rename(columns={'File' : 'genome_id'}).set_index('genome_id', drop=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f827dc61-3771-4b21-b13d-2e69467de166",
   "metadata": {},
   "outputs": [],
   "source": [
    "# adds NCBI assembly level column in main table\n",
    "for i in df.index:\n",
    "    try:\n",
    "        df.loc[i, \"assembly_level\"] = df_ncbi.loc[i, \"assembly_level\"]\n",
    "    except KeyError:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3549045b-4dc8-4e34-98c7-1be646790e36",
   "metadata": {},
   "outputs": [],
   "source": [
    "count_cutoff = 50\n",
    "df_filterd = df[df.Count < count_cutoff]\n",
    "df_filterd_2 = df[df.N50 > 5000000]\n",
    "\n",
    "for genome in df.index:\n",
    "    if genome in df_filterd_2.index:\n",
    "        df.loc[genome, \"sequence_quality\"] = \"HQ\"\n",
    "    elif genome in df_filterd.index:\n",
    "        df.loc[genome, \"sequence_quality\"] = \"MQ\"\n",
    "    else:\n",
    "        df.loc[genome, \"sequence_quality\"] = \"LQ\"\n",
    "df.loc[:, [\"Total\", \"Count\", \"N50\", \"assembly_level\", \"sequence_quality\"]].sort_values(by=\"Count\")\n",
    "\n",
    "Path(\"assets/tables\").mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "df.to_csv(f\"assets/tables/{FIGURE}a_df_seqfu_annotated.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ff31eb5-8246-482f-a76e-40301a7e5278",
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile = Path(f\"assets/tables/{FIGURE}b_mash_hcluster.csv\")\n",
    "outfile.parent.mkdir(parents=True, exist_ok=True)\n",
    "df_mash_cluster = pd.read_csv(outfile, index_col=0)\n",
    "source = pd.concat([df, df_mash_cluster], axis=1)\n",
    "source.hcluster = source.hcluster.fillna(99)\n",
    "source[\"N50_ratio\"] = source[\"N50\"] / source[\"Total\"]\n",
    "source.phylogroup.fillna(\"NA\",inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "608195cb-5404-4b9a-9011-75252206a385",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "source[\"Genome length (Mbp)\"] = source[\"Total\"] / 1000000\n",
    "source = source.rename(columns={\"sequence_quality\" : \"Sequence Quality\"})\n",
    "x_col = 'gc'\n",
    "#y_col = 'Genome length (Mbp)'\n",
    "y_col = 'Total'\n",
    "\n",
    "\n",
    "# set up manual color\n",
    "domain = [\"P1\", \"P2\", \"P3\", \"P4\", \"P5\", \"P6\", \"P7\", \"P8\", \"NA\"]\n",
    "range_ = ['#264653','#e9c46a','#808080', '#808080','#f4a261','#808080', \"#e76f51\", \"#2a9d8f\", \"#FFFFFF\"]\n",
    "\n",
    "chart = alt.Chart(source).mark_point().encode(\n",
    "    alt.X(x_col, \n",
    "          scale=alt.Scale(domain=(round(source[x_col].min() - 0.01, 2), \n",
    "                                  round(source[x_col].max() + 0.01, 2))),\n",
    "          axis=alt.Axis(format='p', \n",
    "                        title=\"GC Content\")\n",
    "         ),\n",
    "    alt.Y(y_col,\n",
    "          scale=alt.Scale(domain=(round(source[y_col].min() - 0.2, 2), \n",
    "                                  round(source[y_col].max() + 0.2, 2))),\n",
    "          axis=alt.Axis(format='.2s',\n",
    "                        title=\"Genome Length (bp)\")\n",
    "         ),\n",
    "    shape=alt.Shape('Sequence Quality', legend=alt.Legend(orient='top')),\n",
    "    #size=\"N50_ratio\",\n",
    "    color=alt.Color(\"phylogroup:N\", scale=alt.Scale(domain=domain, range=range_), legend=None),\n",
    "    tooltip=['genome_id', 'Organism', 'Count', 'Total', 'gc', 'N50', 'AuN', 'Min', 'Max', \"phylogroup\"],\n",
    ").mark_point(\n",
    "    filled=True,\n",
    "    stroke='black',\n",
    "    strokeWidth=0.5,\n",
    "    opacity=0.8,\n",
    "    size=100\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",
    ")\n",
    "\n",
    "chart.interactive()\n",
    "chart.save(f\"assets/figures/{FIGURE}/a.svg\")\n",
    "chart.interactive().save(f\"assets/figures/{FIGURE}/a.html\")\n",
    "chart"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00188bd5-9cab-4268-8fa5-afbd321cfd34",
   "metadata": {},
   "source": [
    "[Download Table]({{ project().file_server() }}/tables/df_seqfu_stats.csv){:target=\"_blank\" .md-button}\n",
    "\n",
    "## References\n",
    "<font size=\"2\">\n",
    "{% for i in project().rule_used['seqfu']['references'] %}\n",
    "- *{{ i }}*\n",
    "{% endfor %}\n",
    "</font>"
   ]
  }
 ],
 "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
}