[PYTHON] Essayez d'utiliser PHATE, une méthode de réduction et de visualisation des données biologiques

Essayez d'utiliser PHATE, une méthode de réduction et de visualisation des données biologiques

PHATE (Moon, KR, van Dijk, D., Wang, Z. et al. Nature Biotechnology 37, 1482-1492 (2019)) Essayez d'utiliser -3). Il existe des tonnes de méthodes de réduction de dimension utilisées dans les articles de biologie, mais chaque méthode a ses avantages et ses inconvénients. Par exemple, les méthodes suivantes sont souvent utilisées comme méthodes typiques.

1. PCA: analyse en composantes principales. Retirez l'axe avec la dispersion maximale. Les autres méthodes n'ont pas divers avantages, mais surtout lorsqu'elles sont utilisées pour la visualisation, les entités non linéaires ne peuvent pas être capturées et la structure locale qui reflète les caractéristiques globales de la distribution a tendance à être écrasée sous forme de bruit. Il y a des inconvénients.
2. t-SNE (et UMAP): recherche des emplacements de faible dimension afin que la relation de distance avec le voisinage local soit préservée pour chaque point. Par conséquent, la structure locale de la distribution des données devient une visualisation bien préservée. Si les données sont composées de plusieurs clusters, la visualisation reflètera bien ces clusters, mais les «trajectoires» continues ont tendance à être divisées de manière appropriée. De plus, comme nous voyons rarement des relations à distance, les relations de distance relative entre plusieurs clusters changent avec des nombres aléatoires.
3. MDS: méthode de construction à échelle multidimensionnelle. Conçu pour stocker les relations de distance de tout à tout entre les points de données Il est vulnérable au bruit car il ne regarde pas directement la «structure» des données (il ne fait pas de distinction entre le bruit et les structures caractéristiques telles que les clusters et les orbites).
4. Carte de diffusion: Construisez une matrice de probabilité de transition à partir de la distance entre les points de données et opérez le processus de diffusion (marche aléatoire sur la structure graphique des données). La distribution après t étapes est obtenue par la t-ième puissance de la matrice de probabilité de transition. Je veux obtenir des coordonnées de faible dimension qui reflètent la différence de distribution comme la distance euclidienne, mais qui peuvent être obtenues comme vecteur propre de la matrice de probabilité de transition à la t-ième puissance. Il présente divers avantages tels que la résistance au bruit, la capacité de gérer des caractéristiques locales aux structures globales par le nombre de pas t, et la facilité de trouver des structures telles que des "orbites", mais il a tendance à comprimer différentes orbites à des dimensions différentes. En d'autres termes, il ne convient pas à la visualisation en petites dimensions telles que 2D et 3D.
5. Inférence de trajectoire telle que Monocle2: Comme il s'agit d'une méthode d'estimation de la structure arborescente, il n'est pas clair à l'avance si la prémisse tient. L'estimation est également instable.

etc. Une méthode appelée PHATE (Potential of Heat diffusion for Affinity-based Transition Embedding) a été développée pour éliminer ces inconvénients. De plus, M-PHATE, qui est une extension de PHATE, a été développé comme une incorporation plus générale de tenseurs. Il s'agit d'une méthode de visualisation de faible dimension pour les données qui inclut l'évolution temporelle du même groupe (comme les informations de poids pour chaque époque d'un réseau neuronal).

L'installation est facile, il suffit de pip install phate '' `` Voir PHATE Repository pour la version R, la version Matlab, etc.

Visualisation des données de l'arborescence

Expérimentons avec des données d'arborescence à 100 dimensions (+ bruit) fournies par PHATE. Il se compose de 20 "branches" et compte au total 2 000 points de données. En ce qui concerne la structure, les données montrent que les 19 branches restantes se développent de manière aléatoire à partir de diverses parties de la première branche. Chaque branche fait face à différentes directions dans un espace de 100 dimensions.

import phate

tree_data, tree_clusters = phate.tree.gen_dla(seed=108)

