{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c65b4b47-c99c-4673-971d-c198e15b5a5d",
   "metadata": {
    "tags": []
   },
   "source": [
    "## __Get some basic mapping stats with samtools flagstat__\n",
    "\n",
    "As we will mainly launch python code to analyse mapping results, the kernel of this jupyter book is Python3. So, to run linux commands, we have to add :\n",
    "* `%` to execute a single linux command \n",
    "* `%%bash`  to execute a serie of linux commands\n",
    " \n",
    "### Merge individual flagstat files into an unique file `all_stat.csv` with python code\n",
    "\n",
    "* Add the path of the directory that contains flagstat files  in the following cell (line 8)\n",
    "```flagstat_dir = PATH_2_FLAGSTAT_DIRECTORY```\n",
    "* Then, execute the two following cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3365afa-1224-4a91-82f0-605406771dda",
   "metadata": {},
   "outputs": [],
   "source": [
    "# IMPORT PYTHON PACKAGE USED BY THE CODE\n",
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "# VARIABLE INITIALIZATION\n",
    "\n",
    "## NAME OF THE DIRECTORY THAT CONTAINS FLAGSTAT FILES\n",
    "flagstat_dir =  PATH_2_FLAGSTAT_DIRECTORY\n",
    "\n",
    "## NAME OF THE FILE THAT WILL CONTAIN ALL THE FLAGSTAT RESULTATS\n",
    "stat_file = f\"{flagstat_dir}/all_stat.csv\"\n",
    "\n",
    "# PRINT THE CONTENT OF 2 PREVIOUS VARIABLES INITIALIZED\n",
    "print(\"DIRECTORY : \",flagstat_dir)\n",
    "print(\"FINAL STAT FILE : \",stat_file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8da626d7-fc37-409a-9308-9d81baf089d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# OPEN THE FINAL FILE IN WHICH WE PRINT SOME STATS EXTRACTED \n",
    "# FROM EACH INDIVIDUAL FILE GENERATED BY SAMTOOLS FLAGSTAT\n",
    "with open(stat_file, 'w') as stat: \n",
    "    # WRITE A HEADER LINE IN OUR STAT FILE\n",
    "    header_line = \"sample,mapped,paired,unmapped\"\n",
    "    stat.write(header_line)\n",
    "    \n",
    "    # READING EACH FILE OF THE FLGSTAT DIRECTORY\n",
    "    for file in os.listdir(flagstat_dir):\n",
    "        \n",
    "        filen = flagstat_dir + \"/\" + file\n",
    "        \n",
    "        # If the word \"flagstat\" is in name of file\n",
    "        if \"flagstat\" in file:\n",
    "        \n",
    "            # Extract sample name and save into a new variable newLine \n",
    "            new_line = f\"\\n{file.split('.')[0]},\"\n",
    "            \n",
    "            # OPEN AND READS FLAGSTAT FILE\n",
    "            with open(filen, \"r\") as flagstat:\n",
    "            \n",
    "                # read file line by line\n",
    "                for line in flagstat:\n",
    "                    # remove the line skipper at the endo of the line\n",
    "                    line = line.rstrip()              \n",
    "                    # Keep only line mapped, paired or singleton word\n",
    "                    if ('mapped (' in line and not 'primary' in line) or 'paired (' in line or 'singleton' in line:\n",
    "                        # get percentage value and save it into the varaible called perc\n",
    "                        perc = f\"{line.split('(')[1].split('%')[0]}\"\n",
    "                        new_line += f\"{perc},\"\n",
    "                # WRITE THE LINE ONCE THE FLAGSTAT FILE COMPLETELY READ\n",
    "                stat.write(new_line.strip(\",\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eccefc94-20ef-4508-9ccf-d080f5f14288",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Display the content of the final stat file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f4a2294-5f5c-4ffb-9c0d-0fbea8dfc7e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%ls -lt /home/jovyan/SV_DATA/RESULTS/2_bam_flagstat/ | head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aff90309-8f07-49b2-88ce-8fe0c47fde20",
   "metadata": {},
   "outputs": [],
   "source": [
    "%cat $stat_file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a512002-2fe5-4172-9bc1-cab76a00c76d",
   "metadata": {},
   "source": [
    "#### __Plot mapping rate per sample__  \n",
    "\n",
    "###### Load csv file into a panda dataframe\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b87b0b6-e76c-4950-9b78-3ecd418f6704",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_bam_stat = pd.read_csv(stat_file, index_col=False, sep=\",\")\n",
    "df_bam_stat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb1b6397-f6d2-4d27-9f85-1ee363a4b499",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Basic stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d111bd5-3966-4d48-8a49-8131654ca1f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# only print the values of the \"mapped\" column\n",
    "print(df_bam_stat['mapped'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da8d2754-548d-494e-931e-d1e87d6334ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display only the minimum, maximum and average rate of mapped reads\n",
    "minimum = df_bam_stat[\"mapped\"].min()\n",
    "maximun = df_bam_stat[\"mapped\"].max()\n",
    "mean_flag = df_bam_stat[\"mapped\"].mean()\n",
    "\n",
    "print(\"\\n######## BASIC STATS\\n MAPPED\")       \n",
    "print(f\"\\t%min : {minimum}\\t %max : {maximun}\\t %mean : {mean_flag}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34836741-81ac-4b8f-8067-f0e6192eb3c2",
   "metadata": {},
   "source": [
    "#### Sort by sample name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87b44f61-bf83-4a7f-8938-7de01f5b6c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_bam_stat_sorted=df_bam_stat.sort_values(by=['sample'])\n",
    "df_bam_stat_sorted"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a3b978c-fa7e-45f2-91c5-003d15b722c3",
   "metadata": {},
   "source": [
    "#### Mapping Rate Plots with Python (Seaborn)\n",
    "\n",
    "* Plot the mapping rate (Y-axis) for each sample (X-axis) using the `mapped`and the `sample` columns of the  `df_bam_sorted`dataframe "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "213d0cd9-2c76-40dd-8843-486b99af7c9e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot with seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "plt.figure(figsize = (15,8))\n",
    "sns.scatterplot(x=\"sample\",y=\"paired\", data=df_bam_stat_sorted)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e89f1ba9-105e-433f-b34c-b9927e433347",
   "metadata": {},
   "source": [
    "* Plot the values of the three columns (mapped, paired, unmapped) for each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90db87ba-e830-4c7c-bd3c-5827cb43d637",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot with seaborn\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "plt.figure(figsize = (15,8))\n",
    "ax=sns.scatterplot(x=\"sample\",y=\"value\", hue='variable', data=pd.melt(df_bam_stat_sorted, 'sample'))\n",
    "ax.set_title(\"ADD_TITLE \")\n",
    "ax.set_xlabel(\"ADD_X_LABEL\")\n",
    "ax.set_ylabel(\"ADD_Y_LABEL\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa486934-f381-49a7-9e82-18ff3011ee5c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}