L'objectif est d'implémenter une méthode de réduction de dimension par Diffusion Maps pour des données non linéaires en grande dimension, et de comparer cette méthode à d'autres méthodes communément utilisées.
Pour le projet, vous utiliserez un jeu de données chez l'embryon de souris sur des cellules uniques du mesoderme ayant le potentiel de différencier en cellules du sang, à quatre stades de développement entre E7.0 et E8.5, publié dans
Moignard et al. (2015) Decoding the regulatory network of early blood development from single-cell gene expression measurements, Nature Biotechnology 33, 269–276
Par exemple, les cinq première lignes du fichier
Cell Cbfa2t3h Cdh1 Cdh5 ...
HFA1_001 -2,598855758 -14 1,815101273 ...
HFA1_002 -14 -14 -14 ...
HFA1_003 -14 -14 -14 ...
HFA1_004 -14 -14 -14 ...
HFA1_005 -3,907257698 -14 0,019497022 ...
...
Chaque ligne correspond à une cellule dont l'identifiant est donné dans la première colonne. Les autres colonnes donnent l'expression des gènes pour chacunes de ces cellules (une valeur de -14 indique une expression non-détectée ou non mesurée). L'identifiant de chaque cellule contient l'étape de développement d'où provient la cellule: jour embryonnaire 7, E7.0 (primitive streak, PS), E7.5 (neural plate, NP), E7.5 (head fold, HF). En plus, à E8.5, deux subpopulations de cellules sanguines GFP${}^+$ (four somite, 4SG) et endothéliales Flk${}^+$GFP${}^-$ (4SFG${}^-$).
Télécharger le jeu de données dCt_values.tab.
Lancez un notebook python (où une session python en ligne de commande). Vous pourrez avoir besoin des packages suivants:
import csv
import numpy as np
%matplotlib inline # pour ipython notebook seulement. Pour python en ligne de commande, utilisez pl.show()
import matplotlib as plt
import pylab as pl Importez les fichiers.
Sûrement pas la façon la plus élégante, mais voici néanmoins une méthode pour importer les données.
Les identifiants des cellules sont stockées dans cell, et les expressions dans l'array
d. La liste head contient les noms des gènes.
fileName = "dCt_values.tab"
head = []
data = []
cell = []
with open(fileName, 'rb') as f:
reader = csv.DictReader(f, delimiter = '\t')
# create a dict out of reader, converting all values to integers
for i, row in enumerate(list(reader)):
data.append([])
for key, value in row.iteritems():
if key == "Cell":
cell.append(value)
else:
data[i].append(float(value))
for key, value in row.iteritems():
if key != "Cell":
head.append(key)
d = np.array(data)
Toutes les méthodes de réduction de dimension cherchent à établir un système de coordonnées dans lequel les données peuvent être visualisées et caractérisées en petite dimension.
La diffusion map est une méthode de réduction de dimension qui exploite les liens entre la diffusion de la chaleur et les chaînes de Markov décrivant une marche aléatoire. Une chaîne de Markov est un processus aléatoire en temps discret $x(t)$, $t=0, 1, ...$, dont l'état à $t+1$ ne dépend que de l'état à $t$, et pas des états passés. Ici on i considère L'idée est qu'un jeu de données peut être décrit par des particules diffusant (effectuant une marche aléatoire) dans l'espace. Deux points de l'espace seront "proches" si la probabilité de diffuser d'un point à l'autre est élevé.
On considère un ensemble de $n$ points $\{ x_i \}_{i = 1,...,n}$ avec chaque point $x_i \in \mathbb{R}^p$ (en général, $n \gg p$).
On définit un noyau de similarité $k(x,y)$ pour $x, y \in \{x_i\}_i$ comme une fonction de deux variables à valeur réelle, positive et symétrique:
Le noyau de similarité $k$ n'est pas une distance, en particulier $k(x,x)$ est en général strictement positif, et ne respecte pas l'inégalité du triangle. Plutôt, le noyau définit la géométrie local du jeu de donnée. Un choix populaire est le noyau gaussien $$k(x,y) = \exp \Bigl( - \frac{||x-y||^2}{\varepsilon} \Bigr).$$ On définit la matrice de similarité $L \in \mathbb{R}^{n \times n}$ par $$L_{i,j} = k(x_i, x_j).$$ Par les propriétés du noyau $k$, on a que $L$ est une matrice réelle symétrique à coefficients non-négatifs.
On normalise ensuite la matrice de similarité de façon à ce que la somme de chaque ligne soit égale à 1. On note par $M$ cette matrice normalisée $$M = D^{-1} L,$$ où $D$ est une matrice diagonale positive, $$ D_{i,i} = \sum_{j=1}^{n} L_{i,j}.$$ et $$D_{i,j} = 0, i \neq j.$$
La matrice $M$ possède les caractéristique suivantes
On appelle une matrice qui possède ces deux propriétés une matrice stochastique. Cependant, $M$ n'est pas symétrique. La normalisation de $M$ permet d'interpréter chaque coefficient comme la probabilité de faire un saut aléatoire d'un point vers un autre. La probabilité de sauter du point $x_i$ au point $x_j$ est $M_{i,j}$. La probabilité de sauter de $x_i$ vers n'importe quel point est la somme de la ligne $j$ de $M$, et est égale à 1. On note que si $\mathbf{1}$ est le vecteur avec coefficients 1 partout, $$ M \mathbf{1} = \mathbf{1}.$$ Donc $M$, comme toute matrice stochastique, possède une valeur propre égale à 1 associée à un vecteur propre unité $\mathbf{1}$. Cette propriété sera importante pour la suite.
Soit une position initiale $x_i$, la probabilité de se trouver à la position $x_j$ après $t \in \mathbb{N}$ sauts est donnée par $(M^t)_{i,j}$. En effet, si $M$ est une matrice stochastique, $M^t$ est aussi une matrice stochastique. La probabilité diffuser d'un point du jeu de données à un autre est donc déterminée par les puissances de $M$.
Montrez que si $M$ est une matrice stochastique, alors pour tout entier $t \geq 0$, $M^t$ est une matrice stochastique. ($M^0 = I$ la matrice identité. Montrez ce résultat par récurrence, avec base $t=1$.)
Plus généralement, pour une vecteur ligne $\mathbf{v}^*(t) \in \mathbb{R}^n$ où $\mathbf{v}_i$ est la probabilité d'être au point $x_i$ au temps $t$, on a que la probabilité d'être au point $x_j$ à $t+1$ est donnée par le vecteur $$\mathbf{v}^*(t+1) = \mathbf{v}^*(t) M.$$
Un vecteur de probabilité ou un vecteur stochastique est un vecteur ligne $\mathbf{v}^*$ dont tous les coefficients sont non-négatifs et $$\sum_i \mathbf{v}_i = 1.$$ Montrez que si $\mathbf{v}^*$ est un vecteur de probabilité et $M$ une matrice stochastique alors $\mathbf{v}^* M$ est un vecteur de probabilité
Cette équation discrete définit un système de récurrence sur la distribution des probabilités de la chaîne de Markov des $x_i$. Pour les chaînes de Markov, il est usuel de considérer le produit vecteur ligne à gauche de la matrice. Le résultat du produit vecteur-matrice est un vecteur ligne (noté $*$). Une distribution de probabilité $\mathbf{\pi}$ est dite stationnaire si $$ \mathbf{\pi}^* M = \mathbf{\pi}^*.$$ Une distribution stationnaire est un vecteur propre à gauche de la matrice $M$, associé à la valeur propre 1. On peut noter que $\mathbf{\pi}^* M = ( M^* \mathbf{\pi} )^* = (\mathbf{\pi})^*$, et donc que $\mathbf{\pi}$ est un vecteur propre de la matrice $M^*$. On a noté plus haut que la valeur propre 1 est aussi valeur propre à droitei. En fait, pour toute matrice carrée, les valeurs propres à gauche et à droite sont les mêmes. Les vecteurs ne le sont pas par contre. Comme on sait qu'une matrice stochastique possède une valeur propre 1, il existe au moins un vecteur propre à gauche avec valeur propre 1. De plus, aucune valeur propre ne peut être plus grande que 1. Enfin, on peut montrer qu'il existe au moins un vecteur propre à gauche, avec valeur propre 1, qui est aussi vecteur de probabilité. Ce vecteur n'est pas nécessairement unique (pensez à la matrice identité: tous les vecteurs de probabilité sont vecteurs stationnaires). Cependant, si pour tout $i$ et $j$, $M_{i,j} > 0$, on a l'existence d'un vecteur stationnaire unique. Pour le projet, nous allons construire des matrices stochastiques possédant cette propriété (en prenant le noyau gaussien par exemple). La distribution stationnaire est alors donnée explicitement par $$\pi = \frac{D_{i,i}}{\sum_j D_{j,j}}.$$ (Vous pouvez vérifier que $\pi$ satisfait $\pi^* = \pi^* M$.)
Pour récapituler, notre matrice stochastique $M$ possède une plus grande valeur propre 1, avec vecteur propre à droite le vecteur unité $\mathbb{1}$ et vecteur propre à gauche $\mathbf{\pi}^*$, le vecteur stationnaire. On supposera pour le projet que le vecteur stationnaire est unique.
Pour caractériser plus finement la matrice stochastique $M$, on note que $M$ est semblable à une matrice réelle symétrique à coefficients positifs: $$M_s = D^{1/2} M D^{-1/2}.$$ La matrice $D^{1/2}$ est la matrice diagonale inversible avec coefficients diagonaux $\sqrt{D_{i,i}}$. De plus les matrices $M$ et $M_s$ partagent les mêmes valeurs propres.
Montrez que $M_s$ est réelle et symétrique.
La matrice $M_s$ est diagonalisable et possède une décomposition en valeur propres $M_s = V \Lambda V^{-1}$ avec vecteurs propres orthogonaux et valeurs propres réelles.
En général, la matrice $M_s$ n'aura que des valeurs propres positives.
Montrez que $M_s$ possède
Si $M_s$ possède des valeurs propres multiples, on peut choisir une base orthogonale pour le sous-espace propre. On peut donc former une matrice de passage $V$ unitaire, pour obtenir la décomposition $M_s = V \Lambda V^*$. Comme $M$ est semblable à $M_s$, on obtient aussi une décomposition pour $M$, $$M = D^{-1/2} V \Lambda V^* D^{1/2}.$$ La décomposition de $M$ se fait par deux familles de vecteurs bi-orthogonaux. Notons $v_j$ la $j$eme colone de $V$, $$ \psi_j = \frac{v_j}{\sqrt{D_{j,j}}},$$ $$ \phi_j^* = \sqrt{D_{j,j}} v_j^*.$$ (16/01) erreur sur la définition de $\phi_j$ corrigée. Les deux familles de vecteurs $\psi$ et $\phi$ sont bi-orthogonales: $\phi_i^* \psi_j = 0$ si $i \neq j$ et $\phi_i^* \psi_i = 1$.
Montrez que les familles de vecteurs $\psi$ et $\phi^*$ sont bi-orthogonales.
Les vecteurs $\phi^*$ sont les vecteurs propres à gauche de $M$ et $\psi$ sont les vecteurs propres à droite de $M$. Si on classe les valeurs propres en ordre décroissant de valeurs absolues, on obtient que le premier vecteur propre à gauche est nul autre que la distribution stationnaire $\mathbf{\pi}^*$. Le premier vecteur propre à droite est le vecteur $\mathbb{1}$. On peut exprimer la probabilité d'être au point $x_j$ au temps $t$ sachant qu'on est parti de $x_i$ au temps 0. $$p(x_j,t|x_i,0) = (M^t)_{i,j} = \sum_{\ell=0}^{n-1} \lambda_\ell^t \psi_{i,\ell} \phi_{j,\ell}.$$ (On indexe les vecteurs de 0 à $n-1$, la raison sera évidente bientôt.) On sait que $\psi_0 = \mathbb{1}, \lambda_0 = 1$ et on peut réorganiser les termes comme $$p(x_j,t|x_i,0) = \phi_{0,j} + \sum_{\ell=1}^{n-1} \lambda_\ell^t \psi_{i,\ell} \phi_{j,\ell}.$$ Comme les valeurs propres sont décroissantes, il est raisonnable d'approximer cette somme en ne gardant que les $k$ premiers termes. En fait, comme on considère la une marche aléatoire sur les points $x_i$ on peut considérer la vitesse à laquelle deux marches aléatoires divergeront en fonction de leur état initiale. Pour deux points $x_0, x_1 \in \mathbb{R}^p$, la distance de diffusion est $$D_t^2(x_0, x_1) = \sum_{j=1}^n \frac{\bigl( p(x_j,t|x_0) - p(x_j,t|x_1) \bigr)^2}{\phi_{0,j}}.$$
On dénote la diffusion map comme l'application qui envoie les $n$ points $\{x_i\}_{i=1,...,n}, x_i \in \mathbb{R}^p$ vers les $k$ premières composante de l'espace propre (exluant le premier vecteur propre trivial) $$\Psi_t(x) = \bigl(\lambda_1^t \psi_1(x), \lambda_2^t \psi_2(x),..., \lambda_k^t \psi_k(x) \bigr).$$ avec $\psi_j(x)$ le $i$eme coefficient du vecteur $\psi_j$ si $x = x_i$. La diffusion map est justifié par le résultat suivant (pas à démontrer), qui relie la distance de diffusion et la diffusion map
Thm. La distance de diffusion est égale à la norme euclidienne de dans l'espace de la diffusion map avec tous les $k=(n-1)$ vecteurs propres $$D_t^2(x_0, x_1) = \sum_{j=1}^n \lambda_j^{2t} \bigl( \psi_j(x_0) - \psi_j(x_1) \bigr)^2 = ||\Psi_t(x_0) - \Psi_t(x_1||^2.$$
Implémentez un algorithme pour calculer la diffusion map d'un jeu de données
$\{x_i\}_{i=1,...,n}, x_i \in \mathbb{R}^p$.
Utilisez le noyau gaussien tel que définit plus haut.
Utilisez par exemple python pour définir une fonction qui prend en entrée
un jeu de donnée x, un paramètre $\varepsilon$, et qui renvoie
les vecteurs propres $\psi$, le vecteur propre $\phi_0$, et les valeurs propres $\lambda$
ordonnées en valeurs absolues décroissantes. Définissez une deuxième fonction
qui prend en entrée les vecteur propres, les valeurs propres, un entier $k$ et un entier
$t$, et qui renvoie les $k$ vecteur de la diffusion map au temps $t$.
Le étapes vont comme suit:
Calculez la diffusion map pour le jeux de données de cellules sanguines. Essayez de trouver $\varepsilon$ et $t$ qui reproduisent le mieux les résultats de Moignard et al. dans la figure 2
Comparez la réduction de dimension par diffusion map et par analyse en composantes principales. (Plus de détails à venir)