Hydra t-SNE Distortion Analysis Revisited¶

This notebook revisits the hydra example using exactly the same setup as those in scDEED's Figure 2.

In [1]:

import pandas as pd
import numpy as np
import anndata as ad
import scanpy as sc
from distortions.geometry import Geometry, local_distortions, bind_metric, neighborhoods
from distortions.visualization.interactive import dplot
from sklearn.neighbors import NearestNeighbors
import warnings
warnings.filterwarnings('ignore')

import pandas as pd import numpy as np import anndata as ad import scanpy as sc from distortions.geometry import Geometry, local_distortions, bind_metric, neighborhoods from distortions.visualization.interactive import dplot from sklearn.neighbors import NearestNeighbors import warnings warnings.filterwarnings('ignore')

Configuration¶

All data comes from the two CSV files exported by run_hydra_tsne.R. The cell_type column is the exact cluster.names mapping from the paper's Fig2.R.

In [ ]:

BASE = '/Users/krissankaran/Desktop/collaborations/manifold_learning' # modify as needed
CSV40  = f'{BASE}/scdeed/repro/output/hydra_tsne_perplexity40.csv'
CSV230 = f'{BASE}/scdeed/repro/output/hydra_tsne_perplexity230.csv'
OUTDIR = f'{BASE}/distortion_analysis'

BASE = '/Users/krissankaran/Desktop/collaborations/manifold_learning' # modify as needed CSV40 = f'{BASE}/scdeed/repro/output/hydra_tsne_perplexity40.csv' CSV230 = f'{BASE}/scdeed/repro/output/hydra_tsne_perplexity230.csv' OUTDIR = f'{BASE}/distortion_analysis'

build_edge_links¶

Build neighborhood connections from the embedding object.

In [3]:

def build_edge_links(embedding_df, X_high, perplexity):
    """Build neighbourhood edge links from an embedding DataFrame.

    The DataFrame must have columns ``tSNE_1``, ``tSNE_2``, ``cell_type``
    (the last coming from the paper's ``cluster.names`` vector).
    """
    emb = embedding_df[['tSNE_1','tSNE_2']].values
    ctypes = embedding_df['cell_type'].values
    nn = perplexity

    adata = ad.AnnData(X=X_high)
    adata.obsm['X_tSNE'] = emb
    adata.obs['cell_type'] = ctypes
    sc.pp.neighbors(adata, n_neighbors=nn, n_pcs=5, method='gauss')
    N = neighborhoods(adata, threshold=0.2, outlier_factor=3, embed_key='X_tSNE')
    print(f'neighborhoods() -> {len(N)} connections')
    return N

def build_edge_links(embedding_df, X_high, perplexity): """Build neighbourhood edge links from an embedding DataFrame. The DataFrame must have columns ``tSNE_1``, ``tSNE_2``, ``cell_type`` (the last coming from the paper's ``cluster.names`` vector). """ emb = embedding_df[['tSNE_1','tSNE_2']].values ctypes = embedding_df['cell_type'].values nn = perplexity adata = ad.AnnData(X=X_high) adata.obsm['X_tSNE'] = emb adata.obs['cell_type'] = ctypes sc.pp.neighbors(adata, n_neighbors=nn, n_pcs=5, method='gauss') N = neighborhoods(adata, threshold=0.2, outlier_factor=3, embed_key='X_tSNE') print(f'neighborhoods() -> {len(N)} connections') return N

1. Load data & compute distortions¶

In [4]:

def _radius(emb, nn):
    nbrs = NearestNeighbors(n_neighbors=nn).fit(emb)
    d, _ = nbrs.kneighbors(emb)
    return 3 * float(np.mean(d[:, 1:]))

results = {}  # perplexity -> {'metric': DataFrame, 'neighbors': DataFrame}