L'utilisation est la même que le modèle scikit-learn. Tout d'abord, définissez les paramètres, créez une instance du modèle et transformez les données avec fit_transform ''. Il existe trois paramètres principaux: knn```, decay, `` `t. t``` peut également être déterminé automatiquement à partir des données. Les détails seront décrits plus tard.

phate_operator = phate.PHATE(knn=15, decay=40, t=100)

tree_phated = phate_operator.fit_transform(tree_data)

    Calculating PHATE...
      Running PHATE on 2000 cells and 100 genes.
      Calculating graph and diffusion operator...
        Calculating KNN search...
        Calculated KNN search in 0.50 seconds.
        Calculating affinities...
        Calculated affinities in 0.05 seconds.
      Calculated graph and diffusion operator in 0.56 seconds.
      Calculating diffusion potential...
      Calculated diffusion potential in 0.41 seconds.
      Calculating metric MDS...
      Calculated metric MDS in 10.43 seconds.
    Calculated PHATE in 11.41 seconds.

Comprimé en deux dimensions. Essayez de tracer les résultats.

import matplotlib.pyplot as plt
import seaborn as sns

def plot_2d(coords, clusters, ax):
    for c, color in zip(sorted(list(set(clusters))), plt.cm.tab20.colors):
        samples = clusters == c
        ax.scatter(coords[samples, 0], coords[samples, 1], \
                   label=c, color=color)
    ax.legend(loc='upper right', bbox_to_anchor=(1.1, 1))
    sns.despine()

fig, ax = plt.subplots(figsize=(12, 9))
plot_2d(tree_phated, tree_clusters, ax=ax)
plt.show()

De cette manière, la structure des données peut être maintenue et chaque branche peut être correctement séparée et visualisée. Comparons la visualisation de PCA, t-SNE, UMAP, PHATE avec les mêmes données.

from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
from umap import UMAP

def compare_vis(data, clusters, pca_init_dim=100):
    print('Running PCA...')
    pca_op = PCA(n_components=2)
    data_pca = pca_op.fit_transform(data)
    print('Running t-SNE...')
    pca_init_op = PCA(n_components=pca_init_dim)
    tsne_op = TSNE(n_components=2)
    data_tsne = tsne_op.fit_transform(pca_init_op.fit_transform(data))
    print('Running UMAP...')
    pca_init_op = PCA(n_components=pca_init_dim)
    umap_op = UMAP(n_components=2)
    data_umap = umap_op.fit_transform(pca_init_op.fit_transform(data))
    print('Running PHATE...')
    phate_op = phate.PHATE(n_components=2)
    data_phated = phate_op.fit_transform(data)
    fig, axes = plt.subplots(figsize=(12, 36), nrows=4)
    plot_2d(data_pca, clusters, ax=axes[0])
    plot_2d(data_tsne, clusters, ax=axes[1])
    plot_2d(data_umap, clusters, ax=axes[2])
    plot_2d(data_phated, clusters, ax=axes[3])
    axes[0].set_title('PCA'); axes[1].set_title('t-SNE')
    axes[2].set_title('UMAP'); axes[3].set_title('PHATE')
    plt.show()    

compare_vis(tree_data, tree_clusters)

Dans PCA, la structure arborescente n'est pas bien comprise (en raison de la grande influence du bruit), et dans t-SNE et UMAP, les clusters sont divisés de manière appropriée. Regardons également l'intégration en 3D. Comme cela sera décrit plus tard, dans l'algorithme de PHATE, la dimension intégrée finale n'est liée qu'au calcul MDS final, il n'est donc pas nécessaire de recalculer le tout, et le paramètre (n_components '') est mis à jour en ''. Si vous transformez '', seule la dernière étape sera exécutée.

%matplotlib notebook
from mpl_toolkits.mplot3d import axes3d

def plot_3d(coords, clusters, fig):
    ax = fig.gca(projection='3d')
    for c, color in zip(sorted(list(set(clusters))), plt.cm.tab20.colors):
        samples = clusters == c
        ax.scatter(coords[samples, 0], coords[samples, 1], coords[samples, 2], \
                   label=c, color=color)
    ax.legend()

phate_operator.set_params(n_components=3)
tree_phated = phate_operator.transform()

fig = plt.figure(figsize=(12, 9))
plot_3d(tree_phated, tree_clusters, fig)

Visualisation des données scRNA-seq

Essayez-le avec les données embryonnaires scRNA-seq utilisées pour la démonstration dans l'article. Les données sont un peu trop grandes, et également pour évaluer la stabilité du résultat lors d'un sous-échantillonnage, 10% du total des cellules sont échantillonnées au hasard et le résultat de l'exécution d'un prétraitement ordinaire est effectué à l'avance. C'était. Pour plus de détails sur le prétraitement des données scRNA-seq et les algorithmes tels que PCA, t-SNE et UMAP, voir Course Video.

import numpy as np
import pandas as pd
import scipy

cells = np.genfromtxt('./data/cell_names.tsv', dtype='str', delimiter='\t')
genes = np.genfromtxt('./data/gene_names.tsv', dtype='str', delimiter='\t')
arr = scipy.io.mmread('./data/matrix.mtx')

EBData = pd.DataFrame(arr.todense(), index=cells, columns=genes)
EBData.head()

Tableau des niveaux d'expression génique de 1 679 cellules et 12 993 gènes. Chaque cellule est étiquetée comme suit. Puisqu'il est pris tous les 3 jours jusqu'au 27ème jour, la série chronologique de différenciation devrait être observable.

sample_labels = EBData.index.map(lambda x:x.split('_')[1]).values
print(sample_labels)

    array(['Day 00-03', 'Day 00-03', 'Day 00-03', ..., 'Day 24-27',
           'Day 24-27', 'Day 24-27'], dtype=object)

phate_op = phate.PHATE()
EBData_phated = phate_op.fit_transform(EBData)

%matplotlib inline
fig, ax = plt.subplots(figsize=(12, 9))
plot_2d(EBData_phated, sample_labels, ax=ax)
plt.show()

Le processus de différenciation peut être clairement compris par réduction de dimension + visualisation par PHATE. Surtout, la diffusion de la diversité du 18e au 27e jour. Avec ce genre de sentiment, le bon point de PHATE est que lorsqu'il y a une orbite, c'est une orbite appropriée, quand il y a un cluster, c'est comme un cluster, et la propagation de l'amas est également visualisée afin qu'ils puissent être comparés les uns aux autres. De plus, bien que seulement 10% de toutes les données soient utilisées, la figure du papier qui est dimensionnellement compressée à l'aide de dizaines de milliers de cellules (Paper Le dessin étant presque le même que celui de la figure 1c) en -3), on peut voir que la taille de l'échantillon est également robuste dans une certaine mesure. Je vais essayer la 3D.

phate_op.set_params(n_components=3)
EBData_phated = phate_op.transform()
fig = plt.figure(figsize=(9, 9))
plot_3d(EBData_phated, sample_labels, fig)
plt.show()

Comparez les mêmes données scRNA-seq avec PCA, t-SNE, UMAP, etc.

compare_vis(EBData.values, sample_labels)

Après tout, t-SNE et UMAP ont un fort sens du cluster, et il est difficile de comprendre la trajectoire.

Mise en œuvre de PHATE

Essayez de mettre en œuvre PHATE selon la méthode du papier. Le tout étant composé en combinant différentes méthodes, on ne sait pas quelle partie du traitement apporte une contribution décisive à l'amélioration de la visualisation. Alors, divisons-le en plusieurs étapes et implémentons-le. Pour être honnête, je ne suis pas sûr de l'importance de chaque étape ou de la nécessité, mais je ne sais pas si c'est cette implémentation ... Il comprend les quatre étapes suivantes. Si vous connaissez la carte de diffusion, elle peut être plus facile à avaler car elle est presque la même jusqu'à l'étape 2.

1. Calcul de l'opérateur de diffusion

• Convertir en une matrice de distance, calculer l'affinité locale en laissant le noyau mordre, et configurer une matrice de probabilité de transition pour un pas = opérateur de diffusion

2. Détermination de l'échelle de temps de diffusion "t"

• Déterminez le nombre d'étapes de la diffusion à partir de la quantité d'informations contenues dans la t-ième puissance de l'opérateur de diffusion.

3. Exécutez le processus de diffusion et calculez les "distances potentielles"

• Au lieu de mesurer la distance après diffusion avec la distance de diffusion comme dans la carte de diffusion, calculez avec la distance potentielle nouvellement proposée.

4. Intégration de faible dimension de la distance potentielle avec MDS

• Utilisez MDS (méthode de construction à échelle multidimensionnelle) pour transférer la relation de distance de la distance potentielle à la relation de distance de coordonnées arbitraires de faible dimension.

1. Calcul de l'opérateur de diffusion

Tout d'abord, la matrice de distance euclidienne pour l'ensemble des données est calculée. Pour chaque point de données, appliquez un noyau adaptatif similaire à t-SNE ou UMAP afin que les points de données proches les uns des autres aient un poids élevé (= affinité) et les points de données éloignés un poids faible. Créez une structure graphique. Voici les paramètres knnetalpha (appelés` `` decay dans l'implémentation originale de phate) qui affectent grandement le résultat global. .. Le noyau est un noyau gaussien ordinaire, mais comme la perplexité dans t-SNE et n_neighbours dans UMAP, la bande passante est modifiée de manière adaptative en fonction de la densité à proximité des données. Dans le cas de PHATE, cette bande passante adaptative est simplement fixée à la distance $ \ epsilon_k $ de chaque point de données au point de données le plus proche du knn th. Un autre périphérique s'appelle $ \ alpha $ -decaying kernel. Si seule la largeur de bande est modifiée, l'affinité reste à un point éloigné lorsque la largeur est large (queue lourde). Ensuite, les points pris dans la région de faible densité auront presque la même affinité avec tous les autres points. Pour éviter cela, la diminution de l'affinité est contrôlée par le paramètre $ \ alpha $ en utilisant le noyau suivant. $K_{k, \alpha}(x,y)=\frac{1}{2}exp\left( -\left(\frac{\|x-y\|_2}{\epsilon_k(x)}\right)^{\alpha}\right) + \frac{1}{2}exp\left( -\left(\frac{\|x-y\|_2}{\epsilon_k(y)}\right)^{\alpha}\right)$
Dans le cas de $ \ alpha = 2 $, c'est un noyau gaussien normal, mais si $ \ alpha $ est grand, l'affinité diminue plus rapidement. Les deux knn et $ \ alpha $ affectent la connexion du graphe et doivent être déterminés avec soin. Si knn est trop petit, ou si $ \ alpha $ est trop grand et que le graphe est trop clairsemé, le processus de diffusion ne se déroulera pas correctement. Par contre, si knn est trop grand ou si $ \ alpha $ est trop petit, le tout sera connecté avec un poids égal et la structure ne sera pas capturée. Pour le moment, PHATE est réputé robuste en raison de la différence de paramètres, mais je pense qu'il est préférable de corriger l'un ou l'autre et de l'expérimenter.

