{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2035d3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "import statistics\n",
    "plt.rc('figure', facecolor='w', figsize=[4,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "901095dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "ksizes = list(map(int,\"11 15 20 25 50 100 200 500\".split()))\n",
    "\n",
    "data_jellyfish = []\n",
    "for k in ksizes:\n",
    "    nb_items = 0\n",
    "    nb_missing = 0\n",
    "    for line in open(f\"jellyfish-k{k}.reads_stats\"):\n",
    "        median,mean,stdev = map(float,line.split())\n",
    "        data_jellyfish += [(k,median,mean,stdev)]\n",
    "\n",
    "ksizes = list(map(int,\"2 4 6 8 10 12 15 20\".split()))\n",
    "\n",
    "data_kminmers = []\n",
    "for k in ksizes:\n",
    "    nb_items = 0\n",
    "    nb_missing = 0\n",
    "    for line in open(f\"HG002-kminmer-abundance-in-CHM13v2-k{k}.50k.read_stats\"):\n",
    "        median,mean,stdev = map(float,line.split())\n",
    "        data_kminmers  += [(k,median,mean,stdev)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53b8843b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame([x for x in data_jellyfish],columns=(\"k\",\"median\",\"mean\",\"stdev\"))\n",
    "print(df)\n",
    "sns_plot = sns.lineplot(data=df,x=\"k\", y=\"median\", errorbar=('se',10))\n",
    "sns_plot.set(yscale=\"log\")\n",
    "#ylabels = ['{:}%'.format(int(x*100))  for x in sns_plot.get_yticks()]\n",
    "#sns_plot.set_yticklabels(ylabels)\n",
    "#plt.xticks(sns_plot.get_xticks(), list(map(lambda x:\"%d\"%x,sns_plot.get_xticks())))\n",
    "sns_plot.set(xlabel='k value', ylabel='')\n",
    "plt.title('Number of occurrences in CHM13 k-mers per read in the HG002 DeepConsensus dataset')\n",
    "sns_plot.axhline(1)\n",
    "fig = sns_plot.get_figure()\n",
    "plt.show()\n",
    "\n",
    "fig.savefig(\"seeds_abundances_ref_kmers.pdf\", bbox_inches='tight', dpi =300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "953651a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame([x for x in data_kminmers],columns=(\"k\",\"median\",\"mean\",\"stdev\"))\n",
    "print(df)\n",
    "sns_plot = sns.lineplot(data=df,x=\"k\", y=\"median\", errorbar=('se',10))\n",
    "sns_plot.set(yscale=\"log\")\n",
    "#ylabels = ['{:}%'.format(int(x*100))  for x in sns_plot.get_yticks()]\n",
    "#sns_plot.set_yticklabels(ylabels)\n",
    "#plt.xticks(sns_plot.get_xticks(), list(map(lambda x:\"%d\"%x,sns_plot.get_xticks())))\n",
    "sns_plot.set(xlabel='k value', ylabel='')\n",
    "plt.title('Number of occurrences in CHM13 k-min-mers per read in the HG002 DeepConsensus dataset')\n",
    "sns_plot.axhline(1)\n",
    "fig = sns_plot.get_figure()\n",
    "plt.show()\n",
    "\n",
    "fig.savefig(\"seeds_abundances_ref_kminmers.pdf\", bbox_inches='tight', dpi =300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f152652",
   "metadata": {},
   "outputs": [],
   "source": [
    "#attempt to merge both plots\n",
    "\n",
    "def k_to_kspan(k):\n",
    "    density=0.01\n",
    "    return k/density\n",
    "\n",
    "def k_to_kweight(k):\n",
    "    l=31\n",
    "    return k*l\n",
    "\n",
    "data = []\n",
    "for (k,median,mean,stdev) in data_jellyfish:\n",
    "    data += [(k,k,median,\"k-mers\")]\n",
    "for (k,median,mean,stdev) in data_kminmers:\n",
    "    kspan = k_to_kspan(k)\n",
    "    kweight = k_to_kweight(k)\n",
    "    data += [(kspan,kweight,median,\"k-min-mers\")]\n",
    "\n",
    "ktype = \"kspan\"\n",
    "ktype = \"kweight\"\n",
    "df = pd.DataFrame(data,columns=(\"kspan\",\"kweight\",\"median\",\"type\"))\n",
    "  \n",
    "sns_plot = sns.lineplot(df,x=ktype, y=\"median\",hue=\"type\", errorbar=('se',10))\n",
    "sns_plot.set(yscale=\"log\")\n",
    "sns_plot.set(xscale=\"log\")\n",
    "plt.ylim(0.15, 10000)\n",
    "from matplotlib.ticker import FuncFormatter\n",
    "f = lambda x, pos: f'{int(x)}'\n",
    "sns_plot.yaxis.set_major_formatter(FuncFormatter(f))\n",
    "sns_plot.xaxis.set_major_formatter(FuncFormatter(f))\n",
    "sns_plot.set(xticks=[10,100,500])\n",
    "#sns_plot.set(xticks=[x/100 for x in range(0,10+1,2)])\n",
    "#sns_plot.set_xticklabels([str(x-2) for x in range(0,10+3,2)])\n",
    "#plt.xticks(sns_plot.get_xticks(), sns_plot.get_xticks() * 100)\n",
    "if ktype == \"kspan\":\n",
    "    sns_plot.set(xlabel='Equivalent k-mer length', ylabel='')\n",
    "else:\n",
    "    sns_plot.set(xlabel='Seed weight', ylabel='')\n",
    "sns_plot.get_legend().remove()\n",
    "plt.legend(loc='upper right')\n",
    "plt.title('Median number of occurrences of seeds \\n from the reads in the reference genome')\n",
    "sns_plot.axhline(1, color='k', linestyle='--')\n",
    "sns_plot.axvline(19, color='steelblue', linestyle='--')\n",
    "if ktype == \"kspan\":\n",
    "    sns_plot.axvline(k_to_kspan(7), color='orange', linestyle='--')\n",
    "else:\n",
    "    sns_plot.axvline(k_to_kweight(7), color='orange', linestyle='--')\n",
    "if ktype == \"kspan\":\n",
    "    plt.text(85, 225, \"minimap2's k\", horizontalalignment='right', size='small', color='steelblue')\n",
    "else:\n",
    "    plt.text(65, 225, \"minimap2's k\", horizontalalignment='right', size='small', color='steelblue')\n",
    "if ktype == \"kspan\":\n",
    "    plt.text(2500, 5, \"mapquik's k\", horizontalalignment='right', size='small', color='darkgoldenrod')\n",
    "else:\n",
    "    plt.text(650, 5, \"mapquik's k\", horizontalalignment='right', size='small', color='darkgoldenrod')\n",
    "\n",
    "fig = sns_plot.get_figure()\n",
    "if ktype == \"kspan\":\n",
    "    # create second Axes. Note the 0.0 height\n",
    "    ax2 = fig.add_axes((0.55,-0.1,0.32,0.0))  # i didn't determine those params automatically, these are just empirically set\n",
    "    ax2.yaxis.set_visible(False) # hide the yaxis\n",
    "    def tick_function(X):\n",
    "        ksizes = list(map(int,\"2 5 10 20\".split()))\n",
    "        return [\"%d\" % (ksizes[i]) for i,z in enumerate(X)]\n",
    "    import math\n",
    "    new_tick_locations = list(map(lambda x:math.log(x)/math.log(1.01),np.array([.1,0.2, .5,1.01])))\n",
    "    ax2.set_xticks(new_tick_locations)\n",
    "    ax2.set_xticklabels(tick_function(new_tick_locations))\n",
    "    ax2.set_xlabel(\"k-min-mer's k value (density=0.01)\")\n",
    "plt.show()\n",
    "\n",
    "fig.savefig(\"seeds_abundances_to_ref.pdf\", bbox_inches='tight', dpi =300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d405308e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80f3a927",
   "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
