A multimodal approach for visualizing and identifying electrophysiological cell types in vivo.

Code ↔ Paper

8 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 8 matches

[1] § Methods › Weighted nearest neighbor (WNN) algorithm › β-variational autoencoder and multi-layer perceptron identification architectures ↔ PairOfVAEs_cuda.ipynb, lines 1020–1117 · score 0.96 · Adam optimizer, MLPClassifier, ReLU, hidden layer, latent dim, latent embedding
[2] § Methods › Weighted nearest neighbor (WNN) algorithm › Identification comparison ↔ PairOfVAEs_cuda.ipynb, lines 1119–1226 · score 0.64 · cross validation, confusion matrices, iteration, split, transformation, VAE
[3] § Methods › Weighted nearest neighbor (WNN) algorithm › Distributional shift in waveform shape ↔ juxtacellular/processJianing.ipynb, lines 795–870 · score 0.62 · max_depth, n_estimators, classifier accuracy, trained, cell
[4] § Results › Interpretable representations from PhysMAP enable biological insight ↔ juxtacellular/processJianing.ipynb, lines 127–165 · score 0.61 · inter spike interval, peri stimulus, histograms, latency, onset, mouse
[5] § Methods › Statistics & reproducibility › Extracellular Mouse A1 Dataset ↔ InvivoA1/processSantiago.r, lines 1–71 · score 0.56 · spike quality, ISI violation, amplitude, waveform, cells
[6] § Methods › Weighted nearest neighbor (WNN) algorithm › Identification algorithm ↔ juxtacellular/processJianing.ipynb, lines 795–870 · score 0.55 · GBM classifier, balanced accuracy, split, classes, train, embedded
[7] § Methods › Statistics & reproducibility › Visual Behavior dataset (Allen Institute) ↔ InvivoA1/processSantiago.r, lines 1–71 · score 0.52 · ISI violation, spike waveform, amplitude, quality
[8] § Methods › Weighted nearest neighbor (WNN) algorithm › Joint embedding and label transfer of Ultras and Visual Behavior ↔ juxtacellular/processJianing.r, lines 152–199 · score 0.50 · log ratio, UMAP embedding, identities, PSTH, waveform, cell

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 1,152 lines · 37 KB · no license · 3 matches

