Day3PM_GORG_classifier.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f8450dc9-4557-4926-9187-11c0e8ed73b5",
   "metadata": {},
   "source": [
    "# GORG classifier\n",
    "\n",
    "We developed the GORG (Global Oceans Reference Genomes) classifier to recruit reads to a database made up of single cell genomes and assign taxonomy and function to reads. \n",
    "\n",
    "We have already gone ahead and recruited the same reads from metagenomes that we used this morning to the GORG database, [described here](https://www.sciencedirect.com/science/article/pii/S0092867419312735).\n",
    "\n",
    "The outputs from the GORG-classifier_setup notebook are in this folder storage/lesson_analyses/gorg-classifier\n",
    "\n",
    "Lets take a look in that folder now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ebff344-0e0a-406b-a139-a3c9f17b5502",
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls /mnt/storage/lesson_analyses/gorg-classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a17170c3-c997-4030-8199-b02baccb91a6",
   "metadata": {},
   "source": [
    "Let's first examine the summaries folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c19bbf1-9115-4c97-aadb-9c279f9fc53a",
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls /mnt/storage/lesson_analyses/gorg-classifier/summaries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "435cb902-fc57-491f-ad19-0e57cb804e4f",
   "metadata": {},
   "source": [
    "There is a summary for each metagenome run through the GORG-classifier.  These are short files meant to give a very quick assessment of how the GORG-classifier did in assigning reads to taxonomic and functional categories.\n",
    "\n",
    "The metagenome SRR6507279 is from the same location and timepoint the SAG plate we've been looking at in previous lessons, AG-910. Let's look at this one first.  \n",
    "\n",
    "Note that we are recruiting to a subset of 5,536 SAGs from the GORG-tropics dataset and not just the AG-910 plate that we have been using in previous lessons. The summaries file produces a nice short text file that is only 12 lines long so by increasing the number of lines that we use in the head command we can look at the entire file at once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45dba11b-24d5-4d03-bcff-b92a5c7857de",
   "metadata": {},
   "outputs": [],
   "source": [
    "! head /mnt/storage/lesson_analyses/gorg-classifier/summaries/SRR6507279_summary.txt -n 15"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a74bd670-18e1-4ce4-9960-ffb4c728c299",
   "metadata": {
    "tags": []
   },
   "source": [
    "Now lets look at the annotations file to see what we have in there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3450af3d-2c3f-410a-887c-6985ce1dc8e3",
   "metadata": {},
   "source": [
    "Wow! we are recruiting nearly 77% of all of the reads just to gene sequences and most of these reads (61% of the total) have functional annotations to Non-hypothetical genes. \n",
    "\n",
    "If you recall from this morning, we were only able to recruit a maximum of around 30% of metagenomic reads to SAGs from plate AG-910.  Expanding the reference database to GORG-Tropics, which uses a curated collection of 5,536 SAGs from the 12,715 cell collection, allowed us to recruit a much greater number of metagenomic reads overall.  \n",
    "\n",
    "The real meat of the results are in the **annotations** output directory.  Let's explore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c331af93-c97a-451a-99be-e4b5c1e35715",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "! ls /mnt/storage/lesson_analyses/gorg-classifier/annotations/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7917cbd1-01fb-42e1-b0c3-5143864be98c",
   "metadata": {},
   "source": [
    "We can look at the first 10 rows of these compressed files like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae5da55-a1c0-4603-b7e8-bc6fe1e78eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "!zcat /mnt/storage/lesson_analyses/gorg-classifier/annotations/SRR5720233_annotated.txt.gz | head -n10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ab65355-6217-489b-8cf7-810c7711050b",
   "metadata": {},
   "source": [
    "Alright there is a lot going on there and the formatting is a bit strange so lets go ahead and load all of our libraries and load this file into a dataframe to make it a bit easier to interpret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54be7172-6b74-402e-87f6-138308dccb37",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import os\n",
    "import sys\n",
    "import csv\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import glob\n",
    "import seaborn as sns\n",
    "import time\n",
    "\n",
    "# modify pandas variabls to only show 10 rows.\n",
    "pd.set_option(\"display.max_rows\", 10)\n",
    "\n",
    "# change default figure size \n",
    "plt.rcParams[\"figure.figsize\"] = (8,6)\n",
    "\n",
    "# create output directories\n",
    "outdir = \"./tbls/\"\n",
    "figdir = \"./figs/\"\n",
    "\n",
    "if not os.path.exists(outdir):\n",
    "    os.mkdir(outdir)\n",
    "if not os.path.exists(figdir):\n",
    "    os.mkdir(figdir)\n",
    "\n",
    "\n",
    "annotations_file = '/mnt/storage/lesson_analyses/gorg-classifier/annotations/SRR6507279_annotated.txt.gz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52717af1-e504-40e0-a4d2-27c12efcc80c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=pd.read_csv(annotations_file, compression='gzip', sep='\\t')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "168521c1-4e90-4399-a744-43b9a1825e83",
   "metadata": {},
   "source": [
    "This file contains a row for **every single read** within a metagenome.  You can imagine how these can get pretty massive, and manipulating them in a notebook could get unweildy.  For this lesson, the metagenomes we're using were subsampled to 1 million reads, so these tables are manageable.\n",
    "\n",
    "OK there is a lot going on here lets break this down a little bit by looking at each of the important columns again\n",
    " - status: Tells us if the read was Classified (C) or Unclassified (U)\n",
    " - sequence_id: read name\n",
    " - taxonomy_id: NCBI tax_ID number. Indicates the NCBI taxonomy\n",
    " - taxonomy_ids_lca: NCBI tax_ID number of the last common ancestor if there are multiple matches. This is used to determine the taxonomic assignment\n",
    " - taxonomy_lineage: NCBI formatted taxonomy represented by the tax_ID number.\n",
    " - The rest of the columns tell us information about the protein using the Prokka annotations\n",
    "\n",
    "\n",
    "### Important Note:\n",
    "The taxonomic id here uses the NCBI based taxonomy and not the GTDB based taxonomy that we discussed yesterday. As a result these taxonomies are not directly compatable so it is challenging to try and get them to work together. We are in the process of updating the GORG classifier database so that the GTDB taxonomy will be available as a database as well as some other modifications to make it more user friendly, flexible, and provide a more informative output. Be on the lookout for an update to this program sometime soon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53e2e9f5-b481-40ec-9d0c-a4b4e9ee615c",
   "metadata": {},
   "source": [
    "## Functional characterization of metagenomic reads.\n",
    "\n",
    "Let's look at the functional profiles of this metagenome by taxonomic group. I am going to group this dataframe by the ```taxonomic_lineage``` and ```prokka_product``` columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70e0d653-f555-4fde-a432-9dfc0b5d1fb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# I choose to group on the lineage and the prokka_product but you can pick any columns that you want to. The count column (sequence_id)\n",
    "# needs to be one of the columns that you are not grouping by but other than that it can be anything.\n",
    "grouped=data[data['status']=='C'].groupby(['taxonomic_lineage', 'prokka_product'], as_index= False)['sequence_id'].count()\n",
    "\n",
    "grouped=grouped.rename(columns={'sequence_id':'number_of_reads'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d7b4495-b8f6-44fc-a6b5-f9d59f6d0563",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets start by looking at the most abundant functions by taxnomic group that are in this dataset\n",
    "\n",
    "grouped.sort_values(['number_of_reads'], inplace=True, ascending=False)\n",
    "grouped.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cac1b11-b39c-4520-b4b9-e09680404e8c",
   "metadata": {},
   "source": [
    "Not surprisingly the most abundant proteins are hypothetical proteins since we grouped all hypothetical proteins together in our groupby command. Lets get rid of those proteins so that we can look at the proteins that have a functional annotation and plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c0c890a-e2fb-436f-b7db-b4d628d2ceab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select the 20 most abundant proteins that are not hypothetical proteins\n",
    "sub_group=grouped[grouped['prokka_product']!= 'hypothetical protein'][:20]\n",
    "\n",
    "# genrate the plot \n",
    "sns.barplot(data=sub_group, x='prokka_product', y='number_of_reads', hue='taxonomic_lineage', dodge=False)\n",
    "\n",
    "# create axis labels and rotate axis\n",
    "plt.xlabel('gene')\n",
    "plt.ylabel('Number of reads recruited')\n",
    "plt.title('Most abundant genes in SRR6507279')\n",
    "plt.xticks(rotation=90)\n",
    "\n",
    "# move legend outside of chart\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "plt.savefig(os.path.join(figdir,\"read_recruitment_to_most_abundant_functions.pdf\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "564f2503-7124-45e7-a28f-029d7c7b41a5",
   "metadata": {},
   "source": [
    "Most of the most abundant genes are from Pelagibacteraceae. Not too surprising, but there are a couple of genes from either Flavobacteriaceae or Gammaproteobacteria.\n",
    "\n",
    "Let's next look at the most abundant gene functions overall regardless of taxnomic assignment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd4e3ca1-5606-40a6-b03d-2dbf0de5896e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# group by just the functional annotation to look at functional abundance\n",
    "functional_df=data[data['status']=='C'].groupby('prokka_product', as_index= False)['sequence_id'].count()\n",
    "\n",
    "#rename the sequence_id and sort\n",
    "functional_df=functional_df.rename(columns={'sequence_id':'number_of_reads'})\n",
    "functional_df.sort_values(['number_of_reads'], inplace=True, ascending=False)\n",
    "\n",
    "# extract the first 20 rows for the plot\n",
    "sub_functional_df=functional_df[~functional_df['prokka_product'].isin(['hypothetical protein', 'putative protein'])][:20]\n",
    "\n",
    "# generate the plot\n",
    "sub_functional_df.plot(kind='bar', x='prokka_product', y='number_of_reads')\n",
    "\n",
    "# generate labels\n",
    "plt.xlabel('functional annotation')\n",
    "plt.ylabel('number of reads recruited')\n",
    "plt.title('number of reads from SRR6507279 recruited to functional annotation')\n",
    "\n",
    "# move legend\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "# save plot\n",
    "plt.savefig(os.path.join(figdir,\"most_abundant_functions_in_SRR6507279.pdf\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04b27ece-2d58-4d16-80f1-51922f9a02ac",
   "metadata": {},
   "source": [
    "Cool so it looks like DNA directed RNA polymerases are the most common genes in this metagenome. \n",
    "\n",
    "Note that the two most abundant genes are essentially the same thing but that one of them has an ' at the very end and the other one does not. As a result of this the computer is identifying them as two different genes. \n",
    "\n",
    "It looks like there most of the most common genes are involved in imformation processing (RNA transcriptin, DNA replication, tRNA ligases, DNA repair or DNA synthesis) in some way. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "335a86bd-44da-4cee-b48e-a29086ee3a8e",
   "metadata": {},
   "source": [
    "## Taxonomic characterization of metagenomic reads\n",
    "\n",
    "Next, I am curious what taxonomic lineages are represented within this metagenome, so I will group by taxonomic lineage, and count reads belonging to each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d80d2892-ab94-4207-beb4-b3e9d031b9d9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#group by just the taxonomic lineage to look at lineage abundance\n",
    "tax_df=data[data['status']=='C'].groupby('taxonomic_lineage', as_index= False)['sequence_id'].count()\n",
    "\n",
    "# rename the sequence id to reflect that it is the number of reads and sort high to low\n",
    "tax_df=tax_df.rename(columns={'sequence_id':'number_of_reads'})\n",
    "tax_df.sort_values(['number_of_reads'], inplace=True, ascending=False)\n",
    "\n",
    "# extract the first 20 rows for the plot\n",
    "sub_tax_df=tax_df[:20]\n",
    "# split the classification so that we can use taxonomic levels to color the barchart\n",
    "sub_tax_df[['domain', 'phyla', 'class', 'order', 'family', 'genus', 'species', '']]=sub_tax_df.taxonomic_lineage.str.split(';', expand=True)\n",
    "\n",
    "\n",
    "# generate the plot. I will plot the log10 transformation of the number of reads in order to better visualize these values\n",
    "ax = sns.barplot(data= sub_tax_df, x='taxonomic_lineage', y = 'number_of_reads', hue='order', dodge=False)\n",
    "\n",
    "# create axis labels \n",
    "ax.set(xlabel = 'taxonomic lineage',\n",
    "       ylabel = 'Number of reads recruited',\n",
    "       title = 'Number of reads from SRR6507279 recruited to genera', \n",
    "       yscale = 'log')\n",
    "\n",
    "# rotate x-axis labels\n",
    "plt.xticks(rotation=90)\n",
    "\n",
    "\n",
    "# move legend\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "plt.savefig(os.path.join(figdir,\"most_abundant_lineages_in_SRR6507279.pdf\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b682c399-d542-4c66-b64c-bbf9abdb65a4",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "## 30 minute coffee and coding break\n",
    "\n",
    "During this time get some coffee and snacks if you would like. \n",
    "\n",
    "Additionally everyone pick a gene of interest and a taxonomic lineage of interest and determine:\n",
    "\n",
    "For the gene: in what lineage is it most abundant?  \n",
    "\n",
    "For the lineage: what functional genes are most abundant?  \n",
    "\n",
    "I'll leave some space below in this notebook for you to work on this.\n",
    "<br>\n",
    "\n",
    "When we come back we will look into the distribution of different functional categories in different lineages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1e114fd-f052-459a-8e11-c7bbf628c8e5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ceacae7-16c1-436c-8226-9da6143b4675",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c86b178-6e9a-4076-89e5-3da9bb8f688f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab743ebe-8cdd-445f-965a-983ef17725ff",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8adce295-1379-455e-affa-21d0c0c6d6c4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13668eb6-88d5-4c7a-855c-e5225d864676",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "327a80e3-330a-4932-940e-d0d580c580d6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bc06320-85f0-4dd1-8eaf-3b7358dabf52",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c159a88f-aae6-4d24-84b6-a2f147329662",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5c17140c-1758-4c01-a310-319793780a22",
   "metadata": {},
   "source": [
    "## Examining the distribution of COGs within a metagenome\n",
    "\n",
    "The annotation output files provide more than just a taxonomic assignment and prokka function. If a gene is aligned to the database, GORG-classifier provides an array of different functional characterizations.  Those include:\n",
    "\n",
    "* prokka_gene\n",
    "* prokka_EC_number\n",
    "* prokka_product\n",
    "* swissprot_gene\n",
    "* swissprot_EC_number\n",
    "* swissprot_product\n",
    "* swissprot_eggNOG\n",
    "* swissprot_KO\n",
    "* swissprot_Pfam\n",
    "* swissprot_CAZy\n",
    "* swissprot_TIGRFAMs\n",
    "\n",
    "These may look familiar to you as they are all flavors of annotaton provided in the standard SCGC SAG deliverable that Ramunas and Greg reviewed on day 1.\n",
    "\n",
    "Next, I'm going to look into the funcitonal profiles of this metagenome using the COG IDs provided in the ```swissprot_eggNOG``` column.\n",
    "\n",
    "\n",
    "To do this, I've manipulated the file ```tbls/cog-20.def.csv``` available on the [COG website](https://ftp.ncbi.nih.gov/pub/COG/COG2020/data/).  My manipulated table: ```tbls/cog-20_unambiguous.csv``` has eliminated COGs that are assigned to multiple functional categories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46ae0b51-a2bc-4efa-bf84-10e701563a6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# I generated this file from the cog-20.def.csv file that is available for download on the COG website. \n",
    "# the major difference is that I eliminated COGs that are assigned to multiple functional categories\n",
    "COG_summary=pd.read_csv('tbls/cog-20_unambiguous.csv', index_col=0)\n",
    "# Only keep the columns that we are interested in.\n",
    "COG_summary=COG_summary[['COG_ID', 'Functional_category', 'COG_name', 'Functional category description']].copy()\n",
    "COG_summary.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26fd8b30-bc6c-4482-8592-fb8aecfb97e9",
   "metadata": {},
   "source": [
    "Next, I'm going to group my annotation table by ```taxonomic_lineage``` and ```swissprot_eggNOG``` columns, and count how many reads fell into each grouping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23a6f8ad-fc72-4282-96c9-711950377f1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "COG_df=data[data['status']=='C'].groupby(['taxonomic_lineage', 'swissprot_eggNOG'], as_index= False)['sequence_id'].count()\n",
    "\n",
    "COG_df=COG_df.rename(columns={'sequence_id':'number_of_reads'})\n",
    "COG_df.sort_values(['number_of_reads'], inplace=True, ascending=False)\n",
    "COG_df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3669e4d-d036-46ab-967d-4b237ca6b4cc",
   "metadata": {},
   "source": [
    "Above you can see that the ```swissprot_eggNOG``` column contains two codes, an EGGNOG id and a CCOG id, separated by a ','.  I'm going to split that column into three columns, as some have more than one COG listed. I'm going to only consider the first listed COG.\n",
    "\n",
    "I'm also going to correct the COG ID by removing the first 'C'.  Not sure why swissprot has added that additional 'C'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "675f10be-c125-450e-bcd0-1ce7b972d2d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# split the eggNOG/COG column so that we can use the COG. There are several genes that have multiple COGs assigned to them. \n",
    "# For the purposes of this exercise we are just going to use the first one\n",
    "COG_df[['eggNOG', 'COG_ID', 'addl_COGs']]=COG_df['swissprot_eggNOG'].str.split(',', 2, expand=True)\n",
    "\n",
    "#for some reason eggNOG adds an extra \"C\" to the beginning of their COG designation so we need to remove that\n",
    "COG_df['COG_ID']=COG_df['COG_ID'].str[1:]\n",
    "COG_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fcb4cae-913c-4f1e-b76b-58bbd7623bb7",
   "metadata": {},
   "source": [
    "Finallly, I'm going only keep the COG_ID column, and merge it with the COG_summary table we'd loaded at the beginning of this section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2d578b0-6e1e-4bd8-b5d0-53ca878c989b",
   "metadata": {},
   "outputs": [],
   "source": [
    "COG_df = COG_df[['COG_ID','taxonomic_lineage', 'number_of_reads']].merge(COG_summary, on='COG_ID', how='left')\n",
    "\n",
    "COG_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5433bc23-8d1c-4df8-a1b4-1e7e55eeb3af",
   "metadata": {},
   "source": [
    "Now the fun part... let's group the COGs by their taxonomic lineage and the COG Functional category.  This time we need to sum the number of reads instead of counting them, since they'd already been counted within the previous set of transofrmations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48bb10eb-aee8-481b-8c40-167c5f009245",
   "metadata": {},
   "outputs": [],
   "source": [
    "func_cat_df=COG_df.groupby(['taxonomic_lineage', 'Functional category description', 'Functional_category'], as_index= False)['number_of_reads'].sum()\n",
    "\n",
    "func_cat_df.sort_values(['number_of_reads'], ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8ff3e9d-ba1d-41bc-8ed9-c60cc82d58fd",
   "metadata": {},
   "source": [
    "Now lets use the taxonomic abundance dataframe that we generated previously to extract the 5 most abundant lineages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8ebe40a-ebc6-45dd-9912-b9c98d992dc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we are extracting the five first five lineages and we already sorted this table so we know that they are the lineages with the most reads\n",
    "common_lineage=sub_tax_df['taxonomic_lineage'][:5].to_list()\n",
    "# check if the taxonomic_lineage is one of the five we extracted previously and if it is copy to a new df\n",
    "common_func_cat_df=func_cat_df[func_cat_df['taxonomic_lineage'].isin(common_lineage)]\n",
    "\n",
    "common_func_cat_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00f8f785-4e35-401b-8fb8-f2f12ea2884f",
   "metadata": {},
   "source": [
    "And now to plot the abundance of different COG categories by the 5 most abundant genera"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "573f5f5c-0e9c-4dc3-a73b-05c45125c41c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set plot size\n",
    "plt.rcParams['figure.figsize'] = [15, 10]\n",
    "\n",
    "# and now to plot it\n",
    "ax = sns.barplot(data= common_func_cat_df, x='Functional category description', y='number_of_reads', hue='taxonomic_lineage')\n",
    "\n",
    "# generate labels\n",
    "ax.set(xlabel = 'COG functional category',\n",
    "       ylabel = 'number of reads recruited')\n",
    "\n",
    "\n",
    "plt.xticks(rotation=90)\n",
    "\n",
    "# move legend\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "# save plot\n",
    "plt.savefig(os.path.join(figdir,\"COG_functional_categories.pdf\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dbd906e-f1fe-4e6f-a1ba-755ff953ffc2",
   "metadata": {},
   "source": [
    "Alright this really does not tell us all that much since the Pelagibacteraceae have the most reads being recruited to every functional category, because they are just so abundant in the ocean. Instead of looking at the total number of reads per category lets calculate the percent of reads that are mapping to each category for each taxonomic lineage. This will give us some idea of differences in the genome structure and gene content of each lineage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e248868-c95d-474e-b7ed-405206be8421",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#use our old df to extract the total number of reads that were recruited to each tax lineage\n",
    "total_reads=sub_tax_df[['taxonomic_lineage', 'number_of_reads']][:5].copy()\n",
    "total_reads=total_reads.rename(columns={'number_of_reads':'total_number_of_reads'})\n",
    "\n",
    "#merge the total number of reads with our dataframe and use it to calculate the percent of reads in each functional category\n",
    "common_func_cat_df=common_func_cat_df.merge(total_reads, on='taxonomic_lineage')\n",
    "common_func_cat_df['percent_of_tax_reads']=common_func_cat_df['number_of_reads']/common_func_cat_df['total_number_of_reads']*100\n",
    "\n",
    "\n",
    "# create the barplot\n",
    "sns.barplot(data= common_func_cat_df, x='Functional category description', y='percent_of_tax_reads', hue='taxonomic_lineage')\n",
    "\n",
    "# generate labels\n",
    "plt.xlabel('COG functional category')\n",
    "plt.xticks(rotation=90)\n",
    "plt.ylabel('percent of reads recruited per lineage')\n",
    "plt.title('percent of reads per COG functional category per lineage')\n",
    "\n",
    "# move legend\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "# make it a bit bigger\n",
    "plt.rcParams['figure.figsize'] = [15, 10]\n",
    "\n",
    "# save plot\n",
    "plt.savefig(os.path.join(figdir,\"COG_functional_categories_percent.pdf\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f265382-c092-4069-a884-703e77871461",
   "metadata": {},
   "source": [
    "Now we can start to see some difference between the lienages. The Alphaproteobacteria lineages Pelagibacteraceae and Rhodospirillaceae seem to have a higher percentage of reads (fraction of genome) related to Amino acid transport and metabolism and Energy production and conversion than the other lineages, while the Gammaproteobacteria lineage have the highest fraction of reads recruited to Cell/wall and replication compared to other lineages"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08cebe35-544d-4d10-80e7-988689ceca0f",
   "metadata": {},
   "source": [
    "The cell below is entirely to demonstrate how to do this for multiple metagenomes, imagine that you have a time series and you want to compare abundance/functions across the time series. You can do this in a notebook and with pandas but it will quickly get out of hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89fda29a-5142-485f-b294-4bff3663e6e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# start timer to check how long this takes \n",
    "start = time.time()\n",
    "\n",
    "# create an empty dataframe that we will populate in this cell\n",
    "example_master_df=pd.DataFrame(columns=['taxonomic_lineage', 'prokka_product'])\n",
    "\n",
    "# define the files that I want to summarize I am limiting it to ones where the number ends with a 2 to reduce time\n",
    "# If I wanted to do it for all files it would be *.txt.gz\n",
    "infiles=glob.glob('/mnt/storage/lesson_analyses/gorg-classifier/annotations/*2_annotated.txt.gz')\n",
    "# loop through the files\n",
    "for file in infiles:\n",
    "    # extract the base name so that we can use it to rename the column with the metagenome name\n",
    "    name=os.path.basename(file)\n",
    "    name,ext=name.split('_')\n",
    "    # read the file\n",
    "    data=pd.read_csv(file, compression='gzip', sep='\\t')\n",
    "    # count the number of reads per taxonomic_lineage and product and rename the \n",
    "    counts=data[data['status']=='C'].groupby(['taxonomic_lineage', 'prokka_product'], as_index= False)['sequence_id'].count()\n",
    "    counts=counts.rename(columns={'sequence_id':name})\n",
    "    \n",
    "    #merge the new dataframe with the master dataframe\n",
    "    example_master_df=example_master_df.merge(counts, on=['taxonomic_lineage', 'prokka_product'], how='outer')\n",
    "    \n",
    "#timer end    \n",
    "end = time.time()\n",
    "print('it took', end - start, 'seconds to load 4 txt files, count the number of reads per function and lineage and \\nmerge the counts for each metagenome into a single dataframe')\n",
    "\n",
    "example_master_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "368dd135-c414-4379-86aa-f65bb94dc61a",
   "metadata": {},
   "outputs": [],
   "source": [
    "example_master_df.to_csv(os.path.join(outdir,'example_master_metagenome_read_recruitment_summary.csv'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08dca66f-135d-4689-be36-0225c8a20706",
   "metadata": {},
   "source": [
    "So it took ~21 seconds to generate this df for 4 metagenomes where we recruited 1 million reads per metagenome. While it was doable for this course you can easily see how this would rapidly get out of hand if you tried to increase either \n",
    "1. the number of reads (1 billion+ read metagenomes are not uncommon) processing 1 billion reads from each of these four metagenomes would take ~10 hours\n",
    "1. the number of metagenomes (Imagine you are comparing hundreds of metagenomes from different locations)\n",
    "\n",
    "or both. \n",
    "\n",
    "In order to get around such issues we generally will write a script that will take care of doing the above summation and then once we have that we can just load this final dataframe into our notebook. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "321950dd-1014-4689-8662-46484a5a4368",
   "metadata": {},
   "source": [
    "## Further exploration ideas\n",
    "1. What species has the most reads recruiting to Photosystem I? how about beta-xylanase? or some other gene you are intersted in?\n",
    "1. What percent of reads recruit to Archaea compared to Bacteria? How about Alphaproteobacteria?\n",
    "1. As you saw this morning summer metagenomes had higher recruitment to AG-910 than winter metagenomes. We started with a summer sample. Get a winter metagenome and compare  taxonomic and functional differences between the two.\n",
    "     1. A winter metagenome SRR5720233\n",
    "         - Hint think about how you are going to ensure that you can differentiate between the two metagenomes\n",
    "     1. Metadata about other metagenomes can be found in this file (/mnt/storage/lessons/Day3PM_GORG-classifier/bats_metagenomes_of_interest.csv)\n",
    "         \n",
    "1. The first number of an EC number indicates the class of an enzyme. Use the EC number to look at the prevelance of different enzyme classes for a single species, between species, or between metagenomes\n",
    "     - Hint the dictionary below contains the enzyme classes associated with the first digit of the EC number\n",
    "```\n",
    "            EC_nums={1: 'Oxidoreductase', 2 : \"Transferases\", 3 : \"Hydrolases\", 4 : \"Lysases\", 5 : \"Isomerases\", 6 : \"Ligases\", 7 : \"Translocases\"}\n",
    "```\n",
    "\n",
    "There are also lots of other columns that you can use to investigate these results. Consider using the KO numbers instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7352cde6-435c-4dd4-899d-153239d7dc25",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "biopy",
   "language": "python",
   "name": "biopy"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}