for perp, csv_path in [(40, CSV40), (230, CSV230)]:
    df = pd.read_csv(csv_path).sample(frac=0.1)
    X_high = df[['PC_1','PC_2','PC_3','PC_4','PC_5']].values
    emb = df[['tSNE_1','tSNE_2']].values
    nn = perp

    # geometry & distortions
    geom = Geometry('brute', laplacian_method='geometric',
                   laplacian_kwds={'scaling_epps':5},
                   adjacency_kwds={'n_neighbors':nn},
                   affinity_kwds={'radius':_radius(emb, nn)})
    H, Hvv, Hs = local_distortions(emb, X_high, geom)
    Hs[Hs > 2.5] = 2.5
    Hs /= Hs.mean()
    for i in range(len(H)):
        H[i] = Hvv[i] @ np.diag(Hs[i]) @ Hvv[i].T

    df_metric = bind_metric(emb, Hvv, Hs)
    df_metric['cell_type'] = df['cell_type'].values   # from R script's cluster.names

    # edge links
    N = build_edge_links(df, X_high, perp)

    results[perp] = {'metric': df_metric, 'neighbors': N}
    print(f'  Perplexity {perp}: mean Hs = {Hs.mean():.3f}, {len(N)} edges')

def _radius(emb, nn): nbrs = NearestNeighbors(n_neighbors=nn).fit(emb) d, _ = nbrs.kneighbors(emb) return 3 * float(np.mean(d[:, 1:])) results = {} # perplexity -> {'metric': DataFrame, 'neighbors': DataFrame} for perp, csv_path in [(40, CSV40), (230, CSV230)]: df = pd.read_csv(csv_path).sample(frac=0.1) X_high = df[['PC_1','PC_2','PC_3','PC_4','PC_5']].values emb = df[['tSNE_1','tSNE_2']].values nn = perp # geometry & distortions geom = Geometry('brute', laplacian_method='geometric', laplacian_kwds={'scaling_epps':5}, adjacency_kwds={'n_neighbors':nn}, affinity_kwds={'radius':_radius(emb, nn)}) H, Hvv, Hs = local_distortions(emb, X_high, geom) Hs[Hs > 2.5] = 2.5 Hs /= Hs.mean() for i in range(len(H)): H[i] = Hvv[i] @ np.diag(Hs[i]) @ Hvv[i].T df_metric = bind_metric(emb, Hvv, Hs) df_metric['cell_type'] = df['cell_type'].values # from R script's cluster.names # edge links N = build_edge_links(df, X_high, perp) results[perp] = {'metric': df_metric, 'neighbors': N} print(f' Perplexity {perp}: mean Hs = {Hs.mean():.3f}, {len(N)} edges')

neighborhoods() -> 133 connections
  Perplexity 40: mean Hs = 1.000, 133 edges
neighborhoods() -> 27 connections
  Perplexity 230: mean Hs = 1.000, 27 edges

3. Interactive plot — Perplexity 40¶

Hover over an edge to see the true-vs-embedding distance ratio.

In [5]:

color_mapping = {
   'stem cell': 'Stem Cells',
   'stem cell/progenitor': 'Stem Cells',
   'nb1': 'Nematoblasts', 'nb2': 'Nematoblasts', 'nb3': 'Nematoblasts',
   'nb4': 'Nematoblasts', 'nb5': 'Nematoblasts', 'nematocyte': 'Nematocytes',
    'neuron ec1': 'Neurons', 'neuron ec2': 'Neurons', 'neuron ec3': 'Neurons',
    'neuron ec4': 'Neurons', 'neuron ec5': 'Neurons', 'neuron en1': 'Neurons',
    'neuron en2': 'Neurons', 'neuron en3': 'Neurons', 'spumous mucous gland cell':
    'Gland Cells', 'granular mucous gland cell': 'Gland Cells',
    'zymogen gland cell': 'Gland Cells', 'male germline': 'Germline',
    'female germline 1': 'Germline', 'female germline 2 nurse cell': 'Germline',
    'ecEp-nem(id)': 'Ectodermal Ep.', 'ecEp-nb(pd)': 'Ectodermal Ep.',
    'battery cell 1(mp) ': 'Ectodermal Ep.', 'battery cell 2(mp) ': 'Ectodermal Ep.',
    'enEp tent-nem(pd)': 'Endodermal Ep.', 'enEp-nem(pd)': 'Endodermal Ep.',
    'enEp-nb(pd)': 'Endodermal Ep.', 'enEp-tent-nem(pd)': 'Endodermal Ep.',
    'neuron/gland cell progenitor': 'Progenitors',
    'neuron progenitor': 'Progenitors', 'head': 'Misc/Structural', 'tentacle': 'Misc/Structural',
    'foot': 'Misc/Structural', 'basal disk': 'Misc/Structural', 'db': 'Misc/Structural'
}
# apply mapping to create a new column
results[230]['metric']['group'] = results[230]['metric']['cell_type'].map(color_mapping)
results[40]['metric']['group'] = results[40]['metric']['cell_type'].map(color_mapping)