# %% [markdown]
# # PhysMAP: Process Jianing Data
#
# This notebook ports the R analysis from `processJianing.r` to Python using Scanpy and Muon.
#
# **Original R pipeline:**
# - Seurat v4 for single-cell analysis
# - Multi-modal integration via FindMultiModalNeighbors (WNN)
#
# **Python implementation:**
# - Scanpy for single-modality processing
# - Muon for multi-modal data handling and integration
# %% [markdown]
# ## 1. Imports and Configuration
# %%
import scanpy as sc
import muon as mu
import numpy as np
import pandas as pd
from scipy.io import loadmat
import matplotlib.pyplot as plt
import seaborn as sns
from anndata import AnnData
from muon import MuData
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import warnings
warnings.filterwarnings('ignore')
# Set random seed for reproducibility
np.random.seed(42)
# Configure scanpy
sc.settings.verbosity = 1
sc.settings.set_figure_params(dpi=100, frameon=False, figsize=(6, 5))
# %%
# Constants (from constants.R)
RESOLUTION = 2
UMAP_SEED = 42
UMAP_NEIGHBORS = 20
UMAP_MIN_DIST = 0.2
UMAP_METRIC = "cosine"
UMAP_COMPONENTS = 10
N_PCS = 100
# %% [markdown]
# ## 2. Data Loading
#
# Load the MATLAB data file and create a MuData object with multiple modalities.
# %%
def clr_normalize(X):
"""
Centered Log-Ratio (CLR) normalization per cell (row).
Equivalent to Seurat's NormalizeData(method='CLR', margin=4).
For each cell, computes: log(x + 1) - mean(log(x + 1))
This centers the log-transformed data around zero for each cell.
"""
# Add pseudocount and log transform
X_pseudo = X + 1
log_X = np.log(X_pseudo)
# Compute geometric mean per cell (row) and center
geometric_mean = np.mean(log_X, axis=1, keepdims=True)
X_clr = log_X - geometric_mean
return X_clr
# %%
def extract_matlab_strings(arr):
"""
Extract strings from MATLAB cell arrays loaded by scipy.io.loadmat.
MATLAB cell arrays come as nested numpy arrays of objects.
"""
arr = np.asarray(arr)
# Flatten if needed
if arr.ndim > 1:
arr = arr.flatten()
# If it's a single element containing the whole array as string, something is wrong
# This happens when squeeze_me=True mangles cell arrays
if len(arr) == 1 and isinstance(arr[0], str) and '[' in arr[0]:
# This is a string representation of the array, not the actual array
raise ValueError("MATLAB cell array was converted to string representation. Reload with squeeze_me=False")
result = []
for item in arr:
if isinstance(item, np.ndarray):
# Nested array - extract the string
if item.size > 0:
val = item.flat[0]
result.append(str(val).strip())
else:
result.append('')
elif isinstance(item, str):
result.append(item.strip())
else:
result.append(str(item).strip())
return np.array(result)
def load_rda_file(filepath):
"""
Load an R .Rda file and return the first object as a numpy array.
"""
try:
import pyreadr
result = pyreadr.read_r(filepath)
# Get the first (and usually only) object
key = list(result.keys())[0]
return result[key].values.flatten()
except ImportError:
print(f"Warning: pyreadr not installed. Cannot load {filepath}")
print("Install with: pip install pyreadr")
return None
except Exception as e:
print(f"Warning: Could not load {filepath}: {e}")
return None
def load_jianing_data(mat_file_path):
"""
Load the Jianing mouse S1 dataset from MATLAB file.
Creates a MuData object with the following modalities:
- ISI: Inter-spike interval data (100 features)
- WF: Waveform data (351 features)
- PSTH: Peri-stimulus time histogram (151 features)
- features: Hand-crafted features (7 features)
Returns:
mdata: MuData object with all modalities
metadata_df: DataFrame with cell metadata
"""
import os
# Load MATLAB file - use squeeze_me=False to preserve array structure
mat_data = loadmat(mat_file_path, squeeze_me=False, struct_as_record=False)
# Extract dimensions
depth_arr = np.asarray(mat_data['depth']).flatten()
n_cells = len(depth_arr)
cell_ids = [f"cell_{i}" for i in range(n_cells)]
# Extract continuous variables
depth = depth_arr.astype(np.float64)
latency = np.asarray(mat_data['latency']).flatten().astype(np.float64)
latency = np.nan_to_num(latency, nan=0.0)
onset_rate = np.asarray(mat_data['onset_rate']).flatten().astype(np.float64)
onset_rate = np.nan_to_num(onset_rate, nan=0.0)
# Extract categorical variables (MATLAB cell arrays need special handling)
layer = extract_matlab_strings(mat_data['layer'])
cell_type = extract_matlab_strings(mat_data['true_cell_type'])
# Create layerCellType (combined label)
layer_cell_type = np.array([f"{ct}-{l}" for ct, l in zip(cell_type, layer)])
# Extract feature matrices
ISI = np.asarray(mat_data['ISI']).astype(np.float64)
WF = np.asarray(mat_data['waveformV']).astype(np.float64)
PSTH = np.asarray(mat_data['AllResp']).astype(np.float64)
# Ensure correct orientation (cells x features)
if ISI.shape[0] != n_cells:
ISI = ISI.T
if WF.shape[0] != n_cells:
WF = WF.T
if PSTH.shape[0] != n_cells:
PSTH = PSTH.T
# Extract additional features from 'features' matrix
features_mat = np.asarray(mat_data['features']).astype(np.float64)
if features_mat.shape[0] != n_cells:
features_mat = features_mat.T
features_mat = np.nan_to_num(features_mat, nan=0.0)
# Load width and ratio_p2t from .Rda files (like the R code does)
base_dir = os.path.dirname(mat_file_path)
width_path = os.path.join(base_dir, 'width.Rda')
width = load_rda_file(width_path)
if width is None:
print("Warning: Could not load width.Rda, using zeros")
width = np.zeros(n_cells)
else:
width = width.astype(np.float64)
print(f"Loaded width from .Rda: min={width.min():.4f}, max={width.max():.4f}")
ratio_p2t_path = os.path.join(base_dir, 'ratio_p2t.Rda')
ratio_p2t = load_rda_file(ratio_p2t_path)
if ratio_p2t is None:
print("Warning: Could not load ratio_p2t.Rda, using zeros")
ratio_p2t = np.zeros(n_cells)
else:
ratio_p2t = ratio_p2t.astype(np.float64)
print(f"Loaded ratio_p2t from .Rda: min={ratio_p2t.min():.4f}, max={ratio_p2t.max():.4f}")
# Create metadata DataFrame
metadata_df = pd.DataFrame({
'Layer': layer,
'CellType': cell_type,
'layerCellType': layer_cell_type,
'Depth': depth,
'latency': latency,
'onset': onset_rate,
'width': width,
'ratio_p2t': ratio_p2t
}, index=cell_ids)
# Create AnnData objects for each modality
adata_isi = AnnData(
X=ISI.astype(np.float32),
obs=metadata_df.copy()
)
adata_isi.var_names = [f"ISI_{i}" for i in range(ISI.shape[1])]
adata_wf = AnnData(
X=WF.astype(np.float32),
obs=metadata_df.copy()
)
adata_wf.var_names = [f"WF_{i}" for i in range(WF.shape[1])]
adata_psth = AnnData(
X=PSTH.astype(np.float32),
obs=metadata_df.copy()
)
adata_psth.var_names = [f"PSTH_{i}" for i in range(PSTH.shape[1])]
# Create features AnnData (hand-crafted features)
feature_names = ['Depth', 'OnsetRate', 'Latency', 'StimFR', 'StimFF']
if features_mat.shape[1] >= 7:
feature_names.extend(['width', 'ratio_p2t'])
else:
# Add width and ratio_p2t to features matrix
features_mat = np.column_stack([features_mat, width, ratio_p2t])
feature_names.extend(['width', 'ratio_p2t'])
adata_features = AnnData(
X=features_mat.astype(np.float32),
obs=metadata_df.copy()
)
adata_features.var_names = feature_names[:features_mat.shape[1]]
# Create concatenated modality (WF + ISI + PSTH)
concat_data = np.concatenate([WF, ISI, PSTH], axis=1)
adata_concat = AnnData(
X=concat_data.astype(np.float32),
obs=metadata_df.copy()
)
concat_var_names = ([f"WF_{i}" for i in range(WF.shape[1])] +
[f"ISI_{i}" for i in range(ISI.shape[1])] +
[f"PSTH_{i}" for i in range(PSTH.shape[1])])
adata_concat.var_names = concat_var_names
# Create MuData object
mdata = MuData({
'ISI': adata_isi,
'WF': adata_wf,
'PSTH': adata_psth,
'features': adata_features,
'concat': adata_concat
})
print(f"\nLoaded {n_cells} cells")
print(f"Modalities:")
print(f" - ISI: {adata_isi.shape}")
print(f" - WF: {adata_wf.shape}")
print(f" - PSTH: {adata_psth.shape}")
print(f" - features: {adata_features.shape}")
print(f" - concat: {adata_concat.shape}")
print(f"\nCell types: {np.unique(cell_type)}")
print(f"Layers: {np.unique(layer)}")
return mdata, metadata_df
# %%
# Load the data
mat_file = "JianingData/MergedData.mat"
mdata, metadata_df = load_jianing_data(mat_file)
# %%
# Inspect the MuData structure
print(mdata)
# %% [markdown]
# ## 3. Single-Modality Processing
#
# Define a function to process each modality individually, equivalent to `calcRepresentation()` in the R code.
# %%
def process_modality(adata, modality_name, n_pcs=30, normalize=True,
n_neighbors=20, metric='cosine', random_state=42):
"""
Process a single modality through the full Seurat-equivalent pipeline.
Steps:
1. CLR normalization (optional)
2. Find variable features (using all features for electrophysiology data)
3. Scale data (z-score normalization)
4. PCA
5. Find neighbors
6. UMAP (2D)
7. Leiden clustering
Args:
adata: AnnData object
modality_name: Name of the modality (for storing results)
n_pcs: Number of principal components
normalize: Whether to apply CLR normalization
n_neighbors: Number of neighbors for kNN graph
metric: Distance metric for neighbors
random_state: Random seed for reproducibility
Returns:
adata: Processed AnnData object with embeddings and clusters
"""
adata = adata.copy()
# Store raw data
adata.layers['raw'] = adata.X.copy()
# 1. CLR normalization (if requested)
if normalize:
adata.X = clr_normalize(adata.X)
adata.layers['normalized'] = adata.X.copy()
# 2. For electrophysiology data, use all features as variable
# (unlike scRNA-seq where we select highly variable genes)
adata.var['highly_variable'] = True
# 3. Scale data (z-score normalization)
sc.pp.scale(adata, max_value=None)
# 4. PCA
n_pcs_actual = min(n_pcs, min(adata.n_obs, adata.n_vars) - 1)
sc.tl.pca(adata, n_comps=n_pcs_actual, svd_solver='arpack', random_state=random_state)
# 5. Find neighbors
sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs_actual,
metric=metric, random_state=random_state)
# 6. UMAP
sc.tl.umap(adata, random_state=random_state, min_dist=UMAP_MIN_DIST)
# 7. Leiden clustering
sc.tl.leiden(adata, resolution=RESOLUTION, random_state=random_state)
# Rename clustering column for this modality
adata.obs[f'{modality_name}_leiden'] = adata.obs['leiden']
print(f"Processed {modality_name}:")
print(f" - PCA components: {n_pcs_actual}")
print(f" - Clusters found: {adata.obs['leiden'].nunique()}")
return adata
# %% [markdown]
# ## 4. Process Each Modality
# %%
# Process ISI modality (with CLR normalization)
adata_isi = process_modality(
mdata['ISI'],
'ISI',
n_pcs=N_PCS,
normalize=True,
metric=UMAP_METRIC
)
# %%
# Process WF modality (without normalization, per R code)
adata_wf = process_modality(
mdata['WF'],
'WF',
n_pcs=N_PCS,
normalize=False, # WF is not normalized in the R code
metric=UMAP_METRIC
)
# %%
# Process PSTH modality (with CLR normalization)
adata_psth = process_modality(
mdata['PSTH'],
'PSTH',
n_pcs=N_PCS,
normalize=True,
metric=UMAP_METRIC
)
# %%
# Process concatenated modality (with CLR normalization)
# adata_concat = process_modality(
# mdata['concat'],
# 'concat',
# n_pcs=N_PCS,
# normalize=True,
# metric=UMAP_METRIC
# )
# %% [markdown]
# ## 5. Individual Modality Visualizations
# %%
# Create UMAP plots for each modality colored by layerCellType
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
modalities = [
(adata_wf, 'WF'),
(adata_isi, 'ISI'),
(adata_psth, 'PSTH'),
]
for ax, (adata, name) in zip(axes, modalities):
sc.pl.umap(adata, color='layerCellType', ax=ax, show=False,
title=name, legend_loc='none', frameon=False, size=100)
plt.tight_layout()
plt.show()
# %%
# Create UMAP plots colored by cluster ID for each modality
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
modalities_for_clustering = [
(adata_wf, 'WF', 'WF_leiden'),
(adata_isi, 'ISI', 'ISI_leiden'),
(adata_psth, 'PSTH', 'PSTH_leiden')
]
for ax, (adata, name, cluster_key) in zip(axes, modalities_for_clustering):
sc.pl.umap(adata, color=cluster_key, ax=ax, show=False,
title=name, legend_loc='on data', frameon=False,size=100)
plt.tight_layout()
plt.show()
# %% [markdown]
# ## 6. Multi-Modal Integration
#
# Integrate the three modalities (WF, ISI, PSTH) using a weighted nearest neighbors approach similar to Seurat's FindMultiModalNeighbors.
# %%
# Update the MuData object with processed modalities
mdata.mod['ISI'] = adata_isi
mdata.mod['WF'] = adata_wf
mdata.mod['PSTH'] = adata_psth
# mdata.mod['concat'] = adata_concat
# Update the MuData object
mdata.update()
# %%
# Muon's multi-modal weighted nearest neighbors integration
# This computes cell-specific weights for each modality based on how well
# each modality captures the local neighborhood structure
# First, ensure each modality has its PCA stored with a unique key
# We need to copy the PCA embeddings to the MuData level
for mod_name in ['WF', 'ISI', 'PSTH']:
mod = mdata.mod[mod_name]
# Store PCA in obsm with modality-specific key
key = f'X_pca_{mod_name}'
mdata.obsm[key] = mod.obsm['X_pca'][:, :N_PCS]
print(f"Stored {key}: {mdata.obsm[key].shape}")
print(f"\nAvailable representations in mdata.obsm: {list(mdata.obsm.keys())}")
# %%
# Compute multi-modal weighted nearest neighbors using Muon
# This is similar to Seurat's FindMultiModalNeighbors
# We only integrate WF, ISI, and PSTH (not features or concat)
print("Computing multi-modal weighted nearest neighbors...")
# Create a subset MuData with only the modalities we want to integrate
mdata_wnn = MuData({
'WF': mdata.mod['WF'],
'ISI': mdata.mod['ISI'],
'PSTH': mdata.mod['PSTH']
})
# Run Muon's WNN on the subset
mu.pp.neighbors(
mdata_wnn,
key_added='wnn',
n_neighbors=UMAP_NEIGHBORS,
n_bandwidth_neighbors=20,
n_multineighbors=200,
metric='euclidean',
low_memory=False
)
print("WNN graph computed successfully!")
# %%
# Compute UMAP on the WNN graph
mu.tl.umap(mdata_wnn, neighbors_key='wnn', random_state=UMAP_SEED, min_dist=UMAP_MIN_DIST)
# Copy the WNN UMAP back to the main mdata object
mdata.obsm['X_wnn_umap'] = mdata_wnn.obsm['X_umap'].copy()
print(f"WNN UMAP computed: {mdata.obsm['X_wnn_umap'].shape}")
# Check if modality weights were computed
weight_keys = [k for k in mdata_wnn.obs.columns if 'weight' in k.lower()]
if weight_keys:
print(f"\nModality weights available: {weight_keys}")
# Copy weights to main mdata
for key in weight_keys:
mdata.obs[key] = mdata_wnn.obs[key].values
# %% [markdown]
# ## 7. Final Clustering on Integrated Data
# %%
# Run Leiden clustering on the WNN neighbor graph
# The WNN graph is stored at the MuData level in obsp, so we use scanpy directly
# First, let's see what's available
print("Available in mdata_wnn.obsp:", list(mdata_wnn.obsp.keys()))
print("Available in mdata_wnn.uns:", list(mdata_wnn.uns.keys()))
# The WNN connectivities should be in obsp with a key like 'wnn_connectivities'
# We need to run leiden using scanpy on an AnnData that has the WNN graph
# Create a temporary AnnData to run leiden on
import anndata
adata_for_clustering = anndata.AnnData(
X=np.zeros((mdata_wnn.n_obs, 1)), # Dummy X matrix
obs=mdata_wnn.obs.copy()
)
# Copy the WNN graph to the temporary AnnData
connectivities_key = 'wnn_connectivities' if 'wnn_connectivities' in mdata_wnn.obsp else 'connectivities'
distances_key = 'wnn_distances' if 'wnn_distances' in mdata_wnn.obsp else 'distances'
adata_for_clustering.obsp['connectivities'] = mdata_wnn.obsp[connectivities_key]
if distances_key in mdata_wnn.obsp:
adata_for_clustering.obsp['distances'] = mdata_wnn.obsp[distances_key]
# Add the neighbors info to uns
adata_for_clustering.uns['neighbors'] = {
'connectivities_key': 'connectivities',
'distances_key': 'distances',
'params': {'n_neighbors': UMAP_NEIGHBORS, 'method': 'wnn'}
}
# Run Leiden clustering
sc.tl.leiden(adata_for_clustering, resolution=RESOLUTION, random_state=UMAP_SEED)
# Copy the clustering results back to main mdata
mdata.obs['wnn_leiden'] = adata_for_clustering.obs['leiden'].values
mdata_wnn.obs['leiden'] = adata_for_clustering.obs['leiden'].values
print(f"Found {mdata.obs['wnn_leiden'].nunique()} clusters at resolution {RESOLUTION}")
# %% [markdown]
# ## 8. Final Visualizations
# %%
# Prepare data for plotting
# Get UMAP coordinates
umap_df = pd.DataFrame(
mdata.obsm['X_wnn_umap'],
columns=['UMAP1', 'UMAP2'],
index=mdata.obs_names
)
# Get metadata from one of the modalities (they all have the same metadata)
mod_obs = mdata.mod['WF'].obs
umap_df['layerCellType'] = mod_obs['layerCellType'].values
umap_df['CellType'] = mod_obs['CellType'].values
umap_df['cluster'] = mdata.obs['wnn_leiden'].values
# Copy additional metadata
for col in ['Depth', 'latency', 'onset', 'width', 'ratio_p2t']:
if col in mod_obs.columns:
umap_df[col] = mod_obs[col].values
# %%
# Create combined visualization (similar to pComb in R)
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# Plot 1: WNN UMAP colored by layerCellType
ax = axes[0]
for label in umap_df['layerCellType'].unique():
mask = umap_df['layerCellType'] == label
ax.scatter(umap_df.loc[mask, 'UMAP1'], umap_df.loc[mask, 'UMAP2'],
label=label, s=30, alpha=0.7)
ax.set_xlabel('UMAP1')
ax.set_ylabel('UMAP2')
ax.set_title('WNN UMAP - Layer+CellType')
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=8)
ax.set_aspect('equal')
# Plot 2: WNN UMAP colored by CellType
ax = axes[1]
for label in umap_df['CellType'].unique():
mask = umap_df['CellType'] == label
ax.scatter(umap_df.loc[mask, 'UMAP1'], umap_df.loc[mask, 'UMAP2'],
label=label, s=30, alpha=0.7)
ax.set_xlabel('UMAP1')
ax.set_ylabel('UMAP2')
ax.set_title('WNN UMAP - CellType')
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=8)
ax.set_aspect('equal')
# Plot 3: WNN UMAP colored by cluster
ax = axes[2]
scatter = ax.scatter(umap_df['UMAP1'], umap_df['UMAP2'],
c=umap_df['cluster'].astype(int), cmap='tab20', s=30, alpha=0.7)
ax.set_xlabel('UMAP1')
ax.set_ylabel('UMAP2')
ax.set_title(f'WNN UMAP - Clusters (res={RESOLUTION})')
ax.set_aspect('equal')
# Add cluster labels at centroids
for cluster in umap_df['cluster'].unique():
mask = umap_df['cluster'] == cluster
centroid_x = umap_df.loc[mask, 'UMAP1'].mean()
centroid_y = umap_df.loc[mask, 'UMAP2'].mean()
ax.annotate(str(cluster), (centroid_x, centroid_y), fontsize=10, fontweight='bold')
plt.tight_layout()
plt.show()
# %%
# Feature plots: Width and Peak-to-Trough ratio
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Width plot
ax = axes[0]
# Convert width to milliseconds if needed
width_ms = umap_df['width'] * 1000 if umap_df['width'].max() < 1 else umap_df['width']
scatter = ax.scatter(umap_df['UMAP1'], umap_df['UMAP2'],
c=width_ms, cmap='viridis', s=30, alpha=0.7)
plt.colorbar(scatter, ax=ax, label='Width (ms)')
ax.set_xlabel('UMAP1')
ax.set_ylabel('UMAP2')
ax.set_title('Width of Waveform')
ax.set_aspect('equal')
# Peak-to-trough ratio plot
ax = axes[1]
scatter = ax.scatter(umap_df['UMAP1'], umap_df['UMAP2'],
c=umap_df['ratio_p2t'], cmap='viridis', s=30, alpha=0.7)
plt.colorbar(scatter, ax=ax, label='Peak to Trough Ratio')
ax.set_xlabel('UMAP1')
ax.set_ylabel('UMAP2')
ax.set_title('Peak to Trough')
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
# %%
# Latency plot
fig, ax = plt.subplots(figsize=(8, 6))
# Use log scale for latency
log_latency = np.log1p(umap_df['latency'])
scatter = ax.scatter(umap_df['UMAP1'], umap_df['UMAP2'],
c=log_latency, cmap='viridis', s=30, alpha=0.7)
plt.colorbar(scatter, ax=ax, label='Log(Latency + 1)')
ax.set_xlabel('UMAP1')
ax.set_ylabel('UMAP2')
ax.set_title('Latency')
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
# %%
# Width vs Peak-to-Trough scatter plot with marginal histograms
# (equivalent to ggscatterhist in R)
# Prepare data with log transformations
plot_df = umap_df.copy()
plot_df['log_width'] = np.log(plot_df['width'] + 1e-6)
plot_df['log_ratio_p2t'] = np.log(plot_df['ratio_p2t'] + 1e-6)
# Create joint plot
g = sns.jointplot(
data=plot_df,
x='log_width',
y='log_ratio_p2t',
hue='layerCellType',
kind='scatter',
height=8,
alpha=0.7
)
g.ax_joint.set_xlabel('Log(Width)')
g.ax_joint.set_ylabel('Log(Peak-to-Trough Ratio)')
plt.suptitle('Width vs Peak-to-Trough Ratio by Cell Type', y=1.02)
plt.show()
# %% [markdown]
# ## 9. Modal Weight Visualization
#
# Since we used equal weights for the simple concatenation approach, this section shows placeholder visualizations. For true WNN-style adaptive weights, consider implementing Seurat's WNN algorithm or using muon's multi-modal analysis tools.
# %%
# Use the modality weights computed by Muon's WNN
# Muon stores them as 'WF:mod_weight', 'ISI:mod_weight', 'PSTH:mod_weight'
print("Modality weights from Muon WNN:")
print(f"Available weight columns: {[col for col in mdata_wnn.obs.columns if 'weight' in col.lower()]}")
# Rename weights for easier access
weight_mapping = {
'WF:mod_weight': 'WF_weight',
'ISI:mod_weight': 'ISI_weight',
'PSTH:mod_weight': 'PSTH_weight'
}
for old_name, new_name in weight_mapping.items():
if old_name in mdata_wnn.obs.columns:
mdata_wnn.obs[new_name] = mdata_wnn.obs[old_name]
mdata.obs[new_name] = mdata_wnn.obs[old_name].values
print(f"{new_name} - mean: {mdata_wnn.obs[old_name].mean():.3f}, std: {mdata_wnn.obs[old_name].std():.3f}")
# Add weights to umap_df for plotting
for weight_col in ['WF_weight', 'ISI_weight', 'PSTH_weight']:
if weight_col in mdata.obs.columns:
umap_df[weight_col] = mdata.obs[weight_col].values
# %%
# Visualize modality weights
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# Plot 1: Weight histograms (like p.hist in R)
ax = axes[0]
colors = {'WF_weight': '#F39922', 'ISI_weight': '#12A84B', 'PSTH_weight': '#0B78BE'}
labels = {'WF_weight': 'WF', 'ISI_weight': 'ISI', 'PSTH_weight': 'PSTH'}
for weight_col in ['WF_weight', 'ISI_weight', 'PSTH_weight']:
if weight_col in umap_df.columns:
ax.hist(umap_df[weight_col], bins=20, alpha=0.6,
color=colors[weight_col], label=labels[weight_col])
ax.set_xlabel('Weight')
ax.set_ylabel('Count')
ax.set_title('Modality Weight Distributions')
ax.legend()
# Plot 2: UMAP colored by dominant modality
ax = axes[1]
if all(col in umap_df.columns for col in ['WF_weight', 'ISI_weight', 'PSTH_weight']):
# Determine dominant modality for each cell
weight_cols = ['WF_weight', 'ISI_weight', 'PSTH_weight']
dominant = umap_df[weight_cols].idxmax(axis=1)
dominant_colors = dominant.map({'WF_weight': colors['WF_weight'],
'ISI_weight': colors['ISI_weight'],
'PSTH_weight': colors['PSTH_weight']})
for mod, color in [('WF_weight', colors['WF_weight']),
('ISI_weight', colors['ISI_weight']),
('PSTH_weight', colors['PSTH_weight'])]:
mask = dominant == mod
ax.scatter(umap_df.loc[mask, 'UMAP1'], umap_df.loc[mask, 'UMAP2'],
c=color, label=labels[mod], s=30, alpha=0.7)
ax.set_xlabel('UMAP1')
ax.set_ylabel('UMAP2')
ax.set_title('Dominant Modality per Cell')
ax.legend()
ax.set_aspect('equal')
# Plot 3: Stacked bar showing weight composition by cell type
ax = axes[2]
if 'CellType' in umap_df.columns:
weight_by_type = umap_df.groupby('CellType')[['WF_weight', 'ISI_weight', 'PSTH_weight']].mean()
weight_by_type.plot(kind='bar', stacked=True, ax=ax,
color=[colors['WF_weight'], colors['ISI_weight'], colors['PSTH_weight']])
ax.set_xlabel('Cell Type')
ax.set_ylabel('Mean Weight')
ax.set_title('Average Modality Weights by Cell Type')
ax.legend(['WF', 'ISI', 'PSTH'], bbox_to_anchor=(1.02, 1), loc='upper left')
ax.set_xticklabels(ax.get_xticklabels(), rotation=45, ha='right')
plt.tight_layout()
plt.show()
# %% [markdown]
# ## 10. Classification Analysis
#
# Compare classification accuracy between multimodal (WNN) and unimodal (WF, ISI, PSTH) representations for the five most prevalent cell types: E-4, E-5, FS-4, FS-5, SOM-nan.
#
# This mirrors the analysis in `classifyData.r` using scikit-learn's GradientBoostingClassifier.
# %%
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import balanced_accuracy_score
def classify_embeddings(embeddings, labels, n_repeats=20, test_size=0.2):
"""
Perform classification using GradientBoostingClassifier.
Mirrors the R code in classifyData.r which uses:
- method = "gbm" (Gradient Boosting Machine)
- 80/20 train/test split
- Multiple random seeds for robust evaluation
Args:
embeddings: numpy array of shape (n_samples, n_features)
labels: numpy array of class labels
n_repeats: number of times to repeat classification with different splits
test_size: fraction of data to use for testing
Returns:
per_class_accuracies: dict mapping class -> list of balanced accuracies
overall_accuracies: list of overall balanced accuracies (mean of per-class)
classes: list of unique class names
"""
classes = np.unique(labels)
per_class_acc = {cls: [] for cls in classes}
overall_acc = []
for seed in range(1, n_repeats + 1):
# Split data
X_train, X_test, y_train, y_test = train_test_split(
embeddings, labels, test_size=test_size,
random_state=seed, stratify=labels
)
# Train GBM classifier
clf = GradientBoostingClassifier(
n_estimators=100,
learning_rate=0.1,
max_depth=3,
random_state=seed
)
clf.fit(X_train, y_train)
# Predict
y_pred = clf.predict(X_test)
# Calculate per-class balanced accuracy
run_per_class_acc = []
for cls in classes:
# For each class, calculate balanced accuracy treating it as binary
cls_mask = (y_test == cls)
if cls_mask.sum() > 0:
# True positive rate (sensitivity) for this class
tp_rate = ((y_pred == cls) & (y_test == cls)).sum() / (y_test == cls).sum()
# True negative rate (specificity) for this class
tn_rate = ((y_pred != cls) & (y_test != cls)).sum() / (y_test != cls).sum()
balanced_acc = (tp_rate + tn_rate) / 2 * 100
per_class_acc[cls].append(balanced_acc)
run_per_class_acc.append(balanced_acc)
# Overall balanced accuracy = mean of per-class balanced accuracies
# This matches R's confusionMatrix approach
overall_bal_acc = np.mean(run_per_class_acc)
overall_acc.append(overall_bal_acc)
return per_class_acc, overall_acc, classes
# %%
# Filter for the 5 prevalent cell types used in the R analysis
# These are: E-4, E-5, FS-4, FS-5, SOM-nan
target_cell_types = ['E-4', 'E-5', 'FS-4', 'FS-5', 'SOM-nan']
# Get layerCellType labels from modality metadata
layer_cell_types = mdata.mod['WF'].obs['layerCellType'].values
# Create boolean mask for cells matching target types
mask = np.isin(layer_cell_types, target_cell_types)
# Get filtered labels
filtered_labels = layer_cell_types[mask]
print(f"Total cells: {len(layer_cell_types)}")
print(f"Cells in target types: {mask.sum()}")
print(f"\nCell type distribution:")
for ct in target_cell_types:
count = (filtered_labels == ct).sum()
print(f" {ct}: {count}")
# %%
# Compute higher-dimensional UMAP embeddings for classification
# This matches R's UMAP.components in constants.R
n_umap_dims = UMAP_COMPONENTS # Currently 10
# Dynamic key name based on dimensions
umap_key = f'X_umap_{n_umap_dims}d'
wnn_umap_key = f'X_wnn_umap_{n_umap_dims}d'
print(f"Computing {n_umap_dims}D UMAP embeddings for classification...")
print(f"Storing in keys: '{umap_key}' and '{wnn_umap_key}'")
# WF
sc.tl.umap(adata_wf, n_components=n_umap_dims, random_state=UMAP_SEED, min_dist=UMAP_MIN_DIST)
adata_wf.obsm[umap_key] = adata_wf.obsm['X_umap'].copy()
print(f" WF: {adata_wf.obsm[umap_key].shape}")
# ISI
sc.tl.umap(adata_isi, n_components=n_umap_dims, random_state=UMAP_SEED, min_dist=UMAP_MIN_DIST)
adata_isi.obsm[umap_key] = adata_isi.obsm['X_umap'].copy()
print(f" ISI: {adata_isi.obsm[umap_key].shape}")
# PSTH
sc.tl.umap(adata_psth, n_components=n_umap_dims, random_state=UMAP_SEED, min_dist=UMAP_MIN_DIST)
adata_psth.obsm[umap_key] = adata_psth.obsm['X_umap'].copy()
print(f" PSTH: {adata_psth.obsm[umap_key].shape}")
# WNN - on the combined neighbor graph
mu.tl.umap(mdata_wnn, neighbors_key='wnn', n_components=n_umap_dims,
random_state=UMAP_SEED, min_dist=UMAP_MIN_DIST)
mdata.obsm[wnn_umap_key] = mdata_wnn.obsm['X_umap'].copy()
print(f" WNN: {mdata.obsm[wnn_umap_key].shape}")
print("\nUMAP embeddings computed!")
# %%
# Run classification on each modality's UMAP embeddings
n_repeats = 20
# Prepare UMAP embeddings for each modality (using dynamic keys from previous cell)
embeddings_dict = {
'WNN': mdata.obsm[wnn_umap_key][mask],
'WF': adata_wf.obsm[umap_key][mask],
'ISI': adata_isi.obsm[umap_key][mask],
'PSTH': adata_psth.obsm[umap_key][mask],
}
# Add features (width, ratio_p2t) if available
if 'width' in mod_obs.columns and 'ratio_p2t' in mod_obs.columns:
width_vals = mod_obs['width'].values[mask]
ratio_vals = mod_obs['ratio_p2t'].values[mask]
# Replace NaN with 0
width_vals = np.nan_to_num(width_vals, nan=0.0)
ratio_vals = np.nan_to_num(ratio_vals, nan=0.0)
embeddings_dict['Features'] = np.column_stack([width_vals, ratio_vals])
print(f"Using {n_umap_dims}D UMAP embeddings for classification")
print(f"Embedding dimensions:")
for name, emb in embeddings_dict.items():
print(f" {name}: {emb.shape}")
# Run classification for each modality
results = {}
print("\nRunning classification analysis (this may take a few minutes)...\n")
for mod_name, embeddings in embeddings_dict.items():
print(f"Classifying {mod_name}...", end=" ")
per_class_acc, overall_acc, classes = classify_embeddings(
embeddings, filtered_labels, n_repeats=n_repeats
)
results[mod_name] = {
'per_class': per_class_acc,
'overall': overall_acc,
'classes': classes
}
print(f"Mean accuracy: {np.mean(overall_acc):.1f}% ± {np.std(overall_acc):.1f}%")
print("\nClassification complete!")
# %%
# Create summary DataFrame for plotting
summary_data = []
for mod_name, res in results.items():
for cls in res['classes']:
acc_values = res['per_class'][cls]
summary_data.append({
'Modality': mod_name,
'CellType': cls,
'Accuracy': np.mean(acc_values),
'Std': np.std(acc_values),
'SE': np.std(acc_values) / np.sqrt(len(acc_values))
})
summary_df = pd.DataFrame(summary_data)
# Also create raw data for statistical tests
raw_data = []
for mod_name, res in results.items():
for cls in res['classes']:
for run_idx, acc in enumerate(res['per_class'][cls]):
raw_data.append({
'Run': run_idx + 1,
'Modality': mod_name,
'CellType': cls,
'Accuracy': acc
})
raw_df = pd.DataFrame(raw_data)
print("Summary statistics:")
print(summary_df.pivot_table(index='CellType', columns='Modality', values='Accuracy').round(1))
# %%
# Plot classification accuracy by cell type and modality
# Similar to the R plot in classifyData.r
fig, ax = plt.subplots(figsize=(12, 6))
# Define colors and markers for each modality
modality_colors = {
'WNN': '#E41A1C', # Red
'WF': '#377EB8', # Blue
'ISI': '#4DAF4A', # Green
'PSTH': '#984EA3', # Purple
'Features': '#FF7F00' # Orange
}
modality_order = ['WNN', 'WF', 'ISI', 'PSTH', 'Features']
cell_types = sorted(summary_df['CellType'].unique())
# Plot with dodge (offset) for each modality
width = 0.15
x = np.arange(len(cell_types))
for i, mod in enumerate(modality_order):
if mod not in results:
continue
mod_data = summary_df[summary_df['Modality'] == mod].set_index('CellType')
# Reorder to match cell_types
means = [mod_data.loc[ct, 'Accuracy'] if ct in mod_data.index else 0 for ct in cell_types]
sems = [mod_data.loc[ct, 'SE'] if ct in mod_data.index else 0 for ct in cell_types]
offset = (i - len(modality_order)/2 + 0.5) * width
ax.errorbar(x + offset, means, yerr=sems,
fmt='o-', color=modality_colors[mod],
label=mod, capsize=3, markersize=8, linewidth=2)
ax.set_xticks(x)
ax.set_xticklabels(cell_types, rotation=45, ha='right')
ax.set_xlabel('Cell Type')
ax.set_ylabel('Classification Accuracy (%)')
ax.set_title(f'Classification Accuracy at Embedding-D = {n_umap_dims}')
ax.legend(title='Modality', bbox_to_anchor=(1.02, 1), loc='upper left')
ax.set_ylim(0, 105)
ax.axhline(y=100, color='gray', linestyle='--', alpha=0.3)
ax.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.show()
# %%
# Overall accuracy comparison (barplot with error bars)
fig, ax = plt.subplots(figsize=(8, 5))
modalities = []
means = []
sems = []
for mod in modality_order:
if mod in results:
modalities.append(mod)
means.append(np.mean(results[mod]['overall']))
sems.append(np.std(results[mod]['overall']) / np.sqrt(len(results[mod]['overall'])))
x = np.arange(len(modalities))
colors = [modality_colors[m] for m in modalities]
bars = ax.bar(x, means, yerr=sems, capsize=5, color=colors, alpha=0.8, edgecolor='black')
ax.set_xticks(x)
ax.set_xticklabels(modalities)
ax.set_ylabel('Overall Balanced Accuracy (%)')
ax.set_title('Overall Classification Accuracy by Modality')
ax.set_ylim(0, 100)
ax.axhline(y=20, color='gray', linestyle='--', alpha=0.5, label='Chance level (5 classes)')
# Add value labels on bars
for bar, mean, sem in zip(bars, means, sems):
ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + sem + 2,
f'{mean:.1f}%', ha='center', va='bottom', fontsize=10)
plt.tight_layout()
plt.show()
# Print summary
print("\nOverall Accuracy Summary:")
print("-" * 40)
for mod in modality_order:
if mod in results:
acc = results[mod]['overall']
print(f"{mod:10s}: {np.mean(acc):5.1f}% ± {np.std(acc):4.1f}%")
# %% [markdown]
# # Print summary
# print("=" * 50)
# print("Analysis Summary")
# print("=" * 50)
# print(f"\nTotal cells: {mdata.n_obs}")
# print(f"\nModalities processed:")
# print(f" - ISI: {mdata.mod['ISI'].n_vars} features")
# print(f" - WF: {mdata.mod['WF'].n_vars} features")
# print(f" - PSTH: {mdata.mod['PSTH'].n_vars} features")
# print(f"\nWNN UMAP dimensions: {mdata.obsm['X_wnn_umap'].shape}")
# print(f"\nClusters found: {mdata.obs['wnn_leiden'].nunique()}")
#
# # Get metadata from modality
# mod_obs = mdata.mod['WF'].obs
# print(f"\nCell types:")
# print(mod_obs['CellType'].value_counts())
# print(f"\nLayer distribution:")
# print(mod_obs['Layer'].value_counts())
#
# # Print weight summary if available
# if 'WF_weight' in mdata.obs.columns:
# print(f"\nModality weight summary:")
# for mod in ['WF', 'ISI', 'PSTH']:
# col = f'{mod}_weight'
# if col in mdata.obs.columns:
# print(f" - {mod}: mean={mdata.obs[col].mean():.3f}, std={mdata.obs[col].std():.3f}")
# %%
# Print summary
print("=" * 50)
print("Analysis Summary")
print("=" * 50)
print(f"\nTotal cells: {mdata.n_obs}")
print(f"\nModalities processed:")
print(f" - ISI: {mdata.mod['ISI'].n_vars} features")
print(f" - WF: {mdata.mod['WF'].n_vars} features")
print(f" - PSTH: {mdata.mod['PSTH'].n_vars} features")
# print(f"\nJoint embedding dimensions: {mdata.obsm['X_joint_pca'].shape[1]}")
print(f"\nClusters found: {mdata.obs['wnn_leiden'].nunique()}")
# print(f"\nCell types:")
# print(mdata.obs['CellType'].value_counts())
# print(f"\nLayer distribution:")
# print(mdata.obs['Layer'].value_counts())
# %%
# Save processed data
# mdata.write('processed_jianing.h5mu')
# print("Saved processed data to processed_jianing.h5mu")
# %%
# %%

