Enjeux du single-cell

• données de haut-débit en biologie cellulaire (en complément du haut-débit en biologie moléculaire)
• jusqu’à présent, on avait accès à la diversité de forme, de location, de fonction
• accès à la description moléculaire de la variabilité de la cellule, donc des tissus

\(\rightsquigarrow\) accès à de multiples technologies de la biologie moléculaire au niveau single-cell. Très développés en immunologie et cancerologie.

Questions statistiques

• grand \(n\) (jusqu’à 1,000,000), grand \(p\) (\(approx 20,000\) chez l’homme).
• modèles de comptages (GLM) plutôt que gaussiens mais
  • sur-dispersion
  • excès de zéros (soit pas d’expression, soit pas atteint par l’échantillonage)

Méthode de réduction de dimension pour données de comtptage

• Présentation des variantes de l’ACP et de la contribution de Franck et Ghislain.
• La NMF fonctionne aussi bien mais il est plus compliqué d’ajouter la sur-dispersion et l’excès de 0.
• Rappel sur le Stochastic Neighbor Embedding (SNE) et sa variante tSNE

Modèle

Soient les paramètres suivants permettant de décrire le processus de transcription:

• \(p_{1cg}\) décrit la probabilité de reverse transcription pour le gène \(g\) dans la cellule \(c\),
• \(p_{2cg}\) décrit la probabilité pour qu’un brin transcrit soit amplifié et séquencé pour le gène \(g\) dans la cellule \(c\),

Le modèle proposé est alors le suivant :

• \(\lambda_{cg}\) la vraie expression génétique, i.e. la quantité d’ARNm associée au gène \(g\) dans la cellule \(c\), telle que $_{cg}H_g $;
• la qualité de l’ADNc est décrit par \(W_{cg} \sim b(\mathrm{p}{1cg}, \lambda_{cg})\),
• le nombre de lectures (read) observés \(Y_{cg} \sim b(\mathrm{p}{2cg}, W_{cg})\),

Ainsi,

\[ Y_{cg} \sim b(\lambda_{cg}, \underbrace{\alpha_{cg}}_{p_{1cg} \times p_{2cg}}) \xrightarrow[ \alpha_{cg}\to 0]{} \mathcal{P}(\alpha_{cg}\lambda_{cg})\]

Inférence

La distribution \(H_g\) est estimée par un modèle non-paramétrique (mélange de spline cubic et d’un dirac en 0).

Notations et rappels

Soient

• \(i \in \{1,\dots,n\}\) un ensemble de mines
• \(j \in \{1,\dots,m\}\) un ensemble d’usines
• \(\Sigma_n = \{\mathbf{a} \in \mathbb{R}_+^n | \sum_i a_i = 1\}\)

On se donne

• \(a \in \Sigma_n\) et \(b\in\Sigma_m\)
• \(\cup(a,b)\) un ensemble de couplage tels que

\[\cup(a,b) = \left\{\mathbf{P} \in \mathbb{R}_+^{n\times m } \ | \ \mathbf{P} \mathbb{1}_m = a \ \mathrm{ et } \ \mathbf{P}^\top \mathbb{1}_n = b \right\}\]

• \(\mathcal{C} = (C_{ij}) \in \mathbb{R}^{n\times m}\), avec \(C_{ij}\) le coût pour déplacer une unité de masse de \(i\) à \(j\).

Le problème de Kantorovich est le suivant :

\[\mathrm{L}_\mathcal{C}(a,b) = \min_{\mathbb{P}\in \cup(a,b)} \langle \mathcal{C}, \mathbb{B} \rangle \]

Régularisation entropique

Motivation: Il s’agit de rendre le problème

• Calculable
• Parallélisable (GPU)
• différentiable en \((a,b)\)
  • + convexe
  • + régulier

On considère

\[ H(\mathbb{P}) = - \sum_{i,j} P_{ij} (\log P_{ij} - 1), \quad P_{ij}$ \geq 0\]

avec la convention \(0 \log 0 = 0\). Soit \(\varepsilon > 0\) un paramètre. Le problème régularisé s’écrit

\[ \mathrm{L}^\varepsilon_\mathcal{C}(a,b) = \min_{\mathbb{P}\in \cup(a,b)} \left\{ \langle \mathcal{C}, \mathbb{P} \rangle - \varepsilon H(\mathbb{P}) \right\} \]

Proposition Gominetti - San Martin ’94

\[ \mathbb{P}^\varepsilon \xrightarrow[\varepsilon\to 0]{} \arg \min_{\substack{\mathbb{P}\in\cup(a,b) \\ \langle \mathbf{C}, \mathbb{P}\rangle}} - H(\mathbb{P}) \] En particulier, \[ \mathrm{L}_{\mathcal{C}}^\varepsilon \xrightarrow[\varepsilon\to 0]{} \mathrm{L}_{\mathcal{C}} \]

Algorithme de Sinkhorn

Notation \(K = (K_{ij}) = \left(e^{-C_{ij}/\varepsilon}\right)_{ij}\)

Proposition Il existe \(u\in\mathbb{R}_+^n\) et \(v\in\mathbb{R}_+^m\) tels que

