OSCR

Endothelial <i>Adgrl2</i> Expression and Alternative Splicing Controls the Cerebrovasculature.

Code ↔ Paper

2 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 2 matches
  1. [1] § Materials and Methods › Quantification and statistical analysis › Vascular morphology analysis ↔ adaptive_vascular_analysis.py, lines 356–421 · score 0.90 · vessel mask, vessel diameter, Frangi filter, vessel area, Adaptive, Equalization
  2. [2] § Results › Endothelial Adgrl2 alternative splicing program prevents ectopic synaptogenesis ↔ adaptive_vascular_analysis.py, lines 226–294 · score 0.54 · vessel widths, point densities, branch points, CD31, vasculature

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 authors' code

Python · 785 lines · 42 KB · CC-BY-4.0 · 2 matches

  1. import os
  2. import geojson
  3. from shapely.geometry import shape
  4. from shapely import affinity
  5. import rasterio
  6. from rasterio.transform import from_origin
  7. from rasterio.features import rasterize
  8. from rasterio.transform import Affine # Import Affine directly
  9. import numpy as np
  10. import scipy.ndimage as ndi
  11. from skimage import filters, morphology, measure, segmentation, io, restoration
  12. from skimage.filters import frangi, threshold_local
  13. from skimage.feature import peak_local_max
  14. from skimage.color import label2rgb
  15. from skimage.exposure import equalize_adapthist, rescale_intensity
  16. import matplotlib.pyplot as plt
  17. import pandas as pd
  18. from aicsimageio import AICSImage # This library works for TIFF, CZI, and more!
  19. import warnings
  20. import time
  21. import itertools # To iterate through parameter combinations
  22. import concurrent.futures # **** ADD FOR PARALLELISM ****
  23. import traceback # For detailed error logging in parallel tasks
  26. # --- Configuration ---
  27. # Input/Output Folders
  28. DATA_FOLDER = 'images'
  29. OUTPUT_CSV_FILENAME = 'choroid_plexus_analysis_param_sweep_um_params.csv'
  30. OUTPUT_PLOT_FOLDER = os.path.join(DATA_FOLDER, 'analysis_plots_param_sweep_um_params')
  32. # --- Global Physical Constant ---
  33. # This is the assumed physical size of a pixel if not found in image metadata.
  34. # ALL micrometer (_UM) parameters below are converted to pixels using this value.
  35. ASSUMED_PIXEL_SIZE_UM = 0.638148 # µm/pixel
  36. FORCE_ASSUMED_PIXEL_SIZE = True
  38. SAVE_PLOTS = True
  39. if SAVE_PLOTS and not os.path.exists(OUTPUT_PLOT_FOLDER):
  40. os.makedirs(OUTPUT_PLOT_FOLDER)
  42. # Channel indices
  43. DAPI_CHANNEL_INDEX = 0
  44. CD31_CHANNEL_INDEX = 0
  46. # --- Preprocessing Parameters (ALL SPATIAL PARAMS NOW IN MICROMETERS) ---
  47. # Note: These values are converted to pixels internally for image processing functions.
  48. CD31_MEDIAN_RADIUS_UM = 0.22 # In µm. Original: 1 pixel -> 1 * 0.2196 ≈ 0.22 µm
  49. CD31_ROLLING_BALL_RADIUS = 0 # Disabled
  50. DAPI_MEDIAN_RADIUS_UM = 0.44 # In µm. Original: 2 pixels -> 2 * 0.2196 ≈ 0.44 µm
  52. APPLY_CLAHE_TO_CD31 = True
  53. CLAHE_KERNEL_SIZE_DIVISOR = 8 # Unitless: divides image dims to get tile size.
  54. CLAHE_CLIP_LIMIT = 0.01 # Unitless: clipping limit for contrast.
  56. # --- Segmentation Parameters (ALL SPATIAL PARAMS NOW IN MICROMETERS) ---
  57. NUCLEI_THRESHOLD_METHOD = 'otsu'
  58. NUCLEI_MIN_DISTANCE_PEAKS_UM = 0.44 # In µm. Original: 2 pixels -> 2 * 0.2196 ≈ 0.44 µm
  60. # --- VESSELNESS SIGMA RANGES TO TEST (DEFINED IN MICROMETERS) ---
  61. # Each element is a tuple: (start_sigma_um, stop_sigma_um, step_sigma_um)
  62. # Original pixel range(20, 70, 1) -> (20*0.2196, 70*0.2196, 1*0.2196) -> (4.39, 15.37, 0.22)
  63. VESSELNESS_SIGMA_RANGES_TO_TEST_UM = [
  64. (4.4, 15.4, 0.22)
  65. # e.g., (2.0, 8.0, 0.5) to test sigmas from 2µm to 8µm in 0.5µm steps.
  66. ]
  68. # --- Post-processing Filter Parameters (ALREADY IN MICROMETERS) ---
  69. MIN_VESSEL_WIDTH_UM = 2.0 # In µm. Remove vessels narrower than this.
  70. MIN_VESSEL_SEGMENT_LENGTH_UM = 20.0 # In µm. Remove skeleton segments shorter than this.
  71. MIN_BRANCH_POINT_DISTANCE_UM = 10 # <<<--- Minimum distance between branch points in micrometers.
  73. # --- PARAMETER SWEEP DEFINITION (Adaptive Threshold in MICROMETERS) ---
  74. # Original pixel block_size 801 -> 801 * 0.2196 ≈ 175.9 µm. Let's use 176 µm.
  75. BLOCK_SIZES_TO_TEST_UM = [176.0]
  76. OFFSETS_TO_TEST = [-0.0001] # This is a unitless intensity offset, not a spatial parameter.
  77. ADAPTIVE_METHOD = 'gaussian'
  79. # --- PARALLELISM CONFIGURATION ---
  80. NUM_WORKERS = os.cpu_count() - 2 # A slightly more conservative default
  81. if NUM_WORKERS < 1:
  82. NUM_WORKERS = 1
  84. # --- Configuration Printout (Updated for Micrometer Parameters) ---
  85. print(f"--- Configuration ---")
  86. print(f"ASSUMED PIXEL SIZE: {ASSUMED_PIXEL_SIZE_UM} µm/pixel (used for all conversions)")
  87. print(f"Data Folder: {DATA_FOLDER}")
  88. print(f"Output CSV: {OUTPUT_CSV_FILENAME}")
  89. print(f"Save Plots: {SAVE_PLOTS}")
  90. if SAVE_PLOTS: print(f"Plot Folder: {OUTPUT_PLOT_FOLDER}")
  91. print(f"DAPI Channel: {DAPI_CHANNEL_INDEX}, CD31 Channel: {CD31_CHANNEL_INDEX}")
  92. print(f"CD31 Median Radius: {CD31_MEDIAN_RADIUS_UM} µm")
  93. print(f"DAPI Median Radius: {DAPI_MEDIAN_RADIUS_UM} µm")
  94. print(f"Apply CLAHE to CD31: {APPLY_CLAHE_TO_CD31}")
  95. if APPLY_CLAHE_TO_CD31:
  96. print(f" CLAHE Kernel Size Divisor: {CLAHE_KERNEL_SIZE_DIVISOR}")
  97. print(f" CLAHE Clip Limit: {CLAHE_CLIP_LIMIT}")
  98. print(f"Nuclei Threshold: {NUCLEI_THRESHOLD_METHOD}, Min Distance: {NUCLEI_MIN_DISTANCE_PEAKS_UM} µm")
  100. print(f"Vesselness Sigma Ranges to Test (in µm):")
  101. for i, (start, stop, step) in enumerate(VESSELNESS_SIGMA_RANGES_TO_TEST_UM):
  102. print(f" - Config {i+1}: start={start}µm, stop={stop}µm, step={step}µm")
  104. print(f"Minimum Vessel Width Filter: {MIN_VESSEL_WIDTH_UM} µm")
  105. print(f"Minimum Vessel Segment Length: {MIN_VESSEL_SEGMENT_LENGTH_UM} µm")
  106. print(f"--- Parameter Sweep Configuration (Adaptive Threshold) ---")
  107. print(f"Block Sizes to Test: {BLOCK_SIZES_TO_TEST_UM} µm")
  108. print(f"Offsets to Test: {OFFSETS_TO_TEST}")
  109. print(f"Adaptive Method: {ADAPTIVE_METHOD}")
  110. print(f"--- Parallelism Configuration ---")
  111. print(f"Using up to {NUM_WORKERS} worker processes.")
  112. print(f"---------------------------------")
  115. # --- Helper Functions ---
  116. ### CHANGED ### - Renamed function from load_czi_image to load_image
  117. def load_image(path):
  118. """Loads a bio-image (TIFF, CZI, etc.) using aicsimageio and extracts pixel size."""
  119. try:
  120. img_aics = AICSImage(path)
  121. dims_order = img_aics.dims.order
  122. has_c = 'C' in dims_order
  123. has_y = 'Y' in dims_order
  124. has_x = 'X' in dims_order
  126. if has_c and has_y and has_x:
  127. image_stack = img_aics.get_image_data("CYX", T=0, Z=0)
  128. elif has_y and has_x:
  129. print(f"Warning: Image {os.path.basename(path)} seems 2D (missing C dimension). Loading as YX.")
  130. image_stack = img_aics.get_image_data("YX", T=0, Z=0)
  131. else:
  132. print(f"Error: Image {os.path.basename(path)} missing Y or X dimension. Cannot process.")
  133. return None, None
  135. pixel_size_x_um = img_aics.physical_pixel_sizes.X
  136. pixel_size_y_um = img_aics.physical_pixel_sizes.Y
  137. # Use the global constant as the fallback
  138. provided_pixel_size = ASSUMED_PIXEL_SIZE_UM
  140. if pixel_size_x_um is None or pixel_size_y_um is None:
  141. # ### MODIFIED ### - Added 'TIFF' to warning for clarity. This is the expected path for many TIFFs.
  142. print(f"Warning: Could not read pixel size from image metadata (e.g., common for standard TIFF files). Using assumed value: {provided_pixel_size} µm.")
  143. pixel_size_xy_um = provided_pixel_size
  144. elif not np.isclose(pixel_size_x_um, pixel_size_y_um):
  145. print(f"Warning: Anisotropic pixel size in {os.path.basename(path)} metadata (X: {pixel_size_x_um}, Y: {pixel_size_y_um}). Using assumed value: {provided_pixel_size} µm.")
  146. pixel_size_xy_um = provided_pixel_size
  147. else:
  148. pixel_size_xy_um = pixel_size_x_um
  150. print(f"Image loaded: {os.path.basename(path)}, shape={image_stack.shape}, dtype={image_stack.dtype}, using pixel_size={pixel_size_xy_um:.4f} µm")
  151. return image_stack, pixel_size_xy_um
  152. except FileNotFoundError:
  153. print(f"Error: Image file not found at {path}")
  154. return None, None
  155. except Exception as e:
  156. print(f"Error loading image {path}: {e}")
  157. traceback.print_exc()
  158. return None, None
  160. def load_roi_polygon(geojson_path):
  161. """Loads the first Polygon feature from a GeoJSON file."""
  162. try:
  163. with open(geojson_path) as f:
  164. gj = geojson.load(f)
  165. roi_polygon = None
  166. if isinstance(gj, geojson.FeatureCollection):
  167. for feature in gj.features:
  168. if isinstance(feature.geometry, geojson.Polygon):
  169. roi_polygon = shape(feature.geometry)
  170. break
  171. elif isinstance(gj, geojson.Feature):
  172. if isinstance(gj.geometry, geojson.Polygon):
  173. roi_polygon = shape(gj.geometry)
  174. elif isinstance(gj, geojson.Polygon):
  175. roi_polygon = shape(gj)
  177. if roi_polygon:
  178. if not roi_polygon.is_valid:
  179. print(f"Warning: Loaded polygon from {geojson_path} is invalid, attempting to buffer(0).")
  180. roi_polygon = roi_polygon.buffer(0)
  181. if not roi_polygon.is_valid or roi_polygon.is_empty:
  182. print(f"Error: Polygon from {geojson_path} remains invalid or empty after buffer(0).")
  183. return None
  184. return roi_polygon
  185. else:
  186. print(f"Error: No Polygon geometry found in {geojson_path}")
  187. return None
  188. except FileNotFoundError:
  189. print(f"Error: GeoJSON file not found at {geojson_path}")
  190. return None
  191. except Exception as e:
  192. print(f"Error loading GeoJSON {geojson_path}: {e}")
  193. return None
  195. def create_roi_mask(image_shape_2d, roi_polygon, pixel_size_um):
  196. """Creates a 2D boolean mask from a Shapely polygon in µm coordinates."""
  197. if roi_polygon is None:
  198. print("Error: create_roi_mask called with None polygon.")
  199. return None
  200. # This transform correctly maps micron-space coordinates of the polygon to pixel-space of the image.
  201. transform = Affine(pixel_size_um, 0.0, 0.0, 0.0, pixel_size_um, 0.0)
  202. try:
  203. mask = rasterize(
  204. [(roi_polygon, 1)],
  205. out_shape=image_shape_2d,
  206. transform=transform,
  207. fill=0,
  208. dtype=np.uint8,
  209. all_touched=False
  210. )
  211. if not np.any(mask):
  212. img_extent_x_um = image_shape_2d[1] * pixel_size_um
  213. img_extent_y_um = image_shape_2d[0] * pixel_size_um
  214. print(f"Warning: Rasterization produced an empty mask for polygon.")
  215. print(f"-> Polygon bounds (µm): {roi_polygon.bounds}")
  216. print(f"-> Image expected extent (µm): X=[0, {img_extent_x_um:.2f}], Y=[0, {img_extent_y_um:.2f}]")
  217. print(f"-> Check GeoJSON coordinates, their scaling to µm, and pixel size ({pixel_size_um:.4f} µm).")
  218. except Exception as e:
  219. print(f"Error during rasterization: {e}")
  220. traceback.print_exc()
  221. raise
  222. return mask.astype(bool)
  226. def analyze_image_roi(image_path, geojson_path, current_sigma_range_um, adaptive_block_size_um, adaptive_offset):
  227. """
  228. Performs the full analysis pipeline for a single image/GeoJSON pair
  229. using specified parameters defined in micrometers.
  230. """
  231. # (The top part of the function remains the same)
  232. (param_sigma_start_um, param_sigma_stop_um, param_sigma_step_um) = current_sigma_range_um
  233. results = {
  234. 'filename_image': os.path.basename(image_path),
  235. 'filename_geojson': os.path.basename(geojson_path) if geojson_path else 'N/A (Full Image)',
  236. 'vesselness_sigma_start_um': param_sigma_start_um, 'vesselness_sigma_stop_um': param_sigma_stop_um, 'vesselness_sigma_step_um': param_sigma_step_um,
  237. 'adaptive_block_size_um': adaptive_block_size_um, 'adaptive_offset': adaptive_offset,
  238. 'min_vessel_width_um_param': MIN_VESSEL_WIDTH_UM, 'min_vessel_segment_length_um_param': MIN_VESSEL_SEGMENT_LENGTH_UM,
  239. 'pixel_size_um': None, 'roi_area_pixels': None, 'roi_area_um2': None, 'roi_area_mm2': None,
  240. 'vessel_area_pixels': None, 'vessel_area_um2': None, 'vaf_percent': None,
  241. 'mean_vessel_diameter_um': None, 'raw_skeleton_pixels': None, 'raw_skeleton_length_um': None,
  242. 'filtered_skeleton_pixels': None, 'filtered_skeleton_length_um': None, 'segments_analyzed': None, 'segments_kept': None,
  243. 'filtered_branch_points_count': None, 'filtered_branch_point_density_mm2': None, 'nuclei_count': None, 'nuclear_density_mm2': None,
  244. 'error_message': None
  245. }
  246. run_id_params = (f"Sigmas={param_sigma_start_um}-{param_sigma_stop_um}s{param_sigma_step_um}um, "
  247. f"Block={adaptive_block_size_um}um, Offset={adaptive_offset}, "
  248. f"MinWidth={MIN_VESSEL_WIDTH_UM}µm, MinLen={MIN_VESSEL_SEGMENT_LENGTH_UM}µm")
  249. run_id = f"{os.path.basename(image_path)} ({run_id_params})"
  250. print(f"--- Starting Task: {run_id} ---")
  251. start_time = time.time()
  253. try:
  254. image_stack, pixel_size_um = load_image(image_path)
  255. if image_stack is None or pixel_size_um is None or pixel_size_um <= 0:
  256. raise ValueError(f"Failed to load image or determine valid pixel size for {image_path}")
  258. if image_stack.ndim == 3 and image_stack.shape[0] == 1:
  259. print(f" -> Info for {run_id}: Image has shape (1, Y, X). Squeezing to 2D.")
  260. image_stack = image_stack.squeeze(axis=0)
  261. # ################################################
  263. results['pixel_size_um'] = pixel_size_um
  264. AREA_UNIT_FACTOR = pixel_size_um * pixel_size_um
  265. LENGTH_UNIT_FACTOR = pixel_size_um
  266. DENSITY_AREA_MM2 = 1e6
  268. cd31_median_radius_px = int(round(CD31_MEDIAN_RADIUS_UM / pixel_size_um))
  269. dapi_median_radius_px = int(round(DAPI_MEDIAN_RADIUS_UM / pixel_size_um))
  270. nuclei_min_distance_peaks_px = int(round(NUCLEI_MIN_DISTANCE_PEAKS_UM / pixel_size_um))
  271. if nuclei_min_distance_peaks_px < 1: nuclei_min_distance_peaks_px = 1
  272. adaptive_block_size_px = int(round(adaptive_block_size_um / pixel_size_um))
  273. if adaptive_block_size_px % 2 == 0: adaptive_block_size_px += 1
  274. if adaptive_block_size_px < 3:
  275. print(f"Warning for {run_id}: Calculated block size ({adaptive_block_size_px}px) is too small. Setting to 3.")
  276. adaptive_block_size_px = 3
  277. sigma_start_px = int(round(param_sigma_start_um / pixel_size_um))
  278. sigma_stop_px = int(round(param_sigma_stop_um / pixel_size_um))
  279. sigma_step_px = int(round(param_sigma_step_um / pixel_size_um))
  280. if sigma_start_px < 1: sigma_start_px = 1
  281. if sigma_stop_px <= sigma_start_px: sigma_stop_px = sigma_start_px + 1
  282. if sigma_step_px < 1: sigma_step_px = 1
  283. frangi_sigma_range_px = range(sigma_start_px, sigma_stop_px, sigma_step_px)
  284. if not list(frangi_sigma_range_px):
  285. error_msg = (f"Vesselness sigma range in pixels is empty. "
  286. f"Params(um): start={param_sigma_start_um}, stop={param_sigma_stop_um}, step={param_sigma_step_um} "
  287. f"-> Converted(px): start={sigma_start_px}, stop={sigma_stop_px}, step={sigma_step_px}")
  288. raise ValueError(error_msg)
  289. print(f" -> Converted Params for {run_id}:")
  290. print(f" CD31 Median Radius: {CD31_MEDIAN_RADIUS_UM}µm -> {cd31_median_radius_px}px")
  291. print(f" DAPI Median Radius: {DAPI_MEDIAN_RADIUS_UM}µm -> {dapi_median_radius_px}px")
  292. print(f" Nuclei Min Dist: {NUCLEI_MIN_DISTANCE_PEAKS_UM}µm -> {nuclei_min_distance_peaks_px}px")
  293. print(f" Adaptive Block Size: {adaptive_block_size_um}µm -> {adaptive_block_size_px}px")
  294. print(f" Frangi Sigmas: {current_sigma_range_um}µm -> range({frangi_sigma_range_px.start}, {frangi_sigma_range_px.stop}, {frangi_sigma_range_px.step})px")
  296. img_shape = image_stack.shape
  297. dapi_img, cd31_img = None, None
  298. img_shape_2d = None
  299. if len(img_shape) == 3 and img_shape[0] >= max(DAPI_CHANNEL_INDEX, CD31_CHANNEL_INDEX) + 1:
  300. dapi_img = image_stack[DAPI_CHANNEL_INDEX, :, :]
  301. cd31_img = image_stack[CD31_CHANNEL_INDEX, :, :]
  302. img_shape_2d = (img_shape[1], img_shape[2])
  303. elif len(img_shape) == 2:
  304. print(f"Info for {run_id}: Handling as 2D image (shape {img_shape}). Assigning channels based on index.")
  305. img_shape_2d = img_shape
  306. if CD31_CHANNEL_INDEX == 0 and DAPI_CHANNEL_INDEX != 0:
  307. cd31_img = image_stack
  308. elif DAPI_CHANNEL_INDEX == 0 and CD31_CHANNEL_INDEX != 0:
  309. dapi_img = image_stack
  310. else:
  311. # Default or ambiguous case: assume it's the primary channel (CD31)
  312. cd31_img = image_stack
  313. if DAPI_CHANNEL_INDEX == CD31_CHANNEL_INDEX:
  314. print(f"Warning for {run_id}: DAPI and CD31 channels are the same. Only CD31 will be processed from the 2D image.")
  315. else:
  316. print(f"Info for {run_id}: Assuming 2D image is CD31 channel (index {CD31_CHANNEL_INDEX}).")
  317. else:
  318. raise ValueError(f"Image {run_id}: Unexpected image dimensions: {img_shape}.")
  322. if cd31_img is None: raise ValueError(f"Image {run_id}: CD31 image could not be extracted. Check channel indices.")
  323. if img_shape_2d is None: raise ValueError(f"Image {run_id}: Could not determine 2D shape.")
  325. roi_mask = None
  326. if geojson_path:
  327. print(f" -> Loading ROI from {os.path.basename(geojson_path)}...")
  328. roi_polygon_pixels = load_roi_polygon(geojson_path)
  329. if roi_polygon_pixels is None:
  330. raise ValueError(f"Failed to load ROI polygon from {geojson_path}")
  331. roi_polygon_um = affinity.scale(roi_polygon_pixels, xfact=pixel_size_um, yfact=pixel_size_um, origin=(0, 0))
  332. roi_mask = create_roi_mask(img_shape_2d, roi_polygon_um, pixel_size_um)
  333. if roi_mask is None: raise ValueError(f"Image {run_id}: Failed to create ROI mask from GeoJSON.")
  334. else:
  335. print(f" -> No GeoJSON file provided. Using entire image as the Region of Interest.")
  336. roi_mask = np.ones(img_shape_2d, dtype=bool)
  338. total_roi_area_pixels = np.sum(roi_mask)
  339. if total_roi_area_pixels == 0:
  340. raise ValueError(f"ROI mask for {run_id} is empty. Check GeoJSON coordinates and image size.")
  344. total_roi_area_um2 = total_roi_area_pixels * AREA_UNIT_FACTOR
  345. total_roi_area_mm2 = total_roi_area_um2 / DENSITY_AREA_MM2
  346. results.update({'roi_area_pixels': total_roi_area_pixels, 'roi_area_um2': total_roi_area_um2, 'roi_area_mm2': total_roi_area_mm2})
  348. cd31_for_filter = cd31_img.astype(np.float32)
  349. if cd31_median_radius_px > 0:
  350. cd31_for_filter = filters.median(cd31_for_filter, morphology.disk(cd31_median_radius_px))
  352. if APPLY_CLAHE_TO_CD31:
  353. actual_ntiles_for_clahe = CLAHE_KERNEL_SIZE_DIVISOR
  354. print(f" -> Applying CLAHE to CD31 for {run_id} (ntiles={actual_ntiles_for_clahe}, clip_limit={CLAHE_CLIP_LIMIT})")
  356. normalized_cd31_for_clahe = rescale_intensity(
  357. cd31_for_filter,
  358. in_range='image',
  359. out_range=(0.0, 1.0)
  360. )
  362. cd31_clahe_output = equalize_adapthist(
  363. normalized_cd31_for_clahe,
  364. kernel_size=actual_ntiles_for_clahe,
  365. clip_limit=CLAHE_CLIP_LIMIT,
  366. nbins=256
  367. )
  368. cd31_for_filter = cd31_clahe_output.astype(np.float32)
  369. print(f" -> CLAHE applied. Output range: {cd31_for_filter.min():.3f} - {cd31_for_filter.max():.3f}, dtype: {cd31_for_filter.dtype}")
  371. frangi_output_full = np.zeros_like(cd31_for_filter)
  372. try:
  373. with warnings.catch_warnings():
  374. warnings.simplefilter("ignore", category=UserWarning)
  375. frangi_output_full = frangi(cd31_for_filter.astype(np.float32), sigmas=frangi_sigma_range_px, black_ridges=False)
  376. except Exception as frangi_err:
  377. raise ValueError(f"Frangi filtering failed for {run_id}: {frangi_err}")
  379. vesselness_image_roi = np.zeros_like(frangi_output_full)
  380. vesselness_image_roi[roi_mask] = frangi_output_full[roi_mask]
  382. vessel_mask_initial = np.zeros_like(vesselness_image_roi, dtype=bool)
  383. if np.any(vesselness_image_roi > 1e-9):
  384. try:
  385. local_thresh_map = threshold_local(
  386. frangi_output_full,
  387. block_size=adaptive_block_size_px,
  388. method=ADAPTIVE_METHOD,
  389. offset=adaptive_offset
  390. )
  391. vessel_mask_initial = (frangi_output_full > local_thresh_map) & roi_mask
  392. if not np.any(vessel_mask_initial):
  393. print(f"Warning for {run_id}: Initial vessel mask is empty after adaptive thresholding.")
  394. except Exception as e:
  395. raise ValueError(f"Adaptive thresholding failed for {run_id}: {e}")
  396. else:
  397. print(f"Warning for {run_id}: No significant Frangi signal in ROI. Initial vessel mask empty.")
  400. mean_diameter_um = 0.0
  401. if np.any(vessel_mask_initial):
  402. distance_map_on_initial_mask_um = ndi.distance_transform_edt(vessel_mask_initial, sampling=[pixel_size_um, pixel_size_um])
  403. skeleton_of_initial_mask = morphology.skeletonize(vessel_mask_initial) & roi_mask
  404. if np.any(skeleton_of_initial_mask):
  405. radii_at_skeleton_points_um = distance_map_on_initial_mask_um[skeleton_of_initial_mask]
  406. valid_radii = radii_at_skeleton_points_um[radii_at_skeleton_points_um > 1e-6]
  407. if len(valid_radii) > 0:
  408. mean_diameter_um = 2 * np.mean(valid_radii)
  409. results['mean_vessel_diameter_um'] = mean_diameter_um
  410. vessel_mask = vessel_mask_initial.copy()
  411. if MIN_VESSEL_WIDTH_UM is not None and MIN_VESSEL_WIDTH_UM > 0 and pixel_size_um > 0:
  412. if np.any(vessel_mask):
  413. min_width_px = MIN_VESSEL_WIDTH_UM / pixel_size_um
  414. if min_width_px >= 1:
  415. radius_px_se = int(round((min_width_px - 1) / 2))
  416. if radius_px_se >= 1:
  417. vessel_mask = morphology.binary_opening(vessel_mask, footprint=morphology.disk(radius_px_se))
  418. vessel_area_pixels = np.sum(vessel_mask)
  419. vessel_area_um2 = vessel_area_pixels * AREA_UNIT_FACTOR
  420. vaf = (vessel_area_pixels / total_roi_area_pixels) * 100 if total_roi_area_pixels > 0 else 0
  421. results.update({'vessel_area_pixels': vessel_area_pixels, 'vessel_area_um2': vessel_area_um2, 'vaf_percent': vaf})
  425. raw_skeleton_pixels = 0; raw_skeleton_length_um = 0
  426. filtered_skeleton_pixels = 0; filtered_skeleton_length_um = 0
  427. segments_analyzed = 0; segments_kept = 0
  428. num_branch_points = 0; branch_point_density = 0
  429. skeleton_raw = None; skeleton_filtered = None; branch_points_mask = None
  431. if vessel_area_pixels > 5:
  432. skeleton_raw = morphology.skeletonize(vessel_mask) & roi_mask
  433. raw_skeleton_pixels = np.sum(skeleton_raw)
  434. raw_skeleton_length_um = raw_skeleton_pixels * LENGTH_UNIT_FACTOR
  435. results.update({'raw_skeleton_pixels': raw_skeleton_pixels, 'raw_skeleton_length_um': raw_skeleton_length_um})
  437. skeleton_filtered = skeleton_raw.copy()
  438. if MIN_VESSEL_SEGMENT_LENGTH_UM > 0:
  439. min_len_px = MIN_VESSEL_SEGMENT_LENGTH_UM / pixel_size_um
  440. if raw_skeleton_pixels > 0:
  441. labels, num_labels = measure.label(skeleton_raw, connectivity=2, return_num=True)
  442. segments_analyzed = num_labels
  443. kept_labels = [lbl for lbl in range(1, num_labels + 1) if np.sum(labels == lbl) >= min_len_px]
  444. segments_kept = len(kept_labels)
  445. if segments_kept < segments_analyzed:
  446. skeleton_filtered = skeleton_raw & np.isin(labels, kept_labels)
  448. filtered_skeleton_pixels = np.sum(skeleton_filtered)
  449. filtered_skeleton_length_um = filtered_skeleton_pixels * LENGTH_UNIT_FACTOR
  451. if filtered_skeleton_pixels > 0:
  452. # ### NEW/MODIFIED LOGIC FOR BRANCH POINT DETECTION ###
  454. # 1. Find all candidate branch points (original method)
  455. neighbor_kernel = np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]], dtype=np.uint8)
  456. num_neighbors = ndi.convolve(skeleton_filtered.astype(np.uint8), neighbor_kernel, mode='constant', cval=0)
  457. branch_points_candidates_mask = (num_neighbors > 2) & skeleton_filtered
  459. # 2. Consolidate candidate clusters into single points
  460. num_branch_points = 0
  461. branch_points_mask = np.zeros_like(skeleton_filtered, dtype=bool)
  463. if np.any(branch_points_candidates_mask):
  464. # Convert the minimum distance from micrometers to pixels
  465. min_dist_px = max(1, int(round(MIN_BRANCH_POINT_DISTANCE_UM / pixel_size_um)))
  466. print(f" -> Merging branch points within {MIN_BRANCH_POINT_DISTANCE_UM}µm ({min_dist_px}px).")
  468. # Create a distance map where clusters of candidates become "hills"
  469. distance_from_candidates = ndi.distance_transform_edt(branch_points_candidates_mask)
  471. # Find the peak of each "hill", enforcing the minimum distance
  472. coords = peak_local_max(distance_from_candidates,
  473. min_distance=min_dist_px,
  474. exclude_border=False)
  476. # The number of coordinates found is our new branch point count
  477. num_branch_points = coords.shape[0]
  479. # Create the final mask with only the consolidated points for plotting
  480. if num_branch_points > 0:
  481. branch_points_mask[tuple(coords.T)] = True
  483. branch_point_density = num_branch_points / total_roi_area_mm2 if total_roi_area_mm2 > 0 else 0
  490. results.update({'filtered_skeleton_pixels': filtered_skeleton_pixels, 'filtered_skeleton_length_um': filtered_skeleton_length_um,
  491. 'segments_analyzed': segments_analyzed, 'segments_kept': segments_kept, 'filtered_branch_points_count': num_branch_points,
  492. 'filtered_branch_point_density_mm2': branch_point_density})
  493. num_nuclei = 0; nuclear_density = 0; nuclei_labels = None
  494. if dapi_img is not None:
  495. dapi_roi = dapi_img.copy().astype(np.float32)
  496. dapi_roi[~roi_mask] = 0
  497. if np.any(dapi_roi):
  498. dapi_processed = dapi_roi
  499. if dapi_median_radius_px > 0:
  500. dapi_processed = filters.median(dapi_processed, morphology.disk(dapi_median_radius_px))
  501. dapi_processed[~roi_mask] = 0
  502. dapi_pixels_in_roi = dapi_processed[roi_mask]
  503. if np.any(dapi_pixels_in_roi > 1e-6):
  504. significant_dapi_pixels = dapi_pixels_in_roi[dapi_pixels_in_roi > np.percentile(dapi_pixels_in_roi, 5)]
  505. if len(np.unique(significant_dapi_pixels)) > 1:
  506. nuclei_thresh_val = filters.threshold_otsu(significant_dapi_pixels)
  507. nuclei_thresholded_mask = (dapi_processed > nuclei_thresh_val) & roi_mask
  508. if np.any(nuclei_thresholded_mask):
  509. distance_dapi = ndi.distance_transform_edt(nuclei_thresholded_mask)
  510. coords = peak_local_max(distance_dapi, min_distance=nuclei_min_distance_peaks_px,
  511. labels=measure.label(nuclei_thresholded_mask), exclude_border=False)
  512. if coords.size > 0:
  513. mask_peaks = np.zeros(distance_dapi.shape, dtype=bool)
  514. mask_peaks[tuple(coords.T)] = True
  515. markers, _ = ndi.label(mask_peaks)
  516. nuclei_labels = segmentation.watershed(-distance_dapi, markers, mask=nuclei_thresholded_mask)
  517. nuclei_labels[~roi_mask] = 0
  518. num_nuclei = len(np.unique(nuclei_labels)) -1
  519. nuclear_density = num_nuclei / total_roi_area_mm2 if total_roi_area_mm2 > 0 else 0
  520. results.update({'nuclei_count': num_nuclei, 'nuclear_density_mm2': nuclear_density})
  521. # Plotting code remains the same...
  522. if SAVE_PLOTS:
  523. try:
  524. fig, axes = plt.subplots(2, 3, figsize=(18, 10))
  525. plot_title = (f"{os.path.basename(image_path)} | Sig={param_sigma_start_um}-{param_sigma_stop_um}s{param_sigma_step_um}µm, "
  526. f"B={adaptive_block_size_um}µm, O={adaptive_offset}, "
  527. f"MinWidth={MIN_VESSEL_WIDTH_UM}µm, MinLen={MIN_VESSEL_SEGMENT_LENGTH_UM}µm\n"
  528. f"ROI Area: {total_roi_area_um2:.1f} µm² ({total_roi_area_mm2:.3f} mm²)")
  529. fig.suptitle(plot_title, fontsize=10)
  530. ax = axes.ravel()
  531. def plot_img(ax_plot, img_data, title, cmap='gray', vmin=None, vmax=None, roi_m=None):
  532. if img_data is None:
  533. ax_plot.set_title(title + ' (N/A)'); ax_plot.axis('off'); return None
  534. effective_roi = roi_m if roi_m is not None else np.ones_like(img_data, dtype=bool)
  535. img_to_plot = img_data.copy().astype(float)
  536. img_to_plot[~effective_roi] = np.nan
  537. valid_pixels = img_data[effective_roi & np.isfinite(img_data)]
  538. vmin_calc, vmax_calc = (np.percentile(valid_pixels, [1, 99]) if valid_pixels.size > 10 else (0, 1))
  539. ax_plot.imshow(img_to_plot, cmap=cmap, vmin=vmin_calc, vmax=vmax_calc, interpolation='nearest')
  540. ax_plot.set_title(title); ax_plot.axis('off')
  542. plot_img(ax[0], cd31_img, 'CD31 (Full)', cmap='viridis', roi_m=roi_mask)
  543. plot_img(ax[1], dapi_img, 'DAPI (Full)', cmap='magma', roi_m=roi_mask)
  544. plot_img(ax[2], vesselness_image_roi, 'Frangi Vesselness (ROI)', cmap='hot', roi_m=roi_mask)
  546. overlay = label2rgb(vessel_mask.astype(int), image=rescale_intensity(cd31_img, out_range=(0,1)), bg_label=0, alpha=0.3, colors=['red'])
  547. ax[3].imshow(overlay); ax[3].set_title(f'Vessel Mask (VAF={vaf:.1f}%)'); ax[3].axis('off')
  550. # --- Define visualization parameters ---
  551. SKELETON_LINE_THICKNESS_PX = 2 # <-- Adjust for skeleton thickness
  552. BRANCH_POINT_MARKER_RADIUS_PX = 4 # <-- Adjust for branch point size
  553. # print(skeleton_raw)
  554. # --- Plot 4: Raw Skeleton Overlay (Thickened) ---
  555. background_raw = rescale_intensity(cd31_img, out_range=(0,1))
  556. overlay_raw = np.stack([background_raw]*3, axis=-1)
  557. if skeleton_raw is not None:
  558. dilated_skeleton_raw = morphology.binary_dilation(skeleton_raw, footprint=morphology.disk(SKELETON_LINE_THICKNESS_PX))
  559. overlay_raw[dilated_skeleton_raw] = [0.0, 1.0, 1.0] # Cyan
  560. ax[4].imshow(overlay_raw)
  561. ax[4].set_title(f'RAW Skel (L={raw_skeleton_length_um:.0f}µm)')
  562. ax[4].axis('off')
  564. # --- Plot 5: Filtered Skeleton & Branch Points (Thickened) ---
  565. background_filt = rescale_intensity(cd31_img, out_range=(0,1))
  566. overlay_filt = np.stack([background_filt]*3, axis=-1)
  567. # Draw skeleton first (lime green)
  568. if skeleton_filtered is not None:
  569. dilated_skeleton_filt = morphology.binary_dilation(skeleton_filtered, footprint=morphology.disk(SKELETON_LINE_THICKNESS_PX))
  570. overlay_filt[dilated_skeleton_filt] = [0.0, 1.0, 0.0] # Lime Green
  571. # Draw branch points on top (red)
  572. if branch_points_mask is not None and np.any(branch_points_mask):
  573. dilated_branch_mask = morphology.binary_dilation(branch_points_mask, footprint=morphology.disk(BRANCH_POINT_MARKER_RADIUS_PX))
  574. overlay_filt[dilated_branch_mask] = [1.0, 0.0, 0.0] # Red
  575. ax[5].imshow(overlay_filt)
  576. ax[5].set_title(f'FILTERED Skel (L={filtered_skeleton_length_um:.0f}µm) & Br ({num_branch_points})')
  577. ax[5].axis('off')
  580. # if branch_points_mask is not None and np.any(branch_points_mask): skel_overlay_filt[branch_points_mask] = (1,0,0)
  581. # ax[5].imshow(skel_overlay_filt); ax[5].set_title(f'FILTERED Skel (L={filtered_skeleton_length_um:.0f}µm) & Br ({num_branch_points})'); ax[5].axis('off')
  585. plt.tight_layout(rect=[0, 0.03, 1, 0.95])
  586. offset_str = f"{adaptive_offset:.4f}".replace('.', 'p').replace('-', 'neg')
  587. min_len_str = f"minLen{MIN_VESSEL_SEGMENT_LENGTH_UM:.1f}um".replace('.','p') if MIN_VESSEL_SEGMENT_LENGTH_UM is not None and MIN_VESSEL_SEGMENT_LENGTH_UM > 0 else "minLenOff"
  588. min_width_str = f"minWidth{MIN_VESSEL_WIDTH_UM:.1f}um".replace('.','p') if MIN_VESSEL_WIDTH_UM is not None and MIN_VESSEL_WIDTH_UM > 0 else "minWidthOff"
  589. sigma_fn_part = f"sigmas_{param_sigma_start_um}-{param_sigma_stop_um}s{param_sigma_step_um}um".replace('.','p')
  590. block_fn_part = f"block{adaptive_block_size_um}um".replace('.', 'p')
  591. output_fig_path = os.path.join(OUTPUT_PLOT_FOLDER, f"{os.path.splitext(os.path.basename(image_path))[0]}_{sigma_fn_part}_{block_fn_part}_offset{offset_str}_{min_width_str}_{min_len_str}_analysis.png")
  592. plt.savefig(output_fig_path, dpi=450)
  593. plt.close(fig)
  594. except Exception as plot_err:
  595. print(f"!!!! Warning: Could not generate plot for {run_id}. Error: {plot_err} !!!!")
  596. traceback.print_exc()
  597. if plt.gcf().get_axes(): plt.close(plt.gcf())
  598. except Exception as e:
  599. print(f"!!!! TASK FAILED: {run_id}. Error: {e} !!!!")
  600. traceback.print_exc()
  601. results['error_message'] = f"Task Failed: {str(e)}"
  602. end_time = time.time()
  603. status = "FAILED" if results['error_message'] else "finished"
  604. print(f"--- Task {run_id} {status} in {end_time - start_time:.2f} seconds ---")
  605. return results
  610. def load_image(path):
  611. """Loads a bio-image (TIFF, CZI, etc.) using aicsimageio and extracts pixel size."""
  612. try:
  613. img_aics = AICSImage(path)
  614. dims_order = img_aics.dims.order
  615. has_c = 'C' in dims_order
  616. has_y = 'Y' in dims_order
  617. has_x = 'X' in dims_order
  619. if has_c and has_y and has_x:
  620. image_stack = img_aics.get_image_data("CYX", T=0, Z=0)
  621. elif has_y and has_x:
  622. print(f"Warning: Image {os.path.basename(path)} seems 2D (missing C dimension). Loading as YX.")
  623. image_stack = img_aics.get_image_data("YX", T=0, Z=0)
  624. else:
  625. print(f"Error: Image {os.path.basename(path)} missing Y or X dimension. Cannot process.")
  626. return None, None
  628. pixel_size_xy_um = ASSUMED_PIXEL_SIZE_UM
  630. if FORCE_ASSUMED_PIXEL_SIZE:
  631. print(f"Info: FORCE_ASSUMED_PIXEL_SIZE is True. Using assumed value: {ASSUMED_PIXEL_SIZE_UM} µm.")
  632. else:
  633. pixel_size_x_um = img_aics.physical_pixel_sizes.X
  634. pixel_size_y_um = img_aics.physical_pixel_sizes.Y
  636. if pixel_size_x_um is None or pixel_size_y_um is None:
  637. print(f"Warning: Could not read pixel size from image metadata (e.g., common for standard TIFF files). Using assumed value: {pixel_size_xy_um} µm.")
  638. elif not np.isclose(pixel_size_x_um, pixel_size_y_um):
  639. print(f"Warning: Anisotropic pixel size in {os.path.basename(path)} metadata (X: {pixel_size_x_um}, Y: {pixel_size_y_um}). Using assumed value: {pixel_size_xy_um} µm.")
  640. else:
  641. # Only use the file's pixel size if it's not being forced and it's valid.
  642. pixel_size_xy_um = pixel_size_x_um
  644. print(f"Image loaded: {os.path.basename(path)}, shape={image_stack.shape}, dtype={image_stack.dtype}, using pixel_size={pixel_size_xy_um:.4f} µm")
  645. return image_stack, pixel_size_xy_um
  646. except FileNotFoundError:
  647. print(f"Error: Image file not found at {path}")
  648. return None, None
  649. except Exception as e:
  650. print(f"Error loading image {path}: {e}")
  651. traceback.print_exc()
  652. return None, None
  656. # --- Main Execution Logic (Updated for Micrometer Parameters) ---
  657. if __name__ == "__main__":
  658. all_results = []
  659. master_start_time = time.time()
  661. # ### CHANGED ### - Search for .tif and .tiff files instead of .czi
  662. image_files = sorted([f for f in os.listdir(DATA_FOLDER) if f.lower().endswith(('.tif', '.tiff'))])
  663. if not image_files:
  664. print(f"Error: No .tif or .tiff files found in {DATA_FOLDER}")
  665. exit()
  667. # ### CHANGED ### - Updated print statement
  668. print(f"Found {len(image_files)} image files. Generating analysis tasks...")
  670. tasks_to_run = []
  671. # ### CHANGED ### - Loop over `image_files` instead of `czi_files` and handle missing ROIs
  672. for image_file in image_files:
  673. base_name = os.path.splitext(image_file)[0]
  674. image_path = os.path.join(DATA_FOLDER, image_file)
  676. # Check for case-insensitive geojson file
  677. geojson_path_options = [
  678. os.path.join(DATA_FOLDER, base_name + '.geojson'),
  679. os.path.join(DATA_FOLDER, base_name + '.GEOJSON')
  680. ]
  681. geojson_path = None
  682. for path_option in geojson_path_options:
  683. if os.path.exists(path_option):
  684. geojson_path = path_option
  685. break
  687. # ### CHANGED ### - This logic now adds a task REGARDLESS of whether an ROI is found.
  688. # If geojson_path is None, the analysis function will use the full image.
  689. if geojson_path is None:
  690. print(f"--- NOTE: No matching GeoJSON for '{base_name}'. Will analyze full image. ---")
  692. # Iterate over the new micrometer-based parameter lists
  693. param_combinations = list(itertools.product(
  694. VESSELNESS_SIGMA_RANGES_TO_TEST_UM,
  695. BLOCK_SIZES_TO_TEST_UM,
  696. OFFSETS_TO_TEST
  697. ))
  698. for sigma_range_um, block_size_um, offset in param_combinations:
  699. # ### CHANGED ### - Append `image_path` to tasks, and `geojson_path` can now be None
  700. tasks_to_run.append((image_path, geojson_path, sigma_range_um, block_size_um, offset))
  702. if not tasks_to_run:
  703. print("Error: No analysis tasks were generated. Check for images in the folder and parameter settings.")
  704. else:
  705. num_tasks = len(tasks_to_run)
  706. num_images = len(image_files)
  707. print(f"Generated {num_tasks} analysis tasks from {num_images} image file(s).")
  708. print(f"Starting parallel execution with {NUM_WORKERS} workers...")
  710. completed_count = 0
  711. with concurrent.futures.ProcessPoolExecutor(max_workers=NUM_WORKERS) as executor:
  712. # ### CHANGED ### - Call the renamed `analyze_image_roi` function
  713. future_to_task_args = {executor.submit(analyze_image_roi, *args): args for args in tasks_to_run}
  714. for future in concurrent.futures.as_completed(future_to_task_args):
  715. task_args = future_to_task_args[future]
  716. try:
  717. result_data = future.result()
  718. all_results.append(result_data)
  719. except Exception as exc:
  720. print(f'\n!!!! Task for {os.path.basename(task_args[0])} with params (sigma_range_um={task_args[2]}, block_um={task_args[3]}, offset={task_args[4]}) failed unexpectedly at executor level: {exc} !!!!')
  721. traceback.print_exc()
  722. (failed_sigma_start, failed_sigma_stop, failed_sigma_step) = task_args[2]
  723. # ### CHANGED ### - Update dictionary key for failed tasks
  724. all_results.append({
  725. 'filename_image': os.path.basename(task_args[0]),
  726. 'filename_geojson': os.path.basename(task_args[1]) if task_args[1] else 'N/A (Full Image)',
  727. 'vesselness_sigma_start_um': failed_sigma_start,
  728. 'vesselness_sigma_stop_um': failed_sigma_stop,
  729. 'vesselness_sigma_step_um': failed_sigma_step,
  730. 'adaptive_block_size_um': task_args[3],
  731. 'adaptive_offset': task_args[4],
  732. 'min_vessel_width_um_param': MIN_VESSEL_WIDTH_UM,
  733. 'min_vessel_segment_length_um_param': MIN_VESSEL_SEGMENT_LENGTH_UM,
  734. 'error_message': f'Executor level failure: {exc}'
  735. })
  736. finally:
  737. completed_count += 1
  738. print(f"Progress: {completed_count}/{num_tasks} tasks completed.", end='\r' if completed_count < num_tasks else '\n')
  739. print(f"Parallel execution finished. Processed {completed_count} tasks.")
  741. if all_results:
  742. print(f"\n--- Consolidating {len(all_results)} results for CSV output. ---")
  743. results_df = pd.DataFrame(all_results)
  745. # ### CHANGED ### - Update column name for output CSV
  746. column_order = [
  747. 'filename_image', 'filename_geojson',
  748. 'vesselness_sigma_start_um', 'vesselness_sigma_stop_um', 'vesselness_sigma_step_um',
  749. 'adaptive_block_size_um', 'adaptive_offset',
  750. 'min_vessel_width_um_param', 'min_vessel_segment_length_um_param',
  751. 'pixel_size_um',
  752. 'roi_area_pixels', 'roi_area_um2', 'roi_area_mm2',
  753. 'vessel_area_pixels', 'vessel_area_um2', 'vaf_percent',
  754. 'mean_vessel_diameter_um',
  755. 'raw_skeleton_pixels', 'raw_skeleton_length_um',
  756. 'segments_analyzed', 'segments_kept',
  757. 'filtered_skeleton_pixels', 'filtered_skeleton_length_um',
  758. 'filtered_branch_points_count', 'filtered_branch_point_density_mm2',
  759. 'nuclei_count', 'nuclear_density_mm2',
  760. 'error_message'
  761. ]
  763. for col in column_order:
  764. if col not in results_df.columns:
  765. results_df[col] = pd.NA
  767. results_df = results_df.reindex(columns=column_order)
  768. output_csv_path = os.path.join(DATA_FOLDER, OUTPUT_CSV_FILENAME)
  769. try:
  770. # ### CHANGED ### - Update sort key
  771. results_df.sort_values(
  772. by=['filename_image', 'vesselness_sigma_start_um', 'vesselness_sigma_step_um', 'adaptive_block_size_um', 'adaptive_offset'],
  773. inplace=True, na_position='last'
  774. )
  775. results_df.to_csv(output_csv_path, index=False, float_format='%.5f')
  776. print(f"Results successfully saved to: {output_csv_path}")
  777. except Exception as e:
  778. print(f"!!!! Error saving CSV file to {output_csv_path}: {e} !!!!")
  779. traceback.print_exc()
  780. else:
  781. print("\nNo results were generated or collected.")
  783. master_end_time = time.time()
  784. print(f"\nTotal execution time: {master_end_time - master_start_time:.2f} seconds.")
  785. print("--- Script Finished ---")