import numpy as np
from scipy.spatial.distance import pdist, squareform

def diffusion_operator(data, knn=5, alpha=40):
    #Construction d'une matrice de distance euclidienne
    dist = squareform(pdist(data))
    #Calcul de la bande passante adaptative epsilon. Triez la matrice de distance à la valeur knnth.
    epsilons = np.sort(dist, axis=1)[:, knn]
    # alpha-decaying kernel
    kernel_mtx = np.exp(-1.* np.power(dist / epsilons[:, np.newaxis], alpha))
    # L1-Normalisé avec la norme. Cela devient une distribution de probabilité pour chaque ligne.
    l1norm_kernel_mtx = kernel_mtx / np.abs(kernel_mtx).sum(axis=1)[:, np.newaxis]
    #Symétrie.
    transition_mtx = 0.5 * l1norm_kernel_mtx + 0.5 * l1norm_kernel_mtx.T
    return transition_mtx

2. Détermination de l'échelle de temps de diffusion "t"

Le dernier paramètre important est "t" qui détermine l'échelle de temps de diffusion. Décidez du nombre d'étapes à suivre avec la marche aléatoire en fonction de l'affinité (= combien de fois l'opérateur de diffusion doit être élevé). Le processus de diffusion sert également à débruiter les données. Les points de données connectés par des chemins courts à forte affinité sont plus susceptibles de se visiter lors d'une promenade aléatoire, et les points aberrants avec uniquement des chemins à faible affinité sont moins susceptibles d'être visités, de sorte que les premiers se déplacent à chaque fois que l'étape de diffusion progresse. Soyez souligné. Par conséquent, l'échelle de temps de diffusion peut être considérée comme le degré d'élimination du bruit. Si elle est trop petite, beaucoup de bruit restera, et si elle est trop grande, des signaux importants peuvent être supprimés en tant que bruit. Dans PHATE, comme méthode de détermination automatique de l'échelle de temps de diffusion t à partir des données, le bruit a été complètement éliminé en examinant la «quantité d'informations» contenue dans la t-ième puissance de l'opérateur de diffusion à divers t, tandis que les données d'origine. Déterminez t en chronométrant afin que la structure de ne soit pas écrasée. Plus précisément, la valeur propre est calculée pour chaque (conjugué symétrique) de l'opérateur de diffusion à la t-ième puissance pour calculer entropie de von Neumann (VNE). À mesure que t augmente, VNE diminue. La réduction initiale de VNE est due à l'élimination du bruit, et les valeurs propres non essentielles passent rapidement à zéro. Les valeurs propres qui reflètent la structure des données diminuent plus lentement. Par conséquent, si vous sélectionnez t au moment où le taux de réduction change, vous pouvez sélectionner automatiquement t qui supprime le bruit et ne réduit pas la structure. Ici, après avoir calculé le VNE avec différents t, effectuez deux régressions linéaires sur les côtés gauche et droit de chaque t pour détecter le «point de coude» du tracé en choisissant le point avec le moins d'erreur totale. ..