1. \(P_{ij}^\star = u_i K_{ij} v_k\)
2. \(u \odot (Kv) = a\); \(v \odot (K^\top u) = b\)

Algorithme

• \(v^{(0)} = \mathbb{1}_m\)
• \(u^{(t+1)} = \frac{a}{K v^{(t)}}\) (division élément par élément)
• \(v^{(t+1)} = \frac{b}{K^\top u^{(t+1)}}\)

Remarque: Soit \(T\) le nombre d’itérations. \(\mathbb{P}^{(T)} \triangleq \mathrm{diag}(u^{(T)}) K \mathrm{diag}(v^{(T)}\) n’appartient pas nécessairement à \(\cup(a,b)\). On peut trouver \(\hat{\mathbb{P}}^{(T)} \in \cup(a,b)\) proche de \(\mathbb{P}^{(T)}\) (il s’agit de l’étape de “Rounding Step”).

Théorème (Franklin et Lorentz, ’89). On note

\[\eta(K) = \max_{i,j,k,\ell} \frac{K_{ik} K_{j\ell}}{K_{jk}K_{i\ell}}, \qquad \lambda(K) = \frac{\sqrt{\eta(K)} - 1}{\sqrt{\eta(K)} + 1} < 1 \] alors \[ \left\| \log P^{(t)} - \log P^\star \right \|_\infty = \mathcal{O}(\lambda(K)^{2t}) \]

Remarque (Calcul de \(\mathrm{L}_\mathcal{C}(a,b)\): complexité, cas \(n=m\)). Soit \(\tau>0\). On peut obtenir \(\hat{\mathbb{P}}\in \cup(a,b)\) tel que

\[ \langle \hat{\mathbb{P}}, \mathcal{C}\rangle \leq \mathrm{L}_\mathcal{C}(a,b) + \tau \] à l’aide de \(\varepsilon = \frac{\tau}{4\log(n)}\) en \[ \underbrace{\mathcal{O}(\|\mathcal{C}\|_\infty^3)\log(n)\tau^{-3}}_{\text{itération de Sinkhorn}} + \underbrace{\mathcal{O}(n^2)}_{\text{rounding step}} \]

Remarque: (Sinkhorn parallèle/GPU). Soient - \(a_1,\dots,a_N \in \Sigma_n\) - \(b_1,\dots,b_N \in \Sigma_m\)

On souhaite calculer \(\mathrm{L}^\varepsilon_{\mathcal{C}}(a_i,b_i)\) pour \(i=1,\dots,N\). En notant \(A = (a_1|\dots|a_N)\) et \(B = (b_1|\dots|b_N)\), idem pour \(U^{(t)},V^{(t)}\), on peut écrire \[ U^{(t)} = \frac{A}{K V^{(t)}}, \quad V^{(t)} = \frac{B}{K^\top U^{(t)}}\]

Dualité

Proposition Le problème dual s’écrit

\[\mathrm{L}^\varepsilon_{\mathcal{C}} \max_{\substack{f\in\mathbb{R}^n \\ g \in \mathbb{R}^m}} Q(f,g), \quad Q(f,g) = \left\{ \langle f,a\rangle + \langle g, b\rangle - \varepsilon \langle e^{f/\varepsilon}, e^{g/\varepsilon}\rangle\right\} \] Si \((f,g)\) est une solution, alors: \[ \mathbb{P}^\star = \mathrm{diag}(\underbrace{e^{f/\varepsilon}}_{u}) K \mathrm{diag}(\underbrace{e^{g/\varepsilon}}_{v})\]

Proposition \((a,b) \rightarrow \mathrm{L}_{\mathcal{C}}^\varepsilon (a,b)\) est convexe et différentiable de gradient \(\nabla \mathrm{L}_{\mathcal{C}}^\varepsilon (a,b) = (f,g)\) où \((f,g)\) est la seule solution du problème dual tel que \(\sum_i f_i = \sum_j g_j = 0\).

Remarque Sinkhorn \(\Leftrightarrow\) Block coordinate descent dans le dual

Les gradients en \(f\) et \(g\) s’écrivent

\[ \nabla_f Q(f,g) = a - e^{f\varepsilon} \odot K e^{g\varepsilon} \]

\[ \nabla_g Q(f,g) = b - e^{g\varepsilon} \odot K^\top e^{f\varepsilon} \] Ainsi \[\left\{\begin{array}{l} f^{(t+1)} : \nabla_f Q(f^{(t)},g^{(t)}) = 0\\ g^{(t+1)} : \nabla_f Q(f^{(t)},g^{(t)}) = 0\\ \end{array}\right. \Leftrightarrow \left\{\begin{array}{l} f^{(t+1)} = \varepsilon \log (a) - \varepsilon \log (K e^{g^{(t)}/\varepsilon})\\ g^{(t+1)} = \varepsilon \log (b) - \varepsilon \log (K^\top e^{f^{(t+1)}/\varepsilon})\\ \end{array}\right. \] Et après un peu de touaillage on retrouve les valeurs de \(u^{(t)}, v^{(t)}\) de l’algorithme de Sinkhorn.