{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbcaae60",
   "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": "0be1234d",
   "metadata": {},
   "outputs": [],
   "source": [
    "kmers = list(map(int,\"11 13 15 17 19 21 25 31 41 51 61 81 101 151 201 251 301 501\".split()))\n",
    "\n",
    "data_kmers = []\n",
    "for k in kmers:\n",
    "    abundances = {}\n",
    "    mean = 0\n",
    "    nb_items = 0\n",
    "    genome_abundance_1_total = 0\n",
    "    for x,line in enumerate(open(\"chm13v2.0-\"+str(k)+\".histo2D\")):\n",
    "        l = list(map(int,line.split()[1:]))\n",
    "        assert(int(line.split()[0].strip(\":\")) == x)\n",
    "        abundances[x]=l\n",
    "        genome_abundance_1_total += l[1]\n",
    "        if x == 0:\n",
    "            genome_abundance_1_missing = l[1]\n",
    "    data_kmers += [(k,genome_abundance_1_missing/genome_abundance_1_total)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b7479b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "df = pd.DataFrame([x for x in data_kmers],columns=(\"k\",\"nb_kmers_abundance_1_in_genome\"))\n",
    "\n",
    "print(df)\n",
    "sns_plot = sns.lineplot(x=df[\"k\"], y=df[\"nb_kmers_abundance_1_in_genome\"])\n",
    "\n",
    "#sns_plot.set(xscale=\"log\")\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",
    "\n",
    "ylabels = ['{:}%'.format(int(x*100))  for x in sns_plot.get_yticks()]\n",
    "sns_plot.set_yticklabels(ylabels)\n",
    "\n",
    "#plt.xticks(sns_plot.get_xticks(), sns_plot.get_xticks() * 100)\n",
    "sns_plot.set(xlabel='k value', ylabel='')\n",
    "#sns_plot.get_legend().remove()\n",
    "plt.title('Uniquely occurring reference k-mers \\n missing from HiFi DeepConsensus reads')\n",
    "\n",
    "fig = sns_plot.get_figure()\n",
    "fig.savefig(\"seeds_kmers.pdf\", bbox_inches='tight', dpi =300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44f380ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "ksizes = list(map(int,\"2 4 6 8 10\".split()))\n",
    "\n",
    "data_kminmers = []\n",
    "for k in ksizes:\n",
    "    abundances = {}\n",
    "    mean = 0\n",
    "    nb_items = 0\n",
    "    genome_abundance_1_total = 0\n",
    "    for x,line in enumerate(open(f\"2DHisto-k{k}-d0.01-l31.hist2D\")):\n",
    "        l = list(map(int,line.split()))\n",
    "        abundances[x]=l\n",
    "        genome_abundance_1_total += l[1]\n",
    "        if x == 0:\n",
    "            genome_abundance_1_missing = l[1]\n",
    "    data_kminmers += [(k,genome_abundance_1_missing/genome_abundance_1_total)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1437142b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame([x for x in data_kminmers],columns=(\"k\",\"nb_kminmers_abundance_1_in_genome\"))\n",
    "\n",
    "print(df)\n",
    "sns_plot = sns.lineplot(x=df[\"k\"], y=df[\"nb_kminmers_abundance_1_in_genome\"])\n",
    "\n",
    "#sns_plot.set(xscale=\"log\")\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",
    "\n",
    "ylabels = ['{:}%'.format(int(x*100))  for x in sns_plot.get_yticks()]\n",
    "sns_plot.set_yticklabels(ylabels)\n",
    "\n",
    "#plt.xticks(sns_plot.get_xticks(), sns_plot.get_xticks() * 100)\n",
    "sns_plot.set(xlabel='k value', ylabel='')\n",
    "#sns_plot.get_legend().remove()\n",
    "plt.title('Uniquely occurring reference k-min-mers \\n missing from HiFi DeepConsensus reads')\n",
    "\n",
    "\n",
    "fig = sns_plot.get_figure()\n",
    "\n",
    "# create second Axes. Note the 0.0 height\n",
    "ax2 = fig.add_axes((0.15,-0.1,0.72,0.0))\n",
    "ax2.yaxis.set_visible(False) # hide the yaxis\n",
    "new_tick_locations = np.array([0,.25, .5,.75,1])\n",
    "def tick_function(X):\n",
    "    ksizes = list(map(int,\" 2 4 6 8 10\".split()))\n",
    "    density=0.01\n",
    "    return [\"%d\" % (ksizes[i]/density) for i,z in enumerate(X)]\n",
    "ax2.set_xticks(new_tick_locations)\n",
    "ax2.set_xticklabels(tick_function(new_tick_locations))\n",
    "ax2.set_xlabel(\"Equivalent k span in bases\")\n",
    "plt.show()\n",
    "\n",
    "fig.savefig(\"seeds_kminers.pdf\", bbox_inches='tight', dpi =300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbe9da56",
   "metadata": {
    "scrolled": true
   },
   "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 x,y in data_kmers:\n",
    "    data[x]=(x,y,\"k-mers\")\n",
    "for x,y in data_kminmers:\n",
    "    kspan = k_to_kspan(x)\n",
    "    kweight = k_to_kweight(x)\n",
    "    data[kspan] = (kweight,y,\"k-min-mers\")\n",
    "\n",
    "ktype = \"kspan\"\n",
    "ktype = \"kweight\"\n",
    "\n",
    "print(data)\n",
    "df = pd.DataFrame([(x,data[x][0],data[x][1],data[x][2]) for x in data],columns=(\"kspan\",\"kweight\",\"perc_kmers_abundance_1_in_genome\",\"type\"))\n",
    "  \n",
    "sns_plot = sns.lineplot(df,x=ktype,y=\"perc_kmers_abundance_1_in_genome\",hue=\"type\")\n",
    "\n",
    "#sns_plot.set(xscale=\"log\")\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",
    "\n",
    "ylabels = ['{:}%'.format(int(x*100))  for x in sns_plot.get_yticks()]\n",
    "sns_plot.set_yticklabels(ylabels)\n",
    "\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.title('Percentage of seeds that occur uniquely \\n in the genome yet are absent from reads')\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.3,-0.1,0.56,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 4 6 8 10\".split()))\n",
    "        return [\"%d\" % (ksizes[i]) for i,z in enumerate(X)]\n",
    "    new_tick_locations = np.array([0,.25, .5,.75,1])\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_2D.pdf\", bbox_inches='tight', dpi =300)\n"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
