# Import
import pandas as pd
import numpy as np
from typing import Optional, Tuple, List, Dict, Any
from tqdm import tqdm
import time
import sqlite3
import os
from ..core.find_cell_neighbors_3d import(
create_neighbor_edge_table_database_3d,
get_anndata_from_database,
query_cell_type_pairs,
get_graph_statistics,
save_surfaces_to_pickle,
load_surfaces_from_pickle,
save_graph_state_to_pickle,
load_graph_state_from_pickle
)
from ..core.find_cell_neighbors_centroid_3d import create_neighbor_edge_table_database_centroid_3d
from ..core.compute_interscellar_volumes_3d_absolute import(
build_interscellar_volume_database_from_neighbors,
create_global_interscellar_mesh_zarr,
create_global_cell_only_volumes_zarr,
export_interscellar_volumes_to_duckdb,
get_anndata_from_interscellar_database,
export_interscellar_volumes_to_anndata,
ANNDATA_AVAILABLE
)
# API: Wrapper functions
[docs]
def find_cell_neighbors_3d(
ome_zarr_path: str,
metadata_csv_path: str,
max_distance_um: float = 0.5,
voxel_size_um: tuple = (0.56, 0.28, 0.28),
centroid_prefilter_radius_um: float = 75.0,
cell_id: str = 'CellID',
cell_type: str = 'phenotype',
centroid_x: str = 'X_centroid',
centroid_y: str = 'Y_centroid',
centroid_z: str = 'Z_centroid',
db_path: Optional[str] = None,
output_csv: Optional[str] = None,
output_anndata: Optional[str] = None,
n_jobs: int = 1,
return_connection: bool = False,
save_surfaces_pickle: Optional[str] = None,
load_surfaces_pickle: Optional[str] = None,
save_graph_state_pickle: Optional[str] = None
) -> Tuple[Optional[pd.DataFrame], Optional[object], Optional[object]]:
print("=" * 60)
print("InterSCellar: Surface-based Cell Neighbor Detection - 3D")
print("=" * 60)
overall_start_time = time.time()
print(f"\n1. Loading metadata from: {metadata_csv_path}...")
step1_start = time.time()
try:
metadata_df = pd.read_csv(metadata_csv_path)
print(f"Loaded {len(metadata_df)} cells")
except Exception as e:
raise ValueError(f"Error loading metadata CSV: {e}")
required_cols = [cell_id, cell_type, centroid_x, centroid_y, centroid_z]
missing_cols = [col for col in required_cols if col not in metadata_df.columns]
if missing_cols:
raise ValueError(f"Missing required columns in metadata: {missing_cols}")
step1_time = time.time() - step1_start
print(f"Step 1 completed in {step1_time:.2f} seconds")
metadata_dir = os.path.dirname(metadata_csv_path) if os.path.dirname(metadata_csv_path) else "."
base_name = os.path.splitext(os.path.basename(metadata_csv_path))[0]
if db_path is None:
db_path = os.path.join(metadata_dir, f"{base_name}_neighbor_graph.db")
print(f"db_path: {db_path}")
if output_csv is None:
output_csv = os.path.join(metadata_dir, f"{base_name}_neighbors_3d.csv")
print(f"output_csv: {output_csv}")
if output_anndata is None:
output_anndata = os.path.join(metadata_dir, f"{base_name}_neighbors_3d.h5ad")
print(f"output_anndata: {output_anndata}")
if save_surfaces_pickle is None:
base_name = os.path.splitext(db_path)[0]
save_surfaces_pickle = f"{base_name}_surfaces.pkl"
print(f"Surfaces pickle path: {save_surfaces_pickle}")
if save_graph_state_pickle is None:
base_name = os.path.splitext(db_path)[0]
save_graph_state_pickle = f"{base_name}_graph_state.pkl"
print(f"Graph state pickle path: {save_graph_state_pickle}")
print(f"\n2. Building neighbor graph...")
print(f"Parameters: max_distance={max_distance_um}μm, n_jobs={n_jobs}")
if max_distance_um == 0.0:
print(f"Mode: Touching cells only")
else:
print(f"Mode: Touching cells + near-neighbors")
step2_start = time.time()
try:
conn = create_neighbor_edge_table_database_3d(
ome_zarr_path=ome_zarr_path,
metadata_df=metadata_df,
max_distance_um=max_distance_um,
voxel_size_um=voxel_size_um,
centroid_prefilter_radius_um=centroid_prefilter_radius_um,
db_path=db_path,
cell_id=cell_id,
cell_type=cell_type,
centroid_x=centroid_x,
centroid_y=centroid_y,
centroid_z=centroid_z,
output_csv=output_csv,
output_anndata=output_anndata,
n_jobs=n_jobs,
save_surfaces_pickle=save_surfaces_pickle,
load_surfaces_pickle=load_surfaces_pickle,
save_graph_state_pickle=save_graph_state_pickle
)
print(f"Saving global surface to: {save_surfaces_pickle}")
print(f"Saving graph state to: {save_graph_state_pickle}")
step2_time = time.time() - step2_start
print(f"Neighbor graph created successfully")
print(f"Step 2 completed in {step2_time:.2f} seconds")
except Exception as e:
raise RuntimeError(f"Error in neighbor detection pipeline: {e}")
print(f"\n3. Retrieving results...")
step3_start = time.time()
neighbor_table_df = None
if output_csv:
try:
neighbor_table_df = pd.read_sql_query("SELECT * FROM neighbors ORDER BY pair_id", conn)
print(f"Neighbor table: {len(neighbor_table_df)} pairs")
except Exception as e:
print(f"Warning: Could not retrieve neighbor table: {e}")
adata = None
try:
adata = get_anndata_from_database(conn)
if adata is not None:
print(f"AnnData object created: {adata.shape}")
else:
print(f"Warning: AnnData not available (install with: pip install anndata)")
except Exception as e:
print(f"Warning: Could not create AnnData object: {e}")
try:
stats = get_graph_statistics(conn)
print(f"Graph statistics: {stats['total_cells']} cells, {stats['total_edges']} pairs")
except Exception as e:
print(f"Warning: Could not retrieve statistics: {e}")
step3_time = time.time() - step3_start
print(f"Step 3 completed in {step3_time:.2f} seconds")
overall_time = time.time() - overall_start_time
print(f"\n4. Pipeline completed successfully!")
print(f"Total execution time: {overall_time:.2f} seconds")
print(f"Database: {db_path}")
if output_csv:
print(f"CSV output: {output_csv}")
if output_anndata:
print(f"AnnData output: {output_anndata}")
print("=" * 60)
if return_connection:
return neighbor_table_df, adata, conn
else:
conn.close()
return neighbor_table_df, adata, None
[docs]
def find_cell_neighbors_centroid_3d(
metadata_csv_path: str,
radius_um: float,
voxel_size_um: tuple = (0.56, 0.28, 0.28),
db_path: Optional[str] = None,
output_csv: Optional[str] = None,
output_cells_csv: Optional[str] = None,
output_anndata: Optional[str] = None,
cell_id: str = "CellID",
cell_type: str = "phenotype",
centroid_x: str = "X_centroid",
centroid_y: str = "Y_centroid",
centroid_z: str = "Z_centroid",
return_connection: bool = False,
) -> Tuple[Optional[pd.DataFrame], Optional[object], Optional[object]]:
print("=" * 60)
print("InterSCellar: Centroid-Sphere Cell Neighbor Detection - 3D")
print("=" * 60)
overall_start_time = time.time()
print(f"\n1. Loading metadata from: {metadata_csv_path}...")
step1_start = time.time()
try:
metadata_df = pd.read_csv(metadata_csv_path)
print(f"Loaded {len(metadata_df)} cells")
except Exception as e:
raise ValueError(f"Error loading metadata CSV: {e}")
required_cols = [cell_id, cell_type, centroid_x, centroid_y, centroid_z]
missing_cols = [col for col in required_cols if col not in metadata_df.columns]
if missing_cols:
raise ValueError(f"Missing required columns in metadata: {missing_cols}")
step1_time = time.time() - step1_start
print(f"Step 1 completed in {step1_time:.2f} seconds")
metadata_dir = os.path.dirname(metadata_csv_path) if os.path.dirname(metadata_csv_path) else "."
base_name = os.path.splitext(os.path.basename(metadata_csv_path))[0]
if db_path is None:
db_path = os.path.join(metadata_dir, f"{base_name}_centroid_neighbor_graph_3d.db")
print(f"db_path: {db_path}")
if output_csv is None:
output_csv = os.path.join(metadata_dir, f"{base_name}_centroid_neighbors_3d.csv")
print(f"output_csv: {output_csv}")
if output_anndata is None:
output_anndata = os.path.join(metadata_dir, f"{base_name}_centroid_neighbors_3d.h5ad")
print(f"output_anndata: {output_anndata}")
if output_cells_csv:
print(f"output_cells_csv: {output_cells_csv}")
print(f"\n2. Building centroid-based neighbor graph...")
print(f"Parameters: radius={radius_um}μm")
step2_start = time.time()
conn = create_neighbor_edge_table_database_centroid_3d(
metadata_df=metadata_df,
radius_um=radius_um,
voxel_size_um=voxel_size_um,
db_path=db_path,
output_csv=output_csv,
output_anndata=output_anndata,
output_cells_csv=output_cells_csv,
cell_id=cell_id,
cell_type=cell_type,
centroid_x=centroid_x,
centroid_y=centroid_y,
centroid_z=centroid_z,
)
step2_time = time.time() - step2_start
print(f"Step 2 completed in {step2_time:.2f} seconds")
print("Centroid neighbor graph created successfully")
print(f"\n3. Retrieving results...")
neighbor_table_df = pd.read_sql_query("SELECT * FROM neighbors ORDER BY pair_id", conn)
print(f"Neighbor table: {len(neighbor_table_df)} pairs")
adata = get_anndata_from_database(conn)
if adata is not None:
print(f"AnnData object created: {adata.shape}")
else:
print("Warning: AnnData not available (install with: pip install anndata)")
try:
stats = get_graph_statistics(conn)
print(f"Graph statistics: {stats['total_cells']} cells, {stats['total_edges']} pairs")
except Exception as e:
print(f"Warning: Could not retrieve statistics: {e}")
overall_time = time.time() - overall_start_time
print(f"\n4. Pipeline completed successfully!")
print(f"Total execution time: {overall_time:.2f} seconds")
print(f"Database: {db_path}")
if output_csv:
print(f"CSV output: {output_csv}")
if output_anndata:
print(f"AnnData output: {output_anndata}")
print("=" * 60)
if return_connection:
return neighbor_table_df, adata, conn
conn.close()
return neighbor_table_df, adata, None
[docs]
def compute_interscellar_volumes_3d(
ome_zarr_path: str,
neighbor_pairs_csv: str,
global_surface_pickle: Optional[str] = None,
halo_bboxes_pickle: Optional[str] = None,
neighbor_db_path: Optional[str] = None,
voxel_size_um: tuple = (0.56, 0.28, 0.28),
db_path: Optional[str] = None,
output_csv: Optional[str] = None,
output_anndata: Optional[str] = None,
output_mesh_zarr: Optional[str] = None,
output_cell_only_zarr: Optional[str] = None,
max_distance_um: float = 3.0,
intracellular_threshold_um: float = 1.0,
n_jobs: int = 4,
return_connection: bool = False,
intermediate_results_dir: str = "intermediate_interscellar_results",
resume: Optional[bool] = None,
output_name_tag: str = "absolute",
) -> Tuple[Optional[pd.DataFrame], Optional[object], Optional[object]]:
print("=" * 60)
print("InterSCellar: Volume Computation - 3D")
print("=" * 60)
overall_start_time = time.time()
csv_dir = os.path.dirname(neighbor_pairs_csv) if os.path.dirname(neighbor_pairs_csv) else "."
csv_base_name = os.path.splitext(os.path.basename(neighbor_pairs_csv))[0]
csv_base_name = csv_base_name.replace("_neighbors_3d", "").replace("_neighbors", "").replace("neighbors", "")
if neighbor_db_path is None:
possible_db = os.path.join(csv_dir, f"{csv_base_name}_neighbor_graph.db")
if os.path.exists(possible_db):
neighbor_db_path = possible_db
print(f"Detected neighbor_db_path: {neighbor_db_path}")
else:
alt_db = os.path.join(csv_dir, os.path.basename(neighbor_pairs_csv).replace("_neighbors_3d.csv", "_neighbor_graph.db").replace("_neighbors.csv", "_neighbor_graph.db").replace(".csv", "_neighbor_graph.db"))
if os.path.exists(alt_db):
neighbor_db_path = alt_db
print(f"Detected neighbor_db_path: {neighbor_db_path}")
# Outputs are tagged to distinguish this pipeline from the adaptive one, which
# writes <stem>_adaptive_*. Inputs (neighbor graph db, surfaces/halo pickles) keep
# the untagged stem, so their auto-detection below is unaffected.
output_base_name = f"{csv_base_name}_{output_name_tag}" if output_name_tag else csv_base_name
if db_path is None:
db_path = os.path.join(csv_dir, f"{output_base_name}_interscellar_volumes.db")
print(f"db_path: {db_path}")
if output_csv is None:
output_csv = os.path.join(csv_dir, f"{output_base_name}_volumes.csv")
print(f"output_csv: {output_csv}")
if output_anndata is None:
output_anndata = os.path.join(csv_dir, f"{output_base_name}_volumes.h5ad")
print(f"output_anndata: {output_anndata}")
if global_surface_pickle is None or halo_bboxes_pickle is None:
if neighbor_db_path:
base_name = os.path.splitext(neighbor_db_path)[0]
search_dir = os.path.dirname(neighbor_db_path) if os.path.dirname(neighbor_db_path) else "."
else:
csv_dir = os.path.dirname(neighbor_pairs_csv) if os.path.dirname(neighbor_pairs_csv) else "."
csv_basename = os.path.basename(neighbor_pairs_csv)
base_name_stem = csv_basename.replace("_neighbors_3d.csv", "").replace("_neighbors.csv", "").replace("neighbors.csv", "").replace(".csv", "")
possible_db = os.path.join(csv_dir, f"{base_name_stem}_neighbor_graph.db")
if os.path.exists(possible_db):
base_name = os.path.splitext(possible_db)[0]
search_dir = csv_dir
else:
base_name = os.path.join(csv_dir, base_name_stem)
search_dir = csv_dir
if global_surface_pickle is None:
possible_surfaces = [
f"{base_name}_surfaces.pkl",
f"{base_name}_graph_surfaces.pkl",
]
if os.path.exists(search_dir):
for f in os.listdir(search_dir):
if f.endswith('.pkl') and 'surface' in f.lower() and os.path.basename(base_name) in f:
possible_surfaces.append(os.path.join(search_dir, f))
for path in possible_surfaces:
if os.path.exists(path):
global_surface_pickle = path
print(f"Detected global_surface_pickle: {global_surface_pickle}")
break
else:
global_surface_pickle = f"{base_name}_surfaces.pkl"
if halo_bboxes_pickle is None:
possible_halo = [
f"{base_name}_halo_bboxes.pkl",
f"{base_name.replace('_surfaces', '_halo_bboxes')}.pkl",
]
if os.path.exists(search_dir):
for f in os.listdir(search_dir):
if f.endswith('.pkl') and 'halo' in f.lower() and os.path.basename(base_name) in f:
possible_halo.append(os.path.join(search_dir, f))
for path in possible_halo:
if os.path.exists(path):
halo_bboxes_pickle = path
print(f"Detected halo_bboxes_pickle: {halo_bboxes_pickle}")
break
else:
halo_bboxes_pickle = f"{base_name}_halo_bboxes.pkl"
if output_mesh_zarr is None:
base_name = os.path.splitext(db_path)[0]
if base_name.endswith('_interscellar_volumes'):
base_name = base_name[:-len('_interscellar_volumes')]
output_mesh_zarr = f"{base_name}_interscellar_volumes.zarr"
print(f"output_mesh_zarr: {output_mesh_zarr}")
if output_cell_only_zarr is None:
base_name = os.path.splitext(db_path)[0]
if base_name.endswith('_interscellar_volumes'):
base_name = base_name[:-len('_interscellar_volumes')]
output_cell_only_zarr = f"{base_name}_cell_only_volumes.zarr"
print(f"output_cell_only_zarr: {output_cell_only_zarr}")
print(f"\n1. Validating input files...")
step1_start = time.time()
required_files = [ome_zarr_path]
if not neighbor_pairs_csv and not neighbor_db_path:
raise ValueError("Must provide either neighbor_pairs_csv or neighbor_db_path")
if neighbor_pairs_csv and not os.path.exists(neighbor_pairs_csv):
if neighbor_db_path and os.path.exists(neighbor_db_path):
print(f"Warning: neighbor_pairs_csv not found: {neighbor_pairs_csv}")
print(f"Use db instead: {neighbor_db_path}")
else:
raise FileNotFoundError(f"neighbor_pairs_csv not found: {neighbor_pairs_csv}")
elif neighbor_pairs_csv:
required_files.append(neighbor_pairs_csv)
if neighbor_db_path and not os.path.exists(neighbor_db_path):
if neighbor_pairs_csv and os.path.exists(neighbor_pairs_csv):
print(f"Warning: neighbor_db_path not found: {neighbor_db_path}")
print(f"Use CSV instead: {neighbor_pairs_csv}")
else:
raise FileNotFoundError(f"neighbor_db_path not found: {neighbor_db_path}")
for file_path in required_files:
if not os.path.exists(file_path):
raise FileNotFoundError(f"Required file not found: {file_path}")
if global_surface_pickle and not os.path.exists(global_surface_pickle):
print(f"Warning: global_surface_pickle not found: {global_surface_pickle}")
print(f"Loading from graph_state.pkl as fallback")
elif global_surface_pickle:
print(f"Using global_surface_pickle: {global_surface_pickle}")
if halo_bboxes_pickle and not os.path.exists(halo_bboxes_pickle):
print(f"Warning: halo_bboxes_pickle not found: {halo_bboxes_pickle}")
elif halo_bboxes_pickle:
print(f"Using halo_bboxes_pickle: {halo_bboxes_pickle}")
print(f"All input files found")
step1_time = time.time() - step1_start
print(f"Step 1 completed in {step1_time:.2f} seconds")
print(f"\n2. Computing interscellar volumes...")
print(f"Parameters: max_distance={max_distance_um}μm, intracellular_threshold={intracellular_threshold_um}μm")
print(f"Voxel size: {voxel_size_um} μm")
step2_start = time.time()
try:
import zarr
print(f"Loading segmentation mask from: {ome_zarr_path}")
zarr_group = zarr.open(ome_zarr_path, mode='r')
if 'labels' in zarr_group:
mask_3d = zarr_group['labels'][0, 0]
elif '0' in zarr_group and '0' in zarr_group['0']:
mask_3d = zarr_group['0']['0'][0, 0]
else:
mask_3d = None
for key in zarr_group.keys():
if hasattr(zarr_group[key], 'shape') and len(zarr_group[key].shape) >= 3:
mask_3d = zarr_group[key]
break
if mask_3d is None:
raise ValueError("Could not find 3D segmentation mask in ome-zarr file")
print(f"Mask shape: {mask_3d.shape}, dtype: {mask_3d.dtype}")
if mask_3d.dtype.byteorder == '>':
mask_3d = mask_3d.astype(mask_3d.dtype.newbyteorder('='))
if output_anndata is None:
base_name = os.path.splitext(db_path)[0]
output_anndata = f"{base_name}.h5ad"
print(f"Auto-setting output_anndata: {output_anndata}")
conn, volume_results = build_interscellar_volume_database_from_neighbors(
mask_3d=mask_3d,
neighbor_pairs_csv=neighbor_pairs_csv,
neighbor_db_path=neighbor_db_path, # Use database if available (more efficient)
global_surface_pickle=global_surface_pickle,
halo_bboxes_pickle=halo_bboxes_pickle,
voxel_size_um=voxel_size_um,
db_path=db_path,
output_csv=output_csv,
output_anndata=output_anndata,
output_mesh_zarr=output_mesh_zarr,
max_distance_um=max_distance_um,
intracellular_threshold_um=intracellular_threshold_um,
n_jobs=n_jobs,
intermediate_results_dir=intermediate_results_dir,
resume=resume,
)
print(f"\n3. Verifying mesh zarr completion...")
import zarr
mesh_zarr_exists = False
if output_mesh_zarr and os.path.exists(output_mesh_zarr) and os.path.isdir(output_mesh_zarr):
try:
zarr_group = zarr.open(output_mesh_zarr, mode='r')
if 'interscellar_meshes' in zarr_group:
final_pairs = zarr_group.attrs.get('num_pairs', 0)
zarr_shape = zarr_group['interscellar_meshes'].shape
max_pair_id = np.asarray(zarr_group['interscellar_meshes']).max()
print(f"Mesh zarr verified: {output_mesh_zarr}")
print(f"Shape: {zarr_shape}")
print(f"Total pairs written: {final_pairs}")
print(f"Max pair ID: {max_pair_id}")
mesh_zarr_exists = True
else:
print(f"Mesh zarr exists but missing 'interscellar_meshes' key")
except Exception as e:
print(f"Error verifying mesh zarr: {e}")
else:
print(f"Mesh zarr not found at: {output_mesh_zarr}")
if output_cell_only_zarr:
print(f"\n4. Creating cell-only volumes zarr...")
if not mesh_zarr_exists:
print(f"Warning: Interscellar mesh zarr not available. Skipping cell-only zarr creation.")
print(f"Run the cell-only zarr creation separately after the interscellar zarr is ready.")
else:
try:
create_global_cell_only_volumes_zarr(
original_segmentation_zarr=ome_zarr_path,
interscellar_volumes_zarr=output_mesh_zarr,
output_zarr_path=output_cell_only_zarr
)
if os.path.exists(output_cell_only_zarr):
print(f"Cell-only volumes zarr created and verified: {output_cell_only_zarr}")
else:
print(f"Warning: Cell-only zarr creation reported success but file not found")
except Exception as e:
print(f"Error: Failed to create cell-only volumes zarr: {e}")
print(f"Other outputs (CSV, DB, Zarr) still available")
import traceback
traceback.print_exc()
step2_time = time.time() - step2_start
print(f"Interscellar volumes computed successfully")
print(f"Step 2 completed in {step2_time:.2f} seconds")
except Exception as e:
raise RuntimeError(f"Error in interscellar volume computation pipeline: {e}")
print(f"\n5. Retrieving results")
step3_start = time.time()
volume_results_df = None
if output_csv:
try:
volume_results_df = pd.read_sql_query("SELECT * FROM interscellar_volumes", conn)
print(f"Volume results table: {len(volume_results_df)} pairs")
except Exception as e:
print(f"Warning: Could not retrieve volume results table: {e}")
adata = None
if output_anndata:
if os.path.exists(output_anndata):
try:
if ANNDATA_AVAILABLE:
try:
import anndata as ad
adata = ad.read_h5ad(output_anndata)
print(f"AnnData file verified: {output_anndata}")
print(f"Shape: {adata.shape}")
print(f"Weighted adjacency matrix with interscellar volumes")
print(f"Component volumes in layers: edt_volume, intracellular_volume, touching_surface_area")
except Exception as e:
print(f"Warning: AnnData file exists but could not be loaded: {e}")
else:
print(f"AnnData file created: {output_anndata}")
print(f"Warning: AnnData package not available for verification (install with: pip install anndata)")
except Exception as e:
print(f"Warning: Could not verify AnnData file: {e}")
else:
print(f"Warning: AnnData file not found at: {output_anndata}")
step3_time = time.time() - step3_start
print(f"Step 3 completed in {step3_time:.2f} seconds")
print(f"\n6. Exporting to DuckDB format...")
step4_start = time.time()
try:
duckdb_output = db_path.replace('.db', '.duckdb')
export_interscellar_volumes_to_duckdb(conn, duckdb_output)
step4_time = time.time() - step4_start
print(f"DuckDB export completed in {step4_time:.2f} seconds")
except Exception as e:
print(f"Warning: DuckDB export failed: {e}")
overall_time = time.time() - overall_start_time
print(f"\n7. Pipeline completed successfully!")
print(f"Total execution time: {overall_time:.2f} seconds")
print(f"Database: {db_path}")
print(f"DuckDB: {db_path.replace('.db', '.duckdb')}")
if output_csv:
print(f"CSV output: {output_csv}")
if output_anndata:
print(f"AnnData output: {output_anndata}")
print(f"Interscellar volumes zarr: {output_mesh_zarr}")
print(f"Cell-only volumes zarr: {output_cell_only_zarr}")
try:
from ..core.compute_interscellar_volumes_3d_absolute import _cleanup_intermediate_results
_cleanup_intermediate_results(intermediate_results_dir)
except Exception as e:
print(f"Warning: Could not clean up intermediate results: {e}")
print(f"Intermediate results directory: {intermediate_results_dir}")
print("=" * 60)
if return_connection:
return volume_results_df, adata, conn
else:
conn.close()
return volume_results_df, adata, None
[docs]
def compute_interscellar_volumes_3d_adaptive(
ome_zarr_path: str,
neighbor_pairs_csv: str,
voxel_size_um: tuple = (0.56, 0.28, 0.28),
max_distance_um: float = 3.0,
surface_distance_um: float = 0.5,
rho_threshold: float = 0.5,
contact_rim_um: float = 0.0,
max_inward_um: Optional[float] = None,
output_csv: Optional[str] = None,
output_mesh_zarr: Optional[str] = None,
rejected_csv: Optional[str] = None,
output_dir: Optional[str] = None,
output_name_tag: str = "adaptive",
global_surface_pickle: Optional[str] = None,
halo_bboxes_pickle: Optional[str] = None,
exclude_truncated: bool = False,
reject_unbridged: bool = True,
preview: bool = True,
chunk_size: int = 500,
workers_load_volume: bool = False,
resume: bool = False,
n_jobs: int = 1
) -> Optional[pd.DataFrame]:
"""Compute interscellar volumes with the adaptive pathway.
Each volume is the extracellular corridor between two neighboring cells plus the
intracellular rind reaching into each cell. Unlike compute_interscellar_volumes_3d,
which reaches a fixed distance into each cell, the adaptive pathway dilates inwards
from the cell surface by a depth ratio relative to the cell's maximum reach.
Args:
ome_zarr_path: 3D segmentation label volume (.zarr or .npy).
neighbor_pairs_csv: Neighbor pairs CSV or neighbor graph .db from
find_cell_neighbors_3d.
voxel_size_um: (z, y, x) voxel size in micrometers.
max_distance_um: Maximum surface-to-surface distance the corridor may span.
surface_distance_um: Distance from the surface that defines each cell's rim.
rho_threshold: Depth ratio, relative to the cell's maximum reach, that the
volume extends into each cell.
contact_rim_um: Extra rim width around direct-contact interfaces.
max_inward_um: Optional cap on how far the volume reaches into a cell.
output_csv: Per-pair volumes CSV. Defaults to <stem>_adaptive_volumes.csv.
output_mesh_zarr: Pair-labeled volumes zarr. Defaults to
<stem>_adaptive_interscellar_volumes.zarr.
rejected_csv: CSV of rejected pairs. Defaults to
<stem>_adaptive_rejected_pairs.csv.
output_dir: Directory for auto-named outputs. Defaults to the directory of
neighbor_pairs_csv.
output_name_tag: Tag inserted into auto-named outputs.
global_surface_pickle: Optional surfaces .pkl; must match the mask shape.
halo_bboxes_pickle: Optional halo bounding boxes .pkl; must lie inside the mask.
exclude_truncated: Skip pairs where either cell touches a face of the volume.
reject_unbridged: Reject pairs with no extracellular corridor reaching both
cells. Disable only for comparison.
preview: Write the dense interscellar_meshes/overlap_count arrays alongside the
lossless per-pair archive.
chunk_size: Pairs processed per chunk.
workers_load_volume: Each worker loads the whole volume instead of reading
pair crops lazily.
resume: Skip pairs already in output_csv and append to existing outputs.
n_jobs: Number of worker processes.
Returns:
DataFrame of per-pair volume measurements read back from output_csv.
Note:
Interscellar volumes may overlap; a voxel shared by several pairs belongs to
each of them in the per-pair archive, which calculate_interscellar_scores_3d
reads. The dense interscellar_meshes array can show only one pair per voxel.
"""
print("=" * 60)
print("InterSCellar: Adaptive Volume Computation - 3D")
print("=" * 60)
from ..core.compute_interscellar_volumes_3d_adaptive import (
main as _adaptive_main,
default_output_paths
)
if output_csv is None:
output_csv = default_output_paths(neighbor_pairs_csv, output_name_tag, output_dir)[0]
z, y, x = voxel_size_um
argv = [
"--mask", ome_zarr_path,
"--pairs", neighbor_pairs_csv,
"--out", output_csv,
"--name-tag", output_name_tag,
"--max-distance-um", str(max_distance_um),
"--surface-distance-um", str(surface_distance_um),
"--rho-threshold", str(rho_threshold),
"--contact-rim-um", str(contact_rim_um),
"--z", str(z), "--y", str(y), "--x", str(x),
"--chunk-size", str(chunk_size),
"--preview", "both" if preview else "none",
"--n-jobs", str(n_jobs),
]
optional = {
"--max-inward-um": max_inward_um,
"--out-zarr": output_mesh_zarr,
"--rejected-out": rejected_csv,
"--output-dir": output_dir,
"--global-surface": global_surface_pickle,
"--halo-bboxes": halo_bboxes_pickle,
}
for flag, value in optional.items():
if value is not None:
argv += [flag, str(value)]
flags = {
"--exclude-truncated": exclude_truncated,
"--no-reject-unbridged": not reject_unbridged,
"--workers-load-volume": workers_load_volume,
"--resume": resume,
}
argv += [flag for flag, enabled in flags.items() if enabled]
# The pipeline reports invalid inputs through argparse, which exits; surface that
# as an exception so a notebook or calling script is not terminated.
try:
_adaptive_main(argv)
except SystemExit as exc:
if exc.code not in (None, 0):
raise RuntimeError(
"Adaptive interscellar volume computation failed (see the error above)"
) from None
volume_results_df = None
if os.path.exists(output_csv):
volume_results_df = pd.read_csv(output_csv)
print(f"Volume results table: {len(volume_results_df)} pairs")
print("=" * 60)
return volume_results_df
[docs]
def compute_cell_only_volumes_3d(
ome_zarr_path: str,
interscellar_volumes_zarr: str,
output_zarr_path: Optional[str] = None,
neighbor_db_path: Optional[str] = None,
nuclei_segmentation_zarr: Optional[str] = None,
slab: int = 16,
n_jobs: int = 1
) -> pd.DataFrame:
"""Compute cell-only volumes: each cell minus every interscellar volume.
A voxel keeps its cell label unless it lies in any interscellar volume, or, when
nuclei_segmentation_zarr is given, in any nucleus. Works with the output of either
pathway: an absolute store (label array '0') or an adaptive store (the
interscellar_meshes grid, or the lossless pair archive when the store was written
with preview=False). The volume is processed in Z-slabs, so memory stays at a few
slabs per worker regardless of dataset size.
Args:
ome_zarr_path: Cell segmentation zarr (cell IDs as labels).
interscellar_volumes_zarr: Interscellar volumes zarr from
compute_interscellar_volumes_3d or compute_interscellar_volumes_3d_adaptive.
output_zarr_path: Output zarr. Defaults to the interscellar zarr's name with
_interscellar_volumes replaced by _cell_only_volumes (plus _nuclei_excluded
when nuclei are given), in the same directory.
neighbor_db_path: Optional neighbor graph .db; supplies the cell_type column.
nuclei_segmentation_zarr: Optional nuclei segmentation zarr on the same grid;
every voxel inside a nucleus is also removed.
slab: Z slices processed at a time.
n_jobs: Number of worker processes, each handling one Z-slab at a time.
Returns:
DataFrame with one row per cell that keeps any voxels: cell_id,
cell_only_volume_voxels, cell_only_volume_um3 and cell_type.
Note:
Only the union of the interscellar volumes is subtracted, so overlapping
adaptive volumes remove each shared voxel once. Cells lying entirely inside
interscellar volumes (or nuclei) have no row. Volumes in um3 use the voxel size
recorded in the interscellar zarr, or (0.56, 0.28, 0.28) when it has none.
"""
print("=" * 60)
print("InterSCellar: Cell-Only Volumes Computation - 3D")
print("=" * 60)
from ..core.compute_cell_only_volumes_3d import compute_cell_only_volumes
summary = compute_cell_only_volumes(
cell_segmentation_zarr=ome_zarr_path,
interscellar_zarr=interscellar_volumes_zarr,
output_zarr_path=output_zarr_path,
slab=slab,
nuclei_segmentation_zarr=nuclei_segmentation_zarr,
n_jobs=n_jobs,
)
voxel_size_um = summary["voxel_size_um"]
if voxel_size_um is None:
voxel_size_um = (0.56, 0.28, 0.28)
print(f"Warning: voxel_size_um not found in the interscellar zarr, using default: {voxel_size_um}")
voxel_volume_um3 = float(np.prod(voxel_size_um))
cell_types = {}
if neighbor_db_path and os.path.exists(neighbor_db_path):
try:
conn = sqlite3.connect(neighbor_db_path)
try:
cells_df = pd.read_sql_query("SELECT cell_id, cell_type FROM cells", conn)
cell_types = dict(zip(cells_df['cell_id'], cells_df['cell_type']))
print(f"Loaded cell types for {len(cell_types)} cells from neighbor database")
except Exception as e:
print(f"Warning: Could not load cell types from database: {e}")
finally:
conn.close()
except Exception as e:
print(f"Warning: Could not access neighbor database: {e}")
remaining = summary["cell_voxels_after"]
cell_ids = np.asarray(sorted(remaining), dtype=np.int64)
voxels = np.asarray([remaining[c] for c in cell_ids.tolist()], dtype=np.int64)
cell_only_df = pd.DataFrame({
'cell_id': cell_ids.astype(int),
'cell_only_volume_voxels': voxels.astype(int),
'cell_only_volume_um3': (voxels * voxel_volume_um3).astype(float),
'cell_type': [cell_types.get(int(cid), 'unknown') for cid in cell_ids],
})
print(f"Computed cell-only volumes for {len(cell_only_df)} cells")
if len(cell_only_df):
print(f"Total cell-only volume: {cell_only_df['cell_only_volume_um3'].sum():.2f} μm³")
print(f"Mean cell-only volume: {cell_only_df['cell_only_volume_um3'].mean():.2f} μm³")
print("=" * 60)
return cell_only_df
[docs]
def calculate_interscellar_scores_3d(
interscellar_volumes_zarr: str,
spot_zarrs: Any,
biomarker: Optional[str] = None,
volumes_csv: Optional[str] = None,
output_csv: Optional[str] = None,
voxel_size_um: Optional[Tuple[float, float, float]] = None,
n_jobs: int = 1,
decay_power: float = 1.0,
reference_distance_um: Any = None,
pair_batch: int = 32,
per_point_csv: Optional[str] = None
) -> pd.DataFrame:
print("=" * 60)
print("InterSCellar: Interscellar Score Computation - 3D")
print("=" * 60)
from ..core.calculate_interscellar_scores_3d_point import (
calculate_interscellar_scores_3d as _calculate_scores
)
written_csv = _calculate_scores(
interscellar_zarr=interscellar_volumes_zarr,
spot_zarrs=spot_zarrs,
biomarker=biomarker,
volumes_csv=volumes_csv,
output_csv=output_csv,
voxel_size_um=voxel_size_um,
n_jobs=n_jobs,
decay_power=decay_power,
reference_distance_um=reference_distance_um,
pair_batch=pair_batch,
per_point_csv=per_point_csv
)
scores_df = pd.read_csv(written_csv)
print("=" * 60)
return scores_df