{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cluster Differentially Methylated Genes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-16T03:38:08.524584Z",
     "start_time": "2022-02-16T03:38:06.124950Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from ALLCools.clustering import one_vs_rest_dmg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-16T03:39:54.199362Z",
     "start_time": "2022-02-16T03:39:54.195979Z"
    }
   },
   "outputs": [],
   "source": [
    "mcds_paths = 'geneslop2k_frac.mcds'\n",
    "cell_meta_path = '../step_by_step/100kb/L1.ClusteringResults.csv.gz'\n",
    "cluster_col = 'L1'\n",
    "\n",
    "obs_dim = 'cell'\n",
    "var_dim = 'geneslop2k'\n",
    "mc_type = 'CHN'\n",
    "\n",
    "top_n = 1000\n",
    "auroc_cutoff = 0.8\n",
    "adj_p_cutoff = 0.001\n",
    "fc_cutoff = 0.8\n",
    "max_cluster_cells = 2000\n",
    "max_other_fold = 5\n",
    "cpu = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-16T03:38:10.251633Z",
     "start_time": "2022-02-16T03:38:10.100034Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AllcPath</th>\n",
       "      <th>mCCCFrac</th>\n",
       "      <th>mCGFrac</th>\n",
       "      <th>mCGFracAdj</th>\n",
       "      <th>mCHFrac</th>\n",
       "      <th>mCHFracAdj</th>\n",
       "      <th>FinalReads</th>\n",
       "      <th>InputReads</th>\n",
       "      <th>MappedReads</th>\n",
       "      <th>DissectionRegion</th>\n",
       "      <th>...</th>\n",
       "      <th>Sample</th>\n",
       "      <th>leiden</th>\n",
       "      <th>mCHFrac.1</th>\n",
       "      <th>tsne_0</th>\n",
       "      <th>tsne_1</th>\n",
       "      <th>L1</th>\n",
       "      <th>L1_proba</th>\n",
       "      <th>CellTypeAnno</th>\n",
       "      <th>umap_0</th>\n",
       "      <th>umap_1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>10E_M_0</th>\n",
       "      <td>/gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...</td>\n",
       "      <td>0.008198</td>\n",
       "      <td>0.822633</td>\n",
       "      <td>0.821166</td>\n",
       "      <td>0.041640</td>\n",
       "      <td>0.033718</td>\n",
       "      <td>1626504.0</td>\n",
       "      <td>4407752</td>\n",
       "      <td>2892347.0</td>\n",
       "      <td>10E</td>\n",
       "      <td>...</td>\n",
       "      <td>10E_190625</td>\n",
       "      <td>13</td>\n",
       "      <td>0.041640</td>\n",
       "      <td>57.602540</td>\n",
       "      <td>-5.024663</td>\n",
       "      <td>c11</td>\n",
       "      <td>0.864367</td>\n",
       "      <td>MGE-Sst</td>\n",
       "      <td>5.288734</td>\n",
       "      <td>9.726882</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10E_M_1</th>\n",
       "      <td>/gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...</td>\n",
       "      <td>0.006019</td>\n",
       "      <td>0.743035</td>\n",
       "      <td>0.741479</td>\n",
       "      <td>0.024127</td>\n",
       "      <td>0.018218</td>\n",
       "      <td>2009998.0</td>\n",
       "      <td>5524084</td>\n",
       "      <td>3657352.0</td>\n",
       "      <td>10E</td>\n",
       "      <td>...</td>\n",
       "      <td>10E_190625</td>\n",
       "      <td>11</td>\n",
       "      <td>0.024127</td>\n",
       "      <td>-45.191850</td>\n",
       "      <td>-11.135287</td>\n",
       "      <td>c7</td>\n",
       "      <td>0.669400</td>\n",
       "      <td>CA3</td>\n",
       "      <td>-3.702348</td>\n",
       "      <td>7.514084</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10E_M_10</th>\n",
       "      <td>/gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...</td>\n",
       "      <td>0.006569</td>\n",
       "      <td>0.750172</td>\n",
       "      <td>0.748520</td>\n",
       "      <td>0.027665</td>\n",
       "      <td>0.021235</td>\n",
       "      <td>1383636.0</td>\n",
       "      <td>3455260</td>\n",
       "      <td>2172987.0</td>\n",
       "      <td>10E</td>\n",
       "      <td>...</td>\n",
       "      <td>10E_190625</td>\n",
       "      <td>11</td>\n",
       "      <td>0.027665</td>\n",
       "      <td>-46.905564</td>\n",
       "      <td>-8.491459</td>\n",
       "      <td>c7</td>\n",
       "      <td>0.787267</td>\n",
       "      <td>CA3</td>\n",
       "      <td>-2.797569</td>\n",
       "      <td>7.604081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10E_M_101</th>\n",
       "      <td>/gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...</td>\n",
       "      <td>0.006353</td>\n",
       "      <td>0.760898</td>\n",
       "      <td>0.759369</td>\n",
       "      <td>0.026547</td>\n",
       "      <td>0.020323</td>\n",
       "      <td>2474670.0</td>\n",
       "      <td>7245482</td>\n",
       "      <td>4778768.0</td>\n",
       "      <td>10E</td>\n",
       "      <td>...</td>\n",
       "      <td>10E_190625</td>\n",
       "      <td>11</td>\n",
       "      <td>0.026547</td>\n",
       "      <td>-53.480022</td>\n",
       "      <td>-1.604433</td>\n",
       "      <td>c7</td>\n",
       "      <td>0.526933</td>\n",
       "      <td>CA3</td>\n",
       "      <td>-0.310848</td>\n",
       "      <td>8.465321</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10E_M_102</th>\n",
       "      <td>/gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...</td>\n",
       "      <td>0.005409</td>\n",
       "      <td>0.752980</td>\n",
       "      <td>0.751637</td>\n",
       "      <td>0.019497</td>\n",
       "      <td>0.014164</td>\n",
       "      <td>2430290.0</td>\n",
       "      <td>7004754</td>\n",
       "      <td>4609570.0</td>\n",
       "      <td>10E</td>\n",
       "      <td>...</td>\n",
       "      <td>10E_190625</td>\n",
       "      <td>7</td>\n",
       "      <td>0.019497</td>\n",
       "      <td>-25.967990</td>\n",
       "      <td>13.813133</td>\n",
       "      <td>c30</td>\n",
       "      <td>0.924000</td>\n",
       "      <td>CA1</td>\n",
       "      <td>0.252257</td>\n",
       "      <td>-3.450731</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 27 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                    AllcPath  mCCCFrac  \\\n",
       "10E_M_0    /gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...  0.008198   \n",
       "10E_M_1    /gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...  0.006019   \n",
       "10E_M_10   /gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...  0.006569   \n",
       "10E_M_101  /gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...  0.006353   \n",
       "10E_M_102  /gale/raidix/rdx-4/mapping/10E/CEMBA190625-10E...  0.005409   \n",
       "\n",
       "            mCGFrac  mCGFracAdj   mCHFrac  mCHFracAdj  FinalReads  InputReads  \\\n",
       "10E_M_0    0.822633    0.821166  0.041640    0.033718   1626504.0     4407752   \n",
       "10E_M_1    0.743035    0.741479  0.024127    0.018218   2009998.0     5524084   \n",
       "10E_M_10   0.750172    0.748520  0.027665    0.021235   1383636.0     3455260   \n",
       "10E_M_101  0.760898    0.759369  0.026547    0.020323   2474670.0     7245482   \n",
       "10E_M_102  0.752980    0.751637  0.019497    0.014164   2430290.0     7004754   \n",
       "\n",
       "           MappedReads DissectionRegion  ...      Sample  leiden mCHFrac.1  \\\n",
       "10E_M_0      2892347.0              10E  ...  10E_190625      13  0.041640   \n",
       "10E_M_1      3657352.0              10E  ...  10E_190625      11  0.024127   \n",
       "10E_M_10     2172987.0              10E  ...  10E_190625      11  0.027665   \n",
       "10E_M_101    4778768.0              10E  ...  10E_190625      11  0.026547   \n",
       "10E_M_102    4609570.0              10E  ...  10E_190625       7  0.019497   \n",
       "\n",
       "              tsne_0     tsne_1   L1  L1_proba CellTypeAnno    umap_0  \\\n",
       "10E_M_0    57.602540  -5.024663  c11  0.864367      MGE-Sst  5.288734   \n",
       "10E_M_1   -45.191850 -11.135287   c7  0.669400          CA3 -3.702348   \n",
       "10E_M_10  -46.905564  -8.491459   c7  0.787267          CA3 -2.797569   \n",
       "10E_M_101 -53.480022  -1.604433   c7  0.526933          CA3 -0.310848   \n",
       "10E_M_102 -25.967990  13.813133  c30  0.924000          CA1  0.252257   \n",
       "\n",
       "             umap_1  \n",
       "10E_M_0    9.726882  \n",
       "10E_M_1    7.514084  \n",
       "10E_M_10   7.604081  \n",
       "10E_M_101  8.465321  \n",
       "10E_M_102 -3.450731  \n",
       "\n",
       "[5 rows x 27 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cell_meta = pd.read_csv('../../cell_level/step_by_step/100kb/L1.ClusteringResults.csv.gz', index_col=0)\n",
    "cell_meta.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate DMG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-16T03:42:10.661960Z",
     "start_time": "2022-02-16T03:39:56.597735Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating cluster c0 DMGs.\n",
      "Calculating cluster c1 DMGs.\n",
      "Calculating cluster c10 DMGs.\n",
      "Calculating cluster c11 DMGs.\n",
      "Calculating cluster c12 DMGs.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/hanliu/miniconda3/envs/allcools_new/lib/python3.8/site-packages/xarray/core/indexing.py:1227: PerformanceWarning: Slicing is producing a large chunk. To accept the large\n",
      "chunk and silence this warning, set the option\n",
      "    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):\n",
      "    ...     array[indexer]\n",
      "\n",
      "To avoid creating the large chunks, set the option\n",
      "    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):\n",
      "    ...     array[indexer]\n",
      "  return self.array[key]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Calculating cluster c13 DMGs.Calculating cluster c14 DMGs.\n",
      "Calculating cluster c15 DMGs.\n",
      "Calculating cluster c16 DMGs.\n",
      "Calculating cluster c17 DMGs.\n",
      "Calculating cluster c18 DMGs.\n",
      "c17 Finished.\n",
      "Calculating cluster c19 DMGs.\n",
      "c16 Finished.\n",
      "Calculating cluster c2 DMGs.\n",
      "c15 Finished.\n",
      "Calculating cluster c20 DMGs.\n",
      "c14 Finished.\n",
      "Calculating cluster c21 DMGs.\n",
      "c13 Finished.\n",
      "Calculating cluster c22 DMGs.\n",
      "c12 Finished.\n",
      "Calculating cluster c23 DMGs.\n",
      "c11 Finished.\n",
      "Calculating cluster c24 DMGs.\n",
      "c10 Finished.\n",
      "Calculating cluster c25 DMGs.\n",
      "c18 Finished.\n",
      "Calculating cluster c26 DMGs.\n",
      "c21 Finished.\n",
      "Calculating cluster c27 DMGs.\n",
      "c19 Finished.\n",
      "Calculating cluster c28 DMGs.\n",
      "c22 Finished.\n",
      "Calculating cluster c29 DMGs.\n",
      "c20 Finished.\n",
      "Calculating cluster c3 DMGs.\n",
      "c24 Finished.\n",
      "Calculating cluster c30 DMGs.\n",
      "c23 Finished.\n",
      "Calculating cluster c31 DMGs.\n",
      "c25 Finished.\n",
      "Calculating cluster c32 DMGs.\n",
      "Calculating cluster c33 DMGs.\n",
      "c27 Finished.\n",
      "c26 Finished.\n",
      "Calculating cluster c34 DMGs.\n",
      "c29 Finished.\n",
      "Calculating cluster c35 DMGs.\n",
      "c28 Finished.\n",
      "Calculating cluster c36 DMGs.\n",
      "c1 Finished.\n",
      "Calculating cluster c37 DMGs.\n",
      "c30 Finished.\n",
      "Calculating cluster c38 DMGs.\n",
      "c31 Finished.\n",
      "Calculating cluster c39 DMGs.\n",
      "c33 Finished.\n",
      "Calculating cluster c4 DMGs.\n",
      "c34 Finished.\n",
      "Calculating cluster c40 DMGs.\n",
      "c32 Finished.\n",
      "Calculating cluster c5 DMGs.\n",
      "c37 Finished.\n",
      "Calculating cluster c6 DMGs.\n",
      "c35 Finished.\n",
      "Calculating cluster c7 DMGs.\n",
      "c36 Finished.\n",
      "Calculating cluster c8 DMGs.\n",
      "c38 Finished.\n",
      "Calculating cluster c9 DMGs.\n",
      "c39 Finished.\n",
      "c40 Finished.\n",
      "c2 Finished.\n",
      "c3 Finished.\n",
      "c0 Finished.\n",
      "c7 Finished.\n",
      "c9 Finished.\n",
      "c8 Finished.\n",
      "c6 Finished.\n",
      "c5 Finished.\n",
      "c4 Finished.\n"
     ]
    }
   ],
   "source": [
    "dmg_table = one_vs_rest_dmg(cell_meta,\n",
    "                            group=cluster_col,\n",
    "                            mcds_paths=mcds_paths,\n",
    "                            obs_dim=obs_dim,\n",
    "                            var_dim=var_dim,\n",
    "                            mc_type=mc_type,\n",
    "                            top_n=top_n,\n",
    "                            adj_p_cutoff=adj_p_cutoff,\n",
    "                            fc_cutoff=fc_cutoff,\n",
    "                            auroc_cutoff=auroc_cutoff,\n",
    "                            max_cluster_cells=max_cluster_cells,\n",
    "                            max_other_fold=max_other_fold,\n",
    "                            cpu=cpu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-02-16T03:47:55.142943Z",
     "start_time": "2022-02-16T03:47:55.043092Z"
    }
   },
   "outputs": [],
   "source": [
    "dmg_table.to_hdf(f'{cluster_col}.OneVsRestDMG.hdf', key='data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.8.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}