from tqdm import tqdm_notebook as tqdm

def von_Neumann_entropy(P):
    # symmetric conjugate (diffusion affinity matrix)Calculs de
    norm_factor = np.sqrt(P.sum(axis=1))
    A = P / norm_factor[:, np.newaxis] / norm_factor
    #Décomposition de valeur unique
    eig, _ = np.linalg.eigh(A)
    eig = eig[eig > 0]
    eig /= eig.sum()
    #VNE est une entropie de Shannon normalisée aux valeurs propres
    VNE = -1.0 * (eig * np.log(eig)).sum()
    return VNE

def find_knee_point(vals):
    y = vals
    x = np.arange(len(y))
    candidates = x[2:len(y) - 2]
    errors = np.zeros(len(candidates))
    errors[:] = np.nan
    #Pour chaque t, régression linéaire sur les côtés gauche et droit pour calculer l'erreur
    #Il semble y avoir un meilleur moyen, mais ici, je calcule honnêtement un par un...
    for i, pnt in enumerate(candidates):
        # left area linear regression (y = ax + b)
        #Calculez la méthode du carré minimum selon le manuel sur le côté gauche du point.
        x_mean = x[:pnt].mean()
        y_mean = y[:pnt].mean()
        a = ((x[:pnt] - x_mean) * (y[:pnt] - y_mean)).sum() / np.power(x[:pnt] - x_mean, 2).sum()
        b = y_mean - a * x_mean
        left_err = (a * x[:pnt] + b) - y[:pnt]
        # right area linear regression
        x_mean = x[pnt:].mean()
        y_mean = y[pnt:].mean()
        a = ((x[pnt:] - x_mean) * (y[pnt:] - y_mean)).sum() / np.power(x[pnt:] - x_mean, 2).sum()
        b = y_mean - a * x_mean
        right_err = (a * x[pnt:] + b) - y[pnt:]
        # total error
        errors[i] = np.abs(left_err).sum() + np.abs(right_err).sum()
    #Soit le point avec la plus petite erreur le point de genou
    knee_ind = candidates[np.nanargmin(errors)]
    return knee_ind

