genotype-graph.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graph of SNPs in the genome\n",
    "The idea here is to build a similarity graph between SNPs. In the raw data, there are many SNPs that are identical, only their names change (see data exploration notebook). Eliminating the duplicates, or grouping them together as a single feature, reduces the analysis time and complexity. Going one step further in that direction, it may be interesting to group together the SNPs that are highly similar as well, while still keeping them distinguishable. This is the idea developed here with a graph approach.\n",
    "\n",
    "In this notebook, SNPs are nodes of the graph. They are connected to their k nearest neighbors. The connections are weighted by the similarity of the SNPs according to a chosen distance. To each SNP is associated a vector encoding its variations over the BXD mouse dataset. Two SNPs are similar if their vectors are close in term of Euclidean distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.sparse import csr_matrix\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import sklearn.metrics\n",
    "import sklearn.neighbors\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Config for accessing the data on the s3 storage\n",
    "storage_options = {'anon':True, 'client_kwargs':{'endpoint_url':'https://os.unil.cloud.switch.ch'}}\n",
    "s3_path = 's3://lts2-graphnex/BXDmice/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data\n",
    "genotype_path = os.path.join(s3_path, 'geno_reduced.csv.gz')\n",
    "#genotype_path = os.path.join(s3_path, 'genotype_BXD.csv.gz')\n",
    "genotype = pd.read_csv(genotype_path, storage_options=storage_options)\n",
    "print('File {} Opened.'.format(genotype_path))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing the distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the data as a numpy array\n",
    "geno_values = genotype.loc[:,'B6D2F1':].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Default distance is Euclidean\n",
    "num_neighbors = 4\n",
    "geno_knn = sklearn.neighbors.kneighbors_graph(geno_values, num_neighbors, mode='distance')\n",
    "# Optionally, one can use the following function to compute all the distances:\n",
    "#geno_distances = sklearn.metrics.pairwise_distances(geno_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Distribution of weights\n",
    "plt.hist(geno_knn.data, bins=20)\n",
    "plt.title('Distribution of distances')\n",
    "plt.xlabel('Distance')\n",
    "plt.ylabel('Nb of edges')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Distance to weight\n",
    "# Modify the non-zero values to turn them into weights instead of distances\n",
    "def distance2weight(d):\n",
    "    sigma = 1\n",
    "    return np.exp(- sigma * d)\n",
    "    \n",
    "M = geno_knn.copy()\n",
    "M.data = distance2weight(geno_knn.data)\n",
    "\n",
    "print('A distance of 1 becomes a weight of {}.'.format(str(distance2weight(1))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Distribution of weights\n",
    "plt.hist(M.data, bins=20)\n",
    "plt.title('Distribution of weights')\n",
    "plt.xlabel('Weight value')\n",
    "plt.ylabel('Nb of edges')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.from_scipy_sparse_matrix(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adding info on the nodes of the graph\n",
    "genoinfo_dic = genotype[['SNP','Chr','Pos']].to_dict()\n",
    "nx.set_node_attributes(G,genoinfo_dic['SNP'],name='SNP') # SNP id\n",
    "nx.set_node_attributes(G,genoinfo_dic['Chr'],name='Chr') # Chromosome\n",
    "nx.set_node_attributes(G,genoinfo_dic['Pos'],name='Pos') # position inside the chromosome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Saving the graph as a gexf file readable with Gephi.\n",
    "nx.write_gexf(G,'SNPgraph.gexf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graph plotted using Gephi, colored by chromosome. SNPs share their similarity mostly with their spatial neighbors inside chromosomes. Most of the subgraph have a unique color (SNPs from the same chromosome).\n",
    "\n",
    "![SNP graph](SNPgraphChr.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications of the graph\n",
    "There are different possible applications of this graph\n",
    "\n",
    "* SNP reduction: reduce the number of SNPs to analyse by grouping similar ones together or keeping only some representatives\n",
    "    * subsample the SNPs regularly over the graph\n",
    "    * coarsen the graph in order to merge the most similar SNPs together\n",
    "* SNP similarity encoding:\n",
    "    * associate a new feature vector to SNPs accounting for their similarities (node2vec on this graph)\n",
    "    * use the graph as an input of a GNN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A clustering can be done using the [\"Graclus\" clustering approach](https://www.cs.utexas.edu/users/inderjit/public_papers/multilevel_pami.pdf) implemented [here](https://github.com/rusty1s/pytorch_cluster)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "venv"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}