ndoc_roc_gpu.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import cudf\n",
    "from numba import jit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   position   5mer  depth_tgt  depth_ref  med_current  mad_current  \\\n",
      "0         2  GAGGA          3          3          0.0          0.0   \n",
      "1         3  AGGAA          3          3          0.0          0.0   \n",
      "2         4  GGAAA          4          4          0.0          0.0   \n",
      "3         5  GAAAG          4          4          0.0          0.0   \n",
      "4         6  AAAGC          4          4          0.0          0.0   \n",
      "\n",
      "   med_currentR  mad_currentR  current_ratio  scoreSide1  scoreSide2  score  \n",
      "0           0.0           0.0            0.0         0.0         0.0    0.0  \n",
      "1           0.0           0.0            0.0         0.0         0.0    0.0  \n",
      "2           0.0           0.0            0.0         0.0         0.0    0.0  \n",
      "3           0.0           0.0            0.0         0.0         0.0    0.0  \n",
      "4           0.0           0.0            0.0         0.0         0.0    0.0  \n",
      "(1936, 12)\n"
     ]
    }
   ],
   "source": [
    "nanodoc_16s = cudf.read_csv(\"/home/bhargavam/Documents/data/ndoc_final/ndoc_16s_1.txt\", sep='\\t')\n",
    "nanodoc_16s.columns = ['position', '5mer', 'depth_tgt', 'depth_ref', 'med_current', 'mad_current', 'med_currentR', 'mad_currentR', 'current_ratio', 'scoreSide1', 'scoreSide2', 'score']\n",
    "\n",
    "print(nanodoc_16s.head())\n",
    "print(nanodoc_16s.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "known_mods_16s = [716, 727, 1166, 1167, 1407, 1602, 1607, 1698, 1716, 1718, 1719]\n",
    "known_mods_16s_pos = [x+w for x in known_mods_16s]\n",
    "known_mods_16s_neg = [x-w for x in known_mods_16s]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def conditions_16s(a):\n",
    "    if a in (known_mods_16s + known_mods_16s_pos + known_mods_16s_neg):\n",
    "        return 1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "user defined function compilation failed.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/site-packages/cudf/core/indexed_frame.py:1821\u001b[0m, in \u001b[0;36mIndexedFrame._apply\u001b[0;34m(self, func, kernel_getter, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1820\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m-> 1821\u001b[0m     kernel, retty \u001b[39m=\u001b[39m _compile_or_get(\n\u001b[1;32m   1822\u001b[0m         \u001b[39mself\u001b[39;49m, func, args, kernel_getter\u001b[39m=\u001b[39;49mkernel_getter\n\u001b[1;32m   1823\u001b[0m     )\n\u001b[1;32m   1824\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mException\u001b[39;00m \u001b[39mas\u001b[39;00m e:\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/contextlib.py:79\u001b[0m, in \u001b[0;36mContextDecorator.__call__.<locals>.inner\u001b[0;34m(*args, **kwds)\u001b[0m\n\u001b[1;32m     78\u001b[0m \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_recreate_cm():\n\u001b[0;32m---> 79\u001b[0m     \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwds)\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/site-packages/cudf/core/udf/utils.py:207\u001b[0m, in \u001b[0;36m_compile_or_get\u001b[0;34m(frame, func, args, kernel_getter)\u001b[0m\n\u001b[1;32m    206\u001b[0m \u001b[39m# check to see if we already compiled this function\u001b[39;00m\n\u001b[0;32m--> 207\u001b[0m cache_key \u001b[39m=\u001b[39m _generate_cache_key(frame, func)\n\u001b[1;32m    208\u001b[0m \u001b[39mif\u001b[39;00m precompiled\u001b[39m.\u001b[39mget(cache_key) \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/site-packages/cudf/core/udf/utils.py:175\u001b[0m, in \u001b[0;36m_generate_cache_key\u001b[0;34m(frame, func)\u001b[0m\n\u001b[1;32m    167\u001b[0m \u001b[39m\"\"\"Create a cache key that uniquely identifies a compilation.\u001b[39;00m\n\u001b[1;32m    168\u001b[0m \n\u001b[1;32m    169\u001b[0m \u001b[39mA new compilation is needed any time any of the following things change:\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    172\u001b[0m \u001b[39m- The existence of the input columns masks\u001b[39;00m\n\u001b[1;32m    173\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m    174\u001b[0m \u001b[39mreturn\u001b[39;00m (\n\u001b[0;32m--> 175\u001b[0m     \u001b[39m*\u001b[39mcudautils\u001b[39m.\u001b[39;49mmake_cache_key(\n\u001b[1;32m    176\u001b[0m         func, \u001b[39mtuple\u001b[39;49m(_all_dtypes_from_frame(frame)\u001b[39m.\u001b[39;49mvalues())\n\u001b[1;32m    177\u001b[0m     ),\n\u001b[1;32m    178\u001b[0m     \u001b[39m*\u001b[39m(col\u001b[39m.\u001b[39mmask \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m \u001b[39mfor\u001b[39;00m col \u001b[39min\u001b[39;00m frame\u001b[39m.\u001b[39m_data\u001b[39m.\u001b[39mvalues()),\n\u001b[1;32m    179\u001b[0m     \u001b[39m*\u001b[39mframe\u001b[39m.\u001b[39m_data\u001b[39m.\u001b[39mkeys(),\n\u001b[1;32m    180\u001b[0m )\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/site-packages/cudf/utils/cudautils.py:202\u001b[0m, in \u001b[0;36mmake_cache_key\u001b[0;34m(udf, sig)\u001b[0m\n\u001b[1;32m    200\u001b[0m names \u001b[39m=\u001b[39m udf\u001b[39m.\u001b[39m\u001b[39m__code__\u001b[39m\u001b[39m.\u001b[39mco_names\n\u001b[0;32m--> 202\u001b[0m \u001b[39mif\u001b[39;00m udf\u001b[39m.\u001b[39;49m\u001b[39m__closure__\u001b[39;49m \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m    203\u001b[0m     cvars \u001b[39m=\u001b[39m \u001b[39mtuple\u001b[39m(x\u001b[39m.\u001b[39mcell_contents \u001b[39mfor\u001b[39;00m x \u001b[39min\u001b[39;00m udf\u001b[39m.\u001b[39m\u001b[39m__closure__\u001b[39m)\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'CPUDispatcher' object has no attribute '__closure__'",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[19], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m nanodoc_16s[\u001b[39m'\u001b[39m\u001b[39mmods\u001b[39m\u001b[39m'\u001b[39m] \u001b[39m=\u001b[39m nanodoc_16s[\u001b[39m'\u001b[39;49m\u001b[39mposition\u001b[39;49m\u001b[39m'\u001b[39;49m]\u001b[39m.\u001b[39;49mapply(conditions_16s)\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/contextlib.py:79\u001b[0m, in \u001b[0;36mContextDecorator.__call__.<locals>.inner\u001b[0;34m(*args, **kwds)\u001b[0m\n\u001b[1;32m     76\u001b[0m \u001b[39m@wraps\u001b[39m(func)\n\u001b[1;32m     77\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39minner\u001b[39m(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwds):\n\u001b[1;32m     78\u001b[0m     \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_recreate_cm():\n\u001b[0;32m---> 79\u001b[0m         \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwds)\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/site-packages/cudf/core/series.py:2382\u001b[0m, in \u001b[0;36mSeries.apply\u001b[0;34m(self, func, convert_dtype, args, **kwargs)\u001b[0m\n\u001b[1;32m   2379\u001b[0m \u001b[39mif\u001b[39;00m convert_dtype \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mTrue\u001b[39;00m:\n\u001b[1;32m   2380\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mSeries.apply only supports convert_dtype=True\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m-> 2382\u001b[0m result \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_apply(func, _get_scalar_kernel, \u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m   2383\u001b[0m result\u001b[39m.\u001b[39mname \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mname\n\u001b[1;32m   2384\u001b[0m \u001b[39mreturn\u001b[39;00m result\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/contextlib.py:79\u001b[0m, in \u001b[0;36mContextDecorator.__call__.<locals>.inner\u001b[0;34m(*args, **kwds)\u001b[0m\n\u001b[1;32m     76\u001b[0m \u001b[39m@wraps\u001b[39m(func)\n\u001b[1;32m     77\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39minner\u001b[39m(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwds):\n\u001b[1;32m     78\u001b[0m     \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_recreate_cm():\n\u001b[0;32m---> 79\u001b[0m         \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwds)\n",
      "File \u001b[0;32m~/mambaforge/lib/python3.9/site-packages/cudf/core/indexed_frame.py:1825\u001b[0m, in \u001b[0;36mIndexedFrame._apply\u001b[0;34m(self, func, kernel_getter, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1821\u001b[0m     kernel, retty \u001b[39m=\u001b[39m _compile_or_get(\n\u001b[1;32m   1822\u001b[0m         \u001b[39mself\u001b[39m, func, args, kernel_getter\u001b[39m=\u001b[39mkernel_getter\n\u001b[1;32m   1823\u001b[0m     )\n\u001b[1;32m   1824\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mException\u001b[39;00m \u001b[39mas\u001b[39;00m e:\n\u001b[0;32m-> 1825\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\n\u001b[1;32m   1826\u001b[0m         \u001b[39m\"\u001b[39m\u001b[39muser defined function compilation failed.\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m   1827\u001b[0m     ) \u001b[39mfrom\u001b[39;00m \u001b[39me\u001b[39;00m\n\u001b[1;32m   1829\u001b[0m \u001b[39m# Mask and data column preallocated\u001b[39;00m\n\u001b[1;32m   1830\u001b[0m ans_col \u001b[39m=\u001b[39m _return_arr_from_dtype(retty, \u001b[39mlen\u001b[39m(\u001b[39mself\u001b[39m))\n",
      "\u001b[0;31mValueError\u001b[0m: user defined function compilation failed."
     ]
    }
   ],
   "source": [
    "nanodoc_16s['mods'] = nanodoc_16s['position'].apply(conditions_16s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[31mERROR: \u001b[0mThe function received no value for the required argument: w\n",
      "Usage: ipykernel_launcher.py FILENAME W\n",
      "\n",
      "For detailed information on this command, run:\n",
      "  ipykernel_launcher.py --help\n"
     ]
    },
    {
     "ename": "FireExit",
     "evalue": "2",
     "output_type": "error",
     "traceback": [
      "An exception has occurred, use %tb to see the full traceback.\n",
      "\u001b[0;31mFireExit\u001b[0m\u001b[0;31m:\u001b[0m 2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/bhargavam/mambaforge/lib/python3.9/site-packages/IPython/core/interactiveshell.py:3441: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.\n",
      "  warn(\"To exit: use 'exit', 'quit', or Ctrl-D.\", stacklevel=1)\n"
     ]
    }
   ],
   "source": [
    "## define function to generalize the CSV generation for ROC curve generation for 16s rRNA\n",
    "\n",
    "def roc(filename, w):\n",
    "    # read in the output file from nanodoc analysis for 16s rRNA\n",
    "    \n",
    "    nanodoc_16s = cudf.read_csv(filename, sep='\\t')\n",
    "    nanodoc_16s.columns = ['position', '5mer', 'depth_tgt', 'depth_ref', 'med_current', 'mad_current', 'med_currentR', 'mad_currentR', 'current_ratio', 'scoreSide1', 'scoreSide2', 'score']\n",
    "\n",
    "    ## calculate values for ROC for 16s rRNA - add mods column with either yes for known modifications or no for all other positions\n",
    "\n",
    "    known_mods_16s = [716, 727, 1166, 1167, 1407, 1602, 1607, 1698, 1716, 1718, 1719]\n",
    "    known_mods_16s_pos = [x+w for x in known_mods_16s]\n",
    "    known_mods_16s_neg = [x-w for x in known_mods_16s]\n",
    "\n",
    "    def conditions_16s(a):\n",
    "        if a in (known_mods_16s + known_mods_16s_pos + known_mods_16s_neg):\n",
    "            return 1\n",
    "        else:\n",
    "            return 0\n",
    "        \n",
    "    nanodoc_16s['mods'] = nanodoc_16s['position'].apply(conditions_16s)\n",
    "\n",
    "    ## sort the rows based on the log values\n",
    "\n",
    "    sorted_16s = nanodoc_16s.sort_values('score')\n",
    "\n",
    "    ## calculate variables required for ROC visualization\n",
    "\n",
    "    thresholds = list(abs(sorted_16s['score']))\n",
    "\n",
    "    roc_point = []\n",
    "\n",
    "    for threshold in thresholds:\n",
    "            \n",
    "        tp = 0; fp = 0; fn = 0; tn = 0\n",
    "\n",
    "        for index, instance in sorted_16s.iterrows():\n",
    "            actual = instance[\"mods\"]\n",
    "            prediction = abs(instance[\"score\"])\n",
    "\n",
    "            #print(actual, prediction)\n",
    "            \n",
    "            if prediction >= threshold:\n",
    "                prediction_class = 1\n",
    "            else:\n",
    "                prediction_class = 0\n",
    "\n",
    "            if prediction_class == 1 and actual == 1:\n",
    "                tp = tp + 1\n",
    "            elif actual == 1 and prediction_class == 0:\n",
    "                fn = fn + 1\n",
    "            elif actual == 0 and prediction_class == 1: \n",
    "                fp = fp + 1\n",
    "            elif actual == 0 and prediction_class == 0:\n",
    "                tn = tn + 1\n",
    "\n",
    "        #print(tp, fp, fn, tn)\n",
    "\n",
    "        tpr = tp / (tp + fn)\n",
    "        fpr = fp / (tn + fp)\n",
    "\n",
    "        #print(tpr, fpr)\n",
    "\n",
    "        roc_point.append([tpr, fpr])\n",
    "    \n",
    "    ## write out ROC values into a dataframe\n",
    "\n",
    "    pivot = cudf.DataFrame(roc_point, columns=[\"tpr\", \"fpr\"])\n",
    "    pivot[\"threshold\"] = thresholds\n",
    "    pivot.to_csv(filename + \"_roc_\" + str(w) + \"nt_window.tsv\")\n",
    "    \n",
    "    return\n",
    "\n",
    "if __name__ == '__main__':\n",
    "  fire.Fire(roc)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e267e0a521d3dc2ae2033b479444b4e5aa0019554de3a258d9d3d1db7a2e250a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}