def find_optimal_timescale(P, max_t=100):
    VNEs = np.zeros(max_t)
    for t in tqdm(range(1, max_t+1)):
        #Calculez la t-ième puissance et la VNE de l'opérateur de diffusion pour chaque t.
        Pt = np.linalg.matrix_power(P, t)
        VNEs[t-1] = von_Neumann_entropy(Pt)
    knee_point_ind = find_knee_point(VNEs)
    optimal_t = knee_point_ind + 1
    fig, ax = plt.subplots()
    ax.plot(range(len(VNEs)), VNEs)
    ax.scatter(knee_point_ind, VNEs[knee_point_ind], marker='o', color='r')
    plt.show()
    return optimal_t

3. Exécutez le processus de diffusion et calculez les "distances potentielles"

Dans le cas de la carte de diffusion, si l'opérateur de diffusion est décomposé en valeurs propres, les coordonnées qui reflètent la distance de diffusion (la distance au carré de chaque point de la distribution de probabilité après t étapes) peuvent être obtenues telles quelles. Cependant, dans ce calcul de distance, la sensibilité des valeurs de faible probabilité est faible, il est donc difficile de refléter des informations sur la distance entre des points distants (structure globale des données). Par conséquent, la distance euclidienne est calculée après conversion logarithmique de l'opérateur de diffusion après diffusion en t. C'est ce qu'on appelle la «distance potentielle». L'article explique en détail ce que signifie physiquement cet indice de distance nouvellement introduit (mais je n'en étais pas sûr).

