context for analysis notebooks

07-analysis/README.md

+# Analysis
+
+The code in this folder concerns downstream analysis after genome assembly and annotation.
+
+### _Pycnogonum litorale_ in the chelicerate/arthropod context
+
+To gain a first understanding of the _P. litorale_ genome we compared it to publicly available
+genome assemblies for arthropods in general and chelicerates in particular using common metrics.
+This notebook produced most of the material for figures 2 and 3 in the manuscript.
+
+- [The genome assembly of _P. litorale_ in the broader arthropod/chelicerate context](genomic_context.ipynb)
+
+### Hox genes
+
+We used manually curated CDS's for _P. litorale_ Hox genes to locate them in the genome. The
+putative CDS's overlap with predicted gene models and form a single cluster on pseudochromosome 56.
+All Hox genes are present, except for _Hox9/abdominalA_, for which contradicting data exists for
+pycnogonids. This notebook produced material for figure 4 and suppl. figure 3 of the manuscript.
+
+- [Exploring Hox genes](hoxfinder.ipynb)
+
+### EggNOG-mapper output
+
+A brief look at the functional annotation provided by EggNOG-mapper for the predicted pycnogonid
+peptide sequences, broken down by prediction quality (also see how the annotation [was
+produced](../06-annotation/README.md)). We find that annotation quality correlates with annotation
+confidence (Iso-seq > BRAKER round 1 > BRAKER round 2 > de-novo).
+
+- [functional annotation overview](emapper_output.ipynb)
\ No newline at end of file

07-analysis/emapper_output.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Parsing functional annotation\n",
+    "\n",
+    "We used EggNOG-mapper to annotate the functional annotation of the proteins. The output of\n",
+    "EggNOG-mapper is a tab-separated file that contains information about the seed ortholog of each\n",
+    "query, as well as information about the corresponding orthogroup and its phylogenetic context. Here\n",
+    "we will parse this output to reduce it to one description per gene, and see how well the different\n",
+    "annotation sources are represented."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pandas as pd\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First read the data. Remember that emapper output has four comment lines in the header and another\n",
+    "three in the footer:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "emapper_loc = \"/Volumes/scratch/pycnogonum/genome/draft/annot_merge/out.emapper.annotations\"\n",
+    "\n",
+    "emapper = pd.read_csv(emapper_loc, sep=\"\\t\", skiprows=4, skipfooter=3, engine=\"python\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We need to extract the gene ID from the protein ID. This is mostly dot-separated, but PacBio IsoSeq \n",
+    "genes have more dots in their names, so we need to parse those ddifferently."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def parse_gene_id(x):\n",
+    "    if \"PB\" in x:\n",
+    "        parts = x.split(\".\")\n",
+    "        return \".\".join(parts[:2])\n",
+    "    elif x.startswith(\"r2\") or x.startswith(\"g\") or x.startswith(\"at\"):\n",
+    "        return x.split(\".\")[0]\n",
+    "    else:\n",
+    "        return ValueError(\"Unknown gene ID format\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "apply the parsing on each row..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "emapper[\"gene\"] = emapper[\"#query\"].apply(parse_gene_id)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "...and now group by gene ID and keep the best-scoring isoform/CDS for each gene."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "summary = emapper.groupby(\"gene\").apply(lambda x: x.sort_values(\"score\", ascending=False).head(1), include_groups=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we need to extract the annotation type from the gene IDs, something we can easily do by the \n",
+    "start of the gene ID.\n",
+    "\n",
+    "* P: PacBio IsoSeq (round 1)\n",
+    "* g: BRAKER3 (round 1)\n",
+    "* r: BRAKER3 (round 2)\n",
+    "* a: assembled transcriptomes (round 3)\n",
+    "\n",
+    "We will keep track of these in a dictionary to make calculations easier. First, let's find out how\n",
+    "many genes we have from each source.\n",
+    "\n",
+    "```bash\n",
+    "$ grep \"      gene    \" merged_sorted.gff3 | grep PB | wc -l\n",
+    "8904\n",
+    "$ grep \"      gene    \" merged_sorted.gff3 | grep r2 | wc -l\n",
+    "2223\n",
+    "$ grep \"      gene    \" merged_sorted.gff3 | grep at | wc -l\n",
+    "774\n",
+    "$ grep \"      gene    \" merged_sorted.gff3 | grep \"=g\" | wc -l\n",
+    "3596\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "translate = {\n",
+    "    \"P\": \"PacBio\",\n",
+    "    \"g\": \"BRAKER3_r1\",\n",
+    "    \"r\": \"BRAKER3_r2\",\n",
+    "    \"a\": \"Trinity\",\n",
+    "}\n",
+    "\n",
+    "totals = {\n",
+    "    \"PacBio\": 8904,\n",
+    "    \"BRAKER3_r1\": 3596,\n",
+    "    \"BRAKER3_r2\": 2223,\n",
+    "    \"Trinity\": 774\n",
+    "}\n",
+    "\n",
+    "rounds = {\n",
+    "    \"PacBio\": 1,\n",
+    "    \"BRAKER3_r1\": 1,\n",
+    "    \"BRAKER3_r2\": 2,\n",
+    "    \"Trinity\": 3\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "annotated = pd.DataFrame(summary[\"#query\"].str[:1].value_counts())\n",
+    "annotated.index = [translate[x] for x in annotated.index]\n",
+    "annotated[\"totals\"] = [totals[x] for x in annotated.index]\n",
+    "annotated[\"round\"] = [rounds[x] for x in annotated.index]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can now see what percentage of gene models find support for each annotation source."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "PacBio        73.69\n",
+       "BRAKER3_r1    57.40\n",
+       "BRAKER3_r2    51.73\n",
+       "Trinity       31.27\n",
+       "dtype: float64"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(annotated[\"count\"] / annotated[\"totals\"] * 100).round(2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As expected, the more trustworthy annotations find more orthologs. The same pattern holds if we\n",
+    "aggregate the round 1 sources:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "round\n",
+       "1    69.00\n",
+       "2    51.73\n",
+       "3    31.27\n",
+       "dtype: float64"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "tmp = annotated.groupby(\"round\").sum()\n",
+    "(tmp[\"count\"] / tmp[\"totals\"] * 100).round(2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Write out the reduced annotation file for downstream analysis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "summary_loc = \"/Volumes/scratch/pycnogonum/genome/draft/annot_merge/out.emapper.best.annotations\"\n",
+    "summary.to_csv(summary_loc, sep=\"\\t\", index=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "np.int64(10017)"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "annotated[\"count\"].sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "np.int64(15497)"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "annotated[\"totals\"].sum()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "genbank",
+   "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.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}