group_order = ['Stem Cells', 'Nematoblasts', 'Nematocytes', 'Gland Cells', 'Neurons',
               'Germline', 'Ectodermal Ep.', 'Endodermal Ep.', 'Progenitors', 'Misc/Structural']

for p in (40, 230):
    df = results[p]['metric'].copy()
    df['group'] = pd.Categorical(df['group'], categories=group_order, ordered=True)
    df['_orig_pos'] = np.arange(len(df))          # position N's indices refer to
    df_sorted = df.sort_values('group', kind='stable').reset_index(drop=True)

    old_to_new = dict(zip(df_sorted['_orig_pos'], df_sorted.index))
    results[p]['metric'] = df_sorted.drop(columns='_orig_pos')

    N_old = results[p]['neighbors']
    results[p]['neighbors'] = {
        old_to_new[center]: [old_to_new[nb] for nb in neighbors]
        for center, neighbors in N_old.items()
    }

color_mapping = { 'stem cell': 'Stem Cells', 'stem cell/progenitor': 'Stem Cells', 'nb1': 'Nematoblasts', 'nb2': 'Nematoblasts', 'nb3': 'Nematoblasts', 'nb4': 'Nematoblasts', 'nb5': 'Nematoblasts', 'nematocyte': 'Nematocytes', 'neuron ec1': 'Neurons', 'neuron ec2': 'Neurons', 'neuron ec3': 'Neurons', 'neuron ec4': 'Neurons', 'neuron ec5': 'Neurons', 'neuron en1': 'Neurons', 'neuron en2': 'Neurons', 'neuron en3': 'Neurons', 'spumous mucous gland cell': 'Gland Cells', 'granular mucous gland cell': 'Gland Cells', 'zymogen gland cell': 'Gland Cells', 'male germline': 'Germline', 'female germline 1': 'Germline', 'female germline 2 nurse cell': 'Germline', 'ecEp-nem(id)': 'Ectodermal Ep.', 'ecEp-nb(pd)': 'Ectodermal Ep.', 'battery cell 1(mp) ': 'Ectodermal Ep.', 'battery cell 2(mp) ': 'Ectodermal Ep.', 'enEp tent-nem(pd)': 'Endodermal Ep.', 'enEp-nem(pd)': 'Endodermal Ep.', 'enEp-nb(pd)': 'Endodermal Ep.', 'enEp-tent-nem(pd)': 'Endodermal Ep.', 'neuron/gland cell progenitor': 'Progenitors', 'neuron progenitor': 'Progenitors', 'head': 'Misc/Structural', 'tentacle': 'Misc/Structural', 'foot': 'Misc/Structural', 'basal disk': 'Misc/Structural', 'db': 'Misc/Structural' } # apply mapping to create a new column results[230]['metric']['group'] = results[230]['metric']['cell_type'].map(color_mapping) results[40]['metric']['group'] = results[40]['metric']['cell_type'].map(color_mapping) group_order = ['Stem Cells', 'Nematoblasts', 'Nematocytes', 'Gland Cells', 'Neurons', 'Germline', 'Ectodermal Ep.', 'Endodermal Ep.', 'Progenitors', 'Misc/Structural'] for p in (40, 230): df = results[p]['metric'].copy() df['group'] = pd.Categorical(df['group'], categories=group_order, ordered=True) df['_orig_pos'] = np.arange(len(df)) # position N's indices refer to df_sorted = df.sort_values('group', kind='stable').reset_index(drop=True) old_to_new = dict(zip(df_sorted['_orig_pos'], df_sorted.index)) results[p]['metric'] = df_sorted.drop(columns='_orig_pos') N_old = results[p]['neighbors'] results[p]['neighbors'] = { old_to_new[center]: [old_to_new[nb] for nb in neighbors] for center, neighbors in N_old.items() }

In [30]:

pal = ['#8C1F28', '#D91438', '#F2AEBB', '#A6375F', '#4BA695', "#135832", '#F2CF66', '#D9AE79', '#F2E0D5', '#F25A38']

