Meeting in Lyon, 2019, February the 2nd

Meeting du 2 février 2019

Statistics for Single-cell data analysis

Présentation générale de Franck sur l’analyse de données single-cell.

Introduction au Single-Cell

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

Gene expression distribution deconvolution in single-cell RNA sequencing

Bertrand présente le papier de Wang et al. (2018)

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).

Régularisation entropique

Présentation de Joon sur le chapitre 4 de Peyré, Cuturi, and others (2017)

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)\)\((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.

Reference

Peyré, Gabriel, Marco Cuturi, and others. 2017. “Computational Optimal Transport.”

Wang, Jingshu, Mo Huang, Eduardo Torre, Hannah Dueck, Sydney Shaffer, John Murray, Arjun Raj, Mingyao Li, and Nancy R Zhang. 2018. “Gene Expression Distribution Deconvolution in Single-Cell Rna Sequencing.” Proceedings of the National Academy of Sciences 115 (28). National Acad Sciences: E6437–E6446.

Avatar
SingleStatOmics
ANR project 2019 - 2022