{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replicate averaging\n",
    "====================\n",
    "\n",
    "Every `Profile` subclass (`rnav.data.Profile`, `SHAPEMaP`, `DanceMaP`, `RNPMaP`) optionally accepts a list of file names or rnav.Data objects as `input_data`.\n",
    "\n",
    "1. the active metric column is averaged across replicates\n",
    "2. the standard error of the mean is stored as `<metric_column>_sem` and registered as the `error_column`, so error bars appear automatically on profile and skyline plots\n",
    "\n",
    "This applies to the `shapemap`, `dancemap`, `rnpmap`, `profile` and other related [standard data keywords](../get_started/loading_data).\n",
    "\n",
    "1. [Averaging replicate profiles](#Averaging-replicate-profiles)\n",
    "2. [Loading replicates from files](#Loading-replicates-from-files)\n",
    "3. [Normalizing before averaging](#Normalizing-before-averaging)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import rnavigate as rnav"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Averaging replicate profiles\n",
    "-----------------------------\n",
    "\n",
    "Below, 3 synthetic replicate measurements are built around a shared \"true\" signal using `rnav.data.Profile.from_array`.\n",
    "Each behaves exactly like a `Profile` loaded from a file. Passing a list of them to `rnav.data.Profile` averages them the same way a list of `shapemap` `profile.txt` paths would."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "sequence = \"AUGC\" * 40  # 160 nt shared sequence\n",
    "true_signal = rng.random(160) * 2\n",
    "\n",
    "replicates = [\n",
    "    rnav.data.Profile.from_array(\n",
    "        input_data=np.clip(true_signal + rng.normal(scale=0.15, size=160), 0, None),\n",
    "        sequence=sequence,\n",
    "        metric=\"Profile\",\n",
    "    )\n",
    "    for _ in range(3)\n",
    "]\n",
    "\n",
    "combined = rnav.data.Profile(input_data=replicates, metric=\"Profile\")\n",
    "combined.data[[\"Nucleotide\", \"Profile\", \"Profile_sem\"]].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`combined.error_column` is now `\"Profile_sem\"`, so error bars are drawn automatically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_sample = rnav.Sample(sample=\"combined replicates\", my_profile=combined)\n",
    "\n",
    "plot = rnav.plot_profile(\n",
    "    samples=[my_sample],\n",
    "    profile=\"my_profile\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading replicates from files\n",
    "\n",
    "In practice, replicates usually come from files. Any profile-creating data keyword accepts a list of paths in place of a single file. RNAvigate loads each one, checks that they share a sequence, and averages them:\n",
    "\n",
    "```python\n",
    "my_sample = rnav.Sample(\n",
    "    sample=\"example\",\n",
    "    shapemap=[\n",
    "        \"replicate_1_profile.txt\",\n",
    "        \"replicate_2_profile.txt\",\n",
    "        \"replicate_3_profile.txt\",\n",
    "    ],\n",
    ")\n",
    "```\n",
    "\n",
    "This is equivalent to constructing each `SHAPEMaP` replicate individually and combining them with `SHAPEMaP.from_replicates([...])`. The same list syntax works for `dancemap` and `rnpmap`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalizing before averaging\n",
    "-------------------------------\n",
    "\n",
    "Replicates are sometimes on different scales, e.g. differing sequencing depth or overall reactivity.\n",
    "Pass `normalize_kwargs` to renormalize each replicate with `Profile.normalize(**normalize_kwargs)` *before* averaging, rather than only normalizing the final combined profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaled_replicates = [\n",
    "    rnav.data.Profile.from_array(\n",
    "        input_data=np.clip((true_signal + rng.normal(scale=0.15, size=160)) * scale, 0, None),\n",
    "        sequence=sequence,\n",
    "        metric=\"Profile\",\n",
    "    )\n",
    "    for scale in [0.5, 1.0, 1.5]  # 3 replicates on very different scales\n",
    "]\n",
    "\n",
    "raw_mean = rnav.data.Profile(\n",
    "    input_data=scaled_replicates, metric=\"Profile\"\n",
    ").data[\"Profile\"].mean()\n",
    "\n",
    "rescaled_mean = rnav.data.Profile(\n",
    "    input_data=scaled_replicates,\n",
    "    metric=\"Profile\",\n",
    "    normalize_kwargs={\"norm_method\": \"boxplot\", \"new_profile\": \"Profile\"},\n",
    ").data[\"Profile\"].mean()\n",
    "\n",
    "print(f\"averaged without renormalizing: {raw_mean:.3f}\")\n",
    "print(f\"averaged after per-replicate renormalizing: {rescaled_mean:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: If the metric has a `color_column` distinct from its metric column, RNAvigate warns\n",
    "that coloring falls back to the first replicate. Averaging does not currently combine colors.\n",
    "\n",
    "See also: [Custom profiles](custom_profiles) for other ways to construct `Profile` objects, and [standard data keywords](../get_started/loading_data) for the `shapemap`, `dancemap`, `rnpmap`, and `profile` keywords."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rnavigate (3.12.13.final.0)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}