21_FigS6a_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 warnings\n",
    "import yaml\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8beb514-435e-45a3-8d0a-23a45c8898b7",
   "metadata": {},
   "source": [
    "## File Configurations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4426b2b8-2f93-4ad1-8eb4-9c83bf496fd2",
   "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": "60bb08cb-33bf-41eb-bf4a-3e7e182e9ecb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "bgcflow_dir = Path(notebook_configuration[\"bgcflow_dir\"])\n",
    "project_name = \"qc_saccharopolyspora\"\n",
    "report_dir = bgcflow_dir / f\"data/processed/{project_name}\"\n",
    "FIGURE = \"Figure_S6\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb9c9dce-20fe-4e6c-99fd-ce6731167069",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Seqfu result\n",
    "seqfu_table = report_dir / \"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",
    "gtdbtk_table = report_dir / \"tables/gtdbtk.bac120.summary.tsv\"\n",
    "df_gtdb = pd.read_csv(gtdb_table).set_index('genome_id')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10cc0cdd-e07c-4a9a-aff9-ac19f76b1f4d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_gtdbtk = pd.read_csv(gtdbtk_table, sep=\"\\t\").rename(columns={\"user_genome\" : \"genome_id\"}).set_index(\"genome_id\")\n",
    "tax_mapping = {\"d\" : \"Domain\",\n",
    "              \"p\" : \"Phylum\",\n",
    "              \"c\" : \"Class\",\n",
    "              \"o\" : \"Order\",\n",
    "              \"f\" : \"Family\",\n",
    "              \"g\" : \"Genus\",\n",
    "              \"s\" : \"Organism\"}\n",
    "\n",
    "for index in df_gtdbtk.index:\n",
    "    tax = [i for i in df_gtdbtk.loc[index, \"classification\"].split(\";\")]\n",
    "    for level in tax:\n",
    "        key = level.split(\"__\")[0]\n",
    "        if key == \"g\":\n",
    "            genus = level.split(\"__\")[-1]\n",
    "        if level == \"s__\":\n",
    "            level = f\"s__{genus} sp.\"\n",
    "        df_gtdb.loc[index, tax_mapping[key]] = level\n",
    "df_gtdb.Species = [i.split()[-1] for i in df_gtdb.Organism]\n",
    "df_gtdb.to_csv(f\"assets/tables/{FIGURE}_df_gtdb.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",
    "\n",
    "df = pd.concat([df_seqfu, df_gtdb], join=\"inner\", axis=1).reset_index()\n",
    "df = df.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}_df_seqfu_annotated.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d8fe1ca-30f9-472a-9e01-4e9db61825b8",
   "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
}