p40e = (dplot(results[40]['metric'], height=500, width=600)
         .mapping(x='embedding_0', y='embedding_1', color='group')
         .geom_ellipse(radiusMin=0.01, radiusMax =  3)
         .inter_edge_link(N=results[40]['neighbors'], backgroundOpacity=0.1, strokeWidth=0.1, strokeOpacity=0.1, highlightStrokeWidth=1, highlightColor="#000000", stroke="#000000")
         .scale_color(scheme=pal)
         .labs(title='Hydra t-SNE – Perplexity 40'))
p40e

pal = ['#8C1F28', '#D91438', '#F2AEBB', '#A6375F', '#4BA695', "#135832", '#F2CF66', '#D9AE79', '#F2E0D5', '#F25A38'] p40e = (dplot(results[40]['metric'], height=500, width=600) .mapping(x='embedding_0', y='embedding_1', color='group') .geom_ellipse(radiusMin=0.01, radiusMax = 3) .inter_edge_link(N=results[40]['neighbors'], backgroundOpacity=0.1, strokeWidth=0.1, strokeOpacity=0.1, highlightStrokeWidth=1, highlightColor="#000000", stroke="#000000") .scale_color(scheme=pal) .labs(title='Hydra t-SNE – Perplexity 40')) p40e

Out[30]:

dplot(dataset=[{'embedding_0': 25.290912945276, 'embedding_1': -16.3388301634052, 'x0': -0.8209345032702061, '…

In [ ]:

p40e.save(f'hydra_perplexity40_distortion.svg')

p40e.save(f'hydra_perplexity40_distortion.svg')

5. Interactive plot — Perplexity 230¶

In [20]:

p230e = (dplot(results[230]['metric'], aspect_ratio=1.0, height=500, width=600)
         .mapping(x='embedding_0', y='embedding_1', color='group')
         .geom_ellipse(radiusMin=1, radiusMax =  4)
         .inter_edge_link(N=results[230]['neighbors'], backgroundOpacity=0.1, strokeWidth=0.1, strokeOpacity=.1, highlightStrokeWidth=1, highlightColor="#000000", stroke="#000000")
         .scale_color(scheme=pal)
         .labs(title='Hydra t-SNE – Perplexity 230'))

p230e

p230e = (dplot(results[230]['metric'], aspect_ratio=1.0, height=500, width=600) .mapping(x='embedding_0', y='embedding_1', color='group') .geom_ellipse(radiusMin=1, radiusMax = 4) .inter_edge_link(N=results[230]['neighbors'], backgroundOpacity=0.1, strokeWidth=0.1, strokeOpacity=.1, highlightStrokeWidth=1, highlightColor="#000000", stroke="#000000") .scale_color(scheme=pal) .labs(title='Hydra t-SNE – Perplexity 230')) p230e

Out[20]:

dplot(dataset=[{'embedding_0': 29.4832281806549, 'embedding_1': 6.08555783581568, 'x0': 0.3629299071214185, 'y…

In [ ]:

p230e.save(f'hydra_perplexity230_distortion_hover_neighbor.svg')

p230e.save(f'hydra_perplexity230_distortion_hover_neighbor.svg')

In [ ]:

df = results[230]['metric']
unmapped = df[df['group'].isna()]['cell_type'].unique()
if len(unmapped) > 0:
    print(f"Warning: The following types were not mapped: {unmapped}")

df = results[230]['metric'] unmapped = df[df['group'].isna()]['cell_type'].unique() if len(unmapped) > 0: print(f"Warning: The following types were not mapped: {unmapped}")

6. (Optional) Save SVGs and CSVs to disk¶

!mkdir -p "{OUTDIR}" p40.save(f'{OUTDIR}/hydra_perplexity40_distortion.svg') p40e.save(f'{OUTDIR}/hydra_perplexity40_distortion_edges.svg') p230.save(f'{OUTDIR}/hydra_perplexity230_distortion.svg') p230e.save(f'{OUTDIR}/hydra_perplexity230_distortion_edges.svg') results[40]['metric'].to_csv(f'{OUTDIR}/hydra_perplexity40_with_distortion.csv', index=False) results[230]['metric'].to_csv(f'{OUTDIR}/hydra_perplexity230_with_distortion.csv', index=False) results[40]['neighbors'].to_csv(f'{OUTDIR}/hydra_perplexity40_neighborhoods.csv', index=False) results[230]['neighbors'].to_csv(f'{OUTDIR}/hydra_perplexity230_neighborhoods.csv', index=False) print('All outputs saved.')