processJianing.ipynb at commit 08819f4, no license · at the source

Overview

Authors: Eric Kenji Lee¹, Asım E Gül², Greggory Heller³, Anna Lakunina^4,5,6, Han Yu^7,8, Andrew Shelton⁶, Shawn Olsen⁶, Nicholas A Steinmetz⁹, Cole Hurwitz⁸, Santiago Jaramillo^4,5, Pawel F Przytycki¹⁰, Chandramouli Chandrasekaran^1,11,12,13

ORCID iDs: Eric Kenji Lee, Greggory Heller, Anna Lakunina, Han Yu, Andrew Shelton, Shawn Olsen, Nicholas A Steinmetz, Cole Hurwitz, Santiago Jaramillo, Pawel F Przytycki, Chandramouli Chandrasekaran

13 affiliations

Department of Psychological and Brain Sciences, Boston University, Boston, MA USA
Department of Psychology, Boğaziçi University, Beşiktaş, Istanbul Turkey
Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA USA
Department of Biology, University of Oregon, Eugene, OR USA
Institute of Neuroscience, University of Oregon, Eugene, OR USA
Allen Institute for Neural Dynamics, Seattle, WA USA
Department of Electrical Engineering, Columbia University, New York City, NY USA
Zuckerman Institute, Columbia University, New York City, NY USA
Department of Neurobiology and Biophysics, University of Washington, Seattle, WA USA
Faculty of Computing & Data Sciences, Boston University, Boston, MA USA
Center for Systems Neuroscience, Boston University, Boston, MA USA
Department of Biomedical Engineering, Boston University, Boston, MA USA
Department of Anatomy & Neurobiology, Boston University, Boston, MA USA

