18_Fig7b_d_Ranthipeptide_gggenomes.R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b0f08259-b159-4c4a-8948-2bce9fcad16c",
   "metadata": {},
   "source": [
    "## Multiple Sequence Alignment of Mycofactocin candidate clusters\n",
    "This analysis requires project `TIGR03997` containing 17 BGCs that have match with ARTS model TIGR03997. The project was run through the `BGC` comparison subworkflow by running `bgcflow run --snakefile workflow/BGC`"
   ]
  },
  {
   "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",
    "import yaml\n",
    "\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "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\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6dcec69-0334-4065-87b4-bd85b4e97f0b",
   "metadata": {},
   "source": [
    "## File configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a72f31f-672e-4f34-8f17-a6bad8f3ddd6",
   "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": "65ad65bd-ad30-4ee4-8537-e4c04360eae8",
   "metadata": {},
   "outputs": [],
   "source": [
    "bgcflow_dir = Path(notebook_configuration[\"bgcflow_dir\"])\n",
    "PROJECT_NAME = \"TIGR03997_with_reference\"\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_7\"\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": [
    "G = nx.read_edgelist(bgcflow_dir / f\"data/interim/mmseqs2/{PROJECT_NAME}/{PROJECT_NAME}_{ANTISMASH_VERSION}_cluster.tsv\")\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 = bgcflow_dir / 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": {
    "scrolled": true
   },
   "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": "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()"
   ]
  },
  {
   "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(\"\")"
   ]
  },
  {
   "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",
    "number_of_bgcs = len(df.seq_id.unique())\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 [number_of_bgcs]:\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": "326ba854-12ec-454a-a44e-bf474e8099e9",
   "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": "33793fa3-fc9c-41ce-b129-8a50da97688f",
   "metadata": {},
   "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": "e7686ee2-42ae-4b1e-8bd5-f2af574e028f",
   "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": "ab3b17f7-1c96-4426-8919-7ad57bb24ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare data tables\n",
    "gtdtbk = bgcflow_dir / f\"data/processed/{PROJECT_NAME}/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": "49cb852e-456e-4f42-9d00-9a60666a1f5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "treefile = f\"assets/data/{FIGURE}/msa/{PROJECT_NAME}/core_BGC_concat_clean.fasta.treefile\""
   ]
  },
  {
   "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', 'M. tuberculosis'='red')) + \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,2))\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/processed/Mycobacterium/antismash/6.1.1/GCF_000195955.2/NC_000962.3.region003.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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26b305b1-850f-43d3-b20a-9866990a445b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_query_mapping = pd.read_parquet(bgcflow_dir / \"data/processed/Mycobacterium/data_warehouse/6.1.1/cdss.parquet\").set_index(\"locus_tag\").loc[df_cblaster_hits.Query.unique()]\n",
    "q_map = df_query_mapping.loc[:, ['name', 'region_id']].T.to_dict()\n",
    "q_map['JDNMFAPF_07350']['name'] = 'mftA'\n",
    "q_map['JDNMFAPF_07360']['name'] = 'mftB'\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": "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",
    "                        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",
    "                cblaster_cds_dict[locus_tag] = {\"seq_id\" : seq_id, \"start\" : start, \"end\" : end, \"strand\" : strand, \"translation\":translation}"
   ]
  },
  {
   "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\"] = cds_dict[query_self]['start']\n",
    "    df_cblaster_hits.loc[ctr, \"End\"] = cds_dict[query_self]['end'] - 10\n",
    "    if 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\"] = cds_dict[query_self]['seq_id']\n",
    "    df_cblaster_hits.loc[ctr, \"query_description\"] = description_map[query_self]\n",
    "    ctr = ctr + 1 \n",
    "df_cblaster_hits"
   ]
  },
  {
   "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']]\n",
    "df_cog"
   ]
  },
  {
   "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(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": "028b78a8-d150-4103-87a7-2663d82cdc74",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R -w 800 -h 600\n",
    "library(\"svglite\")\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'][region_seq['seq_id'] == 'NZ_CP040605.1.region009'] <- 25000\n",
    "#region_seq['length'] <- 5000\n",
    "region_seq\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(1,3,7,9,11,12,15)\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', 'grey', 'white', 'white')) +\n",
    "    #scale_color_viridis_d(\"Blast hits & Features\", direction = -1) +\n",
    "    scale_linetype(\"Graphs\") +  theme(legend.position=\"none\")\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=\"assets/figures/Figure_7/b.svg\", plot=fig, device=svglite, width=2800, height=1800, units=\"px\")\n",
    "ggsave(plot=fig, width=1600, height=1200, units=\"px\", dpi=300, filename=\"assets/figures/Figure_7/b.pdf\", useDingbats=FALSE)\n",
    "fig"
   ]
  },
  {
   "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_precursor_peptide = df_cds[df_cds.seq_id.isin(df_cog.seq_id.unique())]"
   ]
  },
  {
   "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 = f'assets/data/{FIGURE}/msa/{PROJECT_NAME}_precursor.fasta'\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": "ed6747d2-31be-4dd3-a665-e3847446cc75",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "if (!requireNamespace(\"devtools\", quietly=TRUE))\n",
    "    install.packages(\"devtools\")\n",
    "devtools::install_github(\"YuLab-SMU/ggmsa\")\n",
    "\n",
    "library(ggmsa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2de54661-ce1f-4bc4-a65f-56101b97ce70",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "#install.packages(\"svglite\")\n",
    "library(\"svglite\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5b1f84e-9375-4404-ad15-79b548aca5d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure_svg = f\"assets/figures/{FIGURE}/{FIGURE}c.svg\"\n",
    "figure_pdf = f\"assets/figures/{FIGURE}/{FIGURE}c.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,figure_svg,figure_pdf\n",
    "protein_sequences <- precursor_sequence_aligned\n",
    "fig <- ggmsa(protein_sequences, start = 0, end = 50, char_width = 0.5, seq_name = T, \n",
    "             color = \"Clustal\", by_conservation=F, show.legend = F) + geom_msaBar()\n",
    "ggsave(filename=figure_svg, plot=fig, device=svglite, width=2800, height=1600, units=\"px\")\n",
    "ggsave(plot=fig, width=2800, height=1600, units=\"px\", dpi=300, filename=figure_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
}