02_bio_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3485b63b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "import gseapy as gp\n",
    "import matplotlib_venn as mvenn\n",
    "\n",
    "import utils_description as utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d97a246",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(style=\"white\", font='Arial', font_scale=1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc0222ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "infile_edges    = 'results/edges_scores.csv'\n",
    "infile_clusters = 'results/results_clustering.csv'\n",
    "infile_metadata = 'data/ADNI/ADNIMERGE_processed.csv'\n",
    "infile_network  = 'data/networks/PPI_SNAP_brain_False.edgelist'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d911a3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "edges_scores = pd.read_csv(infile_edges, index_col=0).T\n",
    "clusters     = pd.read_csv(infile_clusters, index_col=0)['cluster_3'].rename('cluster')\n",
    "metadata     = pd.read_csv(infile_metadata, index_col=0)\n",
    "network      = nx.read_edgelist(infile_network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12c6a6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.concat([metadata, clusters], axis=1, join='inner')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd6165e2",
   "metadata": {},
   "source": [
    "# Representative graphs of each cluster\n",
    "Features of each representative graph defining each cluster.\n",
    "\n",
    "### Obtain representative graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "196412da",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "utils.obtain_cluster_graphs_edges_scores(edges_scores, clusters, network)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83a6b263",
   "metadata": {},
   "source": [
    "### Compare the similarity between graphs\n",
    "Against the original network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8aca8760",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for k in range(len(data['cluster'].unique())):\n",
    "    \n",
    "    net_cluster_file = f'results/cluster{k}.edgelist'\n",
    "    cluster_network  = nx.read_edgelist(net_cluster_file)\n",
    "    \n",
    "    sim   = utils.jaccard_similarity(network, cluster_network)\n",
    "    sim_w = utils.weighted_jaccard(network, cluster_network)\n",
    "    \n",
    "    print(f'Cluster {k} vs. Original: ', sim, sim_w)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97cae874",
   "metadata": {},
   "source": [
    "Between clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cedfe111",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(data['cluster'].unique())):\n",
    "    \n",
    "    net_cluster_file_i = f'results/cluster{i}.edgelist'\n",
    "    cluster_network_i  = nx.read_edgelist(net_cluster_file_i)\n",
    "    \n",
    "    for j in range(i, len(data['cluster'].unique())):\n",
    "        \n",
    "        if i != j:\n",
    "    \n",
    "            net_cluster_file_j = f'results/cluster{j}.edgelist'\n",
    "            cluster_network_j  = nx.read_edgelist(net_cluster_file_j)\n",
    "\n",
    "            sim   = utils.jaccard_similarity(cluster_network_i, cluster_network_j)\n",
    "            sim_w = utils.weighted_jaccard(cluster_network_i, cluster_network_j)\n",
    "            \n",
    "            print(f'Cluster {j} vs. Cluster {i}: ', sim, sim_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54aa3659",
   "metadata": {},
   "source": [
    "### Global graphs metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3c615b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add unit weights to original network\n",
    "nx.set_edge_attributes(network, values = 1, name = 'weight')\n",
    "print('Original network')\n",
    "utils.graph_metrics(network)\n",
    "print()\n",
    "\n",
    "for k in [0, 1, 2]:\n",
    "\n",
    "    net_cluster_file = f'results/cluster{k}.edgelist'\n",
    "    cluster_network  = nx.read_edgelist(net_cluster_file)\n",
    "    print(f'Cluster {k} network')\n",
    "    utils.graph_metrics(cluster_network)\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9546aaf5",
   "metadata": {},
   "source": [
    "# Significantly different edge scores between clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1965fbf8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "edges_scores = edges_scores.loc[data.index]\n",
    "statistics_edges = utils.get_significant_edges(edges_scores, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "903ca562",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(statistics_edges[['c0_mean', 'c1_mean', 'c2_mean', 'F', 'pvalue', 'significant', 'tukey']].to_latex())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "378a395d",
   "metadata": {},
   "source": [
    "### Visualization of the significant edge scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "323e98bd",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "columns = statistics_edges.index.to_list()\n",
    "columns.append('cluster')\n",
    "\n",
    "heatmap_data = pd.concat([edges_scores, clusters], axis=1, join='inner')\n",
    "heatmap_data = heatmap_data[columns]\n",
    "heatmap_data = heatmap_data.sort_values(by='cluster')\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "\n",
    "# Edges\n",
    "data1 = heatmap_data.copy()\n",
    "data1['cluster'] = float('nan')\n",
    "ax = sns.heatmap(data1.T, cmap=\"rocket\")\n",
    "\n",
    "# Clusters\n",
    "data2 = heatmap_data.copy()\n",
    "data2[statistics_edges.index] = float('nan')\n",
    "print(data2['cluster'].value_counts())\n",
    "my_cmap = sns.hls_palette(n_colors=3)\n",
    "sns.heatmap(data2.T, cmap=my_cmap, cbar=False)\n",
    "\n",
    "plt.tick_params(axis='x', which='both', bottom=False, top=False, labelbottom=False)\n",
    "plt.tight_layout()\n",
    "# plt.savefig('figures/heatmap_edges.png', dpi=500)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c84a724",
   "metadata": {},
   "source": [
    "### Pathway Enrichment Analysis\n",
    "To show biological differences between clusters. Make a list of nodes (genes) in edges that (i) are significantly different between two clusters, and (ii) have the lowest mean. That is, edges that are significantly \"worse\" (more affected) between clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb2c1448",
   "metadata": {},
   "outputs": [],
   "source": [
    "c0_set, c1_set, c2_set, original_set = utils.make_gene_list(statistics_edges, network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fcc2178",
   "metadata": {},
   "outputs": [],
   "source": [
    "source = 'Reactome_2022' # Reactome_2022\n",
    "pea = utils.enrichment_analysis([source], c0_set, c1_set, c2_set, original_set, statistics_edges.index)\n",
    "pea['Adjusted P-value'] = pea['Adjusted P-value'].round(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35631c1d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# print(pea[['Term', 'Adjusted P-value', 'Genes', 'Edges', 'cluster']].to_latex(index=None))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aa009c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "c0_terms = set(pea.loc[pea['cluster'] == 'Cluster 0']['Term'].values.tolist())\n",
    "c1_terms = set(pea.loc[pea['cluster'] == 'Cluster 1']['Term'].values.tolist())\n",
    "c2_terms = set(pea.loc[pea['cluster'] == 'Cluster 2']['Term'].values.tolist())\n",
    "\n",
    "plt.figure(figsize=(15, 15))\n",
    "venn = mvenn.venn3([c0_terms, c1_terms, c2_terms], ('Cluster 1', 'Cluster 2', 'Cluster 3'))\n",
    "\n",
    "venn.get_label_by_id('100').set_text('\\n'.join(c0_terms-c1_terms-c2_terms))\n",
    "venn.get_label_by_id('010').set_text('\\n'.join(c1_terms-c2_terms-c0_terms))\n",
    "venn.get_label_by_id('001').set_text('\\n'.join(c2_terms-c1_terms-c0_terms))\n",
    "\n",
    "venn.get_label_by_id('101').set_text('\\n'.join(c0_terms&c2_terms-c1_terms))\n",
    "venn.get_label_by_id('111').set_text('\\n'.join(c0_terms&c2_terms&c1_terms))\n",
    "\n",
    "plt.tight_layout()\n",
    "# plt.savefig('figures/venn3_reactome.png', dpi=500)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27dc24df",
   "metadata": {},
   "source": [
    "### Variants analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7976b7ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "variants = pd.read_csv('results/processed_variants_ADNI_WGS.csv', index_col=0)\n",
    "\n",
    "nodes_file = open(f'results/nodes_significant.txt', 'r')\n",
    "genes = nodes_file.read().split('\\n')\n",
    "nodes_file.close()\n",
    "\n",
    "variants_sel = variants.loc[variants['SYMBOL'].isin(genes)]\n",
    "variants_sel = variants_sel.sort_values(by='SYMBOL', ascending=True)\n",
    "variants_sel.drop(columns='SYMBOL', inplace=True)\n",
    "variants_sel = variants_sel.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89b39b15",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.metrics as metrics\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split, StratifiedKFold\n",
    "from imblearn.ensemble import BalancedRandomForestClassifier\n",
    "import numpy as np\n",
    "# import shap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28e72d78",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "data_variants = variants.drop(columns='SYMBOL')\n",
    "data_variants = data_variants.drop_duplicates(keep=False)\n",
    "# dataset = pd.concat([data_variants.T, clusters], axis=1, join='inner')\n",
    "dataset = pd.concat([variants_sel, clusters], axis=1, join='inner')\n",
    "\n",
    "y = dataset['cluster']\n",
    "x = dataset.drop(columns=['cluster'])\n",
    "\n",
    "for task in ['C1_vs_All', 'C2_vs_All', 'C3_vs_All']:\n",
    "    \n",
    "    if task == 'C1_vs_All':\n",
    "        y_tmp = y.replace({0:'Positive', 1:'Negative', 2:'Negative'})\n",
    "        title = '(a) Cluster 1 vs. All\\n'\n",
    "        lim = 0.15\n",
    "        \n",
    "    elif task == 'C2_vs_All':\n",
    "        y_tmp = y.replace({0:'Negative', 1:'Positive', 2:'Negative'})\n",
    "        title = '(b) Cluster 2 vs. All\\n'\n",
    "        lim = 0.10\n",
    "        \n",
    "    elif task == 'C3_vs_All':\n",
    "        y_tmp = y.replace({0:'Negative', 1:'Negative', 2:'Positive'})\n",
    "        title = '(c) Cluster 3 vs. All\\n'\n",
    "        lim = 0.26\n",
    "\n",
    "    x_train, x_test, y_train, y_test = train_test_split(x, y_tmp, stratify=y_tmp, test_size=0.2, random_state=42)\n",
    "\n",
    "    model = BalancedRandomForestClassifier(random_state=42)\n",
    "    model.fit(x_train, y_train)\n",
    "    y_pred = model.predict(x_test)\n",
    "    \n",
    "    print(task)\n",
    "    print()\n",
    "\n",
    "    cm = pd.DataFrame(metrics.confusion_matrix(y_test, y_pred))\n",
    "    plt.figure(figsize=(3.5, 2.5))\n",
    "    ax = sns.heatmap(cm, annot=True, cmap='crest', cbar=False, fmt='g')\n",
    "    ax.tick_params(tick2On=False, labelsize=False)\n",
    "    # plt.suptitle(title)\n",
    "    plt.title('Confusion matrix')\n",
    "    plt.xlabel('True')\n",
    "    plt.ylabel('Predicted')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(f'figures/conf_matrix_{task}.png', dpi=300)\n",
    "\n",
    "    data_symbols = variants.groupby(variants.index)['SYMBOL'].apply(lambda x: ',\\n'.join(x.astype(str)))\n",
    "\n",
    "    importances = model.feature_importances_\n",
    "    forest_importances = pd.DataFrame(importances, index=x.columns)\n",
    "    forest_importances = pd.concat([forest_importances, data_symbols], axis=1, join='inner')\n",
    "    forest_importances.reset_index(inplace=True)\n",
    "    forest_importances['variant'] = forest_importances['index'] + ' (' + forest_importances['SYMBOL'] + ')'\n",
    "\n",
    "    plt.figure(figsize=(4, 6))\n",
    "    sns.barplot(forest_importances.sort_values(by=0, ascending=False).head(10), y='variant', x=0, palette='viridis')\n",
    "    plt.title('Feature importances')\n",
    "    plt.ylabel('')\n",
    "    plt.xlim(0, lim)\n",
    "    plt.tight_layout()\n",
    "    # plt.savefig(f'figures/importances_{task}.png', dpi=300)\n",
    "\n",
    "    for v in forest_importances.sort_values(by=0, ascending=False).head(10)['index']:\n",
    "        print(v)\n",
    "    print()"
   ]
  }
 ],
 "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}