15_Fig6cd_Lanthipeptide_erythreapeptin_gggenomes.R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b0f08259-b159-4c4a-8948-2bce9fcad16c",
   "metadata": {},
   "source": [
    "## Multiple Sequence Alignment of Lanthipeptide candidate clusters\n",
    "This analysis requires project `` containing  BGCs that . The project was run through the `BGC` comparison subworkflow by running `bgcflow run --snakefile workflow/BGC`"
   ]
  },
  {
   "cell_type": "raw",
   "id": "00ba0bcc-a0b4-48b4-81c4-f5a4e614f32c",
   "metadata": {},
   "source": [
    "%%bash\n",
    "pip install pyarrow\n",
    "pip install networkx\n",
    "pip install biopython\n",
    "mamba install -c conda-forge r-ggalt -y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dba257e-d2a2-4a58-b7b7-c80c9340f77a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from pathlib import Path\n",
    "from Bio import SeqIO\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "from Bio.Seq import Seq\n",
    "import yaml\n",
    "\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12591861-c6b7-4b20-801c-8e193a033de8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%R #Run this to install R dependencies\n",
    "#if (!require(\"BiocManager\", quietly = TRUE))\n",
    "#    install.packages(\"BiocManager\")\n",
    "    \n",
    "if (!requireNamespace(\"devtools\", quietly=TRUE))\n",
    "    install.packages(\"devtools\")\n",
    " \n",
    "#BiocManager::install(\"IRanges\")\n",
    "devtools::install_github(\"thackl/thacklr\")\n",
    "devtools::install_github(\"thackl/gggenomes\")\n",
    "devtools::install_github(\"YuLab-SMU/ggmsa\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdee64d3-e8de-472d-99c4-40d0d7db2026",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "library(tidyverse)\n",
    "library(gggenomes)\n",
    "library(patchwork)  # arrange multiple plots\n",
    "library(ggtree) \n",
    "library(\"phangorn\")\n",
    "library(IRanges)\n",
    "library(ggmsa)\n",
    "\n",
    "library(\"svglite\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ad2b2f9-637f-4587-a8a3-f11e015d30a6",
   "metadata": {},
   "source": [
    "## File Configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faee9184-86c7-4aca-b3bf-2ee0b1f20fa2",
   "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": "8806aa52-f729-4e15-8384-05b1f56a1141",
   "metadata": {},
   "outputs": [],
   "source": [
    "bgcflow_dir = Path(notebook_configuration[\"bgcflow_dir\"])\n",
    "PROJECT_NAME = \"Erythreapeptin\"\n",
    "PROJECT_NAME2 = \"mq_saccharopolyspora\"\n",
    "\n",
    "report_dir = bgcflow_dir / f\"data/processed/{PROJECT_NAME2}\"\n",
    "PROJECT_CONFIG_DIR = bgcflow_dir / f\"config/{PROJECT_NAME}\"\n",
    "\n",
    "SAMPLE_FILE_NAME = \"samples.csv\"\n",
    "ANTISMASH_VERSION = \"6.1.1\"\n",
    "FIGURE = \"Figure_6\"\n",
    "FIGURE_N = 'b'\n",
    "FIGURE_MASH = \"Figure_3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eac7414b-a92b-44ff-ba27-204a85cf7a88",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_config = pd.read_csv(PROJECT_CONFIG_DIR / SAMPLE_FILE_NAME)\n",
    "samples_config"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4af49c6-930b-4e88-bfc0-65bd61cee868",
   "metadata": {},
   "source": [
    "### Extract COG information from mmseqs2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f666754b-487e-4ad0-907a-dffb39b23f33",
   "metadata": {},
   "outputs": [],
   "source": [
    "mmseq_cluster = bgcflow_dir / f\"data/interim/mmseqs2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}_cluster.tsv\"\n",
    "assert mmseq_cluster.is_file(), f\"Make sure to run the subworkflow for project {PROJECT_NAME} by running: bgcflow run --snakefile workflow/BGC\"\n",
    "G = nx.read_edgelist(mmseq_cluster)\n",
    "# get the connected components of G\n",
    "components = nx.connected_components(G)\n",
    "\n",
    "# sort the components based on their size\n",
    "sorted_components = sorted(components, key=lambda x: len(x), reverse=True)\n",
    "\n",
    "result = {\"feature_id\": [], \"cluster_id\" : [], \"cluster_n\" : [], \"locus_tag\" : []}\n",
    "for num, g in enumerate(sorted_components):\n",
    "    size = len(g)\n",
    "    for item in g:\n",
    "        result['feature_id'].append(f\"cds-{item}\")\n",
    "        result['cluster_id'].append(f\"cog_{num+1:02d}\")\n",
    "        result['cluster_n'].append(size)\n",
    "        result['locus_tag'].append(item)\n",
    "df_mmseqs2 = pd.DataFrame.from_dict(result)\n",
    "df_mmseqs2.to_csv(f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_mmseqs_cog.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1018b815-461a-43b0-b51b-1bc633d8ac92",
   "metadata": {},
   "source": [
    "### Preparation of the dataset using minimap2 and mmseqs\n",
    "This analysis refers to the example shown here: https://thackl.github.io/gggenomes/articles/emales.html\n",
    "\n",
    "Here, COGs were defined using mmseqs2 and links between genes are create using minimap2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96c41fd6-8307-4aad-a818-2c30622d2bcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "region_seq_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.fna\")\n",
    "region_feat_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.gbk\")\n",
    "region_minimap2_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.paf\")\n",
    "cogs_raw = f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_mmseqs_cog.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c11bc6c1-b666-4d1f-8739-70ccad830df3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i region_seq_raw,region_feat_raw,region_minimap2_raw,cogs_raw\n",
    "region_seq <- read_seqs(region_seq_raw)\n",
    "region_feat <- read_feats(region_feat_raw)\n",
    "region_link <- read_paf(region_minimap2_raw)\n",
    "cogs <- read_csv(cogs_raw)\n",
    "\n",
    "cogs %<>% mutate(\n",
    "  cluster_label = paste0(cluster_id, \" (\", cluster_n, \")\"),\n",
    "  cluster_label = fct_lump_min(cluster_label, 5, other_level = \"rare\"),\n",
    "  cluster_label = fct_lump_min(cluster_label, 15, other_level = \"medium\"),\n",
    "  cluster_label = fct_relevel(cluster_label, \"rare\", after=Inf))\n",
    "\n",
    "merged_feat <- merge(region_feat, cogs, by.x = \"feat_id\", by.y = \"feature_id\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "250465a4-8b36-49a8-91ee-1fbf650afe78",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggenomes_feat_table = f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_gggenomes_feat.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f51a328f-a722-4450-950d-b7f8b5df2f76",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i ggenomes_feat_table\n",
    "write_csv(merged_feat, ggenomes_feat_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9460d435-4aaa-4983-9582-70f679b2598b",
   "metadata": {},
   "source": [
    "### GG Genomes Vis\n",
    "This is the default GG Genomes visualization using mmseqs2 as links, limited to the first 25000 bp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff33aa3f-95bd-48ab-af1d-25a45fc812af",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R # Manipulate length of the gene clusters\n",
    "#region_seq['length'][region_seq['seq_id'] == merged_feat[353+1,'seq_id']] <- merged_feat[353+1,'end']\n",
    "#region_seq['length'][region_seq['seq_id'] == merged_feat[86,'seq_id']] <- merged_feat[86,'end']\n",
    "#region_seq['length'][region_seq['seq_id'] == 'NZ_CP061007.1.region011'] <- 25000\n",
    "#region_seq['length'][region_seq['seq_id'] == 'NZ_GL877879.1.region006'] <- 25000\n",
    "#region_seq['length'][region_seq['seq_id'] == 'NZ_JACHIW010000001.1.region012'] <- 25000\n",
    "#region_seq['length'][region_seq['seq_id'] == 'NZ_CP040605.1.region009'] <- 25000\n",
    "region_seq['length'] <- 25000\n",
    "region_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "123184d1-cd87-4ee2-ab55-9daa5fb1d49c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 1200 -h 600\n",
    "p <- gggenomes(\n",
    "    genes=merged_feat,  # a gene track, added as first feat track\n",
    "    seqs=region_seq) %>% add_links(region_link) %>% flip_seqs(1,3,7,9,11,12,15)\n",
    "\n",
    "p2_mmseqs2 <- p +\n",
    "    geom_seq(aes(color=strand), arrow=TRUE) + \n",
    "    geom_seq() + geom_bin_label() +\n",
    "    geom_link(offset = c(0.2, 0.2), color=\"white\", alpha=.8) +  # the first link track\n",
    "    geom_gene(position = \"identity\", aes(fill=cluster_label)) +  # the first feat track filtered for geneish feats: CDS, mRNA, ..\n",
    "    labs(caption=\"Nucleotide Size\") + \n",
    "    scale_fill_manual(\"Conserved genes\", \n",
    "                      labels=sort(unique(merged_feat$cluster_label)),\n",
    "                      values=c(\"#f94144\",\"#f3722c\",\"#f8961e\",\"#f9844a\",\"#f9c74f\",\"#90be6d\",\"#43aa8b\",\"#4d908e\",\"#577590\",\"#277da1\", \n",
    "                               'white', 'white', 'white', 'white'))\n",
    "p2_mmseqs2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5910c3ca-632e-4523-a1f9-8f000de51919",
   "metadata": {},
   "source": [
    "## Utilize Clinker for Links\n",
    "Now, we will use clinker instead of mmseqs2 for the links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2764358d-d05e-44bd-861d-8e447942737a",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(bgcflow_dir / f\"data/processed/{PROJECT_NAME}/clinker/{ANTISMASH_VERSION}/clinker.txt\", \"r\") as f:\n",
    "    clinker_data = f.readlines()\n",
    "\n",
    "clinker_dict = {}\n",
    "container = {}\n",
    "for num, line in enumerate(clinker_data):\n",
    "    if line.startswith(\"----\"):\n",
    "        header = clinker_data[num - 1].strip(\"\\n\")\n",
    "        if header not in clinker_dict.keys():\n",
    "            clinker_dict[header] = []\n",
    "    else:\n",
    "        if line.startswith(\"----\") or line.startswith(\"Query\") or \"vs\" in line:\n",
    "            pass\n",
    "        else:\n",
    "            data = line.strip(\"\\n\").split()\n",
    "            seq_id = header.split(\" vs \")\n",
    "            if len(data) == 0:\n",
    "                data = [None, None, None, None]\n",
    "            data = seq_id + data\n",
    "            assert len(data) == 6, data\n",
    "            clinker_dict[header].append(data)\n",
    "            container[num] = data\n",
    "            \n",
    "\n",
    "df_clinker = pd.DataFrame.from_dict(container, orient=\"index\", columns=[\"seq_id\", \"seq_id2\", \"locus_tag\", \"locus_tag2\", \"identity\", \"similarity\"])\n",
    "df_clinker.identity = df_clinker.identity.astype(float)\n",
    "df_clinker = df_clinker[df_clinker.identity >= 0.5]\n",
    "df_clinker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "694b6895-e736-4809-8914-76c6552dfdc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "cds_dict = {}\n",
    "for i in samples_config.index:\n",
    "    gbk = Path(\"../../../../\") / samples_config.loc[i, \"gbk_path\"]\n",
    "    seq_id = samples_config.loc[i, \"bgc_id\"]\n",
    "    print(seq_id)\n",
    "    with open(gbk, \"r\") as input_handle:\n",
    "        for record in SeqIO.parse(input_handle, \"genbank\"):\n",
    "            antismash_data = record.annotations['structured_comment']['antiSMASH-Data']\n",
    "            for feature in record.features:\n",
    "                if feature.type == \"CDS\":\n",
    "                    translation = feature.qualifiers[\"translation\"][0]  \n",
    "                    strand = feature.location.strand\n",
    "                    if strand > 0:\n",
    "                        start, end = int(feature.location.start), int(feature.location.end)\n",
    "                    elif strand < 0:\n",
    "                        end, start = int(feature.location.start), int(feature.location.end)\n",
    "                    try:\n",
    "                        locus_tag = feature.qualifiers[\"locus_tag\"][0]\n",
    "                    except KeyError:\n",
    "                        print(seq_id, feature.qualifiers.keys())\n",
    "                        locus_tag = feature.qualifiers[\"protein_id\"][0]\n",
    "                        print(locus_tag)\n",
    "                    if seq_id.startswith(\"BGC\"):\n",
    "                        locus_tag = feature.qualifiers[\"protein_id\"][0]\n",
    "                    cds_dict[locus_tag] = {\"seq_id\" : seq_id, \"start\" : start, \"end\" : end, \"strand\" : strand, \"translation\":translation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f559ec6-10a8-41bf-9b06-6e62fdad855f",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in df_clinker.index:\n",
    "    locus_tag1 = df_clinker.loc[i, \"locus_tag\"]\n",
    "    locus_tag2 = df_clinker.loc[i, \"locus_tag2\"]\n",
    "    try:\n",
    "        if locus_tag1 != None:\n",
    "            df_clinker.loc[i, \"start\"] = int(cds_dict[locus_tag1]['start'])\n",
    "            df_clinker.loc[i, \"end\"] = int(cds_dict[locus_tag1]['end'])\n",
    "    except KeyError as e:\n",
    "        print(e, df_clinker.loc[i, \"seq_id\"])\n",
    "    try:\n",
    "        if locus_tag2 != None:\n",
    "            df_clinker.loc[i, \"start2\"] = int(cds_dict[locus_tag2]['start'])\n",
    "            df_clinker.loc[i, \"end2\"] = int(cds_dict[locus_tag2]['end'])\n",
    "    except KeyError as e:\n",
    "        print(e, df_clinker.loc[i, \"seq_id2\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71c0f578-6cbe-427f-8598-c9e1362da649",
   "metadata": {},
   "outputs": [],
   "source": [
    "region_clinker_raw = f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_clinker_links.csv\"\n",
    "df_clinker.dropna().to_csv(region_clinker_raw, index=False)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "c1358ca6-965b-44a7-b6b7-cbc17033fa0e",
   "metadata": {},
   "source": [
    "region_seq_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.fna\")\n",
    "region_feat_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.gbk\")\n",
    "region_minimap2_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.paf\")\n",
    "cogs_raw = f\"assets/tables/{PROJECT_NAME}_mmseqs_cog.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec4d5b0e-ed69-416b-a1b3-1eb5c1af08ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i region_seq_raw,region_feat_raw,region_clinker_raw,cogs_raw\n",
    "region_seq <- read_seqs(region_seq_raw)\n",
    "region_feat <- read_feats(region_feat_raw)\n",
    "region_link <- read_csv(region_clinker_raw)\n",
    "cogs <- read_csv(cogs_raw)\n",
    "\n",
    "cogs %<>% mutate(\n",
    "  cluster_label = paste0(cluster_id, \" (\", cluster_n, \")\"),\n",
    "  cluster_label = fct_lump_min(cluster_label, 5, other_level = \"rare\"),\n",
    "  cluster_label = fct_lump_min(cluster_label, 15, other_level = \"medium\"),\n",
    "  cluster_label = fct_relevel(cluster_label, \"rare\", after=Inf))\n",
    "\n",
    "merged_feat <- merge(region_feat, cogs, by.x = \"feat_id\", by.y = \"feature_id\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0285722d-7044-461a-9686-cc77435fffc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 1200 -h 600\n",
    "region_seq['length'] <- 20000\n",
    "\n",
    "p <- gggenomes(\n",
    "    genes=merged_feat,  # a gene track, added as first feat track\n",
    "    seqs=region_seq) %>% add_links(region_link) #%>% flip_seqs(1,3,7,9,11,12,15)\n",
    "\n",
    "p2_clinker <- p + geom_seq(aes(color=strand), arrow=TRUE) + \n",
    "    geom_seq() + geom_bin_label() +\n",
    "    geom_link() +  # the first link track\n",
    "    geom_gene(position = \"identity\", aes(fill=cluster_label)) +  # the first feat track filtered for geneish feats: CDS, mRNA, ..\n",
    "    labs(caption=\"Nucleotide Size\") +\n",
    "    scale_fill_manual(\"Conserved genes\", \n",
    "                      labels=sort(unique(merged_feat$cluster_label)),\n",
    "                      values=c(\"#f94144\",\"#f3722c\",\"#f8961e\",\"#f9844a\",\"#f9c74f\",\"#90be6d\",\"#43aa8b\",\"#4d908e\",\"#577590\",\"#277da1\", \n",
    "                               'white', 'white', 'white', 'white'))\n",
    "p2_clinker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd4fdcee-77d8-4232-aa5a-6556e7a737ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 1800 -h 600\n",
    "p2_clinker + p2_mmseqs2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9050232d-f648-4b7d-9557-f5140ee74d25",
   "metadata": {},
   "source": [
    "### Adding BGC tree based on multiple alignment of selected COGs"
   ]
  },
  {
   "cell_type": "raw",
   "id": "3b827726-a370-4136-bf72-357b2a34ceba",
   "metadata": {},
   "source": [
    "%%bash\n",
    "mamba install -c bioconda clustalo iqtree seqkit trimal -y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccda521e-ffb3-4d5e-ba37-3f31a7635df4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_gggenomes_feat.csv\")\n",
    "correct_naming = df.set_index('locus_tag.y').loc[:, [\"seq_id\", 'cluster_id']].T.to_dict()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "456b6228-6e70-42cc-8097-677e7564c176",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate tree\n",
    "for cog in df.cluster_label.unique():\n",
    "    if cog not in ['rare', 'medium']:\n",
    "        subset = df[df.cluster_label == cog]\n",
    "        product = subset['product'].unique()\n",
    "        print(cog, \"\\n* \".join(product))\n",
    "        print(\"\")\n",
    "        \n",
    "print(\"Total Seq:\", len(df.seq_id.unique()))\n",
    "number_of_bgcs = len(df.seq_id.unique()) - 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4e4f6d7-6c5c-47e9-adde-f23f3c0b16fb",
   "metadata": {},
   "source": [
    "From COG distribution, we can see that only COG_02 match criteria to build MSA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7796a74e-ff8d-44cf-acde-dc4697199510",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(bgcflow_dir / f\"data/interim/mmseqs2/{PROJECT_NAME}/msa_{PROJECT_NAME}_{ANTISMASH_VERSION}.db\") as f:\n",
    "    data = f.readlines()\n",
    "    for num, i in enumerate(data):\n",
    "        if i.startswith(\">\"):\n",
    "            rename = i.split()\n",
    "            locus_tag = rename[0].replace(\">\", \"\")\n",
    "            seq_id = correct_naming[locus_tag]['seq_id']\n",
    "            data[num] = \" \".join([f'>{seq_id}', locus_tag, \"\\n\"])\n",
    "\n",
    "            \n",
    "\n",
    "for num, i in enumerate(data):\n",
    "    if (\"#cl\" in i):\n",
    "        count = [int(x.strip(\"n=\")) for x in i.split() if \"n=\" in x][0]\n",
    "        if count in [15]:\n",
    "            print(i)\n",
    "            msa = data[(num+2):(num+2)+(count*2)]\n",
    "            filename = i.split(\"|\")[0].replace(\"#\", \"\").replace('\\0', '')\n",
    "            filename = Path(f'assets/data/{FIGURE}/msa/{PROJECT_NAME}/{filename}_aligned.fasta')\n",
    "            filename.parent.mkdir(exist_ok=True, parents=True)\n",
    "            print(filename)\n",
    "            with open(filename, \"w\") as outfile:\n",
    "                outfile.writelines(msa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49cb852e-456e-4f42-9d00-9a60666a1f5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "treefile = f\"assets/data/{FIGURE}/msa/{PROJECT_NAME}/core_BGC_concat_clean.fasta.treefile\"\n",
    "Path(treefile).parent.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6679a10-d12d-457c-b575-1544c740a1a4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%bash -s \"$PROJECT_NAME\" \"$FIGURE\" \"$treefile\" # USE FOR MULTIPLE GENES\n",
    "\n",
    "if [ -f $3 ]; then\n",
    "    echo \"File $3 already exists. Skipping command.\"\n",
    "else\n",
    "    echo \"Building Tree...\"\n",
    "    seqkit concat -o assets/data/$2/msa/$1/core_BGC_concat.fasta assets/data/$2/msa/$1/*_aligned.fasta\n",
    "    trimal -in assets/data/$2/msa/$1/core_BGC_concat.fasta -out assets/data/$2/msa/$1/core_BGC_concat_clean.fasta\n",
    "    iqtree -s assets/data/$2/msa/$1/core_BGC_concat_clean.fasta --redo\n",
    "fi"
   ]
  },
  {
   "cell_type": "raw",
   "id": "407cd3de-a085-4348-abd6-b819ce8b9537",
   "metadata": {},
   "source": [
    "%%bash -s \"$PROJECT_NAME\" # USE FOR SINGLE GENE\n",
    "#seqkit concat -o assets/data/$2/msa/$1/core_BGC_concat.fasta assets/data/$2/msa/$1/*_aligned.fasta\n",
    "cp assets/data/$2/msa/$1/*_aligned.fasta assets/data/$2/msa/$1/core_BGC_concat.fasta\n",
    "trimal -in assets/data/$2/msa/$1/core_BGC_concat.fasta -out assets/data/$2/msa/$1/core_BGC_concat_clean.fasta\n",
    "iqtree -s assets/data/$2/msa/$1/core_BGC_concat_clean.fasta --redo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a47397b1-087f-46a2-a2cf-cdb6049199e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare data tables\n",
    "gtdtbk = bgcflow_dir / f\"data/processed/{PROJECT_NAME2}/tables/df_gtdb_meta.csv\"\n",
    "mash = Path(f\"assets/tables/{FIGURE_MASH}b_mash_hcluster.csv\")\n",
    "ncbi = report_dir / \"tables/df_ncbi_meta.csv\"\n",
    "\n",
    "# merge data tables for tree\n",
    "df_mash = pd.read_csv(mash, index_col=0)\n",
    "df_gtdbtk = pd.read_csv(gtdtbk, index_col=0)\n",
    "df_gtdbtk = df_gtdbtk.loc[:, [\"Domain\", \"Phylum\", \"Class\", \"Order\", \"Family\", \"Genus\", \"Species\", \"Organism\"]]\n",
    "df_tree_annotation = pd.concat([df_mash, df_gtdbtk], axis=1, join=\"inner\")\n",
    "df_ncbi = pd.read_csv(ncbi, index_col=0).loc[:, \"strain\"]\n",
    "df_tree_annotation = pd.concat([df_tree_annotation, df_ncbi], axis=1, join=\"inner\")\n",
    "\n",
    "## manual annotation for tree\n",
    "### Rename MASH based Species Phylogroup\n",
    "# df[\"Species_Phylogroup\"] = [f\"P{i + 1}\" for i in df.hcluster]\n",
    "df_tree_annotation = df_tree_annotation#.reset_index(drop=False).rename(columns={\"index\" : \"genome_id\"})\n",
    "tree_mapping = df_tree_annotation.T.to_dict()\n",
    "df_tree_annotation.columns\n",
    "df_tree_annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9666f748-ddcc-4a52-a7ba-360d7dc04463",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_config = pd.read_csv(PROJECT_CONFIG_DIR / SAMPLE_FILE_NAME)\n",
    "for i in samples_config.index:\n",
    "    genome_id = samples_config.loc[i, \"genome_id\"]\n",
    "    for c in df_tree_annotation.columns:\n",
    "        if genome_id in tree_mapping.keys():\n",
    "            samples_config.loc[i, c] = tree_mapping[genome_id][c]"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b7003176-c412-4044-9ad0-64d9bfeb16e5",
   "metadata": {},
   "source": [
    "# Renaming of some strains based on redundancy\n",
    "df.loc[9,\"strain\"] = \"NRRL 30141\"\n",
    "df.loc[8,\"strain\"] = \"DSM 44228\"\n",
    "df.loc[22,\"strain\"] = \"DSM 40517\"\n",
    "df.loc[24,\"strain\"] = \"DSM 40517\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc9667eb-4d82-4921-917a-61d130635283",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Rename MIBIG entries\n",
    "for i in samples_config.index:\n",
    "    if samples_config.loc[i, \"bgc_id\"].startswith(\"BGC\"):\n",
    "        samples_config.loc[i, \"bgc_id\"] = samples_config.loc[i, \"bgc_id\"] + \".region001\"\n",
    "        samples_config.loc[i, \"phylogroup\"] = \"MIBIG\"\n",
    "    elif samples_config.loc[i, \"bgc_id\"] == \"NC_000962.3.region003\":\n",
    "        samples_config.loc[i, \"phylogroup\"] = \"M. tuberculosis\"\n",
    "        samples_config.loc[i, \"color_code\"] = \"red\"\n",
    "### Set up tip labels to show\n",
    "#samples_config['tip_label'] = [f\"{samples_config.loc[i, 'phylogroup']} | {samples_config.loc[i, 'bgc_id']} | {samples_config.loc[i, 'Organism'].replace('s__Saccharopolyspora','S')} ({samples_config.loc[i, 'strain']})\" for i in samples_config.index]\n",
    "samples_config['tip_label'] = [f\"{samples_config.loc[i, 'phylogroup']} | {samples_config.loc[i, 'bgc_id']}\" for i in samples_config.index]\n",
    "\n",
    "# save to intermediate file\n",
    "tree_annotation_file = f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_ggtree.csv\"\n",
    "samples_config.to_csv(tree_annotation_file, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9fe5a45-53bc-47df-bd1d-532605ae6ae1",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3551927-0179-422a-a76e-0b0e778c2c5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 800 -h 600 -i treefile,tree_annotation_file\n",
    "tree <- read.tree(treefile)\n",
    "tree_data <- read.csv(tree_annotation_file)\n",
    "\n",
    "# midpoint root\n",
    "tree <- phangorn::midpoint(tree)\n",
    "tree <- ladderize(reorder(tree))\n",
    "\n",
    "t <- ggtree(tree)\n",
    "t <- t %<+% tree_data\n",
    "\n",
    "t2 = t + geom_tiplab(aes(label=tip_label, fill=phylogroup),\n",
    "                     size=3, hjust=-0.03, family='sans', #align=T,\n",
    "                    linetype = \"dotted\", linesize = 1) + hexpand(1.5)\n",
    "    \n",
    "t3 = t2 + geom_tippoint(size=2.6, alpha=0.8, aes(color=phylogroup, stroke=1)) + \n",
    "        scale_color_manual(values = c('P1'='#264653', 'P2'='#e9c46a', 'P3'='#808080', \n",
    "                           'P4'='#808080', 'P5'='#f4a261', 'P6'='#808080', \n",
    "                           'P7'='#e76f51', 'P8'='#2a9d8f', 'MIBIG'='blue')) + \n",
    "        theme(legend.position='none')\n",
    "\n",
    "fig <- t3 + scale_y_continuous(expand=c(0.01,0.7,0.01,0.7)) + p2_clinker %>% pick_by_tree(t2) + plot_layout(widths = c(2,3))\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62d3df39-3db6-4329-8a07-ab144aa6c4ff",
   "metadata": {},
   "source": [
    "### Add mapping to known reference using CBlaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a38b92e1-6d93-4f76-b132-cb9d2bfd3e54",
   "metadata": {},
   "outputs": [],
   "source": [
    "cblaster_db_path = str(bgcflow_dir / \"data/processed/mq_saccharopolyspora/cblaster/bgcs/6.1.1/cblaster_bgc_db.dmnd\")\n",
    "query_gbk = bgcflow_dir / \"data/external/bgc_selection/MIBIG/BGC0000513.region001.gbk\"\n",
    "cblaster_out = Path(f\"assets/tables/{FIGURE}{FIGURE_N}_cblaster_{PROJECT_NAME}_{query_gbk.stem}.csv\")\n",
    "cblaster_out.parent.mkdir(parents=True, exist_ok=True)\n",
    "query_gbk = str(query_gbk)\n",
    "cblaster_out = str(cblaster_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46ffdd31-245d-4f3f-8748-31306168b297",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$cblaster_db_path\" \"$query_gbk\" \"$cblaster_out\"\n",
    "# first, create env using mamba \n",
    "# mamba create -f <your bgcflow dir>/workflow/envs/cblaster.yaml\n",
    "mamba run -n cblaster cblaster search -m local -db $1 -qf $2 -o $3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "523a7d90-2c61-4cf8-89ec-28a324371920",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(cblaster_out, \"r\") as f:\n",
    "    cblaster = f.readlines()\n",
    "\n",
    "cluster_marker = []    \n",
    "for num, line in enumerate(cblaster):\n",
    "    if line.startswith(\"Cluster\"):\n",
    "        cluster_marker.append(num)\n",
    "container = {'Query' : [], 'Subject' : [], 'Identity' : [], 'Coverage' : [], 'E-value' : [], 'Bitscore' : [], 'Start' : [], 'End' : [], 'Strand' : []}\n",
    "for num, item in enumerate(cluster_marker):\n",
    "    if num == len(cluster_marker)-1:\n",
    "        data = [x.split() for x in cblaster[item+2:]]\n",
    "    else:\n",
    "        data = [x.split() for x in cblaster[item+2:cluster_marker[num+1]-6]]\n",
    "    for v in data:\n",
    "        for n, c in enumerate(container.keys()):\n",
    "            container[c].append(v[n])\n",
    "df_feat = pd.read_parquet(report_dir / \"data_warehouse/6.1.1/cdss.parquet\")\n",
    "seq_id_mapping = df_feat.set_index(\"locus_tag\").region_id.to_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa1e6289-c16c-4d59-bc2d-4b1256fe1446",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cblaster_hits = pd.DataFrame.from_dict(container)\n",
    "df_cblaster_hits.Query.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72c3e820-ea2e-41a6-ba5c-36f906d387be",
   "metadata": {},
   "outputs": [],
   "source": [
    "cblaster_cds_dict = {}\n",
    "seq_id = Path(query_gbk).stem\n",
    "with open(query_gbk, \"r\") as input_handle:\n",
    "    for record in SeqIO.parse(input_handle, \"genbank\"):\n",
    "        antismash_data = record.annotations['structured_comment']['antiSMASH-Data']\n",
    "        for feature in record.features:\n",
    "            if feature.type == \"CDS\":\n",
    "                start, end = int(feature.location.start), int(feature.location.end)\n",
    "                strand = feature.location.strand\n",
    "                translation = feature.qualifiers[\"translation\"][0]\n",
    "                try:\n",
    "                    gene = feature.qualifiers['gene'][0]\n",
    "                except KeyError:\n",
    "                    #print(feature.qualifiers)\n",
    "                    gene = feature.qualifiers['product'][0]\n",
    "                try:\n",
    "                    locus_tag = feature.qualifiers[\"locus_tag\"][0]\n",
    "                    #print(feature.qualifiers)\n",
    "                except KeyError:\n",
    "                        #print(seq_id, feature.qualifiers.keys())\n",
    "                        locus_tag = feature.qualifiers[\"protein_id\"][0]\n",
    "                        print(locus_tag)\n",
    "                #if seq_id.startswith(\"BGC\"):\n",
    "                #        locus_tag = feature.qualifiers[\"protein_id\"][0]\n",
    "                cblaster_cds_dict[locus_tag] = {\"seq_id\" : seq_id, \"start\" : start, \"end\" : end, \n",
    "                                                \"strand\" : strand, \"translation\":translation,\n",
    "                                               \"name\" : gene}\n",
    "df_cds_cblaster = pd.DataFrame.from_dict(cblaster_cds_dict).T\n",
    "df_cds_cblaster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddd59f35-30ae-4730-aba1-cba0bff9e1f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_query_mapping = df_cds_cblaster.loc[df_cblaster_hits.Query.unique()]\n",
    "q_map = df_query_mapping.loc[:, ['name', 'seq_id']].T.to_dict()\n",
    "q_map['SACE_4230']['name'] = 'EryKC'\n",
    "q_map['SACE_4231']['name'] = 'EryS'\n",
    "df_cblaster_hits['query_description'] = [q_map[i]['name'] for i in df_cblaster_hits.Query]\n",
    "df_cblaster_hits['seq_id'] = [seq_id_mapping[i] for i in df_cblaster_hits.Subject]\n",
    "#df_cblaster_hits['feat_id'] = [f\"cds-{i}\" for i in df_cblaster_hits.Subject]\n",
    "df_cblaster_hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68a82797-e9c8-4ada-9df1-d1ffeee29e15",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cog = pd.read_csv(f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_mmseqs_cog.csv\")\n",
    "df_cog[df_cog.locus_tag.isin(df_cblaster_hits[\"Query\"].unique())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba8e343f-c2dc-4cd3-a3a0-5989d68fb24f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ctr = len(df_cblaster_hits)\n",
    "description_map = df_cblaster_hits.set_index(\"Query\").loc[:, \"query_description\"].to_dict()\n",
    "for query_self in df_cblaster_hits.Query.unique():\n",
    "    print(query_self)\n",
    "    df_cblaster_hits.loc[ctr, \"Query\"] = query_self\n",
    "    df_cblaster_hits.loc[ctr, \"Subject\"] = query_self\n",
    "    df_cblaster_hits.loc[ctr, \"Identity\"] = 100\n",
    "    df_cblaster_hits.loc[ctr, \"Coverage\"] = 100\n",
    "    df_cblaster_hits.loc[ctr, \"E-value\"] = 0\n",
    "    df_cblaster_hits.loc[ctr, \"Bitscore\"] = 999\n",
    "    df_cblaster_hits.loc[ctr, \"Start\"] = cblaster_cds_dict[query_self]['start']\n",
    "    df_cblaster_hits.loc[ctr, \"End\"] = cblaster_cds_dict[query_self]['end'] - 10\n",
    "    if cblaster_cds_dict[query_self]['strand'] > 0:\n",
    "        df_cblaster_hits.loc[ctr, \"Strand\"] = \"+\"\n",
    "    else:\n",
    "        df_cblaster_hits.loc[ctr, \"Strand\"] = \"-\"\n",
    "    df_cblaster_hits.loc[ctr, \"seq_id\"] = cblaster_cds_dict[query_self]['seq_id']\n",
    "    df_cblaster_hits.loc[ctr, \"query_description\"] = description_map[query_self]\n",
    "    ctr = ctr + 1 \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "617abbc9-7ce2-48e4-9e6b-b0744c0b5b1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cog = pd.read_csv(f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_mmseqs_cog.csv\")\n",
    "df_cog = df_cog.merge(df_cblaster_hits, left_on=\"locus_tag\", right_on=\"Subject\", how=\"outer\")\n",
    "df_cog = df_cog.rename(columns={'feature_id':'feat_id', 'Start':'start', 'End':'end', 'region_id':'seq_id', 'Query':'blast_hit', 'query_description':'blast_desc'})  \n",
    "df_cog = df_cog.dropna().loc[:, ['feat_id', 'start', 'end', 'seq_id', 'blast_hit', 'blast_desc', 'cluster_id', 'Identity', 'Coverage', 'E-value', 'Bitscore']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0f0325e-3aba-4378-bca7-741e230413ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cog['end'] = df_cog['end'].astype(int) - df_cog['start'].astype(int)\n",
    "df_cog['start'] = 0\n",
    "blast_hits = f\"assets/tables/{FIGURE}{FIGURE_N}_cblaster_{PROJECT_NAME}_blast_hits.csv\"\n",
    "df_cog.to_csv(blast_hits, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4cfb6b9-8a1b-4272-a2c5-7bd8c9772f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i blast_hits\n",
    "blast <- read_csv(blast_hits)\n",
    "blast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9f5cfe5-68d7-4d83-8c93-e599dad8af9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "region_seq_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.fna\")\n",
    "region_feat_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.gbk\")\n",
    "region_minimap2_raw = str(bgcflow_dir / f\"data/interim/minimap2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}.paf\")\n",
    "cogs_raw = f\"assets/tables/{FIGURE}{FIGURE_N}_{PROJECT_NAME}_mmseqs_cog.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d464b885-4a3e-40b2-b65d-39371699a847",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -i region_seq_raw,region_feat_raw,region_minimap2_raw,cogs_raw\n",
    "region_seq <- read_seqs(region_seq_raw)\n",
    "region_feat <- read_feats(region_feat_raw)\n",
    "region_link <- read_paf(region_minimap2_raw)\n",
    "region_link <- read_csv(region_clinker_raw)\n",
    "cogs <- read_csv(cogs_raw)\n",
    "\n",
    "cogs %<>% mutate(\n",
    "  cluster_label = paste0(cluster_id, \" (\", cluster_n, \")\"),\n",
    "  cluster_label = fct_lump_min(cluster_label, 5, other_level = \"rare\"),\n",
    "  cluster_label = fct_lump_min(cluster_label, 15, other_level = \"medium\"),\n",
    "  cluster_label = fct_relevel(cluster_label, \"rare\", after=Inf))\n",
    "\n",
    "merged_feat <- merge(region_feat, cogs, by.x = \"feat_id\", by.y = \"feature_id\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be0b0250-135e-412d-b8c1-85a823e84971",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 1200 -h 600\n",
    "#region_seq['length'] <- 30000\n",
    "\n",
    "p <- gggenomes(\n",
    "    genes=merged_feat,  # a gene track, added as first feat track\n",
    "    seqs=region_seq) %>% add_links(region_link) %>% add_subfeats(blast, .transform = \"none\") %>% flip_seqs(3,5,12)\n",
    "\n",
    "p2_clinker <- p + geom_seq(aes(color=strand), arrow=TRUE) + \n",
    "    geom_seq() + #geom_bin_label() +\n",
    "    geom_link() +  # the first link track\n",
    "    geom_gene(position = \"identity\", aes(fill=cluster_label)) +  # the first feat track filtered for geneish feats: CDS, mRNA, ..\n",
    "    labs(caption=\"Nucleotide Size\") +\n",
    "    scale_fill_manual(\"Conserved genes\", \n",
    "                      labels=sort(unique(merged_feat$cluster_label)),\n",
    "                      values=c(\"#f94144\",\"#f3722c\",\"#f8961e\",\"#f9844a\",\"#f9c74f\",\"#90be6d\",\"#43aa8b\",\"#4d908e\",\"#577590\",\"#277da1\", \n",
    "                               'white', 'white', 'white', 'white'))\n",
    "p2_clinker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e34ca467-54d7-42ab-9a30-a0129db03afb",
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile_svg = f\"assets/figures/{FIGURE}/{FIGURE}c.svg\"\n",
    "outfile_pdf = f\"assets/figures/{FIGURE}/{FIGURE}c.pdf\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "028b78a8-d150-4103-87a7-2663d82cdc74",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 800 -h 600 -i treefile,tree_annotation_file,outfile_svg,outfile_pdf\n",
    "tree <- read.tree(treefile)\n",
    "tree_data <- read.csv(tree_annotation_file)\n",
    "\n",
    "# midpoint root\n",
    "tree <- phangorn::midpoint(tree)\n",
    "tree <- ladderize(reorder(tree))\n",
    "\n",
    "t <- ggtree(tree)\n",
    "t <- t %<+% tree_data\n",
    "\n",
    "t2 = t + geom_tiplab(aes(label=tip_label, fill=phylogroup),\n",
    "                     size=2.55, hjust=-0.03, family='sans', #align=T,\n",
    "                    linetype = \"dotted\", linesize = 1) + hexpand(2.5)\n",
    "    \n",
    "t3 = t2 + geom_tippoint(size=2.5, alpha=0.8, aes(color=phylogroup, stroke=1)) + \n",
    "        scale_color_manual(values = c('P1'='#264653', 'P2'='#e9c46a', 'P3'='#808080', \n",
    "                           'P4'='#808080', 'P5'='#f4a261', 'P6'='#808080', \n",
    "                           'P7'='#e76f51', 'P8'='#2a9d8f', 'MIBIG'='blue')) + \n",
    "        theme(legend.position='none')\n",
    "\n",
    "tree_ids <- t3$data %>% filter(isTip) %>% arrange(-y)\n",
    "label <- tree_ids$label\n",
    "label_with_missing <- c(label, \"NZ_JADDUE010000019.1.region001\")        \n",
    "        \n",
    "region_seq['length'][region_seq['seq_id'] == 'NZ_GL877879.1.region013'] <- 40000\n",
    "p <- gggenomes(\n",
    "    genes=merged_feat,  # a gene track, added as first feat track\n",
    "    seqs=region_seq) %>% add_links(region_link) %>% add_subfeats(blast, .transform = \"none\")\n",
    "\n",
    "\n",
    "p2 <- p +\n",
    "    geom_seq(aes(color=strand), arrow=TRUE) + \n",
    "    geom_seq() + #geom_bin_label(size = 3.5, expand_left =.7) + \n",
    "    geom_link(offset = c(0.2, 0.2), color=\"white\", alpha=.8) +  # the first link track\n",
    "    geom_gene(position = \"identity\", aes(fill=cluster_label)) +  # the first feat track filtered for geneish feats: CDS, mRNA, ..\n",
    "    labs(caption=\"Size (bp)\") + \n",
    "    geom_feat(aes(color=blast_desc), size=1, position=\"strandpile\") +\n",
    "    scale_fill_manual(\"Conserved genes\", \n",
    "                      labels=sort(unique(merged_feat$cluster_label)),\n",
    "                      values=c(\n",
    "                               \"#f94144\",\"#f3722c\",\"#f8961e\",#\"#f9844a\",\"#f9c74f\", \n",
    "                               #\"#90be6d\",\"#43aa8b\",\"#4d908e\",\"#577590\",\"#277da1\", \n",
    "                               'grey', 'white', 'white')) +\n",
    "    #scale_color_viridis_d(\"Blast hits & Features\", direction = -1) +\n",
    "    scale_linetype(\"Graphs\") +  theme(legend.position=\"none\")\n",
    "\n",
    "\n",
    "\n",
    "fig <- t3 + scale_y_continuous(expand=c(0.08,0.7,0.02,0.7)) + p2 %>% pick(label_with_missing) %>% flip_seqs(6,11,14) + plot_layout(widths = c(1,2))\n",
    "    \n",
    "#fig <- t3 + scale_y_continuous(expand=c(0.01,0.7,0.01,0.7)) + p2 %>% pick_by_tree(t2) + plot_layout(widths = c(2,2))\n",
    "ggsave(file=outfile_svg, plot=fig, device=svglite, width=2000, height=1600, units=\"px\")\n",
    "ggsave(plot=fig, width=1600, height=1200, units=\"px\", dpi=300, filename=outfile_pdf, useDingbats=FALSE)\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1752e898-0ee1-4e61-bc9f-40100e5b2ca6",
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile_svg = f\"assets/figures/{FIGURE}/{FIGURE}c_legend.svg\"\n",
    "outfile_pdf = f\"assets/figures/{FIGURE}/{FIGURE}c_legend.pdf\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d86b00d-9571-4e36-a91c-3bb0dd4afeb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 800 -h 600 -i treefile,tree_annotation_file,outfile_svg,outfile_pdf\n",
    "fig <- fig +  theme(legend.position=\"right\")\n",
    "ggsave(file=outfile_svg, plot=fig, device=svglite, width=2000, height=1600, units=\"px\")\n",
    "ggsave(plot=fig, width=1600, height=1200, units=\"px\", dpi=300, filename=outfile_pdf, useDingbats=FALSE)\n",
    "fig"
   ]
  },
  {
   "cell_type": "raw",
   "id": "a516b4b0-1aa7-41cf-96db-035a3949e441",
   "metadata": {},
   "source": [
    "• NZ_CP069353.1.region023\n",
    "• NZ_JADDUE010000019.1.region001\n",
    "• NZ_PDBV01000001.1.region025"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bde6969f-131e-440d-8107-3ec4355614c1",
   "metadata": {},
   "source": [
    "### MSA of precursor peptide\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc4df49-4b39-46b7-bf39-cafc23212553",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cds = pd.DataFrame.from_dict(cds_dict).T\n",
    "df_cds = pd.concat([df_cds, df_cds_cblaster])\n",
    "df_precursor_hits = df_cblaster_hits[df_cblaster_hits.query_description == \"EryS\"]\n",
    "df_precursor_peptide = df_cds[df_cds.index.isin(df_precursor_hits.Subject.unique())]\n",
    "df_precursor_peptide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80961d68-5726-445b-bd9e-72776d0690d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_precursor_peptide = df_precursor_peptide[df_precursor_peptide.translation#.str.contains(\"IDGMCGVY\")]\n",
    "fasta = []\n",
    "for i in df_precursor_peptide.index:\n",
    "    header = f'>{df_precursor_peptide.loc[i, \"seq_id\"]}'# {df_precursor.loc[i, \"locus_tag\"]} MftA'\n",
    "    fasta.append(header)\n",
    "    print(header, i)\n",
    "    fasta.append(df_precursor_peptide.loc[i, \"translation\"])\n",
    "    \n",
    "precursor_sequence = Path(f'assets/data/{FIGURE}/msa/{PROJECT_NAME}_precursor.fasta')\n",
    "precursor_sequence.parent.mkdir(parents=True, exist_ok=True)\n",
    "with open(precursor_sequence,'w') as file:\n",
    "    file.write('\\n'.join(fasta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bd6f7bf-cf69-4efe-aa03-0bc482d536b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "precursor_sequence_aligned = f'assets/data/{FIGURE}/msa/{PROJECT_NAME}_precursor_aligned.fasta'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52b5480d-401d-48d1-a605-ba3fcbd74410",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$precursor_sequence\" \"$precursor_sequence_aligned\"\n",
    "clustalo -i $1 -o $2 --force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f6f0425-571c-4b09-aabd-47389c31f1eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile_svg = f\"assets/figures/{FIGURE}/{FIGURE}d.svg\"\n",
    "outfile_pdf = f\"assets/figures/{FIGURE}/{FIGURE}d.pdf\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d54ce923-d1d4-42f6-9533-f43113497ef2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 800 -h 400 -i precursor_sequence_aligned -i outfile_pdf -i outfile_svg \n",
    "protein_sequences <- precursor_sequence_aligned\n",
    "fig <- ggmsa(protein_sequences, start = 10, end = 55, char_width = 0.5, seq_name = T, \n",
    "             color = \"Clustal\", by_conservation=F, show.legend = F) + geom_msaBar()\n",
    "ggsave(filename=outfile_svg, plot=fig, device=svglite, width=2400, height=1600, units=\"px\")\n",
    "ggsave(plot=fig, width=2400, height=1600, units=\"px\", dpi=300, filename=outfile_pdf, useDingbats=FALSE)\n",
    "fig"
   ]
  }
 ],
 "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
}