Versuchen Sie es mit PHATE, einer Methode zur Reduzierung und Visualisierung biologischer Daten

PHATE (Moon, KR, D. van Dijk, Z. Wang et al., Nature Biotechnology 37, 1482–1492 (2019)) Versuchen Sie es mit -3). In Biologiepapieren werden unzählige Methoden zur Dimensionsreduzierung verwendet, aber jede Methode hat ihre Vor- und Nachteile. Beispielsweise werden die folgenden Methoden häufig als typische Methoden verwendet.

PCA: Hauptkomponentenanalyse. Nehmen Sie die Achse mit der maximalen Streuung heraus. Es gibt verschiedene Vorteile, die andere Methoden nicht haben, aber insbesondere wenn sie zur Visualisierung verwendet werden, können nichtlineare Merkmale nicht erfasst werden, und die lokale Struktur, die die globalen Merkmale der Verteilung widerspiegelt, wird tendenziell als Rauschen zerstört. Es gibt Nachteile.
t-SNE (und UMAP): Suchen Sie nach niedrigdimensionalen Platzierungen, damit die Entfernungsbeziehung zur lokalen Nachbarschaft für jeden Punkt erhalten bleibt. Daher wird die lokale Struktur der Datenverteilung zu einer gut erhaltenen Visualisierung. Wenn Daten aus mehreren Clustern bestehen, spiegelt die Visualisierung diese Cluster gut wider, aber kontinuierliche "Trajektorien" werden in der Regel entsprechend aufgeteilt. Da wir selten Fernbeziehungen sehen, ändern sich auch die relativen Entfernungsbeziehungen zwischen mehreren Clustern mit Zufallszahlen.
MDB: Konstruktionsmethode im mehrdimensionalen Maßstab. Angelegt, um alle Entfernungsbeziehungen zwischen Datenpunkten zu speichern. Es ist anfällig für Rauschen, da es die "Struktur" der Daten nicht direkt betrachtet (es unterscheidet nicht zwischen Rauschen und charakteristischen Strukturen wie Clustern und Umlaufbahnen).
Diffusionskarte: Erstellen Sie eine Übergangswahrscheinlichkeitsmatrix aus dem Abstand zwischen Datenpunkten und führen Sie den Diffusionsprozess aus (Random Walk auf der Diagrammstruktur der Daten). Die Verteilung nach t Schritten ergibt sich aus der t-ten Potenz der Übergangswahrscheinlichkeitsmatrix. Ich möchte niedrigdimensionale Koordinaten erhalten, die den Verteilungsunterschied als euklidischen Abstand widerspiegeln, aber als Eigenvektor der Übergangswahrscheinlichkeitsmatrix zur t-ten Potenz erhalten werden können. Es hat verschiedene Vorteile, wie z. B. Rauschbeständigkeit, Fähigkeit, von lokalen Merkmalen zu globalen Strukturen mit der Anzahl der Schritte t zu gelangen, und leicht zu findende Strukturen wie "Umlaufbahnen", aber es neigt dazu, unterschiedliche Umlaufbahnen auf unterschiedliche Dimensionen zu komprimieren. Mit anderen Worten, es ist nicht für die Visualisierung in kleinen Dimensionen wie 2D und 3D geeignet.
Trajektorieninferenz wie Monocle2: Da es sich um eine Methode zur Schätzung der Baumstruktur handelt, ist im Voraus nicht klar, ob die Prämisse gilt. Die Schätzung ist ebenfalls instabil.

etc. Um diese Nachteile zu beseitigen, wurde eine Methode namens PHATE (Potential of Heat Diffusion for Affinity-based Transition Embedding) entwickelt. Darüber hinaus wurde M-PHATE, eine Erweiterung von PHATE, als allgemeinere Einbettung von Tensoren entwickelt. Dies ist eine niedrigdimensionale Visualisierungsmethode für Daten, die die zeitliche Entwicklung derselben Gruppe umfasst (z. B. Gewichtsinformationen für jede Epoche eines neuronalen Netzwerks).

