E - Seurat comparison.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "16dca5dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "from statannotations.Annotator import Annotator\n",
    "from statannot import add_stat_annotation\n",
    "from sklearn.metrics import adjusted_mutual_info_score as ANMI\n",
    "from matplotlib.patches import Patch\n",
    "\n",
    "\n",
    "import statannot\n",
    "import json\n",
    "import glob\n",
    "import random\n",
    "import nmi\n",
    "import string\n",
    "\n",
    "from pathlib import Path\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "from helps import *\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "plt.style.use('seaborn')\n",
    "plt.rcParams['font.size'] = 15\n",
    "plt.rcParams[\"xtick.labelsize\"]=13\n",
    "plt.rcParams[\"ytick.labelsize\"]=13\n",
    "plt.rcParams[\"axes.titlesize\"]=15\n",
    "plt.rcParams[\"figure.dpi\"]=600\n",
    "plt.rcParams[\"savefig.format\"]=\"pdf\"\n",
    "plt.rcParams[\"savefig.bbox\"]=\"tight\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7406d0e1",
   "metadata": {},
   "source": [
    "# NMI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c9689b4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3723, 26)\n"
     ]
    },
    {
     "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>orig.ident</th>\n",
       "      <th>nCount_mRNA</th>\n",
       "      <th>nFeature_mRNA</th>\n",
       "      <th>Group</th>\n",
       "      <th>mRNA_snn_res.0.5</th>\n",
       "      <th>seurat_clusters</th>\n",
       "      <th>seurat_clusters_mRNA</th>\n",
       "      <th>lncRNA_snn_res.0.5</th>\n",
       "      <th>seurat_clusters_lncRNA</th>\n",
       "      <th>mRNA.weight</th>\n",
       "      <th>...</th>\n",
       "      <th>hSBM.mRNA_2</th>\n",
       "      <th>hSBM.lncRNA_0</th>\n",
       "      <th>hSBM.lncRNA_1</th>\n",
       "      <th>hSBM.lncRNA_2</th>\n",
       "      <th>hSBM.mRNA.lncRNA_0</th>\n",
       "      <th>hSBM.mRNA.lncRNA_1</th>\n",
       "      <th>hSBM.mRNA.lncRNA_2</th>\n",
       "      <th>triSBM.mRNA.lncRNA_0</th>\n",
       "      <th>triSBM.mRNA.lncRNA_1</th>\n",
       "      <th>triSBM.mRNA.lncRNA_2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>hs.LumR.CGTGAGCTCACATACG.1</th>\n",
       "      <td>SeuratProject</td>\n",
       "      <td>441.838876</td>\n",
       "      <td>351</td>\n",
       "      <td>LumR</td>\n",
       "      <td>7</td>\n",
       "      <td>2</td>\n",
       "      <td>7</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0.739659</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>133</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>50</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>27</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hs.LumR.AACTGGTTCAGCTGGC.1</th>\n",
       "      <td>SeuratProject</td>\n",
       "      <td>375.232822</td>\n",
       "      <td>194</td>\n",
       "      <td>LumR</td>\n",
       "      <td>7</td>\n",
       "      <td>3</td>\n",
       "      <td>7</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>0.585476</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>74</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hs.Lum.CACAGGCCATTAACCG.1</th>\n",
       "      <td>SeuratProject</td>\n",
       "      <td>308.878196</td>\n",
       "      <td>254</td>\n",
       "      <td>Lum</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.528964</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>124</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>13</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>38</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hs.Basal.TTTCCTCAGCCAGTTT.1</th>\n",
       "      <td>SeuratProject</td>\n",
       "      <td>361.110331</td>\n",
       "      <td>295</td>\n",
       "      <td>Basal</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>0.718288</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>39</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>33</td>\n",
       "      <td>10</td>\n",
       "      <td>1</td>\n",
       "      <td>17</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hs.Basal.TCTATTGCAAGAAGAG.1</th>\n",
       "      <td>SeuratProject</td>\n",
       "      <td>328.966235</td>\n",
       "      <td>257</td>\n",
       "      <td>Basal</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.828428</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>39</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>16</td>\n",
       "      <td>10</td>\n",
       "      <td>1</td>\n",
       "      <td>17</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                orig.ident  nCount_mRNA  nFeature_mRNA  Group  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1   SeuratProject   441.838876            351   LumR   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1   SeuratProject   375.232822            194   LumR   \n",
       "hs.Lum.CACAGGCCATTAACCG.1    SeuratProject   308.878196            254    Lum   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1  SeuratProject   361.110331            295  Basal   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1  SeuratProject   328.966235            257  Basal   \n",
       "\n",
       "                             mRNA_snn_res.0.5  seurat_clusters  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1                  7                2   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1                  7                3   \n",
       "hs.Lum.CACAGGCCATTAACCG.1                   4                1   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1                 2                4   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1                 2                0   \n",
       "\n",
       "                             seurat_clusters_mRNA  lncRNA_snn_res.0.5  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1                      7                   2   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1                      7                   3   \n",
       "hs.Lum.CACAGGCCATTAACCG.1                       4                   1   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1                     2                   4   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1                     2                   0   \n",
       "\n",
       "                             seurat_clusters_lncRNA  mRNA.weight  ...  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1                        2     0.739659  ...   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1                        3     0.585476  ...   \n",
       "hs.Lum.CACAGGCCATTAACCG.1                         1     0.528964  ...   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1                       4     0.718288  ...   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1                       0     0.828428  ...   \n",
       "\n",
       "                             hSBM.mRNA_2  hSBM.lncRNA_0  hSBM.lncRNA_1  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1             0            133              1   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1             0             74              1   \n",
       "hs.Lum.CACAGGCCATTAACCG.1              0            124              1   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1            1             39              5   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1            1             39              5   \n",
       "\n",
       "                            hSBM.lncRNA_2  hSBM.mRNA.lncRNA_0  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1              0                  50   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1              0                   2   \n",
       "hs.Lum.CACAGGCCATTAACCG.1               0                  13   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1             1                  33   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1             1                  16   \n",
       "\n",
       "                             hSBM.mRNA.lncRNA_1  hSBM.mRNA.lncRNA_2  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1                    8                   0   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1                    2                   0   \n",
       "hs.Lum.CACAGGCCATTAACCG.1                     4                   0   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1                  10                   1   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1                  10                   1   \n",
       "\n",
       "                             triSBM.mRNA.lncRNA_0  triSBM.mRNA.lncRNA_1  \\\n",
       "hs.LumR.CGTGAGCTCACATACG.1                     27                     0   \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1                      7                     0   \n",
       "hs.Lum.CACAGGCCATTAACCG.1                      38                     2   \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1                    17                     3   \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1                    17                     3   \n",
       "\n",
       "                             triSBM.mRNA.lncRNA_2  \n",
       "hs.LumR.CGTGAGCTCACATACG.1                      1  \n",
       "hs.LumR.AACTGGTTCAGCTGGC.1                      1  \n",
       "hs.Lum.CACAGGCCATTAACCG.1                       1  \n",
       "hs.Basal.TTTCCTCAGCCAGTTT.1                     0  \n",
       "hs.Basal.TCTATTGCAAGAAGAG.1                     0  \n",
       "\n",
       "[5 rows x 26 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df=pd.read_csv(\"Results/Seurat/seurat_hsbm_metadata.txt\", sep=\"\\t\")\n",
    "print(df.shape)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9542f400",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d664041b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NMI seurat_clusters_wnn vs seurat_clusters_mRNA --> 0.855\n",
      "ANMI seurat_clusters_wnn vs seurat_clusters_mRNA --> 0.855\n",
      "\n",
      "\n",
      "NMI seurat_clusters_wnn vs seurat_clusters_lncRNA --> 0.57\n",
      "ANMI seurat_clusters_wnn vs seurat_clusters_lncRNA --> 0.569\n",
      "\n",
      "\n",
      "NMI triSBM.mRNA.lncRNA_1 vs hSBM.mRNA.lncRNA_1 --> 0.494\n",
      "ANMI triSBM.mRNA.lncRNA_1 vs hSBM.mRNA.lncRNA_1 --> 0.492\n",
      "\n",
      "\n",
      "NMI triSBM.mRNA.lncRNA_1 vs hSBM.lncRNA_1 --> 0.433\n",
      "ANMI triSBM.mRNA.lncRNA_1 vs hSBM.lncRNA_1 --> 0.431\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "couples = [(\"seurat_clusters_wnn\",\"seurat_clusters_mRNA\"), \n",
    "           (\"seurat_clusters_wnn\",\"seurat_clusters_lncRNA\"),\n",
    "           (\"triSBM.mRNA.lncRNA_1\",\"hSBM.mRNA.lncRNA_1\"),\n",
    "           (\"triSBM.mRNA.lncRNA_1\",\"hSBM.lncRNA_1\")]\n",
    "for coup in couples:\n",
    "    NMI=np.around(nmi.compute_normalised_mutual_information(df[coup[0]], df[coup[1]]), decimals=3)\n",
    "    print(\"NMI\", coup[0], \"vs\", coup[1], \"-->\", NMI)\n",
    "    print(\"ANMI\", coup[0], \"vs\", coup[1], \"-->\", round(ANMI(df[coup[0]],  df[coup[1]]), 3))\n",
    "    print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "872b080e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seurat_clusters_mRNA vs Group NMI: 0.784    NMI/NMI* 464.48469547160386\n",
      "seurat_clusters_lncRNA vs Group NMI: 0.672    NMI/NMI* 614.1855921717237\n",
      "seurat_clusters_wnn vs Group NMI: 0.78    NMI/NMI* 450.95479195259486\n"
     ]
    }
   ],
   "source": [
    "couples=[(\"seurat_clusters_mRNA\",\"Group\"),\n",
    "           (\"seurat_clusters_lncRNA\",\"Group\"),\n",
    "           (\"seurat_clusters_wnn\",\"Group\")]\n",
    "for coup in couples:\n",
    "    NMI=np.around(nmi.compute_normalised_mutual_information(df[coup[0]], df[coup[1]]),decimals=3)\n",
    "    nmi_rand=0\n",
    "    for k in range(1000):\n",
    "        a=df[coup[0]].to_list()\n",
    "        np.random.shuffle(a)\n",
    "        nmi_rand+=nmi.compute_normalised_mutual_information(df[coup[1]],a)/1000\n",
    "    print(coup[0], \"vs\", coup[1], \"NMI:\", NMI, \"   NMI/NMI*\", NMI/nmi_rand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f68ece3c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8976faad",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.replace(\"LumR\", \"Lum-R\", inplace=True)\n",
    "df.replace(\"BasalR\", \"Basal-R\", inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9253f711",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=pd.DataFrame(index=df.index)\n",
    "data[\"mRNA weight\"]=df[\"mRNA.weight\"]\n",
    "data[\"lncRNA weight\"]=df[\"lncRNA.weight\"]\n",
    "data=data.melt()\n",
    "data.columns=[\"Gene family\",\"WNN weight\"]\n",
    "data.index=list(df.index)+list(df.index)\n",
    "data[\"typehisto\"]=[df.loc[ind][\"Group\"] for ind in data.index]\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dbbd5ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "subtypes=list(set(data[\"typehisto\"]))\n",
    "fig, axs = plt.subplots(1,2, figsize=(14,6), dpi=600)\n",
    "\n",
    "sns.boxplot(data=data, x=\"Gene family\",y=\"WNN weight\", ax=axs[0],\n",
    "            palette=[\"green\",\"violet\"], saturation=0.65)\n",
    "test_results = add_stat_annotation(data=data, x=\"Gene family\",y=\"WNN weight\", ax=axs[0],\n",
    "                                   box_pairs=[(\"mRNA weight\",\"lncRNA weight\")],\n",
    "                                   test='Mann-Whitney', text_format='simple',\n",
    "                                   comparisons_correction=None,\n",
    "                                   text_annot_custom=[\"Mann-Whitney U-Test p-value < 1$^{-10}$\"],\n",
    "                                   loc='inside', verbose=2)\n",
    "\n",
    "\n",
    "sns.boxplot(data=data, x=\"Gene family\",y=\"WNN weight\", hue=\"typehisto\", ax=axs[1],\n",
    "           palette=dict(zip(subtypes, nmi.set_colors(subtypes))),  saturation=0.65)\n",
    "pairs=[((\"mRNA weight\", \"Basal\"),   (\"lncRNA weight\", \"Basal\")),\n",
    "       ((\"mRNA weight\", \"Basal-R\"), (\"lncRNA weight\", \"Basal-R\")),\n",
    "       ((\"mRNA weight\", \"Lum\"),   (\"lncRNA weight\", \"Lum\")),\n",
    "       ((\"mRNA weight\", \"Lum-R\"),   (\"lncRNA weight\", \"Lum-R\"))]\n",
    "test_results = add_stat_annotation(data=data, x=\"Gene family\",y=\"WNN weight\",\n",
    "                                   hue=\"typehisto\", ax=axs[1], box_pairs=pairs,\n",
    "                                   test='Mann-Whitney', text_format='simple',\n",
    "                                   comparisons_correction=None,\n",
    "                                   text_annot_custom=[\"MW U p < 1$^{-10}$\" for i in range(4)],\n",
    "                                   loc='inside', verbose=1)\n",
    "\n",
    "\n",
    "axs[0].set_ylabel(\"WNN weight\", weight=\"bold\", size=20)\n",
    "axs[1].set_ylabel(\"WNN weight\", weight=\"bold\", size=20)\n",
    "\n",
    "axs[0].set_xlabel(\"Gene family\", weight=\"bold\", size=20)\n",
    "axs[1].set_xlabel(\"Gene family\", weight=\"bold\", size=20)\n",
    "\n",
    "axs[0].set_xticks([0,0.25,0.5,0.75,1])\n",
    "axs[1].set_yticks([0,0.25,0.5,0.75,1])\n",
    "\n",
    "title=f\"Fig_7\"\n",
    "plt.savefig(f\"Results/Figures/{title}.png\", dpi=600)\n",
    "#plt.savefig(f\"Results/Figures/{title}.pdf\", dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e915591f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "485098e1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e9be05f",
   "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.10.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "197px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}