Institutions: Boston University (United States); Boğaziçi University (Türkiye); Massachusetts Institute of Technology (United States); University of Oregon (United States); Allen Institute for Neural Dynamics (United States); Mortimer B. Zuckerman Mind Brain Behavior Institute (United States); Columbia University (United States); University of Washington (United States)

Journal: Nature communications, volume 17, issue 1, article 5241

Dates: received 29 September 2025; accepted 19 March 2026; published online 15 April 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s41467-026-71331-0 · PMID 41986318 · PMCID PMC13260831 · OpenAlex W7154465353

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), extracellular electrophysiology (units, LFP) (modality), mouse (organism)

Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Machine learning, Evoked potentials, Graphs, fMRI & imaging, Single-unit activity, calcium imaging

Keywords: Machine learning, Computational neuroscience, Extracellular recording, Functional clustering

MeSH: Electrophysiological Phenomena*, Neurons*, Action Potentials, Animals, Mice, Multiomics (* major topic)

Topic: Cardiac electrophysiology and arrhythmias (Cardiology and Cardiovascular Medicine, Medicine), according to OpenAlex

Funding: NINDS NIH HHS (R01 NS131993, RF1 NS131993); U.S. Department of Health & Human Services | NIH | National Institute of Neurological Disorders and Stroke (U01NS113252, R01NS122969, RF1NS131993, R00NS092972, R01NS121409, F31NS131018, R21NS135361); U.S. Department of Health & Human Services | NIH | National Institute of Neurological Disorders and Stroke (NINDS) (R01NS122969, RF1NS131993, R21NS135361, R01NS121409, R00NS092972, U01NS113252, F31NS131018); Whitehall Foundation (2019-12-77); Brain and Behavior Research Foundation (27923); U.S. Department of Health & Human Services | NIH | National Institute on Deafness and Other Communication Disorders (NIDCD) (R01DC01553); U.S. Department of Health & Human Services | NIH | National Institute on Deafness and Other Communication Disorders (R01DC01553)