Die Installation ist einfach, nur `` `pip install phate``` Siehe PHATE Repository für R-Version, Matlab-Version usw.

Visualisierung von Baumstrukturdaten

Experimentieren wir mit 100-dimensionalen Baumstrukturdaten (+ Rauschen), die von PHATE bereitgestellt werden. Es besteht aus 20 "Zweigen" und hat insgesamt 2.000 Datenpunkte. Was die Struktur betrifft, zeigen die Daten, dass die verbleibenden 19 Zweige zufällig aus verschiedenen Teilen des ersten Zweigs wachsen. Jeder Zweig zeigt im 100-dimensionalen Raum in verschiedene Richtungen.

import phate

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

Die Verwendung entspricht dem Scikit-Lernmodell. Stellen Sie zuerst die Parameter ein, erstellen Sie eine Instanz des Modells und transformieren Sie die Daten mit `fit_transform```. Es gibt drei Hauptparameter: knn```, zerfall```, t```. `` t``` kann auch automatisch aus den Daten ermittelt werden. Details werden später beschrieben.

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.

Zweidimensional komprimiert. Versuchen Sie, die Ergebnisse zu zeichnen.

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

Auf diese Weise kann die Struktur der Daten beibehalten und jeder Zweig ordnungsgemäß getrennt und visualisiert werden. Vergleichen wir die Visualisierung von PCA, t-SNE, UMAP, PHATE mit denselben Daten.

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)

