project/misc_illustrations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figures for Illustration of concepts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "FIGUREFOLDER = Path('figures')\n",
    "FIGUREFOLDER.mkdir(exist_ok=True, parents=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'xtick.labelsize': 'xx-large',\n",
    "                     'ytick.labelsize': 'xx-large',\n",
    "                     'axes.titlesize': 'xx-large',\n",
    "                     'axes.labelsize': 'xx-large',\n",
    "                     })\n",
    "# {k:v for k,v in plt.rcParams.items() if 'tick' in k and 'size' in k}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imputation by random draw from normal distribution\n",
    "\n",
    "- currently commenly used approach at NNF CPR for downstream experimentation task\n",
    "- see also illustrations in [Lazar 2016, Figure 1](https://pubs.acs.org/doi/10.1021/acs.jproteome.5b00981#fig1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 25.0\n",
    "stddev = 1.0\n",
    "\n",
    "x = np.linspace(mu - 3, mu + 3, num=101)\n",
    "\n",
    "y_normal = scipy.stats.norm.pdf(x, loc=mu, scale=stddev)\n",
    "\n",
    "mu_shifted = mu - (1.8 * stddev)\n",
    "stddev_shifted = 0.3 * stddev\n",
    "print(f\"Downshifted: {mu_shifted = }, {stddev_shifted = }\")\n",
    "y_impute = scipy.stats.norm.pdf(x, loc=mu - (1.8 * stddev), scale=0.3 * stddev)\n",
    "\n",
    "colors = plt.cm.viridis([0.25, 0.75])\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(5, 4))\n",
    "\n",
    "for y, c in zip([y_normal, y_impute], colors):\n",
    "    ax.plot(x, y, color=c,)\n",
    "    ax.fill_between(x, y, color=c)\n",
    "ax.set_xlabel('log2 intensity')\n",
    "ax.set_ylabel('density')\n",
    "ax.set_label(\"test\")\n",
    "ax.legend([\"original\", \"down shifted\"])\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_normal_imputation')\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_normal_imputation.pdf')\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_normal_imputation_highres', dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3646ec23",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'xtick.labelsize': 'large',\n",
    "                     'ytick.labelsize': 'large',\n",
    "                     'axes.titlesize': 'large',\n",
    "                     'axes.labelsize': 'large',\n",
    "                     })\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3, 2))\n",
    "\n",
    "for y, c in zip([y_normal], colors):\n",
    "    ax.plot(x, y, color=c,)\n",
    "    # ax.fill_between(x, y, color=c)\n",
    "    ax.set_xlabel('log2 intensity')\n",
    "    ax.set_ylabel('density')\n",
    "    ax.set_label(\"test\")\n",
    "    # ax.legend([\"original\", \"down shifted\"])\n",
    "fig.tight_layout()\n",
    "\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_normal')\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_normal.pdf')\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_normal_highres', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log transformations and errors\n",
    "\n",
    "- what does log2 transformation mean for the error\n",
    "\n",
    "If the error is calculated in log2 space, the larger values have to be\n",
    "predicted with higher precision (in comparison to the original space)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_original_error_log2(x: float, error_log2: float):\n",
    "    return 2 ** (np.log2(x) + error_log2) - x\n",
    "\n",
    "\n",
    "print(\n",
    "    f\"{get_original_error_log2(1e9, 0.5) = :,.1f}\",\n",
    "    f\"{get_original_error_log2(1e8, 0.5) = :,.1f}\",\n",
    "    sep='\\n'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we try to find the rel log2 error equalling the original error, this can be done by\n",
    "equating:\n",
    "\n",
    "$$ \\exp(\\ln(a)+e) - a = \\exp(\\ln(a)+e^*) - b $$\n",
    "\n",
    "Setting $a$, $e$ and $b$ we want to solve for $e^*$, which gives\n",
    "\n",
    "$$ e^* = \\ln \\left(\\frac{\\exp\\big(\\ln(a)+e\\big) - a + b}{a}  \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rel_error(measurment, log_error, other_measurment):\n",
    "    numerator = 2 ** (np.log2(measurment) + log_error)\n",
    "    numerator -= measurment\n",
    "    numerator += other_measurment\n",
    "\n",
    "    denominator = other_measurment\n",
    "    return np.log2(numerator / denominator)\n",
    "\n",
    "\n",
    "rel_error = rel_error(1.e9, 0.5, 1e8)\n",
    "print(f\"{rel_error = :.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f\"0.500 rel to 1e9: {get_original_error_log2(1e9, 0.5) :,.1f}\",\n",
    "    f\"{rel_error:.3f} rel to 1e8: {get_original_error_log2(1e8, rel_error) :,.1f}\",\n",
    "    sep='\\n'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the relative error of 0.5 for $10^9$ is five times larger for $10^8$ in the logspace,\n",
    "whereas the error in the original space is the same"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Volcano plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7883109",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Sample data for the volcano plot\n",
    "np.random.seed(42)\n",
    "fold_change = np.random.default_rng().normal(0, 1, 1000)\n",
    "p_value = np.random.default_rng().uniform(0, 1, 1000)\n",
    "\n",
    "# Volcano plot\n",
    "# Assuming you have two arrays, fold_change and p_value, containing the fold change values and p-values respectively\n",
    "\n",
    "# Set the significance threshold for p-value\n",
    "significance_threshold = 0.05\n",
    "\n",
    "# Set the fold change threshold\n",
    "fold_change_threshold = 2\n",
    "\n",
    "# Create a boolean mask for significant points\n",
    "significant_mask = (p_value < significance_threshold)\n",
    "\n",
    "# Create a boolean mask for points that meet the fold change threshold\n",
    "fold_change_mask = (abs(fold_change) > fold_change_threshold)\n",
    "\n",
    "# Combine the masks to get the final mask for significant points\n",
    "final_mask = significant_mask & fold_change_mask\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3, 3))\n",
    "# Plot the volcano plot\n",
    "_ = ax.scatter(fold_change, -np.log10(p_value), c='gray', alpha=0.5, s=10)\n",
    "_ = ax.scatter(fold_change[final_mask], -np.log10(p_value[final_mask]), c='red', alpha=0.7, s=20)\n",
    "\n",
    "# Add labels and title\n",
    "_ = ax.set_xlabel('Log2 fold change')\n",
    "_ = ax.set_ylabel('-log10(p-value)')\n",
    "_ = ax.set_title('Volcano Plot')\n",
    "\n",
    "# Add significance threshold lines\n",
    "_ = ax.axhline(-np.log10(significance_threshold), color='black', linestyle='--')\n",
    "_ = ax.axvline(fold_change_threshold, color='black', linestyle='--')\n",
    "_ = ax.axvline(-fold_change_threshold, color='black', linestyle='--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b65b4e01",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "fig.tight_layout()\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_volcano.png', dpi=300)\n",
    "fig.savefig(FIGUREFOLDER / 'illustration_volcano.pdf')  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c324f9d8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a3c193f136b3677f7a0d8b8f2344336d9d9fbcf8449c7b9fd96c2b5d44d3d77c"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}