Citations: cited by 2 papers (Europe PMC); 105 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 8 matches between paragraphs and lines of code.

EricKenjiLee/PhysMAP_Manuscript

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 08819f45c7742001a760096c64bb00109bd8062e, 11 February 2026

Languages: R (17), Python (3), Jupyter (2), MATLAB (1)

Size: 77 files, 23 scripts

Software Heritage: not archived

Found in: “Code availability”

Holds: README, environment (poetry.lock, pyproject.toml, juxtacellular/requirements.txt), 2 notebooks

Not found: license file, CITATION.cff, tests, continuous integration, documentation

Tools: caret (9 files), reshape2 (9 files), tidyverse (8 files), NumPy (5 files), pandas (5 files), Matplotlib (4 files), seaborn (3 files), anndata (2 files), ggplot2 (2 files), ggpubr (2 files), reticulate (2 files), Scanpy (2 files), scikit-learn (2 files), SciPy (2 files), Seurat (2 files), cowplot (1 file), PyTorch (1 file), UMAP (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

24 files

CellExplorerv2/AllenOnly.r, R, 245 lines
CellExplorerv2/classifyCellExp.r, R, 118 lines
CellExplorerv2/processCellExp.m, MATLAB, 110 lines
CellExplorerv2/processCellExp.r, R, 248 lines
ConfMatPlotting.R, R, 39 lines
InvivoA1/processSantiago.r, R, 331 lines, 2 matches
PairOfVAEs_cuda.ipynb, Jupyter, 1,278 lines, 2 matches
constants.R, R, 17 lines
juxtacellular/Generate_Traces.R, R, 47 lines
juxtacellular/classifyData.r, R, 264 lines
juxtacellular/classifyDataTwoModality.r, R, 235 lines
juxtacellular/classifyRaw.R, R, 125 lines
juxtacellular/generate_example_data.py, Python, 130 lines
juxtacellular/helperFunctions.r, R, 156 lines
juxtacellular/physmap_app.py, Python, 471 lines
juxtacellular/physmap_utils.py, Python, 490 lines
juxtacellular/processJianing.ipynb, Jupyter, 1,152 lines, 3 matches
juxtacellular/processJianing.r, R, 199 lines, 1 match
juxtacellular/processJianingTwoModality.R, R, 189 lines
lookupTable/lookupTable.r, R, 38 lines
lookupTable/mapReference.r, R, 231 lines
lookupTable/mapReferenceAllen.r, R, 229 lines
lookupTable/mapReferenceAllen_Kenji.r, R, 242 lines
README.md, Text, 70 lines

At the source: github.com/EricKenjiLee/PhysMAP_Manuscript

Zenodo 18239113

License: CC-BY-4.0

State: the link answers, verified on 29 September 2026

Evidence: files inventoried

Size: 1 file

Software Heritage: not checked

Found in: “Code availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)