In PCA ist die Baumstruktur (aufgrund des großen Einflusses von Rauschen) nicht gut verstanden, und in t-SNE und UMAP sind die Cluster entsprechend unterteilt. Schauen wir uns auch die Einbettung in 3D an. Wie später beschrieben wird, bezieht sich im PHATE-Algorithmus die endgültige eingebettete Dimension nur auf die endgültige MDB-Berechnung, sodass das Ganze nicht neu berechnet werden muss, und der Parameter (`n_components```) wird auf` aktualisiert. Wenn Sie "" "transformieren, wird nur der letzte Schritt ausgeführt.

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

Visualisierung von scRNA-seq-Daten

Probieren Sie es mit den embryonalen scRNA-seq-Daten aus, die für die Demonstration in diesem Artikel verwendet wurden. Die Daten sind etwas zu groß, und um die Stabilität des Ergebnisses beim Downsampling zu bewerten, werden 10% der gesamten Zellen zufällig abgetastet, und das Ergebnis einer gewöhnlichen Vorbehandlung wird im Voraus durchgeführt. Es war. Einzelheiten zur Vorverarbeitung von scRNA-seq-Daten und zu Algorithmen wie PCA, t-SNE und UMAP finden Sie unter Kursvideo.

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

Genexpressionsniveau-Tabelle von 1.679 Zellen und 12.993 Genen. Jede Zelle ist wie folgt gekennzeichnet. Da es alle 3 Tage bis zum 27. Tag eingenommen wird, sollte die Zeitreihe der Differenzierung beobachtbar sein.

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

Der Differenzierungsprozess kann durch Dimensionsreduktion + Visualisierung durch PHATE klar verstanden werden. Insbesondere die Verbreitung der Vielfalt vom 18. bis zum 27. Tag. Mit dieser Art von Gefühl ist der gute Punkt von PHATE, dass wenn es eine Umlaufbahn gibt, es eine richtige Umlaufbahn ist, wenn es einen Cluster gibt, es wie ein Cluster ist und die Ausbreitung des Clusters auch visualisiert wird, damit sie miteinander verglichen werden können. Obwohl nur 10% aller Daten verwendet werden, wird die Abbildung des Papiers mit Zehntausenden von Zellen dimensional komprimiert (Papier). Da die Zeichnung fast dieselbe ist wie in Abb. 1c) in -3), ist ersichtlich, dass die Stichprobengröße in gewissem Maße auch robust ist. Ich werde 3D versuchen.

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

Vergleichen Sie dieselben scRNA-seq-Daten mit PCA, t-SNE, UMAP usw.

compare_vis(EBData.values, sample_labels)

Schließlich haben t-SNE und UMAP einen starken Sinn für Cluster, und es ist schwierig, die Flugbahn zu verstehen.

Implementierung von PHATE

Versuchen Sie, PHATE gemäß der in diesem Dokument beschriebenen Methode zu implementieren. Da das Ganze durch die Kombination verschiedener Methoden zusammengesetzt ist, ist nicht klar, welcher Teil der Verarbeitung einen entscheidenden Beitrag zur Verbesserung der Visualisierung leistet. Teilen wir es also in mehrere Schritte auf und implementieren es. Um ehrlich zu sein, bin ich mir nicht sicher, wie wichtig jeder Schritt oder welche Notwendigkeit ist, aber ich bin mir nicht sicher, ob es sich um diese Implementierung handelt ... Es besteht aus den folgenden vier Schritten. Wenn Sie die Diffusionskarte kennen, ist das Schlucken möglicherweise einfacher, da sie bis zu Schritt 2 fast gleich ist.

Berechnung des Diffusionsoperators

In eine Distanzmatrix konvertieren, die lokale Affinität durch Beißen des Kernels berechnen und eine Übergangswahrscheinlichkeitsmatrix für einen Schritt = Diffusionsoperator konfigurieren

Bestimmung der Diffusionszeitskala "t"

Bestimmen Sie anhand der Informationsmenge, die in der t-ten Potenz des Diffusionsoperators enthalten ist, wie viele Schritte die Diffusion abläuft.

Führen Sie den Diffusionsprozess durch und berechnen Sie "mögliche Entfernungen"

Anstatt den Abstand nach der Diffusion mit dem Diffusionsabstand wie in der Diffusionskarte zu messen, berechnen Sie mit dem neu vorgeschlagenen potenziellen Abstand.

Niedrigdimensionale Einbettung des potenziellen Abstands mit MDS

Verwenden Sie MDS (Multidimensional Scale Construction Method), um die Abstandsbeziehung der potenziellen Entfernung auf die Abstandsbeziehung beliebiger niedrigdimensionaler Koordinaten zu übertragen.

1. Berechnung des Diffusionsoperators

Zunächst wird die euklidische Distanzmatrix für die gesamten Daten berechnet. Wenden Sie für jeden Datenpunkt einen adaptiven Kernel an, der t-SNE oder UMAP ähnelt, sodass Datenpunkte, die nahe beieinander liegen, ein großes Gewicht (= Affinität) und Datenpunkte, die weit entfernt sind, ein kleines Gewicht haben. Erstellen Sie eine Diagrammstruktur. Hier kommen die Parameter `knn``` und` alpha (in der ursprünglichen Implementierung von phate als `` `zerfall bezeichnet), die einen signifikanten Einfluss auf das Gesamtergebnis haben. .. Der Kernel ist ein gewöhnlicher Gaußscher Kernel, aber wie die Ratlosigkeit in t-SNE und n_Nachbarn in UMAP wird die Bandbreite entsprechend der Dichte in der Nähe der Daten adaptiv geändert. Im Fall von PHATE wird diese adaptive Bandbreite einfach auf den Abstand $ \ epsilon_k $ von jedem Datenpunkt zu dem Datenpunkt eingestellt, der dem knn-ten am nächsten liegt. Ein anderes Gerät heißt $ \ alpha $ -decaying Kernel. Wenn nur die Bandbreite geändert wird, bleibt die Affinität bis zu einem weit entfernten Punkt, wenn die Breite breit ist (schwerer Schwanz). Dann haben die Punkte, die aus dem Bereich niedriger Dichte entnommen wurden, fast die gleiche Affinität zu allen anderen Punkten. Um dies zu verhindern, wird die Abnahme der Affinität durch den Parameter $ \ alpha $ unter Verwendung des folgenden Kernels gesteuert. $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)$
Im Fall von $ \ alpha = 2 $ handelt es sich um einen normalen Gaußschen Kernel. Wenn jedoch $ \ alpha $ groß ist, nimmt die Affinität schneller ab. Sowohl knn als auch $ \ alpha $ beeinflussen die Verbindung des Graphen und sollten sorgfältig bestimmt werden. Wenn knn zu klein oder $ \ alpha $ zu groß und der Graph zu dünn ist, läuft der Diffusionsprozess nicht gut ab. Wenn andererseits knn zu groß oder $ \ alpha $ zu klein ist, wird das Ganze mit gleichem Gewicht verbunden und die Struktur wird nicht erfasst. Derzeit wird behauptet, dass PHATE aufgrund der unterschiedlichen Parameter robust ist, aber ich denke, es ist besser, einen der beiden zu reparieren und damit zu experimentieren.

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

def diffusion_operator(data, knn=5, alpha=40):
    #Euklidische Distanzmatrixkonstruktion
    dist = squareform(pdist(data))
    #Berechnung des adaptiven Bandbreiten-Epsilons. Sortieren Sie die Abstandsmatrix auf den knnth-Wert.
    epsilons = np.sort(dist, axis=1)[:, knn]
    # alpha-decaying kernel
    kernel_mtx = np.exp(-1.* np.power(dist / epsilons[:, np.newaxis], alpha))
    # L1-Normalisiert mit Norm. Es wird eine Wahrscheinlichkeitsverteilung für jede Zeile.
    l1norm_kernel_mtx = kernel_mtx / np.abs(kernel_mtx).sum(axis=1)[:, np.newaxis]
    #Symmetrie.
    transition_mtx = 0.5 * l1norm_kernel_mtx + 0.5 * l1norm_kernel_mtx.T
    return transition_mtx

2. Bestimmung der Diffusionszeitskala "t"

Der letzte wichtige Parameter ist "t", der die Zeitskala der Diffusion bestimmt. Entscheiden Sie anhand der Affinität, wie viele Schritte Sie mit dem Random Walk ausführen möchten (= wie oft der Diffusionsoperator angehoben werden soll). Der Diffusionsprozess dient auch dazu, die Daten zu entrauschen. Datenpunkte, die durch kurze Pfade mit hoher Affinität verbunden sind, besuchen sich eher auf einem zufälligen Spaziergang, und Ausreißerpunkte mit nur Pfaden mit niedriger Affinität werden seltener besucht, sodass erstere mit jedem Fortschreiten des Diffusionsschritts häufiger besucht werden. Hervorgehoben werden. Daher kann die Diffusionszeitskala als Grad der Rauschentfernung angesehen werden. Wenn sie zu klein ist, bleibt viel Rauschen zurück, und wenn sie zu groß ist, können wichtige Signale als Rauschen entfernt werden. In PHATE wurde als ein Verfahren zum automatischen Bestimmen der Diffusionszeitskala t aus den Daten das Rauschen vollständig entfernt, indem die "Informationsmenge" untersucht wurde, die in der t-ten Potenz des Diffusionsoperators bei verschiedenen t enthalten ist, während die Originaldaten. Bestimmen Sie t durch Timing, damit die Struktur von nicht gequetscht wird. Insbesondere wird der Eigenwert für jedes (symmetrische Konjugat) des Diffusionsoperators zur t-ten Potenz berechnet, um [von Neumann-Entropie] zu berechnen (https://en.wikipedia.org/wiki/Von_Neumann_entropy) (VNE). Wenn t zunimmt, nimmt VNE ab. Die anfängliche VNE-Reduzierung ist auf die Rauschentfernung zurückzuführen, und nicht wesentliche Eigenwerte gehen schnell auf Null. Die Eigenwerte, die die Struktur der Daten widerspiegeln, nehmen langsamer ab. Wenn Sie also t zum Zeitpunkt des Umschaltens der Reduktionsrate auswählen, können Sie automatisch t auswählen, das Rauschen entfernt und die Struktur nicht zusammenbricht. Führen Sie hier nach der Berechnung des VNE mit verschiedenen t zwei lineare Regressionen auf der linken und rechten Seite jedes t durch, um den "Kniepunkt" des Diagramms zu ermitteln, indem Sie den Punkt mit dem geringsten Gesamtfehler auswählen. ..

from tqdm import tqdm_notebook as tqdm

def von_Neumann_entropy(P):
    # symmetric conjugate (diffusion affinity matrix)Berechnung von
    norm_factor = np.sqrt(P.sum(axis=1))
    A = P / norm_factor[:, np.newaxis] / norm_factor
    #Einzigartige Wertzerlegung
    eig, _ = np.linalg.eigh(A)
    eig = eig[eig > 0]
    eig /= eig.sum()
    #VNE ist eine normalisierte Eigenwert-Shannon-Entropie
    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
    #Für jedes t eine lineare Regression auf der linken und rechten Seite, um den Fehler zu berechnen
    #Es scheint einen besseren Weg zu geben, aber hier berechne ich ehrlich eins nach dem anderen...
    for i, pnt in enumerate(candidates):
        # left area linear regression (y = ax + b)
        #Berechnen Sie die Methode des minimalen Quadrats gemäß dem Lehrbuch auf der linken Seite des Punkts.
        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()
    #Der Punkt mit dem kleinsten Fehler sei der Kniepunkt
    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)):
        #Berechnen Sie die t-te Potenz und VNE des Diffusionsoperators für jedes 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. Führen Sie den Diffusionsprozess durch und berechnen Sie "mögliche Entfernungen"

Im Fall der Diffusionskarte können, wenn der Diffusionsoperator in Eigenwerte zerlegt wird, die Koordinaten erhalten werden, die den Diffusionsabstand (den quadratischen Abstand jedes Punktes der Wahrscheinlichkeitsverteilung nach t Schritten) widerspiegeln, wie er ist. Bei dieser Entfernungsberechnung ist die Empfindlichkeit von Werten mit geringer Wahrscheinlichkeit jedoch gering, so dass es schwierig ist, Informationen über die Entfernung zwischen entfernten Punkten (globale Datenstruktur) wiederzugeben. Daher wird der euklidische Abstand nach logarithmischer Umwandlung des Diffusionsoperators nach t-stufiger Diffusion berechnet. Dies wird als "potentielle Entfernung" bezeichnet. Das Papier erklärt ausführlich, was dieser neu eingeführte Entfernungsindex physikalisch bedeutet (aber ich war mir nicht sicher).

def potential_distances(P, t):
    #multiplizieren Sie den Diffusionsoperator mit t
    Pt = np.linalg.matrix_power(P, t)
    # -log(P^t)Konvertieren und berechnen Sie die euklidische Entfernung
    return squareform(pdist(-1.0 * np.log(Pt + 1e-7)))

4. Niedrigdimensionale Einbettung des potenziellen Abstands mit MDS

Anstelle einer Eigenwertzerlegung wie der Diffusionskarte wird die in 3 berechnete potenzielle Entfernung in eine niedrigere Dimension verschoben, die für die Visualisierung durch MDS (mehrdimensionale Skalenkonstruktionsmethode) geeignet ist. Zunächst werden niedrigdimensionale Koordinaten unter Verwendung des klassischen MDS (= PCoA) berechnet, und das metrische MDS wird unter Verwendung dieses als Anfangswert berechnet. Das metrische MDS wird vom SMACOF-Algorithmus von scikit-learn ausgeführt.

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

Lassen Sie uns die obigen Schritte tatsächlich mit Baumstrukturdaten ausführen.

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

Es wurden ähnliche Ergebnisse wie bei der ursprünglichen Implementierung erzielt. Tatsächlich wird der obige Algorithmus selbst nicht berechnet, indem Wege entwickelt werden, um die Stabilität der numerischen Berechnung bei jedem Schritt zu verbessern, oder indem verschiedene Maßnahmen getroffen werden, um Speicher zu sparen und mit hoher Geschwindigkeit zu berechnen, sondern er ist ungefähr. Es scheint, dass es Berechnungen durchführt, aber es sieht in der Regel so aus.

[PYTHON] Versuchen Sie es mit PHATE, einer Methode zur Reduzierung und Visualisierung biologischer Daten