2.3.StatsCircle.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stats Circle\n",
    "\n",
    "This notebook replicates the stats circle plots provided by BGI in their HTML reports using the outputs from STARsolo.\n",
    "\n",
    "This is not a perfect replication as the analysis pipelines are different, but it is a very close approximation.\n",
    "\n",
    "Notably, we do not have a count for `non-relevant` HDMIs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from math import pi\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import pandas as pd\n",
    "import scipy.io as sio\n",
    "\n",
    "from matplotlib.patches import Patch\n",
    "from matplotlib.lines import Line2D\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.facecolor'] = 'white'\n",
    "plt.rcParams['axes.facecolor'] = 'white'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples = [] # A list of the samples to process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "group_colours = [\n",
    "    [(255, 97, 79)],\n",
    "    [(255, 189, 132), (255, 206, 162)],\n",
    "    [(162, 210, 135)],\n",
    "    [(152, 193, 221), (194, 226, 249), (215, 237, 251)],\n",
    "    [(148, 158, 206), (177, 184, 215)],\n",
    "    [(247, 188, 255), (237, 210, 240)],\n",
    "]\n",
    "\n",
    "group_colours = [[tuple([c / 255 for c in rgb]) for rgb in group] for group in group_colours]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 1, figsize=(10, 10), subplot_kw=dict(projection=\"polar\"))\n",
    "for n, sample in enumerate(samples):\n",
    "    # ax = axs[n // 3, n % 3]\n",
    "    ax = axs\n",
    "    sample_folder = (\n",
    "        f\"output/{sample}/\"\n",
    "    )\n",
    "    main_log = {}\n",
    "    with open(f\"{sample_folder}/{sample}_Log.final.out\") as f:\n",
    "        for line in f:\n",
    "            if \"|\" in line:\n",
    "                line = line.strip().split(\"|\")\n",
    "                main_log[line[0].strip()] = line[1].strip()\n",
    "    bc_stats = {}\n",
    "    with open(f\"{sample_folder}/{sample}_Solo.out/Barcodes.stats\") as f:\n",
    "        for line in f:\n",
    "            line = line.strip().split()\n",
    "            bc_stats[line[0]] = int(line[1])\n",
    "    feature_stats = {}\n",
    "    with open(f\"{sample_folder}/{sample}_Solo.out/GeneFull/Features.stats\") as f:\n",
    "        for line in f:\n",
    "            line = line.strip().split()\n",
    "            feature_stats[line[0]] = int(line[1])\n",
    "    summary = pd.read_csv(f\"{sample_folder}/{sample}_Solo.out/GeneFull/Summary.csv\", sep=\",\", index_col=0, header=None)\n",
    "    mtx = sio.mmread(f\"{sample_folder}/{sample}_Solo.out/GeneFull/raw/matrix.mtx\")\n",
    "    dedup_umis = mtx.sum()\n",
    "    data = {\n",
    "        \"total_reads\": int(summary.loc[\"Number of Reads\"].values[0]),\n",
    "        \"valid_hdmis\": bc_stats[\"yesWLmatchExact\"] + bc_stats[\"yesOneWLmatchWithMM\"],\n",
    "        \"invalid_hdmis\": bc_stats[\"noNoWLmatch\"] + bc_stats[\"noUMIhomopolymer\"] + bc_stats[\"noTooManyWLmatches\"],\n",
    "        \"discarded_hdmis\": 0,\n",
    "        \"not_relevant\": 0,\n",
    "        \"mapped\": feature_stats[\"yesWLmatch\"] + feature_stats[\"noNoFeature\"],\n",
    "        \"unmapped\": feature_stats[\"noUnmapped\"],\n",
    "        \"multimap\": feature_stats[\"MultiFeature\"],\n",
    "        \"annotated\": feature_stats[\"yessubWLmatch_UniqueFeature\"],\n",
    "        \"unannotated\": feature_stats[\"noNoFeature\"],\n",
    "        \"dedup_umis\": dedup_umis,\n",
    "        \"dup_umis\": feature_stats[\"yessubWLmatch_UniqueFeature\"] - dedup_umis,\n",
    "    }\n",
    "    data[\"relevant\"] = data[\"valid_hdmis\"]\n",
    "    data_norm = {k: v / data[\"total_reads\"] * 100 for k, v in data.items()}\n",
    "    groups = [\n",
    "        [\"total_reads\"],\n",
    "        [\"valid_hdmis\", \"invalid_hdmis\"],\n",
    "        [\"relevant\"],\n",
    "        [\"mapped\", \"unmapped\", \"multimap\"],\n",
    "        [\"annotated\", \"unannotated\"],\n",
    "        [\"dedup_umis\", \"dup_umis\"],\n",
    "    ]\n",
    "    group_heights = [\n",
    "        -0.2,\n",
    "        1,\n",
    "        2.2,\n",
    "        3.4,\n",
    "        4.6,\n",
    "        5.8,\n",
    "    ]\n",
    "\n",
    "    for group, g_colours, height in zip(groups, group_colours, group_heights):\n",
    "        group_data = [data_norm[k] for k in group]\n",
    "        group_labels = [k for k in group]\n",
    "        startangle = 0\n",
    "        xs = [(i * pi * 2) / 100 for i in group_data]\n",
    "        for n, (dp, c) in enumerate(zip(xs, g_colours)):\n",
    "            left = (startangle * pi * 2) / 360  # this is to control where the bar starts\n",
    "            if n == 0:\n",
    "                startangle += dp * 360 / (2 * pi)\n",
    "            else:\n",
    "                startangle += dp * 360 / (2 * pi)\n",
    "            dp = dp - ((0.5 * pi * 2) / 100)\n",
    "            ax.barh(height, dp, left=left, height=1, color=c)\n",
    "\n",
    "\n",
    "    ax.spines.clear()\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(f'{sample}\\nSaturation {float(summary.loc[\"Sequencing Saturation\"]) * 100:.4f}%')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.9.5 - BiomedBundle",
   "language": "Python",
   "name": "biomedbundle"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}