B1 - hSBM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbb7d4d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import graph_tool.all as gt\n",
    "\n",
    "import json\n",
    "import random\n",
    "import nmi\n",
    "import glob\n",
    "import string\n",
    "\n",
    "from pathlib import Path\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "from time import localtime, strftime\n",
    "from sbmtm import sbmtm\n",
    "from nsbm import nsbm\n",
    "\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": "66eeac1f",
   "metadata": {},
   "source": [
    "# hSBM: experiment run\n",
    "\n",
    "Here we show how to reproduce all the experiments with hSBM. To select the desidered experiment change the value of variable dataset choosing one of the three following values: hSBM-mRNA, hSBM-lncRNA, hSBM-mRNA-lncRNA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "168add40",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset=\"hSBM-mRNA\"\n",
    "print(dataset)\n",
    "df=pd.read_csv(f\"Results/{dataset}/{dataset}.csv.gz\",index_col=0)\n",
    "df=df.sample(frac=1, axis=1)\n",
    "df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d9e4d8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "model=sbmtm()\n",
    "model.make_graph_from_BoW_df(df)\n",
    "V=model.get_V()\n",
    "D=model.get_D()\n",
    "print(V/(D*(D-1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb82b105",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3 minutes for the hSBM-mRNA experiment\n",
    "# 1 minute for the hSBM-lncRNA experiment\n",
    "# 3 minutes for the hSBM-mRNA-lncRNA experiment\n",
    "#with a i5-8265U 4 cores 1.60 GHz laptop\n",
    "\n",
    "print(strftime(\"%Y-%m-%d %H:%M:%S\", localtime()))\n",
    "model.fit(n_init=1)\n",
    "print(strftime(\"%Y-%m-%d %H:%M:%S\", localtime()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd30e619",
   "metadata": {},
   "outputs": [],
   "source": [
    "Path(f\"Results/{dataset}/hSBM\").mkdir(parents=True, exist_ok=True)\n",
    "model.save_graph(filename=f\"Results/{dataset}/hSBM/{dataset}-graph.xml.gz\")\n",
    "model.save_levels(f\"Results/{dataset}/hSBM/{dataset}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa045e56",
   "metadata": {},
   "source": [
    "# hSBM: experiment analysis\n",
    "Here you can find all the steps we followed to analyse the outcome of all the hSBM's experiments.\n",
    "You can reproduce all the experiments changing the variable \"dataset\" that corresponds to the experiment's name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d366ef5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset=\"hSBM-mRNA\"\n",
    "print(dataset)\n",
    "info=pd.read_csv(\"HelperFiles/ENS-Info.txt\",index_col=0,sep=\"\\t\")\n",
    "df=pd.read_csv(f\"Results/{dataset}/{dataset}.csv.gz\",index_col=0)\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "024d9ced",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "labels=pd.read_csv(\"HelperFiles/All-datasets-labels.csv\",index_col=0)\n",
    "labels=labels.loc[df.columns]\n",
    "labels.index=labels.index.astype(str)\n",
    "subtypes=list(sorted(set(labels.typehisto)))\n",
    "print(df.shape)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2963b8c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "performances={}\n",
    "for level in [0,1,2]:\n",
    "    with open(f\"Results/{dataset}/hSBM/{dataset}-cluster-level-{level}.txt\") as f:\n",
    "        clusters=json.load(f)\n",
    "    cluster_df=pd.DataFrame.from_dict(clusters,orient=\"index\")\n",
    "    labels[\"hSBM\"]=\"--\"\n",
    "    for i in range(len(clusters)):\n",
    "        labels[\"hSBM\"].loc[np.asarray(np.asarray(clusters[str(i)])[:,0])]=i\n",
    "    \n",
    "    NMI=np.around(nmi.compute_normalised_mutual_information(labels.typehisto,labels[\"hSBM\"]),decimals=3)\n",
    "    nmi_rand=0\n",
    "    for k in range(1000):\n",
    "        a=labels[\"hSBM\"].to_list()\n",
    "        np.random.shuffle(a)\n",
    "        nmi_rand+=nmi.compute_normalised_mutual_information(labels[\"typehisto\"],a)/1000\n",
    "\n",
    "    performances[f\"Level {level}\"]=[NMI,NMI/nmi_rand]\n",
    "with open(f\"Results/{dataset}/hSBM/{dataset}_NMI.json\", 'w') as fp:\n",
    "    json.dump(performances, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "787355f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(f\"Results/{dataset}/hSBM/{dataset}_NMI.json\") as f:\n",
    "        performances=json.load(f)\n",
    "for key in performances.keys():\n",
    "    print(f\"{key} NMI: {performances[key][0]}, NMI/NMI*: {int(performances[key][1])}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5370b6d8",
   "metadata": {},
   "source": [
    "# NMI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf4d0df5",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(f\"Results/{dataset}/hSBM/{dataset}_NMI.json\") as f:\n",
    "        perfs=json.load(f)\n",
    "perfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4605db3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,3,figsize=(25,8), dpi=600)\n",
    "axs=axs.flatten()\n",
    "ax=0\n",
    "for level in [0,1,2]:\n",
    "    with open(f\"Results/{dataset}/hSBM/{dataset}-cluster-level-{level}.txt\") as f:\n",
    "        clusters=json.load(f)\n",
    "    cluster_df=pd.DataFrame.from_dict(clusters,orient=\"index\")\n",
    "    labels[\"hSBM\"]=\"--\"\n",
    "    for i in range(len(clusters)):\n",
    "        labels[\"hSBM\"].loc[np.asarray(np.asarray(clusters[str(i)])[:,0])]=i\n",
    "     \n",
    "    labels[\"typehisto_1\"]=pd.Series(list(labels[\"typehisto\"])).astype('category').cat.codes.values    \n",
    "    fraction_sites = pd.DataFrame(index=labels[\"hSBM\"].unique(), columns=sorted(labels[\"typehisto_1\"].unique())[::-1]).fillna(0)\n",
    "    for sample in labels[[\"hSBM\",\"typehisto_1\"]].values:\n",
    "        fraction_sites.at[sample[0],sample[1]] += 1\n",
    "\n",
    "    fraction_sites = fraction_sites.sort_values(by=list(fraction_sites.columns), ascending=True)\n",
    "    fraction_sites.columns=subtypes[::-1]\n",
    "    fraction_sites.index=[i for i in range(len(fraction_sites))]\n",
    "    fraction_sites.plot.bar(stacked=True, color=dict(zip(subtypes, nmi.set_colors(subtypes))),\n",
    "                           width=1, alpha=0.65, ax=axs[level])    \n",
    "       \n",
    "    axs[level].set_xlabel(f\"\\nClustering level {level}\", size=25, weight='bold')\n",
    "    axs[level].set_ylabel(\"Number of cells\", size=25, weight='bold')\n",
    "    axs[level].yaxis.set_major_formatter(FormatStrFormatter('%.0f'))\n",
    "    axs[level].tick_params(axis='both', which='major', labelsize=25, rotation=0)\n",
    "    \n",
    "    legend_properties = {'weight':'bold', \"size\":\"x-large\"}\n",
    "    if level==1:\n",
    "        axs[level].legend(loc=(0.55,0.65), prop=legend_properties)\n",
    "    else:\n",
    "        axs[level].get_legend().remove()        \n",
    "    \n",
    "    axs[level].text(-0.055, 1.1, string.ascii_uppercase[level],\n",
    "                 transform=axs[ax].transAxes, size=35, weight='bold',rotation=0)\n",
    "    axs[level].xaxis.set_major_locator(plt.MaxNLocator(min(12, len(set(labels.hSBM))+1)))\n",
    "    ax+=1   \n",
    "    \n",
    "fig.tight_layout(pad=1)    \n",
    "title=f\"Fig_1\"\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": "4e0c475d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "16a677a3",
   "metadata": {},
   "source": [
    "# Data reading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f897ab44",
   "metadata": {},
   "outputs": [],
   "source": [
    "level=1\n",
    "data=hSBM_data(dataset=dataset,labels=labels, lev=level, info=info)\n",
    "print(data.keys())\n",
    "\n",
    "labels=data[\"sample_cluster\"]\n",
    "topic_gene_prob=data[\"topic_gene_prob\"]\n",
    "topic_gene=data[\"topic_gene\"]\n",
    "topic_gene_genename=data[\"topic_gene_genename\"]\n",
    "topic_gene_raw=data[\"topic_gene_raw\"]\n",
    "p_topic_sample_cent=data[\"p_topic_sample_cent\"]\n",
    "p_topic_sample=data[\"p_topic_sample\"]\n",
    "\n",
    "\n",
    "Path(f\"Results/{dataset}/hSBM/Data\").mkdir(parents=True, exist_ok=True)\n",
    "for key in data.keys():\n",
    "    data[key].to_csv(f\"Results/{dataset}/hSBM/Data/{key}-level-{level}.csv\")\n",
    "\n",
    "\n",
    "print(\"N clusters: \", len(set(labels.hSBM)))\n",
    "print(\"N topics:\", p_topic_sample_cent.shape[1])\n",
    "\n",
    "clusters=sorted(list(set(labels.hSBM)))\n",
    "\n",
    "p_topic_class_cent=pd.DataFrame(index=clusters, columns=p_topic_sample.columns)    \n",
    "for cla in clusters:\n",
    "    in_sample=labels[labels[\"hSBM\"]==cla].index\n",
    "    p_topic_class_cent.loc[cla]=p_topic_sample_cent.loc[in_sample].mean() \n",
    "p_topic_class_cent.index.names=[\"clusters\"]\n",
    "data[\"p_topic_class_cent\"]=p_topic_class_cent\n",
    "\n",
    "path_to_save=f\"Results/{dataset}/hSBM/Outcome analysis/Level {level}/\"\n",
    "Path(f\"{path_to_save}\").mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01ca7044",
   "metadata": {},
   "source": [
    "# Topic - cluster assignment\n",
    "Here we assign to each cluster the most appropriate topics following the method explained in the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edb8369b",
   "metadata": {},
   "outputs": [],
   "source": [
    "topic_arr, threshold = loop_topics(clusters, p_topic_class_cent, direction=\"up\")\n",
    "with open(f\"{path_to_save}/{dataset} hSBM level {level} Topic up clusters threshold {np.around(threshold, decimals=2)}.json\", 'w') as fp:\n",
    "    json.dump(topic_arr, fp)\n",
    "topic_arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "376809ba",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "assigned_topic=flat_list([topic_arr[k] for k in topic_arr.keys()])\n",
    "assigned_topic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ba48996",
   "metadata": {},
   "source": [
    "## Enrichement test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "468a3a09",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Here you need to choose the appropriate file between \"MSigDB.json\" for the hSBM-mRNA experiment \n",
    "#and \"lncSEA_red.json\" for the hSBM-lncRNA experiment. \n",
    "\n",
    "path=f\"Results/{dataset}/hSBM/Outcome analysis/Level {level}\"\n",
    "Path(f\"{path}/Enrichment test topics\").mkdir(parents=True, exist_ok=True)\n",
    "with open('HelperFiles/MSigDB.json') as f:\n",
    "    all_lists=json.load(f)\n",
    "print(len(all_lists.keys()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34d58c25",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Please note that this is the slowest step of the procedure: it takes about 20-30 seconds to test \n",
    "#each topic with a i5-8265U 4 cores 1.60 GHZ laptop\n",
    "enrichment_test(all_lists, topic_gene_genename, topic_gene_raw, info, dataset[5:],\n",
    "               f\"{path}/Enrichment test topics\", dataset, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26588c83",
   "metadata": {},
   "outputs": [],
   "source": [
    "path=f\"Results/{dataset}/hSBM/Outcome analysis/Level {level}\"\n",
    "dfs_hgt={}\n",
    "for topic in topic_gene.columns:\n",
    "    dfs_hgt[topic]=pd.read_csv(f\"{path}/Enrichment test topics/{dataset} level {level} Enrichment Test {dataset[5:]}-topic {topic}.csv\",index_col=0)\n",
    "    dfs_hgt[topic].drop(dfs_hgt[topic][dfs_hgt[topic].fdr<3].dropna().index,inplace=True)\n",
    "    dfs_hgt[topic].sort_values(by=\"fdr\",inplace=True,ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d412cfec",
   "metadata": {},
   "source": [
    "Before calling the function topics_names check if in the first 10-20 positions of each dataframe in the dictionary dfs_hgt there are gene sets that totally don't matter with your data.\n",
    "After that topics_names has run, you should check the assigment makes sense, otherwise you can manually change it. \n",
    "\n",
    "Please not that the developed procedure doesn't take in account the biological meaning of gene sets in the database, it's just a fully reproducibile method that assign the best \"name\" to each topic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "409171cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "topic_name=topics_names(topics=assigned_topic,enr_test_outcome=dfs_hgt,\n",
    "                        database=all_lists)\n",
    "topic_name.sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11ddc07a",
   "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": "231px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}