adaptive_vascular_analysis.py, under CC-BY-4.0 · at the source

Overview

Authors: Alexander King1,2, Catherine Garcia1, Crisylle Blanton1,2, Anna Chen1, Amna Ahmad1, David Lukacsovich3, Csaba Földy3, Takako Makita4, Garret R Anderson1
ORCID iDs: Alexander King, Garret R Anderson
  1. Department of Molecular, Cell, and Systems Biology, University of California – Riverside, Riverside, California 92521
  2. Neuroscience Graduate Program, University of California – Riverside, Riverside, California 92521
  3. Laboratory of Neural Connectivity, Brain Research Institute, Faculties of Medicine and Science, University of Zürich, Zurich, Switzerland
  4. Department of Regenerative Medicine and Cell Biology, Medical University of South Carolina, Charleston, South Carolina 29425
Institutions: University of California, Riverside (United States); University of Zurich (Switzerland); Medical University of South Carolina (United States)
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience, volume 46, issue 16, article e0019262026
Dates: received 5 January 2026; accepted 12 March 2026; published online 24 March 2026; in print 22 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1523/jneurosci.0019-26.2026 · PMID 41876233 · PMCID PMC13108381 · OpenAlex W7140218929
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism), stroke (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Evoked potentials, fMRI & imaging
Keywords: adhesion G-protein coupled receptors, alternative splicing, cell–cell recognition, cerebrovasculature, latrophilin
MeSH: Alternative Splicing*, Cerebrovascular Circulation*, Endothelial Cells*, Receptors, G-Protein-Coupled*, Animals, Female, Male, Mice, Mice, Inbred C57BL, Mice, Knockout, Neurons (* major topic)
Topic: Single-cell and spatial transcriptomics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Whitehall Foundation (2018-08-01)
Citations: not cited yet (Europe PMC); 73 references in the paper

Abstract

Central nervous system development requires parallel but interrelated processes of neural circuit assembly and vascularization. Intersecting between these two processes is the cell-adhesion G-protein coupled receptor Adgrl2. In select neuronal populations, Adgrl2 is localized and control the assembly of specific synaptic sites. In non-neuronal brain cells, Adgrl2 is restricted in expression to endothelial cells. Testing for Adgrl2 function in these cells in mice (of either sex), here we find that endothelial cell-specific Adgrl2 deletion results in an impairment in cerebrovascular integrity. To understand how it might be possible for Adgrl2 to function independently in neuronal and endothelial contexts, we surveyed Adgrl2 transcripts within these cell classes. By analyzing single-cell RNA sequencing datasets, we find that Adgrl2 mRNA is subject to robust cell type-specific alternative splicing that results in distinct isoforms being produced in neurons compared with endothelial cells. To probe the functional significance of this alternative splicing, we forced expression of the neuronal isoform of Adgrl2 in endothelial cells. This resulted in altered cerebrovascular properties including the formation of ectopic glutamatergic synaptic contacts onto endothelial cells, indicating alterations in the cell–cell recognition process. Functionally, in direct contrast to endothelial Adgrl2 deletion, this genetic expression switch instead enhances blood–brain barrier integrity. This overly restrictive cerebrovascular function results in dysregulation of blood to cerebrospinal fluid homeostasis, enlargement of brain ventricles, and a higher risk of hydrocephalus. Thus, alternative splicing serves as a cell type-specific mechanism that provides isoform-specific Adgrl2 for discerning functions controlling neural circuit assembly and cerebrovascular homeostasis.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

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

Zenodo 15802164

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Languages: Python (1)
Size: 1 file, 1 script
Software Heritage: not checked
Found in: “Data Availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), scikit-image (1 file), SciPy (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
1 file
At the source: doi.org/10.5281/zenodo.15802164

AllenCellModeling/aicsimageio

License: BSD-3-Clause
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 8dfaaba9893c0ee5931938a0def3a4d93a5d13eb, 1 December 2025
Languages: Python (76), Jupyter (2)
Size: 118 files, 78 scripts
Software Heritage: not archived
Found in: the text, “Vascular morphology analysis”
Holds: README, license file, environment (setup.cfg, setup.py, presentations/2021-dask-life-sciences/environment.yml), tests, continuous integration, documentation, 2 notebooks
Not found: CITATION.cff
Tools: NumPy (39 files), xarray (22 files), tifffile (6 files), imageio (3 files), pandas (2 files), Matplotlib (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
80 files
At the source: github.com/AllenCellModeling/aicsimageio

shapely/shapely

License: BSD-3-Clause
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 30c0a51039d4b4b0b3b74ba117039e2c5849906b, 29 September 2026
Languages: Python (158), C/C++ (32), C (31), Shell (1)
Size: 293 files, 222 scripts
Software Heritage: archived
Found in: the text, “Vascular morphology analysis”
Holds: README, license file, CITATION.cff, environment (pyproject.toml, ci/Dockerfile, docker/Dockerfile.arm64, docker/Dockerfile.valgrind, docs/environment.yml), tests, continuous integration, documentation
Tools: NumPy (64 files), Matplotlib (29 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
224 files
At the source: github.com/shapely/shapely

Zenodo 15802163

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Languages: Python (1)
Size: 1 file, 1 script
Software Heritage: not checked
Found in: DataCite
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), scikit-image (1 file), SciPy (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
1 file
At the source: doi.org/10.5281/zenodo.15802163

The paper's code and data availability statement is in the Data section.

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:

  • 4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 302 scripts, each with its path and the digest of its content;
  • 2 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

Datasets cited

Data Availability

Raw sequencing data used in this study are accessible from NCBI Gene Expression Omnibus database through accession numbers GSE185862 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE185862), GSE133291 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133291), GSE121653 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121653), GSE98816 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98816), GSE99235 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99235), and GSE159812 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159812); or Synapse.org (syn51015750). Raw and analyzed data will be provided upon request. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request (). Publicly available code used in this study are accessible as described. Unique code generated in this study for vascular analysis is available at Zenodo (https://doi.org/10.5281/zenodo.15802164).

Reproduced under the paper's license (CC BY), from the paper cited above.

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, 9 authors, 5 keywords, 11 MeSH terms, 1 funder, 72 references.

Cite

This paper

King, A., Garcia, C., Blanton, C., Chen, A., Ahmad, A., Lukacsovich, D., Földy, C., Makita, T., & Anderson, G. R. (2026). Endothelial &lt;i&gt;Adgrl2&lt;/i&gt; Expression and Alternative Splicing Controls the Cerebrovasculature. The Journal of neuroscience : the official journal of the Society for Neuroscience, 46(16), e0019262026. https://doi.org/10.1523/jneurosci.0019-26.2026

BibTeX

@article{king2026endothelial,
author = {King, Alexander and Garcia, Catherine and Blanton, Crisylle and Chen, Anna and Ahmad, Amna and Lukacsovich, David and Földy, Csaba and Makita, Takako and Anderson, Garret R},
title = {{Endothelial \&lt;i\&gt;Adgrl2\&lt;/i\&gt; Expression and Alternative Splicing Controls the Cerebrovasculature}},
journal = {The Journal of neuroscience : the official journal of the Society for Neuroscience},
year = {2026},
month = apr,
volume = {46},
number = {16},
pages = {e0019262026},
publisher = {Society for Neuroscience},
issn = {0270-6474},
doi = {10.1523/jneurosci.0019-26.2026},
url = {https://doi.org/10.1523/jneurosci.0019-26.2026},
pmid = {41876233},
pmcid = {PMC13108381}
}

RIS

TY - JOUR
AU - King, Alexander
AU - Garcia, Catherine
AU - Blanton, Crisylle
AU - Chen, Anna
AU - Ahmad, Amna
AU - Lukacsovich, David
AU - Földy, Csaba
AU - Makita, Takako
AU - Anderson, Garret R
TI - Endothelial &lt;i&gt;Adgrl2&lt;/i&gt; Expression and Alternative Splicing Controls the Cerebrovasculature
T2 - The Journal of neuroscience : the official journal of the Society for Neuroscience
J2 - J Neurosci
PY - 2026
DA - 2026/04/22
VL - 46
IS - 16
SP - e0019262026
SN - 0270-6474
PB - Society for Neuroscience
DO - 10.1523/jneurosci.0019-26.2026
UR - https://doi.org/10.1523/jneurosci.0019-26.2026
LA - en
ER -

CSL-JSON

{
"id": "10.1523/jneurosci.0019-26.2026",
"type": "article-journal",
"title": "Endothelial &lt;i&gt;Adgrl2&lt;/i&gt; Expression and Alternative Splicing Controls the Cerebrovasculature",
"container-title": "The Journal of neuroscience : the official journal of the Society for Neuroscience",
"author": [
{
"family": "King",
"given": "Alexander"
},
{
"family": "Garcia",
"given": "Catherine"
},
{
"family": "Blanton",
"given": "Crisylle"
},
{
"family": "Chen",
"given": "Anna"
},
{
"family": "Ahmad",
"given": "Amna"
},
{
"family": "Lukacsovich",
"given": "David"
},
{
"family": "Földy",
"given": "Csaba"
},
{
"family": "Makita",
"given": "Takako"
},
{
"family": "Anderson",
"given": "Garret R"
}
],
"container-title-short": "J Neurosci",
"volume": "46",
"issue": "16",
"page": "e0019262026",
"DOI": "10.1523/jneurosci.0019-26.2026",
"PMID": "41876233",
"PMCID": "PMC13108381",
"ISSN": "0270-6474",
"publisher": "Society for Neuroscience",
"URL": "https://doi.org/10.1523/jneurosci.0019-26.2026",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
22
]
]
}
}

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-71619-1
Structurally exclusive Teneurin complexes orchestrate divergent programs in early cortical development.
Journal: Nature communications
In common: mouse, 10 references
[2] doi:10.1038/s41467-026-73373-w [code]
Mapping neuro-vascular unit communications reveals distinct angiogenic programs across developing mouse brain regions.
Journal: Nature communications
In common: imageio, tifffile, scikit-image, 4 other tools, mouse, 3 references
[3] doi:10.1016/j.cub.2026.04.038 [code]
Ten3-Lphn2-mediated target selection across the extended hippocampal network demonstrates a repeated strategy for circuit assembly.
Journal: Current biology : CB
In common: mouse, 7 references
[4] doi:10.1016/j.isci.2026.116213 [code]
Opioid receptor distribution in the claustrum-dorsal endopiriform complex.
Journal: iScience
In common: xarray, imageio, tifffile, 4 other tools, cellular / molecular, 1 reference
[5] doi:10.1038/s41586-026-10679-1 [code]
Cortical development dynamics across autism spectrum disorder mouse models.
Journal: Nature
In common: imageio, tifffile, scikit-image, 4 other tools, mouse, cellular / molecular, 2 references
[6] doi:10.1371/journal.pcbi.1014571 [code]
SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: imageio, tifffile, scikit-image, 4 other tools, 2 references
[7] doi:10.1016/j.celrep.2026.117255
The atypical adhesion GPCR ADGRA1 controls hippocampal inhibitory circuit function.
Journal: Cell reports
In common: cellular / molecular, 6 references
[8] doi:10.1016/j.isci.2026.115689 [code]
Unperturbed dye-based imaging of spontaneous synchronized calcium activity in iPSC-derived neuronal cultures.
Journal: iScience
In common: xarray, imageio, tifffile, 4 other tools, cellular / molecular
[9] doi:10.7554/elife.109717 [code]
Retrosplenial cortex enables context-dependent goal-directed sensorimotor transformation.
Journal: eLife
In common: imageio, tifffile, scikit-image, 4 other tools, mouse, 1 reference
[10] doi:10.3389/fnbeh.2026.1895371 [code]
Discrete and differentiated encoding of distinct components of spatial experience occurs within the proximodistal subfields of the dorsoventral axis of the hippocampus.
Journal: Frontiers in behavioral neuroscience
In common: xarray, imageio, tifffile, 4 other tools

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.