def potential_distances(P, t):
    #multiplier l'opérateur de diffusion par t
    Pt = np.linalg.matrix_power(P, t)
    # -log(P^t)Convertir et calculer la distance euclidienne
    return squareform(pdist(-1.0 * np.log(Pt + 1e-7)))

4. Intégration de faible dimension de la distance potentielle avec MDS

Au lieu d'une décomposition en valeurs propres comme la carte de diffusion, la distance potentielle calculée en 3 est déplacée vers une dimension inférieure adaptée à la visualisation par MDS (méthode de construction à échelle multidimensionnelle). Premièrement, les coordonnées de faible dimension sont calculées en utilisant le MDS classique (= PCoA), et le MDS métrique est calculé en utilisant cela comme valeur initiale. metric MDS est exécuté par l'algorithme SMACOF de scikit-learn.

from sklearn.manifold import smacof

def mds(dist, n_components=2):
    # classical multidimensional scaling (= PCoA)
    n = len(dist)
    J = np.eye(n) - np.ones((n, n))/n
    B = -J.dot(dist**2).dot(J)/2
    evals, evecs = np.linalg.eigh(B)
    idx   = np.argsort(evals)[::-1]
    evals = evals[idx]
    evecs = evecs[:,idx]
    pos = np.where(evals > 0)[0]
    initial_coords  = evecs[:,pos].dot(np.diag(np.sqrt(evals[pos])))[:, :n_components]
    # metric multidimensional scaling (SMACOF algorithm)
    coords = smacof(dist, n_components=n_components, init=initial_coords, n_init=1)[0]
    return coords

Exécutons en fait les étapes ci-dessus avec les données de l'arborescence.

diffop = diffusion_operator(tree_data, knn=15, alpha=40)
optimal_t = find_optimal_timescale(diffop)
print(optimal_t)

potential_dist = potential_distances(diffop, optimal_t)
coords = mds(potential_dist)

fig, ax = plt.subplots(figsize=(12, 9))
plot_2d(coords, tree_clusters, ax=ax)
plt.show()

Des résultats similaires à l'implémentation originale ont été obtenus. En fait, l'algorithme ci-dessus lui-même n'est pas calculé en concevant des moyens d'améliorer la stabilité du calcul numérique à chaque étape, ou en prenant diverses mesures pour économiser de la mémoire et calculer à grande vitesse, mais il est approximatif. Il semble qu'il fasse des calculs, mais cela ressemble à ceci en règle générale.