Punctate analysis methods

Author: mogres

.python-version

@@ -0,0 +1 @@
+3.9

__init__.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 
 """Top-level package for cellpack_analysis."""
 

cellpack_analysis/analysis/occupancy_analysis/calculate_available_space.py

@@ -9,63 +9,65 @@
 distances are saved in a grid directory for each cellid.
 Distances are normalized by the cell diameter and saved in the grid directory.
 """
+import logging
 from concurrent.futures import ProcessPoolExecutor, as_completed
-from pathlib import Path
 
 import pandas as pd
 from tqdm import tqdm
 
+from cellpack_analysis.lib.file_io import get_project_root
 from cellpack_analysis.lib.mesh_tools import calculate_grid_distances
 
+log = logging.getLogger(__name__)
 # %% set structure id
-STRUCTURE_ID = "RAB5A"  # SLC25A17: peroxisomes, RAB5A: early endosomes
+STRUCTURE_ID = "mean"  # SLC25A17: peroxisomes, RAB5A: early endosomes
 # %% set up parameters for grid
 SPACING = 2
 # %% set file paths and setup parameters
-base_datadir = Path(__file__).parents[3] / "data"
-print(f"Data directory: {base_datadir}")
+base_datadir = get_project_root() / "data"
+log.info(f"Data directory: {base_datadir}")
 
 # %% select cellids to use
 use_mean_shape = True
 if use_mean_shape:
-    mesh_folder = base_datadir / "average_shape_meshes"
     cellids_to_use = ["mean"]
 else:
-    mesh_folder = base_datadir / f"structure_data/{STRUCTURE_ID}/meshes/"
     df_cellid = pd.read_csv("s3://cellpack-analysis-data/all_cellids.csv")
     df_struct = df_cellid.loc[df_cellid["structure_name"] == STRUCTURE_ID]
-    cellids_to_use = df_struct.loc[df_struct["8dsphere"], "CellId"]
-print(f"Using {len(cellids_to_use)} cellids")
+    cellids_to_use = df_struct.loc[df_struct["8dsphere"], "CellId"].tolist()
+mesh_folder = base_datadir / f"structure_data/{STRUCTURE_ID}/meshes/"
+log.info(f"Using {len(cellids_to_use)} cellids")
 # %% get meshes for cellids used
+# cellids_to_use = [cellids_to_use[0]]
 cellid_list = []
 nuc_meshes_to_use = []
 mem_meshes_to_use = []
 for cellid in cellids_to_use:
-    nuc_mesh = mesh_folder / f"nuc_mesh_{cellid}.obj"
-    mem_mesh = mesh_folder / f"mem_mesh_{cellid}.obj"
-    if nuc_mesh.exists() and mem_mesh.exists():
+    nuc_mesh_path = mesh_folder / f"nuc_mesh_{cellid}.obj"
+    mem_mesh_path = mesh_folder / f"mem_mesh_{cellid}.obj"
+    if nuc_mesh_path.exists() and mem_mesh_path.exists():
         cellid_list.append(cellid)
-        nuc_meshes_to_use.append(nuc_mesh)
-        mem_meshes_to_use.append(mem_mesh)
-print(f"Found {len(nuc_meshes_to_use)} meshes")
+        nuc_meshes_to_use.append(nuc_mesh_path)
+        mem_meshes_to_use.append(mem_mesh_path)
+log.info(f"Found {len(nuc_meshes_to_use)} meshes")
 # %% set up grid results directory
-grid_dir = base_datadir / f"structure_data/{STRUCTURE_ID}/grid_distances"
+grid_dir = base_datadir / f"structure_data/{STRUCTURE_ID}/grid_distances/"
 grid_dir.mkdir(exist_ok=True, parents=True)
-# %% run in parallel
-PARALLEL = True
-recalculate = True
+# %% run the workflow
+save_dir = None
+PARALLEL = False
+recalculate = False
 calc_nuc_distances = True
 calc_mem_distances = True
 calc_z_distances = True
-calc_inside_mem = True
-chunk_size = 50000
+calc_scaled_nuc_distances = True
+chunk_size = 20000
 if PARALLEL:
-    num_cores = 4
+    num_cores = 8
     results = []
     with ProcessPoolExecutor(max_workers=num_cores) as executor:
-        futures = []
         for i in range(len(nuc_meshes_to_use)):
-            futures.append(
+            results.append(
                 executor.submit(
                     calculate_grid_distances,
                     nuc_meshes_to_use[i],
@@ -77,13 +79,15 @@
                     calc_nuc_distances,
                     calc_mem_distances,
                     calc_z_distances,
+                    calc_scaled_nuc_distances,
                     chunk_size,
                 )
             )
 
-        with tqdm(total=len(futures), desc="CellIDs done") as pbar:
-            for future in as_completed(futures):
-                pbar.update(1)
+        with tqdm(total=len(results), desc="CellIDs done") as pbar:
+            for result in as_completed(results):
+                if result.result:
+                    pbar.update(1)
 else:
     results = []
     for i in tqdm(range(len(nuc_meshes_to_use)), desc="CellIDs done"):
@@ -98,6 +102,9 @@
                 calc_nuc_distances=calc_nuc_distances,
                 calc_mem_distances=calc_mem_distances,
                 calc_z_distances=calc_z_distances,
+                calc_scaled_nuc_distances=calc_scaled_nuc_distances,
                 chunk_size=chunk_size,
             )
         )
+
+# %%

cellpack_analysis/analysis/occupancy_analysis/visualize_available_space_distribution.py

@@ -1,16 +1,13 @@
 # %% [markdown]
-"""
-## Analyze and visualize available space distribution
-"""
+# # Analyze and visualize available space distribution
+import logging
 import pickle
 from pathlib import Path
 
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 import seaborn as sns
-
-# %%
 import trimesh
 from tqdm import tqdm
 
@@ -19,16 +16,23 @@
     round_away_from_zero,
 )
 
-# %% set pixel size
+log = logging.getLogger(__name__)
+
+# %% [markdown]
+# ## Set parameters and file paths
+MEAN_SHAPE = False
+
 PIX_SIZE = 0.108  # um per pixel
 
-# %% set structure id
-STRUCTURE_ID = "SLC25A17"  # SLC25A17: peroxisomes, RAB5A: early endosomes
-# %% set file paths and setup parameters
+STRUCTURE_ID = "SEC61B"  # SLC25A17: peroxisomes, RAB5A: early endosomes
+STRUCTURE_NAME = "ER_peroxisome"
+
 base_datadir = Path(__file__).parents[3] / "data"
 base_results_dir = Path(__file__).parents[3] / "results"
 
-results_dir = base_results_dir / f"stochastic_variation_analysis/{STRUCTURE_ID}/rules/"
+results_dir = (
+    base_results_dir / f"stochastic_variation_analysis/{STRUCTURE_NAME}/rules/"
+)
 results_dir.mkdir(exist_ok=True, parents=True)
 
 figures_dir = results_dir / "figures"
@@ -37,24 +41,23 @@
 grid_dir = base_datadir / f"structure_data/{STRUCTURE_ID}/grid_distances"
 grid_dir.mkdir(exist_ok=True, parents=True)
 
-print(f"Results directory: {results_dir}")
-print(f"Figures directory: {figures_dir}")
-print(f"Grid directory: {grid_dir}")
+log.info(f"Results directory: {results_dir}")
+log.info(f"Figures directory: {figures_dir}")
+log.info(f"Grid directory: {grid_dir}")
 
-# %%
-use_mean_shape = False
-
-# %% select cellids to use
-if use_mean_shape:
+# %% [markdown]
+# ## Select cellids to use
+if MEAN_SHAPE:
     mesh_folder = base_datadir / "average_shape_meshes"
     cellids_to_use = ["mean"]
 else:
     mesh_folder = base_datadir / f"structure_data/{STRUCTURE_ID}/meshes/"
     df_cellid = pd.read_csv("s3://cellpack-analysis-data/all_cellids.csv")
     df_struct = df_cellid.loc[df_cellid["structure_name"] == STRUCTURE_ID]
-    cellids_to_use = df_struct.loc[df_struct["8dsphere"], "CellId"]
-print(f"Using {len(cellids_to_use)} cellids")
-# %% get meshes for cellids used
+    cellids_to_use = df_struct.loc[df_struct["8dsphere"], "CellId"].tolist()
+log.info(f"Using {len(cellids_to_use)} cellids")
+# %% [markdown]
+# ## Get meshes for cellids used
 cellid_list = []
 nuc_meshes_to_use = []
 mem_meshes_to_use = []
@@ -65,10 +68,12 @@
         cellid_list.append(cellid)
         nuc_meshes_to_use.append(nuc_mesh)
         mem_meshes_to_use.append(mem_mesh)
-print(f"Found {len(nuc_meshes_to_use)} meshes")
+log.info(f"Found {len(nuc_meshes_to_use)} meshes")
 
-# %% load meshes
-if use_mean_shape:
+# %% [markdown]
+# # Grid spacing illustration
+# ## Load meshes
+if MEAN_SHAPE:
     nuc_mesh = trimesh.load_mesh(
         base_datadir / "average_shape_meshes/nuc_mesh_mean.obj"
     )
@@ -79,86 +84,119 @@
     nuc_mesh = trimesh.load_mesh(nuc_meshes_to_use[0])
     mem_mesh = trimesh.load_mesh(mem_meshes_to_use[0])
 
-# %% set up grid
-SPACING = 2
+# %% [markdown]
+# ## Get grid points
+SPACING = 5
 bounds = mem_mesh.bounds
 bounding_box = round_away_from_zero(bounds)
 all_points = get_list_of_grid_points(bounding_box, SPACING)
-# %% explicit inside-outside check
-print("Calculating nuc inside check")
+
+# %% [markdown]
+# ## Run inside-outside check
+log.info("Calculating nuc inside check")
 inside_nuc = nuc_mesh.contains(all_points)
-print("Calculating mem inside check")
+log.info("Calculating mem inside check")
 inside_mem = mem_mesh.contains(all_points)
-
-# %% find points inside mem but outside nuc
+# %% [markdown]
+# ## Plot grid points
 inside_mem_outside_nuc = inside_mem & ~inside_nuc
-# %% plot grid point scatter plot
+
 fig, ax = plt.subplots(dpi=300)
 all_points_scaled = all_points * PIX_SIZE
+centroid = np.mean(all_points_scaled, axis=0)
 ax.scatter(
-    all_points_scaled[inside_mem_outside_nuc, 0],
-    all_points_scaled[inside_mem_outside_nuc, 1],
+    all_points_scaled[inside_mem_outside_nuc, 0] - centroid[0],
+    all_points_scaled[inside_mem_outside_nuc, 1] - centroid[1],
     c="magenta",
-    label="inside mem outside nuc",
-    s=0.1,
-    alpha=0.7,
+    label="Available points",
+    s=0.5,
+    alpha=1,
 )
 ax.scatter(
-    all_points_scaled[inside_nuc, 0],
-    all_points_scaled[inside_nuc, 1],
+    all_points_scaled[inside_nuc, 0] - centroid[0],
+    all_points_scaled[inside_nuc, 1] - centroid[1],
     c="cyan",
     label="inside nuc",
-    s=0.1,
-    alpha=0.7,
+    s=0.5,
+    alpha=1,
 )
-ax.set_xlabel("x (\u03BCm)")
-ax.set_ylabel("y (\u03BCm)")
-ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1))
+ax.set_xlabel("x (\u03bcm)")
+ax.set_ylabel("y (\u03bcm)")
+# ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1))
 ax.set_aspect("equal")
+# ax.set_aspect(1.3)
 plt.show()
-fig.savefig(figures_dir / "grid_points.png", bbox_inches="tight")
+file_name = "grid_points"
+if MEAN_SHAPE:
+    file_name += "_mean"
+fig.savefig(figures_dir / f"{file_name}.png", bbox_inches="tight")
 
-# %% load mesh information
+# %% [markdown]
+# # Plot distance from nucleus for all shapes
+# ## Load mesh information
 file_path = grid_dir.parent / "mesh_information.dat"
 with open(file_path, "rb") as f:
     mesh_information_dict = pickle.load(f)
-# %% load saved distances
+# %% [markdown]
+# ## Load distances
 # normalization = "cell_diameter"
 normalization = None
 nuc_distances = []
 mem_distances = []
 for cellid in tqdm(cellids_to_use):
-    normalization_factor = mesh_information_dict[str(cellid)].get(normalization, 1)
+    if normalization is not None:
+        normalization_factor = mesh_information_dict[str(cellid)].get(normalization, 1)
+    else:
+        normalization_factor = 1
     nuc_distances.append(
         np.load(grid_dir / f"nuc_distances_{cellid}.npy") / normalization_factor
     )
     mem_distances.append(
         np.load(grid_dir / f"mem_distances_{cellid}.npy") / normalization_factor
     )
 
-# %% plot distance distribution kdeplot
+# %% [markdown]
+# ## Plot distance distribution as kde
 fig, ax = plt.subplots(dpi=300)
-cmap = plt.get_cmap("jet", len(nuc_distances))
+cmap = plt.get_cmap("viridis", len(nuc_distances))
+color_inds = np.random.permutation(len(nuc_distances))
 all_nuc_distances = []
 for i in tqdm(range(len(nuc_distances))):
     distances_to_plot = nuc_distances[i] * PIX_SIZE
     distances_to_plot = distances_to_plot[distances_to_plot > 0]
-    sns.kdeplot(distances_to_plot, ax=ax, color=cmap(i + 1), alpha=0.3)
+    sns.kdeplot(
+        distances_to_plot, ax=ax, color=cmap(color_inds[i]), alpha=0.4, linewidth=1
+    )
     all_nuc_distances.append(distances_to_plot)
+    # break
 all_nuc_distances = np.concatenate(all_nuc_distances)
-sns.kdeplot(all_nuc_distances, ax=ax, color=cmap(0), linewidth=2)
-mean_distance = np.mean(all_nuc_distances)
+sns.kdeplot(all_nuc_distances, ax=ax, color="k", linewidth=3)
+mean_distance = np.mean(all_nuc_distances).item()
 ax.axvline(mean_distance, color="black", linestyle="--")
-ax.set_title(f"Distance to nucleus\nMean: {mean_distance:.2f}\u03BCm")
-ax.set_xlabel("Distance (\u03BCm)")
+# %% [markdown]
+# Adjust plot
+# xlim = ax.get_xlim()
+plt.rcdefaults()
+ax.set_xlim((-1, 10))
+ax.set_title(f"Distance to nucleus\nMean: {mean_distance:.2f}\u03bcm")
+ax.set_xlabel("Distance (\u03bcm)")
 ax.set_ylabel("Probability density")
-plt.show()
-# %% plot distance distribution histogram
+file_name = "nuc_distance_kde"
+fig.savefig(figures_dir / f"{file_name}.png", bbox_inches="tight")
+fig.savefig(figures_dir / f"{file_name}.svg")
+# plt.show()
+fig
+# %% [markdown]
+# ## Plot distance distribution histogram
 fig, ax = plt.subplots(dpi=300)
 nuc_distances = np.array(nuc_distances)
 distances_to_plot = nuc_distances[0][nuc_distances[0] > 0] * PIX_SIZE
 ax.hist(distances_to_plot, bins=100)
-ax.set_xlabel("Distance (\u03BCm)")
+ax.set_xlabel("Distance (\u03bcm)")
 ax.set_ylabel("Number of points")
 ax.set_title("Distance to nucleus")
+file_name = "nuc_distance_hist"
+fig.savefig(figures_dir / f"{file_name}.png", bbox_inches="tight")
 plt.show()
+
+# %%