29 September 2026: the link answers (HTTP 200)

At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: EricKenjiLee/PhysMAP_Manuscript, Zenodo 18239113

Read it in the paper: doi.org/10.1038/s41467-026-71331-0.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
23 scripts, each with its path and the digest of its content;
8 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: EricKenjiLee/PhysMAP_Manuscript, Zenodo 18239113

Read it in the paper: doi.org/10.1038/s41467-026-71331-0.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 4 keywords, 6 MeSH terms, 7 funders, 92 references.

Cite

This paper

Lee, E. K., Gül, A. E., Heller, G., Lakunina, A., Yu, H., Shelton, A., Olsen, S., Steinmetz, N. A., Hurwitz, C., Jaramillo, S., Przytycki, P. F., & Chandrasekaran, C. (2026). A multimodal approach for visualizing and identifying electrophysiological cell types in vivo. Nature communications, 17(1), 5241. https://doi.org/10.1038/s41467-026-71331-0

BibTeX

@article{lee2026multimodal, author = {Lee, Eric Kenji and Gül, Asım E and Heller, Greggory and Lakunina, Anna and Yu, Han and Shelton, Andrew and Olsen, Shawn and Steinmetz, Nicholas A and Hurwitz, Cole and Jaramillo, Santiago and Przytycki, Pawel F and Chandrasekaran, Chandramouli}, title = {{A multimodal approach for visualizing and identifying electrophysiological cell types in vivo}}, journal = {Nature communications}, year = {2026}, month = apr, volume = {17}, number = {1}, pages = {5241}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-71331-0}, url = {https://doi.org/10.1038/s41467-026-71331-0}, pmid = {41986318}, pmcid = {PMC13260831} }

RIS

TY - JOUR AU - Lee, Eric Kenji AU - Gül, Asım E AU - Heller, Greggory AU - Lakunina, Anna AU - Yu, Han AU - Shelton, Andrew AU - Olsen, Shawn AU - Steinmetz, Nicholas A AU - Hurwitz, Cole AU - Jaramillo, Santiago AU - Przytycki, Pawel F AU - Chandrasekaran, Chandramouli TI - A multimodal approach for visualizing and identifying electrophysiological cell types in vivo T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/04/15 VL - 17 IS - 1 SP - 5241 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-71331-0 UR - https://doi.org/10.1038/s41467-026-71331-0 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-71331-0", "type": "article-journal", "title": "A multimodal approach for visualizing and identifying electrophysiological cell types in vivo", "container-title": "Nature communications", "author": [ { "family": "Lee", "given": "Eric Kenji" }, { "family": "Gül", "given": "Asım E" }, { "family": "Heller", "given": "Greggory" }, { "family": "Lakunina", "given": "Anna" }, { "family": "Yu", "given": "Han" }, { "family": "Shelton", "given": "Andrew" }, { "family": "Olsen", "given": "Shawn" }, { "family": "Steinmetz", "given": "Nicholas A" }, { "family": "Hurwitz", "given": "Cole" }, { "family": "Jaramillo", "given": "Santiago" }, { "family": "Przytycki", "given": "Pawel F" }, { "family": "Chandrasekaran", "given": "Chandramouli" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "5241", "DOI": "10.1038/s41467-026-71331-0", "PMID": "41986318", "PMCID": "PMC13260831", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-71331-0", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 15 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-76939-w [code]: HIPPIE: a generative model for electrophysiological analysis across species, technologies, and modalities.
Journal: Nature communications
In common: caret, reticulate, UMAP, 15 other tools, mouse, 10 references
[2] doi:10.1038/s41592-026-03076-z [code]: Neuropixels Opto: combining high-resolution electrophysiology and optogenetics.
Journal: Nature methods
In common: PyTorch, pandas, SciPy, 2 other tools, extracellular electrophysiology (units, LFP), mouse, 5 references, 2 authors
[3] doi:10.1038/s41586-026-10629-x [code]: Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: reticulate, UMAP, anndata, 13 other tools, genetics / omics, mouse
[4] doi:10.1016/j.xcrm.2026.102651 [code]: Integrative CSF profiling identifies disease-specific immune responses in leptomeningeal disease.
Journal: Cell reports. Medicine
In common: reticulate, UMAP, anndata, 13 other tools, genetics / omics
[5] doi:10.1016/j.xcrm.2026.102766 [code]: A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: UMAP, anndata, Scanpy, 12 other tools, genetics / omics, 2 references
[6] doi:10.1186/s13059-026-04177-w [code]: Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: reticulate, UMAP, anndata, 12 other tools, genetics / omics, mouse, 1 reference
[7] doi:10.1038/s41586-026-10490-y [code]: Lineage and organ signals sequentially build organ intrinsic nervous systems.
Journal: Nature
In common: UMAP, anndata, Scanpy, 12 other tools, mouse, 2 references
[8] doi:10.1186/s13073-026-01704-z [code]: Gene expression profiling enables refined parcellation of cortical layers in the heterogeneous human cerebral cortex.
Journal: Genome medicine
In common: UMAP, anndata, Scanpy, 13 other tools, genetics / omics, mouse
[9] doi:10.1038/s44320-026-00208-7 [code]: Interpretable deep generative ensemble learning for single-cell omics with Hydra.
Journal: Molecular systems biology
In common: reticulate, UMAP, anndata, 12 other tools, 1 reference
[10] doi:10.21203/rs.3.rs-9676637/v1 [code]: A Comprehensive Benchmarking of Spatial Deconvolution and Domain Detection Methods across Diverse Tissues and Spatial Transcriptomic Technologies
Journal: Research Square (preprint)
In common: reticulate, UMAP, anndata, 12 other tools, genetics / omics, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.