12b - ASRAnnotationAnalysis_Processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c6f1dbd7-530e-4b3e-b059-b571c0a2211b",
   "metadata": {},
   "source": [
    "This part of the pipeline processes the results from the Ancestral State Reconstruction using the Count tool, and annotates the exchanged genes using COG and BRITE."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "108cc388-7a8c-49e6-be2d-84c61c5e0f08",
   "metadata": {},
   "source": [
    "### Paths and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7416196b-0870-438e-b39c-87cfac4d4b60",
   "metadata": {},
   "source": [
    "#### Pipeline input folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4f9b652-e25c-4f5f-95e0-e0dad7608c75",
   "metadata": {},
   "outputs": [],
   "source": [
    "apa.file = \"families.tsv\"\n",
    "tree.file = \"input_ready.tree\"\n",
    "indices.path = \"02-QC/indices\"\n",
    "cog.file = \"matrix_annotations.tsv\"\n",
    "kog.file = \"07-PangenomeAnnotation/KOG/processed_output/merged_full_annot\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fcd289d-e2a4-4bd1-9243-d2225ee402d0",
   "metadata": {},
   "source": [
    "#### Pipeline output folders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c69fd64-81a3-414d-8bc0-deaa986d6611",
   "metadata": {},
   "outputs": [],
   "source": [
    "task_root = \"12-ASR-analysis\"\n",
    "ASR_output = paste(task_root, \"output\", sep = \"/\")\n",
    "\n",
    "system(paste('mkdir -p', task_root, ASR_output), intern = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcda2af3-2b63-43f1-bf24-4098f52f50f5",
   "metadata": {},
   "source": [
    "#### Tool pointers and parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d326062f-3640-489d-b89b-03b46bdc943e",
   "metadata": {},
   "outputs": [],
   "source": [
    "cog.cats.file = \"utils/COG_cats.tsv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98991973-8db8-43ee-b6f3-0023ab1b1a53",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "library(data.table)\n",
    "library(dplyr)\n",
    "library(tidyr)\n",
    "library(Matrix)\n",
    "library(ggplot2)\n",
    "library(ggh4x)\n",
    "library(RColorBrewer)\n",
    "library(pals)\n",
    "library(IRdisplay)\n",
    "library(ape)\n",
    "library(ggtree)\n",
    "library(stringr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7efb5f17-f54e-404a-b5b7-58be7c7768eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "getPalette = colorRampPalette(rev(stepped()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70ccff45-c93a-4620-b3be-c98f016dbce2",
   "metadata": {},
   "source": [
    "## Reading files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e568620-3006-41a7-b54f-ef79bb433885",
   "metadata": {},
   "outputs": [],
   "source": [
    "root = getwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53bbfcd2-cbc5-4d8c-9817-f8395c252077",
   "metadata": {},
   "outputs": [],
   "source": [
    "setwd(task_root)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69482100-a34c-44ea-8ea8-b504a5401718",
   "metadata": {},
   "source": [
    "#### Ancestral P/A data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74809cd4-ab98-4829-bb7b-9c4d2ab5d274",
   "metadata": {},
   "outputs": [],
   "source": [
    "apa = as.data.frame(fread(apa.file))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c9b0b00-03a1-42a5-89ef-048392bf3310",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(apa)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2d77d1f-ca6c-4fae-a206-64334535fe44",
   "metadata": {},
   "source": [
    "Genes that are gained in the root node (with label 1) are the LCA genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba8eca6c-de4a-448d-95d9-f1f4da43d4f4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "common = apa$'1' != 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df3399d3-ce88-458d-af01-43e0a0e8c061",
   "metadata": {},
   "source": [
    "#### Phylogeny"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "818725c6-6fa1-45d3-a95c-c1179fba0449",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = read.tree(tree.file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7343dd2-4d6a-4557-bfe6-27a5919369c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggtree(tree) + geom_text(aes(label = label), nudge_x = -0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f12e6755-08e0-4cd7-973c-bb210945f8e9",
   "metadata": {},
   "source": [
    "#### Genome indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36a7c15b-d296-4a57-b785-9a7b690651cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "index.1 = read.table(paste(paste(root, indices.path, sep = \"/\"), \"group1\", sep = \"/\"))$V1\n",
    "index.4 = read.table(paste(paste(root, indices.path, sep = \"/\"), \"group4\", sep = \"/\"))$V1\n",
    "index.14a = read.table(paste(paste(root, indices.path, sep = \"/\"), \"group14a\", sep = \"/\"))$V1\n",
    "index.14b = read.table(paste(paste(root, indices.path, sep = \"/\"), \"group14b\", sep = \"/\"))$V1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "553645a3-c6a3-4931-982c-7ba1d4d49d92",
   "metadata": {},
   "outputs": [],
   "source": [
    "mrca.1 = getMRCA(tree, tip = index.1)\n",
    "mrca.4 = getMRCA(tree, tip = index.4)\n",
    "mrca.14a = getMRCA(tree, tip = index.14a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bafd2fae-8406-46b6-a310-85f9fde5e469",
   "metadata": {},
   "source": [
    "## Identifying the exchanged genes per node"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80899e1b-7671-4865-96e7-bfcb959661f4",
   "metadata": {},
   "source": [
    "### Auxiliary functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9f0c92b-2e3e-4984-abac-a1c11374d499",
   "metadata": {},
   "source": [
    "The conversion between node IDs (number assigned internally by `ape`) and node labels (number assigned externally by `Count`) has been pushed in such a way by the conversion script for the inputs for `Count` (`convert_ASR_inputs.R`) that they can be easily interconverted by adding or subtracting the number of terminal leafs (i.e. the leafs that had a non-numerical label before input conversion)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea3cd6f8-bcaf-4eae-b284-fa33f60d2239",
   "metadata": {},
   "outputs": [],
   "source": [
    "nodeID_to_nodeLabel = function(nodeID, tree) {\n",
    "    return(nodeID - length(tree$tip.label))}\n",
    "nodeLabel_to_nodeID = function(nodeLabel, tree) {\n",
    "    return(nodeLabel + length(tree$tip.label))}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4183bc4a-0046-42b0-8066-cb44e6821adc",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Returns the node ID of a given label in a given phylogeny object.\n",
    "##\n",
    "## PARAMS\n",
    "## label      node label to get the ID of\n",
    "## tree       phylogeny in the form of an ape tree object\n",
    "##\n",
    "## OUTPUT\n",
    "## node ID of the tree node with the given label\n",
    "##\n",
    "findID = function(label, tree) {\n",
    "    lab = str_replace_all(label, ' ', '_')\n",
    "    ID.tip = which(tree$tip.label == lab)\n",
    "    if (length(ID.tip) == 0) {\n",
    "        ID.node = which(tree$node.label == lab)\n",
    "        ID = nodeLabel_to_nodeID(ID.node, tree)\n",
    "    }\n",
    "    else {\n",
    "        ID = ID.tip\n",
    "    }\n",
    "    return(ID)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25dfbfa9-3043-4e5b-ae40-1d40c6cf6aa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Sorts the gene exchange events into gene gain and loss events, neglecting the number of genes exchanged.\n",
    "## Also returns the full list of tree nodes and gene families.\n",
    "##\n",
    "## PARAMS\n",
    "## apa      the gene presence table of the ASR produced by Count (presences.tsv)\n",
    "## tree     the phylogeny in the form of an ape tree object\n",
    "##\n",
    "## OUTPUT\n",
    "## a list of the binary gained and lost sparse matrices, and the array of tree nodes and gene families\n",
    "##\n",
    "identify_type = function(apa, tree) {\n",
    "    families = apa$name\n",
    "    nodes = colnames(apa)[str_detect(colnames(apa), \"GCF [0-9]+\\\\.[0-9]|[1-9][0-9]+|[2-9]\")]\n",
    "\n",
    "    # Initialise two sparse matrices to represent whether a gene gain or loss, respectively, has taken place in a certain node\n",
    "    gained.c = Matrix(FALSE, length(families), length(nodes)+1, sparse = TRUE)\n",
    "    lost.c = gained.c\n",
    "\n",
    "    # We'll compare the presence of gene families of each node with the one of its ancestor node\n",
    "    for (node in nodes) {\n",
    "        # Find the node label of the ancestor node using the IDs of the current node and its ancestor in the phylogeny object\n",
    "        node.ID = findID(node, tree)\n",
    "        ancestor.ID = tree$edge[tree$edge[,2] == node.ID, 1]\n",
    "        ancestor = nodeID_to_nodeLabel(ancestor.ID, tree)\n",
    "\n",
    "        # Get the presence of gene families in both tree nodes from the ASR presence table from Count\n",
    "        pa.comp = subset(apa, select = c(ancestor, node))\n",
    "\n",
    "        # Gained genes are present in this node but not in the ancestor node\n",
    "        gained = pa.comp[,1] == 0 & pa.comp[,2] > 0\n",
    "        gained.c[,node.ID] = gained\n",
    "\n",
    "        # Lost genes are present in the ancestor node but not in this node\n",
    "        lost = pa.comp[,1] > 0 & pa.comp[,2] == 0\n",
    "        lost.c[,node.ID] = lost\n",
    "    }\n",
    "    \n",
    "    res = list('gained' = gained.c, 'lost' = lost.c, 'nodes' = nodes, 'families' = families)\n",
    "    return(res)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bfc6973-37cb-4731-a07e-0bfe6348cd7c",
   "metadata": {},
   "source": [
    "### Determining the exchange type of all gene families"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d4b29d0-5054-4d33-88a2-2a58684c8353",
   "metadata": {},
   "outputs": [],
   "source": [
    "identified.c = identify_type(apa, tree)\n",
    "gained.c = identified.c$gained\n",
    "lost.c = identified.c$lost\n",
    "nodes = identified.c$nodes\n",
    "families = identified.c$families"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95f4b1b1-61de-40b1-8ddb-989d233dfdec",
   "metadata": {},
   "outputs": [],
   "source": [
    "rm(apa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db736cb1-6314-4098-ab5a-eaccd73fe81a",
   "metadata": {},
   "outputs": [],
   "source": [
    "save(gained.c, lost.c, nodes, families, file = 'output/sorted_genes.RData')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49205eab-fa8f-4475-a2ff-c2f2415147ae",
   "metadata": {},
   "source": [
    "### Overview trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02406dce-8bce-4a13-bf76-e3fcd658f263",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Returns a ape tree object supplemented with the number of gene exchange events by tree node\n",
    "##\n",
    "## PARAMS\n",
    "## tree       phylogeny in the form an ape tree object\n",
    "## ggl_data   sparse binary matrix indicating whether a certain type of gene exchange event has taken place in a certain tree node;\n",
    "##            produced by identify_type()\n",
    "##\n",
    "## OUTPUT\n",
    "## an ape tree object joined with a gene exchange metadata column\n",
    "##\n",
    "extend_tree_with_ggl = function(tree, ggl_data) {\n",
    "    dt = data.frame(node = 1:(length(tree$node.label)+length(tree$tip.label)), trait = colSums(ggl_data))\n",
    "    gt = full_join(fortify(tree), dt, by = \"node\")\n",
    "    return(gt)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7a32f74-e206-49bf-9636-1d27bab14cb0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gt.gained = extend_tree_with_ggl(tree, gained.c)\n",
    "gt.gained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b428794f-12f0-42d3-9585-1542e79c0fea",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir.create('output/gained', recursive = TRUE)\n",
    "svg('output/gained/overview_tree.svg', width = 8, height = 12)\n",
    "ggtree(gt.gained, aes(color=.data$trait), size=0.5) +\n",
    "    labs(colour='Total genes gained') +\n",
    "    scale_color_gradientn(colours = magma(12), transform = \"log10\", na.value = \"black\")\n",
    "dev.off()\n",
    "display_svg(file = 'output/gained/overview_tree.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a8aa61e-df88-49b7-b043-a910915f7192",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gt.lost = extend_tree_with_ggl(tree, lost.c)\n",
    "gt.lost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86e72159-ef2e-4540-aad5-159491900529",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir.create('output/lost', recursive = TRUE)\n",
    "svg('output/lost/overview_tree.svg', width = 8, height = 12)\n",
    "ggtree(gt.lost, aes(color=.data$trait), size=0.5) +\n",
    "    labs(colour='Total genes lost') +\n",
    "    scale_colour_gradientn(colours = magma(12), transform = \"log10\", na.value = \"black\")\n",
    "dev.off()\n",
    "display_svg(file = 'output/lost/overview_tree.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a600459-6a4d-4c41-b9c9-e1c884b84c24",
   "metadata": {},
   "source": [
    "## Characterising the exchanged genes of some nodes in particular"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c57a86e-be0a-43d5-980d-3efba5a9af67",
   "metadata": {},
   "source": [
    "### Reading the annotation files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f75495b7-50bf-486f-a682-7fff19d42c0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cog = as.data.frame(fread(cog.file, header = FALSE, col.names = c('family', 'annotation', 'category')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e1d13f9-b016-4e38-b730-ee21595a5fbe",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bdd2265-d6d6-46ac-86c9-48a1d0cfced6",
   "metadata": {},
   "outputs": [],
   "source": [
    "kog = as.data.frame(fread(paste(root, kog.file, sep = \"/\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb9df1d5-1684-4643-ba3e-46b88599ca64",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7c23fda-9c2b-4ef2-9f8c-4bc613d3e97b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cog.cats = read.table(paste(root, cog.cats.file, sep = \"/\"), sep = \"\\t\", row.names = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54c60811-2a53-4947-be8b-78901ff8ed05",
   "metadata": {},
   "source": [
    "### Auxiliary functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef1c3978-5da7-465d-98b0-0847118327d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Gathers the gene family labels that were exchanged in a certain tree node\n",
    "##\n",
    "## PARAMS\n",
    "## var.c      sparse binary matrix produced by identify_type() indicating in which tree node a certain gene family was exchanged\n",
    "## families   array of gene family labels produced by identify_type()\n",
    "## mrca       tree node to get the exchange gene family labels for\n",
    "##\n",
    "## OUTPUT\n",
    "## a single-column dataframe listing the gene families that were exchanged in the given tree node\n",
    "##\n",
    "gather_genes = function(var.c, families, mrca) {\n",
    "    var.mrca = as.data.frame(families[var.c[,mrca]])\n",
    "    colnames(var.mrca) = c('family')\n",
    "    return(var.mrca)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6382da0-71fc-46cf-bf93-fa47a74a7b17",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Aggregates and counts the COG annotations of an array of gene families\n",
    "##\n",
    "## PARAMS\n",
    "## var        a single-column dataframe with a list of gene families; produced by gather_genes()\n",
    "## full_cog   the COG category annotations listed by gene family\n",
    "##\n",
    "## OUTPUT\n",
    "## a dataframe with relative COG category frequencies\n",
    "##\n",
    "aggregate_and_count_cog = function(var, full_cog) {\n",
    "    cog = left_join(var, full_cog, by = 'family') # left join to preserve the unknown gene families\n",
    "    cog.freq = cog %>% arrange(family) %>% count(category)\n",
    "\n",
    "    # Redistribute the plural annotations (e.g. 'BE')\n",
    "    for (c in cog.freq$category) {\n",
    "        # Case for known gene families\n",
    "        if (!is.na(c)) {\n",
    "            # Only do something for plural annotations\n",
    "            if (nchar(c) > 1) {\n",
    "                idx = which(cog.freq$category == c)\n",
    "                c.split = str_split(c, '')[[1]]\n",
    "                # Redistributing plural annotations to the separate single annotation categories, creating a new category if non-existent\n",
    "                for (cs in c.split) {\n",
    "                    where = which(cog.freq$category == cs)\n",
    "                    if (!length(where) == 0) {\n",
    "                        cog.freq[where,\"n\"] = cog.freq[where,\"n\"] + cog.freq[idx,\"n\"]\n",
    "                    }\n",
    "                    else {\n",
    "                        cog.freq[nrow(cog.freq)+1,] = list(cs, cog.freq[idx,\"n\"])\n",
    "                    }\n",
    "                }\n",
    "                cog.freq = cog.freq[-c(idx),]\n",
    "            }\n",
    "        }\n",
    "        # Unknown gene families end up with a NA label, so redistributing those to the unannotated category ('-')\n",
    "        else {\n",
    "            idx = which(is.na(cog.freq$category))\n",
    "            cog.freq[which(cog.freq$category == '-'),'n'] = cog.freq[which(cog.freq$category == '-'),'n'] + cog.freq[idx,'n']\n",
    "            cog.freq = cog.freq[-c(idx),]\n",
    "        }\n",
    "    }\n",
    "    # Convert to relative frequencies\n",
    "    cog.freq = arrange(cog.freq, category) %>% mutate(freq = n/sum(n)*100) %>% arrange(desc(n))\n",
    "    return(cog.freq)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec6017fc-d5d4-4830-9712-b16686273d67",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Aggregates and counts the KOG annotations at a given BRITE level of an array of gene families\n",
    "##\n",
    "## PARAMS\n",
    "## var        a single-column dataframe with a list of gene families; produced by gather_genes()\n",
    "## full_kog   the KOG category annotations listed by gene family and BRITE level\n",
    "## level      BRITE level to filter the KOG annotation dataframe for\n",
    "##\n",
    "## OUTPUT\n",
    "## a dataframe with relative KOG category frequencies of a given BRITE level\n",
    "##\n",
    "aggregate_and_count_kog = function(var, full_kog, level) {\n",
    "    kog = left_join(var, full_kog, by = join_by('family' == 'ID')) %>% arrange(family) # left join to preserve the unknown gene families\n",
    "    kog.freq = subset(kog, select = c('family', level)) # filter\n",
    "        %>% distinct() # deduplicate\n",
    "        %>% arrange(family) %>% count(.data[[level]]) # count\n",
    "        %>% mutate(freq = n/sum(n)*100) %>% arrange(desc(n)) # convert to relative frequencies\n",
    "    return(kog.freq)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44d34b27-c845-453b-9822-46afdf73c810",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Gathers and counts both COG and KOG annotations of both gene exchange events for a certain tree node\n",
    "##\n",
    "## PARAMS\n",
    "## gained.c      sparse binary matrix indicating whether a gene gain event has taken place in a tree node; produced by identify_type()\n",
    "## lost.c        sparse binary matrix indicating whether a gene loss event has taken place in a tree node; produced by identify_type()\n",
    "## families      array of gene family labels; produced by identify_type()\n",
    "## mrca          tree node for which the annotations of the exchanged genes need to be examined\n",
    "## cog           full COG annotation table for this genome set; expects the converted COG annotation list that was the input for Count\n",
    "## kog           full KOG annotation table for this genome set; expects the processed KOG-BRITE table produced by notebook 07b2\n",
    "## write_file    flag indicating whether the frequency tables should be written away (default = FALSE)\n",
    "## output        output directory in which the frequency tables will be saved as tsv files; ignored if write_file is FALSE\n",
    "## prefix        prefix for the filename of the frequency tables to distinguish different genome sets; ignored if write_file is FALSE\n",
    "##\n",
    "## OUTPUT\n",
    "## a nested list of both gained and lost gene families with their COG and KOG annotations at all levels\n",
    "##\n",
    "gather_and_count = function(gained.c, lost.c, families, mrca, cog, kog, write_file = FALSE, output = NULL, prefix = NULL) {\n",
    "    ## Gained genes\n",
    "    gained = gather_genes(gained.c, families, mrca)\n",
    "    writeLines(c(\"Number of genes gained\", nrow(gained)))\n",
    "    cog.gained.freq = aggregate_and_count_cog(gained, cog)\n",
    "    kog.gained.freq.A = aggregate_and_count_kog(gained, kog, 'A')\n",
    "    kog.gained.freq.B = aggregate_and_count_kog(gained, kog, 'B')\n",
    "    kog.gained.freq.C = aggregate_and_count_kog(gained, kog, 'C')\n",
    "    kog.gained.freq.D = aggregate_and_count_kog(gained, kog, 'D')\n",
    "\n",
    "    # Saving results\n",
    "    if write_file {\n",
    "        dir.create(paste(output, 'gained', sep = \"/\"), recursive = TRUE)\n",
    "        path = paste(output, 'gained', prefix, sep = \"/\")\n",
    "        write.table(cog.gained.freq, paste(path, 'cog', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.gained.freq.A, paste(path, 'kog.A', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.gained.freq.B, paste(path, 'kog.B', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.gained.freq.C, paste(path, 'kog.C', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.gained.freq.D, paste(path, 'kog.D', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        writeLines('Results for gained genes saved!')\n",
    "    }\n",
    "\n",
    "    ## Lost genes\n",
    "    lost = gather_genes(lost.c, families, mrca)\n",
    "    writeLines(c(\"Number of genes lost\", nrow(lost)))\n",
    "    cog.lost.freq = aggregate_and_count_cog(lost, cog)\n",
    "    kog.lost.freq.A = aggregate_and_count_kog(lost, kog, 'A')\n",
    "    kog.lost.freq.B = aggregate_and_count_kog(lost, kog, 'B')\n",
    "    kog.lost.freq.C = aggregate_and_count_kog(lost, kog, 'C')\n",
    "    kog.lost.freq.D = aggregate_and_count_kog(lost, kog, 'D')\n",
    "\n",
    "    # Saving results\n",
    "    if write_file {\n",
    "        dir.create(paste(output, 'lost', sep = \"/\"), recursive = TRUE)\n",
    "        path = paste(output, 'lost', prefix, sep = \"/\")\n",
    "        write.table(cog.lost.freq, paste(path, 'cog', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.lost.freq.A, paste(path, 'kog.A', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.lost.freq.B, paste(path, 'kog.B', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.lost.freq.C, paste(path, 'kog.C', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        write.table(kog.lost.freq.D, paste(path, 'kog.D', sep = \".\"), sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "        writeLines('Results for lost genes saved!')\n",
    "    }\n",
    "\n",
    "    res = list('gained' = \n",
    "               list('listed' = gained,\n",
    "                    'cog' = cog.gained.freq, \n",
    "                    'kog' = list('A' = kog.gained.freq.A, 'B' = kog.gained.freq.B, 'C' = kog.gained.freq.C, 'D' = kog.gained.freq.D)),\n",
    "               'lost' =\n",
    "               list('listed' = lost,\n",
    "                    'cog' = cog.lost.freq,\n",
    "                    'kog' = list('A' = kog.lost.freq.A, 'B' = kog.lost.freq.B, 'C' = kog.lost.freq.C, 'D' = kog.lost.freq.D))\n",
    "               )\n",
    "    return(res)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7012c0d-d4ee-41ca-b35f-e04d9a7d4892",
   "metadata": {},
   "source": [
    "### MRCA of cluster 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c461781-bad7-4b32-b98d-4f86c265273f",
   "metadata": {},
   "outputs": [],
   "source": [
    "counted = gather_and_count(gained.c, lost.c, families, mrca.1, cog, kog, TRUE, 'output', 'mrca.1')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8495111-6636-4561-9cc6-6b50c795fa84",
   "metadata": {},
   "source": [
    "#### Gained genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3a49be9-39ee-4c9f-967b-fe504243b722",
   "metadata": {},
   "outputs": [],
   "source": [
    "gained.1 = counted$gained$listed\n",
    "cog.1.gained.freq = counted$gained$cog\n",
    "kog.1.gained.freq.A = counted$gained$kog$A\n",
    "kog.1.gained.freq.B = counted$gained$kog$B\n",
    "kog.1.gained.freq.C = counted$gained$kog$C\n",
    "kog.1.gained.freq.D = counted$gained$kog$D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a558e998-41bf-47dc-9d35-3499c4720fa2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gained.1 %>% left_join(cog, by = \"family\") %>% arrange(annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aca19a4-4fe4-4512-b63d-c0975f267d83",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kog.1.gained.freq.C %>% arrange(desc(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc94e3cc-ba32-4d24-8286-576b0bab14ba",
   "metadata": {},
   "source": [
    "#### Lost genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c45ded70-5838-45b1-98b4-94063f958f50",
   "metadata": {},
   "outputs": [],
   "source": [
    "lost.1 = counted$lost$listed\n",
    "cog.1.lost.freq = counted$lost$cog\n",
    "kog.1.lost.freq.A = counted$lost$kog$A\n",
    "kog.1.lost.freq.B = counted$lost$kog$B\n",
    "kog.1.lost.freq.C = counted$lost$kog$C\n",
    "kog.1.lost.freq.D = counted$lost$kog$D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c86adde6-82d3-4840-aa8a-22bb13827b81",
   "metadata": {},
   "outputs": [],
   "source": [
    "lost.1 %>% left_join(cog, by = \"family\") %>% arrange(annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af3d5f9d-4b70-419a-a198-a65d3465862d",
   "metadata": {},
   "outputs": [],
   "source": [
    "kog.1.lost.freq.C %>% arrange(desc(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b3e28fe-a4bb-4008-9be6-f794690f4083",
   "metadata": {},
   "source": [
    "### MRCA of cluster 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95537adc-9138-4af2-aeb8-034b1ffacd86",
   "metadata": {},
   "outputs": [],
   "source": [
    "counted = gather_and_count(gained.c, lost.c, families, mrca.4, cog, kog, TRUE, 'output', 'mrca.4')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac4e0e9a-2438-480c-bc1d-9f4ebcaa00f0",
   "metadata": {},
   "source": [
    "#### Gained genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8e65280-530e-4dfe-adc7-b81b2835eede",
   "metadata": {},
   "outputs": [],
   "source": [
    "gained.4 = counted$gained$listed\n",
    "cog.4.gained.freq = counted$gained$cog\n",
    "kog.4.gained.freq.A = counted$gained$kog$A\n",
    "kog.4.gained.freq.B = counted$gained$kog$B\n",
    "kog.4.gained.freq.C = counted$gained$kog$C\n",
    "kog.4.gained.freq.D = counted$gained$kog$D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb882350-a349-4d0d-ae22-76c856c1444b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gained.4 %>% left_join(cog, by = \"family\") %>% arrange(annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41e9fe17-b7ab-46a2-a884-9b9260bef19a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kog.4.gained.freq.C %>% arrange(desc(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb0b73e4-61d8-4c17-86f0-e196ae3f0443",
   "metadata": {},
   "source": [
    "#### Lost genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70b1644d-c87b-4750-bb4a-10ad4d918889",
   "metadata": {},
   "outputs": [],
   "source": [
    "lost.4 = counted$lost$listed\n",
    "cog.4.lost.freq = counted$lost$cog\n",
    "kog.4.lost.freq.A = counted$lost$kog$A\n",
    "kog.4.lost.freq.B = counted$lost$kog$B\n",
    "kog.4.lost.freq.C = counted$lost$kog$C\n",
    "kog.4.lost.freq.D = counted$lost$kog$D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34b44cc5-d79a-4640-8bd6-220d13e49877",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lost.4 %>% left_join(cog, by = \"family\") %>% arrange(annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3849d5f4-0f65-4a9d-a38a-a750ae62ef8f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kog.4.lost.freq.C %>% arrange(desc(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5664daf4-807e-42a0-a5d0-953042bb47ad",
   "metadata": {},
   "source": [
    "### MRCA of cluster 14a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b68498a-ca58-4dfe-8d7a-f0826309cf65",
   "metadata": {},
   "outputs": [],
   "source": [
    "counted = gather_and_count(gained.c, lost.c, families, mrca.14a, cog, kog, TRUE, 'output', 'mrca.14a')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f9904fd-b9f3-4b7f-9bab-100cc12dfa6c",
   "metadata": {},
   "source": [
    "#### Gained genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88630e83-8c9c-4a0e-a85a-05eb7882d8dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "gained.14a = counted$gained$listed\n",
    "cog.14a.gained.freq = counted$gained$cog\n",
    "kog.14a.gained.freq.A = counted$gained$kog$A\n",
    "kog.14a.gained.freq.B = counted$gained$kog$B\n",
    "kog.14a.gained.freq.C = counted$gained$kog$C\n",
    "kog.14a.gained.freq.D = counted$gained$kog$D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "022fe176-503d-47a5-8397-778af52f9e8c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "gained.14a %>% left_join(cog, by = \"family\") %>% arrange(annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dda6cca8-dbc8-4441-bbc7-fb50cb93052a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kog.14a.gained.freq.C"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccc1efe1-4cc5-4e26-a8c2-ca489fa88766",
   "metadata": {},
   "source": [
    "#### Lost genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51020f6e-8531-445c-a011-dba94866da0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lost.14a = counted$lost$listed\n",
    "cog.14a.lost.freq = counted$lost$cog\n",
    "kog.14a.lost.freq.A = counted$lost$kog$A\n",
    "kog.14a.lost.freq.B = counted$lost$kog$B\n",
    "kog.14a.lost.freq.C = counted$lost$kog$C\n",
    "kog.14a.lost.freq.D = counted$lost$kog$D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "990a4432-d7da-4c1e-aea1-60c45f985b1c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lost.14a %>% left_join(cog, by = \"family\") %>% arrange(annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a661beb-2864-40c0-90a9-ab2c332d631b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kog.14a.lost.freq.C"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "435105b2-02ba-496b-81e1-a22a2e01361b",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccd1c50b-67dc-4cef-8579-92b59a82ebc0",
   "metadata": {},
   "source": [
    "#### Gained genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40090fa0-4ecf-4e2b-ba38-ee10e252de83",
   "metadata": {},
   "outputs": [],
   "source": [
    "gained.cog.toplot = bind_rows(list('1' = cog.1.gained.freq, '4' = cog.4.gained.freq, '14a' = cog.14a.gained.freq), .id = 'cluster')\n",
    "n.colors = length(unique(gained.cog.toplot$category))\n",
    "present.cogs = sort(intersect(unique(gained.cog.toplot$category), rownames(cog.cats)))\n",
    "\n",
    "svg('output/gained/COG_Freqs_Gained.svg', height = 4.5, width = 6)\n",
    "ggplot(gained.cog.toplot, aes(x = cluster, y = freq, fill = category)) +\n",
    "  geom_bar(stat = \"identity\", position = \"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4),\n",
    "                    labels = factor(paste0(\"(\", present.cogs, \") \", cog.cats[present.cogs,]))) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster\", y = \"% COGs\", fill = \"COG category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/gained/COG_Freqs_Gained.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d06ed5b-9c71-4a1b-9a3d-c757339bca2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "gained.kog.A.toplot = bind_rows(list('1' = kog.1.gained.freq.A, '4' = kog.4.gained.freq.A, '14a' = kog.14a.gained.freq.A), .id = 'cluster')\n",
    "n.colors = length(unique(gained.kog.A.toplot$A))\n",
    "\n",
    "svg('output/gained/KOG_A_Freqs_Gained.svg', height = 2, width = 4.5)\n",
    "ggplot(gained.kog.A.toplot, aes(x = cluster, y = freq, fill = A)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4)) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster\", y = \"% KOGs\", fill = \"BRITE-A category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/gained/KOG_A_Freqs_Gained.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53fbfc02-f83e-407c-8310-bbfd95c58051",
   "metadata": {},
   "outputs": [],
   "source": [
    "gained.kog.B.toplot = bind_rows(list('1' = kog.1.gained.freq.B, '4' = kog.4.gained.freq.B, '14a' = kog.14a.gained.freq.B), .id = 'cluster')\n",
    "n.colors = length(unique(gained.kog.B.toplot$B))\n",
    "\n",
    "svg('output/gained/KOG_B_Freqs_Gained.svg', height = 5, width = 5)\n",
    "ggplot(gained.kog.B.toplot, aes(x = cluster, y = freq, fill = B)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4)) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster\", y = \"% KOGs\", fill = \"BRITE-B category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/gained/KOG_B_Freqs_Gained.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be80237c-7b50-490e-b616-bb79be30b9f7",
   "metadata": {},
   "source": [
    "#### Lost genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "178aa3db-adfa-4fda-8266-171dc4f9d67c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lost.cog.toplot = bind_rows(list('1' = cog.1.lost.freq, '4' = cog.4.lost.freq, '14a' = cog.14a.lost.freq), .id = 'cluster')\n",
    "n.colors = length(unique(lost.cog.toplot$category))\n",
    "present.cogs = sort(intersect(unique(lost.cog.toplot$category), rownames(cog.cats)))\n",
    "\n",
    "svg('output/lost/COG_Freqs_Lost.svg', height = 4, width = 5.5)\n",
    "ggplot(lost.cog.toplot, aes(x = cluster, y = freq, fill = category)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4),\n",
    "                    labels = factor(paste0(\"(\", present.cogs, \") \", cog.cats[present.cogs,]))) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster\", y = \"% COGs\", fill = \"COG category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/lost/COG_Freqs_Lost.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdbd163f-98cf-447a-8ac4-a3e038c0d630",
   "metadata": {},
   "outputs": [],
   "source": [
    "lost.kog.A.toplot = bind_rows(list('1' = kog.1.lost.freq.A, '4' = kog.4.lost.freq.A, '14a' = kog.14a.lost.freq.A), .id = 'cluster')\n",
    "n.colors = length(unique(lost.kog.A.toplot$A))\n",
    "\n",
    "svg('output/lost/KOG_A_Freqs_Lost.svg', height = 2, width = 4)\n",
    "ggplot(lost.kog.A.toplot, aes(x = cluster, y = freq, fill = A)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4)) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster\", y = \"% KOGs\", fill = \"BRITE-A category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/lost/KOG_A_Freqs_Lost.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "991ec28a-8529-4b2c-a968-07ac9b3b071d",
   "metadata": {},
   "outputs": [],
   "source": [
    "lost.kog.B.toplot = bind_rows(list('1' = kog.1.lost.freq.B, '4' = kog.4.lost.freq.B, '14a' = kog.14a.lost.freq.B), .id = 'cluster')\n",
    "n.colors = length(unique(lost.kog.B.toplot$B))\n",
    "\n",
    "svg('output/lost/KOG_B_Freqs_Lost.svg', height = 5, width = 4.5)\n",
    "ggplot(lost.kog.B.toplot, aes(x = cluster, y = freq, fill = B)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4)) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster\", y = \"% KOGs\", fill = \"BRITE-B category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/lost/KOG_B_Freqs_Lost.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "909ef6a8-f5a1-4026-96f5-6c94ec78531e",
   "metadata": {},
   "source": [
    "### Overview plot of the number of exchanged genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5126fe96-30e7-4116-85b8-80876d32864c",
   "metadata": {},
   "outputs": [],
   "source": [
    "numbers.exchanged.toplot = bind_rows(list(\n",
    "    '1' = bind_cols(list('gained' = nrow(gained.1), 'lost' = nrow(lost.1))),\n",
    "    '4' = bind_cols(list('gained' = nrow(gained.4), 'lost' = nrow(lost.4))),\n",
    "    '14a' = bind_cols(list('gained' = nrow(gained.14a), 'lost' = nrow(lost.14a)))\n",
    "    ), .id = 'cluster') %>% pivot_longer(cols = c('gained', 'lost'), names_to = \"type\", values_to = \"counts\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9857965a-c2b9-429b-8acf-e7632c9a9a28",
   "metadata": {},
   "outputs": [],
   "source": [
    "svg('output/exchanged.svg', height = 2, width = 3)\n",
    "ggplot(numbers.exchanged.toplot, aes(x = cluster, y = counts, fill = cluster)) +\n",
    "    geom_bar(stat = \"identity\") +\n",
    "    facet_grid2(~type, scales = \"free\", independent = 'y') +\n",
    "    scale_x_discrete(expand = c(0,0)) +\n",
    "    theme(panel.spacing = unit(1, \"lines\"), legend.position = \"none\") +\n",
    "    labs(x = 'rRNA cluster MRCA', y = 'Number of genes')\n",
    "dev.off()\n",
    "display_svg(file = \"output/exchanged.svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8fac683-05aa-4dca-bdf1-9bd0d0796e74",
   "metadata": {},
   "source": [
    "## Characterising the genes in the LCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dde115e9-5f9b-46b9-8db9-77d367a29d24",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca.genes = as.data.frame(families[common])\n",
    "colnames(lca.genes) = c(\"family\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9fb71f3-b3ce-4ed0-8e77-168ae25c48e1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lca.genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d07cb98-abbe-49b5-891e-ad0b236e6c93",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lca.cog.freq = aggregate_and_count_cog(lca.genes, cog)\n",
    "lca.kog.freq.A = aggregate_and_count_kog(lca.genes, kog, 'A')\n",
    "lca.kog.freq.B = aggregate_and_count_kog(lca.genes, kog, 'B')\n",
    "lca.kog.freq.C = aggregate_and_count_kog(lca.genes, kog, 'C')\n",
    "lca.kog.freq.D = aggregate_and_count_kog(lca.genes, kog, 'D')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c59da3c-43c5-4d7d-8d02-7b149e030626",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lca.kog.freq.C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a053938-f28c-403a-965f-7cae625833f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir.create('output/common', recursive = TRUE)\n",
    "write.table(lca.cog.freq, 'output/common/lca.cog', sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "write.table(lca.kog.freq.A, 'output/common/lca.kog.A', sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "write.table(lca.kog.freq.B, 'output/common/lca.kog.B', sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "write.table(lca.kog.freq.C, 'output/common/lca.kog.C', sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)\n",
    "write.table(lca.kog.freq.D, 'output/common/lca.kog.D', sep = \"\\t\", quote = FALSE, col.names = TRUE, row.names = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2e95734-c00b-422d-8e90-f49fb8e33669",
   "metadata": {},
   "source": [
    "#### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dcc22c9-dde8-48f6-b725-5dec3a43859b",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca.cog.toplot = bind_rows(list('LCA' = lca.cog.freq), .id = 'cluster')\n",
    "n.colors = length(unique(lca.cog.toplot$category))\n",
    "present.cogs = sort(intersect(unique(lca.cog.toplot$category), rownames(cog.cats)))\n",
    "\n",
    "svg('output/common/COG_Freqs.svg', height = 4, width = 5)\n",
    "ggplot(lca.cog.toplot, aes(x = cluster, y = freq, fill = category)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4),\n",
    "                    labels = factor(paste0(\"(\", present.cogs, \") \", cog.cats[present.cogs,]))) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"\", y = \"% COGs\", fill = \"COG category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/common/COG_Freqs.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "252dc0a9-a241-46dc-845f-38b9882fed7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca.kog.A.toplot = bind_rows(list('LCA' = lca.kog.freq.A), .id = 'cluster')\n",
    "n.colors = length(unique(lca.kog.A.toplot$A))\n",
    "\n",
    "svg('output/common/KOG_A_Freqs.svg', height = 2, width = 3.5)\n",
    "ggplot(lca.kog.A.toplot, aes(x = cluster, y = freq, fill = A)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4)) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"\", y = \"% KOGs\", fill = \"BRITE-A category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/common/KOG_A_Freqs.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7696e30-ca04-4254-9537-959c08195597",
   "metadata": {},
   "outputs": [],
   "source": [
    "lca.kog.B.toplot = bind_rows(list('LCA' = lca.kog.freq.B), .id = 'cluster')\n",
    "n.colors = length(unique(lca.kog.B.toplot$B))\n",
    "\n",
    "svg('output/common/KOG_B_Freqs.svg', height = 5, width = 4)\n",
    "ggplot(lca.kog.B.toplot, aes(x = cluster, y = freq, fill = B)) +\n",
    "  geom_bar(stat = \"identity\", position=\"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4)) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"\", y = \"% KOGs\", fill = \"BRITE-B category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"))\n",
    "dev.off()\n",
    "display_svg(file = 'output/common/KOG_B_Freqs.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8db82449-aba5-4f97-b207-cacc865cb508",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## Overview figure for supplementary figures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9f7186b-07af-42a1-9554-243cfa51d267",
   "metadata": {},
   "source": [
    "#### COG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "908ef05b-c76e-49ea-9132-fba0c579247b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cog.toplot = bind_rows(list('Gained' = gained.cog.toplot, 'Lost' = lost.cog.toplot, 'LCA' = lca.cog.toplot), .id = 'type')\n",
    "cog.toplot$type = factor(cog.toplot$type, levels = c('Gained', 'Lost', 'LCA'))\n",
    "cog.toplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0537eccf-6542-4169-8a70-9d690681d378",
   "metadata": {},
   "outputs": [],
   "source": [
    "n.colors = length(unique(cog.toplot$category))\n",
    "present.cogs = sort(intersect(unique(cog.toplot$category), rownames(cog.cats)))\n",
    "\n",
    "svg('output/ASR_COGannotation_overview.svg', height = 4.5, width = 8.5)\n",
    "ggplot(cog.toplot, aes(x = cluster, y = freq, fill = category)) +\n",
    "  facet_grid2(~type, scales = 'free', space = 'free') +\n",
    "  geom_bar(stat = \"identity\", position = \"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4),\n",
    "                    labels = factor(paste0(\"(\", present.cogs, \") \", cog.cats[present.cogs,]))) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster MRCA\", y = \"% COGs\", fill = \"COG category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"), text = element_text(size = 13))\n",
    "dev.off()\n",
    "display_svg(file = 'output/ASR_COGannotation_overview.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e63213e-bf89-4020-a66e-8ee0f80a86e9",
   "metadata": {},
   "source": [
    "#### KOG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "685b138a-01ff-4d34-bd83-2e3c785f5f16",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "kog.B.toplot = bind_rows(list('Gained' = gained.kog.B.toplot, 'Lost' = lost.kog.B.toplot, 'LCA' = lca.kog.B.toplot), .id = 'type')\n",
    "kog.B.toplot$type = factor(kog.B.toplot$type, levels = c('Gained', 'Lost', 'LCA'))\n",
    "kog.B.toplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b02ef2de-cb9b-4c2d-ab97-9a48e5e6426e",
   "metadata": {},
   "outputs": [],
   "source": [
    "n.colors = length(unique(kog.B.toplot$B))\n",
    "\n",
    "svg('output/ASR_KOGBannotation_overview.svg', height = 5.25, width = 7.5)\n",
    "ggplot(kog.B.toplot, aes(x = cluster, y = freq, fill = B)) +\n",
    "  facet_grid2(~type, scales = 'free', space = 'free') +\n",
    "  geom_bar(stat = \"identity\", position = \"stack\") +\n",
    "  scale_fill_manual(values = getPalette(n.colors), \n",
    "                    guide = guide_legend(ncol = 1, keyheight = 0.8, keywidth = 0.4)) +\n",
    "  scale_y_continuous(expand = c(0,0)) +\n",
    "  scale_x_discrete(expand = c(0,0)) +\n",
    "  labs(x = \"rRNA cluster MRCA\", y = \"% KOGs\", fill = \"BRITE-B category\") +\n",
    "  theme(panel.spacing = unit(1, \"lines\"), text = element_text(size = 13))\n",
    "dev.off()\n",
    "display_svg(file = 'output/ASR_KOGBannotation_overview.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd6f7635-8030-4bea-ad83-dfddbfeb29ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "gained.kog.B.toplot\n",
    "lost.kog.B.toplot\n",
    "lca.kog.B.toplot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29b6070a-96cc-479a-8616-afe245b95fd5",
   "metadata": {},
   "source": [
    "## Saving and session info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3827644-27ff-4b92-aa90-fe79ee8b0af6",
   "metadata": {},
   "outputs": [],
   "source": [
    "save.image(file = 'output/env_output.RData')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9b67cddb-329e-4b76-ad65-45518a8002f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "R version 4.1.2 (2021-11-01)\n",
       "Platform: x86_64-pc-linux-gnu (64-bit)\n",
       "Running under: Linux Mint 21.2\n",
       "\n",
       "Matrix products: default\n",
       "BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3\n",
       "LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so\n",
       "\n",
       "locale:\n",
       " [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              \n",
       " [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    \n",
       " [5] LC_MONETARY=nl_BE.UTF-8    LC_MESSAGES=en_US.UTF-8   \n",
       " [7] LC_PAPER=nl_BE.UTF-8       LC_NAME=C                 \n",
       " [9] LC_ADDRESS=C               LC_TELEPHONE=C            \n",
       "[11] LC_MEASUREMENT=nl_BE.UTF-8 LC_IDENTIFICATION=C       \n",
       "\n",
       "attached base packages:\n",
       "[1] stats     graphics  grDevices utils     datasets  methods   base     \n",
       "\n",
       "other attached packages:\n",
       " [1] stringr_1.5.1      ggtree_3.11.1      ape_5.8            IRdisplay_1.1     \n",
       " [5] pals_1.8           RColorBrewer_1.1-3 ggh4x_0.2.8        ggplot2_3.5.1     \n",
       " [9] Matrix_1.6-5       tidyr_1.3.1        dplyr_1.1.4        data.table_1.15.4 \n",
       "\n",
       "loaded via a namespace (and not attached):\n",
       " [1] treeio_1.27.0.002  pbdZMQ_0.3-11      tidyselect_1.2.0   repr_1.1.7        \n",
       " [5] purrr_1.0.2        lattice_0.22-5     ggfun_0.1.4        colorspace_2.1-0  \n",
       " [9] vctrs_0.6.5        generics_0.1.3     htmltools_0.5.8.1  base64enc_0.1-3   \n",
       "[13] gridGraphics_0.5-1 utf8_1.2.4         rlang_1.1.3        pillar_1.9.0      \n",
       "[17] glue_1.7.0         withr_3.0.0        uuid_1.2-0         lifecycle_1.0.4   \n",
       "[21] munsell_0.5.0      gtable_0.3.4       mapproj_1.2.11     evaluate_0.23     \n",
       "[25] memoise_2.0.1      fastmap_1.1.1      parallel_4.1.2     fansi_1.0.6       \n",
       "[29] Rcpp_1.0.12        scales_1.3.0       cachem_1.0.8       IRkernel_1.3.2    \n",
       "[33] jsonlite_1.8.8     fs_1.6.3           digest_0.6.35      aplot_0.2.2       \n",
       "[37] stringi_1.8.3      grid_4.1.2         cli_3.6.2          tools_4.1.2       \n",
       "[41] yulab.utils_0.1.4  magrittr_2.0.3     maps_3.4.2         patchwork_1.2.0   \n",
       "[45] lazyeval_0.2.2     tibble_3.2.1       dichromat_2.0-0.1  crayon_1.5.2      \n",
       "[49] pkgconfig_2.0.3    tidytree_0.4.6     ggplotify_0.1.2    R6_2.5.1          \n",
       "[53] nlme_3.1-155       compiler_4.1.2    "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sessionInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1546c0bd-ee87-4fab-82c6-2c9ea544d885",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}